!  DUSTBUNNY.F90 -- Add-in module for PC-ORD version 6
!  This simultaneous population process model was designed  to be as simple as possible 
!      to develop dust bunny distributions by community development in unoccupied space. 
!      Populations are described by N, the number of individuals. The model is discrete 
!      with respect to time. 
!  The main components of the model are stochastic immigration, community-limited growth, 
!      species-specific competitive ability, and variable species performance along an 
!      environmental gradient.


      SUBROUTINE DUSTBUNNY

! FILES
!      FILE CONTENTS      READ      WRITE     FILE NAMES
!  result file              ---     write=21  RESULT.TXT
!  options                read=35    ---      DUSTBUNNYOPT.TXT
!  PC-ORD options file    read=10    ---      PCORD6.OPT
!  *.wk1 spreadsheet, main matrix with species abundances or log abundances
!  *.wk1 spreadsheet, second matrix with sample unit total abundance,
!                     number of generations, and positions on environmental gradients

                                  ! The following compiler directive is specific
                                  ! to Intel Visual Fortran
!DEC$ ATTRIBUTES DLLEXPORT:: DUSTBUNNY

      USE IFWIN                    ! Requests Intel Visual Fortran windows library module
                                   ! Intel Fortran Windows applications must use the IFWIN module 
      USE IFPORT
      INCLUDE 'ThermInt6.f90'        ! Include interfaces to therm4.dll
      PARAMETER (MAXENV=10)

      CHARACTER*1, ALLOCATABLE :: VTYPE(:),VTYPE2(:)
      CHARACTER*1 DUMMY,COMMA,BLANK1
      CHARACTER*8 AROWS,ACOLS,AENV,ATEMP,BLANK8,LABEL
      CHARACTER*8, ALLOCATABLE :: ENT(:), ATT(:), ENV(:)  ! row and column names
      CHARACTER*9 OUTPUTFILE
      CHARACTER*10 YESNO(0:1),SEEDSOURCE(0:1),EWEIGHTING(0:2)
      CHARACTER*16 BLANK16
      CHARACTER*20 VERSION                    ! PC-ORD version number
      CHARACTER*70 PROCESS                    ! Process identifier used in therm4.dll
      REAL*4, ALLOCATABLE :: X(:,:), E(:,:)   ! data matrices
      REAL*4, ALLOCATABLE :: SPMAX(:), SPMAXLOG(:), ROWTOT(:), ROWTOTLOG(:), TEMPVEC(:), WORKVEC(:)
      REAL*4, ALLOCATABLE :: ENVOPT(:,:), COMPCOEF(:), IMMPRESSURE(:), POISSONPARAM(:), RINTRINSIC(:)
      REAL*4 AVEGEN,RANGEGEN,STICK(MAXENV),ETARGET(MAXENV)
      INTEGER*4 IMMIGRANTS
      INTEGER*2 IHR,IMIN,ISEC,IHUN,IMON,IDAY,HANDL
      
      LOGICAL FILEOPENED,FILEEXIST

      COMMON/VERS /VERSION

INCLUDE 'ThermCom.F90'  ! include commons needed for thermometer
INCLUDE 'StartDLL.F90'  ! initialize dll

      DATA SEEDSOURCE/'  Use time','User-suppl'/
      DATA YESNO/'        No','       Yes'/
      DATA EWEIGHTING/'   Uniform','Descen-Lin','Descend-BS'/

      CKESC=0
      BLANK1=' '
      BLANK8='        '
      BLANK16='                '
      LABEL=BLANK8
      COMMA=','
      INUMBER=0
      IMATRIX=0
      ICHECK=0
      ITBLHEADER=1

      NROWSTOMAKE=20
      NCOLSTOMAKE=10
      MONTESOURCE=0
      ISEED1=1
      RMAX=2.
      NGENERATIONS=20
      AVEGEN=20.    ! Average number of generations to run model
      RANGEGEN=20.  ! Range in number of generations (half above, half below the mean)
      CARRYINGCAPACITY=1000.  ! Community-level carrying capacity
      TOLERANCE=30.           ! Niche width, spread of species response function to environment
      DISPLIMITATION=1.
      NGRADIENTS=1
      AROWS='SUs     '
      ACOLS='Species '
      AENV='Environ '
      XMISSING=-9999.
      IWRITELOGX=1
      LMAX=15


      LOOP1BEG=1  ! Initialize optional loop parameters as if none are being used
      LOOP1END=1
      LOOP1INCR=1
      LOOP2BEG=0
      LOOP2END=0
      LOOP2INCR=1
      LOOP3BEG=0
      LOOP3END=0
      LOOP3INCR=1
      LOOP4BEG=0
      LOOP4END=0
      LOOP4INCR=1
      LOOPCOUNTER=0

      OPEN (21, FILE='RESULT.TXT', STATUS='UNKNOWN')  ! Open result file

! Read PC-ORD version number from PCORD6.OPT.  You can also use
! this file to read the currently selected options for other
! PC-ORD procedures.  If you wish to use the title line in your
! output, read the second line in this file.
      OPEN (10,FILE='PCORD6.OPT',STATUS='OLD',IOSTAT=IERR)
      IF (IERR.NE.0) THEN
         IERR=11
         GOTO 999
      ENDIF
      READ (10,'(A20)',ERR=888) VERSION
      GOTO 889
888   WRITE (21,'(//A)') ' Error reading PC-ORD options file'
      CLOSE (10,STATUS='KEEP')
      GOTO 999
