1 |
jmc |
1.16 |
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 |
jmc |
1.10 |
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 |
jmc |
1.1 |
C !INTERFACE: |
17 |
jmc |
1.10 |
SUBROUTINE CHEAPAML_COARE3_FLUX( |
18 |
jmc |
1.16 |
I i,j,bi,bj, iceOrNot, |
19 |
|
|
I tSurf, windSq, |
20 |
jmc |
1.15 |
O hf, ef, evap, Rnl, ssqt, q100, cdq, cdu, |
21 |
jmc |
1.16 |
O dSensdTs, dEvapdTs, dLWdTs, dQAdTs, |
22 |
jmc |
1.13 |
I myIter, myThid ) |
23 |
jmc |
1.10 |
|
24 |
|
|
C !DESCRIPTION: |
25 |
|
|
|
26 |
|
|
C !USES: |
27 |
|
|
IMPLICIT NONE |
28 |
jmc |
1.1 |
C === Global variables === |
29 |
|
|
#include "SIZE.h" |
30 |
|
|
#include "EEPARAMS.h" |
31 |
|
|
#include "PARAMS.h" |
32 |
|
|
#include "CHEAPAML.h" |
33 |
|
|
|
34 |
jmc |
1.10 |
C !INPUT PARAMETERS: |
35 |
jmc |
1.16 |
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 |
|
|
C windSq :: relative wind (vs surface motion) speed square |
40 |
|
|
C myIter :: Current iteration number in simulation |
41 |
|
|
C myThid :: My Thread Id number |
42 |
jmc |
1.10 |
INTEGER i,j,bi,bj |
43 |
jmc |
1.16 |
INTEGER iceOrNot |
44 |
|
|
_RL tSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
45 |
|
|
_RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
46 |
jmc |
1.13 |
INTEGER myIter, myThid |
47 |
jmc |
1.10 |
C !OUTPUT PARAMETERS: |
48 |
jmc |
1.15 |
C cdu :: surface drag coeff (for wind stress) |
49 |
jmc |
1.16 |
_RL hf, ef, evap, Rnl, ssqt, q100, cdq |
50 |
|
|
_RL cdu |
51 |
|
|
C derivative vs surf. temp of Sensible, Evap, LW, q100 |
52 |
|
|
_RL dSensdTs, dEvapdTs, dLWdTs, dQAdTs |
53 |
jmc |
1.10 |
CEOP |
54 |
|
|
|
55 |
|
|
C !LOCAL VARIABLES: |
56 |
|
|
INTEGER iter,nits |
57 |
jmc |
1.15 |
_RL tau,L,psu,pst,Bf |
58 |
|
|
_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 |
wienders |
1.6 |
|
63 |
jmc |
1.10 |
C various constants |
64 |
jmc |
1.15 |
_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.15 |
_RL Ct,zetu,L10,charn |
69 |
wienders |
1.4 |
|
70 |
jmc |
1.10 |
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 |
jmc |
1.1 |
|
77 |
jmc |
1.10 |
C sea surface temperature without skin correction |
78 |
jmc |
1.16 |
c tsw=theta(i,j,1,bi,bj) |
79 |
|
|
tsw = tSurf(i,j) |
80 |
|
|
Tas = Tair(i,j,bi,bj) |
81 |
jmc |
1.1 |
|
82 |
jmc |
1.10 |
C net upward long wave |
83 |
|
|
Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +). |
84 |
wienders |
1.2 |
|
85 |
jmc |
1.10 |
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.15 |
C- convert from mb to spec. humidity kg/kg |
90 |
jmc |
1.10 |
qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es) |
91 |
jmc |
1.16 |
|
92 |
jmc |
1.10 |
tta = Tas+celsius2K |
93 |
jmc |
1.15 |
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 |
wienders |
1.4 |
C saturation no more at the top: |
104 |
jmc |
1.10 |
ssqt=ssq0*EXP( lath*(ssq1-ssq2/tta) ) / p0 |
105 |
wienders |
1.4 |
|
106 |
jmc |
1.10 |
IF (useFreshWaterFlux) THEN |
107 |
|
|
q=qair(i,j,bi,bj) |
108 |
|
|
ELSE |
109 |
|
|
q=QaR*ssqt |
110 |
|
|
ENDIF |
111 |
jmc |
1.1 |
|
112 |
jmc |
1.10 |
C Wave parameters |
113 |
|
|
cwave=gravity*wavesp(i,j,bi,bj)/twoPI |
114 |
|
|
lwave=cwave*wavesp(i,j,bi,bj) |
115 |
wienders |
1.2 |
|
116 |
jmc |
1.10 |
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 |
jmc |
1.16 |
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 |
|
|
u = windSq(i,j) |
124 |
jmc |
1.14 |
Du= SQRT(u + Wg**2 ) !include gustiness in wind spd. difference |
125 |
jmc |
1.10 |
u = SQRT(u) |
126 |
wienders |
1.6 |
Dt=tsw-Tas-gamma_blk*zt !potential temperature difference. |
127 |
jmc |
1.7 |
Dq=qs-q |
128 |
jmc |
1.1 |
|
129 |
jmc |
1.10 |
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 |
jmc |
1.1 |
zi=600. _d 0 |
142 |
|
|
|
143 |
jmc |
1.10 |
C ************* Grachev and Fairall (JAM, 1997) ********** |
144 |
|
|
|
145 |
|
|
Ct=xkar/LOG(zt/zot10) ! Temperature transfer coefficient |
146 |
wienders |
1.6 |
CC=xkar*Ct/Cd ! z/L vs Rib linear coefficient |
147 |
jmc |
1.7 |
Ribcu=-zu/(zi*0.004 _d 0*xBeta**3) ! Saturation or plateau Rib |
148 |
jmc |
1.15 |
Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2) |
149 |
jmc |
1.10 |
IF (Ribu.LT.0. _d 0) THEN |
150 |
jmc |
1.1 |
zetu=CC*Ribu/(1. _d 0+Ribu/Ribcu) ! Unstable G and F |
151 |
jmc |
1.10 |
ELSE |
152 |
jmc |
1.1 |
zetu=CC*Ribu*(1. _d 0 +27. _d 0/9. _d 0*Ribu/CC) ! Stable |
153 |
jmc |
1.10 |
ENDIF |
154 |
jmc |
1.1 |
L10=zu/zetu ! MO length |
155 |
jmc |
1.10 |
IF (zetu.GT.50. _d 0) THEN |
156 |
jmc |
1.1 |
nits=1 |
157 |
jmc |
1.10 |
ELSE |
158 |
jmc |
1.1 |
nits=3 ! number of iterations |
159 |
jmc |
1.10 |
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 |
jmc |
1.7 |
charn=0.011 _d 0 !then modify Charnock for high wind speeds Chris data |
169 |
jmc |
1.10 |
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 |
jmc |
1.1 |
zo=charn*usr*usr/gravity + 0.11 _d 0*visa/usr !after Smith 1988 |
178 |
jmc |
1.10 |
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 |
jmc |
1.1 |
zo=1200. _d 0*wavesh(i,j,bi,bj)*(wavesh(i,j,bi,bj)/lwave)**4.5 |
183 |
jmc |
1.10 |
& + 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 |
jmc |
1.13 |
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.15 |
& ( 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.13 |
WRITE(errorMessageUnit,'(A,1P4E17.9)') |
202 |
jmc |
1.15 |
& ' tsr,tta,q,qsr =', tsr, tta, q, qsr |
203 |
jmc |
1.13 |
CALL MDS_FLUSH( errorMessageUnit, myThid ) |
204 |
|
|
CALL MDS_FLUSH( standardMessageUnit, myThid ) |
205 |
|
|
ENDIF |
206 |
jmc |
1.10 |
zoq = MIN( 1.15 _d -4, 5.5 _d -5/rr**0.6 _d 0 ) |
207 |
|
|
zot = zoq |
208 |
|
|
|
209 |
jmc |
1.15 |
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.10 |
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.15 |
Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr) |
218 |
jmc |
1.10 |
IF (Bf.GT.0. _d 0) THEN |
219 |
jmc |
1.1 |
Wg=xBeta*(Bf*zi)**.333 _d 0 |
220 |
jmc |
1.10 |
ELSE |
221 |
jmc |
1.1 |
Wg=0.2 _d 0 |
222 |
jmc |
1.10 |
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 |
jmc |
1.12 |
tau=rhoa*usr*usr !stress N/m2 |
228 |
jmc |
1.15 |
hf=-cpair*rhoa*usr*tsr !sensible W/m2 |
229 |
|
|
ef=-lath*rhoa*usr*qsr !latent W/m2 |
230 |
jmc |
1.10 |
evap=-rhoa*usr*qsr |
231 |
|
|
cdq = evap/Dq |
232 |
jmc |
1.14 |
cdu = tau/Du |
233 |
jmc |
1.16 |
|
234 |
jmc |
1.10 |
q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L)) |
235 |
|
|
|
236 |
jmc |
1.16 |
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 |
|
|
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 |
|
|
C-- for now, ignores derivative of q100 relative to Tsurf: |
262 |
|
|
dQAdTs = 0. |
263 |
|
|
|
264 |
jmc |
1.10 |
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 |
jmc |
1.7 |
|
346 |
jmc |
1.10 |
RETURN |
347 |
|
|
END |