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 U10UCH=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