      PROGRAM MTCLIM
CM  MTCLIM
CM  VERSION 1.0 6-1-89
CM  COPYRIGHT 1987, 1988, 1989, 1990  Steven W. Running and Joseph C. Coughlan
CD  THIS PROGRAM PREDICTS MICROCLIMATE CONDITIONS FOR ANY SITE
CD  ON MOUNTAINOUS TERRAIN GIVEN BASE STATION METEOROLOGICAL DATA,
CD  SITE AND BASE STATION CHARACTERISTICS. -RRN
CD
CD  ORIGINAL CODE WRITTEN BY R.R.NEMANI
CD  REWRITTEN BY J.C.COUGHLAN ON 4-1-89, 6-1-89
CD  MICROSOFT FORTRAN VERSION 4.1 IN STANDARD, TRANSPORTABLE FORTRAN
CD
CD  IF MTCLIM RUNS SLOW, TOO SLOW, THEN SEE SUBROUTINE RAD FOR
CD  SUGGESTIONS ON OPTIMIZING THE MODEL FOR SPEED.  IT CAN BE OPTIMIZED
CD  TO RUN IN HALF THE TIME IT NOW TAKES WITH A FEW CHANGES TO 1 DO
CD  LOOP.  SEE SUBROUTINE RAD.  -JCC

CVF MTCLIM FILES/VARIABLES
CVF  'IFILE':FILE CONTAINING BASE AND SITE DESCRIPTIONS - SEE 'README.DOC'
CVF          CODE BELOW DESCRIBES SOME OF THE FILES FORMAT.
CVF   BFILE: INPUT- BASE STATION CLIMATE DATA
CVF          FILE STRUCTURE: 1 LINE PER RECORD IN FREE FORMAT
CVF          THE ORDER OF VARIABLES IS FOUND IN SUBROUTINE BREAD.
CVF   SFILE: OUTPUT- SIMULATED SITE CLIMATE FILE
CVF          WRITTEN IN FORMATTED OUTPUT AT END OF PROGRAM MTCLIM
CVF   INFILE: INTEGER FILE UNIT NUMBER CONNECTED TO BFILE
CVF
CVI SITE VARIABLES:
CVI   SLAT :  SITE LATITUDE
CVI   SELEV : SITE ELEVATION IN METERS
CVI   SASPCT: SITE ASPECT DEGREES
CVI   SSLOPE: SITE SLOPE  %
CVI   SLAI:   SITE LEAF AREA INDEX (ALL SIDES)
CVI   SISOH:  SITE ISOHYET
CVI   SALBDO: SITE ALBEDO %
CVI   SEHORZ: EAST-HORIZON
CVI   SWHORZ: WEST HORIZON
CVI   SFILE:  FILE NAME FOR OUTPUT OF MICROCLIMATE OF THE SITE
CVI
CVI BASE VARIABLES:
CVI   BELEV: BASE ELEVATION IN METERS
CVI   BISOH: BASE ISOHYET
CVI   BMAX:  BASE MAX TEMPERATURE, CELSIUS
CVI   BMIN:  BASE MIN TEMPERATURE, CELSIUS
CVI   BPPT:  BASE PRECIPITATION OF STATION(S), 1 OR 2 ARE SUPPORTED
CVI   BFILE: NAME OF THE FILE CONTAINING BASE STATION MET DATA
CVI
CVI PARAMETERS:
CVI   TLAPSE: LAPSE RATE FOR AIR TEMPERATURE
CVI   DLAPSE: LAPSE RATE FOR DEW POINT TEMPERATURE
CVI   MAXLAP: LAPSE RATE FOR MAX TEMPERATURE
CVI   MINLAP: LAPSE RATE FOR MIN TEMPERATURE
CVI   NDAYS:  # OF DAYS TO BE SIMULATED
CVI   MLAI:   MAXIMUM LEAF AREA INDEX, ASSUMED TO BE 10
CVI
CVO OUTPUT VARIABLES:
CVO   SRAD : SITE INCIDENT RADIATION
CVO   STEMP: DAYLIGHT AVERAGE AIR TEMPERATURE FOR THE SITE
CVO   SMIN:  NIGHT MINIMUM TEMPERATURE FOR THE SITE
CVO   SMAX:  MAXIMUM TEMPERATURE FOR THE SITE
CVO   SHUMD: DAYLIGHT AVERAGE RELATIVE HUMIDITY FOR THE SITE
CVO   SPPT:  PRECIPITATION AT THE SITE, CM


      REAL BMAX,BMIN,BDEW,SELEV,SLAT,SASPCT,SSLOPE
      REAL BISOH(2),BPPT(2),TRAN(365)
      REAL SISOH,SEHORZ,SWHORZ,SLAI,MLAI,SALBDO
      REAL SARAD,FARAD,TLAPSE,DLAPSE,SRAD,FRAD,MAXLAP,MINLAP
      REAL BELEV, TRANCF,TEMPCF,SMIN,SMAX,SHUMD,SPPT,STEMP
      REAL X, C2F, F2C, F2M
      CHARACTER*1 ANS0,ANS1,ANS2,ANS3,ANS4
      CHARACTER *12 SFILE,BFILE,IFILE
      INTEGER INFILE
      INTEGER NDAYS,NPPT,J,JDAY,FLAG
      LOGICAL USEJD, USEENG, USEDEW, USETOT, USETHR
C     MAX LAI IS SET TO 10.0
      DATA MLAI/10.0/

C  OPTIONAL CONVERSION FUNCTIONS FOR ENGLISH UNIT INPUT FILES ONLY!
C       TO CONVERT INPUT FILES FROM ENGLISH UNITS INTO METRIC
C       CONVERSIONS ALSO NEEDED IN SUBROUTINE TRANS.
C       OUTPUT CAN BE CONVERTED TO ENGLISH UNITS VIA 'INIT.DAT' FILE
C       SEE 'README.DOC' ON FLOPPY DISK FOR DETAILS
C      FEET TO METERS
       F2M(X) = X*0.3048
