C     MODEL SUBROUTINE FOR SCHIZ SIMULATION 
C     OCTOBER 19, 1971 
C     BOB STOUT, PROGRAMMER 
C     D. MALIN, SIMULATION DESIGNER 
C     MODIFIED AUGUST, 1972 BY BOB STOUT 
C     MODIFIED JANUARY, 1973 BY BOB STOUT
C                 
C                 
      SUBROUTINE MODEL(N,NCOND,*) 
C                 
      INTEGER*4 SI(12),SF(12), OVNAM 
      INTEGER*4 IV(12), SINK 
      LOGICAL*4 COSTPT, DATOUT
      REAL*4 FV(12) 
      REAL*4 EXCOST(9),ECOSTF(2),SCOSTB(9),SCOSTI(2),MCOSTF(2),ACOSTF(3)
      REAL*4 OHDRT(3),SREIM,TCOST(2),TOCST(2)
      INTEGER*4 YES1,YESY,ERSP
!!      DATA YES1/'1'/,  YESY/'Y'/
      INTEGER*4 INST,SAMP,NT,EXTER,PCOST,IX,IR
      LOGICAL*4 LEAVE,SKIZ,MAXPT,PNTIT
      COMMON /VARS1/IV,FV,SI,SF 
      COMMON /IO1/NOV,SINK,OVNAM(8,6),COSTPT,DATOUT
      COMMON /IO/IDEV1,IDEV2,IDEV3,IDEV4
C                 
      INTEGER*4 REL,ADOPT,DIAG,SEXI,SEXS 
      EQUIVALENCE (INDEX,IV(1)),(REL,IV(2)),(ADOPT,IV(3)), 
     1 (DIAG,IV(4)),(SEXI,IV(5)),(SEXS,IV(6)),(METHOD,IV(7)),
     2 (INST,IV(8)),(SAMP,IV(9)),(EXTER,IV(10)),(PCOST,IV(11))
C                 
      REAL*4 METH(2), NOISE 
C     METH CONTAINS THE MAIN EFFECT PARAMETERS FOR METHOD 
C     BASE IS THE BASE RATE OF SCHIZ IN THE GENERAL POPULATION 
C     DMINT IS THE COEFFICIENT FOR THE DIAG-METHOD INTERACTION 
C                 
C     SEXEFF CONTAINS PARAMETERS FOR THE SEX EFFECT.  EFFECT DEPENDS
C     ONLY ON WHETHER SEXI(INDEX S)=SEXS(ACTUAL S).  FOR SEXI=SEXS 
C     SEXEFF(2) IS USED, OTHERWISE SEXEFF(1). 
      REAL*4 SEXEFF(2) 
C                 
C     RATES CONTAINS THE SCHIZ RATES FOR VARIOUS RELATIVES: 1=UNREL,
C     2=AUNTUNCL, 3=COUSIN, 4=SIBLING, 5=FRATTWIN, 6=IDTWIN, 7=PARENT,
C     8=CHILD, 9=NIECNEPH
      REAL*4 RATES(9) 
C                 
C     DIAGEF AND ADOEFF CONTAIN THE MAIN EFFECT PARAMETERS FOR STRICY-
C     NESS OF DIAGNOSIS AND ADOPT 
      REAL*4 DIAGEF(2), ADOEFF(3) 
C                 
C------------------------------ 
C                 
      DOUBLE PRECISION MODNM1,MODNM2,MODNM3
      COMMON /DATANM/ MODNM1,MODNM2,MODNM3
      DATA MODNM1/'SHGM.DAT'/,MODNM2/'SHGM3.DAT'/,MODNM3/'SHGM.BIN'/
      DATA YES1/'1'/,  YESY/'Y'/
      METHOD=IV(7)+SI(7)
      TOHRT=OHDRT(INST)+TOCST(SAMP)
	NOISE=NOISE!DEC-10 REQUIES THIS
	NCOND=NCOND!DEC-10 BUG MAKES THIS NECESSARY
      XN=N        
      SEX=0.      
C     COMPUTE OVERALL SEX EFFECT 
      IF(SEXI.EQ.0 .OR. SEXS.EQ.0)GO TO 2 
C                 
C     IT ONLY MATTERS WHETHER SEXI=SEXS, NOT WHAT THE SEXES ARE 
      SEX=SEXEFF(1) 
      IF(SEXI.EQ.SEXS)SEX=SEXEFF(2) 
C     IDENTICAL TWINS ALWAYS HAVE SAME SEX 
 2    IF(REL.EQ.6)SEX=SEXEFF(2) 
      IF(SEX.GT.0.01)GO TO 3 
