C	THIS IS THE SUBROUTINE FOR DOPE
C	WRITTEN BY HOWARD EICHENBAUM AND TRUDY VILLARS
C	AND PROGRAMMED BY RICH NUSSLOCH ONE DREARY
C	DAY IN N0VEMBER 1972
C
C
C      VERSION 2--PROGRAMMED IN SUMMER 1973
C
C	THIS INCLUDES REVISED DOSE RESOPONSE CURVES
C	AND REVISIONS TO MAKE THE SIMULATION PERFORM
C	IN A MANNER MORE AKIN TO THE EVERYDAY FIASCOS(??)
C	OF RESEARCH
C
C      VERSION 2 AUTHORS:  HOWARD EICHENBAUM, TRUDY VILLARS, AND
C      RICHARD NUSSLOCH--PROGRAMMED BY R. NUSSLOCH
C
C
      SUBROUTINE MODEL(N,NCONDS,*)
C
C>>>>>>>>>> SECTION 1
C
C	VARIABLE SET UP--SEE ALSO SECTION 13
C
C@@@@@@@@@@  SECTION 1.1
C	VARIABLE TYPE DECLARATION
C
      REAL*4 MK,I,STRDOS(3,3),DV(6),K,FV(12),DMISS(6),S(15,20)
      REAL*4 DOSE,TR,R,A,B,SMALD,MK1,R1,TR1
      REAL*4 Y,TOLV,LDH(4,3,2),VARI,FA1,REFAC,F1,F2,F3,F4
      REAL*4 A1,B1,SCKFAC,T(6,10),Y3,PDEAD,DH,DL
      REAL*4 Y1,B2,RETY,RETY9,Y9,COST,X2,X1,ERR,PRTI,STDE,STDE2,RND1,RND2
      INTEGER*4 TASK,DRUG,TEST,STRAIN,TOL,RET,SINK,INACQ1(3,3,4)
      INTEGER*4 IV(12),DRUG1,NDEAD,XDRUG,INRET1(4,2,4)
      INTEGER*4 IIDX,ND,DRM(9),I1,I2,I3,M,M1
      LOGICAL*4 FLGS(6),DEAD1,ODRET
C
C@@@@@@@@@@  SECTION 1.2
C      INITIALIZATION FOR TIME-OF-DAY EFFECT
C
      INTEGER*4 IHRS
      REAL*4 TDAY
      EQUIVALENCE (TDAY,FV(7))
C
C      END OF SECTION 1.2 
C
C@@@@@@@@@@  SECTION 1.3
C	DATA, COMMON, AND EQUALVALENCE DECLARATIONS
C
!!      DATA DMISS /6*0./
!!      DATA FLGS /.TRUE.,5*.FALSE./
      COMMON /VARS1/IV,FV
      COMMON /IO1/NOV,SINK
      COMMON /IO/IDEV1,IDEV2,IDEV3
      EQUIVALENCE (TASK,IV(1)),(DRUG1,IV(2)),(TEST,IV(3)),
     1(STRAIN,IV(4)),(MK,FV(1)),(TOL,IV(5)),(I,FV(2)),
     2(TR,FV(3)),(R,FV(4)),(RET,IV(6)),(VARI,FV(6)),
     3(XDRUG,IV(7))
C  
      DOUBLE PRECISION MODNM1,MODNM2,MODNM3
      COMMON /DATANM/ MODNM1,MODNM2,MODNM3
      DATA DMISS /6*0./
      DATA FLGS /.TRUE.,5*.FALSE./
      DATA MODNM1/'DOPE.DAT'/, MODNM2/'DOPE3.DAT'/,MODNM3/'DOPE.BIN'/
C
C---------------------
C  
C  
	RND=RND!DEC-10 COMPILER BUG REQUIRES THIS
	NCONDS=NCONDS!DITTO
C@@@@@@@@@@  SECTION 1.4
C	IF TRIALS (TR) = 0 OR
C	IF REINFORCEMENT < 1, SIMULATE ALL SUBJECTS WITH MISSING DATA
C
      IF(R.LT.1 .OR. TR.EQ.0) GO TO 887
C
C@@@@@@@@@@  SECTION 1.5
C      INITIALIZE VALUES FOR PROGRAM FOR EACH GROUP
C
      DRUG=DRUG1
      TOLV=TOL
      MK1=MK
      TR1=TR
      FLGS(1)=.TRUE.
      FLGS(2)=TEST.GE.2
      IIDX=1
      IF(I.GE.0) IIDX=2
C
C>>>>>>>>>> SECTION 2
C	XDRUG ROUTINE
C
C@@@@@@@@@@  SECTION 2.1
C	PASS IF XDRUG=0
C
      IF(XDRUG.EQ.0) GO TO 727
