PROGRAM ZELIG
C
C     ZELIG version 1.0
C     D.L. Urban, June 1990
C
C          ZELIG is a transect version of FORET.
C          This version simulates canopies by distributing
C     each tree's leaf area vertically along it's bole (the
C     approach is taken from the FORSKA model of Rik Leemans).
C     In default mode, shading effects are implemented by averaging
C     available light over each tree's canopy (in 1-m increments).
C     In interactive mode, the light regime is partitioned into
C     diffuse and direct-beam radiation, and shading is based
C     on sun angle and tree height.
C          Below-ground constraints include soil fertility (as
C     maximum biomass increment) and soil moisture (as drought-
C     days per growing season, based on a Thornthwaite-Mather
C     routine with a single-layer soil, from the Pastor-Post model).
C
C     The program is modular and is designed for easy modification.
C
C
C     MAIN:
C
C     PARAMETERS TO DIMENSION ZELIG:
C     MP=MAX PLOTS; MT=MAX TREES PER PLOT; MS=MAX SPECIES;
C     MH=MAX HEIGHT OF FOLIAGE PROFILES.
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
C     MAIN LOOP:
      CALL INITL
      DO 1 KYR=1,NYRS
      CALL WEATHR
      CALL MORTAL
      CALL GROW
      CALL REGEN
      CALL BOOKS
      IF (MOD(KYR,IPRT).EQ.0) CALL PRINT
      IF (MOD(KYR,IPCH).EQ.0) CALL PUNCH
    1 CONTINUE
      STOP
      END
C
C     SUBROUTINES:
C
C     INITL INITIALIZES PARAMETERS AND READS INPUT DATA.
C     INPUT DATA ARE ECHO-PRINTED TO THE SCREEN.
C
      SUBROUTINE INITL
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      CHARACTER*72 TITLE
      CHARACTER*24 LOCALE
C     I/O FILES:
C     SOME INFO IS ECHO-PRINTED TO UNIT 3, THE SCREEN;
      OPEN(3,FILE='ECHO')
C     UNIT 4 IS AN ANNUAL TRACER ON BASAL AREA PER SPECIES:
      OPEN(4,FILE='TRACER')
C     UNIT 5 IS INPUT DRIVER DATA:
      OPEN(5,FILE='DRIVER')
C     UNIT 6 IS IBM RUN-TIME DIAGNOSTICS AND SYSTEM PRINT.
C     UNIT 7 IS PUNCHED PLOT-LEVEL OUTPUT:
      OPEN(7,FILE='PUNCH')
C     UNIT 8 IS PRINTED STAND-LEVEL OUTPUT:
      OPEN(8,FILE='PRINT')
C     UNIT 9 IS INITIAL-CONDITION DATA, IF USED:
      OPEN(9,FILE='INDATA')
C     UNIT 10 IS OUTPUT LEAF-AREA PROFILE:
      OPEN(10,FILE='PROFILE')
C     UNIT 11 IS USER-PROVIDED DIAGNOSTICS, IF USED:
      OPEN(11,FILE='LOG')
C     READ RUN TITLE AND CONTROL PARAMETERS
      READ(5,5001) TITLE, MODE, INDATA, NPLOTS, NSPP,
     2  NYRS, IPRT, IPCH, ILOG
 5001 FORMAT(A72/8I5)
C     READ SITE DATA
      READ(5,5002) LOCALE, PLAT, PLON, ELEV, THETA, PHIB, PHID,
     2  XCHT, FC, WP, SF
 5002 FORMAT(A24/2F6.1,F7.1/3F6.3,2X,F2.0/3F6.2)
C     COMPUTE PLOT SIZE FROM SUN ANGLE AND EXPECTED CANOPY HEIGHT:
C     SHADOW LENGTH FOR CANOPY TREE:
      SL=XCHT/TAN(THETA)
C     NO. OF CELLS TO AGGREGATE:
      NAC=(SL/10.0+0.5)+1
C     MXS = HT (M) OF DIAGONAL X.S. THRU A PLOT AT THIS THETA:
      MXS=10.0*TAN(THETA)+0.5
C     ADJUST FOR INDEPENDENT VS INTERACTIVE-PLOT MODE:
      IF (MODE.EQ.0) THEN
        AREA=100.0*NAC
        XAREA=FLOAT(NPLOTS)*AREA/10000.0
        NMAX=100*NAC
        WRITE(8,8001) AREA
        WRITE(3,8001) AREA
        ELSE
        AREA=100.0
        NMAX=100
        XAREA=0.01*NPLOTS
        AA=100.0*NAC
        WRITE(8,8002) AA, NAC, MXS
        WRITE(3,8002) AA, NAC, MXS
        ENDIF
 8001 FORMAT(/5X,'ZELIG IS IN INDEPENDENT-PLOT MODE',
     2  /10X,'PLOT AREA:',F8.1,' SQ. M')
 8002 FORMAT(/5X,'ZELIG IS IN INTERACTIVE-TRANSECT MODE',
     2  /10X,'AGGREGATE AREA:',F8.1,' SQ. M (',I2,' CELLS)',
     3  /10X,'VERTICAL STEP SIZE THRU LEAF PROFILE:',I3,' M')
