| 1 |
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 |
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 |
I tSurf, windSq, |
| 20 |
O hf, ef, evap, Rnl, ssqt, q100, cdq, cdu, |
| 21 |
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 |
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 |
INTEGER i,j,bi,bj |
| 43 |
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 |
INTEGER myIter, myThid |
| 47 |
C !OUTPUT PARAMETERS: |
| 48 |
C cdu :: surface drag coeff (for wind stress) |
| 49 |
_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 |
CEOP |
| 54 |
|
| 55 |
C !LOCAL VARIABLES: |
| 56 |
INTEGER iter,nits |
| 57 |
_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 |
|
| 63 |
C various constants |
| 64 |
_RL u,q,Tas,tta,zi,es,qs,tsw |
| 65 |
_RL psiu,psit,zot10,Ct10,CC,Ribu |
| 66 |
_RL Du,Wg,Dt,Dq,u10,zo10,Cd10,Ch10 |
| 67 |
_RL xBeta,visa,Ribcu,QaR |
| 68 |
_RL Ct,zetu,L10,charn |
| 69 |
|
| 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 |
C- convert from mb to spec. humidity kg/kg |
| 90 |
qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es) |
| 91 |
|
| 92 |
tta = Tas+celsius2K |
| 93 |
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 |
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 |
u = windSq(i,j) |
| 124 |
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 |
Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2) |
| 149 |
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 |
& ( 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 |
WRITE(errorMessageUnit,'(A,1P4E17.9)') |
| 202 |
& ' tsr,tta,q,qsr =', tsr, tta, q, qsr |
| 203 |
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 |
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 |
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 |
Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr) |
| 218 |
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 |
hf=-cpair*rhoa*usr*tsr !sensible W/m2 |
| 229 |
ef=-lath*rhoa*usr*qsr !latent W/m2 |
| 230 |
evap=-rhoa*usr*qsr |
| 231 |
cdq = evap/Dq |
| 232 |
cdu = tau/Du |
| 233 |
|
| 234 |
q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L)) |
| 235 |
|
| 236 |
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 |
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