C*****FILE PVAX05.F C SIM APPLICATION FOR STATISTICS OF POPULATION TRAJECTORY C.....INCOMPLETE PVA PROJECTION (NO DEMOGRAPHIC STOCHASTICITY, C SAMPLING JOINT DISTRIBUTION OF PARAMETER UNCERTAINTY), WITH C POPULATION MODELED AS C * NO SPATIAL STRUCTURE C * NO AGE STRUCTURE C * NO EXPLICIT ENVIRONMENTAL COVARIATE C * DENSITY DEPENDENCE C ** HARD CEILNG C * DEMOGRAPHIC STOCHASTICITY C ** NONE C * UNATTRIBUTED ENVIRONMENTAL VARIATION C ** LOGNORMAL DISTRIBUTION OF FACTOR OF INCREASE (LAMBDA) C * PARAMETER UNCERTAINTY C ** JOINT (POSTERIOR) DISTRIBUTION FOR MEAN AND STANDARD C DEVIATION OF ENVIRONMENTAL VARIATION (LOGNORMAL) C OF FACTOR OF INCREASE C * QUASI-EXTINCTION C ** FIRST PASSAGE TIME REPORTED AND TERMINATES TRAJECTORY C * SUB TIME-STEP PROCESSES C ** NONE; FACTOR OF INCREASE IS NOT DECOMPOSED INTO C BIRTHS AND DEATHS C * NUMERICAL C ** POPULATION STATE IS REPRESENTED AS INTEGER COUNT C ** DYNAMICS ARE ROUNDED (NOT TRUNCATED) BACK TO INTEGER 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) NKH HARD POPULATION CEILING C SPV(5) 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 OTHER SPECIFICATIONS C EXPECTS A FILE NAMED PVAX05.SMP WHICH WILL BE A (POSSIBLY RENAMED) C OUTPUT SMP FILE (FROM MTG RUNNING APPLICATION SPECIFIC FILE C RPVAX05) CONTAINING JOINT POSTERIOR SAMPLE OF C UXF MEAN OF LOGNORMAL DISTRIBUTION OF ENVIRONMENTAL C VARIATION IN LAMBDA C SXF SXB OF LOGNORMAL DISTRIBUTION OF ENVIRONMENTAL C VARIATION IN LAMBDA C WHERE UXF IS STORED AS PARAMETER 1 AND SXF C IS STORED AS PARAMETER 2 C THIS SMP FILE IS NOT NAMED IN THE JOB FILE, AND READING IT C IS NOT MANAGED BY THE SIM SHELL; IT IS MANAGED AND READ BY C CALLS TO SUBROUTINES SMPFI AND SMPFR IN THE CODE BELOW. 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='PVAX05 ' 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 HARD POPULATION CEILNG TO INTEGER, LOAD IN IWV(5) FNKH=SPV(4) NKH=IFIX(FNKH+0.01) IWV(5)=NKH C.....CONVERT SWITCH TO INTEGER, RECODE, AND LOAD IN IWV(4) SW=SPV(5) ISW=IFIX(SW+0.01)-1 IWV(4)=ISW C.....SET UP FOR READING THE SMP FILE LUNC=5 LUNP=8 CALL SMPFI(LUNS,LUNC,LUNF,LUNP,CFNAM,NPT,MPS,NHPS,IWV(11),NSMP, * FWV) IWV(6)=NPT IWV(7)=MPS IWV(8)=NHPS IWV(9)=NSMP 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.....HARD POPULATION CEILING NKH=IWV(5) C.....RETRIAL SWITCH ISW=IWV(4) C.....DOWNLOAD OF JOB SPECIFICATIONS COMPLETE: C UNLOAD SMP FILE CONTROL SPECIFICATIONS C..... LUNP=8 NPT=IWV(6) MPS=IWV(7) NHPS=IWV(8) NSMP=IWV(9) 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.....PARAMETER UNCERTAINTY C SAMPLE FROM THE SMP FILE JOINT VALUES OF THE MEAN AND THE STD OF C THE (LOGNORMAL) DISTRIBUTION OF ENVIRONMENTAL VARIATION IN LAMBA, C FOR THIS TRAJECTORY C..... CALL SMPFR(LUNP,NPT,MPS,NHPS,IWV(11),NSMP,FWV,SMV) IWV(9)=NSMP UXF=SMV(1) SXF=SMV(2) C.....LOG SPACE PARAMETERS FOR THE ENVIRONMENTAL VARIATION IN LAMBDA CALL LNSPC(UXF,SXF,ULF,SLF) C.....LOOP THROUGH TIME STEPS DO 1140 I=1,IT C.....ENVIRONMENTAL STOCHASTICITY C SAMPLE LOGNORMAL FOR LAMBDA THIS TIME STEP C..... CALL LNORG(ULF,SLF,FA,ISD) C.....APPLY FACTOR OF INCREASE IN FLOATING POINT TO ADVANCE POPULATION FN=FLOAT(N) FN=FA*FN C.....ROUND TO INTEGER N=IFIX(FN+0.5) C.....ENFORCE HARD CEILING IF (N-NKH) 1100,1100,1010 1010 N=NKH C.....CHECK FOR FIRST PASSAGE 1100 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*****