C $Header: /u/gcmpack/MITgcm/pkg/cheapaml/cheapaml_coare3_flux.F,v 1.16 2013/06/17 13:45:14 jmc Exp $
C $Name: $
#include "CHEAPAML_OPTIONS.h"
C-- File cheapaml_coare3_flux.F:
C-- Contents:
C-- o CHEAPAML_COARE3_FLUX
C-- o PSIU (Function)
C-- o PSIT (Function)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: CHEAPAML_COARE3_FLUX
C !INTERFACE:
SUBROUTINE CHEAPAML_COARE3_FLUX(
I i,j,bi,bj, iceOrNot,
I tSurf, windSq,
O hf, ef, evap, Rnl, ssqt, q100, cdq, cdu,
O dSensdTs, dEvapdTs, dLWdTs, dQAdTs,
I myIter, myThid )
C !DESCRIPTION:
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "CHEAPAML.h"
C !INPUT PARAMETERS:
C i, j :: local indices of current grid-point
C bi, bj :: current tile indices
C iceOrNot :: 0=open water, 1=ice cover, 2=ice+snow
C tSurf :: surface temperature
C windSq :: relative wind (vs surface motion) speed square
C myIter :: Current iteration number in simulation
C myThid :: My Thread Id number
INTEGER i,j,bi,bj
INTEGER iceOrNot
_RL tSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER myIter, myThid
C !OUTPUT PARAMETERS:
C cdu :: surface drag coeff (for wind stress)
_RL hf, ef, evap, Rnl, ssqt, q100, cdq
_RL cdu
C derivative vs surf. temp of Sensible, Evap, LW, q100
_RL dSensdTs, dEvapdTs, dLWdTs, dQAdTs
CEOP
C !LOCAL VARIABLES:
INTEGER iter,nits
_RL tau,L,psu,pst,Bf
_RL CD,usr,tsr,qsr
c _RL ttas,ttt,ttt2,pt,essqt
_RL zo,zot,zoq,RR,zL
_RL twoPI,cwave,lwave
C various constants
_RL u,q,Tas,tta,zi,es,qs,tsw
_RL psiu,psit,zot10,Ct10,CC,Ribu
_RL Du,Wg,Dt,Dq,u10,zo10,Cd10,Ch10
_RL xBeta,visa,Ribcu,QaR
_RL Ct,zetu,L10,charn
C Constants and coefficients (Stull 1988 p640).
xBeta = 1.2 _d 0 !Given as 1.25 in Fairall et al.(1996)
twoPI = 2. _d 0*PI
visa = 1.326 _d -5
C default relative humidity
QaR = 0.8 _d 0
C sea surface temperature without skin correction
c tsw=theta(i,j,1,bi,bj)
tsw = tSurf(i,j)
Tas = Tair(i,j,bi,bj)
C net upward long wave
Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +).
C Teten''s return s air svp es in mb
es = (1.0007 _d 0 + 3.46 _d -6*p0)*6.1121 _d 0
& *EXP( 17.502 _d 0*tsw/(240.97 _d 0+tsw) )
es = es*0.98 _d 0 !reduced for salinity Kraus 1972 p. 46
C- convert from mb to spec. humidity kg/kg
qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es)
tta = Tas+celsius2K
c ttas=tta+gamma_blk*zt
c ttt=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk
c ttt2=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk-celsius2K
c pt = p0*(1.-gamma_blk*CheapHgrid(i,j,bi,bj)/ttas)
c & **(gravity/gamma_blk/gasR)
c essqt = (1.0007 _d 0 + 3.46 _d -6*pt)*6.1121 _d 0
c & *EXP( 17.502 _d 0*ttt2/(240.97 _d 0+ttt2) )
C- convert from mb to spec. humidity kg/kg
c ssqt = 0.62197 _d 0*essqt/(pt -0.378 _d 0*essqt)
C- LANL formulation
C saturation no more at the top:
ssqt=ssq0*EXP( lath*(ssq1-ssq2/tta) ) / p0
IF (useFreshWaterFlux) THEN
q=qair(i,j,bi,bj)
ELSE
q=QaR*ssqt
ENDIF
C Wave parameters
cwave=gravity*wavesp(i,j,bi,bj)/twoPI
lwave=cwave*wavesp(i,j,bi,bj)
C Initial guesses
zo = 0.0001 _d 0
Wg = 0.5 _d 0 !Gustiness factor initial guess
C Air-sea differences - includes warm layer in Dt and Dq
c u = (uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2
c & + (vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2
u = windSq(i,j)
Du= SQRT(u + Wg**2 ) !include gustiness in wind spd. difference
u = SQRT(u)
Dt=tsw-Tas-gamma_blk*zt !potential temperature difference.
Dq=qs-q
C **************** neutral coefficients ******************
u10 = Du*LOG(10. _d 0/zo)/LOG(zu/zo)
usr = 0.035 _d 0*u10
zo10= 0.011 _d 0*usr*usr/gravity+0.11 _d 0*visa/usr
Cd10= (xkar/LOG(10. _d 0/zo10))**2
Ch10= 0.00115 _d 0
Ct10= Ch10/SQRT(Cd10)
zot10=10. _d 0/EXP(xkar/Ct10)
Cd = (xkar/LOG(zu/zo10))**2
C standard coare3 boundary layer height
zi=600. _d 0
C ************* Grachev and Fairall (JAM, 1997) **********
Ct=xkar/LOG(zt/zot10) ! Temperature transfer coefficient
CC=xkar*Ct/Cd ! z/L vs Rib linear coefficient
Ribcu=-zu/(zi*0.004 _d 0*xBeta**3) ! Saturation or plateau Rib
Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2)
IF (Ribu.LT.0. _d 0) THEN
zetu=CC*Ribu/(1. _d 0+Ribu/Ribcu) ! Unstable G and F
ELSE
zetu=CC*Ribu*(1. _d 0 +27. _d 0/9. _d 0*Ribu/CC) ! Stable
ENDIF
L10=zu/zetu ! MO length
IF (zetu.GT.50. _d 0) THEN
nits=1
ELSE
nits=3 ! number of iterations
ENDIF
C First guess M-O stability dependent
C scaling params.(u*,t*,q*) to estimate zo and z/L
usr= Du*xkar/(LOG(zu/zo10)-psiu(zu/L10))
tsr=-(Dt)*xkar/(LOG(zt/zot10)-psit(zt/L10))
qsr=-(Dq)*xkar/(LOG(zq/zot10)-psit(zq/L10))
charn=0.011 _d 0 !then modify Charnock for high wind speeds Chris data
IF (Du.GT.10. _d 0) charn=0.011 _d 0
& + (0.018 _d 0-0.011 _d 0)*(Du-10.)/(18.-10.)
IF (Du.GT.18. _d 0) charn=0.018 _d 0
C **** Iterate across u*(t*,q*),zo(zot,zoq) and z/L including cool skin ****
DO iter=1,nits
IF (WAVEMODEL.EQ.'Smith') THEN
zo=charn*usr*usr/gravity + 0.11 _d 0*visa/usr !after Smith 1988
ELSEIF (WAVEMODEL.EQ.'Oost') THEN
zo=(50./twoPI)*lwave*(usr/cwave)**4.5 _d 0
& + 0.11 _d 0*visa/usr !Oost et al.
ELSEIF (WAVEMODEL.EQ.'TayYel') THEN
zo=1200. _d 0*wavesh(i,j,bi,bj)*(wavesh(i,j,bi,bj)/lwave)**4.5
& + 0.11 _d 0*visa/usr !Taylor and Yelland
ENDIF
rr=zo*usr/visa
C *** zoq and zot fitted to results from several ETL cruises ************
IF ( rr.LE.0. ) THEN
WRITE(errorMessageUnit,'(A,I8,I4,A,5I4)')
& 'CHEAPAML_COARE3_FLUX: myIter,iter=', myIter, iter,
& ' , in: i,j,bi,bj,thid=', i, j, bi, bj, myThid
WRITE(errorMessageUnit,'(A,1P4E17.9)')
& ' rr,zo,usr,visa=', rr, zo, usr, visa
WRITE(errorMessageUnit,'(A,1P4E17.9)')
& ' L,zu,zL,zt =', L, zu, zL, zt
WRITE(errorMessageUnit,'(A,1P4E16.8)')
& ' ln(zu/zo),psu,diff,zL*=', LOG(zu/zo), psu, LOG(zu/zo)-psu,
& ( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )
& /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) )
WRITE(errorMessageUnit,'(A,1P4E17.9)')
& ' tsr,tta,q,qsr =', tsr, tta, q, qsr
CALL MDS_FLUSH( errorMessageUnit, myThid )
CALL MDS_FLUSH( standardMessageUnit, myThid )
ENDIF
zoq = MIN( 1.15 _d -4, 5.5 _d -5/rr**0.6 _d 0 )
zot = zoq
zL=xkar*gravity*zu*( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )
& /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) )
L=zu/zL
psu=psiu(zu/L)
pst=psit(zt/L)
usr=Du*xkar/(LOG(zu/zo)-psiu(zu/L))
tsr=-(Dt)*xkar/(LOG(zt/zot)-psit(zt/L))
qsr=-(Dq)*xkar/(LOG(zq/zoq)-psit(zq/L))
Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr)
IF (Bf.GT.0. _d 0) THEN
Wg=xBeta*(Bf*zi)**.333 _d 0
ELSE
Wg=0.2 _d 0
ENDIF
Du=SQRT(u**2 + Wg**2) !include gustiness in wind spd.
ENDDO
C compute surface fluxes and other parameters
tau=rhoa*usr*usr !stress N/m2
hf=-cpair*rhoa*usr*tsr !sensible W/m2
ef=-lath*rhoa*usr*qsr !latent W/m2
evap=-rhoa*usr*qsr
cdq = evap/Dq
cdu = tau/Du
q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L))
C-- compute derivative of surface fluxes relatice to Tsurf
C- dSensdTs = -cpair*rhoa*usr*(tsr/Dt)
dSensdTs = cpair*rhoa*usr
& *xkar/(LOG(zt/zot10)-psit(zt/L10))
C- dEvapdTs = -rhoa*usr* d/dTs(qsr)
C d/dTs(qsr)= (-xkar/(LOG(zq/zoq)-psit(zq/L)) )* d/dTs(qs)
C d/dTs(qs) = 0.62197 _d 0*p0/(p0 -0.378 _d 0*es)**2 * d/dTs(es)
C d/dTs(es) = (0.98)* es * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2
C- this simplifies (using qs) to:
dEvapdTs = rhoa*usr*( xkar/(LOG(zq/zoq)-psit(zq/L)) )
& * qs*p0/(p0 -0.378 _d 0*es)
c & *0.98 _d 0
& * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2
if (iceornot.EQ.0) THEN
c dLWdTs = 4. _d 0*ocean_emissivity*stefan*tsr*tsr*tsr
dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
ELSEIF (iceornot.EQ.2) THEN
c dLWdTs = 4. _d 0*snow_emissivity*stefan*tsr*tsr*tsr
dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
ELSEIF (iceornot.EQ.1) THEN
c dLWdTs = 4. _d 0*ice_emissivity*stefan*tsr*tsr*tsr
dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
ENDIF
C-- for now, ignores derivative of q100 relative to Tsurf:
dQAdTs = 0.
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C !ROUTINE: PSIU
C !INTERFACE:
_RL FUNCTION psiu(zL)
C !DESCRIPTION:
C psiu and psit evaluate stability function for wind speed and scalars
C matching Kansas and free convection forms with weighting f
C convective form follows Fairall et al (1996) with profile constants
C from Grachev et al (2000) BLM
C stable form from Beljaars and Holtslag (1991)
C !USES:
IMPLICIT NONE
#include "EEPARAMS.h"
C !INPUT PARAMETERS:
_RL zL
C !LOCAL VARIABLES:
_RL x,y,psik,psic,f,c
CEOP
IF (zL.LT.0.0) THEN
x = (1. - 15.*zL)**.25 !Kansas unstable
psik=2.*LOG((1.+x)/2.)+LOG((1.+x*x)/2.)-2.*ATAN(x)+2.*ATAN(oneRL)
y = (1. - 10.15 _d 0*zL)**.3333 _d 0 !Convective
psic = 1.5*LOG((1.+y+y*y)/3.)
& - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) )
& + 4.*ATAN(oneRL)/SQRT(3. _d 0)
f = zL*zL/(1.+zL*zL)
psiu = (1.-f)*psik+f*psic
ELSE
c = MIN( 50. _d 0, 0.35 _d 0*zL ) !Stable
c psiu=-((1.+1.*zL)**1.+.6667*(zL-14.28)/EXP(c)+8.525)
psiu = -( (1.+zL) + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c)
& + 8.525 _d 0 )
ENDIF
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C !ROUTINE: PSIT
C !INTERFACE:
_RL FUNCTION psit(zL)
C !DESCRIPTION:
C !USES:
IMPLICIT NONE
#include "EEPARAMS.h"
C !INPUT PARAMETERS:
_RL zL
C !LOCAL VARIABLES:
_RL x,y,psik,psic,f,c
CEOP
IF (zL.LT.0.0) THEN
x = (1. - 15.*zL)**.5 !Kansas unstable
psik = 2.*LOG((1.+x)/2.)
y = (1. - 34.15 _d 0*zL)**.3333 _d 0 !Convective
psic = 1.5*LOG((1.+y+y*y)/3.)
& - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) )
& + 4.*ATAN(oneRL)/SQRT(3. _d 0)
f = zL*zL/(1.+zL*zL)
psit = (1.-f)*psik+f*psic
ELSE
c = MIN( 50. _d 0, 0.35 _d 0*zL ) !Stable
psit = -( (1.+2.*zL/3.)**1.5
& + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c)
& + 8.525 _d 0 )
ENDIF
RETURN
END