889   CLOSE (10,STATUS='KEEP')

      WRITE (21,43)   !  Write header into result file
43    FORMAT (//25('*'),' Dust Bunny Generator ',25('*'))
      CALL DATEPCORD(21,1,IMON,IDAY,IHR,IMIN,ISEC)

! Read option file.  If not present give message and take 25% sample.
      OPEN (UNIT=35,IOSTAT=ICHECK,FILE='DUSTBUNNYOPT.TXT',STATUS='OLD')  ! open import file
      IF (ICHECK.NE.0) THEN  ! file missing or some other problem
         GOTO 10
      ELSE
         WRITE (21,'(/A)') 'Options read from DustBunnyOpt.txt'
         READ (35,*,ERR=10,END=10) DUMMY ! skip first line
         READ (35,*,ERR=10,END=10) NROWSTOMAKE ! second line has number of rows (sample units)to make
         WRITE (21,38) NROWSTOMAKE,' = sample units to make'
38       FORMAT (I10,A)    ! multipurpose formats for reading/writing options
         READ (35,*,ERR=10,END=10) DUMMY 
         READ (35,*,ERR=10,END=10) NCOLSTOMAKE ! number of columns (species) to make
         WRITE (21,38) NCOLSTOMAKE,' = species to make'
         IF (NCOLSTOMAKE.LT.2.OR.NCOLSTOMAKE.GT.10000) GOTO 8

         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) MONTESOURCE ! Source of random number seeds (0=use time, 1=user supplied)
         IF (MONTESOURCE.LT.0.OR.MONTESOURCE.GT.1) GOTO 8
         WRITE (21,48,ERR=10) SEEDSOURCE(MONTESOURCE),' = Random number seed source'
48       FORMAT (A10,A)

         READ (35,*,ERR=10,END=10) DUMMY 
         READ (35,*,ERR=10,END=10) ISEED1            ! Seed for random number generator (positive integer, used only if source = user supplied)
         IF (MONTESOURCE.EQ.1) WRITE (21,38,ERR=10) MONTESEED,' = User-supplied seed'
         READ (35,*,ERR=10,END=10) DUMMY 
         READ (35,*,ERR=10,END=10) RMAX              ! maximum intrinsic growth rate
         WRITE (21,58,ERR=10) RMAX,' = Rmax, maximum intrinsic growth rate'
58       FORMAT (F10.3,A)
         IF (RMAX.LT.0..OR.RMAX.GT.10000.) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) AVEGEN            ! Average number of generations to run model
         WRITE (21,58,ERR=10) AVEGEN,' = Average number of generations to run'
         IF (AVEGEN.LT.1.OR.AVEGEN.GT.100000) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) RANGEGEN          ! Range in number of generations (half above, half below the mean)
         WRITE (21,58,ERR=10) RANGEGEN,' = Range in number of generations (half above, half below the mean)'
         IF (RANGEGEN.LT.0.OR.RANGEGEN.GT.10000) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) CARRYINGCAPACITY  ! Community-level carrying capacity
         WRITE (21,58,ERR=10) CARRYINGCAPACITY,' = Community-level carrying capacity'
         IF (CARRYINGCAPACITY.LT.10.OR.CARRYINGCAPACITY.GT.100000) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) TOLERANCE         ! Spread of species response function to environment
         WRITE (21,58,ERR=10) TOLERANCE,' = Niche width (spread of species response function on environmental gradient)'
         IF (TOLERANCE.LT.1..OR.TOLERANCE.GT.1000.) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) DISPLIMITATION    ! Dispersal limitation factor
         WRITE (21,58,ERR=10) DISPLIMITATION,' = Dispersal limitation factor'
         IF (DISPLIMITATION.LE.0..OR.DISPLIMITATION.GT.10000.) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) NGRADIENTS        ! Number of environmental factors
         WRITE (21,38,ERR=10) NGRADIENTS,' = Number of environmental factors'
         IF (NGRADIENTS.GT.MAXENV) THEN
            WRITE (21,*) 'Error: too many environmentl gradients requested'
            GOTO 999
         ENDIF
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) KWEIGHTS         ! environmental factor weighting
         IF (KWEIGHTS.LT.0.OR.KWEIGHTS.GT.2) GOTO 8
         WRITE (21,48,ERR=10) EWEIGHTING(KWEIGHTS),' = Environmental factor weighting'
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) LOOPS             ! Run loops or not
         IF (LOOPS.NE.0.AND.LOOPS.NE.1) GOTO 8
         WRITE (21,48,ERR=10) YESNO(LOOPS),' = Run loops'

         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) LOOPREPS          ! Reps per innermost loop
         IF (LOOPREPS.LT.1.OR.LOOPREPS.GT.10000) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) LOOPGENS          ! Loop through generations (0/1)
         IF (LOOPGENS.NE.0.AND.LOOPGENS.NE.1) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) LOOPBETA          ! Loop through beta diversity (0/1)
         IF (LOOPBETA.NE.0.AND.LOOPBETA.NE.1) GOTO 8
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) LOOPDISP          ! Loop through dispersal limitation factor (0/1)
         IF (LOOPDISP.NE.0.AND.LOOPDISP.NE.1) GOTO 8
         IF (LOOPS.NE.0) THEN
            WRITE (21,38,ERR=10) LOOPREPS,' = Reps per innermost loop'
            WRITE (21,48,ERR=10) YESNO(LOOPGENS),' = Loop through generations'
            WRITE (21,48,ERR=10) YESNO(LOOPBETA),' = Loop through niche widths'
            WRITE (21,48,ERR=10) YESNO(LOOPDISP),' = Loop through dispersal limitation factor'
            IF (LOOPREPS.GT.1) THEN
               LOOP1BEG=1   ! reps
               LOOP1END=LOOPREPS
               LOOP1INCR=1
            ENDIF
            IF (LOOPGENS.GT.0) THEN
               LOOP2BEG=1   ! generations
               LOOP2END=10
               LOOP2INCR=1
               RANGEGEN=1.
            ENDIF
            IF (LOOPBETA.GT.0) THEN
               LOOP3BEG=15   ! tolerance
               LOOP3END=50
               LOOP3INCR=5
            ENDIF
            IF (LOOPDISP.GT.0) THEN
               LOOP4BEG=1   ! dispersal
               LOOP4END=20
               LOOP4INCR=1
            ENDIF
         ELSE
            LOOPGENS=0
            LOOPBETA=0
            LOOPDISP=0
         ENDIF
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) INTERMEDRESULTS   ! Write intermediate results (0/1)
         WRITE (21,48,ERR=10) YESNO(INTERMEDRESULTS),' = Write intermediate results'
         READ (35,*,ERR=10,END=10) DUMMY
         READ (35,*,ERR=10,END=10) IWRITELOGX   ! Write main matrix as log(x+1) (0/1)
         WRITE (21,48,ERR=10) YESNO(IWRITELOGX),' = Write main matrix as log(x+1)'

         CLOSE (35,STATUS='KEEP')        ! close options file
      ENDIF

      GOTO 12