C      DEGREE F TO DEGREE C
       F2C(X) = (X-32)*0.5556
C      DEGREE C TO DEGREE F
       C2F(X) = X *1.8+32.0
C      DEGREE F/FT TO DEGREE C/M LAPSE RATE
       FF2CM(X) = X*1.822


C BEGIN
      WRITE(*,'(////////////////)')
      WRITE(*,*) '                            MT-CLIM V1.0'
      WRITE(*,*) '                    --------------------------'
      WRITE(*,*) ' '
      WRITE(*,*) '         A MOUNTAINOUS TERRAIN MICROCLIMATE SIMULATOR'
      WRITE(*,*) '                         COPYRIGHT 1989.'
      WRITE(*,*) ' '
      WRITE(*,*) ' '
      WRITE(*,*) '  1) MT-CLIM NEEDS A INITIALIZATION FILE DESCRIBING'
      WRITE(*,*) '     THE SITE AND BASE STATION.'
      WRITE(*,*) ' '
      WRITE(*,*) ' '
      WRITE(*,*) ' '
      WRITE(*,*) ' '
      WRITE(*,*) ' '
      WRITE(*,*) ' WHAT IS THE INITILIZATION FILE? (12 CHARACTER LIMIT)'
      WRITE(*,*) ' '
      READ(*,90)IFILE
90    FORMAT(1A12)

C     FORMAT FOR INIT.DAT IS 1 VALUE FOLLOWED BY COMMENTS, PER LINE
C     REAL VALUES IN COLUMNS 1 TO 12, COMMENTS FROM 13 ONWARD.
C     READ 1 VARIABLE PER LINE BY REUSING F12.0 FORMAT STATEMENT

      OPEN(7,FILE=IFILE,STATUS='OLD')
      READ(7,100) BFILE,SFILE
 100  FORMAT(/,/,A12,/,A12)

C     ' ENGLISH(TEMPS IN F & PPT IN INCHES) OR SI(C & CM) UNITS (E OR S) ?')
      READ(7,1) ANS3
C     QUESTIONS NOW ANSWERED IN THE INIT.DAT FILE.
C     ' DO YOU HAVE DEW POINT TEMPERATURE (Y OR N) ?')
      READ(7,1) ANS0
C     ' HOW MANY PRECIP STATIONS (1 OR 2) ?')
      READ(7,*) NPPT
C     ' DO YOU WANT THRESHOLD FOR RADIATION (Y OR N) ?')
      READ(7,1) ANS1
C     ' DO YOU WANT TOTAL RAD OR AVERAGE RAD (T OR A) ?')
      READ(7,1) ANS2
C     ' USE YEARDAY(JULIAN DAY) INPLACE OF MONTH, DAY (Y OR N) ?'
      READ(7,1) ANS4
1     FORMAT(1A1)

      USEENG=(ANS3 .EQ. 'E' .OR. ANS3 .EQ. 'e')
      USEDEW=(ANS0 .EQ. 'Y' .OR. ANS0 .EQ. 'y' )
      USETHR=(ANS1 .EQ. 'Y' .OR. ANS1 .EQ. 'y')
      USETOT=(ANS2 .EQ. 'T' .OR. ANS2 .EQ. 't')
      USEJD=(ANS4 .EQ. 'Y' .OR. ANS4 .EQ. 'y')

      READ(7,*) NDAYS

 101  FORMAT(F12.0)
      READ(7,101) SLAT,SELEV,BELEV,SASPCT,SSLOPE,SLAI
      READ(7,101) SISOH,BISOH(1),BISOH(2),SEHORZ,SWHORZ,SALBDO
      READ(7,101) TRANCF,TEMPCF,TLAPSE,MAXLAP,MINLAP,DLAPSE
      CLOSE(7)

C     CONVERT ENGLISH TO SI
      IF (USEENG) THEN
         SELEV=F2M(SELEV)
         BELEV=F2M(BELEV)
         TLAPSE=FF2CM(TLAPSE)
         MINLAP=FF2CM(MINLAP)
         MAXLAP=FF2CM(MAXLAP)
         DLAPSE=FF2CM(DLAPSE)
      ENDIF

C     COMPUTE ATMOSPHERIC TRANSMISSIVITY
      INFILE = 8
      CALL TRANSM (INFILE,SELEV,TRAN,BFILE,NDAYS,TRANCF,
     &             USEDEW,USEENG,NPPT,USEJD)

C                   INITIALIZE OUTPUT FILE
      OPEN(INFILE,FILE=BFILE,STATUS='OLD')
      OPEN(9,FILE=SFILE)
      CALL INIT9(IFILE,USETOT,USEENG,USEJD)

C     FOR EACH DAY COMPUTE AND PRINT MICROCLIMATE
      DO 11 J=1,NDAYS
C        ECHO TO SCREEN TO INDICATE PROGRAM IS WORKING
         PRINT*,J

         CALL BREAD(INFILE,JDAY,BMAX,BMIN,BDEW,BPPT,NPPT,
     &              USEDEW,USEENG,USEJD)

C        THE VARIABLE FLAG IS SET TO 1 FOR SLOPED TERRAIN; FLAT TO 0
         FLAG=1.0
         CALL RAD (SLAT,SASPCT,SSLOPE,JDAY,SEHORZ,SWHORZ,TRAN,
     1            SALBDO,SRAD,SARAD,FLAG,USETHR,USETOT)
         FLAG=0
         CALL RAD (SLAT,SASPCT,SSLOPE,JDAY,SEHORZ,SWHORZ,TRAN,
     1            SALBDO,FRAD,FARAD,FLAG,USETHR,USETOT)

         CALL TEMP (BMAX,BMIN,SELEV,BELEV,SLAI,MLAI,SARAD,FARAD,
     1            TLAPSE,TEMPCF,MAXLAP,MINLAP,STEMP,SMIN,SMAX)
         CALL HUMD (BDEW,DLAPSE,STEMP,SELEV,BELEV,SHUMD)
         IF(SHUMD .GE. 100.0) SHUMD=99.0

         CALL RAIN (NPPT,BPPT,BISOH,SISOH,SPPT)

         CALL WRITE9(JDAY,SRAD,STEMP,SMAX,SMIN,SHUMD,SPPT,
     &               USEENG, USEJD)

 11   CONTINUE
      STOP
      END
