PROGRAM TPVAX06 C*****TEST TO COMPARE TO SIM APPLICATION PVAX06 C.....NOMINALLY COMPLETE PVA PROJECTION, 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 ** SPECIFIED DIRECTLY AS ADDITIVE COMPONENT OF C TOTAL VARIANCE IN FACTOR OF INCREASE (LAMBDA) C * UNATTRIBUTED ENVIRONMENTAL VARIATION C ** LOGNORMAL DISTRIBUTION OF PER CAPITA BIRTH RATE C * PARAMETER UNCERTAINTY C ** JOINT (POSTERIOR) DISTRIBUTION FOR MEAN AND STANDARD C DEVIATION OF ENVIRONMENTAL VARIATION (LOGNORMAL) C OF BIRTH RATE 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 AFTER C FLOATING POINT APPLICATION OF LAMBDA C ** TRUNCATE POPULATION AT 10E+6 TO AVOID INTEGER OVERFLOW C***** C VERSION OF MAR 16, 2011 C***** CHARACTER*8 SMPNM REAL*8 TBAR,TVAR,ENBAR,ENVAR,DFN DIMENSION ISD(11) DIMENSION IWV(2),FWV(2),SMV(2) C.....JOB SPECIFICTIONS C ISEED SEED FOR RANDOM NUMBER GENERATOR C NTR NUMBER OF TRIALS (MONTE CARLO SAMPLE SIZE) C N1 INITIAL POPULATION C ITH TIME HORIZON C NQ QUASI-EXTINCTION THRESHOLD C NKH HARD CEILING ON POPULATION SIZE C DS STD OF DEMOGRAPHIC STOCHASTICITY IN LAMBDA AT ONE INDIVIDUAL 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..... ISEED=39 NTR=2000000 N1=10 ITH=50 NQ=2 NKH=15 DS=0.2 C.....OPEN OUTPUT FILE OPEN (8,FILE='TPVAX06.OUT') C.....SET OVERFLOW TRAP NTRAP=10000000 C.....INITIALLIZE RANDOM NUMBER GENERATOR ISD(1)=ISEED CALL NRINT(ISD) C.....CONVERT THE DEMOGRAPHIC STOCHASTICITY STD TO VARIANCE DS2=DS*DS C.....SET UP FOR READING THE SMP FILE C LUNC IS DEVICE FOR READING FROM CONSOLE C LUNF IS OUTPUT FILE DEVICE NUMBER C LUNP IS DEVICE NUMBER RESERVED FOR SMP FILE C SMPNM MUST BE DECLARED CHARACTER*8 C IWV IS WORKSPACE AND MUST BE DIMENSIONED AT LEAST 2 C FWV IS WORKSPACE AND MUST BE DIMENSIONED AT LEAST 2 C SMV MUST BE DIMENSIONED AT LEAST 2 C..... SMPNM='PVAX05 ' LUNC=5 LUNF=8 LUNP=9 CALL SMPFI(LUNS,LUNC,LUNF,LUNP,CFNAM,NPT,MPS,NHPS,IWV,NSMP,FWV) C.....FRACTION OF QUASI-EXTINCTIONS C AND DISTRIBUTION OF POPULATION SIZE AT TIME HORIZON C AND DISTRIBUTION OF TIME TO QUASI-EXTINCTION BY TIME HORIZON C..... WRITE (*,*) 'OUTPUT OF PROGRAM TPVAX06' WRITE (*,*) 'NOMINALLY COMPLETE LIST OF TYPES OF STOCHASTICITY' WRITE (*,*) 'APPLY DEMOGRAPHIC STOCHASTICITY TO LAMBDA' WRITE (8,*) 'OUTPUT OF PROGRAM TPVAX06' WRITE (8,*) 'NOMINALLY COMPLETE LIST OF TYPES OF STOCHASTICITY' WRITE (8,*) 'APPLY DEMOGRAPHIC STOCHASTICITY TO LAMBDA' NX=0 NN=0 TBAR=0.0 TVAR=0.0 ENBAR=0.0 ENVAR=0.0 C.....LOOP THROUGH TRIALS DO 1100 I=1,NTR 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,NSMP,FWV,SMV) UXF=SMV(1) SXF=SMV(2) C.....LOOP THROUGH TIME STEPS DO 1070 JT=1,ITH C.....ADD THE DEMOGRAPHIC STOCHASTICITY VARIANCE COMPONENT TO ENVIRONMENTAL FN=FLOAT(N) DV=SXF*SXF+DS2/FN SXF=SQRT(DV) C.....LOG SPACE PARAMETERS FOR THE TOTAL VARIATION IN LAMBDA CALL LNSPC(UXF,SXF,ULF,SLF) C.....SAMPLE LOGNORMAL FOR LAMBDA THIS TIME STEP CALL LNORG(ULF,SLF,FA,ISD) C.....APPLY FACTOR OF INCREASE IN FLOATING POINT TO ADVANCE POPULATION FN=FA*FN C.....ROUND TO INTEGER N=IFIX(FN+0.5) C.....ENFORCE HARD CEILING IF (N.GT.NKH) N=NKH C.....TRAP EXTINCTION IF (N.LE.NQ) THEN NX=NX+1 DFN=DFLOAT(JT) TBAR=TBAR+DFN TVAR=TVAR+DFN*DFN GO TO 1100 ENDIF C.....TRAP OVERFLOW IF (N.GE.NTRAP) THEN WRITE (*,*) 'POPULATION EXCEEDED OVERFLOW TRAP' GO TO 9999 ENDIF 1070 CONTINUE C.....REACHED END OF TRAJECTORY WITHOUT QUASI-EXTINCTION NN=NN+1 DFN=DFLOAT(N) ENBAR=ENBAR+DFN ENVAR=ENVAR+DFN*DFN 1100 CONTINUE C.....FINISHED TRIALS: NORMALIZE RESULTS FX=FLOAT(NX)/FLOAT(NTR) WRITE (*,*) 'PROBABILITY OF QUASI-EXTINCTION:',FX WRITE (8,*) 'PROBABILITY OF QUASI-EXTINCTION:',FX IF (NX.GT.0) THEN DFN=DFLOAT(NX) TBAR=TBAR/DFN TVAR=TVAR/DFN-TBAR*TBAR TVAR=DSQRT(TVAR) WRITE (*,*) * 'MEAN TIME TO QUASI-EXTINCTION FOR QUASI-EXTINCT:',TBAR WRITE (*,*) * 'STANDARD DEVIATION OF TIME TO QUASI-EXTINCTION' WRITE (*,*) * ' FOR NOT QUASI-EXTINCT: ',TVAR WRITE (8,*) * 'MEAN TIME TO QUASI-EXTINCTION FOR QUASI-EXTINCT:',TBAR WRITE (8,*) * 'STANDARD DEVIATION OF TIME TO QUASI-EXTINCTION' WRITE (8,*) * ' FOR NOT QUASI-EXTINCT: ',TVAR ENDIF IF (NN.GT.0) THEN DFN=DFLOAT(NN) ENBAR=ENBAR/DFN ENVAR=ENVAR/DFN-ENBAR*ENBAR ENVAR=DSQRT(ENVAR) WRITE (*,*) * 'MEAN FINAL POPULATION NOT QUASI-EXTINCT: ',ENBAR WRITE (*,*) * 'STANDARD DEVIATION OF FINAL POPULATION NOT QUASI-EXTINCT:' WRITE (*,*) * ' ',ENVAR WRITE (8,*) * 'MEAN FINAL POPULATION NOT QUASI-EXTINCT: ',ENBAR WRITE (8,*) * 'STANDARD DEVIATION OF FINAL POPULATION NOT QUASI-EXTINCT:' WRITE (8,*) * ' ',ENVAR ENDIF C.....PROGRAM COMPLETE 9999 STOP END