C      THIS IS A SIMULATION WHICH SIMULATES 
C      RATS BAR PRESSING FOR A REINFORCEMENT
C      OF ELECTRICAL CURRENT APPLIED DIRECTLY
C      TO THEIR BRAINS VIA CHRONICALLY IMPLANTED
C      ELECTRODES.
C                 
C      THIS SUBROUTINE WAS WRITTEN BY MARK KLEIN
C      AND PROGRAMMED BY MARK KLEIN WITH THE AID
C      OF RICH NUSSLOCH IN JANUARY OF 1974
C                 
C                 
C      THIS IS THE FIRST VERSION SO PLEASE PARDON
C      THE MISTAKES.
      SUBROUTINE MODEL(N,NCONDS,*)
C                 
C                 
C@@@@@@ SECTION 1 HAS ENTIRELY TO DO WITH HOUSE 
C      KEEPING CHORES AND THE LIKE.
C                 
C                 
C      VARIABLE TYPE DECLARATION:
C                 
      INTEGER*4 PHASE,POLES,DIRECT,SHAPE,ELECTR,LCYCLE
      INTEGER*4 AREA,LESION,DL1,DRUG,EXPER,N,NDED,LIVE
      INTEGER*4 OPDED1,OPDED2,NBAD,IPROB1,IPROB2,IPROB3
      INTEGER*4 MICROC,ITEMP,FIRSTT(11),CLOCK,TTIME
      INTEGER*4 EXPER1,EXPER2
      INTEGER*4 NOV,SINK,IDEV1,IDEV2,IDEV3,IV(12)
      REAL*4 MICROA,FREQ,PULSE,TRAIN,TESTH,RTEMP,ANDOSE
      REAL*4 DDOSE,INJECT,FDEPRV,PDEAD,OPER,RATE,T1RATE
      REAL*4 T2RATE,FIRSTE(11),DARK(24),DIV,OVERD(5)
      REAL*4 TEFF,DREFF(4,35),SHOT(4),MORT,LMFB(7,7)
      REAL*4 DMISS(6),DV(6),PRESS(20,2,6),PULL,SEX
      REAL*4 SIGN(5),FV(12)
      LOGICAL*4 FLGS(6),EOVER
C                 
C       COMMON STATEMENTS
C                 
      COMMON /VARS1/IV,FV/IO1/NOV,SINK/IO/IDEV1,IDEV2,IDEV3
C                 
C      DATA STATEMENTS
C                 
!!      DATA DMISS/6*350./,FLGS/6*.FALSE./
!!      DATA EOVER/.FALSE./
C                 
C      EQUIVALENCE STATEMENTS
C                 
      EQUIVALENCE (MICROA,FV(1)),(PHASE,IV(1)),(POLES,IV(2)),
     1(DIRECT,IV(3)),(SHAPE,IV(4)),(FREQ,FV(2)),(PULSE,FV(3)),
     2(TRAIN,FV(4)),(ELECTR,IV(5)),(LCYCLE,IV(6)),(TESTH,FV(5)),
     3(RTEMP,FV(6)),(AREA,IV(7)),(ANDOSE,FV(8)),(LESION,IV(8)),
     4(DL1,IV(9)),(DRUG,IV(10)),(DDOSE,FV(9)),(INJECT,FV(10)),
     5(FDEPRV,FV(11)),(EXPER1,IV(11)),(SEX,FV(7))
C                 
C  ************************MTS SPECIFIC SECTION ********************
C  @@@@@@@@@SECTION TO CONTACT CLOCK FOR TIME OF DAY DEFAULT
C      REMOVE THIS WHOLE SECTION OR MODIFY APPROPRIATELY TO FIT
C      THE RESIDENT SYSTEM. 
C                 
       INTEGER*4 IT
       INTEGER*4 JM
       EQUIVALENCE (IT,JM)
C                 
C  ******************************END MTS SECTION *******************
C                 
C        SET UP FOR DATA FILES
C 
      DOUBLE PRECISION MODNM1,MODNM2,MODNM3
      COMMON /DATANM/ MODNM1,MODNM2,MODNM3
      DATA MODNM1/'RATS.DAT'/,MODNM2/'RATS3.DAT'/,MODNM3/'RATS.BIN'/
      DATA DMISS/6*350./,FLGS/6*.FALSE./
      DATA EOVER/.FALSE./