C END MTCLIM


CM SUBROUTINE TRANSM (INFILE,SELEV,TRAN,BFILE,NDAYS,TRANCF,
CM                    USEDEW,USEENG,NPPT,USEJD)
CD
CD  ORIGINAL CODE WRITTEN BY R. K. NEMANI
CD  REWRITTEN BY J. C. COUGHLAN ON 4-1-89
CD  MICROSOFT FORTRAN VERSION 4.1 IN STANDARD, TRANSPORTABLE FORTRAN
CD
CD    THIS SUBROUTINE COMPUTES TRANSMISSIVITY FOR EACH DAY BASED ON
CD    BRISTOW, K.L. AND G.S. CAMPBELL 1984  On the relationship between
CD    incoming solar radiation and daily maximum and minimum temperature,
CD    Agric. For. Meteorol. 31, 159-66
CD
CV  PARAMETERS
CVF  FILES:
CVF    INFILE :UNIT NUMBER FOR BASE STATION CLIMATIC INPUT FILE
CVI  INPUT:
CVI    SELEV  :SITE ELEVATION IN METERS
CVI    BFILE  :INPUT FILE NAME OF BASE STATION MET. DATA
CVI    NDAYS  :TOTAL NUMBER OF DAYS TO SIMULATE
CVI    USEDEW :FLAG TO DEW POINT PRESENT IN BFILE
CVI    NPPT   :NUMBER OF PREC STATIONS IN BFILE
CVO  OUTPUT:
CVO    TRANCF :CLEAR SKY TRANSMISSIVITY AT SEA LEVEL
CVO    TRAN   :ATMOSPHERIC TRANSMISSIVITY ARRAY
CVI
CVL  LOCAL VARIABLES
CVL    TAMP  :TEMPERATURE AMPLITUDE
CVL    DRAIN :RAINY DAYS
CVL    CLTRAN:CLEAR SKY TRANSMISSIVITY
CVL    XTRANS:ACTUAL TRANSMISSIVITY
CVL    TRANCF:SEA LEVEL CLEAR SKY TRANS (D.M.GATES 1980 BIOPHYSICAL ECOLOGY)
CVL    PCTRAN:% TRANSMITTANCE OF CLEAR SKY POTENTIAL. TEMP AMP OF 20 = 100%
CVL    PPTMIN:MIN PPT FOR REDUCING ATMOSPHERIC TRANS (CM)
CVL    TRANMN:(CONSTANT) MINIMUM TRANSMITTANCE IN %
CVL    JDAY  :YEARDAY FROM INPUT FILE
CVL    JDAY1 :FIRST DAY IN THE INPUT FILE
CVL    JDAYL :LAST DAY IN THE FILE
CVL    KDAY  :POINTER TO ARRAYS REFERENCED BY YEARDAY

      SUBROUTINE TRANSM (INFILE,SELEV,TRAN,BFILE,NDAYS,TRANCF,
     &                   USEDEW,USEENG,NPPT,USEJD)
      REAL AMP(365),DRAIN(365),TRAN(365)
      REAL SELEV,TRANCF
      LOGICAL USEDEW,USEENG
      CHARACTER *12 BFILE
      INTEGER NPPT, INFILE, NDAYS
      LOGICAL USEJD

      REAL BPPTA(2), BPPT,BMAX,BMIN,BDEW,BMAX1,BMIN1,BPPT1,TAMP
      REAL DIFF,CLTRAN, PCTRAN,XTRANS
      REAL PPTMIN, TRANMN
      INTEGER K,KDAY,JDAY,JDAY1,JDAYL,I,M,J

C  CONSTANTS
      DATA PPTMIN/0.254/
      DATA TRANMN/0.1/

C  BEGIN
      K=0
      KDAY=0
      OPEN(INFILE,FILE=BFILE,STATUS='OLD')
      REWIND(INFILE)

C     READ FOR THE 1ST DAY TO INITIALIZE THE LOOP
      CALL BREAD(INFILE,JDAY,BMAX,BMIN,BDEW,BPPTA,NPPT,
     &           USEDEW,USEENG,USEJD)
      IF (NPPT .EQ. 2) THEN
         BPPT = (BPPTA(1) + BPPTA(2)) / 2
      ELSE
         BPPT = BPPTA(1)
      ENDIF

C     MAKE NOTE OF THE FIRST DAY
      JDAY1 = JDAY

C     READ TO OBTAIN TEMP AMPLITUDE
      DO 12 J=1, NDAYS-1

         CALL BREAD(INFILE,KDAY,BMAX1,BMIN1,BDEW,BPPTA,NPPT,
     &              USEDEW,USEENG,USEJD)
         IF (NPPT .EQ. 2) THEN
           BPPT1 = (BPPTA(1) + BPPTA(2)) / 2
         ELSE
           BPPT1 = BPPTA(1)
         ENDIF

C        COMPUTE AMPLITUDE
         TAMP = BMAX - ((BMIN + BMIN1)/2.)

C        RAINY DAY CORRECTIONS
         IF (BPPT .GT. PPTMIN) THEN
           K=K+1
           DRAIN(K) = JDAY
           TAMP =  TAMP* 0.75
         ENDIF