8     WRITE (21,9)
9     FORMAT (//'Error: Input value in option file is out of range.')
      GOTO 999
10    WRITE (21,11)
11    FORMAT (//&
&     'The options file `DustBunnyOpt.txt` was not found or there was an'/&
&     'error reading it.  Number of rows to make was set to 20, number of columns to 10'/&
&     'Other defaults will be used.'//&
&     'The options file must be named DustBunnyOpt.txt and must be'/&
&     'present in the default directory.')
12    CONTINUE

! Set relative strengths of gradients, according to descending series of broken stick model
         IF (NGRADIENTS.EQ.1) THEN
            STICK(1)=1.
         ELSE
            IF (KWEIGHTS.EQ.0) THEN  ! Uniform weights
               DO K=1,NGRADIENTS
                  STICK(K)=1.
               END DO
            ELSEIF (KWEIGHTS.EQ.1) THEN  ! Descending linear
               DO K=1,NGRADIENTS
                  STICK(K)= 1. - FLOAT(K-1)/FLOAT(NGRADIENTS)
               END DO
            ELSE                     ! Broken stick weights (descending)
               XGRADIENTS=FLOAT(NGRADIENTS)
               DO K=1,NGRADIENTS
                  STICK(K)=0.
                  DO J=K,NGRADIENTS
                     STICK(K)=STICK(K)+1./FLOAT(J)
                  END DO
                  STICK(K)=STICK(K)/FLOAT(NGRADIENTS)
               END DO
            ENDIF
         ENDIF

         WRITE (21,'(/A/A)') 'Relative gradient strengths','Gradient  Rel.Strength'
         DO K=1,NGRADIENTS
            WRITE (21,'(I8,F12.6)') K,STICK(K)
         END DO

      NENT=NROWSTOMAKE
      NATT=NCOLSTOMAKE
      NENV=2+NGRADIENTS
! Allocate memory for main matrix, row/col names, etc.
      LGTH=MAX0(NENT,NATT)
      LGTH2=MAX0(NENT,NENV)
      ALLOCATE (X(NENT,NATT),ENT(NENT),ATT(NATT),VTYPE(LGTH),&
&               ENVOPT(NATT,NGRADIENTS),COMPCOEF(NATT),IMMPRESSURE(NATT),POISSONPARAM(NATT),RINTRINSIC(NATT),&
&               E(NENT,NENV),ENV(NENV),VTYPE2(LGTH2),SPMAX(NATT),SPMAXLOG(NATT),ROWTOT(NENT),ROWTOTLOG(NENT),&
&               WORKVEC(NATT),TEMPVEC(NATT),&
&               STAT=IERR)
      IF (IERR.NE.0) THEN
         CALL FINISWIN(9,1,INUMBER,LABEL)
         GOTO 999
      ENDIF

! more initializations
      GRADIENTLENGTH=140.  ! usable part is from 0 to 100, but extends from -20 to 120
      ENV(1)='SUtotal '
      VTYPE2(1)='Q'
      ENV(2)='Time    '
      VTYPE2(2)='Q'
      DO K=1,NGRADIENTS
         ATEMP='Env     '
          WRITE (ATEMP(4:5),'(I2.2)') K  ! Internal write
         ENV(2+K)=ATEMP
         VTYPE2(2+K)='Q'
      END DO

! initialize random number generator
      IF (MONTESOURCE.EQ.0) THEN
         CALL GETTIM (IHR,IMIN,ISEC,IHUN)
         ISEED1=INT4(ISEC)*100+INT4(IHUN)
      ELSE ! seed is user supplied
         CONTINUE
      ENDIF
      IDUM=(-1)*ISEED1  ! setting IDUM to negative number request initialization
      XDUMMY=RAN1(IDUM,IX1,IX2,IX3)! initialize random number generator

      ICURRENT=1
      IF (LOOPS.EQ.1) THEN
         IBEGINRUNS=1
         NRUNSTOTAL=1
         IF (LOOPREPS.GT.0) NRUNSTOTAL=NRUNSTOTAL*LOOPREPS
         IF (LOOPGENS.EQ.1) NRUNSTOTAL=NRUNSTOTAL*LOOP2END
         IF (LOOPBETA.EQ.1) NRUNSTOTAL=NRUNSTOTAL*(1+(LOOP3END-LOOP3BEG)/LOOP3INCR)
         IF (LOOPDISP.EQ.1) NRUNSTOTAL=NRUNSTOTAL*LOOP4END
         WRITE (21,*) 'Requesting total number of data sets =',NRUNSTOTAL
      ENDIF

      DO IREP=LOOP1BEG,LOOP1END,LOOP1INCR  ! reps
      DO LOOP2=LOOP2BEG,LOOP2END,LOOP2INCR ! generations
         IF (LOOPGENS.EQ.1) AVEGEN=FLOAT(LOOP2)
      DO LOOP3=LOOP3BEG,LOOP3END,LOOP3INCR ! tolerance
         IF (LOOPBETA.EQ.1) TOLERANCE=FLOAT(LOOP3)
      DO LOOP4=LOOP4BEG,LOOP4END,LOOP4INCR ! dispersal
         IF (LOOPDISP.EQ.1) DISPLIMITATION=FLOAT(LOOP4)

         LOOPCOUNTER=LOOPCOUNTER+1
         CALL PROGRWIN(IBEGINRUNS,NRUNSTOTAL,LOOPCOUNTER,ICURRENT,CKESC)
         IF (CKESC.EQ.1) THEN
            CALL FINISWIN(50,0,0,BLANK8)
            GOTO 999
         ENDIF




         DO J=1,NATT  ! Define species parameters
            DO K=1,NGRADIENTS
               ENVOPT(J,K)=150.*FLOAT(J-1)/FLOAT(NATT-1) - 25.0   ! Evenly spaced environmental optima range from -25. to +125
                                                                  !   such that some species optima fall outside of the
                                                                  !   sampled range (0,100)
!               ENVOPT(J,K)=100.*FLOAT(J)/FLOAT(NATT)   ! Evenly spaced environmental optima range (100/NATT, 100)
            END DO
         END DO
         
         IF (NGRADIENTS.GT.1) THEN
            DO K=2,NGRADIENTS  ! Shuffle values in gradients 2 through NGRADIENTS; this makes the env'al gradients independent
               DO J=1,NATT     ! Copy data into vector
                  TEMPVEC(J)=ENVOPT(J,K)
               END DO
               CALL SHUFFVAR(TEMPVEC,WORKVEC,NATT,IX1,IX2,IX3)  ! shuffle this vector
               DO J=1,NATT     ! Copy shuffled vector back into array
                  ENVOPT(J,K)=TEMPVEC(J)
               END DO
            END DO
         ENDIF

         DO J=1,NATT  ! Define species parameters
            VTYPE(J)='Q'
            ATEMP=BLANK8                   ! Make column label
            WRITE (ATEMP(3:6),'(I4.4)') J  ! Internal write
            ATEMP(1:2)='Sp'
            ATT(J)=ATEMP
            ATEMP(7:7)=BLANK1
            ATEMP(8:8)=BLANK1

325         COMPCOEF(J)=RAN1(1,IX1,IX2,IX3)        ! Get random number for competitive ability, range (0,1)
            IF (COMPCOEF(J).LE.0.) GOTO 325        ! Avoid later division by zero, get another random number
            POISSONPARAM(J)= -1.*COMPCOEF(J) + 1.  ! Assign poisson parameter; is negative linear function of COMPCOEF
            POISSONPARAM(J)= POISSONPARAM(J)/DISPLIMITATION
            RMIN=1.                                ! Minimum intrinsic growth rate
            RINTRINSIC(J)=-1.*(RMAX-RMIN)*COMPCOEF(J) + RMAX   ! Intrinsic growth rate is negative linear function of COMPCOEF
         END DO

         IF (LOOPS.EQ.0) THEN
            WRITE (21,340)
340         FORMAT (/'Species parameters'/&
&                    '  Sp   CompCoef PoissParam          R')
            DO J=1,NATT
               WRITE (21,345) J,COMPCOEF(J),POISSONPARAM(J),RINTRINSIC(J)
345            FORMAT (I4,3F11.4)
            END DO
            
            WRITE (21,350) (ENV(2+K), K=1,NGRADIENTS)
350         FORMAT (/'Species parameters, environmental optima'/&
&                    '  Sp      ',10(A8,3x))
            DO J=1,NATT
               WRITE (21,352) J,(ENVOPT(J,K),K=1,NGRADIENTS)
352            FORMAT (I4,10F11.4)
            END DO
         ENDIF

         PROCESS='Generating communities...'
         PROCESS(26:26)=CHAR(0)
         CALL UPDATE(PROCESS,CKESC)
         IF (CKESC.EQ.1) THEN
            CALL FINISWIN(50,0,0,BLANK8)
            GOTO 999
         ENDIF

         IF (INTERMEDRESULTS.EQ.1) WRITE (21,355)
355      FORMAT (/'Intermediate Results'&
&                /'  SU   t  Sp    X(I,J)       Z    GROWTH    IMMIGR     Total')

         IBEGIN=1
         IEND=NENT

         DO I=1,NENT  ! Cycle through rows to be generated

            IF (LOOPS.EQ.0) CALL PROGRWIN(IBEGIN,IEND,I,ICURRENT,CKESC)
            IF (CKESC.EQ.1) THEN
               CALL FINISWIN(50,0,0,BLANK8)
               GOTO 999
            ENDIF
            DO K=1,NGRADIENTS
               E(I,K+2)=100.*RAN1(1,IX1,IX2,IX3)   ! Choose a random point on the gradient
               ETARGET(K)=E(I,K+2)
            END DO
            NGENERATIONS=NINT(AVEGEN-RANGEGEN/2. + RANGEGEN*RAN1(1,IX1,IX2,IX3))
            E(I,2)=FLOAT(NGENERATIONS)  !  Record this in second matrix as a temporal variable.

                                         ! Make ROW label
            ATEMP=BLANK8
            WRITE (ATEMP(3:6),'(I4.4)') I  ! Internal write
            ATEMP(1:2)='SU'
            ATEMP(7:7)=BLANK1
            ATEMP(8:8)=BLANK1
            ENT(I)=ATEMP

            DO J=1,NATT
               X(I,J)=0.          ! Initialize each species at zero
            END DO

            DO KG=1,NGENERATIONS   ! Colonize this sample unit
               COMMTOTAL=0.
               DO J=1,NATT        ! Get community total
                  COMMTOTAL=COMMTOTAL+X(I,J)
               END DO

               DO J=1,NATT  ! Loop through species for this time step

                                                 ! Calculate environmental effect on species growth (range: 0 = failure, 1 = optimum)
                  Z=1.
                  DO K=1,NGRADIENTS
                     EVALUE=ETARGET(K)
                     RAWEFFECT=GAUSSRESP(TOLERANCE,ENVOPT(J,K),EVALUE)  ! Effect of evalue on species, according to
                     EFFECT=1.-(1.-RAWEFFECT)*STICK(K)                  ! Moderate that effect by moving it towards 1, 
                                                                        !    moving it farther for weaker env'al variables.
                                                                        ! Env'al factors with full effect have stick(k)=1
                     Z=Z*EFFECT  ! Accumulate environmental effects multiplicatively
                  END DO

                  L=0
                  POISZERO=POISSON(L,LMAX,POISSONPARAM(J))
                  SUMNONZEROP=1.-POISZERO
                  ADJPOISZERO=1.-Z*SUMNONZEROP
                  CUMULPOISSON=ADJPOISZERO
!                  write (21,*) poiszero,'=poiszero ',sumnonzerop,'=sumnonzerop ',adjpoiszero,'=adjpoiszero ',cumulpoisson,'=cumulpoisson '

                  XRAND=RAN1(1,IX1,IX2,IX3)   ! Get random number to determine immigration
                                              ! Find where the random number falls in the Poisson series
                                              ! Get number of immigrants
                  IMMIGRANTS=15
                  IF (XRAND.LT.CUMULPOISSON) THEN ! random number falls in FIRST interval (L=0),
                     IMMIGRANTS=0                 ! so Immigrants=0
                     GOTO 765
                  ENDIF
                  DO L=1,LMAX
                     CUMULPOISSON=CUMULPOISSON + Z * POISSON(L,LMAX,POISSONPARAM(J))
!                     WRITE (21,*) CUMULPOISSON,'=CUMULPOISSON',XRAND,'=XRAND'
                     IF (XRAND.LE.CUMULPOISSON) THEN ! random number falls in interval, so assign number of immigrants as L
                        IMMIGRANTS=L
                        GOTO 765
                     ELSE
                        CYCLE
                     ENDIF
                  END DO
                  WRITE (21,*) 'Warning: Immigrants = ',LMAX,' was reached. Possible problem with routine.',CUMULPOISSON,'=CUM.POIS.',XRAND,'=XRAND'
                  IMMIGRANTS=LMAX
765               CONTINUE

                                                 ! Calculate population growth
                  POPGROWTH= RINTRINSIC(J) * Z * X(I,J) * (1. - COMMTOTAL/(COMPCOEF(J)*CARRYINGCAPACITY))

                  IF (INTERMEDRESULTS.EQ.1) WRITE (21,555,ERR=556) I,KG,J,X(I,J),Z,POPGROWTH,IMMIGRANTS,COMMTOTAL
555               FORMAT (3I4,3F10.3,I10,F10.3)
556               CONTINUE

                  X(I,J) = X(I,J) + POPGROWTH + FLOAT(IMMIGRANTS)  ! Calculate new population size
                  IF (X(I,J).LT.0.) X(I,J)=0.    ! Population cannot be negative

               END DO ! End species loop
            END DO    ! End time step loop
         END DO       ! End loop through sample units.

 ! Get matrix mean, number of zeros, and row totals
         NZEROS=0
         DO J=1,NATT
            SPMAX(J)=0.
            SPMAXLOG(J)=0.
            DO I=1,NENT
               IF (X(I,J).GT.SPMAX(J)) SPMAX(J)=X(I,J)
               IF (X(I,J).LE.0.) NZEROS=NZEROS+1
            END DO
            IF (SPMAX(J).GT.0) SPMAXLOG(J)=ALOG10(SPMAX(J)+1.)
         END DO
         
         NEMPTYSPP=0   ! Count the number of empty spp to later exclude these from the %zeros
         DO J=1,NATT
            IF (SPMAX(J).LE.0) NEMPTYSPP=NEMPTYSPP+1
         END DO
         NZEROS=NZEROS - NEMPTYSPP*NENT
         NATTADJ=NATT-NEMPTYSPP  ! Number of species adjusted after empties removed

         NNONZEROS=NENT*NATTADJ-NZEROS
         RELMATRIXSUM=0.
         RELMATRIXSUMLOG=0.
         NEMPTYSU=0
         
         DO I=1,NENT
            ROWTOT(I)=0.
            ROWTOTLOG(I)=0.
            DO J=1,NATT
               ROWTOT(I)=ROWTOT(I)+X(I,J)
               IF (X(I,J).GE.0) ROWTOTLOG(I)=ROWTOTLOG(I)+ ALOG10(X(I,J)+1.)
               IF (SPMAX(J).GT.0.) THEN
                  RELMATRIXSUM=RELMATRIXSUM+X(I,J)/SPMAX(J)
                  IF (SPMAXLOG(J).GT.0) THEN
                     RELMATRIXSUMLOG=RELMATRIXSUMLOG+ALOG10(X(I,J)+1)/SPMAXLOG(J)
                  ENDIF
               ENDIF
            END DO
            IF (ROWTOT(I).LE.0.) NEMPTYSU=NEMPTYSU+1
            E(I,1)=ROWTOT(I)
         END DO
         RELMATRIXMEAN=1. - RELMATRIXSUM/FLOAT(NENT*NATTADJ)
         RELMATRIXMEANLOG=1. - RELMATRIXSUMLOG/FLOAT(NENT*NATTADJ)
         PCTZEROS=100.*FLOAT(NZEROS)/FLOAT(NENT*NATTADJ)

         IF (IWRITELOGX.NE.1.AND.LOOPS.NE.1) THEN
            CALL WRITEMAIN(NENT,NATT,X,ENT,ATT,VTYPE,AROWS,ACOLS,LGTH,IERR)
            IF (IERR.NE.0) GOTO 999
         ENDIF

! Calculate averaged sorensen distance among rows and beta diversity
         CALL AVESORDIST(X,NENT,NATT,ATT,AVEDISS,IERR)
         IF (IERR.NE.0) THEN
            AVEHALFCH=XMISSING
            AVEDISS=XMISSING
                        ! Delete the error file if there was one
            INQUIRE(FILE='PCORD.ERR',EXIST=FILEEXIST)
            IF (FILEEXIST) THEN
               OPEN (9,FILE='PCORD.ERR',STATUS='UNKNOWN')
               CLOSE (9,STATUS='DELETE')
            ENDIF
         ELSE         ! Calculate average distance in half changes
            IF (AVEDISS.GE.1.) THEN
               AVEHALFCH=0.
            ELSE
               AVEHALFCH=ALOG10(1.-AVEDISS)/ALOG10(0.5)
            ENDIF
            AVEDISS=100.*AVEDISS
         ENDIF
!         BETAW=FLOAT(NATT)/(FLOAT(NNONZEROS)/FLOAT(NENT)) - 1.
         BETAW=FLOAT(NATTADJ*NENT)/FLOAT(NNONZEROS) - 1.

! Repeat the above with log-transformed data
         DO I=1,NENT   ! first log transform
            DO J=1,NATT
               IF (X(I,J).LT.0) THEN
                  X(I,J)=0.
               ELSE
                  X(I,J) = ALOG10(X(I,J)+1.)
               ENDIF
            END DO
         END DO
         CALL AVESORDIST(X,NENT,NATT,ATT,AVEDISSLOG,IERR)
         IF (IERR.NE.0) THEN
            AVEHALFCHLOG=XMISSING
            AVEDISSLOG=XMISSING
                        ! Delete the error file if there was one
            INQUIRE(FILE='PCORD.ERR',EXIST=FILEEXIST)
            IF (FILEEXIST) THEN
               OPEN (9,FILE='PCORD.ERR',STATUS='UNKNOWN')
               CLOSE (9,STATUS='DELETE')
            ENDIF
         ELSE         ! Calculate average distance in half changes
            IF (AVEDISSLOG.GE.1.) THEN
               AVEHALFCHLOG=0.
            ELSE
               AVEHALFCHLOG=ALOG10(1.-AVEDISSLOG)/ALOG10(0.5)
            ENDIF
            AVEDISSLOG=100.*AVEDISSLOG
         ENDIF
         AVERICHNESS= NATTADJ*(1. - PCTZEROS/100.)

        IF (LOOPS.NE.1) THEN
            WRITE (21,'(/F10.4,A)') PCTZEROS,' = percentage of zeros in matrix, excluding empty species'
            WRITE (21,'(F10.4,A)') BETAW,' = beta diversity, Whittaker`s'
            WRITE (21,'(/A)') 'Statistics on raw abundance data'
            WRITE (21,'(F10.4,A)') RELMATRIXMEAN,' = 1- matrix mean after relativization by species maximum, excluding empties'
            WRITE (21,'(F10.4,A)') AVEDISS,' = average Sorensen distance among rows'
            WRITE (21,'(F10.4,A)') AVEHALFCH,' = beta diversity, average half changes'
            WRITE (21,'(/A)') 'Statistics after log(x+1) transformation of data'
            WRITE (21,'(F10.4,A)') RELMATRIXMEANLOG,' = 1- matrix mean after relativization by species maximum, excluding empties'
            WRITE (21,'(F10.4,A)') AVEDISSLOG,' = average Sorensen distance among rows'
            WRITE (21,'(F10.4,A)') AVEHALFCHLOG,' = beta diversity, average half changes'
            WRITE (21,'(F10.4,A)') AVERICHNESS,' = average species richness'

            IF (IWRITELOGX.EQ.1) THEN
               CALL WRITEMAIN(NENT,NATT,X,ENT,ATT,VTYPE,AROWS,ACOLS,LGTH,IERR)
               IF (IERR.NE.0) THEN
                  WRITE (21,*) 'Error writing main matrix.'
                  GOTO 999
               ENDIF
            ENDIF

                                                   ! Save second matrix
            OUTPUTFILE='TEMP2.WK1'
            HANDL=18
            OPEN (HANDL,FILE=OUTPUTFILE,STATUS='UNKNOWN',FORM='BINARY',BLOCKSIZE=2048,BUFFERED='YES',&
&                 BUFFERCOUNT=30,IOSTAT=IERR)
               IF (IERR.NE.0) THEN
               CALL FINISWIN(10,11,INUMBER,LABEL)
               GOTO 999
            ENDIF

            NORMAL=1
            IHEADER=1
            CALL WRWK15 (HANDL,NENT,NENV,E,ENT,ENV,VTYPE2,AROWS,AENV,LGTH2,NORMAL,IHEADER,IERR)
            IF (IERR.NE.0) GOTO 999
            INQUIRE(18,OPENED=FILEOPENED)
            IF (FILEOPENED)  CLOSE (18,STATUS='KEEP')

            WRITE (21,807) NENT,NATT,NEMPTYSU,NEMPTYSPP,NATTADJ,NENV
807         FORMAT (&
&                    /I10,' sample units generated.'&
&                    /I10,' species generated.'&
&                    /I10,' sample units were empty (all zeros).'&
&                    /I10,' species were empty (all zeros).'&
&                    /I10,' species with at least one positive value.'&
&                    /I10,' variables generated in second matrix.'&
&                    //32('*'),' Procedure completed ',32('*'))

         ELSE  ! write line in output file then continue with loops
            IF (ITBLHEADER.EQ.1) THEN
               WRITE (21,847,ERR=990)
847            FORMAT(/'     Run     SUs     Spp SppPres  NGrads    NGen Tolerance DispLimit       DBI    PctZeros AveSorDist  BetaDivHC  BetaDivWh     DBIlog AveDistLog  BetaHClog  NemptySUs')
               ITBLHEADER=0
            ENDIF
            WRITE (21,857,ERR=995) LOOPCOUNTER,NENT,NATT,NATTADJ,NGRADIENTS,NGENERATIONS,TOLERANCE,DISPLIMITATION,&
&                                  RELMATRIXMEAN,PCTZEROS,AVEDISS,AVEHALFCH,BETAW,&
&                                  RELMATRIXMEANLOG,AVEDISSLOG,AVEHALFCHLOG,AVERICHNESS,NEMPTYSU
857         FORMAT (6I8,2F10.1,F11.4,8F11.3,I8)
         ENDIF
      END DO  ! End of dispersal loop
      END DO  ! End of tolerance loop
      END DO  ! End of generations loop
      END DO  ! End of main loop thru reps

      GOTO 999
990   WRITE (21,*) 'ERROR WRITING OUTPUT HEADER'
      GOTO 999
995   WRITE (21,*) 'ERROR WRITING OUTPUT LINE'

999   CONTINUE
      INCLUDE 'FREELIB.F90'

!  Clean up files
      INQUIRE( 9,OPENED=FILEOPENED)
      IF (FILEOPENED) CLOSE ( 9,STATUS='KEEP',IOSTAT=IERR) ! error message file
      INQUIRE(10,OPENED=FILEOPENED)
      IF (FILEOPENED) CLOSE (10,STATUS='KEEP',IOSTAT=IERR) ! options file
      INQUIRE(18,OPENED=FILEOPENED)
      IF (FILEOPENED) CLOSE (18,STATUS='KEEP',IOSTAT=IERR) ! output spreadsheet
      INQUIRE(21,OPENED=FILEOPENED)
      IF (FILEOPENED) CLOSE (21,STATUS='KEEP',IOSTAT=IERR) ! result file
      INQUIRE(26,OPENED=FILEOPENED)
      IF (FILEOPENED) CLOSE (26,STATUS='KEEP',IOSTAT=IERR) ! input matrices
      INQUIRE(35,OPENED=FILEOPENED)
      IF (FILEOPENED) CLOSE (35,STATUS='KEEP',IOSTAT=IERR) ! DustBunnyOpt.txt
      END

!******************************************************************************
      FUNCTION RAN1(IDUM,IX1,IX2,IX3)   ! From Numerical Recipes, pp. 196-197.
! Returns a uniform random deviate between 0.0 and 1.0. 
! Set IDUM to any negative value to initialize or re-initialize the sequence.
! Set IDUM to any dummy positive integer in subsequent calls.
! This routine is based on three linear congruential generators.  One generator
! is used for the most significant part of the output number, the second for 
! the least significant part, and the third to control a shuffling routine.
! This is the best (and slowest) of the random number generators
! given in the Numerical Recipes book (Press, W. H., et al.  1986.  Numerical
! Recipes.  Cambridge Univ. Press, Cambridge.)

      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)  ! These parameters good for
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)  ! machine with integer overflow
      PARAMETER (M3=243000,IA3=4561,IC3=51349)                   ! at 2**30
      DATA IFF /0/               ! Initialize on first call even if IDUM is not negative
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)     ! Seed the first routine
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)          ! and use it to seed the second
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)          ! and third routines
        DO J=1,97                ! Fill the table with sequential uniform deviates
          IX1=MOD(IA1*IX1+IC1,M1)!   generated by the first two routines.
          IX2=MOD(IA2*IX2+IC2,M2)
          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1   !combine low and high order pieces
        END DO
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)    ! Generate the next number for each sequence
      IX2=MOD(IA2*IX2+IC2,M2)