C                 
	NCONDS=NCONDS! A DEC-10 BUG MAKES THIS NECESSARY
C                 
C      NOW A LITTLE ERROR CHECKING!!!
C                 
 1    IF (SHAPE.EQ.1) GO TO 5
      IF (PULSE/1000.LT.1/FREQ) GO TO 5
      WRITE(SINK,800)
      RETURN      
C@@@@@@@@SECTION 2
C      THIS SECTION MAKES THE STUDENT AWARE 
C      OF SOME OF THE TYPICAL PROBLEMS THAT
C      CAN ARRISE IN THE OPERATION DURING
C      WHICH THE ELECTRODES ARE IMPLANTED.
C                 
C                 
 5    NDED=0      
C      SECTION 2.1
C      THE ANESTHETIC DOSE????
      IF (ANDOSE.LE.50) GO TO 10
      PDEAD=(ANDOSE-50)*0.05
      IF (PDEAD.GT.1.) PDEAD=1.
C                 
C      FIND THE NUMBER OF ANIMALS TO DIE FRON O.D.
C                 
      CALL BINOM(N,PDEAD,NDED)
      NDED=MIN0(NDED,N)
C      LET THE PEOPLE KNOW HOW MANY DIED.
      IF (NDED.GT.0) WRITE(SINK,801)NDED
C                 
C      IS THE DOSE TOO STRONG TO CONTINUE
C                 
      IF (NDED.LT.N/2) GO TO 11
      WRITE(SINK,802)
      RETURN      
C      ANDOSE TOO LOW?
 10    IF(ANDOSE.GT.30)GO TO 11
       WRITE(SINK,816)
       RETURN     
 11   LIVE=N-NDED 
C                 
C                 
      EXPER=EXPER1+EXPER2
      OPER=1-EXPER/1000
      OPER=AMAX1(.015,OPER)
      OPER=AMIN1(.55,OPER)
C      THIS PART OF THE SECTION INCORPORATES
C      EXPERIENCE AND CHANCE TO FIND OUT HOW
C      MANY OF THE LIVING ANIMALS WILL MAKE
C      IT THROUGH THE OPERATION.
      CALL BINOM(LIVE,OPER,NDED)
      LIVE=LIVE-NDED
C      NOW DETERMINE HOW MANY SUFFER FROM THE
C      FIRST MALADY
      CALL URAND2(0,NDED,OPDED1)
      NDED=NDED-OPDED1
C                 
C      INFORM THE FOLKS
C                 
      IF(OPDED1.GT.0)WRITE(SINK,803)OPDED1
      IF(NDED.EQ.0)GO TO 50
C      THE NEXT MALADY!!!!!!
      CALL URAND2(0,NDED,OPDED2)
      NDED=NDED-OPDED2
C      INFORM THE FOLKS!!!!!!
      IF(OPDED2.GT.0)WRITE(SINK,804)OPDED2
C      THE LAST MALADY!!!
      IF (NDED.GT.0) WRITE(SINK,805)NDED
C                 
C                 
C      SECTION 2.3
C      THIS SECTION SIMULATES SOME OF THE VARIOUS
C      PROBLEMS THAT THE EXPERIMENTER MIGHT ENCOUNTER
C      ONCE HE IS ABLE TO BEGIN TRAINING THE RATS
C      AND AFTERWARDS WHEN HE IS VARIFYING THE 
C      LOCATION OF THE ELECTRODES
 50   OPER=OPER/5 
C      HOW MANY RATS WILL BE EFFECTED??????
      CALL BINOM(LIVE,OPER,NBAD)
      IF(NBAD.EQ.0) GO TO 100
      LIVE=LIVE-NBAD
C      HOW MANY RATS HAVE THE FIRST PROBLEM???
      CALL URAND2(0,NBAD,IPROB1)
      NBAD=NBAD-IPROB1
C      LET THE FOLKS KNOW.
      IF(IPROB1.GT.0) WRITE(SINK,806)IPROB1
      IF(NBAD.EQ.0) GO TO 100
C      THE SECOND MALADY
      CALL URAND2(0,NBAD,IPROB2)
      NBAD=NBAD-IPROB2
C      LET THE FOLKS KNOW
      IF(IPROB2.GT.0) WRITE(SINK,807)IPROB2
      IF(NBAD.EQ.0) GO TO 100