C        SWITCH THE VALUES FROM J+1 DAY TO J
         AMP(JDAY) = TAMP
         JDAY = KDAY
         BMAX = BMAX1
         BMIN = BMIN1
         BPPT=BPPT1
 12   CONTINUE

      JDAYL = JDAY

C     EXCEPTION FOR THE LAST DAY
 121  TAMP = BMAX - BMIN
      IF (BPPT .GT. PPTMIN ) TAMP = TAMP*0.75
      AMP(JDAY) =TAMP
 
C     CORRECT AMPL VALUES FOR PRE-RAINY DAYS
      DO 30 I=1,K
         KDAY = DRAIN(I)
C        CANNOT CORRECT THE 1ST AND 2ND DAY SO SKIP
         IF ((KDAY-2) .GE. JDAY1) THEN
            DIFF = AMP(KDAY-2) - AMP(KDAY-1)
            IF (DIFF.GE.2.0) AMP(KDAY-1) = AMP(KDAY-1) * 0.75
         ENDIF
30    CONTINUE

C     COMPUTE CLEAR SKY TRANSMITTANCE AT SITE (GATES 1980)
      CLTRAN = TRANCF + SELEV * 8.0E-5
      CLTRAN = AMIN1(CLTRAN,1.0)

C     COMPUTE TRANSMISSIVITY FOR EACH DAY (BRISTOW AND CAMPBELL 1984)
      DO 40 M = JDAY1, JDAYL
	 PCTRAN = (1 - EXP(-0.003*(AMP(M)**2.4)))
	 XTRANS = CLTRAN * PCTRAN
	 TRAN(M) = AMAX1(XTRANS,TRANMN)
 40   CONTINUE

C     REWIND FILE 8 TO USE AGAIN
      REWIND(8)
      CLOSE (8)
      RETURN
      END


CM  SUBROUTINE RAD
C
      SUBROUTINE RAD(SLAT,SASPCT,SSLOPE,JDAY,SEHORZ,SWHORZ,TRAN,
     & ALBDO,RADN,ARAD,FLAG,USETHR,USETOT)
      INTEGER FLAG
      REAL SLAT,SASPCT,SSLOPE,SEHORZ,SWHORZ
      REAL ALBDO,ARAD,RADN
      REAL DEC(46),SOLCON(12),A(21),TRAN(365)
      LOGICAL USETHR,USETOT

CD     VALIDATED TO BUFFO BY JCC ON 4-19-89
CD     RAD COMPUTES INCIDENT SHORTWAVE RADIATION AND
CD     NET SHORTWAVE RADIATION FOR ANY GIVEN DAY BASED ON SURFACE
CD     CHARACTERISTICS, SUN-EARTH GEOMETRY, TRANSMISSIVITY.
CD  ORIGINAL CODE WRITTEN BY R. K. NEMANI
CD  REWRITTEN BY J. C. COUGHLAN ON 4-1-89
CD  MICROSOFT FORTRAN VERSION 4.1 IN STANDARD, TRANSPORTABLE FORTRAN
CD
CD  THIS CODE CAN BE OPTIMIZED TO THE POINT WHERE IT CUTS TOTAL PROGRAM
CD  EXECUTION TIME TO HALF OF THE PRESENT.  REMOVE AS MUCH FROM THE DO
CD  LOOP AS POSSIBLE.  THE DO LOOP INCREMENTS 10 MINUTES FOR A 24 HOUR
CD  DAY.  SOME CALCULATIONS CAN BE DONE ONCE OUTSIDE THE LOOP I.E.
CD  COS(DSLOPE), AND SAVED IN A VARIABLE TO BE USED IN THE LOOP.  WE
CD  HAVE KEPT IT LESS EFFICIENT TO REDUCE COMPLEXITY.  EFFICIENT CODE
CD  IS LESS READABLE, USUALLY.  ALSO, HAVE HARDWARE SUPPORT FOR MATH
CD  REALLY HELPS TOO, ESP. FOR TRIG. FUNCTIONS. -JCC
CD
CV  VARIABLES
CVI  INPUT
CVI     SLAT: SITE LATITUDE   DEGREES
CVI     SASPECT:SITE ASPECT  %
CVI     SSLOPE:SITE SLOPE  %
CVI     JDAY: CURRENT YEARDAY
CVI     SEHORZ:SITE EAST HORIZON IN DEGREES FROM 0
CVI     SWHORZ:SITE WEST HORI IN DEGREES FROM 0
CVI     TRAN:TRANSMISSIVITY ARRAY   %
CVI     ALBDO:SITE ALBEDO   %
CVI     FLAG:1 MEANS SLOPING TERRAIN, 0 MEANS FLAT SURFACE.
CVI     USETHR:THRESHOLD RADIATION OF 70 W/M^2
CVI     USETOT:TOTAL RADIATION IF = TO T   KJ/M^/DAY
CVI          ELSE DAYLIGHT AVERAGE   W/M^2
CVO  OUTPUT:
CVO     ARAD:ABSORBED RADIATION, W/M2
CVO     RADN:INCIDENT RADIATION  W/M2
CV
CV
CVL
CVL   LOCAL VARIABLE LIST:
CVL     SOL   =SOLAR CONSTANT DERIVED FROM SOLCON ARRAY
CVL     AM    =OPTICAL AIR MASS FOR ANGLES > 21 DEGREES
CVL     A     =OPTICAL AIR MASS ARRAY FOR ANGLES BETWEEN
CVL            0 AND 21 DEGREES ABOVE HORIZON
CVL     DECL  =DECLINATION
CVL     JDAY  =DAY OF YEAR
CVL     ASP   =ASPECT IN DEGREES
CVL     DSLOP =SLOPE IN DEGREES
CVL     H     =ANGLE OF SUN FROM SOLAR NOON
CVL     TRAN  =TRANSMISSIVITY CONSTANT
CVL     TRAM  =TRANSMISSIVITY CORRECTED FOR AIR MASS
CVL     NNH   =CALCULATION INTERVAL IN SECONDS. 600 = 10 MINUTES
CVL     NC    =SECONDS IN ONE DAY (24 HOURS)
CVL     N     =NUMBER OF INTERVALS OF LENGTH NNH IN ONE DAY
CVL     DT    =DIRECT SOLAR PERPENDICULAR TO SUN ON THE
CVL            OUTSIDE OF ATMOSPHERE FOR INTERVAL (KJ/M**2)
CVL     ETF   =DIRECT SOLAR ON OUTSIDE OF ATMOSPHERE
CVL            PARALLEL TO EARTHS SURFACE FOR INTERVAL
CVL     GRAD  =TOTAL DAILY RADIATION AT GIVEN LOCATION (KJ/M**2)
CVL     HRAD  =DIRECT SOLAR ON EARTHS SURFACE (FLAT)       "
CVL     TDIF  =TOTAL DAILY DIFFUSE RADIATION               "
CVL     DIFRAD=DIFFUSE ON SLOPE FOR INTERVAL               "
CLL     DRAD  =DIRECT ON SLOPE FOR INTERVAL                "
CVL     CZA   =COSINE ZENITH ANGLE
CVL     CBSA  =COSINE BEAM SLOPE ANGLE
CVL     GLOBF =GLOBAL RADIATION, DETERMINING DIFFUSE
CVL     DIFFL =DIFFUSE ON FLAT FOR INTERVAL
CVL     DAYL  =DAYLENGTH (HOURS)
CVL     ALBDO =ALBEDO
CVL
C  LOCAL VARIABLES
      INTEGER NNH,NC,N,MO,IDEC,NH,K,ML
      REAL CONV,X,ASP,DLAT,SLOPE,DSLOP,XTRAN,DECL,GRAD,TDIF
      REAL DAYL2,DH,CZA,H,AM,TRAM,DT,ETF,HRAD,GLOBF,DIFFL,CBSA
      REAL DRAD,SE,DIFRAD,RADT,AVERAD,DAYL,SOL,TDRAD