C                 
C     SEX IS RANDOMIZED--COMPUTE AVERAGE EFFECT 
      CALL BINOM(N,.5,I) 
      XI=I        
      SEX=XI*SEXEFF(1)+(XN-XI)*SEXEFF(2) 
      SEX=SEX/XN  
C                 
C------------------------------ 
C                 
C     ***** THE MODEL ***** 
C                 
C      FIRST SET UP SOME THINGS
C                 
    3 LEAVE=.FALSE.
      P1=BASE*(1.+DIAGEF(DIAG))
C                 
C      COMPUTE BASIC COST FACTORS
C                 
      ECST=EXCOST(REL)+ECOSTF(INDEX)+TCOST(SAMP)
      SCST1=SCOSTB(REL)+SCOSTI(INDEX)
      SCST2=SCST1+MCOSTF(METHOD)*SCST1
      SCST3=SCST1+ACOSTF(ADOPT)*SCST1
      SCOST=SCST1+SCST2+SCST3+SREIM
C     SKIP THE MODEL IF INDEX=NONSCHIZ OR REL=UNREL 
      IF(INDEX.EQ.2 .OR. REL.EQ.1)GO TO 4 
C                 
C     HEREDITY-REARING FACTORS 
      P1=RATES(REL)*ADOEFF(ADOPT) 
C                 
C     SEX LINKAGE FACTOR 
      P1=P1+(SEX-1.)*P1*(1.-P1) 
C                 
C     ADD IN BASE RATE 
      P1=P1+BASE*(1.-P1) 
C                 
C     PUT IN GROSS RATE INFLATION FACTORS--METHOD AND DIAGNOSIS 
      P2=1.-P1    
      P2=P1*P2*P2 
      P1=P1+(METH(METHOD)+DIAGEF(DIAG)+DMINT*DIAGEF(DIAG)*METH(METHOD))
     1  *P2       
C                 
C     SELECT SAMPLING FUNCTION
C                 
   4  GO TO(600,500,700),SAMP
C                 
C     NO. OF SUBJECTS DIAGNOSED SCHIZ IS BINOMIAL RV 
 600  LEAVE=.TRUE.
      NT=N        
      CALL BINOM(N,P1,NSCHIZ) 
C                 
C     PRINT RESULTS 
      WRITE(SINK,1001)NSCHIZ 
      GO TO 550   
C                 
C     PASCAL SAMPLING IS CHOSEN
  700 LEAVE=.TRUE.
C                 
C     IR= # DESIRED SCHIZ; Q=PR(NONSCHIZ)
C                 
      IR=N        
      Q=1.-P1     
      CALL PASCAL(IR,Q,IX,NT)
C                 
C     PRINT IT ALL
C                 
      WRITE(SINK,3020) IR,NT
 3020 FORMAT(/' ',I8,' SCHIZOPHRENICS WERE FOUND UNDER EXPERIMENTER.' 
     1 ' ',I8,' SUBJECT-PAIRS (N) WERE SAMPLED TO DO SO.')