12    IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3            ! Use the 3rd sequence to get an integer between 1 and 97
      IF(J.GT.97.OR.J.LT.1) GOTO 12  ! goto replace PAUSE in original code
      RAN1=R(J)                  ! Return that table entry
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1   ! and refill it with a new number
      RETURN
      END

!*****************
FUNCTION GAUSSRESP(TOLERANCE,OPTIMUM,ETARGET)
! Get get coefficient for environmental effect (Z) at ETARGET
!    given optimum and tolerance (s.d.)
! GAUSSRESP=1 when ETARGET lies on the optimum

    REAL*4 GAUSSRESP,DEV2,DIFF,TWONEG
    
    TWONEG=-2.0
    IERR=0
    WEIGHT=DONE

    DIFF=ABS(ETARGET-OPTIMUM)  ! Absolute distance between target point and OPTIMUM
    DEV2=((DIFF/TOLERANCE)**2)/TWONEG
!  WRITE (21,*) DEV2,'=DEV2',DIFF,'=DIFF'
    IF (DEV2.GE.1.) THEN
       GAUSSRESP=1.
    ELSE
       GAUSSRESP=EXP(DEV2)
    ENDIF
! WRITE (21,*) GAUSSRESP,TOLERANCE,OPTIMUM,ETARGET,'IN GAUSSRESP'

    RETURN
    END