C  CONSTANTS
      DATA A/2.90,3.05,3.21,3.39,3.69,3.82,4.07,4.37,4.72,5.12
     1,5.60,6.18,6.88,7.77,8.90,10.39,12.44,15.36,19.79,26.96
     2,30.00/
      DATA DEC/-23.,-22.,-21.,-19.,-17.,-15.,-12.,-9.,-6.,-3.,
     10.,3.,6.,9.,12.,14.,17.,19.,21.,22.,23.,23.5,23.5,23.,
     221.5,20.,18.,16.,14.,12.,9.,6.,3.,0.,-3.,-6.,-9.,-12.,
     3-15.,-17.,-19.,-21.,-22.,-23.,-23.5,-23.5/
      DATA SOLCON/1.445,1.431,1.410,1.389,1.368,1.354,1.354
     1,1.375,1.403,1.424,1.438,1.445/
C  FUNCTIONS
C     CONVERSION STATEMENT FUNCTION FOR DEGREES TO RADIANS
      CONV(X)=X/57.296
C  BEGIN
C     SET THE SLOPE & ASPECT VALUES DEPENDING ON THE FLAG VALUE
      IF (FLAG .GT. 0) THEN
         ASP=conv(SASPCT)
         SLOPE=SSLOPE
      ELSE
         ASP=0
         SLOPE=0
      ENDIF

C     CONVERT PERCENT SLOPE TO DEGREES and to radians
      DLAT=conv(SLAT)
      DSLOP=ATAN(SLOPE/100.)*57.29
      DSLOP=conv(DSLOP)
      XTRAN=TRAN(JDAY)
      NNH=600
      NC=86400
      N=IFIX(NC/NNH+1.)
      DAYL=0.
      MO=IFIX(JDAY/30.+1.)
      IF (MO.GT.12)MO=12

C     SOLCON ARRAY IS IN UNITS OF KW/M**2
      SOL=SOLCON(MO)
      IDEC=IFIX(1.+JDAY/8.)
      DECL=conv(DEC(IDEC))
      GRAD=0.
      TDIF=0.
      TDRAD=0.
      NH=0
      DAYL2=0.

C     DO LOOP INCREMENTS 10 MINUTES, STOPS AFTER 24 HOURS
      DO 11 K=1,N
         NH=NH+NNH
C        DETERMINE ANGLE FROM SOLAR NOON
	 DH=(NH-43200)*.0041667
	 h=conv(dh)
	 CZA=COS(DECL)*COS((DLAT))*COS((H))+SIN((DECL))
     1   *SIN((DLAT))
         IF(CZA.GT. 0.) THEN
C           DAYLENGTH BASED ON SOLAR ELEVATION ABOVE A FLAT HORIZON
            DAYL2=DAYL2+(NNH/3600.)
C           NEXT 6 LINES, DETERMINE OPTICAL AIR MASS
            AM=1./(CZA+.0000001)
            IF(AM.GT.2.9) THEN
               ML=IFIX(ACOS(CZA)/.0174533)-69
               IF(ML.LT.1)ML=1
               IF(ML.GT.21)ML=21
               AM=A(ML)
            ENDIF
            TRAM=XTRAN**AM
            DT=SOL*NNH
            ETF=CZA*DT
            HRAD=ETF*TRAM
            DT=DT*TRAM
            GLOBF=SQRT(HRAD*ETF)
            DIFFL=GLOBF*(1.-GLOBF/ETF)
            CBSA=-SIN(DSLOP)*SIN(ASP)*SIN(H)
     1      *COS(DECL)+(-COS(ASP)*SIN(DSLOP)
     2      *SIN(DLAT)+COS(DSLOP)*COS(DLAT))
     3      *COS(DECL)*COS(H)+(COS(ASP)*SIN(DSLOP)
     4      *COS(DLAT)+COS(DSLOP)*SIN(DLAT))
     5      *SIN(DECL)
            IF(CBSA.GE.0.) THEN
               DRAD=CBSA*DT