C     CHECK DIMENSIONS FOR MODE 0 VS 1 ...
      IF (NMAX.GT.MT) THEN
        WRITE(3,9999)
        STOP
        ENDIF
 9999 FORMAT(//5X,'ZELIG IS OUT-OF-BOUNDS:',
     2  '  REDIMENSION AND RESTART'//)
C     CORRECT MAX BIOMASS INCREMENT TO PLOT SIZE (MG/HA->KG/AREA)
      SF=SF*AREA/10.0
      WRITE(8,8003) TITLE, LOCALE, PLAT, PLON, ELEV, THETA, PHIB,
     2  PHID, SF, FC, WP
      WRITE(3,8003) TITLE, LOCALE, PLAT, PLON, ELEV, THETA, PHIB,
     2  PHID, SF, FC, WP
 8003 FORMAT(/5X,A72,
     2  //5X,'LOCATION:  ',A24,/10X,'LAT:',F6.1,5X,'LONG:',F6.1,
     3  5X,'ELEV:',F7.1,' M',/10X,'SOLAR ANGLE:',F7.2,' RAD',5X,
     4  'DIRECT, DIFFUSE (.%):',F6.2,',',F5.2,/10X,'SOIL FERTILITY',
     5  ' (KG):',F8.2,5X,'FC, WP (CM):',F7.2,',',F6.2)
C     READ TEMPERATURE AND PRECIP DATA (MEANS AND S.D.'S):
      READ(5,5003) (XT(M),M=1,12), (VT(M),M=1,12),
     2  (XR(M),M=1,12), (VR(M),M=1,12)
 5003 FORMAT(12F5.1)
C     CORRECT LATITUDE POINTER FOR SR WEATHR
      LAT=(PLAT+0.5)-30
C     IDUM IS THE SEED FOR THE RANDOM NUMBER GENERATOR
      IDUM=-1
C     INITIALIZE ARRAYS
      DO 11 KP=1,MP
      IND(KP)=0
      BIOM(KP)=0.0
C     SOIL FERTILITY FACTORS:
      SFF(KP,1)=1.0
      SFF(KP,2)=1.0
      SFF(KP,3)=1.0
      CLAI(KP)=0.0
      DO 111 KI=1,MT
      DBH(KP,KI)=0.0
      IBC(KP,KI)=0
      ISP(KP,KI)=0
      NOGRO(KP,KI)=0
  111 CONTINUE
      DO 112 M=1,MH
      AL(KP,M)=1.0
  112 CONTINUE
   11 CONTINUE
      WRITE(8,8004) NPLOTS, XAREA, NSPP
      WRITE(3,8004) NPLOTS, XAREA, NSPP
 8004 FORMAT(/10X,'NUMBER OF PLOTS OR CELLS:',I5,
     2  /10X,'OUTPUT SAMPLES ARE',F6.3,'-HA AGGREGATES',
     3  //5X,'NUMBER OF SPECIES IN SIMULATION:',I5)
C     READ AND COMPUTE SPECIES PARAMETERS ...
C     NEGATING THE FIRST SPECIES' SEED RATE SETS THEM ALL TO 1.0
      KSEED=0
C     LSPP IS LOCAL SPECIES INCLUDED (NON-0 SEED)
      LSPP=0
C     SEED(K) IS SET FROM ASSUMPTION OF 1 STEM/SQ.M, AND A 10-YR
C     STOCKING TIME; RANK SEEDING RATES ARE THEN ADJUSTED TO THIS.
      SMAX=AREA/10.0
      STOT=0.0
      DO 12 K=1,NSPP
      READ(5,5004) MSP(K), AGEMAX(K), DMAX(K), HTMAX(K), G(K),
     2  DDMIN(K), DDMAX(K), LIGHT(K), MDRT(K), NUTRI(K), SEED(K)
 5004 FORMAT(A4,1X,F4.0,F4.0,F5.0,2X,F3.0,1X,2F5.0,1X,3I2,1X,F2.0)
      IF (K.EQ.1.AND.SEED(K).LT.0.0) KSEED=1
      IF (KSEED.EQ.1) SEED(K)=1.0
      STOT=STOT+SEED(K)
      IF (SEED(K).GT.0.0) LSPP=LSPP+1
      B2(K)=2.0*((HTMAX(K)-137.0)/DMAX(K))
      B3(K)=(HTMAX(K)-137.0)/DMAX(K)**2
   12 CONTINUE
C     NUMBER OF SPECIES ALLOWED TO ESTABLISH:
      WRITE(3,8005) LSPP
      WRITE(8,8005) LSPP
 8005 FORMAT(/10X,'NUMBER LOCALLY AVAILABLE:',I5/,
     2  /10X,'SPECIES ',
     3  'AGE, DBH, HT; G; GDDS; LIGHT, DROUT, NUTRI; SEEDS:'/)
C     ADJUST SEED(K) TO PLOT SIZE AND RELATIVIZE PER SPECIES:
      DO 13 K=1,NSPP
      IF (SEED(K).EQ.0.0) GO TO 13
      SEED(K)=SMAX*SEED(K)/STOT
      WRITE(8,8006) MSP(K), AGEMAX(K), DMAX(K), HTMAX(K), G(K),
     2  DDMIN(K), DDMAX(K), LIGHT(K), MDRT(K), NUTRI(K), SEED(K)
      WRITE(3,8006) MSP(K), AGEMAX(K), DMAX(K), HTMAX(K), G(K),
     2  DDMIN(K), DDMAX(K), LIGHT(K), MDRT(K), NUTRI(K), SEED(K)
 8006 FORMAT(10X,A4,1X,F6.1,F6.1,F7.1,2X,F5.1,2X,F6.1,F7.1,1X,3I2,F5.1)
   13 CONTINUE
      IF (INDATA.EQ.0) THEN
        WRITE(8,8007)
        WRITE(3,8007)
        RETURN
        ENDIF
C     CONTINUE IF INDATA=1 AND EXTERNAL INITIALS ARE USED:
      WRITE(8,8008)
      WRITE(3,8008)
 8007 FORMAT(/5X,'SIMULATION INITIATED FROM BARE GROUND')
 8008 FORMAT(/5X,'SIMULATION INITIATED WITH DATA FROM UNIT 9')
C     READ EXTERNAL INITIAL PLOT DATA
      DO 14 KP=1,NPLOTS
      READ(9,9001) IND(KP)
 9001 FORMAT(10X,I5)
      DO 15 KI=1,IND(KP)
      READ(9,9002) ISP(KP,KI), DBH(KP,KI), IBC(KP,KI)
 9002 FORMAT(5X,I5,F10.5,I5)
      NOGRO(KP,KI)=0
   15 CONTINUE
   14 CONTINUE
      KYR=0
      CALL BOOKS
      CALL PRINT
      RETURN
      END
C
C     WEATHR GENERATES TEMP., DEGD, PRECIP, AND SOIL MOISTURE.
C     THE SOIL MOISTURE ROUTINE IS MODIFIED FROM PASTOR AND POST,
C     ORNL/TM-9519, 1985.
C
      SUBROUTINE WEATHR
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      REAL DAYS(12), CLAT(12,20), T(12), R(12), PET, DDM
      DATA DAYS/31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./
C     FOLLOWING ARE PET CORRECTION TERMS FOR 31-50 N LATITUDE
      DATA CLAT/.90,.87,1.03,1.08,1.18,1.18,1.20,1.14,1.03,.98,.89,.88,
     2          .89,.86,1.03,1.08,1.19,1.19,1.21,1.15,1.03,.98,.88,.87,
     3          .88,.86,1.03,1.09,1.19,1.20,1.22,1.15,1.03,.97,.88,.86,
     4          .88,.85,1.03,1.09,1.20,1.20,1.22,1.16,1.03,.97,.87,.86,
     5          .87,.85,1.03,1.09,1.21,1.21,1.23,1.16,1.03,.97,.86,.85,
     6          .87,.85,1.03,1.10,1.21,1.22,1.24,1.16,1.03,.97,.86,.84,
     7          .86,.84,1.03,1.10,1.22,1.23,1.25,1.17,1.03,.97,.85,.83,
     8          .85,.84,1.03,1.10,1.23,1.24,1.25,1.17,1.04,.96,.84,.83,
     9          .85,.84,1.03,1.11,1.23,1.24,1.26,1.18,1.04,.96,.84,.82,
     *          .84,.83,1.03,1.11,1.24,1.25,1.27,1.18,1.04,.96,.83,.81,
     1          .83,.83,1.03,1.11,1.25,1.26,1.27,1.19,1.04,.96,.82,.80,
     2          .82,.83,1.03,1.12,1.26,1.27,1.28,1.19,1.04,.95,.82,.79,
     3          .81,.82,1.02,1.12,1.26,1.28,1.29,1.20,1.04,.95,.81,.77,
     4          .81,.82,1.02,1.13,1.27,1.29,1.30,1.20,1.04,.95,.80,.76,
     5          .80,.81,1.02,1.13,1.28,1.29,1.31,1.21,1.04,.94,.79,.75,
     6          .79,.81,1.02,1.13,1.29,1.31,1.32,1.22,1.04,.94,.79,.74,
     7          .77,.80,1.02,1.14,1.30,1.32,1.32,1.22,1.04,.93,.78,.73,
     8          .76,.80,1.02,1.14,1.31,1.33,1.34,1.23,1.05,.93,.77,.72,
     9          .75,.79,1.02,1.14,1.32,1.34,1.35,1.24,1.05,.93,.76,.71,
     *          .74,.78,1.02,1.15,1.33,1.36,1.37,1.25,1.06,.92,.76,.70/
      IF (KYR.EQ.1) WATER=FC
      DDBASE=5.56
C     DDAYS IS DRY-DAYS; DEGD, DEGREE-DAYS; RTOT, TOTAL RAIN FOR
C     YEAR; PETA, ANNUAL PET; TE, TEMPERATURE EFFICIENCY FOR PET;
C     CDAY, CALENDAR DAY; CUMDD, CUMULATIVE DRY-DAYS;
C     TOTRO, TOTAL RUNOFF (WATER-FC>0);
C     BGS, EGS ARE BEGINNING AND END OF GROWING SEASON.
      DDAYS=0.0
      DEGD=0.0
      RTOT=0.0
      PETA=0.0
      TE=0.0
      CDAY=15.0
      CUMDD=0.0
      TOTRO=0.0
      BGS=1.0
      IGS=0
      EGS=365.0
C     GENERATE WEATHER:
      DO 21 M=1,12
      T(M)=XT(M)+GAUSS(IDUM)*VT(M)
      IF (T(M).GT.DDBASE) DEGD=DEGD+(T(M)-DDBASE)*DAYS(M)
      IF (T(M).LE.0.0) T(M)=0.0
      TE=TE+(0.2*T(M))**1.514
      R(M)=XR(M)+GAUSS(IDUM)*VR(M)
      IF (R(M).LT.0.0) R(M)=0.0
      RTOT=RTOT+R(M)
C     FIND BGS, EGS, AND TOTAL GROWING SEASON LENGTH TGS
C     (IGS MAKES SURE THE SEASON'S FIRST WARMING IS BGS)
      LM=M-1
      IF (LM.EQ.0) LM=12
      IF (T(M).GE.DDBASE.AND.T(LM).GE.DDBASE) GO TO 21
      IF (T(LM).LT.DDBASE.AND.T(M).GE.DDBASE.AND.IGS.EQ.0) THEN
          BGS=TLINE(LM,T(LM),M,T(M),DDBASE)
          IGS=1
          ENDIF
      IF (T(LM).GE.DDBASE.AND.T(M).LT.DDBASE)
     2   EGS=TLINE(LM,T(LM),M,T(M),DDBASE)
   21 CONTINUE
      TGS=EGS-BGS+1
        IF (TGS.LE.0.0) TGS=TGS+365.0
C     TE AND A ARE USED IN THORNTHWAITE-MATHER PET EQUATION
      A=(0.675*TE**3-77.1*TE**2+17920.0*TE+492390.0)*0.000001
C     WEATHER TRACER:
      IF (MOD(KYR,ILOG).EQ.0)
     2  WRITE(11,1107) KYR, BGS, EGS, TGS
 1107 FORMAT(//5X,'SIMULATION YEAR:',I5,
     3  //10X,'BGS:',F7.1,' AND EGS:',F7.1,';  TGS:',F7.1,
     4  //10X,'MONTH, TEMP, PET, RAIN, WATER, RUNOFF, & DRY-DAYS:'/)
      DO 22 M=1,12
      OWATER=WATER
      RUNOFF=0.0
      PET=1.6*((10.0*T(M)/TE)**A)*CLAT(M,LAT)
      PETA=PETA+PET
C     SOIL WATER IS OLD WATER + RAIN - PET, TRUNCATED TO FC OR WP
      WATER=OWATER+R(M)-PET
      IF (WATER.LT.WP) WATER=WP
      IF (WATER.GT.FC) THEN
        RUNOFF=WATER-FC
        TOTRO=TOTRO+RUNOFF
        WATER=FC
        ENDIF
C     INTERPOLATE DROUGHT DAYS BETWEEN MONTHS AS NECESSARY
      OCDAY=CDAY
      CDAY=CDAY+DAYS(M)
C     NOT IN GROWING SEASON:
      IF (CDAY.LE.BGS.OR.OCDAY.GE.EGS) DDM=0.0
C     SUFFICIENT SOIL WATER:
      IF (OWATER.GT.WP.AND.WATER.GT.WP) DDM=0.0
C     DRY BOTH MONTHS:
      IF (OWATER.LE.WP.AND.WATER.LE.WP) THEN
        DDM=DAYS(M)
        IF (OCDAY.LT.BGS.AND.CDAY.GT.BGS) DDM=CDAY-BGS
        IF (OCDAY.LT.EGS.AND.CDAY.GT.EGS) DDM=EGS-OCDAY
        ENDIF
C     DRYING PERIOD:
      IF (OWATER.GT.WP.AND.WATER.LE.WP) THEN
        DDM=DAYS(M)*(WP-WATER)/(OWATER-WATER)
        IF (OCDAY.LT.BGS.AND.CDAY.GT.BGS) DDM=AMIN1(DDM,CDAY-BGS)
        IF (OCDAY.LT.EGS.AND.CDAY.GT.EGS) DDM=EGS-CDAY+DDM
        IF (DDM.LT.0.0) DDM=0.0
        ENDIF
C     WETTING PERIOD:
      IF (OWATER.LE.WP.AND.WATER.GT.WP) THEN
        DDM=DAYS(M)*(WP-OWATER)/(WATER-OWATER)
        IF (OCDAY.LT.BGS.AND.CDAY.GT.BGS) DDM=OCDAY+DDM-BGS
        IF (DDM.LT.0.0) DDM=0.0
        IF (OCDAY.LT.EGS.AND.CDAY.GT.EGS) DDM=AMIN1(DDM,EGS-OCDAY)
        ENDIF
      DDAYS=DDAYS+DDM
      CUMDD=CUMDD+DDM
C     ... MORE TRACER:
      IF (MOD(KYR,ILOG).EQ.0) WRITE(11,1109) M, T(M), PET, R(M),
     2  WATER, RUNOFF, DDM
 1109 FORMAT(10X,I3,2X,2F6.1,F7.1,F7.1,F5.1,F7.1)
   22 CONTINUE
C     RELATIVIZE DRY-DAYS TO GROWING SEASON LENGTH
      DDAYS=DDAYS/TGS
C     ... AND A BIT MORE TRACER:
      IF (MOD(KYR,ILOG).EQ.0) WRITE(11,1199)
     2  DEGD, PETA, RTOT, TOTRO, CUMDD, DDAYS
 1199 FORMAT(/10X,'DEGD:',F7.1,5X,'TOTAL PET:',F6.1,5X,'RAIN:',F6.1,
     2  /10X,'TOTAL RUNOFF:',F7.1,
     3  /10X,'SUMMED DRY-DAYS:',F7.1,3X,'RELATIVE DDI:',F7.3)
C     COMPUTE SOIL MOISTURE AND DEGD FACTORS:
      DO 23 KS=1,NSPP
C     SOIL MOISTURE FACTOR:
      SMF(KS)=DRTF(MDRT(KS),DDAYS)
C     TEMPERATURE:
      DDF(KS)=DEGDF(DDMIN(KS),DDMAX(KS),DEGD)
   23 CONTINUE
      RETURN
      END
C
C     BOOKS KEEPS TRACK OF BIOMASS, LEAF AREA, AND AVAILABLE LIGHT
C
      SUBROUTINE BOOKS
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
C     ALA, ACTUAL LEAF AREA; CLA, CUMULATIVE; BLA, ON-BEAM;
C     BMI, BIOMASS INCREMENT PER PLOT;  BA, BASAL AREA (FOR TRACER).
      REAL ALA(MP,MH), CLA(MP,MH), BLA(MP,MH), BMI(MP)
      REAL BA(MS)
C     C3, COMPENSATION POINT FOR 5 TOLERANCE CLASSES
      REAL C3(5)
C     XK, LIGHT EXTINCTION COEFFICIENT FOR BEER'S LAW
      REAL XK
      DATA XK/0.40/
      DATA C3/0.05,0.06,0.07,0.08,0.09/
C     REPORT PROGRESS TO SCREEN:
      IF (MOD(KYR,IPRT).EQ.0) WRITE(3,3000) KYR
 3000 FORMAT(5X,'SIMULATION YEAR:',I5)
C     INITIALIZE TOTAL BIOMASS AND DENSITY, AND BA, FOR TRACER
      SUMB=0.0
      SUMD=0.0
      DO 301 KS=1,NSPP
      BA(KS)=0.0
  301 CONTINUE
      DO 31 KP=1,NPLOTS
C
C     TO DUMP INITIALS FOR FUTURE RUNS, UN-COMMENT THESE WRITES:
C     IF (KYR.EQ.NYRS) WRITE(9,9003) KP, IND(KP)
C9003 FORMAT(I5,5X,I5)
C
      BMI(KP)=0.0
      BIOM(KP)=0.0
      CLAI(KP)=0.0
      DO 311 M=1,MH
      ALA(KP,M)=0.0
      CLA(KP,M)=0.0
      BLA(KP,M)=0.0
  311 CONTINUE
C     TALLY DENSITY OF STEMS
      SUMD=SUMD+FLOAT(IND(KP))
      IF (IND(KP).EQ.0) GO TO 33
      DO 32 KI=1,IND(KP)
C
C     ... MORE WRITES TO FILE INDATA
C     IF (KYR.EQ.NYRS) WRITE(9,9004) KI, ISP(KP,KI), DBH(KP,KI),
C    2  IBC(KP,KI)
C9004 FORMAT(2I5,F10.5,I5)
C
      KS=ISP(KP,KI)
      D=DBH(KP,KI)
C     BIOMASS (KG) AND BASAL AREA (SQ.M):
      BIOM(KP)=BIOM(KP)+WOOD(D)
      BA(KS)=BA(KS)+3.14159*(D/200.0)**2
C     POTENTIAL BIOMASS INCREMENT:
      TBI=WOOD(D+OPTINC(D,G(KS),B2(KS),B3(KS)))-WOOD(D)
      BMI(KP)=BMI(KP)+TBI
C     COMPUTE TREE HEIGHT AND HANG LEAVES:
      HT=HEIGHT(D,B2(KS),B3(KS))
      IHT=HT
      IF (IHT.LT.1) IHT=1
      IF (IHT.GT.MH) IHT=MH
C     HANG LEAVES IN CANOPY:
C     FD=UNIT FOLIAGE DENSITY
      FD=ALEAF(D)/FLOAT(IHT)
      DO 321 M=IHT,1,-1
C     STOP WHEN OUT OF LIGHT OR BELOW OLD CANOPY
      IF (AL(KP,M).LT.C3(LIGHT(KS)).OR.M.LT.IBC(KP,KI)) GO TO 322
      ALA(KP,M)=ALA(KP,M)+FD
  321 CONTINUE
  322 IBC(KP,KI)=M+1
   32 CONTINUE
C     COMPUTE AVAILABLE LIGHT FOR INDEPENDENT-PLOT MODE
   33 CLA(KP,MH)=ALA(KP,MH)
      AL(KP,MH)=1.0
      DO 331 M=(MH-1),1,-1
      CLA(KP,M)=ALA(KP,M)+CLA(KP,M+1)
      AL(KP,M)=EXP(-XK*CLA(KP,M+1)/AREA)
  331 CONTINUE
C     ACCRUE MEAN PLOT BIOMASS (T) FOR TRACER
      B=BIOM(KP)/1000.0
      SUMB=SUMB+B
C     CORRECT LEAF-AREA INDEX FOR PLOT SIZE:
      CLAI(KP)=CLA(KP,1)/AREA
C     COMPUTE SOIL FERTILITY CONSTRAINT IN TERMS OF MAX BMI ...
      IF (BMI(KP).GT.SF) THEN
        RF=SF/BMI(KP)
        ELSE
        RF=1.0
        ENDIF
C     ... AND COMPUTE FERTILITY FACTORS
      SFF(KP,1)=FERTF(1,RF)
      SFF(KP,2)=FERTF(2,RF)
      SFF(KP,3)=FERTF(3,RF)
   31 CONTINUE
C     TRACERS, CONVERTED TO PER-HA AND DUMPED TO UNIT 4:
      SS=FLOAT(NPLOTS)
      XB=(SUMB/SS)*(10000.0/AREA)
      XD=(SUMD/SS)*(10000.0/AREA)
      XBATOT=0.0
      DO 302 KS=1,NSPP
      BA(KS)=(BA(KS)/SS)*(10000.0/AREA)
      XBATOT=XBATOT+BA(KS)
  302 CONTINUE
      WRITE(4,4001) KYR, XD, XB, XBATOT, (BA(KS),KS=1,NSPP)
 4001 FORMAT(I4,2X,3F8.2/6X,10F6.2/6X,10F6.2/6X,10F6.2)
C     DUMP FOLIAGE PROFILE:
      IF (KYR.EQ.NYRS) THEN
        DO 303 M=MH,1,-1
        WRITE(10,1099) (ALA(KP,M)/AREA,KP=1,NPLOTS)
  303   CONTINUE
        ENDIF
 1099 FORMAT(10F7.2)
C     DUMP DIAGNOSTICS IN LOG YEAR, ELSE GET OUT IF MODE=0:
      IF (MODE.EQ.0) THEN
        IF (MOD(KYR,ILOG).EQ.0) GO TO 360
        RETURN
        ENDIF
C     DO AGGREGATE LEAF AREA PROFILE IF IN INTERACTIVE MODE:
  340 DO 34 KP=1,NPLOTS
      DO 35 M=1,(MH-1)
C     DIRECT-BEAM INSOLATION ACCRUES IN STEPS OF MXS LEFTWARD:
      BLA(KP,M)=DLA(KP,M+1,MXS,-1,ALA,NPLOTS)
      BLAI=BLA(KP,M)/AREA
      ALB=PHIB*EXP(-XK*BLAI)
C     DIFFUSE LIGHT IS TALLIED ACROSS A SKY ARC IN 7 SAMPLES ...
      ALD=0.0
C     ... AS 3 DLA'S TO EACH SIDE, WITH HEIGHTS 2, 6, AND 10 M ...
      DO 351 KM=2,10,4
      XLL=DLA(KP,M+1,KM,-1,ALA,NPLOTS)/AREA
      ALD=ALD+(PHID/7.0)*EXP(-XK*XLL)
      XLR=DLA(KP,M+1,KM,1,ALA,NPLOTS)/AREA
      ALD=ALD+(PHID/7.0)*EXP(-XK*XLR)
  351 CONTINUE
C     ... + THE VERTICAL PROFILE:
      XVL=CLA(KP,M+1)/AREA
      ALD=ALD+(PHID/7.0)*EXP(-XK*XVL)
C     TOTAL LIGHT IS DIRECT-BEAM + DIFFUSE:
      AL(KP,M)=ALB+ALD
   35 CONTINUE
   34 CONTINUE
C     DUMP LEAF AREA PROFILES IF LOG YEAR
      IF (MOD(KYR,ILOG).NE.0) RETURN
  360 WRITE(11,1101)
 1101 FORMAT(/5X,'ACTUAL LEAF-AREA INDEX PROFILE FOR 8 CELLS:'/)
      DO 36 M=MH,1,-1
      WRITE(11,1102) M, (ALA(KP,M)/AREA,KP=1,8)
 1102 FORMAT(5X,I5,8F7.2)
   36 CONTINUE
      WRITE(11,1103)
 1103 FORMAT(/5X,'CUMULATIVE LEAF AREA PROFILE:'/)
      DO 37 M=MH,1,-1
      WRITE(11,1102) M, (CLA(KP,M)/AREA,KP=1,8)
   37 CONTINUE
      IF (MODE.EQ.1) THEN
        WRITE(11,1104) MXS
 1104   FORMAT(/5X,'DIAGONAL LEAF AREA PROFILE (MXS=',I2,'):'/)
        DO 38 M=MH,1,-1
        WRITE(11,1102) M, (BLA(KP,M)/AREA,KP=1,8)
   38   CONTINUE
        ENDIF
      WRITE(11,1105)
 1105 FORMAT(/5X,'AVAILABLE LIGHT PROFILE:'/)
      DO 39 M=MH,1,-1
      WRITE(11,1102) M, (AL(KP,M),KP=1,8)
   39 CONTINUE
      RETURN
      END
C
C     MORTAL KILLS TREES
C
      SUBROUTINE MORTAL
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      DIMENSION LIVE(7), NDEAD1(7), NDEAD2(7)
      LOGICAL AMORT, SMORT
      IF (MOD(KYR,ILOG).EQ.0) WRITE(11,1112)
 1112 FORMAT(//5X,'MORTALITY BY 10-CM SIZE CLASSES:')
      DO 41 KP=1,NPLOTS
C     IDEAD AND NDEAD TALLY DEAD TREES FOR DIAGNOSTICS
      IDEAD=0
      DO 411 KC=1,7
      LIVE(KC)=0
      NDEAD1(KC)=0
      NDEAD2(KC)=0
  411 CONTINUE
      IF (IND(KP).EQ.0) GO TO 41
      DO 42 KI=1,IND(KP)
      KS=ISP(KP,KI)
      D=DBH(KP,KI)
C     TALLY STEMS IN 10-CM DBH CLASSES (FOR DIAGNOSTICS)
      KC=(D/10.0)+1
      IF (KC.GT.7) KC=7
      LIVE(KC)=LIVE(KC)+1
C     MORTALITY CAN COME FROM STRESS, OR AMBIENT ...
C     UNHEALTHY TREE DIES:
      IF (SMORT(NOGRO(KP,KI),IDUM)) THEN
C       LABEL DBH AND TALLY 1 DEAD FROM STRESS
        DBH(KP,KI)=-2.0
        NDEAD2(KC)=NDEAD2(KC)+1
        IDEAD=IDEAD+1
        GO TO 42
        ENDIF
C     HEALTHY TREE DIES:
      IF (AMORT(AGEMAX(KS),IDUM)) THEN
C       LABEL DBH AND TALLY 1 DEAD OF NATURAL CAUSES
        DBH(KP,KI)=-1.0
        NDEAD1(KC)=NDEAD1(KC)+1
        IDEAD=IDEAD+1
        ENDIF
   42 CONTINUE
C     TRACER ON MORTALITY:
      IF (KP.EQ.6.AND.MOD(KYR,ILOG).EQ.0) WRITE(11,1113) KP,
     2  (LIVE(KC),KC=1,7), (NDEAD1(KC),KC=1,7),
     3  (NDEAD2(KC),KC=1,7)
 1113 FORMAT(/10X,'PLOT:',I3,//10X,'ALIVE:',7I4,/10X,'NDEAD:',7I4,
     2  /10X,'SDEAD:',7I4)
C     RESHUFFLE TREE ARRAYS
      IN=0
      IF (IDEAD.EQ.0) GO TO 41
      DO 43 IO=1,IND(KP)
      IF (DBH(KP,IO).EQ.0.0) GO TO 44
      IF (DBH(KP,IO).LT.0.0) GO TO 43
      IN=IN+1
      DBH(KP,IN)=DBH(KP,IO)
      NOGRO(KP,IN)=NOGRO(KP,IO)
      ISP(KP,IN)=ISP(KP,IO)
      IBC(KP,IN)=IBC(KP,IO)
   43 CONTINUE
   44 IND(KP)=IN
      NXS=IN+1
      DO 45 IX=NXS,NMAX
      DBH(KP,IX)=0.0
      NOGRO(KP,IX)=0
      ISP(KP,IX)=0
      IBC(KP,IX)=0
   45 CONTINUE
   41 CONTINUE
      RETURN
      END
C
C     REGEN SPROUTS STUMPS AND PLANTS SAPLINGS
C
      SUBROUTINE REGEN
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      INTEGER NUM(MS)
      REAL PROBS(MS), XS(MS)
      DO 51 KP=1,NPLOTS
C     HEADER FOR DIAGNOSTICS:
      IF (KP.EQ.6.AND.MOD(KYR,ILOG).EQ.0) WRITE(11,1110) KP
 1110 FORMAT(/5X,'REGEN IN PLOT:',I5/)
      IF (IND(KP).EQ.NMAX) GO TO 51
C     ZERO TALLIES FOR SAPLINGS:
C     XSTOT IS EXPECTED TOTAL SAPLINGS;
C     NS IS NUMBER OF SAPLINGS TO PLANT, NPOSS THE NUMBER POSSIBLE.
      XSTOT=0.0
      NS=0
      NPOSS=0
      DO 52 KS=1,NSPP
C     AVAILABLE LIGHT REGEN FACTOR:
      ALRF=ALF(LIGHT(KS),AL(KP,1))
C     REGEN FACTOR IS LIGHT*MIN(WATER,NUTRIENTS)*TEMPERATURE
      RF=ALRF*AMIN1(SMF(KS),SFF(KP,NUTRI(KS)))*DDF(KS)
      XS(KS)=SEED(KS)*RF
      XSTOT=XSTOT+XS(KS)
C     ZERO NUMBER-TO-PLANT
      NUM(KS)=0
   52 CONTINUE
C     NUMBER OF SAPLINGS IS MIN(AVAILABLE,POSSIBLE):
      NSTOT=XSTOT+RAN2(IDUM)
      NPOSS=NMAX-IND(KP)
      NS=MIN(NSTOT,NPOSS)
C     IF NO SEEDLINGS, SKIP ON OUT:
      IF (NS.EQ.0) GO TO 51
C     ELSE, SELECT SAPLINGS BASED ON RELATIVE RATES...
C     COMPUTE RELATIVE PROBABILITY ARRAY:
   53 DO 531 KS=1,NSPP
      PROBS(KS)=XS(KS)/XSTOT
      IF (KS.GE.2) PROBS(KS)=PROBS(KS)+PROBS(KS-1)
  531 CONTINUE
C     SELECT SAPLINGS PROBABILISTICALLY (INS IS INSEEDING):
      INS=0
  532 Y=RAN2(IDUM)
      KS=1
      IF (Y.LT.PROBS(KS)) GO TO 535
  533 DO 534 KS=2,NSPP
      IF (Y.GE.PROBS(KS-1).AND.Y.LT.PROBS(KS)) GO TO 535
  534 CONTINUE
C     ADD A SAPLING OF THE SELECTED SPECIES:
  535 NUM(KS)=NUM(KS)+1
      INS=INS+1
C     KEEP GOING UNTIL ENOUGH (NS) HAVE BEEN SELECTED
      IF (INS.LT.NS) GO TO 532
C     REGEN DIAGNOSTICS:
      IF (KP.EQ.6.AND.MOD(KYR,ILOG).EQ.0) WRITE(11,1111)
     2  NPOSS, NSTOT, NS
 1111 FORMAT(/10X,'NPOSS:',I4,3X,'NSTOT:',I4,3X,'NS:',I4,
     2  //10X,'KS, MSP, XS, NUM:'/)
C     PLANT TREES AND ASSIGN INITIAL DIAMETERS
      IN=IND(KP)
      DO 54 KS=1,NSPP
C     TRACER:
      IF (KP.EQ.6.AND.MOD(KYR,ILOG).EQ.0) WRITE(11,1112)
     2  KS, MSP(KS), XS(KS), NUM(KS)
 1112 FORMAT(10X,I2,1X,A4,1X,F6.2,I5)
      IF (NUM(KS).EQ.0) GO TO 54
      DO 55 KI=1,NUM(KS)
      IN=IN+1
      DBH(KP,IN)=2.50+0.25*GAUSS(IDUM)
      ISP(KP,IN)=KS
      IBC(KP,IN)=1
      NOGRO(KP,IN)=0
   55 CONTINUE
   54 CONTINUE
C     UPDATE TALLY OF INDIVIDUALS
      IND(KP)=IN
   51 CONTINUE
      RETURN
      END
C
C     GROW INCREMENTS TREE DIAMETERS
C
      SUBROUTINE GROW
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      DO 61 KP=1,NPLOTS
C     HEADER FOR GROWTH-FACTOR DIAGNOSTICS
      IF (MOD(KYR,ILOG).EQ.0.AND.KP.EQ.6)
     2  WRITE(11,1114) KP, IND(KP)
 1114 FORMAT(/5X,'GROWTH FACTOR TRACE:',//10X,'PLOT:',I3,5X,
     2  'NUMBER OF TREES:',I5,//10X,
     3  'KI, KSP, DBH, HT, BC, ALF, SMF, SFF, DDF, GF, DINC, NOGRO'/)
      IF (IND(KP).EQ.0) GO TO 61
      DO 62 KI=1,IND(KP)
      KS=ISP(KP,KI)
      D=DBH(KP,KI)
C     COMPUTE TREE HEIGHT:
      HT=HEIGHT(D,B2(KS),B3(KS))
      IHT=HT
      IF (IHT.LT.1) IHT=1
      IF (IHT.GT.MH) IHT=MH
C     AVAILABLE LIGHT FACTOR:
      ALGF=0.0
      IF (IBC(KP,KI).GT.IHT) GO TO 622
  620 DO 621 M=IBC(KP,KI),IHT
      ALGF=ALGF+ALF(LIGHT(KS),AL(KP,M))
  621 CONTINUE
      ALGF=ALGF/FLOAT(IHT-IBC(KP,KI)+1)
  622 CONTINUE
C     GROWTH FACTOR IS LIGHT*MIN(WATER,NUTRIENTS)*TEMPERATURE:
      GF=ALGF*AMIN1(SMF(KS),SFF(KP,NUTRI(KS)))*DDF(KS)
C     POTENTIAL GROWTH INCREMENT
      DM=OPTINC(D,G(KS),B2(KS),B3(KS))
C     REALIZED INCREMENT
      DINC=DM*GF
      DBH(KP,KI)=DBH(KP,KI)+DINC
C     MEAN DBH INCREMENT FOR THIS SPECIES:
      XINC=DMAX(KS)/AGEMAX(KS)
C     NOGRO, IF GROWTH<10% OF OPT, OR DINC<10% AVERAGE
      IF (GF.LT.0.10.OR.DINC.LT.(0.1*XINC)) THEN
        NOGRO(KP,KI)=NOGRO(KP,KI)-1
        ELSE
        NOGRO(KP,KI)=0
        ENDIF
C     DIAGNOSTICS:
      IF (KP.EQ.6.AND.MOD(KYR,ILOG).EQ.0) WRITE(11,1116) KI, KS, D,
     2  IHT, IBC(KP,KI), ALGF, SMF(KS), SFF(KP,NUTRI(KS)), DDF(KS),
     3  GF, DINC, NOGRO(KP,KI)
 1116 FORMAT(10X,2I3,F7.2,2I3,1X,4F5.2,F6.2,F7.2,I3)
C     DUMP GROWTH DIAGNOSTICS FOR EXTERNAL STATS:
C     IF (KYR.EQ.50.OR.KYR.EQ.100.OR.KYR.EQ.NYRS)
C    2  WRITE(10,1016) KYR, KP, KI, KS, D, HT,
C    3    ALGF, SMF(KS), SFF(KP,NUTRI(KS)), DDF(KS), GF
c1016 FORMAT(I4,I5,1X,2I3,F7.2,F6.1,2X,4F5.2,F6.2)
   62 CONTINUE
   61 CONTINUE
      RETURN
      END
C
C     PRINT WRITES AGGREGATE STAND STATISTICS FOR ALL PLOTS
C
      SUBROUTINE PRINT
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      REAL SBA(MS), RIV(MS), RDC(10)
      INTEGER ISDC(MS,10), IDC(10), ISD(MS)
      XAREA=FLOAT(NPLOTS)*AREA/10000.0
      WRITE(8,8009) KYR
 8009 FORMAT(//5X,'SIMULATION YEAR:',I5)
C     ITD IS (INTEGER) TOTAL DENSITY (ITDA, STEMS>2.5 CM);
C     ISD, SBA, RIV ARE SPECIES DENSITY, BASAL AREA, AND IMPORTANCE;
C     IDC IS STEM TALLY BY 10-CM DBH CLASS (RDC, PER HA); ISDC, BY SPP.
      ITD=0
      ITDA=0
      TBA=0.0
      TBIOM=0.0
      TLA=0.0
      TDBH=0.0
      XDBH=0.0
      DO 701 KS=1,NSPP
      ISD(KS)=0
      SBA(KS)=0.0
      DO 702 KC=1,10
      ISDC(KS,KC)=0
  702 CONTINUE
  701 CONTINUE
      DO 703 KC=1,10
      IDC(KC)=0
  703 CONTINUE
      DO 71 KP=1,NPLOTS
      IF (IND(KP).EQ.0) GO TO 71
      DO 72 KI=1,IND(KP)
      KS=ISP(KP,KI)
      D=DBH(KP,KI)
      ITDA=ITDA+1
      IF (D.LT.2.5) GO TO 72
C     INCREMENT DIAMETER CLASSES
      KC=(D/10.0)+1
      IF (KC.GT.10) KC=10
      ISDC(KS,KC)=ISDC(KS,KC)+1
      IDC(KC)=IDC(KC)+1
C     INCREMENT STEM DENSITY
      ISD(KS)=ISD(KS)+1
      ITD=ITD+1
C     INCREMENT BASAL AREA
      TDBH=TDBH+D
      DM=D/100.0
      SBA(KS)=SBA(KS)+3.14159*(DM/2.0)**2
      TBA=TBA+3.14159*(DM/2.0)**2
   72 CONTINUE
      TBIOM=TBIOM+BIOM(KP)
      TLA=TLA+CLAI(KP)
   71 CONTINUE
      SUMRIV=0.0
      WRITE(8,8010)
 8010 FORMAT(//5X,'AGGREGATE STATISTICS:',
     2  //5X,'SPP:     DBH DISTRIBUTION (10-CM CLASSES),',
     3  ' DENSITY, B.A., RIV200'/)
      DO 73 KS=1,NSPP
C     CHECK TO MAKE SURE THERE ARE TREES ...
      IF (ITD.GT.0) THEN
        RSD=100.0*FLOAT(ISD(KS))/FLOAT(ITD)
        RBA=100.0*SBA(KS)/TBA
        RIV(KS)=(RSD+RBA)/2.0
      ELSE
        RSD=0.0
        RBA=0.0
        RIV(KS)=0.0
        ENDIF
C     SKIP SPECIES IF UNREPRESENTED
      IF (ISDC(KS,1).EQ.0.AND.RIV(KS).EQ.0.0) GO TO 73
C     CONVERT SPP DBH CLASSES TO PER-HA AND ROUND TO INTEGERS:
      DO 731 L=1,10
      ISDC(KS,L)=INT(FLOAT(ISDC(KS,L))/XAREA +0.5)
  731 CONTINUE
C     CONVERT OTHER TALLIES TO PER-HA
      SD=FLOAT(ISD(KS))/XAREA
      SBA(KS)=SBA(KS)/XAREA
      WRITE(8,8011) MSP(KS), (ISDC(KS,L),L=1,10), SD,
     2  SBA(KS), RIV(KS)
 8011 FORMAT(5X,A4,I6,9I4,F7.1,1X,F6.2,2X,F6.2)
      SUMRIV=SUMRIV+RIV(KS)
   73 CONTINUE
      DO 74 L=1,10
      RDC(L)=FLOAT(IDC(L))/XAREA
   74 CONTINUE
      TBIOM=TBIOM/1000.0
      TBMHA=TBIOM/XAREA
      TDHA=FLOAT(ITD)/XAREA
      TDA=FLOAT(ITDA)/XAREA
      TBAHA=TBA/XAREA
      XLAI=TLA/FLOAT(NPLOTS)
      IF (ITD.GT.0) XDBH=TDBH/FLOAT(ITD)
C     DUMP STAND SUMMARY:
      WRITE(8,8012) (RDC(L),L=1,10), TDHA, TDA, TBAHA, XDBH, TBMHA, XLAI
 8012 FORMAT(/5X,'STAND ESTIMATES:',
     2  //10X,'DBH DISTRIBUTION (TREES/HA):',/10X,F8.1,F7.1,8F6.1,
     3  /10X,'DENSITY (>2.5 CM):',F8.2,'/HA',2X,'ALL STEMS:',F8.1,'/HA',
     4  /10X,'BASAL AREA:',F8.3,' SQ.M/HA',5X,'MEAN DBH:',F7.2,
     5  ' CM',/10X,'TOTAL WOODY BIOMASS:',F9.3,
     6  ' MG/HA',5X,'LEAF-AREA INDEX:',F7.3)
      RETURN
      END
C
C     PUNCH MIMICS SR PRINT, BUT PUNCHES CELL SUMMARIES TO UNIT 7
C
      SUBROUTINE PUNCH
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      COMMON/CONTRL/ MODE, INDATA, NPLOTS, NSPP, NYRS,
     2       KYR, IPRT, ILOG, IPCH, NMAX, IDUM, MXS
      COMMON/SPP/ MSP(MS), AGEMAX(MS), DMAX(MS), HTMAX(MS),
     2       G(MS), B2(MS), B3(MS), DDMIN(MS), DDMAX(MS), DDF(MS),
     3       LIGHT(MS), NUTRI(MS), MDRT(MS), SMF(MS), SEED(MS)
      COMMON/TREES/ ISP(MP,MT), DBH(MP,MT), IBC(MP,MT), NOGRO(MP,MT)
      COMMON/PLOTS/ IND(MP), BIOM(MP), SFF(MP,3),
     2       CLAI(MP), AL(MP,MH)
      COMMON/SITE/ AREA, LAT, THETA, PHIB, PHID, BGS, EGS, TGS,
     2       FC, WP, SF, XT(12), VT(12), XR(12), VR(12)
      REAL SBA(MS), DC(10)
      INTEGER IDC(10)
      XAREA=AREA/10000.0
      DO 81 KP=1,NPLOTS
      ITD=0
      TBA=0.0
      TDHA=0.0
      TBAHA=0.0
      TBMHA=0.0
      DO 811 KS=1,NSPP
      SBA(KS)=0.0
  811 CONTINUE
      DO 812 L=1,10
      IDC(L)=0
      DC(L)=0.0
  812 CONTINUE
      IF (IND(KP).EQ.0) GO TO 85
      DO 82 KI=1,IND(KP)
      KS=ISP(KP,KI)
      D=DBH(KP,KI)
C     INCREMENT DIAMETER CLASSES
      KC=(D/10.0)+1
      IF (KC.GT.10) KC=10
      IDC(KC)=IDC(KC)+1
C     INCREMENT STEM DENSITY
      ITD=ITD+1
C     INCREMENT BASAL AREA
      DM=D/100.0
      SBA(KS)=SBA(KS)+3.14159*(DM/2.0)**2
      TBA=TBA+3.14159*(DM/2.0)**2
   82 CONTINUE
C     CONVERT TO PER-HA
      TBIOM=BIOM(KP)/1000.0
      TBMHA=TBIOM/XAREA
      TDHA=FLOAT(ITD)/XAREA
      TBAHA=TBA/XAREA
      DO 83 L=1,10
      DC(L)=FLOAT(IDC(L))/XAREA
   83 CONTINUE
      DO 84 KS=1,NSPP
      SBA(KS)=SBA(KS)/XAREA
   84 CONTINUE
   85 WRITE(7,7001) KYR, KP, TDHA, TBAHA, TBMHA, CLAI(KP),
     2  (DC(L),L=1,10), (SBA(KS),KS=1,NSPP)
 7001 FORMAT(2I4,2X,F8.2,2F7.2,F6.2,/5X,F7.1,9F6.1,/5X,10F6.2,
     2  /5X,10F6.2/5X,10F6.2)
   81 CONTINUE
      RETURN
      END
C
C     FUNCTIONS:
C
C     DLA RETURNS A DIAGONAL LEAF-AREA PROFILE WITH VERTICAL
C     STEP-SIZE MXS, IN DIRECTION IDIR (-1=LEFT, 1=RIGHT)
C     (CALLED ONLY IF ZELIG IN IS INTERACTIVE-TRANSECT MODE)
C
      FUNCTION DLA(KPLT,MHT,LXS,IDIR,ALA,NPLOTS)
      PARAMETER (MP=100,MT=1000,MS=30,MH=100)
      REAL ALA(MP,MH)
      INTEGER KPLT, MHT, LXS, IDIR, NPLOTS
C     PL=REL. PATH LENGTH OF LIGHT (RATIO OF HEIGHT TO HYPOTENUSE)
      PL=FLOAT(LXS)/SQRT(LXS**2+100.0)
      DLA=0.0
      KK=KPLT
      M1=MHT
      M2=M1+LXS
      DO 101 MM=M1,M2-1
      IF (MM.GT.MH) GO TO 105
  100 DLA=DLA+ALA(KK,MM)/PL
  101 CONTINUE
  102 IF (IDIR.GT.0) THEN
        KK=KK+1
        ELSE
        KK=KK-1
        ENDIF
      IF (KK.GT.NPLOTS) KK=KK-NPLOTS
      IF (KK.LT.1) KK=KK+NPLOTS
      M1=M2
      M2=M1+LXS
      DO 104 MM=M1,M2-1
      IF (MM.GT.MH) GO TO 105
  103 DLA=DLA+ALA(KK,MM)/PL
  104 CONTINUE
      GO TO 102
  105 RETURN
      END
C
C     RANDOM NUMBER GENERATORS:
C
C     RAN2(IDUM) RETURNS A UNIFORM RANDOM NUMBER ON (0,1).
C     THE ROUTINE IS FROM:  PRESS, W.H., B.P. FLANNERY, S.A.
C     TEUKOLSKY, AND W.T. VETTERLING.  1987.  NUMERICAL RECIPES.
C     CAMBRIDGE UNIVERSITY PRESS, CAMBRIDGE.
C
      FUNCTION RAN2(IDUM)
      DIMENSION IR(97)
      DATA M/714025/
      DATA IA/1366/
      DATA IC/150889/
      DATA IFF/0/
      RM=1.0/FLOAT(M)
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
      IFF=1
      IDUM=MOD(IC-IDUM,M)
      DO 11 J=1,97
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
   11 CONTINUE
      IDUM=MOD(IA*IDUM+IC,M)
      IY=IDUM
      ENDIF
      J=1+(37*IY)/M
      IF (J.GT.97.OR.J.LT.1) PAUSE
      IY=IR(J)
      RAN2=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
C
C     GAUSS(IDUM) RETURNS A NORMAL RANDOM NUMBER, X,S=(0,1).
C     THIS IS ROUTINE GASDEV(IDUM) FROM PRESS ET AL., ABOVE.
C
      FUNCTION GAUSS(IDUM)
      DATA ISET/0/
      IF (ISET.EQ.0) THEN
   21 V1=2.0*RAN2(IDUM)-1
      V2=2.0*RAN2(IDUM)-1
      R=V1**2+V2**2
      IF (R.GE.1.) GO TO 21
      FAC=SQRT(-2.*ALOG(R)/R)
      GSET=V1*FAC
      GAUSS=V2*FAC
      ISET=1
      ELSE
      GAUSS=GSET
      ISET=0
      ENDIF
      RETURN
      END
C
C     ALLOMETRIES:
C
C     FUNCTION HEIGHT (M), ALLOMETRIC ON DBH (CM), PER SPECIES
C
      REAL FUNCTION HEIGHT(DBH,B2,B3)
      REAL DBH, B2, B3
      HEIGHT=(137.0+B2*DBH-B3*DBH**2)/100.0
      RETURN
      END
C
C     LEAF AREA (SQ.M), ALLOMETRIC ON DBH (CM)
C
      REAL FUNCTION ALEAF(DBH)
      REAL DBH
      ALEAF=0.160694*DBH**2.129
      RETURN
      END
C
C     ABOVE-GROUND WOODY BIOMASS (KG), ON DBH (CM)
C
      REAL FUNCTION WOOD(DBH)
      REAL DBH
      WOOD=0.1193*DBH**2.393
      RETURN
      END
C
C     POTENTIAL GROWTH AS OPTIMAL DIAMETER INCREMENT:
C
      REAL FUNCTION OPTINC(D,G,B2,B3)
      REAL D, G, B2, B3
C     GROWTH EQUATION FROM JABOWA:
      GR=(137.0+0.25*B2**2/B3)*(0.5*B2/B3)
      DM=G*D*(1.0-(137.0*D+B2*D**2-B3*D**3)/GR)/
     2  (274.0+3.0*B2*D-4.0*B3*D**2)
      OPTINC=DM
      RETURN
      END
C
C     SPECIES RESPONSE FUNCTIONS:
C
C     AVAILABLE LIGHT FACTOR, BY TOLERANCE CLASS AND LIGHT:
C
      REAL FUNCTION ALF(KT,AL)
      REAL C1(5), C2(5), C3(5)
      DATA C1/1.01,1.04,1.11,1.24,1.49/
      DATA C2/4.62,3.44,2.52,1.78,1.23/
      DATA C3/0.05,0.06,0.07,0.08,0.09/
      ALF=C1(KT)*(1.0-EXP(-C2(KT)*(AL-C3(KT))))
      IF (ALF.LT.0.0) ALF=0.0
      RETURN
      END
C
C     FUNCTION DRTF COMPUTES A DRY-DAY CONSTRAINT FACTOR
C
      REAL FUNCTION DRTF(MDT,DDAYS)
      REAL DT, DDAYS, DRT
      INTEGER MDT
      DT=FLOAT(MDT)/10.0
      DRT=AMIN1(DT,DDAYS)
      DRTF=SQRT((DT-DRT)/DT)
      RETURN
      END
C
C     FUNCTION FERTF COMPUTES QUADRATIC NUTRIENT RESPONSE FACTORS
C     BY NUTRIENT RESPONSE CLASS (1=INTOL, 3=TOL) AND SOIL FERTILITY
C
      REAL FUNCTION FERTF(NRC,SF)
      INTEGER NRC
      REAL C1(3), C2(3), C3(3), SF
      DATA C1 /-0.6274,-0.2352,0.2133/
      DATA C2 /3.600,2.770,1.789/
      DATA C3 /-1.994,-1.550,-1.014/
      FERTF=C1(NRC)+C2(NRC)*SF+C3(NRC)*SF**2
      IF (FERTF.LT.0.0) FERTF=0.0
      RETURN
      END
C
C     FUNCTION DEGDF COMPUTES THE DEGREE-DAY PARABOLIC FACTOR
C     ON MIN, MAX, AND ACTUAL DEGREE-DAYS
C
      REAL FUNCTION DEGDF(DDMIN,DDMAX,DEGD)
      REAL DDMIN, DDMAX, DEGD
      DEGDF=4.0*(DEGD-DDMIN)*(DDMAX-DEGD)/(DDMAX-DDMIN)**2
      IF (DEGDF.LT.0.0) DEGDF=0.0
      RETURN
      END
C
C     MORTALITY FUNCTIONS ARE TRUE IF A TREE DIES
C
C     FUNCTION SMORT IS FOR STRESS MORTALITY
C
      LOGICAL FUNCTION SMORT(NOGRO,IDUM)
      INTEGER NOGRO, IDUM
      IF (NOGRO.LT.-1.AND.RAN2(IDUM).LT.(0.369)) THEN
        SMORT=.TRUE.
        ELSE
        SMORT=.FALSE.
        ENDIF
      RETURN
      END
C
C     FUNCTION AMORT IS FOR AMBIENT (NATURAL) MORTALITY
C
      LOGICAL FUNCTION AMORT(AGEMAX,IDUM)
      INTEGER IDUM
      REAL AGEMAX
      IF (RAN2(IDUM).LT.(4.605/AGEMAX)) THEN
        AMORT=.TRUE.
        ELSE
        AMORT=.FALSE.
        ENDIF
      RETURN
      END
C
C     FUNCTION TLINE DOES LINEAR INTERPOLATIONS ON TEMPERATURES,
C     TO FIND BGS AND EGS, GIVEN MONTHLY TEMPS AND THE GDD BASE
C
      REAL FUNCTION TLINE(M1,T1,M2,T2,DDB)
      REAL D1, T1, D2, T2, DDB
      INTEGER M1, M2
      REAL CD(12)
      DATA CD/15.0,46.0,74.0,105.0,135.0,166.0,196.0,227.0,258.0,
     2  288.0,319.0,349.0/
      D1=CD(M1)
      D2=CD(M2)
        IF (M1.EQ.12) D1=D1-365.0
C     REGRESSION SLOPE IS A; INTERCEPT IS B:
      A=(T2-T1)/(D2-D1)
      B=T2-A*D2
C     SOLVE FOR CALENDAR DAY AT WHICH T=DDBASE:
      D=(DDB-B)/A
      TLINE=D
      RETURN
      END

.