!*************************************************************
FUNCTION POISSON(L,LMAX,POISPARAM)
    REAL*4 FACTORIALS(0:15)
    DATA FACTORIALS/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3628800.,399168000.,4.790016E+08,6.227021E+09,&
&                   8.717829E+10,1.307674E+12/

    POISSON=0.
    IF (POISPARAM.EQ.0.OR.L.LT.0) THEN
       POISSON=0.
    ELSEIF (L.GT.LMAX) THEN
       POISSON=0.
    ELSE
       POISSON=(POISPARAM**L) * EXP(-1.*POISPARAM) / FACTORIALS(L)
    ENDIF
    
    RETURN
    END
!*************************************************************
    SUBROUTINE SHUFFVAR(VEC1,VEC2,N,IX1,IX2,IX3) ! Shuffle array VEC1 of length NENT
    REAL VEC1(N),VEC2(N)

    DO I=1,N         ! Clear array used to mark reassigned elements
       VEC2(I)=0.
    END DO

    IPOINT=N+1
    DO I=1,N
       IPOINT=IPOINT-1
80     XRAND=RAN1(1,IX1,IX2,IX3)             ! get random number
       IRND = 1 + IFIX(FLOAT(IPOINT)*XRAND)  ! convert to pointer ranging up to number of remaining items
       IF (IRND.GT.IPOINT) GOTO 80
       VEC2(I)=VEC1(IRND)  ! move random item from GROUP array into next position in IGROUP
       IF (IPOINT.NE.IRND) VEC1(IRND)=VEC1(IPOINT)  ! move last item into hole just created
    END DO

    DO I=1,N  ! now move shuffled numbers back into original array
       VEC1(I)=VEC2(I)
    END DO

    RETURN
    END