C              THE FOLLOWING THREE LINES COMPUTES A TOPOGRAPHIC REDUCTION OF
C              DIRECT RADIATION
C              EHE = EAST HORIZON ELEVATION (DEGREES)
C              WHE = WEST HORIZON ELEVATION (DEGREES)
               SE=1.57-ACOS(CZA)
               IF(DH.LT.0.AND.SE.LT.CONV(SEHORZ))DRAD=0.
               IF(DH.GT.0.AND.SE.LT.CONV(SWHORZ))DRAD=0.
            ELSE
               DRAD=0
            ENDIF

            DIFRAD=DIFFL*(COS((DSLOP*.5))**2.)
            RADT=(DRAD+DIFRAD)/FLOAT(NNH)
C           THE FOLLOWING PROVIDES A MINIMUM RADIATION THRESHOLD OF
C           70 W/M**2 (0.1 LY/MIN) FOR DAYLENGTH AND RADIATION SUMMATION.
            IF(USETHR) THEN
               IF(RADT.GT.0.07) THEN
                  DAYL = DAYL + (NNH/3600.)
                  GRAD = GRAD + DRAD + DIFRAD
                  TDIF = TDIF + DIFRAD
                  TDRAD = TDRAD + DRAD
               ENDIF
            ELSE
               DAYL = DAYL + (NNH/3600.)
               GRAD = GRAD + DRAD + DIFRAD
               TDIF = TDIF + DIFRAD
               TDRAD = TDRAD + DRAD
            ENDIF
         ENDIF
11    CONTINUE
      IF(USETOT) THEN
         RADN = GRAD
         ARAD = RADN * (1- ALBDO)
      ELSE
         AVERAD=(GRAD/(DAYL*3600.))
C        CONVERT KW/M**2 to W/M**2
         RADN = AVERAD * 1000.0
         ARAD = RADN * (1-ALBDO)
      ENDIF
      RETURN
      END

CM
CM      SUBROUTINE TEMP(BMAX,BMIN,SELEV,BELEV,SLAI,MLAI,SARAD,
CM                     FARAD,TLAPSE,TEMPCF,MAXLAP,MINLAP,STEMP,SMIN,SMAX)
CM
CD     TEMP COMPUTES SITE DAYLIGHT AVERAGE TEMPERATURE, AND SITE MAX &
CD     MIN TEMPS BASED ON BASE STATION DATA AND THEN CORRECTS IT FOR
CD     ELEVATION, SLOPE AND ASPECT.
CD
CVI PARAMETER LIST
CVI
CVI INPUT:
CVI   BMAX   :BASE STATION MAX TEMP
CVI   BMIN   :BASE STATION MIN TEMP
CVI   SELEV  :SITE ELEVATION IN METERS
CVI   BELEV  :BASE ELEVATION IN METERS
CVI   SLAI   :SITE LAI (ALL SIDES)
CVI   MLAI   :MAXIMUM LAI (ALL SIDES) FOR MTCLIM ALGORITHMS
CVI   SARAD  :SITE RADIATION KJ/M^2/DAY
CVI   FARAD  :FLAT SURFACE RADIATION AT SITE KJ/M^2/DAY
CVI   TLAPSE :LAPSE RATE FOR DAYLIGHT AVE TEMP  C/1000M
CVI   TEMPCF :CONSTANT FOR DAYLIGHT AVE. TEMP (RUNNING ET.AL. EQ. 1)
CVI   MAXLAP :LAPSE RATE FOR MAX TEMPS C/1000M
CVI   MINLAP :LAPSE RATE FOR MIN TEMPS C/1000M
CVO OUTPUT:
CVO   STEMP  :SITE DAYLIGHT AVE TEMP C
CVO   SMIN   :SITE MIN TEMP C
CVO   SMAX   :SITE MAX TEMP C
CVL
CVL LOCAL VARIABLE LIST
CVL   MLAPSE: LAPSE RATE FOR NIGHT MIN TEMPERATURES
CVL   DAYLT: DAYLIGHT AVERAGE TEMPERATURE
CVL   TSYNOP: SYNOPTIC TEMPERATURE AT THE SITE
CVL   RADRAT: RATIO OF FLAT AND SLOPE RADIATION
CVL   TADD : TEMPERATURE INCREMENT FOR SOUTH SLOPES
CVL   TSUB : TEMPERATURE DECREMENT FOR NORTH SLOPES
CVL   MLAI : MAXIMUM LAI (ALL SIDES)

      SUBROUTINE TEMP (BMAX,BMIN,SELEV,BELEV,SLAI,MLAI,SARAD,
     1      FARAD,TLAPSE,TEMPCF,MAXLAP,MINLAP,STEMP,SMIN,SMAX)

      REAL BMAX, BMIN, SELEV, BELEV, SLAI, MLAI, SARAD, TEMPCF
      REAL FARAD, TLAPSE, MAXLAP, MINLAP,STEMP, SMIN, SMAX

C  LOCAL VARIABLES
      REAL LAI,TMEAN
      REAL DAYLT,DELEV,TSYNOP,RADRAT,TSUB,TADD

C IF SITE LAI > MAX LAI THEN USE MAX LAI IN EQUATIONS BELOW
      IF (SLAI .GT. MLAI) THEN
        LAI = MLAI
      ELSE
        LAI = SLAI
      ENDIF

C     COMPUTE DAYLIGHT AVERAGE TEMP
      TMEAN= (BMAX+BMIN)/2.0
      DAYLT= ((BMAX-TMEAN)*TEMPCF) + TMEAN

C                CORRECT FOR LAPSE RATE
      DELEV = SELEV-BELEV
      TSYNOP = DAYLT - TLAPSE * (DELEV / 1000.0)
      SMIN = BMIN - (DELEV/1000.0 * MINLAP)
      SMAX = BMAX - (DELEV/1000.0 * MAXLAP)