C
C@@@@@@@@@@  SECTION 2.2
C	PARSE THE VALUE OF XDRUG FOR THE THIRD DIGIT (I2)
C
      I1=XDRUG/10
      F1=I1
      F2=F1/10.
      I3=F2
      F3=I3
      F4=((F2-F3)*10.)+.1
      I2=F4
      IF(I2.LT.1) I2=1
      IF(I2.GT.9) I2=9
C
C@@@@@@@@@@  SECTION 2.3
C      SELECT DRUG VALUE AND DOSAGE ON BASIS OF I2
C
      R1=.1
      GO TO(991,992,991,992,991,992,991,992,993),I2
  992 R1=10.
      GO TO 991
  993 R1=0.
  991 MK1=R1*MK
      DRUG=DRM(I2)
C
C>>>>>>>>>> SECTION 3
C	ROUTINE TO DETERMINE IF THERE IS AN O.D.
C
C@@@@@@@@@@  SECTION 3.1
C	FIRST SELECT END-POINTS OF LETHAL DOSE RANGE AND COMPUTE
C	O.D.-IN-RETENTION-ONLY INDICATOR (ODRET)
C
  727 DH=LDH(DRUG,STRAIN,2)
      DL=LDH(DRUG,STRAIN,1)
      ODRET=(RET.EQ.2).AND.(TEST.EQ.2)
C
C@@@@@@@@@@  SECTION 3.2
C	COMPUTE PROBABILITY OF SUBJECT DYING (PEAD)
C	NOTE THAT MINIMUM(PDEAD) = .02
C
      PDEAD=.02
      IF(MK1.GT.DL .AND. MK1.LT.DH) PDEAD=(MK1-DL)/(DH-DL)
      IF(MK1.GE.DH) PDEAD=1.
C
C@@@@@@@@@@  SECTION 3.3
C	IF ALL ANIMALS O.D., SIMULATE SUBJECTS WITHOUT FURTHER ADO AND
C	COMPUTATION; OTHERWISE, PROCEDE
C
      NDEAD=N
      IF(PDEAD.EQ.1.) GO TO 887
      NDEAD=0
C
C>>>>>>>>>> SECTION 4
C	RESCALE DOSAGE
C
C@@@@@@@@@@  SECTION 4.1
C	SELECT DRUG
C
  715 GO TO(750,751,752,753),DRUG
C
C@@@@@@@@@@  SECTION 4.2
C	UPS
C
  750 IF(MK1.LE.3) DOSE=.667*MK1
      IF(MK1.GT.3 .AND. MK1.LE.5) DOSE=MK1-1
      IF(MK1.GT.5) DOSE=4.+.5*(MK1-5.)
      GO TO 754
C
C@@@@@@@@@@  SECTION 4.3
C	DOWNS
C
  751 IF(MK1.LT.12) DOSE=MK1/12.
      IF(MK1.GE.12 .AND. MK1.LT.25) DOSE=1.+.077*(MK1-12.)
      IF(MK1.GE.25 .AND. MK1.LT.40) DOSE=2.+.067*(MK1-25.)
      IF(MK1.GE.40 .AND.MK1.LT.50) DOSE=3.+.1*(MK1-40.)
      IF(MK1.GE.50) DOSE=4.+.2*(MK1-50)
      GO TO 754
C
C@@@@@@@@@@  SECTION 4.4
C	NARC
C
  752 IF(MK1.LT.60) DOSE=MK1/60.
      IF(MK1.GE.60 .AND. MK1.LT.125) DOSE=1.+.015*(MK1-60)
      IF(MK1.GE.125 .AND. MK1.LT.200) DOSE=2.+.013*(MK1-125)
      IF(MK1.GE.200 .AND. MK1.LT.250) DOSE=3.+.02*(MK1-200)
      IF(MK1.GE.250) DOSE=4.+.2*(MK1-250)
      GO TO 754
C
C@@@@@@@@@@  SECTION 4.5
C	ABIO
C
  753 IF(MK1.LT.400) DOSE=MK1/400.
      IF(MK1.GE.400 .AND. MK1.LT.800) DOSE=1.+.005*(MK1-400)
      IF(MK1.GE.800) DOSE=3.+.01*(MK1-800.)
C
C>>>>>>>>>> SECTION 5
C	SET UP RANDOM ERROR DISTRIBUTION
C
C@@@@@@@@@@  SECTION 5.1
C	PARAMETERS OF ERROR DISTRIBUTION DEPEND ON STRAIN AND TOL
C
  754 X2=0.
      IF(STRAIN.EQ.1) X1=.1
      IF(STRAIN.EQ.2) X1=.075
      IF(STRAIN.EQ.3) X1=.05
      IF(TOL.EQ.2) X2=.25
      ERR=6.*(X1+X2)