C      THE THIRD MALADY
      CALL URAND2(0,NBAD,IPROB3)
      NBAD=NBAD-IPROB3
C      LET THE FOLKS KNOW
      IF (IPROB3.GT.0)WRITE(SINK,808)IPROB3
C      THE LAST ONE
      IF(NBAD.GT.0)WRITE(SINK,809)NBAD
C      LET THE FOLKS KNOW HOW MANY THEY HAVE LEFT TO
C      WORK WITH  
 100  WRITE(SINK,810)LIVE
C                 
C      IN THE EVENT THAT ALL OF THE SUBJECTS WERE
C      WIPED OUT THERE IS NO SENSE IN GOING ANY FURTHER
      IF(LIVE.EQ.0)RETURN
C                 
C                 
C@@@@@@SECTION 3A 
C      THIS SECTION PROCESSES THE FIRST PART OF THE WORK 
C      RELATING TO PHASE,POLES,AND DIRECTION. 
      IF(PHASE.EQ.2) WRITE(SINK,811) 
      IF(POLES.EQ.2) GO TO 200 
      IF((PHASE.EQ.1).AND.(DIRECT.EQ.2)) MICROA=MICROA-150 
      IF(MICROA.LT.0)MICROA=1.
C                 
C                 
C@@@@@@SECTION 4  
C      THIS SECTION CALCULATES THE GROSS ENERGY INPUT USED AS 
C      STIMULUS TO THE RAT EACH TIME IT PRESSES THE BAR. 
C      BASED ON THEIR MICROCOULOMB VALUE, THE BASE RATE WILL BE
C      OBTAINED FROM THE BASE RATE DATA MATRIX.
 200  IF(TRAIN.GE.5.) GO TO 210 
       GO TO 220  
C      STIMULUS IS SO LONG THAT IT BECOMES AVERSIVE 
 210   RATE=.001  
      GO TO 300   
C      MICROCOULOMBS BY RECTANGULAR CURRENT 
 220   MICROC=MICROA*FREQ*TRAIN*PULSE*0.001 
 221  IF(MICROC.LT.21) GO TO 222
      GO TO 210   
C        GET BASE RATE FROM THE MATRIX 
 222    ITEMP=MICROC
       IF(AREA.EQ.7) AREA=6
       IF(ITEMP.LE.0) GO TO 210
       RATE= PRESS(ITEMP,1,AREA) 
C      AN EFFECT DUE JUST TO USING SINE WAVES 
 255   IF(SHAPE.EQ.1) RATE=RATE*1.5 
C                 
C                 
C@@@@@@SECTION 3B 
C      FURTHER CALCULATION OF THE EFFECT OF POLES AND PHASE 
      IF(POLES.EQ.1.AND.PHASE.EQ.2) RATE=RATE*.9 
C                 
C                 
C                 
C@@@@@@SECTION 5  
C      THIS SECTION CALCULATES THE TIME OF DAY EFFECT BASED ON
C      THE SPECIFIED LIGHT CYCLE, THE TEST HOUR, AND THE INTENSITY
C      OF THE STIMULATION.
C                 
C                 
C      INITIALIZE TEMPORARY RATE FOR USE IN FIGURING THE EFFECT OF 
C      THE LIGHT CYCLE 
 300  T1RATE= -109 
      ITEMP=(MICROA/100.)+1. 
C      SECTION 5.1
C      CHECK TO SEE IF RATE COMES UNDER THE EFFECT DUE TO INITIAL 
C      EXPERIENCE IN THE APPARATUS OR AFTERWARDS WHEN THE RATE HAS 
C      SETTLED DOWN TO A CYCLIC SORT OF EVENT 
      IF(FIRSTT(ITEMP).LT.TESTH) GO TO 315 
C      IF TESTING OCURRS DURING THAT INITIAL DECLINE 
      T1RATE=RATE*(1.25-FIRSTE(ITEMP)*TESTH**2/100)
C                 
C      IS THE LIGHT CYCLE SET FOR A CONSTANT LIGHT LEVEL? 
C                 
 315   IF(LCYCLE.EQ.13) GO TO 350 
      IF(LCYCLE.EQ.14) GO TO 350 
