| 1 |
jmc |
1.3 |
C $Header: /u/gcmpack/MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml_coare3_flux.F,v 1.2 2013/05/25 18:19:05 jmc Exp $ |
| 2 |
jmc |
1.1 |
C $Name: $ |
| 3 |
|
|
|
| 4 |
|
|
#include "CHEAPAML_OPTIONS.h" |
| 5 |
|
|
|
| 6 |
|
|
C-- File cheapaml_coare3_flux.F: |
| 7 |
|
|
C-- Contents: |
| 8 |
|
|
C-- o CHEAPAML_COARE3_FLUX |
| 9 |
|
|
C-- o PSIU (Function) |
| 10 |
|
|
C-- o PSIT (Function) |
| 11 |
|
|
|
| 12 |
|
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 13 |
|
|
|
| 14 |
|
|
CBOP |
| 15 |
|
|
C !ROUTINE: CHEAPAML_COARE3_FLUX |
| 16 |
|
|
C !INTERFACE: |
| 17 |
|
|
SUBROUTINE CHEAPAML_COARE3_FLUX( |
| 18 |
|
|
I i,j,bi,bj, iceOrNot, |
| 19 |
jmc |
1.2 |
I tSurf, windSq, |
| 20 |
|
|
O hf, ef, evap, Rnl, ssqt, q100, cdq, cdu, |
| 21 |
jmc |
1.1 |
O dSensdTs, dEvapdTs, dLWdTs, dQAdTs, |
| 22 |
|
|
I myIter, myThid ) |
| 23 |
|
|
|
| 24 |
|
|
C !DESCRIPTION: |
| 25 |
|
|
|
| 26 |
|
|
C !USES: |
| 27 |
|
|
IMPLICIT NONE |
| 28 |
|
|
C === Global variables === |
| 29 |
|
|
#include "SIZE.h" |
| 30 |
|
|
#include "EEPARAMS.h" |
| 31 |
|
|
#include "PARAMS.h" |
| 32 |
|
|
#include "CHEAPAML.h" |
| 33 |
|
|
|
| 34 |
|
|
C !INPUT PARAMETERS: |
| 35 |
|
|
C i, j :: local indices of current grid-point |
| 36 |
|
|
C bi, bj :: current tile indices |
| 37 |
|
|
C iceOrNot :: 0=open water, 1=ice cover, 2=ice+snow |
| 38 |
|
|
C tSurf :: surface temperature |
| 39 |
jmc |
1.2 |
C windSq :: relative wind (vs surface motion) speed square |
| 40 |
jmc |
1.1 |
C myIter :: Current iteration number in simulation |
| 41 |
|
|
C myThid :: My Thread Id number |
| 42 |
|
|
INTEGER i,j,bi,bj |
| 43 |
|
|
INTEGER iceOrNot |
| 44 |
jmc |
1.2 |
_RL tSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
| 45 |
|
|
_RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
| 46 |
jmc |
1.1 |
INTEGER myIter, myThid |
| 47 |
|
|
C !OUTPUT PARAMETERS: |
| 48 |
jmc |
1.2 |
C cdu :: surface drag coeff (for wind stress) |
| 49 |
jmc |
1.1 |
_RL hf, ef, evap, Rnl, ssqt, q100, cdq |
| 50 |
jmc |
1.2 |
_RL cdu |
| 51 |
|
|
C derivative vs surf. temp of Sensible, Evap, LW, q100 |
| 52 |
|
|
_RL dSensdTs, dEvapdTs, dLWdTs, dQAdTs |
| 53 |
jmc |
1.1 |
CEOP |
| 54 |
|
|
|
| 55 |
|
|
C !LOCAL VARIABLES: |
| 56 |
|
|
INTEGER iter,nits |
| 57 |
|
|
_RL tau,L,psu,pst,Bf |
| 58 |
jmc |
1.2 |
_RL CD,usr,tsr,qsr |
| 59 |
|
|
c _RL ttas,ttt,ttt2,pt,essqt |
| 60 |
|
|
_RL zo,zot,zoq,RR,zL |
| 61 |
|
|
_RL twoPI,cwave,lwave |
| 62 |
jmc |
1.1 |
|
| 63 |
|
|
C various constants |
| 64 |
jmc |
1.2 |
_RL u,q,Tas,tta,zi,es,qs,tsw |
| 65 |
jmc |
1.1 |
_RL psiu,psit,zot10,Ct10,CC,Ribu |
| 66 |
|
|
_RL Du,Wg,Dt,Dq,u10,zo10,Cd10,Ch10 |
| 67 |
|
|
_RL xBeta,visa,Ribcu,QaR |
| 68 |
jmc |
1.2 |
_RL Ct,zetu,L10,charn |
| 69 |
jmc |
1.1 |
|
| 70 |
|
|
C Constants and coefficients (Stull 1988 p640). |
| 71 |
|
|
xBeta = 1.2 _d 0 !Given as 1.25 in Fairall et al.(1996) |
| 72 |
|
|
twoPI = 2. _d 0*PI |
| 73 |
|
|
visa = 1.326 _d -5 |
| 74 |
|
|
C default relative humidity |
| 75 |
|
|
QaR = 0.8 _d 0 |
| 76 |
|
|
|
| 77 |
|
|
C sea surface temperature without skin correction |
| 78 |
|
|
c tsw=theta(i,j,1,bi,bj) |
| 79 |
|
|
tsw = tSurf(i,j) |
| 80 |
|
|
Tas = Tair(i,j,bi,bj) |
| 81 |
|
|
|
| 82 |
|
|
C net upward long wave |
| 83 |
|
|
Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +). |
| 84 |
|
|
|
| 85 |
|
|
C Teten''s return s air svp es in mb |
| 86 |
|
|
es = (1.0007 _d 0 + 3.46 _d -6*p0)*6.1121 _d 0 |
| 87 |
|
|
& *EXP( 17.502 _d 0*tsw/(240.97 _d 0+tsw) ) |
| 88 |
|
|
es = es*0.98 _d 0 !reduced for salinity Kraus 1972 p. 46 |
| 89 |
jmc |
1.2 |
C- convert from mb to spec. humidity kg/kg |
| 90 |
jmc |
1.1 |
qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es) |
| 91 |
jmc |
1.2 |
|
| 92 |
jmc |
1.1 |
tta = Tas+celsius2K |
| 93 |
jmc |
1.2 |
c ttas=tta+gamma_blk*zt |
| 94 |
|
|
c ttt=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk |
| 95 |
|
|
c ttt2=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk-celsius2K |
| 96 |
|
|
c pt = p0*(1.-gamma_blk*CheapHgrid(i,j,bi,bj)/ttas) |
| 97 |
|
|
c & **(gravity/gamma_blk/gasR) |
| 98 |
|
|
c essqt = (1.0007 _d 0 + 3.46 _d -6*pt)*6.1121 _d 0 |
| 99 |
|
|
c & *EXP( 17.502 _d 0*ttt2/(240.97 _d 0+ttt2) ) |
| 100 |
|
|
C- convert from mb to spec. humidity kg/kg |
| 101 |
|
|
c ssqt = 0.62197 _d 0*essqt/(pt -0.378 _d 0*essqt) |
| 102 |
|
|
C- LANL formulation |
| 103 |
jmc |
1.1 |
C saturation no more at the top: |
| 104 |
|
|
ssqt=ssq0*EXP( lath*(ssq1-ssq2/tta) ) / p0 |
| 105 |
|
|
|
| 106 |
|
|
IF (useFreshWaterFlux) THEN |
| 107 |
|
|
q=qair(i,j,bi,bj) |
| 108 |
|
|
ELSE |
| 109 |
|
|
q=QaR*ssqt |
| 110 |
|
|
ENDIF |
| 111 |
|
|
|
| 112 |
|
|
C Wave parameters |
| 113 |
|
|
cwave=gravity*wavesp(i,j,bi,bj)/twoPI |
| 114 |
|
|
lwave=cwave*wavesp(i,j,bi,bj) |
| 115 |
|
|
|
| 116 |
|
|
C Initial guesses |
| 117 |
|
|
zo = 0.0001 _d 0 |
| 118 |
|
|
Wg = 0.5 _d 0 !Gustiness factor initial guess |
| 119 |
|
|
|
| 120 |
|
|
C Air-sea differences - includes warm layer in Dt and Dq |
| 121 |
|
|
c u = (uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2 |
| 122 |
|
|
c & + (vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2 |
| 123 |
jmc |
1.2 |
u = windSq(i,j) |
| 124 |
jmc |
1.1 |
Du= SQRT(u + Wg**2 ) !include gustiness in wind spd. difference |
| 125 |
|
|
u = SQRT(u) |
| 126 |
|
|
Dt=tsw-Tas-gamma_blk*zt !potential temperature difference. |
| 127 |
|
|
Dq=qs-q |
| 128 |
|
|
|
| 129 |
|
|
C **************** neutral coefficients ****************** |
| 130 |
|
|
|
| 131 |
|
|
u10 = Du*LOG(10. _d 0/zo)/LOG(zu/zo) |
| 132 |
|
|
usr = 0.035 _d 0*u10 |
| 133 |
|
|
zo10= 0.011 _d 0*usr*usr/gravity+0.11 _d 0*visa/usr |
| 134 |
|
|
Cd10= (xkar/LOG(10. _d 0/zo10))**2 |
| 135 |
|
|
Ch10= 0.00115 _d 0 |
| 136 |
|
|
Ct10= Ch10/SQRT(Cd10) |
| 137 |
|
|
zot10=10. _d 0/EXP(xkar/Ct10) |
| 138 |
|
|
Cd = (xkar/LOG(zu/zo10))**2 |
| 139 |
|
|
|
| 140 |
|
|
C standard coare3 boundary layer height |
| 141 |
|
|
zi=600. _d 0 |
| 142 |
|
|
|
| 143 |
|
|
C ************* Grachev and Fairall (JAM, 1997) ********** |
| 144 |
|
|
|
| 145 |
|
|
Ct=xkar/LOG(zt/zot10) ! Temperature transfer coefficient |
| 146 |
|
|
CC=xkar*Ct/Cd ! z/L vs Rib linear coefficient |
| 147 |
|
|
Ribcu=-zu/(zi*0.004 _d 0*xBeta**3) ! Saturation or plateau Rib |
| 148 |
jmc |
1.2 |
Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2) |
| 149 |
jmc |
1.1 |
IF (Ribu.LT.0. _d 0) THEN |
| 150 |
|
|
zetu=CC*Ribu/(1. _d 0+Ribu/Ribcu) ! Unstable G and F |
| 151 |
|
|
ELSE |
| 152 |
|
|
zetu=CC*Ribu*(1. _d 0 +27. _d 0/9. _d 0*Ribu/CC) ! Stable |
| 153 |
|
|
ENDIF |
| 154 |
|
|
L10=zu/zetu ! MO length |
| 155 |
|
|
IF (zetu.GT.50. _d 0) THEN |
| 156 |
|
|
nits=1 |
| 157 |
|
|
ELSE |
| 158 |
|
|
nits=3 ! number of iterations |
| 159 |
|
|
ENDIF |
| 160 |
|
|
|
| 161 |
|
|
C First guess M-O stability dependent |
| 162 |
|
|
C scaling params.(u*,t*,q*) to estimate zo and z/L |
| 163 |
|
|
|
| 164 |
|
|
usr= Du*xkar/(LOG(zu/zo10)-psiu(zu/L10)) |
| 165 |
|
|
tsr=-(Dt)*xkar/(LOG(zt/zot10)-psit(zt/L10)) |
| 166 |
|
|
qsr=-(Dq)*xkar/(LOG(zq/zot10)-psit(zq/L10)) |
| 167 |
|
|
|
| 168 |
|
|
charn=0.011 _d 0 !then modify Charnock for high wind speeds Chris data |
| 169 |
|
|
IF (Du.GT.10. _d 0) charn=0.011 _d 0 |
| 170 |
|
|
& + (0.018 _d 0-0.011 _d 0)*(Du-10.)/(18.-10.) |
| 171 |
|
|
IF (Du.GT.18. _d 0) charn=0.018 _d 0 |
| 172 |
|
|
|
| 173 |
|
|
C **** Iterate across u*(t*,q*),zo(zot,zoq) and z/L including cool skin **** |
| 174 |
|
|
|
| 175 |
|
|
DO iter=1,nits |
| 176 |
|
|
IF (WAVEMODEL.EQ.'Smith') THEN |
| 177 |
|
|
zo=charn*usr*usr/gravity + 0.11 _d 0*visa/usr !after Smith 1988 |
| 178 |
|
|
ELSEIF (WAVEMODEL.EQ.'Oost') THEN |
| 179 |
|
|
zo=(50./twoPI)*lwave*(usr/cwave)**4.5 _d 0 |
| 180 |
|
|
& + 0.11 _d 0*visa/usr !Oost et al. |
| 181 |
|
|
ELSEIF (WAVEMODEL.EQ.'TayYel') THEN |
| 182 |
|
|
zo=1200. _d 0*wavesh(i,j,bi,bj)*(wavesh(i,j,bi,bj)/lwave)**4.5 |
| 183 |
|
|
& + 0.11 _d 0*visa/usr !Taylor and Yelland |
| 184 |
|
|
ENDIF |
| 185 |
|
|
rr=zo*usr/visa |
| 186 |
|
|
|
| 187 |
|
|
C *** zoq and zot fitted to results from several ETL cruises ************ |
| 188 |
|
|
|
| 189 |
|
|
IF ( rr.LE.0. ) THEN |
| 190 |
|
|
WRITE(errorMessageUnit,'(A,I8,I4,A,5I4)') |
| 191 |
|
|
& 'CHEAPAML_COARE3_FLUX: myIter,iter=', myIter, iter, |
| 192 |
|
|
& ' , in: i,j,bi,bj,thid=', i, j, bi, bj, myThid |
| 193 |
|
|
WRITE(errorMessageUnit,'(A,1P4E17.9)') |
| 194 |
|
|
& ' rr,zo,usr,visa=', rr, zo, usr, visa |
| 195 |
|
|
WRITE(errorMessageUnit,'(A,1P4E17.9)') |
| 196 |
|
|
& ' L,zu,zL,zt =', L, zu, zL, zt |
| 197 |
|
|
WRITE(errorMessageUnit,'(A,1P4E16.8)') |
| 198 |
|
|
& ' ln(zu/zo),psu,diff,zL*=', LOG(zu/zo), psu, LOG(zu/zo)-psu, |
| 199 |
jmc |
1.2 |
& ( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr ) |
| 200 |
|
|
& /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) ) |
| 201 |
jmc |
1.1 |
WRITE(errorMessageUnit,'(A,1P4E17.9)') |
| 202 |
jmc |
1.2 |
& ' tsr,tta,q,qsr =', tsr, tta, q, qsr |
| 203 |
jmc |
1.1 |
CALL MDS_FLUSH( errorMessageUnit, myThid ) |
| 204 |
|
|
CALL MDS_FLUSH( standardMessageUnit, myThid ) |
| 205 |
|
|
ENDIF |
| 206 |
|
|
zoq = MIN( 1.15 _d -4, 5.5 _d -5/rr**0.6 _d 0 ) |
| 207 |
|
|
zot = zoq |
| 208 |
|
|
|
| 209 |
jmc |
1.2 |
zL=xkar*gravity*zu*( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr ) |
| 210 |
|
|
& /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) ) |
| 211 |
jmc |
1.1 |
L=zu/zL |
| 212 |
|
|
psu=psiu(zu/L) |
| 213 |
|
|
pst=psit(zt/L) |
| 214 |
|
|
usr=Du*xkar/(LOG(zu/zo)-psiu(zu/L)) |
| 215 |
|
|
tsr=-(Dt)*xkar/(LOG(zt/zot)-psit(zt/L)) |
| 216 |
|
|
qsr=-(Dq)*xkar/(LOG(zq/zoq)-psit(zq/L)) |
| 217 |
jmc |
1.2 |
Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr) |
| 218 |
jmc |
1.1 |
IF (Bf.GT.0. _d 0) THEN |
| 219 |
|
|
Wg=xBeta*(Bf*zi)**.333 _d 0 |
| 220 |
|
|
ELSE |
| 221 |
|
|
Wg=0.2 _d 0 |
| 222 |
|
|
ENDIF |
| 223 |
|
|
Du=SQRT(u**2 + Wg**2) !include gustiness in wind spd. |
| 224 |
|
|
ENDDO |
| 225 |
|
|
|
| 226 |
|
|
C compute surface fluxes and other parameters |
| 227 |
|
|
tau=rhoa*usr*usr !stress N/m2 |
| 228 |
jmc |
1.2 |
hf=-cpair*rhoa*usr*tsr !sensible W/m2 |
| 229 |
|
|
ef=-lath*rhoa*usr*qsr !latent W/m2 |
| 230 |
jmc |
1.1 |
evap=-rhoa*usr*qsr |
| 231 |
|
|
cdq = evap/Dq |
| 232 |
jmc |
1.2 |
cdu = tau/Du |
| 233 |
|
|
|
| 234 |
jmc |
1.1 |
q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L)) |
| 235 |
|
|
|
| 236 |
jmc |
1.2 |
C-- compute derivative of surface fluxes relatice to Tsurf |
| 237 |
|
|
C- dSensdTs = -cpair*rhoa*usr*(tsr/Dt) |
| 238 |
|
|
dSensdTs = cpair*rhoa*usr |
| 239 |
|
|
& *xkar/(LOG(zt/zot10)-psit(zt/L10)) |
| 240 |
|
|
C- dEvapdTs = -rhoa*usr* d/dTs(qsr) |
| 241 |
|
|
C d/dTs(qsr)= (-xkar/(LOG(zq/zoq)-psit(zq/L)) )* d/dTs(qs) |
| 242 |
|
|
C d/dTs(qs) = 0.62197 _d 0*p0/(p0 -0.378 _d 0*es)**2 * d/dTs(es) |
| 243 |
|
|
C d/dTs(es) = (0.98)* es * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2 |
| 244 |
|
|
C- this simplifies (using qs) to: |
| 245 |
jmc |
1.1 |
dEvapdTs = rhoa*usr*( xkar/(LOG(zq/zoq)-psit(zq/L)) ) |
| 246 |
|
|
& * qs*p0/(p0 -0.378 _d 0*es) |
| 247 |
|
|
c & *0.98 _d 0 |
| 248 |
|
|
& * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2 |
| 249 |
|
|
|
| 250 |
|
|
if (iceornot.EQ.0) THEN |
| 251 |
|
|
c dLWdTs = 4. _d 0*ocean_emissivity*stefan*tsr*tsr*tsr |
| 252 |
|
|
dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr |
| 253 |
|
|
ELSEIF (iceornot.EQ.2) THEN |
| 254 |
|
|
c dLWdTs = 4. _d 0*snow_emissivity*stefan*tsr*tsr*tsr |
| 255 |
|
|
dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr |
| 256 |
|
|
ELSEIF (iceornot.EQ.1) THEN |
| 257 |
|
|
c dLWdTs = 4. _d 0*ice_emissivity*stefan*tsr*tsr*tsr |
| 258 |
|
|
dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr |
| 259 |
|
|
ENDIF |
| 260 |
|
|
|
| 261 |
jmc |
1.2 |
C-- for now, ignores derivative of q100 relative to Tsurf: |
| 262 |
|
|
dQAdTs = 0. |
| 263 |
jmc |
1.1 |
|
| 264 |
|
|
RETURN |
| 265 |
|
|
END |
| 266 |
|
|
|
| 267 |
|
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 268 |
|
|
|
| 269 |
|
|
CBOP 0 |
| 270 |
|
|
C !ROUTINE: PSIU |
| 271 |
|
|
|
| 272 |
|
|
C !INTERFACE: |
| 273 |
|
|
_RL FUNCTION psiu(zL) |
| 274 |
|
|
|
| 275 |
|
|
C !DESCRIPTION: |
| 276 |
|
|
C psiu and psit evaluate stability function for wind speed and scalars |
| 277 |
|
|
C matching Kansas and free convection forms with weighting f |
| 278 |
|
|
C convective form follows Fairall et al (1996) with profile constants |
| 279 |
|
|
C from Grachev et al (2000) BLM |
| 280 |
|
|
C stable form from Beljaars and Holtslag (1991) |
| 281 |
|
|
|
| 282 |
|
|
C !USES: |
| 283 |
|
|
IMPLICIT NONE |
| 284 |
|
|
#include "EEPARAMS.h" |
| 285 |
|
|
|
| 286 |
|
|
C !INPUT PARAMETERS: |
| 287 |
|
|
_RL zL |
| 288 |
|
|
C !LOCAL VARIABLES: |
| 289 |
|
|
_RL x,y,psik,psic,f,c |
| 290 |
|
|
CEOP |
| 291 |
|
|
|
| 292 |
|
|
IF (zL.LT.0.0) THEN |
| 293 |
|
|
x = (1. - 15.*zL)**.25 !Kansas unstable |
| 294 |
|
|
psik=2.*LOG((1.+x)/2.)+LOG((1.+x*x)/2.)-2.*ATAN(x)+2.*ATAN(oneRL) |
| 295 |
|
|
y = (1. - 10.15 _d 0*zL)**.3333 _d 0 !Convective |
| 296 |
|
|
psic = 1.5*LOG((1.+y+y*y)/3.) |
| 297 |
|
|
& - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) ) |
| 298 |
|
|
& + 4.*ATAN(oneRL)/SQRT(3. _d 0) |
| 299 |
|
|
f = zL*zL/(1.+zL*zL) |
| 300 |
|
|
psiu = (1.-f)*psik+f*psic |
| 301 |
|
|
ELSE |
| 302 |
|
|
c = MIN( 50. _d 0, 0.35 _d 0*zL ) !Stable |
| 303 |
|
|
c psiu=-((1.+1.*zL)**1.+.6667*(zL-14.28)/EXP(c)+8.525) |
| 304 |
|
|
psiu = -( (1.+zL) + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c) |
| 305 |
|
|
& + 8.525 _d 0 ) |
| 306 |
|
|
ENDIF |
| 307 |
|
|
RETURN |
| 308 |
|
|
END |
| 309 |
|
|
|
| 310 |
|
|
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 311 |
|
|
|
| 312 |
|
|
CBOP 0 |
| 313 |
|
|
C !ROUTINE: PSIT |
| 314 |
|
|
|
| 315 |
|
|
C !INTERFACE: |
| 316 |
|
|
_RL FUNCTION psit(zL) |
| 317 |
|
|
|
| 318 |
|
|
C !DESCRIPTION: |
| 319 |
|
|
|
| 320 |
|
|
C !USES: |
| 321 |
|
|
IMPLICIT NONE |
| 322 |
|
|
#include "EEPARAMS.h" |
| 323 |
|
|
|
| 324 |
|
|
C !INPUT PARAMETERS: |
| 325 |
|
|
_RL zL |
| 326 |
|
|
C !LOCAL VARIABLES: |
| 327 |
|
|
_RL x,y,psik,psic,f,c |
| 328 |
|
|
CEOP |
| 329 |
|
|
|
| 330 |
|
|
IF (zL.LT.0.0) THEN |
| 331 |
|
|
x = (1. - 15.*zL)**.5 !Kansas unstable |
| 332 |
|
|
psik = 2.*LOG((1.+x)/2.) |
| 333 |
|
|
y = (1. - 34.15 _d 0*zL)**.3333 _d 0 !Convective |
| 334 |
|
|
psic = 1.5*LOG((1.+y+y*y)/3.) |
| 335 |
|
|
& - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) ) |
| 336 |
|
|
& + 4.*ATAN(oneRL)/SQRT(3. _d 0) |
| 337 |
|
|
f = zL*zL/(1.+zL*zL) |
| 338 |
|
|
psit = (1.-f)*psik+f*psic |
| 339 |
|
|
ELSE |
| 340 |
|
|
c = MIN( 50. _d 0, 0.35 _d 0*zL ) !Stable |
| 341 |
|
|
psit = -( (1.+2.*zL/3.)**1.5 |
| 342 |
|
|
& + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c) |
| 343 |
|
|
& + 8.525 _d 0 ) |
| 344 |
|
|
ENDIF |
| 345 |
|
|
|
| 346 |
|
|
RETURN |
| 347 |
|
|
END |
Parent Directory
| 
Revision Log
| 
Revision Graph