C
C>>>>>>>>>> SECTION 6
C	SECTION TO FIND TIME OF DAY IN HOURS--
C
	CALL TIME (HRS)
	DECODE(5,483,HRS) IHRS
483	FORMAT(I2,3X)
C
C
C>>>>>>>>>> SECTION 7
C	SET UP THEORETICAL MEAN POINTS FOR ACQUISITION
C
C@@@@@@@@@@  SECTION 7.1
C	SELECT MAIN EFFECT OF STRAIN BY SELECTING 3 PARAMETERS OF PARABOLIC
C	FUNCTION (4TH IS REINFORCEMENT) AND COMPUTE APPROPRIATE ZERO DOSE
C	THEORETICAL GROUP MEAN (Y1)
C
      K=STRDOS(1,STRAIN)
      A=STRDOS(2,STRAIN)
      B=STRDOS(3,STRAIN)
      Y1=K*(R-A)*(R-A)+B
      Y=Y1
C
C@@@@@@@@@@  SECTION 7.2
C	SELECT STRAIN X TASK X DRUG ACQUISITION-A1-B1-FUNCTION.
C
      M=INACQ1(STRAIN,TASK,DRUG)
C
C@@@@@@@@@@  SECTION 7.3
C	COMPUTE DOSE EFFECTIVENESS (SMALD) ON THE BASIS OF RESCALED DOSAGE, AND
C	TIME OF INJECTION ( I IN THE PROGRAM), TEST, STATE (RET), AND DRUG
C
      SMALD=0.
      IF(((TEST.EQ.1).OR.(TEST.EQ.2.AND.RET.EQ.2)).AND.I.GT.0) GO TO 757
      IF(I.LT.(-60.) .OR. I.GT.60) GO TO 757
      IF(DRUG.EQ.4) TOLV=1.
      IF(I.LT.0) PRTI=(ABS(ABS(I)-15.))/60.
      IF(I.GE.0) PRTI=I/60.
      IF(DRUG.EQ.4 .AND. I.GT.-15 .AND. I.LT.0) PRTI=0.
      SMALD=DOSE*(1.-PRTI)/TOLV
C
C@@@@@@@@@@  SECTION 7.4
C	IF DOWNS, PASSIVE AVOIDANCE IS SELECTED (M=21) GO TO
C	NON-PARABOLIC ROUTINE
C
  757 IF(M.EQ.21) GO TO 102
C
C@@@@@@@@@@  SECTION 7.5
C	FILTER OUT POTENTIAL EXPONENTIAL OVERFLOW BOMB
C
C
      IF(DRUG.EQ.3 .AND. SMALD.GE.100) SMALD=100.
      IF(DRUG.NE.3 .AND. SMALD.GE.35) SMALD=35
C
C@@@@@@@@@@  SECTION 7.6
C	RECOMPUTE A1 & B1 PARAMETERS OF PARABOLIC CURVE ON BASIS OF SMALD,
C	DRUG, TASK, AND STRAIN MAIN AND/OR INTERACTION EFFECTS
C	AND COMPUTE AND FILTER THEORETICAL GROUP MEAN WITH EFFECTIVE DRUG DOSE
C	Y IS THE THEORETICAL GROUP MEAN
C
      A1=A-((S(1,M)*SMALD)+S(2,M)*2**(SMALD/S(3,M))+S(4,M))
     1+S(15,M)*(A-R)
      B1=B/(S(5,M)*(SMALD/S(6,M))**S(7,M)+1.)+
     1(S(8,M)*B**(SMALD/S(9,M))+S(10,M))+(SMALD*S(14,M))
      B2=(S(11,M)-(S(12,M)*(S(13,M)-SMALD)*(S(13,M)-SMALD)))
      IF(B2.LT.0) B2=0.
      B1=B1-B2
      Y=K*(R-A1)*(R-A1)+B1
      IF(DRUG.LE.2.AND.TASK.NE.1.AND.A1.GT.0.AND.R.GT.A1) Y=B1
      GO TO 103
C
C@@@@@@@@@@  SECTION 7.7
C	ROUTINE FOR DOWNS-PASSIVE AVOIDANCE--NOT A PARABOLIC CURVE
C	NOTE THAT TRAINING EFFECT IS DISABLED AS PART OF THE EFFECT
C
  102 IF(TASK.NE.1 .AND. R.GT.A) Y=B
      IF(Y.GE.70.) GO TO 485
      Y=70.
      IF(SMALD.GE.5) GO TO 485
      Y=Y+.2*SMALD*(B-Y)
  485 IF(TR.LE.70 .AND. SMALD.LT.5) TR1=.2*SMALD*(70.-TR)+TR
      IF(TR.GT.70 .OR. SMALD.GE.5) TR1=TR
