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