C                COMPUTE THE RATIO OF SLOPE AND FLAT RADIATION
      RADRAT = SARAD/FARAD

C                ADJUST SYNOPTIC TEMP TO OBTAIN SLOPE TEMP
C     ADDITIONS MADE ON 10/88 TO ADJUST MAX TEMPS TO RADIATION.
C     NOT VALIDATED BY LITERATURE.
C     NOT ADJUSTING MINIMUM TEMPERATURE BECAUSE LONGWAVE NIGHT ADJUSTMENT.

      IF (RADRAT.LT.1.0) THEN
         TSUB=((1/RADRAT)*(1+(LAI/MLAI)))
         STEMP = TSYNOP -  TSUB
         SMAX = SMAX - TSUB
      ELSE
         TADD=(RADRAT* (1-(LAI/MLAI)))
         STEMP = TSYNOP + TADD
         SMAX = SMAX + TADD
      ENDIF
      RETURN
      END



      SUBROUTINE HUMD (BDEW,DLAPSE,STEMP,SELEV,BELEV,SHUMD)
      REAL BDEW, DLAPSE,STEMP,SELEV,BELEV,SHUMD
CM  SUBROUTINE HUMD (BDEW,DLAPSE,STEMP,SELEV,BELEV,SHUMD)
CM                   REAL BDEW, DLAPSE,STEMP,SELEV,BELEV,SHUMD
CD  HUMD COMPUTES RELATIVE HUMIDITY BASED ON BASED
CD  ON BASE STATION DEW POINT
CVI  PARAMETERS
CVI INPUT:
CVI    BDEW   :BASE STATION BEW POINT C
CVI    DLAPSE :DEW POINT LAPSE RATE
CVI    STEMP  :DAYLIGHT AVE. TEMP. LAPSE RATE
CVI    SELEV  :SITE ELEVATION IN METERS
CVI    BELEV  :BASE ELEVATION IN METERS
CVO OUTPUT:
CVO    SHUMD  :SITE HUMIDITY %
CVL
CVL  LOCAL VARIABLES
CVL    SDEW   :SITE DEW POINT  C
CVL    ES     :AMBIENT VAPOR DENSITY
CVL    ESD    :SATURATED VAPOR DENSITY
      REAL SDEW,ES,ESD

C     CORRECT DEW PNT FOR LAPSE
      SDEW = BDEW - (DLAPSE * ((SELEV-BELEV)/1000.))
C     COMPUTE RELATIVE HUMIDITY
      ES = 6.1078 * EXP((17.269*SDEW)/(237.3 + SDEW))
      ESD = 6.1078 * EXP((17.269 * STEMP)/(237.3 + STEMP))
      SHUMD=(ES/ESD) * 100.0
      RETURN
      END



      SUBROUTINE RAIN (NPPT,BPPT,BISOH,SISOH,SPPT)
      INTEGER NPPT
      REAL BISOH(2),SISOH,SPPT,BPPT(2)
      REAL RAT1,RAT2

C     NPPT: NUMBER OF BASE PPT STATIONS

C     THIS SUBROUTINE COMPUTES SITE PRECIPITATION BY MULTIPLYING
C     THE BASE STATION PRECIP WITH THE RATIO OF BASE AND SITE
C     ISOHYETS.

C     COMPUTE RATIO OF THE ISOHYTE(S)
      RAT1 = SISOH/BISOH(1)
      IF (NPPT .EQ. 2) THEN
         RAT2 = SISOH/BISOH(2)
         SPPT = ( BPPT(1)*RAT1 + BPPT(2)*RAT2 ) / 2.0
      ELSE
         SPPT = RAT1*BPPT(1)
      ENDIF
      RETURN
      END


       SUBROUTINE MONDAY(JDAY,MON,DAY)
       INTEGER JDAY,MON,DAY
       INTEGER DAYS(12)
       DATA DAYS/31,28,31,30,31,30,31,31,30,31,30,31/
       DAY=JDAY
       MON = 1
10     IF(DAY-DAYS(MON) .GE. DAYS(MON) ) THEN
          DAY=DAY-DAYS(MON)
          MON = MON+1
          GO TO 10
       ENDIF
       RETURN
       END


       REAL FUNCTION YEARDAY(MON,IDATE)
       INTEGER DAYS(12), I,JDAY,MON
       DATA DAYS/31,28,31,30,31,30,31,31,30,31,30,31/
       JDAY = 0
       DO 10 I = 1, (MON-1)
          JDAY= JDAY+DAYS(I)
10     CONTINUE
       JDAY=JDAY+IDATE
       YEARDAY=JDAY
       RETURN
       END

       
       SUBROUTINE WRITE9(JDAY,SRAD,STEMP,SMAX,SMIN,SHUMD,SPPT,
     & USEENG, USEJD)
       INTEGER JDAY,DAY,MON
       REAL SRAD,STEMP,SMAX,SMIN,SHUMD,SPPT,C2F,X
       LOGICAL USEJD
       LOGICAL USEENG
C      DEGREE C TO DEGREE F
       C2F(X) = X *1.8+32.0

       IF (USEENG) THEN
c        CONVERT TEMPERATURES TO DEG F IF INDICATED
            STEMP = C2F(STEMP)
            SMIN = C2F(SMIN)
            SMAX = C2F(SMAX)
            SPPT = SPPT/2.54
       ENDIF

       IF (USEJD) THEN
         WRITE(9,110) JDAY,SRAD,STEMP,SMAX,SMIN,SHUMD,SPPT
 110     FORMAT(1X,I3,2X,F8.2,4(F6.0),F5.2)
       ELSE
         CALL MONDAY(JDAY,MON,DAY)
         WRITE(9,111) MON,DAY,SRAD,STEMP,SMAX,SMIN,SHUMD,SPPT
 111     FORMAT(1X,I2,2X,I3,2X,F8.2, 4(F6.0),F5.2)
       ENDIF
       RETURN
       END


      SUBROUTINE INIT9(IFILE,USETOT,USEENG,USEJD)
      LOGICAL USETOT,USEENG
      CHARACTER*80 LINE1, LINE3
      CHARACTER*17 LINE2
      CHARACTER*12 IFILE
      LOGICAL USEJD

