1 |
jmc |
1.3 |
C $Header: /u/gcmpack/MITgcm/pkg/cheapaml/cheapaml_coare3_flux.F,v 1.2 2011/02/24 16:11:41 wienders Exp $ |
2 |
jmc |
1.1 |
C $Name: $ |
3 |
|
|
|
4 |
|
|
#include "CHEAPAML_OPTIONS.h" |
5 |
|
|
|
6 |
|
|
C !ROUTINE: CHEAPAML_COARE3_flux |
7 |
|
|
C !INTERFACE: |
8 |
|
|
subroutine cheapaml_COARE3_flux |
9 |
|
|
&(i,j,bi,bj,hf,ef,evap,Rnl,ssqt,q100) |
10 |
|
|
c |
11 |
|
|
implicit none |
12 |
|
|
C === Global variables === |
13 |
|
|
#include "SIZE.h" |
14 |
|
|
#include "EEPARAMS.h" |
15 |
|
|
#include "PARAMS.h" |
16 |
|
|
#include "CHEAPAML.h" |
17 |
|
|
#include "DYNVARS.h" |
18 |
|
|
|
19 |
|
|
integer iter,i,j,bi,bj,nits |
20 |
|
|
_RL hf,ef,evap,tau,L,psu,pst,Bf |
21 |
|
|
_RL CD,usr,tsr,qsr,q100,ssqt,ttas,essqt |
22 |
|
|
_RL zo,zot,zoq,RR,zL,pt,ttt,tta |
23 |
|
|
_RL Rnl,es,twopi,cwave,lwave |
24 |
|
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
25 |
|
|
c various constants |
26 |
|
|
c |
27 |
|
|
_RL u,ts,q,zi,qs,tsw |
28 |
|
|
_RL psiu,psit,zot10,Ct10,CC,Ribu |
29 |
|
|
_RL Du,Wg,Dt,Dq,u10,zo10,Cd10,Ch10 |
30 |
|
|
_RL xBeta,visa,Ribcu,QaR |
31 |
|
|
_RL Ct,zetu,L10,Tas,ta,charn |
32 |
|
|
c |
33 |
jmc |
1.3 |
c |
34 |
|
|
c Constants and coefficients (Stull 1988 p640). |
35 |
jmc |
1.1 |
xBeta=1.2 !Given as 1.25 in Fairall et al.(1996) |
36 |
|
|
twopi=2. _d 0*pi |
37 |
|
|
visa=1.326 _d -5 |
38 |
|
|
c default relative humidity |
39 |
|
|
QaR=0.8 _d 0 |
40 |
|
|
|
41 |
|
|
c sea surface temperature without skin correction |
42 |
|
|
ts=theta(i,j,1,bi,bj) |
43 |
|
|
tsw=ts |
44 |
|
|
Tas=Tair(i,j,bi,bj) |
45 |
|
|
|
46 |
wienders |
1.2 |
|
47 |
jmc |
1.1 |
c net upward long wave |
48 |
jmc |
1.3 |
Rnl= 0.97*(stefan*(tsw+Celsius2K)**4) !Net longwave (up = +). |
49 |
jmc |
1.1 |
c |
50 |
wienders |
1.2 |
c Teten''s returns air svp es in mb |
51 |
jmc |
1.1 |
es = (1.0007+3.46e-6*p0)*6.1121*dexp(17.502*tsw/(240.97+tsw)) !mb |
52 |
|
|
es=es*0.98 !reduced for salinity Kraus 1972 p. 46 |
53 |
|
|
qs=.62197*es/(p0-0.378*es) !convert from mb to spec. humidity kg/kg |
54 |
wienders |
1.2 |
tta=Tas+Celsius2K |
55 |
jmc |
1.1 |
ttas=tta+gamma_blk*zt |
56 |
wienders |
1.2 |
ttt=tta-(cheapaml_h - zt)*gamma_blk |
57 |
jmc |
1.1 |
pt=p0*(1-gamma_blk*cheapaml_h/ttas)**(gravity/gamma_blk/gasR) |
58 |
|
|
essqt = (1.0007+3.46e-6*pt)*6.1121*dexp(17.502*tas/(240.97+tas)) !mb |
59 |
|
|
ssqt = .62197*essqt/(pt-0.378*essqt) !convert from mb to spec. humidity kg/kg |
60 |
|
|
|
61 |
|
|
if (useFreshWaterFlux)then |
62 |
wienders |
1.2 |
q=qair(i,j,bi,bj) |
63 |
jmc |
1.1 |
else |
64 |
wienders |
1.2 |
q=QaR*ssqt |
65 |
jmc |
1.1 |
endif |
66 |
|
|
|
67 |
wienders |
1.2 |
|
68 |
jmc |
1.1 |
c Wave parameters |
69 |
|
|
cwave=gravity*wavesp(i,j,bi,bj)/twopi |
70 |
|
|
lwave=cwave*wavesp(i,j,bi,bj) |
71 |
|
|
c |
72 |
|
|
c Initial guesses |
73 |
|
|
zo=0.0001 |
74 |
|
|
Wg=0.5 !Gustiness factor initial guess |
75 |
|
|
c |
76 |
|
|
c Air-sea differences - includes warm layer in Dt and Dq |
77 |
|
|
u=(uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2+ |
78 |
|
|
&(vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2 |
79 |
|
|
u=dsqrt(u) |
80 |
|
|
Du=(u**2.+Wg**2.)**.5 !include gustiness in wind spd. difference |
81 |
|
|
Dt=tsw-Tas-gamma_blk*zt !potential temperature difference. |
82 |
jmc |
1.3 |
Dq=qs-q |
83 |
jmc |
1.1 |
c |
84 |
|
|
c **************** neutral coefficients ****************** |
85 |
|
|
c |
86 |
|
|
u10=Du*dlog(10. _d 0/zo)/dlog(zu/zo) |
87 |
|
|
usr=0.035*u10 |
88 |
|
|
zo10=0.011*usr*usr/gravity+0.11*visa/usr |
89 |
|
|
Cd10=(xkar/dlog(10. _d 0/zo10))**2 |
90 |
|
|
Ch10=0.00115 _d 0 |
91 |
|
|
Ct10=Ch10/sqrt(Cd10) |
92 |
|
|
zot10=10._d 0/dexp(xkar/Ct10) |
93 |
|
|
Cd=(xkar/dlog(zu/zo10))**2 |
94 |
|
|
|
95 |
|
|
c standard coare3 boundary layer height |
96 |
|
|
zi=600. _d 0 |
97 |
|
|
|
98 |
jmc |
1.3 |
c |
99 |
jmc |
1.1 |
c ************* Grachev and Fairall (JAM, 1997) ********** |
100 |
|
|
c |
101 |
wienders |
1.2 |
ta=Tas+Celsius2K |
102 |
jmc |
1.1 |
Ct=xkar/dlog(zt/zot10) ! Temperature transfer coefficient |
103 |
|
|
CC=xkar*Ct/Cd ! z/L vs Rib linear coefficient |
104 |
jmc |
1.3 |
Ribcu=-zu/(zi*0.004 _d 0*xBeta**3) ! Saturation or plateau Rib |
105 |
jmc |
1.1 |
Ribu=-gravity*zu*(Dt+0.61 _d 0*ta*Dq)/(ta*Du**2) |
106 |
|
|
if (Ribu.lt.0. _d 0) then |
107 |
|
|
zetu=CC*Ribu/(1. _d 0+Ribu/Ribcu) ! Unstable G and F |
108 |
|
|
else |
109 |
|
|
zetu=CC*Ribu*(1. _d 0 +27. _d 0/9. _d 0*Ribu/CC) ! Stable |
110 |
|
|
endif |
111 |
|
|
L10=zu/zetu ! MO length |
112 |
|
|
if (zetu.gt.50. _d 0) then |
113 |
|
|
nits=1 |
114 |
|
|
else |
115 |
|
|
nits=3 ! number of iterations |
116 |
|
|
endif |
117 |
|
|
c |
118 |
|
|
c First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate zo and z/L |
119 |
|
|
c |
120 |
|
|
usr= Du*xkar/(dlog(zu/zo10)-psiu(zu/L10)) |
121 |
|
|
tsr=-(Dt)*xkar/(dlog(zt/zot10)-psit(zt/L10)) |
122 |
|
|
qsr=-(Dq)*xkar/(dlog(zq/zot10)-psit(zq/L10)) |
123 |
jmc |
1.3 |
c |
124 |
|
|
charn=0.011 _d 0 !then modify Charnock for high wind speeds Chris s data |
125 |
|
|
if(Du.gt.10. _d 0) charn=0.011 _d 0 |
126 |
|
|
& +(0.018-0.011)*(Du-10.)/(18.0-10.0) |
127 |
jmc |
1.1 |
if(Du.gt.18. _d 0) charn=0.018 _d 0 |
128 |
jmc |
1.3 |
c |
129 |
jmc |
1.1 |
c **** Iterate across u*(t*,q*),zo(zot,zoq) and z/L including cool skin **** |
130 |
|
|
c |
131 |
|
|
do iter=1,nits |
132 |
|
|
if(WAVEMODEL.eq.'Smith') then |
133 |
|
|
zo=charn*usr*usr/gravity + 0.11 _d 0*visa/usr !after Smith 1988 |
134 |
|
|
else if(WAVEMODEL.eq.'Oost') then |
135 |
|
|
zo=(50./twopi)*lwave*(usr/cwave)**4.5 _d 0+0.11*visa/usr !Oost et al. |
136 |
|
|
else if(WAVEMODEL.eq.'TayYel') then |
137 |
|
|
zo=1200. _d 0*wavesh(i,j,bi,bj)*(wavesh(i,j,bi,bj)/lwave)**4.5 |
138 |
jmc |
1.3 |
& +0.11 _d 0*visa/usr !Taylor and Yelland |
139 |
|
|
endif |
140 |
jmc |
1.1 |
rr=zo*usr/visa |
141 |
|
|
c |
142 |
|
|
c *** zoq and zot fitted to results from several ETL cruises ************ |
143 |
|
|
c |
144 |
|
|
zoq=min(1.15 _d -4,5.5 _d -5/rr**0.6 _d 0) |
145 |
|
|
zot=zoq |
146 |
|
|
c |
147 |
|
|
zL=xkar*gravity*zu*(tsr*(1.+0.61*q)+0.61*ta*qsr) |
148 |
|
|
& /(ta*usr*usr*(1. _d 0+0.61 _d 0*q)) |
149 |
|
|
L=zu/zL |
150 |
|
|
psu=psiu(zu/L) |
151 |
|
|
pst=psit(zt/L) |
152 |
|
|
usr=Du*xkar/(dlog(zu/zo)-psiu(zu/L)) |
153 |
|
|
tsr=-(Dt)*xkar/(dlog(zt/zot)-psit(zt/L)) |
154 |
|
|
qsr=-(Dq)*xkar/(dlog(zq/zoq)-psit(zq/L)) |
155 |
|
|
Bf=-gravity/ta*usr*(tsr+0.61 _d 0*ta*qsr) |
156 |
|
|
if (Bf.gt.0. _d 0) then |
157 |
|
|
Wg=xBeta*(Bf*zi)**.333 _d 0 |
158 |
|
|
else |
159 |
|
|
Wg=0.2 _d 0 |
160 |
|
|
endif |
161 |
|
|
Du=sqrt(u**2.+Wg**2.) !include gustiness in wind spd. |
162 |
|
|
enddo |
163 |
jmc |
1.3 |
|
164 |
jmc |
1.1 |
c compute surface fluxes and other parameters |
165 |
|
|
tau=rhoa*usr*usr*u/Du !stress N/m2 |
166 |
|
|
hf=-cpair*rhoa*usr*tsr !sensible W/m2 |
167 |
|
|
ef=-lath*rhoa*usr*qsr !latent W/m2 |
168 |
|
|
evap=-rhoa*usr*qsr |
169 |
|
|
if(.NOT.useStressOption)THEN |
170 |
|
|
ustress(i,j,bi,bj)=tau*(uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))/u |
171 |
|
|
vstress(i,j,bi,bj)=tau*(vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))/u |
172 |
|
|
endif |
173 |
|
|
q100=qs+qsr*(dlog(100. _d 0/zoq)-psit(100. _d 0/L)) |
174 |
|
|
c |
175 |
jmc |
1.3 |
return |
176 |
jmc |
1.1 |
end |
177 |
|
|
c |
178 |
|
|
c------------------------------------------------------------------ |
179 |
|
|
function psiu(zL) |
180 |
|
|
c |
181 |
|
|
c psiu and psit evaluate stability function for wind speed and scalars |
182 |
|
|
c matching Kansas and free convection forms with weighting f |
183 |
|
|
c convective form follows Fairall et al (1996) with profile constants |
184 |
|
|
c from Grachev et al (2000) BLM |
185 |
|
|
c stable form from Beljaars and Holtslag (1991) |
186 |
|
|
c |
187 |
|
|
implicit none |
188 |
|
|
_RL zL,x,y,psik,psic,f,psiu,c |
189 |
|
|
if(zL.lt.0.0) then |
190 |
|
|
x=(1.-15.*zL)**.25 !Kansas unstable |
191 |
|
|
psik=2.*dlog((1.+x)/2.)+dlog((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) |
192 |
|
|
y=(1.-10.15*zL)**.3333 !Convective |
193 |
|
|
psic=1.5*dlog((1.+y+y*y)/3.)-sqrt(3.)*atan((1.+2.*y)/sqrt(3.)) |
194 |
|
|
& +4.*atan(1.)/sqrt(3.) |
195 |
|
|
f=zL*zL/(1.+zL*zL) |
196 |
|
|
psiu=(1.-f)*psik+f*psic |
197 |
|
|
else |
198 |
|
|
c=min(50.,0.35*zL) !Stable |
199 |
|
|
psiu=-((1.+1.*zL)**1.+.6667*(zL-14.28)/dexp(c)+8.525) |
200 |
|
|
endif |
201 |
|
|
return |
202 |
|
|
end |
203 |
jmc |
1.3 |
c-------------------------------------------------------------- |
204 |
jmc |
1.1 |
function psit(zL) |
205 |
jmc |
1.3 |
|
206 |
jmc |
1.1 |
implicit none |
207 |
|
|
_RL zL,x,y,psik,psic,f,psit,c |
208 |
|
|
if(zL.lt.0.0) then |
209 |
|
|
x=(1.-15.*zL)**.5 !Kansas unstable |
210 |
|
|
psik=2.*dlog((1.+x)/2.) |
211 |
|
|
y=(1.-34.15*zL)**.3333 !Convective |
212 |
|
|
psic=1.5*dlog((1.+y+y*y)/3.)-sqrt(3.)*atan((1.+2.*y)/sqrt(3.)) |
213 |
|
|
& +4.*atan(1.)/sqrt(3.) |
214 |
|
|
f=zL*zL/(1.+zL*zL) |
215 |
|
|
psit=(1.-f)*psik+f*psic |
216 |
|
|
else |
217 |
|
|
c=min(50.,0.35*zL) !Stable |
218 |
|
|
psit=-((1.+2.*zL/3.)**1.5+.6667*(zL-14.28)/dexp(c)+8.525) |
219 |
|
|
endif |
220 |
|
|
return |
221 |
|
|
end |
222 |
jmc |
1.3 |
|
223 |
jmc |
1.1 |
c------------------------------------------------------------- |