C
        PROGRAM FLUXES
C
C  This program computes air/sea fluxes and some boundary layer
C  characterisitcs from data entered at the terminal.  The crucial piece of
C  the program is the subroutine BULK, which follows the theory given by
C  Large and Pond, JPO, (11) 324-336, 1981.  The code for BULK was written 
C  by Bill Large, and modified in small ways by Bob Weller and Jim Price.  
C  The main program was written by Jim Price (who is responsible for any and
C  all mistakes !) in January, 1991. Please contact him if you find an error
C  or can suggest an improvement.  
C
C  The program assumes that unit 6 is the terminal for input, unit 5 is
C  the terminal for output, a file FLUXA3.LIS is opened as unit 1 for 
C  output, and if requested, a file FLUXA3.DAT is opened for output of
C  data in a format suitable for plotting.
C
C  Input data are as follows, all assumed to be in MKS units:
C
C  Observation height; the height above the surface at which the atmospheric
C  observations were made.  The surface temperature is presumed to be made
C  at Z = 0. Reasonable values are in the range 1 to 30 m.
C
C  Wind speed; the mean (say 10 min average) wind at the observation height.
C  Reasonable range 1. to 40. m/sec.
C
C  Air temperature; range -10. to 35. C.
C  
C  Humidity can be entered in one of four ways:
C    Relative humidity; range 0. to 100 (percent). A good default
C     is to take relative humidity = 70.
C    Specific humidity; reasonable range 0. to 30. (g/kg)
C    Dew point temperature; range -10. to 30. (C)
C    Wet bulb temperature; range -10. to 30. (C).
C  The program will first ask you to specify which variable, and
C   then ask for the data.
C
C  Sea surface temperature; range say 0. to 30. C.
C
C  If you request to do a radiation estimate then you will be asked
C   to enter the following additional data:
C
C  Year day; range 0 to 365.99, Jan 10 at noon would be 10.5.
C
C  Latitude (deg) of the observation site.
C
C  Pressure of the atmosphere; reasonable range 900 to 1100 mb. 
C
C  Atmospheric transmission coefficient; range say 0.5 to 1.0.
C
C  Cloud cover; range 0 to 1.0.
C
C
        OPEN (1, FILE='fluxa4.lis')
C
C 
 1000   CONTINUE
C
        WRITE (6,1)
 1      FORMAT(1X,//,1X,
     x 'This program computes air/sea fluxes from data',/,
     x 1X,'entered at the terminal. Data should be entered',/,
     x 1X,'in floating point format, (i.e., observation ',/,
     x 1X,'height should be 10.0 rather than 10).',/) 