C IDENTIFY THE OUTPUT FILE
      WRITE(9,*)' MTCLIM OUTPUT FILE'
      WRITE(9,*)' A LISITNG OF INITILIZATION FILE FOLLOWED BY OUTPUT:'
      WRITE(9,*)

C LIST INITALIZATION FILE IN OUTPUT FILE
      OPEN(UNIT=7,FILE=IFILE)
2     READ(7,1,END=10) LINE1
1        FORMAT(1A80)
         WRITE(9,1)LINE1
         GO TO 2
10    CONTINUE
      WRITE(9,*)

C SET UP HEADINGS FOR OUTPUT FILE
      IF (USEJD) THEN 
         LINE1='JDAY  RADIATION  STEMP MAXT  MINT  RH   PPT'
      ELSE
         LINE1='MON  DAY RADIATION  STEMP MAXT  MINT  RH   PPT'
      ENDIF

      IF (USETOT) THEN
        LINE2= '        KJ/M**2 '
      ELSE
        LINE2= '        W/M**2 '
      ENDIF

      IF (USEENG) THEN
        LINE3=LINE2//'F     F     F     %     INCHES'
      ELSE
        LINE3=LINE2//'C     C     C     %     CM '
      ENDIF

      WRITE(9,111)LINE1,LINE3
111   FORMAT(1A80)

      RETURN
      END

CM  SUBROUTINE BREAD: READ BASE STATION DATA.
      SUBROUTINE BREAD(INFILE,JDAY,BMAX,BMIN,BDEW,BPPT,NPPT,
     &                 USEDEW,USEENG,USEJD)
      INTEGER INFILE, JDAY, NPPT
      REAL BMAX,BMIN
      REAL BPPT(2)
      LOGICAL USEDEW, USEENG
      LOGICAL USEJD
CVI VARIABLES
CVI  INPUT:
CVF       INFILE: INTEGER UNIT NUMBER OF BASE STATION CLIMATIC INPUT FILE
CVI       USEDEW: IS DEW POINT SUPPLIED IN FILE
CVI            IF NOT THEN SET IT TO NIGHT MINIMUM TEMPERATURE.
CVI      USEENG: IS INPUT IN ENGLISH OR SI UNITS
CVO  OUTPUT:
CVO       JDAY: YEARDAY
CVO       BMAX: MAXIMUM TEMPERATURE
CVO       BMIN: MINIMUM TEMPERATURE
CVO       BDEW: DEW POINT TEMPERATURE
CVO       NPPT: NUMBER OF BASE STATIONS 1 OR 2
CF    FUNCTION F2C: DEGREE F TO DEGREE C
      REAL F2C,X,YEARDAY
      INTEGER DAY, MON
      F2C(X) = (X-32)*0.5556

C BEGIN
      IF (USEJD) THEN
        IF(USEDEW) THEN
          IF(NPPT .EQ. 1) THEN
             READ(INFILE,*) JDAY, BMAX, BMIN, BDEW, BPPT(1)
             BPPT(2) = 0.0
          ELSE
             READ(8,*) JDAY, BMAX, BMIN, BDEW, BPPT(1), BPPT(2)
          ENDIF
        ELSE
          IF(NPPT .EQ. 1) THEN
             READ(INFILE,*) JDAY, BMAX, BMIN, BPPT(1)
             BPPT(2) = 0.0
             BDEW=BMIN
          ELSE
             READ(INFILE,*) JDAY, BMAX, BMIN, BPPT(1), BPPT(2)
             BDEW = BMIN
          ENDIF
        ENDIF
      ELSE
        IF(USEDEW) THEN
          IF(NPPT .EQ. 1) THEN
             READ(INFILE,*) MON,DAY, BMAX, BMIN, BDEW, BPPT(1)
             BPPT(2) = 0.0
          ELSE
             READ(8,*) MON,DAY, BMAX, BMIN, BDEW, BPPT(1), BPPT(2)
          ENDIF
        ELSE
          IF(NPPT .EQ. 1) THEN
             READ(INFILE,*) MON,DAY, BMAX, BMIN, BPPT(1)
             BPPT(2) = 0.0
             BDEW=BMIN
          ELSE
             READ(INFILE,*) MON,DAY, BMAX, BMIN, BPPT(1), BPPT(2)
             BDEW = BMIN
          ENDIF
        ENDIF
        JDAY=YEARDAY(MON,DAY)
      ENDIF

      IF (USEENG) THEN
C CONVERT ENGLISH UNITS TO SI
         BPPT(1) = BPPT(1) * 2.54
         BPPT(2) = BPPT(2) * 2.54
         BMAX = F2C(BMAX)
         BMIN = F2C(BMIN)
         BDEW = F2C(BDEW)
      ENDIF

      RETURN
      END
C
C ---------------- END LISTING -----------------

CD COMMENT GIUDE:
CD
CD  'C'   ANY OLD COMMENT NOT TO BE EXTRACTED TO A DOCUMENTATION FILE
CD  'CV'  VARIABLE
CD  'CVF' FILE PARAMETER/VARIABLE
CD  'CVI' INPUT PARAMETER/VARIABLE
CD  'CVO' OUTPUT PARAMETER/VARAIBLE
CD  'CVL' LOCAL VARIABLE
CD  'CM'  MODULE NAME AND CALLINF CONVENTION
CD  'CD'  DESCRIPTION OF MODULE
CD  'CF'  STATEMENT FUNCTIONS
