C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lay.F,v 1.2 2006/05/30 22:44:54 mlosch Exp $
C $Name: $
#include "BULK_FORCE_OPTIONS.h"
CBOP
C !ROUTINE: BULKF_FORMULA_LAY
C !INTERFACE:
SUBROUTINE BULKF_FORMULA_LAY(
I uw, vw, ws, Ta, Qa, tsfCel,
O flwupa, flha, fsha, df0dT,
O ust, vst, evp, ssq, dEvdT,
I iceornot, i,j,bi,bj,myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE BULKF_FORMULA_LAY
C | o Calculate bulk formula fluxes over open ocean or seaice
C | Large and Yeager, 2004, NCAR/TN-460+STR.
C *==========================================================*
C \ev
C
C === Turbulent Fluxes :
C * use the approach "B": shift coeff to height & stability of the
C atmosphere state (instead of "C": shift temp & humid to the height
C of wind, then shift the coeff to this height & stability of the atmos).
C * similar to EXF (except over sea-ice) ; default parameter values
C taken from Large & Yeager.
C * assume that Qair & Tair inputs are from the same height (zq=zt)
C * formulae in short:
C wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
C Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
C Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
C = -Evap * Lvap
C with Ws = wind speed = sqrt(del.u^2 +del.v^2) ;
C del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ;
C Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number
C respectively [no-units], function of height & stability
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "BULKF_PARAMS.h"
C !INPUT/OUTPUT PARAMETERS:
C input:
_RL uw ! zonal wind speed (at grid center) [m/s]
_RL vw ! meridional wind speed (at grid center) [m/s]
_RL ws ! wind speed [m/s] at height zwd
_RL Ta ! air temperature [K] at height zth
_RL Qa ! specific humidity [kg/kg] at heigth zth
_RL tsfCel ! sea-ice or sea surface temperature [oC]
INTEGER iceornot ! 0=open water, 1=sea-ice, 2=sea-ice with snow
INTEGER i,j, bi,bj !current grid point indices
INTEGER myThid ! my Thread Id number
C output:
_RL flwupa ! upward long wave radiation (>0 upward) [W/m2]
_RL flha ! latent heat flux (>0 downward) [W/m2]
_RL fsha ! sensible heat flux (>0 downward) [W/m2]
_RL df0dT ! derivative of heat flux with respect to Tsf [W/m2/K]
_RL ust ! zonal wind stress (at grid center) [N/m2]
_RL vst ! meridional wind stress (at grid center)[N/m2]
_RL evp ! evaporation rate (over open water) [kg/m2/s]
_RL ssq ! surface specific humidity [kg/kg]
_RL dEvdT ! derivative of evap. with respect to Tsf [kg/m2/s/K]
CEOP
#ifdef ALLOW_BULK_FORCE
C == Local variables ==
_RL dflhdT ! derivative of latent heat with respect to T
_RL dfshdT ! derivative of sensible heat with respect to T
_RL dflwupdT ! derivative of long wave with respect to T
_RL Tsf ! surface temperature [K]
_RL Ts2 ! surface temperature square [K^2]
c _RL ht ! height for air temperature [m]
c _RL hq ! height for humidity [m]
c _RL hu ! height for wind speed [m]
c _RL zref ! reference height [m]
_RL wsm ! limited wind speed [m/s] (> umin)
_RL usn ! neutral, zref (=10m) wind speed [m/s]
_RL usm ! usn but limited [m/s] (> umin)
c _RL umin ! minimum wind speed used for drag-coeff [m/s]
_RL lath ! latent heat of vaporization or sublimation [J/kg]
_RL t0 ! virtual temperature [K]
_RL delth ! potential temperature diff [K]
_RL delq ! specific humidity difference [kg/kg]
_RL ustar ! friction velocity [m/s]
_RL tstar ! temperature scale [K]
_RL qstar ! humidity scale [kg/kg]
_RL rd ! = sqrt(Cd) [-]
_RL re ! = Ce / sqrt(Cd) [-]
_RL rh ! = Ch / sqrt(Cd) [-]
_RL rdn, ren, rhn ! neutral, zref (=10m) values of rd, re, rh
_RL stable ! = 1 if stable ; = 0 if unstable
_RL huol ! stability parameter at zwd [-] (=z/Monin-Obuklov length)
_RL htol ! stability parameter at zth [-]
_RL x ! stability function [-]
_RL xsq ! = x^2 [-]
_RL psimh ! momentum stability function
_RL psixh ! latent & sensib. stability function
_RL czol ! = zref*Karman_cst*gravity
_RL zwln ! = log(zwd/zref)
_RL ztln ! = log(zth/zref)
c _RL cdalton ! coeff to evaluate Dalton Number
c _RL mixratio
c _RL ea
c _RL psim_fac
_RL tau ! surface stress coef = rhoA * Ws * Cd
_RL csha ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir
_RL clha ! latent heat flx coef = rhoA * Ws * Ce * Lvap
c _RL zice
c _RL ssq0, ssq1, ssq2 ! constant used in saturated specific humidity
c _RL p0 ! reference sea-level atmospheric pressure [mb]
_RL qs1w, qs2w ! above freezing saturated specific humidity
_RL qs1i, qs2i ! below freezing saturated specific humidity
_RL tmpBlk
_RL half, one, two
INTEGER iter
C == external Functions
C-- Constant
DATA half, one, two
& / 0.5 _d 0 , 1. _d 0 , 2. _d 0 /
c DATA ssq0, ssq1, ssq2
c & / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 /
c DATA p0 / 1013. _d 0 /
DATA qs1w, qs2w
& / 640.38 _d 3 , 5107.0 _d -0 /
DATA qs1i, qs2i
& / 11637.80 _d 3 , 5897.8 _d -0 /
C-- Set surface parameters :
c zice = 0.0005 _d 0
zwln = LOG(zwd/zref)
ztln = LOG(zth/zref)
czol = zref*xkar*gravity
C- Surface Temp.
Tsf = tsfCel+Tf0kel
Ts2 = Tsf*Tsf
C- Wind speed
IF (ws.EQ.0. _d 0) THEN
ws = SQRT(uw*uw + vw*vw)
ENDIF
wsm = MAX(ws,umin)
C--- Compute turbulent surface fluxes
C- Pot. Temp and saturated specific humidity
t0 = Ta*(one + humid_fac*Qa)
IF ( iceornot.EQ.0 ) THEN
lath=Lvap
ssq = saltQsFac*qs1w*EXP( -qs2w/Tsf ) / rhoA
dEvdT = qs2w
ELSE
lath = Lvap+Lfresh
ssq = qs1i*EXP( -qs2i/Tsf ) / rhoA
dEvdT = qs2i
ENDIF
c ssq = ssq0*EXP( lath*(ssq1-ssq2/Tsf) ) / p0
c dEvdT = lath*ssq2
delth = Ta + gamma_blk*zth - Tsf
delq = Qa - ssq
C-- initial guess for exchange coefficients:
C take U_N = del.U ; stability from del.Theta ;
stable = half + SIGN(half, delth)
tmpBlk = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
rdn = SQRT(tmpBlk)
rhn = stable*cStantonS + (one-stable)*cStantonU
ren = cDalton
c rdn=xkar/(LOG(zref/zice))
c rhn=rdn
c ren=rdn
C-- calculate turbulent scales
ustar=rdn*wsm
tstar=rhn*delth
qstar=ren*delq
C--- iterate with psi-functions to find transfer coefficients
DO iter=1,blk_nIter
huol = ( tstar/t0
& +qstar/(Qa + one/humid_fac)
& )*czol/(ustar*ustar)
huol = SIGN( MIN(abs(huol),10. _d 0), huol)
stable = half + SIGN(half, huol)
xsq = SQRT( ABS(one - huol*16. _d 0) )
x = SQRT(xsq)
psimh = -5. _d 0*huol*stable
& + (one-stable)*
& ( LOG( (one + two*x + xsq)*(one+xsq)*.125 )
& -two*ATAN(x) + half*pi )
htol = huol*zth/zwd
xsq = SQRT( ABS(one - htol*16. _d 0) )
psixh = -5. _d 0*htol*stable
& + (one-stable)*( two*LOG(half*(one+xsq)) )
C- Shift wind speed using old coefficient
usn = ws/(one + rdn/xkar*(zwln-psimh) )
usm = MAX(usn, umin)
C- Update the 10m, neutral stability transfer coefficients
tmpBlk = cdrag_1/usm + cdrag_2 + cdrag_3*usm
rdn = SQRT(tmpBlk)
rhn = stable*cStantonS + (one-stable)*cStantonU
ren = cDalton
C- Shift all coefficients to the measurement height and stability.
rd = rdn/(1. _d 0 + rdn*(zwln-psimh)/xkar)
rh = rhn/(1. _d 0 + rhn*(ztln-psixh)/xkar)
re = ren/(1. _d 0 + ren*(ztln-psixh)/xkar)
C-- Update ustar, tstar, qstar using updated, shifted coefficients.
ustar = rd*wsm
qstar = re*delq
tstar = rh*delth
ENDDO
C- Coeff:
tau = rhoA*rd*ws
csha = cpAir*tau*rh
clha = lath*tau*re
C- Turbulent Fluxes
fsha = csha*delth
flha = clha*delq
evp = -flha/lath
ust = tau*rd*uw
vst = tau*rd*vw
C- surf.Temp derivative of turbulent Fluxes
dEvdT = tau*re*ssq*dEvdT/Ts2
dflhdT = -lath*dEvdT
dfshdT = -csha
C--- Upward long wave radiation
IF ( iceornot.EQ.0 ) THEN
flwupa = ocean_emissivity*stefan*Ts2*Ts2
dflwupdT= ocean_emissivity*stefan*Ts2*Tsf*4. _d 0
ELSEIF (iceornot.EQ.2) THEN
flwupa = snow_emissivity*stefan*Ts2*Ts2
dflwupdT = snow_emissivity*stefan*Ts2*Tsf*4. _d 0
ELSE
flwupa = ice_emissivity*stefan*Ts2*Ts2
dflwupdT = ice_emissivity*stefan*Ts2*Tsf*4. _d 0
ENDIF
C- Total derivative with respect to surface temperature
df0dT = -dflwupdT+dfshdT+dflhdT
#endif /*ALLOW_BULK_FORCE*/
RETURN
END