C
C>>>>>>>>>> SECTION 8
C	SET UP THEORETICAL POINTS FOR RETENTION
C
C@@@@@@@@@@  SECTION 8.1
C	FIRST SELECT EFFECTS OF DRUG, INTIME, AND STATE (RET) BY SELECTING
C	PARAMETER OPTION (M OR M1)
C
  103 IF(TEST.EQ.1) GO TO 761
      M1=INRET1(DRUG,IIDX,RET)
      M=M1
C
C@@@@@@@@@@  SECTION 8.2
C      NARC FOR DN AND DD IS A SPECIAL OF CASE M=1 AS RETY MUST BE >=Y1
C
      IF(M.EQ.11) M=1
C
C@@@@@@@@@@  SECTION 8.3
C	COMPUTE FA1 FACTOR FOR ABIO'S EFFECT ON RETENTION ON BASIS OF
C	SMALD AND TRAINING
C
      IF(TR.LE.70.) Y3=Y1*(TR/70.)*(TR/70.)
      IF(TR.GT.70.) Y3=Y1+(TR-70.)
      IF(Y3.LT.0.) Y3=0.
      IF(Y3.GT.TR) Y3=TR
      FA1=1.5-((1.-Y3/TR)*SMALD)
      IF(FA1.LT.1) FA1=1.
C
C@@@@@@@@@@  SECTION 8.4
C	COMPUTE THEORETICAL RETENTION GROUP MEAN (RETY) AND SELECT
C	APPROPRIATE ACQUISITION GROUP MEAN (Y OR Y1--SEE SECTION 7)
C	ACCORDING TO STATE (CALLED RET IN PROGRAM)
C
      RETY=(T(1,M)+T(2,M)*SMALD)*Y1+(T(3,M)*FA1+T(4,M))*Y
      Y=T(5,M)*Y+T(6,M)*Y1
      IF(M1.GE.7 .AND. RETY.LT.Y1) RETY=Y1
C
C>>>>>>>>>> SECTION 9
C	COMPUTE AND ADD:
C
C@@@@@@@@@@  SECTION 9.1
C      TRAINING EFFECT
C
  761 IF(TR.LE.70.) Y=(TR1/70.)*(TR1/70.)*Y
      IF(TR.GT.70.) Y=Y+(TR-70.)
      IF(TR.LE.70. .AND. TEST.EQ.2) RETY=(TR1/70.)*(TR1/70.)*RETY
      IF(TR.GT.70. .AND. TEST.EQ.2) RETY=RETY+(TR-70.)
C
C@@@@@@@@@@  SECTION 9.2
C      EXCESSIVE REINFORCEMENT AND ABIO SICKNESS EFFECTS
C
      REFAC=0.
      SCKFAC=0.
      IF((R.GT.5) .AND. (TASK.GE.2 .OR. STRAIN.EQ.3)) REFAC=(R-5.)*15.
      IF(DRUG.EQ.4 .AND. MK1.GT.500) SCKFAC=Y*.002*(MK1-500)
      RETY=RETY-REFAC
      Y=Y-REFAC
      IF(RET.GE.3 .OR. TEST.EQ.1) Y=Y-SCKFAC
      IF(RET.EQ.2 .OR. RET.EQ.4) RETY=RETY-SCKFAC
C
C>>>>>>>>>> SECTION 10
C      SIMULATE SUBJECTS
C
C@@@@@@@@@@  SECTION 10.1
C	FIRST SET UP OUTPUT SEQUENCE AND BEGIN LOOP THROUGH N SUBJECTS
C
  887 CALL FPOUT1(FLGS,DMISS,&999)
      DO 555 M=1,N
C
C@@@@@@@@@@  SECTION 10.2
C	STORE THEORETICAL GROUP MEANS SO THEY WONT BE LOST FOR NEXT SUBJECT
C
      Y9=Y
      RETY9=RETY
C
C@@@@@@@@@@  SECTION 10.3
C	FUDGE BACK UNDER TOTAL TRIALS SPECIFIED IF THEORETICAL MEAN >= TRIALS
C	INPUT
C
      F1=.91*TR
      F2=.06*F1
      IF(Y9.GT.TR) CALL NRAND(F1,F2,Y9)
      IF(RETY9.GT.TR) CALL NRAND(F1,F2,RETY9)