C                 
C     PRINT COSTS 
C                 
  550 TSCST=SCOST*NT
      OHCST=(ECST+TSCST)*TOHRT
      TCST=ECST+TSCST+OHCST
      WRITE(SINK,3001)
 3001 FORMAT(//' COSTS:'/)
      PNTIT=LEAVE .OR. MAXPT
      IF(PNTIT) WRITE(SINK,3002) ECST
 3002 FORMAT(' EXPERIMENTAL CONDITION SET-UP COST:',F10.2)
      IF(PNTIT) WRITE(SINK,3003) NT,TSCST
 3003 FORMAT(' SAMPLED N= ',I5,' YIELDS SUBJECT COSTS OF: ',F10.2)
      IF(PNTIT) WRITE(SINK,3004) OHCST
 3004 FORMAT(' OVERHEAD COSTS ARE: ',F10.2)
      WRITE(SINK,3005) TCST
 3005 FORMAT(' THE TOTAL COST IS: ',F10.2)
      IF(LEAVE) RETURN
      GO TO 502   
C                 
 1001 FORMAT(1H0,5X,I4,' SUBJECTS DIAGNOSED SCHIZOPHRENIC.') 
C                 
C                 
C     NOW FOR THE BERNOULLI-ONE-AT-A-TIME SAMPLING
C                 
  500 NT=0        
      NFS=0       
      P9=1.-P1    
      IF(P9.LT.0) P9=0
      IF(P9.GT.1) P9=1
      MAXPT=EXTER.EQ.1 .OR. PCOST .EQ.1
      IF(EXTER.EQ.1) WRITE(SINK,3011)
 3011 FORMAT(//' WELCOME TO THE SCHIZGAME SEQUENCE WHICH ALLOWS'/
     1' YOU, THE EXPERIMENTER, TO DECIDE TO CONTINUE'/
     2' EXPERIMENTATION OR STOP.  AT THE END OF EACH'/
     3' SUBJECT-PAIR SAMPLED, YOU WILL BE ASKED IF YOU'/
     4' WISH TO CONTINUE THE EXPERIMENT OR NOT.  IF YOU'/
     5' DECIDE TO CONTINUE, ENTER A "1"; IF NOT, A "0".'//)
      IF(PCOST.NE.3 .OR. EXTER.EQ.1) WRITE(SINK,3017) ECST,SCOST,TOHRT
 3017 FORMAT(//' BASE COSTS:'/' EXPERIMENTAL CONDITION SET-UP: ',F10.2/
     1 ' PER SUBJECT-PAIR SAMPLED COST: ',F10.2/' OVERHEAD RATE: ',F5.3)
      GO TO 501   
  502 WRITE(SINK,3012)
 3012 FORMAT(//' CONTINUE?--"1"=YES, OR "0"=NO.')
      READ(IDEV1,3013)ERSP
 3013 FORMAT(A1)  
      LEAVE=ERSP.NE.YES1 .AND. ERSP.NE.YESY
      IF(LEAVE) GO TO 550
  501 NT=NT+1     
      SKIZ=.FALSE.
      CALL CHANCE(P9,&505)
      SKIZ=.TRUE. 
  505 IF(SKIZ) GO TO 506
      WRITE(SINK,3014) NT
 3014 FORMAT(/' SAMPLED SUBJECT-PAIR NUMBER: ',I5/
     1 ' THE ACTUAL SUBJECT WAS NOT FOUND TO BE SCHIZOPHRENIC.')
  507 WRITE(SINK,3015) NFS,NT
 3015 FORMAT(/' NUMBER OF SCHIZOPHRENICS FOUND SO FAR: ',I5/
     1 ' TOTAL NUMBER OF SUBJECT-PAIRS SAMPLED SO FAR: ',I5)
      GO TO 550   
  506 WRITE(SINK,3016) NT
      NFS=NFS+1   
 3016 FORMAT(/' SAMPLED SUBJECT-PAIR NUMBER: ',I5/
     1 ' THE ACTUAL SUBJECT IS A SCHIZOPHRENIC.')
      WRITE(SINK,3015) NFS,NT
      GO TO 550   
C                 
      ENTRY MINIT 
      OPEN(UNIT=IDEV3,ACCESS='SEQIN',FILE=MODNM2,DIRECTORY='4030,11')
C  
      READ(IDEV3,2000) EXCOST,ECOSTF,SCOSTB,SCOSTI,MCOSTF,ACOSTF,OHDRT,
     1 SREIM,TCOST,TOCST,
     1 METH,BASE,DMINT,SEXEFF,RATES,DIAGEF,ADOEFF
 2000 FORMAT(9F8.2/2F8.2/9F8.2/2(2F8.2/),3F8.2/3F8.2/F8.2/2F8.2/2F8.2/,
     1 2F8.2/2(F8.2/),2F8.2/9F8.2/2F8.2/3F8.2)
      RETURN      
      END         
C     SUBROUTINE FOR THE NEGATIVE BINOMIAL OR PASCAL DISTRIBUTION
C     FOR SAMPLING THE N BERNOULLI TRIALS NEED TO OBTAIN
C     R (IR) SUCCESSES SPECIFIED BY THE USER.  THE SUBROUTINE
C     RETURNS BOTH THE VALUE FOR N (IN IN), AND THE VALUE FOR
C     X (IN IX), WHICH IS THE NUMBER OF FAILURES (ZEROS) IN THE
C     N TRIALS.  P IS THE PR OF A FAILURE (ZERO).
C                 
      SUBROUTINE PASCAL(IR,P,IX,IN)
      INTEGER*4 IR,IX,IN
      REAL*4 P,Q,B,GVAL
C                 
C      ERROR CHECKS AND FIX-UP
C                 
      IF(IR.LE.0) STOP 4001
      IF(P .LT. -.00001 .OR. P .GT. 1.00001) STOP 4002
      IF(P.LT. (.0000001)) P=.0000001
      IF(P.GT. (.9999999)) P=.9999999
C                 
C      START COMPUTATION
C                 
      Q=1.-P      
C                 
C      MAKE SURE WE CAN NOT EXPECT OVERFLOW
C                 
      IF(((IR*Q/P)+IR) .GT. 2000000000) STOP 4003
      B=P/Q       
C                 
C      X (OR IX) IS PASCAL MEANS WE CAN USE CONVOLUTION
C      OF THE GAMMA AND POISSON (MIDAS SIMULATION MANUAL)
C                 
      CALL RGAMMA(IR,B,GVAL)
      CALL POISSN(GVAL,IX)
      IN=IX+IR    
      RETURN      
      END         