C      A TWELVE HOUR LIGHT CYCLE SPECIED? 
      IF (LCYCLE.LT.13) GO TO 325 
C      LIGHT CYCLE DEFAULT 
C  ********************** MTS SPECIFIC SECTION ********************
C      SECTION TO FIND THE TIME OF DAY IN HOURS--USES MTS SPECIFIC
C      SUBROUTINE "TIME" TO FILL IN THE IT ARRAY WITH THE TIME OF
C      DAY.       
C      REMOVE THIS SECTION OR MODIFY APPROPRIATELY FOR THE RESIDENT
C      SYSTEM.    
C                 
       CALL TIME(HRS)
       DECODE(5,320,HRS) IT
 320   FORMAT(I2,3X)
C**************END MTS SPECIFIC SECTION*********************
       CLOCK=JM-1
      IF(CLOCK.GT.6) LCYCLE=-CLOCK+19 
      IF (CLOCK.LE.6) LCYCLE=-CLOCK-5 
C      RESCALE TESTH TO CONFORM TO A STANDARD 24 HOUR SCHEME 
 325   TTIME=TESTH+12-LCYCLE 
      TTIME=TTIME-(TTIME/24)*24 
C      IS THE TEST HOUR DURING THE DARK HALF OF THE CYCLE? 
      IF(TTIME.GE.12) GO TO 375 
      GO TO 360   
C      RESCALE TESTH FOR CONSTANT LIGHT LEVELS 
 350     TTIME= TESTH-(TESTH/24)*24 
      IF(LCYCLE.EQ.14) GO TO 375 
C      RESCALE MICROA TO DETERMINE THE INTERACTION 
 360    ITEMP=MICROA/250.+1. 
      TTIME=(TTIME*15)+180 
C      CALCULATE THE CIRCADIAN EFFECT 
      T2RATE=RATE*(1+SIGN(ITEMP)*SIN(TTIME/57.29576)) 
      GO TO 390   
C      FIND THE EFFECT FOR THE HOUR OF DARKNESS 
  375    T2RATE= RATE*DARK(TTIME+1) 
C      DECIDE WHICH EFFECT TO USE:  CHOOSE THE LARGER OF THE TWO 
  390     RATE=T1RATE 
       IF(T1RATE.LT.T2RATE) RATE=T2RATE 
C                 
C                 
C                 
C@@@@@@SECTION 6  
C      THIS SECTION CALCULATES THE EFFECT DUE TO ELECTRODE MATERIAL 
C      DIFFERENCES 
        IF(ELECTR.EQ.2) GO TO 395 
      IF (TESTH.LE.19) GO TO 395 
C      INTERACTION OF PHASE AND ELECTRODE MATERIAL 
      IF(PHASE.EQ.2) RATE=RATE*.80 
C      INTERACTION OF PHASE, DIRECTION, AND MATERIAL 
      IF(PHASE.EQ.1.AND.DIRECT.EQ.2) RATE=RATE*.80 
C      INTERACTION OF PHASE, POLES, AND MATERIAL 
      IF((POLES.EQ.1).AND.(PHASE.EQ.2)) RATE=RATE*.90 
C                 
C                 
C@@@@@@ SECTION 7 
C      THIS SECTION HANDLES THE EFFECT DUE TO THE AMBIENT ROOM TEMPE
 395  CONTINUE    
 396  IF((RTEMP.GT.22.).AND.(RTEMP.LT.28)) GO TO 400 
C      IS THE TEMPERATURE TOO LOW? 
      IF(RTEMP.GE.2) GO TO 398 
      WRITE(SINK,813) 
      RETURN      
 398  IF(AREA.GT.5) GO TO 400 
C      RESCALE TEMPERATURE JUST IN CASE THERE WASN'T ENOUGH TIME 
C      FOR THE RAT TO COOL OFF OR WARM UP COMPLETELY 
      IF(TESTH.LT.1) RTEMP=RTEMP+(25-RTEMP)/2 
      IF(RTEMP.LT.25) GO TO 399 
C      CALCULATE THE EFFECT DUE TO A HIGH ROOM TEMP 
      RATE=RATE*(1-(RTEMP**2-30*RTEMP+225)/300.0) 
      GO TO 400   
