PROGRAM TPVAX00
C*****TEST TO COMPARE TO SIM APPLICATION PVAX00
C.....INCOMPLETE PVA PROJECTION, DEMOGRAPHIC STOCHASTICITY ONLY
C     * NO OTHER STRUCTURE OR STOCHASTICITY
C     * DEMOGRAPHIC STOCHASTICITY
C        ** BIRTH-DEATH PROCESSES INDEPENDENT AMONG INDIVIDUALS
C        ** BINOMIAL REALIZATION OF PER INDIVIDUAL DEATH RESPONSE
C        ** POISSON REALIZATION OF PER INDIVIDUAL BIRTH RESPONSE
C     * QUASI-EXTINCTION
C        ** FIRST PASSAGE TIME REPORTED, 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*****
C     VERSION OF MAR 16, 2011
C*****
      REAL*8 TBAR,TVAR,ENBAR,ENVAR,DFN
      DIMENSION ISD(11)
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      D       PER CAPITA MORTALITY RATE PER TIME STEP (CONSTANT)
C      B       PER CAPITA BIRTH RATE PER TIME STEP (CONSTANT)
C.....
      ISEED=39
      NTR=100000
      N1=10
      ITH=50
      NQ=0
      D=0.5
      B=0.55
C.....OPEN OUTPUT FILE
      OPEN (8,FILE='TPVAX00.OUT')
C.....SET OVERFLOW TRAP
      NTRAP=10000000
C.....INITIALLIZE RANDOM NUMBER GENERATOR
      ISD(1)=ISEED
      CALL NRINT(ISD)
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 TPVAX00'
      WRITE (*,*) 'DEMOGRAPHIC STOCHASTICITY ONLY'
      WRITE (8,*) 'OUTPUT OF PROGRAM TPVAX00'
      WRITE (8,*) 'DEMOGRAPHIC STOCHASTICITY ONLY'
      NX=0
      NN=0
      TBAR=0.0
      TVAR=0.0
      ENBAR=0.0
      ENVAR=0.0
C.....LOOP TROUGH TRIALS
      DO 1100 I=1,NTR
C.....INITIALIZE TAJECTORY
      N=N1
C.....LOOP THROUGH TIME STEPS
      DO 1070 JT=1,ITH
C.....PROJECTION DYNAMICS
      CALL BINGG(N,D,ND,ISD)
      BN=FLOAT(N)*B
      CALL POIGG(NB,BN,ISD)
      N=N+NB-ND
C.....TRAP QUASI-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