C
C@@@@@@@@@@  SECTION 10.4
C	SET OUTPUT DV'S TO MISSING DATA VALUE (0.) AND DELIVER STUDENT
C	MISSING DATA IF REINFOECEMENT=0 OR PDEAD=1.
C
      DV(1)=0.
      DV(2)=0.
      IF(R.EQ.0 .OR. PDEAD.EQ.1) GO TO 716
C
C@@@@@@@@@@  SECTION 10.5
C	DETERMINE IF THIS SUBJECT DIES
C
      CALL BINOM(1,PDEAD,ND)
      DEAD1=ND.GE.1
      IF(DEAD1) NDEAD=NDEAD+1
      IF(DEAD1.AND.(.NOT.ODRET)) GO TO 716
C
C@@@@@@@@@@  SECTION 10.6
C	ADD IN RANDOM ERROR (SECTION 5) FACTOR.  CAN BE LEFT OUT BY
C	BY SPECIFYING VARI (VVVV)=.283
C
  307 IF(VARI.GE.(.282) .AND. VARI.LE.(.284)) GO TO 711
      STDE=ERR*Y9
      STDE2=STDE/3.
      CALL NRAND(0.,STDE,RND1)
      CALL NRAND(0.,STDE2,RND2)
      IF(TEST.EQ.2) RETY9=RND1*.75+RND2+RETY9
      Y9=Y9+RND1
C
C@@@@@@@@@@  SECTION 10.7
C      ADD IN TIME-OF-DAY EFFECT
C
  711 IF(TDAY.GE.(.668).AND.TDAY.LE.(.670)) GO TO 851
      Y9=Y9*(1-(ABS(12.-IHRS)/24.))
      RETY9=RETY9*(1-(ABS(12.-IHRS)/24.))
C
C
C@@@@@@@@@@  SECTION 10.8
C	TRUNCATE AND FILTER ACQ, IF APPROPRIATE DO LIKEWISE WITH RET SCORE
C	SCORE
C
  851 I1=Y9
      Y9=I1
      IF(Y9.GT.TR) Y9=TR-3.
      IF(Y9.LT.0.) Y9=0
      DV(1)=Y9
      IF((.NOT.FLGS(2)) .OR. DEAD1) GO TO 716
      I1=RETY9
      RETY9=I1
      IF(RETY9.GT.TR) RETY9=TR-3.
      IF(RETY9.LT.0) RETY9=0.
      DV(2)=RETY9
C
C@@@@@@@@@@SECTION 10.9
C      PUT OUT DV VALUES AND END LOOP TO SIMULATE SUBJECTS
C
716	MTMP=M
      CALL FPOUT2(MTMP,1,FLGS,DV,&999)
  555 CONTINUE
C
C>>>>>>>>>> SECTION 11
C	COST; DEAD ANIMALS OUTPUT
C
C@@@@@@@@@@  SECTION 11.1
C	COMPUTE COST BASED ON TEST, STRAIN, STATE, DOSE, N
C
      F1=3.
      IF(STRAIN.EQ.1) F1=5.
      F1=((.5)*TEST)*(.02)*TR+F1
      F2=.1*DOSE
      IF(RET.EQ.4) F2=2.*F2
      F1=F1+F2
      IF(TOL.EQ.2) F1=F1-2.
      COST=F1*N
      IF(COST.LT.0) COST=0.
C
C@@@@@@@@@@  SECTION 11.2
C      PUT OUT COST AND NUMBER OF DEAD ANIMALS
C
      CALL FPOUT3(COST)
      IF(NDEAD.GT.0) WRITE(SINK,303) NDEAD
  303 FORMAT(//2X,I5,1X,'ANIMALS DIED OR O.D.ED--WHY DO YOU THINK
     1 THEY CALL IT DOPE???')
C
C>>>>>>>>>> SECTION 12
C      RETURN FROM SIMULATING A GROUP
C
  999 RETURN
C
C>>>>>>>>>> SECTION 13
C	ENTRY TO SET UP ARRAYS FOR OPTION, MAPPING AND THE LIKE
C	(SEE SECTION 1).  INPUT FILE IS DOPEDAT2
C
      ENTRY MINIT
C  
      OPEN(UNIT=IDEV3,ACCESS='SEQIN',FILE=MODNM2,DIRECTORY='4030,11')
C  
      READ(IDEV3,304) DRM,INACQ1,STRDOS,S,T,INRET1,LDH
  304 FORMAT(9I2/2(18I2/),9F5.1,20(/15F5.1),5(/12F5.1),
     1     2(/16I2),2(/12F5.1))
      RETURN
      END