C      CALCULATE THE EFFECT DUE TO A LOW ROOM TEMP 
 399   DIV=833.33 
      IF((AREA.EQ.2).OR.(AREA.EQ.3)) DIV=2500.0 
       RATE=RATE*(1-(RTEMP**2-50*RTEMP+625)/DIV) 
C                 
C                 
C                 
C@@@@@@SECTION 8  
C      THIS SECTION DEALS WITH THE DRUG VARIABLE 
400    IF(DDOSE.EQ.0) GO TO 500 
       IF(DRUG.EQ.6) GO TO 500 
C      ARE THERE ANY OVERDOSES OF THE PARTICULAR DRUG? 
       IF(DDOSE.LE.OVERD(DRUG)) GO TO 410 
C      YES??????? 
       PDEAD=(DDOSE-OVERD(DRUG)/(OVERD(DRUG)/2)) 
       IF(PDEAD.GT.1) PDEAD=1. 
C      DETERMINE THE NUMBER OF SUBJECTS TO DIE BY THAT OVERDOSE 
      CALL BINOM(LIVE,PDEAD,NDED) 
      IF (NDED.EQ.0) GO TO 410 
      LIVE=LIVE-NDED 
      WRITE(SINK,814)NDED 
      IF(LIVE.EQ.0) RETURN
C      SALINE?    
 410   IF(DRUG.EQ.5) GO TO 500 
C      FIGURE TIME DELAY UNTIL TESTING 
       TIM= TESTH-INJECT 
C      THAT TIME DELAY TOO LONG    THE DRUGS NO LONGER EFFECTIVE????
       IF(TIM.GT.15) GO TO 500 
       IF(TIM.GT.2) GO TO 415 
C      EFFECT OF TIME ON THE ACTION OF THE DRUG? 
       TEFF=TIM/2 
       GO TO 416  
 415     TEFF=1-(TIM-2)/13 
C      RESCALE DRUG DOSE SO AS TO USE A MATRIX 
 416  IDOSE=DDOSE*SHOT(DRUG) 
 460   RATE=RATE*(1+TEFF*DREFF(DRUG,IDOSE)) 
C                 
C                 
C                 
C@@@@@@SECTION 9  
C      THIS SECTION PROCESSES THE LESION VARIABLE 
500    IF(LESION.EQ.8) GO TO 600 
C      SET UP THE MORTALITY RATE FOR THE LESION OPERATION 
       MORT=.20   
       IF(DL1.NE.08) MORT=.5 
      CALL BINOM(LIVE,MORT,NDED) 
      IF (NDED.GT.0) WRITE(SINK,815)NDED
      LIVE=LIVE-NDED 
C      DID THEY WAIT TOO LONG TO TEST FOR THE EFFECT OF THE LESION??
       IF(TESTH.LT.168) GO TO 520 
      IF((LESION.EQ.1.AND.AREA.EQ.6).OR.(DL1.EQ.1.AND.AREA.EQ.6))GOTO510
      IF((LESION.EQ.2.AND.AREA.EQ.1).OR.(DL1.EQ.2.AND.AREA.EQ.1))GOTO510
      IF((LESION.EQ.4.AND.AREA.EQ.3).OR.(DL1.EQ.4.AND.AREA.EQ.3))GOTO510
       GO TO 600  
 510   RATE=RATE*.5
       GO TO 600  
520   IF((LESION.EQ.1.AND.AREA.EQ.6).OR.(DL1.EQ.1.AND.AREA.EQ.6))GOTO590
      IF((LESION.EQ.2.AND.AREA.EQ.1).OR.(DL1.EQ.2.AND.AREA.EQ.1))GOTO590
      IF((LESION.EQ.4.AND.AREA.EQ.3).OR.(DL1.EQ.4.AND.AREA.EQ.3))GOTO590
       GO TO(530,540,600,530,600,530,600),AREA 
C      MFB AND RELATED AREA LESIONS 
 530   RATE=RATE*LMFB(LESION,DL1) 
       GO TO 600  
C      SEPTAL AREA LESIONS 
 540    IF((LESION.EQ.4) .OR.(DL1.EQ.4)) RATE=RATE*.5 
       GO TO 600  
C      IF THE TEST AREA IS THE SAME ONE THAT THEY DID THE LESION IN.
 590    RATE=RATE*.1 
