COMPUTER CODES FOR EQUIVALENT NEUTRAL WIND CALCULATION
PROGRAM FLUX
C APPLICATION OF LKB BULK PARAMETERIZATION COMPUTER CODE.
C The bulk parameterization technique of Liu et al. (J. Atmos.
C Sci., 36, 1722-1735, 1979) essentially solves three
C simultaneous equations representing the non-dimensional wind,
C temperature, and humidity profiles by iterations. This is
C achieved by calling Subroutine Asl1, which in turn calls
C Subroutines: Humlow, Lkb1, Psi, Zeta. The input and
C output parameters (through Common statement) are listed in
C the comment statements of Asl1.
C Please note that
C [1]Theoretically, the validity of the Liu et al. model depends on
C the validity of the similarity theory. For example, when
C turbulence is suppressed by stable density stratification
C (bulk Richardson number exceeds a critical value), the
C results may not be valid and the program may fail to
C converge.
C [2]In case not all the input parameters are available,
C substituting with estimated values is better than leaving the
C missing values as zero. For example, if specific humidity is
C not available, assume a 75% or 80% relative humidity.
C [3]Subroutine Humlow can be used to derive surface specific
C humidity Q from three types of observations, and the
C following are examples
C (1) If dew point temperature (TD) is available,
C Call Humlow (TD,TD,P,Q)
C (2)If air temperature (T) and relative humidity (R) are
C available
C Call Humlow (T,T,P,QA)
C Q=QA*R
C (3)If wet-bulb temperature (TW) is available
C Call Humlow (T,TW,P,Q)
C P is air pressure in mb, which can be taken as 1013 if no
C observation is available
REAL LH
COMMON/PIN/U,T,Q,TS,ZU,ZT,ZQ,P
COMMON/POUT/USR,TSR,QSR,ZO,ZL,RR,RT,RQ
P=1013. !SURFACE PRESSURE
CP=1.0E+03 !ISOBARIC SPECIFIC HEAT
RHO=1.2 !SURFACE AIR DENSITY
LH=2.5E+06 !LATENT HEAT OF VAPORIZATION
C EXAMPLE:
PRINT *,' ENTER WIND,TEMP,AND HUMIDITY SENSOR HEIGHTS'
READ(5,*)ZU,ZT,ZQ
PRINT *,' ENTER WIND,SST,TEMP,REL HUMIDITY'
READ(5,*)U,TS,T,R
CALL HUMLOW(T,T,P,QA)
Q=QA*R
print*,'q=',q
CALL ASL1(IER)
print*,'tsr,qsr,usr',tsr,qsr,usr
IF(IER .LT. 0)GO TO 90
C
C WRITE STABILITY AND ROUGHNESS LENGTH
C
WRITE(6,1)ZL
1 FORMAT(' STABILITY PARAMETER=',E12.2)
WRITE(6,2)ZO
2 FORMAT(' ROUGHNESS LENGTH=',E12.2)
C
C COMPUTE SURFACE STRESS TAU, SENSIBLE HEAT FLUX H, AND LATENT HEAT FLUX E
C
TAU=RHO*USR**2
H=-CP*RHO*USR*TSR
E=-LH*RHO*USR*QSR
WRITE(6,3)TAU
3 FORMAT(' SURFACE STRESS=',F10.3,' N/M**2')
WRITE(6,4)H
4 FORMAT(' SURFACE HEAT FLUX=',F10.3,' W/M**2')
WRITE(6,5)E
5 FORMAT(' SURFACE LATENT HEAT FLUX=',F10.3,' W/M**2')
C
C COMPUTE EQUIVALENT NEUTRAL WIND AT 19.5M HEIGHT
C
U20=2.5*USR*ALOG(19.5/ZO)
WRITE(6,6)U20
6 FORMAT(' EQUIVALENT NEUTRAL WIND AT 19.5M='F10.3)
C
C COMPUTE TRANSFER COEFFICIENT
C
CALL HUMLOW(TS,TS,P,QS)
CD=(USR/U)**2
CH=USR*TSR/(U*(T-TS))
CE=USR*QSR/(U*(Q-QS))
WRITE(6,7)CD,CH,CE
7 FORMAT(' THE COEFS (CD,CH,CE)=',3(1PE12.2))
GO TO 50
90 WRITE(6,8)IER
8 FORMAT(' FAILS TO CONVERGE ',I5)
50 CONTINUE
END
SUBROUTINE ASL1(IER)
C
C TO EVALUATE SURFACE FLUXES, SURFACE ROUGHNESS, AND STABILITY OF
C THE ATMOSPHERIC SURFACE LAYER FROM BULK PARAMETERS ACCORDING TO
C LIU ET AL. (79) JAS 36 1722-1735
C WRITTEN BY TIM LIU ON 5/8/79, FIRST REVISION 8/31/94
C
C INPUT:
C U WIND SPEED IN M/S
C T TEMPERATURE IN DEG C
C Q SPECIFIC HUMIDITY DIMENSIONLESS
C TS SURFACE TEMPERATURE IN DEG C
C ZU HEIGHT OF WIND SENSOR (M)
C ZT HEIGHT OF TEMPERATURE SENSOR (M)
C ZQ HEIGHT OF HUMIDITY SENSOR (M)
C PO SURFACE PRESSURE IN MB (DEFAULT TO 1013.25)
C
C DELTA = RELATIVE TOLERANCE FOR CONVERGENCE
C
C OUTPUT:
C USR,TSR,QSR SCALING QUANTITIES FOR U,T,Q
C ZO,ZL ROUGHNESS AND STABILITY PARAMETERS
C RR,RT,RQ ROUGHNESS REYNOLD NUMBERS FOR U,T,Q
C
C DISCARD OUTPUT IF IER GREATER THAN O
C IER=1 FAIL TO CONVERGE
C
COMMON/PIN/U,T,Q,TS,ZU,ZT,ZQ,PO
COMMON/POUT/USR,TSR,QSR,ZO,ZL,RR,RT,RQ
IER=0
VISA=.15E-4
g=9.8
delta=0.0001
ZL=0.
US=0.
IF(PO .EQ. 0.)PO=1013.25
CALL HUMLOW(TS,TS,PO,QS)
print*,'asl qs=',qs
DU=U-US
DT=T-TS
DQ=Q-QS
TA=273.15+T
USR=.04*DU
N=0
30 CONTINUE
zo=0.11*visa/usr + 0.011*usr*usr/g
PUZ=PSI(1,ZL)
USR=DU*0.4/(ALOG(ZU/ZO)-PUZ)
RR=ZO*USR/VISA
CALL LKB1(RR,RT,1)
CALL LKB1(RR,RQ,2)
ZTL=ZL*ZT/ZU
ZQL=ZL*ZQ/ZU
PTZ=PSI(2,ZTL)
PQZ=PSI(2,ZQL)
ZTSR=ZT*USR/VISA
ZQSR=ZQ*USR/VISA
S=2.2*(ALOG(ZTSR/RT)-PTZ)
D=2.2*(ALOG(ZQSR/RQ)-PQZ)
TSR=DT/S
QSR=DQ/D
CALL ZETA(T,Q,USR,TSR,QSR,ZU,ZLN)
TEST=ABS((ZL-ZLN)/(ZL+1.E-8))
IF(TEST.LT.delta) GO TO 39
N=N+1
IF(N.GT.50)GO TO 95
ZL=ZLN
GO TO 30
39 CONTINUE
GO TO 99
95 IER=1
WRITE(6,1)N
1 FORMAT(1X,24HASL FAILS TO CONVERGE,I5)
99 RETURN
END
SUBROUTINE HUMLOW(T,TW,P,Q)
C
C TO EVALUATE SPECIFIC HUMIDITY Q FROM DRY AND WET BULB TEMP
C T AND TW IN DEG C AND PRESSURE P IN MB
C Q IS THE SATURATION SPECIFIC HUMIDITY AT T IF TW=T
C WRITTEN BY TIM LIU ON 5/3/79
C
DIMENSION A(6)
DATA A/4.436519E-1,1.428946E-2,2.650649E-4,3.031240E-6,
& 2.034081E-8,6.136821E-11/
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 lkb1(rr,rt,iflag)
C TO DETERMINE THE LOWER BOUNDARY VALUE RT OF THE LOGARITHMIC
C PROFILES OF TEMPERATURE (IFLAG=1) OR HUMIDITY (IFLAG=2)
C IN THE ATMOSPHERE FROM ROUGHNESS REYNOLD NUMBER RR BETWEEN
C 0 AND 1000. OUT OF RANGE RR INDICATED BY RT=-999.
C BASED ON LIU ET AL. (1979) JAS 36 1722-1723
C WRITTEN BY WENDY TANG 8/31/94
C
dimension a(9,2)
data a/1.78372207e-02,-9.42581262e-02,-6.86854874e-01,
$ 1.23243995e-01, 9.60776184e-02,-3.83157945e-02,
$ -3.61932655e-03, 2.95832855e-03,-3.08498119e-04,
$ 3.19474376e-01,-8.14176320e-02,-5.96190694e-01,
$ 1.06061761e-01, 8.06995259e-02,-3.26535114e-02,
$ -2.90657805e-03, 2.47245084e-03,-2.58972350e-04/
xx=alog(rr)
xi=1.0
yy=a(1,iflag)
do i=2,9
xi=xi*xx
yy=yy+a(i,iflag)*xi
enddo
rt=exp(yy)
return
end
FUNCTION PSI(ID,ZL)
C
C TO EVALUATE THE STABILITY FUNCTION PSI FOR WIND SPEED (IFLAG=1)
C OR FOR TEMPERATURE AND HUMIDITY PROFILES FROM STABILITY PARAMETER ZL
C SEE LIU ET AL. (1979) JAS 36 1722-1723 FOR DETAILS
C WRITTEN BY TIM LIU ON 9/12/71, REVISED FOR VAX ON 2/10/82
C
IF(ZL)10,20,30
10 CHI=(1.-16.*ZL)**0.25
IF(ID.EQ.1)GO TO 11
PSI=2.*ALOG((1.+CHI*CHI)/2.)
GO TO 99
11 PSI=2.*ALOG((1.+CHI)/2.)+ALOG((1.+CHI*CHI)/2.)-2.*ATAN(CHI)
& +2.*ATAN(1.)
GO TO 99
20 PSI=0.
GO TO 99
30 PSI=-6.*ALOG(1.+ZL)
99 RETURN
END
SUBROUTINE ZETA(T,Q,USR,TSR,QSR,Z,ZL)
C
C TO EVALUATE OBUKHOV'S STABILITY PARAMETER Z/L FROM AVERAGE
C TEMP T IN DEG C, AVERAGE HUMIDITY Q IN GM/GM, HEIGHT Z IN M,
C AND FRICTIONAL VEL.,TEMP., HUM. IN MKS UNITS
C SEE LIU ET AL. (1979) JAS 36 1722-1723 FOR DETAILS
C WRITTEN BY TIM LIU ON 10/1/77, REVISED FOR VAX ON 2/10/82
C
VON=0.4
G=9.81
TA=273.16+T
TV=TA*(1.+0.61*Q)
TVSR=TSR*(1.+0.61*Q)+0.61*TA*QSR
IF(TVSR.EQ.0.)GO TO 10
OB=TV*USR*USR/(G*VON*TVSR)
ZL=Z/OB
GO TO 99
10 ZL=0.
99 RETURN
END
