C*****FILE PVAX01.F C SIM APPLICATION FOR STATISTICS OF POPULATION TRAJECTORY C.....INCOMPLETE PVA PROJECTION (ENVIRONMENAL STOCHASTICITY ONLY), WITH C POPULATION MODELED AS C * NO SPATIAL STRUCTURE C * NO AGE STRUCTURE C * NO DENSITY DEPENDENCE C * NO EXPLICIT ENVIRONMENTAL COVARIATE C * DEMOGRAPHIC STOCHASTICITY C ** NONE C * UNATTRIBUTED ENVIRONMENTAL VARIATION C ** LOGNORMAL DISTRIBUTION OF PER CAPITA BIRTH RATE C * PARAMETER UNCERTAINTY C ** NONE C * QUASI-EXTINCTION C ** FIRST PASSAGE TIME REPORTED AND TERMINATES TRAJECTORY C * SUB TIME-STEP PROCESSES C ** BIRTHS OCCUR AT BEGINING OF TIME STEP C ** DEATHS OCCUR AT END OF TIME STEP BUT DO NOT APPLY C TO THAT TIME STEP'S BIRTHS C * NUMERICAL C ** POPULATION STATE IS REPRESENTED AS INTEGER COUNT C ** ANY FLOATING POINT CALCULATIONS ROUND DOWN FRACTIONAL INDIVIDUALS C ** TRUNCATE POPULATION AT 10E+6 TO AVOID INTEGER OVERFLOW C.....EACH CALL TO SIMSM SIMULATES ONE COMPLETE TRAJECTORY C.....REPORTS AFTER IT TIME STEPS C PRV(1) PROBABILITY OF (QUASI-)EXTINCTION WITHIN TIME HORIZON C (REPORTED AS MEAN IN MARGINAL SUMMARY REPORT; C IGNORE HISTOGRAM--SUPPRESS WITH 0 BINS) C PRV(2) TIME TO FIRST PASSAGE, FROM ABOVE, OF QUASI-EXTINCTION C THRESHOLD, FOR TRAJECTORIES WHICH ACHIEVE FIRST PASSAGE C (CALCULATION WILL BE SLOW IF PROBABILITY OF C (QUASI-)EXTINCTION IS LOW; SUPPRESS BY SETTING SWITCH C TO 1, OR 0 TO SUPPRESS BOTH THIS AND TERMINAL POPULATION, C AND THEN SUPPRESS HISTOGRAM BY 0 BINS) C PRV(3) TERMINAL POPULATION SIZE FOR TRAJECTORIES WHICH DID NOT C ACHIEVE FIRST PASSAGE) C (CALCULATION WILL BE SLOW IF PROBABILITY OF C (QUASI-)EXTINCTION IS HIGH; SUPPRESS BY SETTING SWITCH C TO 2, OR 0 TO SUPPRESS BOTH THIS AND TIME TO C (QUASI-)EXTINCTION AND THEN SUPPRESS HISTOGRAM BY 0 BINS) C.....RECEIVE FROM JOB FILE C SPV(1) N1 INITIAL POPULATION C SPV(2) IT TIME HORIZON C SPV(3) NQ (QUASI-)EXTINCTION THRESHOLD (MAY BE 0) C SPV(4) D PER CAPITA MORTALITY RATE PER TIME STEP (CONSTANT) C SPV(5) UXB MEAN OF ENVIRONMENTAL VARIATION IN PER CAPITA BIRTH C RATE PER TIME STEP (DISTRIBUTION IS LOGNORMAL) C SPV(6) SXB STANDARD DEVIATION OF ENVIRONMENTAL VARIATION IN PER C CAPITA BIRTH RATE PER TIME STEP (DISTRIBUTION IS C LOGNORMAL) C SPV(7) ISW CUMULATION SWITCH C 0 SUPPRESSES CUMULATION OF TIME TO (QUASI-)EXTINCTION C AND TERMINAL POPULATION SIZE C 1 SUPPRESSES CUMULATION OF TIME TO (QUASI-)EXTINCTION C BUT NOT TERMINAL POPULATION SIZE C 2 SUPPRESSES CUMULATION OF TERMINAL POPULATION SIZE C BUT NOT TIME TO (QUASI-)EXTINCTION C 3 CUMULATES BOTH TIME TO (QUASI-)EXTINCTION C AND TERMINAL POPULATION SIZE C***** C VERSION OF MAR 16, 2011 C***** SUBROUTINE SIMIN(LUNS,LUNF,NDT,MDT,DTV,SPV,CFNAM,FWV,IWV) C*****JOB-SPECIFIC DEDICATED SUBROUTINE THAT INITIALIZES VALUES C FOR PROGRAM SIM C.....GIVEN C LUNS DEVICE NUMBER FOR SCREEN OUTPUT, PRESUMED OPEN C LUNF DEVICE NUMBER FOR OUTPUT FILE, PRESUMED OPEN C NDT NUMBER OF OBSERVATIONS IN DATA C MDT NUMBER OF VARIABLES PER OBSERVATION IN DATA C DTV FLOATING POINT ARRAY OF DATA ELEMENTS, IN SINGLE SUBSCRIPT C STORAGE, WHERE ELEMENT IJ=(I-1)*MDT+J IS THE J-TH VARIABLE C FROM OBSERVATION I. C SPV VECTOR OF FIXED PARAMETER VALUES READ FROM JOB FILE C.....RETURNS C CFNAM CHARACTER*8 STRING IDENTIFYING THIS SOURCE FILE, TO MATCH C AGAINST CODE NAME SPECIFIED IN JOB FILE C DTV POSSIBLY REARRANGED, TRANSFORMED, SUBSET OR AUGMENTED; WILL C BE COMMUNICATED TO SUBROUTINE SIMSM IN THIS FORM C SPV POSSIBLY REARRANGED, TRANSFORMED, SUBSET OR AUGMENTED; WILL C BE COMMUNICATED TO SUBROUTINE SIMSM IN THIS FORM C FWV VECTOR AVAILABLE FOR FLOATING POINT WORKSPACE AND COMMUNICATING C ADDITIONAL FLOATING POINT VALUES TO SUBROUTINE SIMSM C IWV VECTOR AVAILABLE FOR INTEGER WORKSPACE AND COMMUNICATING C INTEGER VALUES TO SUBROUTINE SIMSM C OPTIONALLY MAY WRITE ANY DESIRED MESSAGE TO SCREEN (UNIT LUNS) OR C THE OUTPUT FILE (UNIT LUNF) C***** CHARACTER*8 CFNAM DIMENSION DTV(1),SPV(1),FWV(1),IWV(1) C========BEGIN BLOCK OF JOB-SPECIFIC CODE=============================== C.....SELF-IDENTIFY WITH CHARACTER STRING OF UP TO 8 CHARACTERS C IN SINGLE QUOTES (MANDATORY). C DO NOT CHANGE THE NAME OF THE VARIABLE CFNAM C..... CFNAM='PVAX01 ' C.....CONVERT INITIAL POPULATION TO INTEGER, LOAD IN IWV(1) EN1=SPV(1) N1=IFIX(EN1+0.01) IWV(1)=N1 C.....CONVERT TIME HORIZON TO INTEGER, LOAD IN IWV(2) T=SPV(2) IT=IFIX(T+0.01) IWV(2)=IT C.....CONVERT THRESHOLD TO INTEGER, LOAD IN IWV(3) ENQ=SPV(3) NQ=IFIX(ENQ+0.01) IWV(3)=NQ C.....CONVERT SWITCH TO INTEGER, RECODE, AND LOAD IN IWV(4) SW=SPV(7) ISW=IFIX(SW+0.01)-1 IWV(4)=ISW C.....CONVERT MEAN AND STD OF LOG LOGNORMAL ENVIRONMENTAL VARIATION IN C PER CAPITA BIRTH RATE TO MEAN AND STD IN THE LOG SPACE, IN PLACE C..... UXB=SPV(5) SXB=SPV(6) CALL LNSPC(UXB,SXB,ULB,SLB) SPV(5)=ULB SPV(6)=SLB C========END THIS BLOCK OF JOB SPECIFIC CODE============================ 9999 RETURN END C----- SUBROUTINE SIMSM(NDT,MDT,DTV,SPV,FWV,IWV,SMV,ISD,PRV) C*****JOB SPECIFIC DEDICATED PROGRAM TO SAMPLE PROCESS FOR PROGRAM SIM C.....GIVEN C NDT NUMBER OF OBSERVATIONS IN DATA C MDT NUMBER OF VARIABLES PER OBSERVATION IN DATA C DTV FLOATING POINT DATA AS POSSIBLY REARRANGED, TRANSFORMED, C SUBSET OR AUGMENTED BY JOB-SPECIFIC SUBROUTINE SIMIN C SPV VECTOR OF FIXED PARAMETER VALUES, AS SPECIFIED IN THE C JOB FILE AND POSSIBLY MODIFIED IN SIMIN C FWV VECTOR OF FLOATING WORKSPACE, AND POSSIBLE FLOATING POINT C VALUES COMMUNICATED FROM SUBROUTINE SIMIN C IWV VECTOR OF INTEGER WORKSPACE AND POSSIBLE INTEGER VALUES C COMMUNICATED FROM SUBROUTINE SIMIN C SMV FLOATING POINT ARRAY OF ONE NEW SAMPLE FROM THE SMP FILE C INDEXED BY THE PARAMETER DESIGNATION IN THE FILE HEADER C (I.E., THE PARAMETER DESIGNATION IN THE PROGRAM THAT GENERATED C THE SMP FILE) C ISD SEED VECTOR FOR RANDOM NUMBER GENERATOR (DO NOT MODIFY) C.....RETURNS C PRV VECTOR OF VALUES (UP TO 3 ELEMENTS) GENERATED ON THIS C REALIZATION FOR CUMULATION OF JOINT DISTRIBUTION C***** DIMENSION DTV(1),SPV(1),FWV(1),IWV(1),SMV(1),ISD(1),PRV(1) C========BEGIN BLOCK OF JOB-SPECIFIC CODE=============================== C.....UNLOAD JOB SPECIFICATIONS: C INITIAL POPULATION C..... N1=IWV(1) C.....TIME HORIZON IT=IWV(2) C.....POPULATION THRESHOLD NQ=IWV(3) C.....CONSTANT MORTALITY PROBABILITY D=SPV(4) C.....MEAN AND STD OF ENVIRONMENTAL VARIATION OF BIRTH RATE IN LOG SPACE ULB=SPV(5) SLB=SPV(6) C.....RETRIAL SWITCH ISW=IWV(4) C.....DOWNLOAD OF SPECIFICATIONS COMPLETE: C.....INITIALIZE TRIAL OUTCOMES FLAG NTRJ=0 IX=0 IH=0 C.....START A TRIAL TRAJECTORY C INITIALIZE FIRST PASSAGE FLAG, AT NOT C..... 1000 IFP=0 C.....INITIALIZE TRAJECTORY N=N1 C.....LOOP THROUGH TIME STEPS DO 1140 I=1,IT C.....ENVIRONMENTAL STOCHASTICITY C SAMPLE LOGNORMAL FOR THE PER CAPITA BIRTH RATE THIS TIME STEP C..... CALL LNORG(ULB,SLB,B,ISD) C.....DETERMINISTIC BIRTHS, ROUNDING DOWN FN=FLOAT(N) FNB=FN*B NB=IFIX(FNB) C.....DETERMINISTIC DEATHS, ROUNDING DOWN FND=FN*D ND=IFIX(FND) C.....ADVANCE POPULATION C (THIS CANNOT GO NEGATIVE BECAUSE B CANNOT BE NEGATIVE UNDER C THE LOGNORMAL, AND D CANNOT EXCEED 1) C..... N=N+NB-ND C.....CHECK FOR FIRST PASSAGE IF (N-NQ) 1110,1110,1120 1110 IFP=I GO TO 1150 C.....CHECK FOR APPROACHING OVERFLOW 10,000,000 1120 IF (N-10000000) 1140,1140,1130 1130 N=10000000 GO TO 1150 C.....TIME STEP COMPLETE 1140 CONTINUE C.....TRAJECTORY COMPLETE: C.....CHECK FOR FIRST TRIAL 1150 IF (NTRJ) 9999,1160,1190 1160 NTRJ=1 IF (IFP) 9999,1170,1180 1170 PRV(1)=0.0 GO TO 1190 1180 PRV(1)=1.0 C.....CHECK RETRIAL SWITCH C RECODED TO: -1 SUPPRESSES TIME AND POPULATION SIZE C 0 SUPPRESSES TIME ONLY C 1 SUPPRESSES POPULATION SIZE ONLY C 2 NONE SUPPRESSED C..... 1190 IF (ISW) 1200,1210,1230 C.....TIME AND POPULATION SIZE SUPPRESSED 1200 PRV(2)=-9.9 PRV(3)=-9.9 GO TO 9999 C.....TIME ONLY SUPPRESSED 1210 IF (IFP) 9999,1220,1000 1220 PRV(2)=-9.9 PRV(3)=FLOAT(N) GO TO 9999 C.....POPULATION OR NEITHER SUPPRESSED 1230 IF (ISW-1) 9999,1240,1260 C.....POPULATION SUPPRESSED 1240 IF (IFP) 9999,1000,1250 1250 PRV(2)=FLOAT(IFP) PRV(3)=-9.9 GO TO 9999 C.....NEITHER SUPPRESSED 1260 IF (IX) 9999,1270,1290 1270 IF (IFP) 9999,1290,1280 1280 IX=1 PRV(2)=FLOAT(IFP) IF (IH) 9999,1300,9999 1290 IF (IH) 9999,1300,1000 1300 IF (IFP) 9999,1310,1000 1310 IH=1 PRV(3)=FLOAT(N) IF (IX) 9999,1000,9999 C========END THIS BLOCK OF JOB SPECIFIC CODE============================ 9999 RETURN END C***** C END OF FILE C*****