C                 
C                 
CCCCCCCSECTION 10 
C      THIS SECTION HANDLES THE FOOD DEPREVATION VARIABLE 
600     IF(FDEPRV.LE.1) GO TO 700 
C      CALCULATE THE EFFECT 
       IF(FDEPRV.GE.36)RATE=RATE*(.25*FDEPRV-4.4) 
       IF((FDEPRV.GE.20.AND.FDEPRV.LT.36)) FDEPRV=20. 
       IF(FDEPRV.LE.20.)RATE=RATE*SQRT(FDEPRV) 
 700    IF(RATE.GT.350) RATE=355
 701   ITEMP=MICROC
       IF(ITEMP.EQ.0) ITEMP=1
      FLGS(1)=.TRUE.
      CALL FPOUT1(FLGS,DMISS,&999)
C      SIMULATE THE DATA FOR THE RATS THAT ACTUALLY MADE IT
C                 
      DO 799 K=1,LIVE
      CALL NRAND(RATE,PRESS(ITEMP,2,AREA),PULL)
      IF(PULL.LT.0.05) PULL=ABS(PULL)
      IF(PULL.GT.350) PULL=350
      EOVER=(PULL.EQ.350).OR.EOVER
      DV(1)=PULL  
C      OUTPUT THAT SCORE!!
	KTMP=K
      CALL FPOUT2(KTMP,1,FLGS,DV,&999)
 799  CONTINUE    
       CALL FPOUT3(1)
       IF(EOVER) WRITE(SINK,817)
C      INCREMENT THE EXPERIENCE FACTOR
C                 
7999  EXPER2=EXPER2+N
       EOVER=.FALSE.
 800  FORMAT(' IMPOSSIBLE PULSE AND FREQUENCY COMBINATION')
 801  FORMAT(' ',I3,'  ANIMALS DIED FROM AN ANESTHETIC OVERDOSE')
 802  FORMAT('  KILLER ANESTHETIC DOSE TERMINATE EXPERIMENT')
 803  FORMAT(' ',I3,'  ANIMALS DIED:YOU SEVERED THE MID-SAGGITAL'
     *   ' ARTERY.')
 804  FORMAT(' ',I3,'  ANIMALS DIED: YOU LET BODY TEMP GET TOO LOW.')
 805  FORMAT(' ',I3,'  ANIMALS JUST DIED: REASON UNKNOWN.')
 806  FORMAT(' ',I3,'  ANIMALS NO DATA:  ELECTRODES KNOCKED OUT.')
 807  FORMAT(' ',I3,'  ANIMALS NO DATA:  MISSED BRAIN STRUCTURE.')
 808  FORMAT(' ',I3,'  ANIMALS NO DATA:  INFECTION AROUND ELECTRODE.')
 809  FORMAT(' ',I3,'  ANIMALS NO DATA:  MISC!!!')
 810  FORMAT(' ',I3,'  ANIMALS PRODUCE DATA.  TADA!')
 811  FORMAT(' WITH BIPHASIC CURRENT CURRENT DIRECTION IGNORED.')
 812  FORMAT(' WITH SINE WAVE CURRENT IS CONTINUOUS:IGNORE PULSE.')
 813  FORMAT(' TEMPERATURE TOO LOW FOR SUBJECTS TO SURVIVE AT ALL.')
 814  FORMAT(' ',I3,'  ANIMALS DIED FROM A DRUG OVERDOSE.')
 815  FORMAT(' ',I3,'  ANIMALS DIED FROM THE LESION PROCEDURE.')
 816  FORMAT(' ANESTHETIC DOSE TOO LOW TO OPERATE.')
 817  FORMAT(/' MISSING DATA DUE TO EQUIPMENT OVERLOAD')
 999   RETURN     
C                 
C                 
C                 
C                 
C                 
 9999  ENTRY MINIT
C    
      OPEN(UNIT=IDEV3,ACCESS='SEQIN',FILE=MODNM2,DIRECTORY='4030,11')
C     
      READ(IDEV3,850)PRESS,SIGN,OVERD,DREFF,DARK,FIRSTT,FIRSTE
     1,LMFB,EXPER2
 850  FORMAT(24(10F5.1/),5F5.4,5F6.0/11(12F5.3/),8F5.3/
     12(12F5.2/),11I3/6F10.7/5F10.7/7(7F5.3/),I5)
       RETURN     
	END               

