C ********************************************************************** C * * C * SOFTWARE LICENSING * C * * C * This program is free software; you can redistribute * C * it and/or modify it under the terms of the GNU * C * General Public License as published by the Free * C * Software Foundation, either Version 2 of the * C * license, or (at your option) any later version. * C * * C * This program is distributed in the hope that it * C * will be useful, but without any warranty; without * C * even the implied warranty of merchantability or * C * fitness for a particular purpose. See the GNU * C * General Public License for more details. * C * * C * A copy of the GNU General Public License is * C * available at http://www.gnu.org/copyleft/gpl.html * C * or by writing to the Free Software Foundation, Inc.,* C * 59 Temple Place - Suite 330, Boston, MA 02111, USA. * C * * C ********************************************************************** CNG SUBROUTINE HITFLX (ALATJ,XLNG0,IYR,IMON,IDAY,IHR,PSS,TSS, SUBROUTINE HITFLX (ALATJ,XLNG0,IYR,IMON,IDAY,CHRS,PSS,TSS, + DRYBT,VMOD,ANGLE,RH,CLOUD,DHR,SWR,RLW,HNOT,WENOT,GMT) C C SUBROUTINE HITFLX (ALATJ,XLNG0,GMT,IYR,IMON,IDAY,IHR,PSS,TSS, C + DRYBT,VMOD,ANGLE,RH,CLOUD,DHR,SWR,RLW,HNOT,WENOT,HFS,DHFS, C +XKS,TEMPE,HEATXL,ESATS,ESATA,FWW,DIFEA) C C THIS SUBROUTINE CALCULATES SOLAR (SHORT WAVE) RADIATION, C REFLECTED (LONG WAVE) RADIATION, SENSIBLE HEAT FLUX, C LATENT HEAT FLUX, TOTAL HEAT FLUX AND PARTIAL DERIVATIVE, C D (HEAT FLUX) /DT, USED FOR LINEAR CORRECTION OF THE TOTAL C HEAT FLUX BASED ON SOME VALUE OF THE SEA SURFACE TEMPERATURE C DIFFERENT FROM THE GIVEN VALUE. ALL UNITS ARE IN SI. C C INPUT DATA: C C ALATJ = LATITUDE, IN DEGREES C XLNG0 = LONGITUDE, IN DEGREES WEST C GMT = SHIFT IN THE START TIME VERSUS GMT C (FOR INSTANCE, IF LST FOR DELAWARE BAY LOCATION GMT=-5.0) c change IYR to 4-digit number (hli, 04/07/2000) C IYR = YEAR (ONLY TWO LAST DIGITS, SUCH AS 84 FOR 1984) c end change C IMON = MONTH (1 FOR JANUARY, 2 FOR FEBRUARY ETC.) C IDAY = JULIAN DAY C IHR = HOUR (0 FOR MIDNIGHT, 12 FOR MIDDAY, ETC.) C PSS = ATMOSPHERIC PRESSURE AT THE SEA SURFACE LEVEL C (NEWTON/M**2; 1000 MB WILL BE EQUAL TO 1.0E5) C TSS = SEA SURFACE TEMPERATURE (DEG. C) C DRYBT = DRY BULB TEMPERATURE (DEG. C) C VMOD = AMPLITUDE OF THE WIND VELOCITY (M/SEC) C ANGLE = WIND DIRECTION (ANGLE) (IN DEGREES) C RH = RELATIVE HUMIDITY (BETWEEN 0.0 AND 1.0) C CLOUD = CLOUDINESS (BETWEEN 0.0 AND 1.0) C DHR = PERIOD OF AVERAGING FOR SOLAR RADIATION (HOURS) C C FOLLOWINGS ARE INTRODUCED BY QUAMRUL AHSAN 10/16/92 C C XKS = TRANSFER COEFFICIENT C TEMPE = EQUIVALENT TEMPERATURE C HEATXL= LINEARIZED HEAT FLUX COMPUTED FROM TRANFER C COEFFICIENT, EQUIVALENT TEMPERATURE AND ATMOSPHEREIC C TEMPERATURE. HEATXL=-XKS*(TEMPE-DRYBTP) C C OUTPUT DATA: C C SWR = SHORT WAVE RADIATION INCLUDING CLOUDINESS (WATT/M**2) C RLW = REFLECTED LONG WAVE RADIATION (WATT/M**2) C HNOT = SENSIBLE HEAT FLUX (WATT/M**2) C WENOT = LATENT HEAT FLUX (WATT/M**2) C HFS = TOTAL HEAT FLUX (WATT/M**2) C DHFS = D (HEAT FLUX) / DT C C PARAMETER (KB=11,KBM1=KB-1,JMT=1) DIMENSION ATTEN(KB),ATTENZ(KB),ATTZ(KB),MON(12),AMON(12) DIMENSION ALB1(20),ZA(20),DZA(19) DATA ALB1/.719,.656,.603,.480,.385,.300,.250,.193,.164,.131 * ,.103,.084,.071,.061,.054,.039,.036,.032,.031,.030/ DATA ZA/90.,88.,86.,84.,82.,80.,78.,76.,74.,70.,66.,62. * ,58.,54.,50.,40.,30.,20.,10.,0.0/ DATA DZA/8*2.0,6*4.0,5*10.0/ DATA P62/.62/,P19/.0019/,P39/.39/,P05/.05/,P8/.8/ DATA P6/.6/,P006/.006/ DATA CH,CE/0.0015,0.0015/ DATA AMON/0,31,60,91,121,152,182,213,244,274,305,335/ DATA MON/0,31,60,91,121,152,182,213,244,274,305,335/ C C SET UP MON ARRAY BASED ON WHETHER OR NOT THE CURRENT YEAR C IS A LEAP YEAR C IF ( MOD(IYR,4) .NE. 0 ) THEN DO 5 I=3,12 MON(I)=AMON(I)-1 5 CONTINUE END IF C C5D9=5.0/9.0 RADIAN=180.0/3.141593 ALATJZ=ALATJ/RADIAN TREF=273.16 RTREF=1.0/TREF CP=1005.0 RD=2.8708E2 RV=4.6157E2 W=2.4957E6 WCP=W/CP WRV=W/RV POW=RD/CP DW=1.0 SC=1353.0 SIGSB=5.67E-8 EMISS=.97 CPCH = CP*CH CEDW = CE*DW C TSSP=TSS DRYBTP=DRYBT TSS=TSS+273.16 DRYBT=DRYBT+273.16 C**** C**** DOWNWARD IRRADIANCE AFTER PAULSON AND SIMPSON (JPO 1977) C**** (SEE ALSO SIMPSON AND DICKEY - JPO 1981) C**** C**** TYPE I C GAM1=-(1.0/150.0) C GAM2=-(1.0/2000.0) C R1=0.62 C**** TYPE III C GAM1=-(1.0/140.0) C GAM2=-(1.0/790.0) C R1=0.78 C R11=1.0-R1 DO 10 K=1,KB ATTEN(K)=0.0 ATTENZ(K)=0.0 10 CONTINUE C ATTEN(*)=R1*EXP(Z(*)*GAM1)+R11*EXP(Z(*)*GAM2) C**** TAKE Z DERIVATIVE TO CALC THE DIVERGENCE OF DOWNWARD IRRADIANCE C**** MULT BY IRRADIANCE AT SURFACE C ATTENZ(1)=(1.0-ATTEN(1))/DZ(1) C DO 1488 K=2,KB C ATTENZ(K)=(ATTEN(K-1)-ATTEN(K))/DZ(K) C1488 CONTINUE C PRINT 1489,ATTENZ C1489 FORMAT(2X,'ATTENUATION',(1X,10E12.4)) C**** C**** C**** SOLAR CALCULATES THE COS OF THE SUN'S ZENITH ANGLE, C**** THE EARTH'S RADIUS VECTOR,AND SOLAR ALTITUDE FOR C**** EACH LATITUDE. C**** C**** R = EARTH'S RADIUS VECTOR C**** TAUDA = TOTAL DAYLIGHT FRACTION C**** COSZEN = COSINE OF ZENITH ANGLE C**** SOLALT = SOLAR ALTITUDE C**** JD = MON(IMON) + IDAY C C JD IS DESIGNED TO START AT JAN 1, 1984 00:00 AM C JD = JD + 2415020 + 25566 + 5115 C C THIS IF-THEN LOOP RESETS THE STARTING JULIAN DAY TO C JAN. 1 00:00 AM OF THE CURRENT YEAR C c change 2-digit year to 4-digit (hli, 04/07/2000) IF (IYR.LT.1984) THEN IDYR=1984-IYR ILEAP=IDYR/4 NDAY=IDYR*365+ILEAP JD=JD-NDAY ELSE IF (IYR.GT.1984) THEN IDYR=IYR-1984 ILEAP=IDYR/4+1 NDAY=IDYR*365+ILEAP JD=JD+NDAY END IF c end change C CNG FJD = FLOAT(IHR)/24.0 FJD = AMOD(CHRS/24,1.0) ! FJD IS HRS AFTER START OF THAT CALENDAR DAY C CALL SOLAR(JD,FJD,R,DLT,ALP,SLAG,1,HANG,TAUDA,COSZEN,SOLALT, 1 ALATJZ) CALL ZENITH(FJD,DLT,SLAG,ALATJZ,HANG,DHR,COSZRS,FRAC,XLNG0, 1 GMT) TAUDA = FRAC COSZEN= COSZRS RR2=R*R SWDTOP=(SC/RR2)*COSZEN*TAUDA IF (SWDTOP.EQ.0.0) THEN SWD=0.0 SWR=SWD GO TO 73 END IF C**** C**** CALCULATE INSOLATION AT OCEAN SURFACE C**** FXA=1.0 C**** TRANSMISSION COEFF= 0.68 AFTER HSUING JGR'86 TRANS=0.68 TSEC=TRANS**(FXA/COSZEN) C**** CALCULATE DIRECT SURFACE INSOLATION (ASSUME TRANS COEFF=.68) SWDSUR=TSEC*SWDTOP C**** CALCULATE DIRECT+DIFFUSE (AFTER SMITHSONIAN APPROX.) ABSOR=0.91 SWDSUR=SWDSUR+(ABSOR*SWDTOP-SWDSUR)*0.5 C**** C**** CALCULATE TOTAL SWD C**** CALCULATE ALBEDO AFTER PAYNE JAS 72 ZEN=RADIAN*ACOS(COSZEN) IF(ZEN.LT.74.0) GO TO 71 JAB=0.5*(90.0-ZEN)+1.0 GO TO 75 71 IF(ZEN.LT.50.0) GO TO 72 JAB=0.23*(74.0-ZEN)+9.0 GO TO 75 72 JAB=0.1*(50.0-ZEN)+15.0 75 CONTINUE DZEN=-(ZEN-ZA(JAB))/DZA(JAB) ALBEDO=ALB1(JAB)+DZEN*(ALB1(JAB+1)-ALB1(JAB)) SWDSUR=SWDSUR*(1.0-ALBEDO) C**** TAKE CLOUDY CONDITIONS INTO ACCOUNT (AFTER REED JPO 77) SWD=SWDSUR*(1.0-P62*CLOUD+P19*SOLALT) SWR=SWD 73 CONTINUE C C CALCULATE SATURATION MIXING RATIO ESTS=611.0*EXP(WRV*(RTREF-1.0/TSS)) ESATS=ESTS QSAT=RD/RV*ESTS/(PSS-ESTS) C**** CALC SAT VAPOR PRESS FOR DRY BULB TEMPERATURE DRYBT ESTA=611.0*EXP(WRV*(RTREF-1.0/DRYBT)) ESATA= ESTA QASAT=RD/RV*ESTA/(PSS-ESTA) C**** CALC MIXING RATIO FROM SAT VAPOR PRESS AND RH QKM=RH*QASAT C CALCULATE VIRTUAL TEMPERATURE AND AIR DENSITY W1=1.0+QKM W2=RD/RV*W1 W3=RD/RV+QKM W4=W3/W2 VTEMP=DRYBT*W4 W1=RD*VTEMP RHOAIR=PSS/W1 RHOA=PSS/(RD*DRYBT*1.025E3) C C COMPUTE SENSIBLE HEAT FLUX C W1 = RHOAIR * CPCH W2 = W1 * VMOD W3=TSS-DRYBT HNOT = W2 * W3 C C COMPUTE LATENT HEAT FLUX C W1 = RHOAIR * CEDW W2 = W1 * VMOD W3=QSAT-QKM ENOT = W2 * W3 WENOT=W*ENOT C Note by Quamrul ******* 1/19/96 C C WENOT is in w/m2 C PSS, QSAT and QKM are in pascal (N/m2), see p-276 of Boundary Layer Met C we use here 611. not .611 in ESTS formula above thats why in pascal C 1 mmHg = 1.33333 mbar C So to get F(w) in w/m2 / (mmHg) as follows C FWW = WENOT/W3/(PSS/0.622) ! f(w) function Lat Heat = f(w)*(Esat -Ea) FWW = FWW*10**2*1.33333 ! N/m2 to mb to mmHg C DIFEA = W3*(PSS/0.622)/1.0E2/1.3333 ! ESTS - ESTA; diff in vapor press OR as follows DIFEA = (ESTS-RH*ESTA)/1.0E2/1.3333 ! ESTS - ESTA; diff in vapor press C C NET LONG WAVE (AFTER BERLIAND:SEE SIMPSON-PAULSON QJRMS 79) C W1 = DRYBT**3 W2=W1*DRYBT W3=EMISS*SIGSB*W2*(P39-P05*SQRT(ESTA*1.0E-2)) W4=EMISS*SIGSB*4.0*W1*(TSS-DRYBT) C**** TAKE CLOUDINESS INTO ACCOUNT RLW =W3*(1.0-(P8*CLOUD))+W4 AR=SWD SWD=SWD*(1.0-ATTEN(1)) C**** C**** COMPUTE TOTAL HEAT FLUX C**** HFS=SWD-RLW-HNOT-WENOT C C CALCULATE PARTIAL DERIVATIVE DQ/DT: ( DELTA(T)=3 DEG.) C TSS=TSS+3.0 C C CALCULATE SATURATION MIXING RATIO ESTS=611.0*EXP(WRV*(RTREF-1.0/TSS)) QSAT=RD/RV*ESTS/(PSS-ESTS) C C COMPUTE SENSIBLE HEAT FLUX C W1 = RHOAIR * CPCH W2 = W1 * VMOD W3=TSS-DRYBT HNOT1 = W2 * W3 C C COMPUTE LATENT HEAT FLUX C W1 = RHOAIR * CEDW W2 = W1 * VMOD W3=QSAT-QKM ENOT1 = W2 * W3 C C NET LONG WAVE (AFTER BERLIAND:SEE SIMPSON-PAULSON QJRMS 79) C W1 = DRYBT**3 W2=W1*DRYBT W3=EMISS*SIGSB*W2*(P39-P05*SQRT(ESTA*1.0E-2)) W4=EMISS*SIGSB*4.0*W1*(TSS-DRYBT) C**** TAKE CLOUDINESS INTO ACCOUNT RLW1 =W3*(1.0-(P8*CLOUD))+W4 AR=SWD SWD=SWD*(1.0-ATTEN(1)) HFS1=SWD-RLW1-HNOT1-W*ENOT1 DHFS=(HFS1-HFS)/3.0 DRYBT=DRYBT-273.16 TSS=TSS-273.16-3.0 C WRITE (6,80) C WRITE (6,90) IYR,IMON,IDAY,IHR,VMOD,ANGLE, C 1 DRYBT,TSS,CLOUD,PSS,HFS,DHFS,SWR 80 FORMAT(//5X,'INPUT AND OUTPUT INFORMATION'// 1 ' YR ',' MON ',' DAY ',' HOUR',' VELOCITY ', 2 ' ANGLE(T) ',' TEMPER. ',' SST ', 3 ' CLOUD ',' SURF PRESS.',' HEAT FLX ',' DQ/DT ', 4 ' SUN RAD. ') 90 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,2X,1PE10.3,2X,E10.3,2X, 1 E10.3,2X,E10.3,2X,E10.3,2X,E10.3, 2 2X,E10.3,2X,E10.3,2X,E10.3) C**** C**** INCLUDE DIVERGENCE OF DOWNWARD IRRADIANCE C**** C W3=AR*C2DTTS/60.0SO C DO 410 K=2,KB C TB(*,*,K)=TB(*,*,K)+(W3*ATTENZ(K))*FSM(*,*) C 410 CONTINUE C C ************************************************************* C C CALCULATES THE EQUILIBRIUM TEMP (TEMPE) AND TRANSFER COEFF. C (XKS), SETTING NET FLUX=0; AND TS=TE C FOR THEORITICAL COMPUTATIONS PLEASE SEE ONONDAGA LAKE FILE-2 C C ************************************************************* C C TRANSFER COEFFICIENT (XKS) C DELTA=273.16 W1=4.0*DELTA**3 W2=W1*EMISS*SIGSB*(P39-P05*SQRT(ESTA*1.0E-2)) XKS1=W2*(1.0-P8*CLOUD) C W1=16.0*DELTA**3*EMISS*SIGSB W2=12.0*EMISS*SIGSB*(DRYBTP+DELTA)*DELTA**2 XKS2=W1-W2 C XKS3=RHOAIR*CPCH*VMOD C XKS=-(XKS1+XKS2+XKS3) C C C CALCULATION OF EQUILIBRIUM TEMPERATURE (TEMPE) C W1=-SWD W2=EMISS*SIGSB*DELTA**4 W3=P39-P8*SQRT(ESTA*1.0E-2) W4=W3*(1.0-P8*CLOUD) W5=W2*W3*W4 W6=-4.0*EMISS*SIGSB*DELTA**3*DRYBTP W7=-RHOAIR*CPCH*VMOD*DRYBTP C TEMPE=(W1+W5+W6+W7)/XKS C HEATXL=-XKS*(TEMPE-TSS) ! LINEARIZED HEAT FLUX FROM XKS & TTS C RETURN END