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