!*************************************************************
SUBROUTINE WRITEMAIN(NENT,NATT,X,ENT,ATT,VTYPE,AROWS,ACOLS,LGTH,IERR) ! Save main matrix
    LOGICAL FILEOPENED
    INTEGER*2 HANDL
    CHARACTER*9 OUTPUTFILE
    CHARACTER*8 ENT(NENT),ATT(NATT),AROWS,ACOLS
    CHARACTER*1 VTYPE(LGTH)
    DIMENSION X(NENT,NATT)

    INUMBER=0
    IERR=0
    HANDL=18
    OUTPUTFILE='TEMP.WK1'
    OPEN (HANDL,FILE=OUTPUTFILE,STATUS='UNKNOWN',FORM='BINARY',BLOCKSIZE=2048,BUFFERED='YES',&
&         BUFFERCOUNT=30,IOSTAT=IERR)
    IF (IERR.NE.0) THEN
       IERR=10
       CALL FINISWIN(IERR,11,INUMBER,'        ')
       RETURN
    ENDIF

    IHEADER=1
    NORMAL=1
    CALL WRWK15 (HANDL,NENT,NATT,X,ENT,ATT,VTYPE,AROWS,ACOLS,LGTH,NORMAL,IHEADER,IERR)
    INQUIRE(18,OPENED=FILEOPENED)
    IF (FILEOPENED)  CLOSE (18,STATUS='KEEP')
    RETURN
    END

!******************** end of file DustBunnyOpt.f90 *************************