C
        WRITE(6,6)
 6      FORMAT(1X,/,1X,'Enter observation height (m)  ', $)
        READ(5,900) ZU
 900    FORMAT (F14.3)
        WRITE (1,980)
 980    FORMAT (1X,//,1X,'Input data are:  ',/)
        WRITE (1,90) ZU
 90     FORMAT (1X,'Observation height (m) =',F6.1)
C
        ZT = ZU
        ZQ = ZU
C
        WRITE (6,2)
 2      FORMAT(1X,/,1X,'Enter wind speed (m/sec)  ', $)
        READ (5,900) UZ
        WRITE (1,91) UZ
 91     FORMAT (1X,'Wind speed (m/sec) =',F8.2)
        WRITE(6,3)
 3      FORMAT(1X,/,1X,'Enter air temperature (C)  ', $)      
        READ(5,900) TZ
        WRITE(6,4)
        WRITE (1,92) TZ
 92     FORMAT (1X,'Air temperature (C) =',F8.2)
 4      FORMAT(1X,/,1X,
     X  'The humidity may be specified by either the ',/,
     X  1X,'  relative humidity (percent) (enter 1), or, '/,
     X  1X,'  specific humidity (g/kg) (enter 2), or, ',/,
     X  1X,'  dewpoint temperature (C) (enter 3), or, ',/,
     X  1X,'  wet bulb temperature (C) (enter 4)  ',  $)
        READ (5,909) IQ
 909    FORMAT (I1)
        IF(IQ.EQ.1) THEN
        WRITE (6,12)
 12     FORMAT(1X,/,1X,'Enter relative humidity (percent) ', $)      
        READ(5,900) QZ
        END IF
        IF(IQ.EQ.2) THEN
        WRITE (6,13)
 13     FORMAT(1X,/,1X,'Enter specific humidity (g/kg)  ', $)      
        READ(5,900) QZ
        END IF
        IF(IQ.EQ.3) THEN
        WRITE (6,14)
 14     FORMAT(1X,/,1X,'Enter dewpoint temperature (C) ', $)      
        READ(5,900) QZ
        END IF
        IF(IQ.EQ.4) THEN
        WRITE (6,15)
 15     FORMAT(1X,/,1X,'Enter wet bulb temperature (C) ', $)      
        READ(5,900) QZ
        END IF
        WRITE (1,93) IQ, QZ
 93     FORMAT (1X,'Humidity flag, and humidity variable =',I5,F8.3)
        WRITE(6,5)
 5      FORMAT(1X,/,1X,'Enter sea surface temperature (C)  ', $)
        READ(5,900) TSFC
        WRITE (1,94) TSFC
 94     FORMAT (1X,'Surface temperature (C) = ',F8.2)
        NOSTAB = 0
        WRITE (6,26) 
 26     FORMAT (1X,/,1X,
     x  'If you do NOT want to use stability-dependent transfer',/,
     x  1X,'  formulae then enter 1 ',  $)
        READ (5,920) NOSTAB 
        WRITE (1,95) NOSTAB
 95     FORMAT (1X,'Stability flag =',I3)
C
        WRITE (6,30)
 30     FORMAT (1X,/,1X,'If you do want to estimate solar and ',/,
     x  1X,'  infrared radiation fluxes then enter 1 ',$)
        READ (5,920) IRAD
        IF(IRAD.EQ.0) GO TO 39
C
        WRITE (6,31) 
 31     FORMAT (1X,/,1X,'Enter year day ',$)
        READ (5,900) YDAY
        WRITE (1,96) YDAY
 96     FORMAT (1X,'Year day = ', F6.2)
C
        WRITE (6,32) 
 32     FORMAT (1X,/,1X,'Enter latitude (deg) ',$)
        READ (5,900) RLAT
        WRITE (1,97) RLAT
 97     FORMAT (1X,'Latitude =', F6.2)
C
        WRITE (6,101)
 101    FORMAT (1X,/,1X,'Enter atmospheric pressure (mb) ', $)
        READ (5,900) PRESS
        WRITE (1,102) PRESS
 102    FORMAT (1X,'Pressure (mb) =',F9.1)
C
        WRITE (6,34)
 34     FORMAT (1X,/,1X,
     x  'Enter atmospheric transmission coefficient (0 to 1.) ',$)
        READ (5,900) ATMOS
        WRITE (1,98) ATMOS
 98     FORMAT (1X,'Atmospheric transmission coefficient =', F7.2)
C
        WRITE (6,33) 
 33     FORMAT (1X,/,1X,'Enter cloud cover (0 to 1.) ', $)
        READ (5,900) CC
        WRITE (1,99) CC
 99     FORMAT (1X,'Cloud cover =' F5.2)
C
        CALL SOLAR (YDAY,RLAT,ATMOS,SRAD0,SRADMX,SRADAV,SUMR,DAYTIM)          
        SRAD = SRAD0*(1. - 0.71*CC)
        SRADAV = SRADAV*(1. - 0.71*CC)
C

C
 39     CONTINUE
C
C
        CALL BULK
     X    (UZ, ZU, TZ, ZT, IQ, QZ, ZQ, TSFC, NOSTAB, TOL, UW, WT, WQ,
     X    PSIM, PSIT, PHIM, PHIT, 
     X    CD, CDN, CT, CTN, CE, CEN, Z0M, Z0T,
     X    Z0Q, QSFC)
C
C  Compute infrared radiation
C
        SIG = 5.67E-8
        QZA = QZ/1000.
        RMIX = QZA/(1. - QZA) 
        E = RMIX*PRESS/(0.61 + RMIX)
        TAABS = 273.15 + TZ
        TSABS = 273.15 + TSFC
C
        EFACT = 0.39 - 0.05*SQRT(E)
        CFACT = 1. - 0.6*(CC**2)
C
        RINFRA = -0.985*SIG*EFACT*CFACT*(TSABS**4)
C
C  Compute the turbulent fluxes
C
        DENAIR = 1.23
        HCPAIR = 1005.0
        HVAP = 2.45E6
C
        STRESS = DENAIR*(-UW) 
        HFSEN = -DENAIR*HCPAIR*WT 
        HFLAT = -DENAIR*HVAP*WQ/1000.
        ROL = 10./TOL
        RLON = ALOG(ZU/Z0M)
C
        WRITE (6,50)
        WRITE (1,50) 
   50   FORMAT (1X,/,1X,'Estimated fluxes and boundary layer',
     X  ' parameters are:',/)
        WRITE (6,51) STRESS  
        WRITE (1,51) STRESS 
   51   FORMAT (1X,'Wind Stress (Pa) =           ', F10.3)
        WRITE (6,52) HFSEN
        WRITE (1,52) HFSEN
   52   FORMAT (1X,'Sensible Heat Flux (W/m**2) =', F10.0)
        WRITE (6,53) HFLAT
        WRITE (1,53) HFLAT
   53   FORMAT (1X,'Latent Heat Flux (W/m**2) =  ', F10.0)
C
        IF(IRAD.EQ.1) THEN
        WRITE (6,36) SRAD
        WRITE (1,36) SRAD
   36   FORMAT (1X,'Solar Insolation (W/m**2) =  ', F10.0)
        WRITE (6,104) SRADAV
        WRITE (1,104) SRADAV
 104    FORMAT (1X,'Daily Average SI (W/m**2) =  ', F10.0)
        WRITE (6,37) RINFRA
        WRITE (1,37) RINFRA
   37   FORMAT (1X,'Infrared Heat Flux (W/m**2) =', F10.0)
C
        HSUMA = RINFRA + SRADAV + HFLAT + HFSEN
        HSUM = RINFRA + SRAD + HFLAT + HFSEN
        ALBEDO = 0.07
        HNET = (1. - ALBEDO)*SRADAV + RINFRA + HFLAT + HFSEN
C
        WRITE (6,248) HSUMA
        WRITE (1,248) HSUMA
 248    FORMAT (1X,'Daily Avg. Net H Flux (W/m**2) = ', F6.0)
C
        WRITE (6,249) HNET
        WRITE (1,249) HNET
 249    FORMAT (1X,'D A Net H F Absorbed (W/m**2) = ', F7.0)
C
        END IF

C
        WRITE (6,38) 
        WRITE (1,38)
  38    FORMAT (1X)
C
C
        WRITE (6,56) QZ
        WRITE (1,56) QZ
   56   FORMAT (1X,'Air humidity (g/kg) =       ',F10.2)
        WRITE (6,63) QSFC
        WRITE (1,63) QSFC
   63   FORMAT (1X,'Surface air humid. (g/kg) = ',F10.2)
        WRITE (6,54) CD
        WRITE (1,54) CD
   54   FORMAT (1X,'10 m drag coefficient =     ',E10.3)
        WRITE (6,59) CDN
        WRITE (1,59) CDN
   59   FORMAT (1X,'Neutral 10 m drag coeff. =  ',E10.3)
        WRITE (6,57) Z0M
        WRITE (1,57) Z0M
   57   FORMAT (1X,'Roughness length, Z0 (m) =  ',E10.2)
        WRITE (6,44) CT
        WRITE (1,44) CT
   44   FORMAT (1X,'10 m Stanton number =       ',E10.3)
        WRITE (6,45) CTN
        WRITE (1,45) CTN
   45   FORMAT (1X,'Neutral 10 m Stanton no. =  ',E10.3)
        WRITE (6,46) CE
        WRITE (1,46) CE
   46   FORMAT (1X,'10 m Dalton number =        ',E10.3)
        WRITE (6,47) CEN
        WRITE (1,47) CEN
   47   FORMAT (1X,'Neutral 10 m Dalton no. =   ',E10.3)
        WRITE (6,60) RLON
        WRITE (1,60) RLON
   60   FORMAT (1X,'Lon(Z/Z0) =                 ',E10.3)
        WRITE (6,55) ROL
        WRITE (1,55) ROL
   55   FORMAT (1X,'Obukhov length (m) =        ',F10.2)
        WRITE (6,58) TOL
        WRITE (1,58) TOL
   58   FORMAT (1X,'Stability parameter, Z/L =  ',F10.2)
        WRITE (6,61) PSIM
        WRITE (1,61) PSIM
   61   FORMAT (1X,'Stability function, PSIM =  ',E10.3)
        WRITE (6,66) PHIM
        WRITE (1,66) PHIM
   66   FORMAT (1X,'Shear form, PHIM =          ',E10.3)
        WRITE (6,62) PSIT
        WRITE (1,62) PSIT
   62   FORMAT (1X,'Stability function, PSIT =  ',E10.3)
        WRITE (6,67) PHIT
        WRITE (1,67) PHIT
   67   FORMAT (1X,'Shear form, PHIT =          ',E10.3,/)
C
C     
        WRITE (6,70) 
   70   FORMAT (1X,'Do you want to list the profiles ? (1=yes)  ', $)
        READ (5,920) ILIST
C
        IF(ILIST.EQ.1) THEN
C
        WRITE (6,75) 
   75   FORMAT (1X,'Do you want to store the profiles ? (1=yes)  ', $)
        READ (5,920) ISTOR
        IF(ISTOR.EQ.1) OPEN (7, FILE='fluxa4.dat')
C
        VC = 0.4
        USTAR = SQRT(-UW)
        ZZ = 0.
        DZ = 0.1
C       
        DO 71 J=1,22
        ZZ = ZZ + DZ
        IF(ZZ.GT.1.0) DZ = 1.0
        IF(ZZ.GT.9.9) DZ = 10.
C
        RLZZ = ALOG(ZZ/Z0M)
        ZOL = ZZ/ROL
        RLZ = ALOG(ZZ)
        CALL FZOL (ZOL,PHIM,PHIT,PSIM,PSIT)
        UZZ = (USTAR/VC)*(ALOG(ZZ/Z0M) - PSIM)
        TZZ = TSFC - (WT/(VC*USTAR))*(ALOG(ZZ/Z0T) - PSIT)
        QZZ = QSFC - (WQ/(VC*USTAR))*(ALOG(ZZ/Z0Q) - PSIT)
C
        AAA = 0.
        IF(J.EQ.1) WRITE (1,134) 
        IF(J.EQ.1) WRITE (6,134)
 134    FORMAT (1X,6X,'Z      LN(Z)  LN(Z/Z0)  U      T       Q')
        IF(J.EQ.1) WRITE (6,73) AAA, AAA, AAA, AAA, TSFC, QSFC
        IF(J.EQ.1) WRITE (1,73) AAA, AAA, AAA, AAA, TSFC, QSFC
        IF(J.EQ.1.AND.ISTOR.EQ.1)
     X     WRITE (7,76) AAA, AAA, AAA, AAA, TSFC, QSFC
   76   FORMAT (6F10.2)
   73   FORMAT (1X,3F8.1,3F8.2)
        WRITE (6,72) ZZ, RLZ, RLZZ, UZZ, TZZ, QZZ
        WRITE (1,72) ZZ, RLZ, RLZZ, UZZ, TZZ, QZZ
        IF(ISTOR.EQ.1) WRITE (7,76) ZZ, RLZ, RLZZ, UZZ, TZZ, QZZ
   72   FORMAT (1X,3F8.1,3F8.2)
   71   CONTINUE
        END IF
C
        WRITE (6,910)
  910   FORMAT (1X,/,
     X ' Do you want another pass through ? (1=yes) ',$)
        READ (5,920) IEND
  920   FORMAT (I1)
        IF (IEND.EQ.1) GO TO 1000
        END
C
      SUBROUTINE
     x  BULK(UZ,ZU,TZ,ZT,IQ,QZ,ZQ,TSFC,NOSTAB,TOL,UW,WT,WQ,
     x  PSIM,PSIT,PHIM, PHIT,
     x  CD,CDN,CT,CTN,CE,CEN,Z0M,Z0T,Z0Q,QSFC)
C
C
CMNT  BULK ESTIMATION OF FLUXES.  FROM BILL LARGE, NCAR, OCT 84
CMNT  ESTIMATES:
CMNT   1. TOL, THE STABILITY PARAMETER Z/L AT Z=10 METERS
CMNT   2.  UW, THE KINEMATIC MOMENTUM FLUX (M/S)**2
CMNT   3.  WT, THE KINEMATIC SENSIBLE HEAT FLUX (DEGREES M/S)
CMNT   4.  WQ, THE KINEMATIC LATENT HEAT (MOISTURE) FLUX (M/S G/M**3)
CMNT   5.  CD, THE DRAG COEFFICIENT FOR THE GIVEN STABILITY, HEIGHT
C      6.  CDN, THE NEUTRAL DRAG COEFFICIENT
C      6a. CT, CE SAME AS ABOVE FOR STANTON, DALTON NOS.
C      7.  PSIM, THE STABILITY FUNCTION FOR MOMENTUM
C      8.  PSIT, THE STABILITY FUNCTION FOR HEAT AND MOISTURE
C      8a. PHIMT, PHIT, ARE SHEAR FORMS OF THE STABILITY FUNCTIONS
C      9.  Z0M, THE ROUGHNESS LENGTH FROM THE NEUTRAL DRAG COEFFICIENT
C     10.  Z0T, ROUGHNESS FOR TEMPERATURE, (M)
C     11.  Z0Q, ROUGHNESS FOR HUMIDITY, (M)
C     12.  QZ, THE SPECIFIC HUMIDITY (G/KG)
C     13.  QSFC, SATRUARTION SPECIFIC HUMIDITY AT T = TSFC (G/KG).
CMNT
CMNT  FROM THESE INPUTS:
CMNT   1.  UZ, THE MEAN WIND SPEED AT HEIGHT ZU
CMNT   2.  TZ, THE MEAN AIR TEMPERATURE (CELSIUS) AT HEIGHT ZT
CMNT   3.  QZ, THE HUMIDITY VARIABLE (RELATIVE HUMIDITY, SPECIFIC
CMNT          HUMIDITY, DEW POINT OR WET BULB TEMPERATURE) 
CMNT   3a. IQ, A FLAG THAT SPECIFIES WHICH HUMIDITY VARIABLE
CMNT   4.  TSFC, THE SEA SURFACE TEMPERATURE (CELSIUS)
C      5.  NOSTAB, A FLAG SET TO 1 IF NEUTRAL CD IS TO BE USED
C           TO ESTIMATE STRESS.
CMNT  
CMNT  OPTIONS:
CMNT   1. IF SET ZQ=0., REL HUM ASSUMED TO BE 75 PERCENT
CMNT   2. IF SET ZT=0., UW IS FOUND ASSUMING NEUTRAL STABILITY
CMNT   3. IF SET ZU=0., 10 M STRESS IS RETURNED
CMNT
CMNT  ASSUME:
CMNT   1. A NEUTRAL 10 M DRAG COEF CDN=0.00115 U10<UCH=11M/S
CMNT                  =0.00049 + 0.000065 U10  U10>UCH=11M/S
CMNT   2. A NEUTRAL 10 M STANTON NUMBER CTN=0.00115 Z/L<0
CMNT                                       =0.00075 Z/L>0
CMNT   3. A NEUTRAL 10 M DALTON NUMBER  CEN=0.00115 ALWAYS
CMNT   4. VON KARMAN'S CONSTANT, VONK=0.40
CMNT   5. THE SMALL WT REGRESSION, WTR BELOW IF BIG=FALSE
CMNT      WITH (A=.0008 IF Z/L>0, A=.0010 IF Z/L<0)
CMNT
      LOGICAL BIG
      WTR(A,U,DTH)=0.002+A*U*DTH
CMNT
CMNT  6. THE SATURATION HUMIDITY OF AIR AT T CELSIUS QSAT(T)
CMNT
      QSAT(T)=640380000./EXP(5107.4/(273.16+T))
CMNT
      DATA CDA,CTS,CTU,CEA,UCH/0.00115,0.00075,0.00115,0.00115,10./
      DATA RHDEFLT,VONK/75.,0.40/
CMNT                         ARE TEMPS RELIABLE??
      IF (ZT .GT. 0.2) GO TO 1
      TOL=0.0
      GO TO 3
CMNT  YES, CALCULATE A BULK STABILITY PARAMETER AT 10 METERS FROM
CMNT    DELTH,SURFACE - AIR POTENTIAL TEMPERATURE DIFFERENCE
CMNT    QSFC, SURFACE ABSOLUTE HUMIDITY
CMNT    TO, LOCAL MEAN VIRTUAL TEMERATURE
1     DELTH=TSFC-TZ-0.01*ZT
CMNT  THIS IS THE SEA MINUS AIR POT. TEMP. DIFF
      IF (ZQ .LE. 0.0) RH=RHDEFLT
C
C  DETERMINE SPECIFIC HUMIDITY IF IQ.NE.2
C
      IF(IQ.EQ.1) THEN
      RH = QZ
      QZ = RH*QSAT(TZ)*.01
      END IF
      IF(IQ.EQ.3) THEN
      TDEW = QZ
      QZ = QSAT(TDEW)
      END IF
      IF(IQ.EQ.4) THEN
      TWET = QZ
      CALL HUMLOW(TZ, TWET, QZ) 
      QZ = QZ*1000.
      END IF
C
C      WRITE (6,889) QZ
  889 FORMAT (1X,'SPECIFIC HUMIDITY IS', F10.4)
C
      QSFC = QSAT(TSFC)
      DELQ=0.98*QSAT(TSFC)-QZ
      TKELV=273.16 + TZ
      F1=0.00000172
      T0=TKELV+QZ*F1*TKELV**2
CMNT  
CMNT  IS /WT/ LARGE ENOUGH TO USE TRANSFER COEFFICIENTS?
      UDT=UZ*DELTH
      IF ((UDT .LT. 10.) .AND. (UDT. GT. -15.)) GO TO 2
CMNT  YES
      BIG=.TRUE.
      F2=1.00
      IF (DELTH .LT. 0.0) F2=0.70
      TOL=0.0-1000.*F2/T0/UZ**2*(DELTH+F1/F2*DELQ*T0**2)
      GO TO 3
CMNT  
CMNT  NO,     DEAL WITH SMALL SENSIBLE HEAT FLUXES
C
2     BIG=.FALSE.
      A=0.0010
      IF (DELTH .LT. 0.0) A=0.0008
      WT=WTR(A,UZ,DELTH)
      WQ=CEA*UZ*DELQ
      U3=(SQRT(CDA)*UZ)**3
      TOL=0.0-39./U3/T0*(WT+F1*WQ*T0**2)
CMNT
CMNT     RETURN WITH ONLY TOL IF ZU NOT GREATER THAN 0.2
3     IF (ZU .LE. 0.2) RETURN
      ZUOL=TOL*ZU/10.0
      ZTOL=TOL*ZT/10.0
      ZQOL=TOL*ZQ/10.0
CMNT     EVALUATE THE STABILITY FUNCTIONS IN FZOL
      TTOL = TOL
      CALL FZOL(TTOL,P10,PG10,S10,SG10)
      TOL = TTOL
      CALL FZOL(ZUOL,PHIM,PHI,PSIM,PSI)
      CALL FZOL(ZTOL,PHI,PHIT,PSIMT,PSIT)
      CALL FZOL(ZQOL,PHI,PHIQ,PSIMQ,PSIQ)
CMNT
CMNT     SOLVE FOR U10 AND CDN ITERATIVELY
C 
      STAB = 1.0
      IF(NOSTAB.EQ.1) STAB = 0.
C
      RK=(ALOG(ZU/10.0) - STAB*PSIM + S10)/VONK
      ITER=1
      UP=UZ/(1.0+SQRT(CDA)*RK)
22    CP=CDA
      IF (UP .GT. UCH) CP=(0.49+0.065*UP)/1000.
      U10=UZ/(1.0+SQRT(CP)*RK)
      DIFF=(UP-U10)/UP
      UP=U10
      IF (DIFF .LE. 0.01) GO TO 4
      ITER=ITER+1
      IF (ITER .LE. 5) GO TO 22
      WRITE(6,66)ITER,U10,DIFF,CP
66    FORMAT('-',5X,'TOO BAD, FAILED TO CONVERGE IN 5 ITERATIONS',I6,
     13F16.4)
      STOP
CMNT
CMNT      SOLVED FOR U10, NOW FIND CDN AND CD
CMNT      THEN GET UW FROM BULK FORMULA
C
4     CDN=CDA
      IF (U10 .GT. UCH) CDN=(0.49+0.065*U10)/1000.
      SQCD=SQRT(CDN)
      RCD=1.0+SQCD/VONK*(ALOG(ZU/10.0)- STAB*PSIM)
      CD=CDN/RCD**2
      UW=0.0-CD*UZ**2
C
C  COMPUTE ROUGHNESS LENGTH FROM NEUTRAL CD
C
      Z0M = 10./EXP((SQRT(.16/CDN)))
C

CMNT
CMNT  SENSIBLE HEAT FLUX IF ZT GREATER THAN 0.02
      IF (ZT .LE. 0.02) GO TO 5
      CTN=CTU
      IF (DELTH .LT. 0.0) CTN=CTS
      RCDT=1.0+SQCD/VONK*(ALOG(ZT/10.0)-STAB*PSIMT)
      CT=CTN/RCDT/(1.0+CTN/VONK/SQCD*(ALOG(ZT/10.)-STAB*PSIT))
      WT=CT*UZ*DELTH
      IF (BIG) GO TO 5
CMNT       CALCULATE SMALL SENSIBLE HEAT FLUX
      RKT=(ALOG(ZT/10.0)-STAB*PSIT+SG10)/VONK
      DTH10=DELTH-CT/SQRT(CD)*RKT*(DELTH)
      WT=WTR(A,U10,DTH10)
C

CMNT
CMNT      LATENT HEAT (MOISTURE) FLUX IF ZQ SET > 0.02
5     IF (ZQ .GT. 0.02) GO TO 6
CMNT  MOD. ORIGINAL RETURNED IF ZQ .LE. 0.02 WITHOUT
CMNT  CALCULATING A NEW WQ
      WQ=CEA*UZ*DELQ
      RETURN
6     CEN=CEA
      RCDQ=1.0+SQCD/VONK*(ALOG(ZQ/10.0)-STAB*PSIMQ)
      CE=CEN/RCDQ/(1.0+CEN/VONK/SQCD*(ALOG(ZQ/10.)-STAB*PSIQ))
      WQ=CE*UZ*DELQ
C
      Z0T = 10./EXP((SQRT(.16/CTN)))
      Z0Q = 10./EXP((SQRT(.16/CEN)))
C
C      WRITE (6,88) CDN,CTN,CEN, Z0M, Z0T, Z0Q
C   88 FORMAT (1X,'CDN,CTN,CEN, Z0, Z0T, Z0Q ARE ',
C     X //,1X,  6E10.3)
C

      RETURN
      END
      SUBROUTINE FZOL(ZOL,PHIM,PHIG,PSIM,PSIG)
C FZOL:    SUBROUTINE TO ESTIMATE THE STABILITY FUNCTIONS
CMNT       USING THE MOST UP TO DATE FLUX-PROFILE RELATIONSHIPS
CMNT
C
      IF (ZOL .GE. 0.0) GO TO 2
C
CMNT             UNSTABLE
      X=(1.0-16.0*ZOL)**0.25
      PHIM=1.0/X
      PHIG=PHIM**2
      PSIG=2.0*ALOG(0.5+0.5*X**2)
      PSIM=0.5*PSIG+2.0*ALOG(0.5*(1.0+X))-2.0*ATAN(X)+1.571
      RETURN
CMNT
CMNT     STABLE
C
2     PHIG=1.0+7.0*ZOL
      PHIM=1.0+7.0*ZOL
      PSIG=-7.0*ZOL
      PSIM=-7.0*ZOL
      RETURN
      END          

        SUBROUTINE HUMLOW(T,TW,Q)
C
C TO EVALUATE SPECIFIC HUMIDITY Q FROM DRY AND WET BULB TEMP
C T AND TW ACCORDING TO LOWE(77) JAM 16 100-103
C WRITTEN BY TIM LIU ON 5/3/79, REVISED FOR VAX ON 2/10/82
C
        DIMENSION A(6)
        DATA A/4.436519E-1,1.428946E-2,2.650649E-4,3.031240E-6,
     $  2.034081E-8,6.136821E-11/
        P=1013.25
        X=0.
        DO 100 I=1,6
        J=7-I
        X=(X+A(J))*TW
 100    CONTINUE
        ES=6.107800+X
        Q=0.622*ES/(P-ES)-4.045E-04*(T-TW)
        RETURN
        END
C
      SUBROUTINE SOLAR(YDAY,RLAT,ATC,SRAD,SRADMX,SRADAV,SUMR,SUMTIM)
C
C
C  This subroutine computes solar insolation from:
C  YDAY, the year day, i.e., 10.5 is Jan 10th at 1200 L
C  RLAT, the latitude (deg)
C  ATC, the atmospheric transmission coefficient, range of
C    1 to 0, nominal clear air is 0.8.
C
C  IALLDY, a flag, if = 0, do data for the exact time YDAY only,
C                  if = 1, compute data for entire day. 
C  and returns:
C
C  SRAD, insolation at the time YDAY, (W/m**2),
C  SRADMX, the maximum insolation on YDAY, (W/m**2),
C  SRADAV, the daily average insolation on YDAY, (W/m**2), 
C  SUMR, the total insolation on that day, (J/m**2),
C  SUMTIM, the duration of insolation, (sec).
C
C  Written by Jim Price, February, 1991.
C
C
      R0 = 1370.
C
      SUMR = 0.
      RMAX = 0.
      SUMTIM = 0.
C
C
      DT = .001*8.64E4
      NSTEP = 1000
      DO 1 J=1,NSTEP
      T = YDAY + .001*FLOAT(J-1)
C
      SA = SINA(T,RLAT)
C
      IF(SA.GT.0.) THEN
      RSS = R0*SA*ATC**(1./SA)
      RDIFF = R0*SA*(0.91 - ATC**(1./SA))/2.
      RSS = RSS + RDIFF
      SUMR = SUMR + DT*RSS
      SUMTIM = SUMTIM + DT
      IF(RSS.GT.RMAX) RMAX = RSS
C
      END IF
C
      IF(J.EQ.1) SRAD = RSS
C
    1 CONTINUE
C
      SRADMX = RMAX
      SRADAV = SUMR/8.64E4
C
      RETURN 
      END
C
      FUNCTION SINA(TIME,RLAT)
      TP = 3.14159*2.0
      D2R = TP/360.
      YRPH = (TP/365.)*(357.0 - TIME)
      THTA = -23.0*COS(YRPH)*D2R
C
C  SINA IS THE SINE OF SUN'S ALTITUDE
C
      SINA = COS(THTA)*COS(RLAT*D2R)*COS((TIME + 0.5)*TP) + 
     C  SIN(THTA)*SIN(RLAT*D2R)
      RETURN
      END
