C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lanl.F,v 1.8 2006/05/25 17:30:54 jmc Exp $ C $Name: $ #include "BULK_FORCE_OPTIONS.h" CBOP C !ROUTINE: BULKF_FORMULA_LANL C !INTERFACE: SUBROUTINE BULKF_FORMULA_LANL( I uw, vw, us, Ta, Qa, nc, tsf_in, O flwupa, flha, fsha, df0dT, O ust, vst, evp, ssq, dEvdT, I iceornot, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE BULKF_FORMULA_LANL c | o Calculate bulk formula fluxes over open ocean or seaice C *==========================================================* C \ev C swd -- bulkf formula used in bulkf and ice pkgs C taken from exf package C 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 = transfer coefficient for momentum, sensible C & latent heat flux [no units] C *==========================================================* 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 us ! wind speed [m/s] at height hu _RL Ta ! air temperature [K] at height ht _RL Qa ! specific humidity [kg/kg] at heigth ht _RL nc ! fraction cloud cover _RL tsf_in ! sea-ice or sea surface temperature [oC] INTEGER iceornot ! 0=open water, 1=sea-ice, 2=sea-ice with snow 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 ht ! height for air temperature [m] c _RL hq ! height for humidity [m] _RL hu ! height for wind speed [m] c _RL zref ! reference height [m] _RL usm ! wind speed limited [m/s] c _RL umin ! minimum wind speed used for drag-coeff [m/s] _RL lath ! latent heat of vaporization or sublimation _RL t0 ! virtual temperature [K] _RL deltap ! 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 ! initial (neutral) values of rd, re, rh _RL stable ! = 1 if stable ; = 0 if unstable _RL huol ! stability parameter [-] _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 aln ! = log(ht/zref) c _RL cdalton ! coeff to evaluate Dalton Number c _RL mixratio c _RL ea c _RL psim_fac _RL tau ! surface stress coef = ? _RL csha ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir _RL clha ! latent heat flx coef = rhoA * Ws * Ce * Lvap _RL zice _RL ssq0, ssq1, ssq2 ! constant used in surface specific humidity _RL p0 ! reference sea-level atmospheric pressure [mb] _RL bulkf_Cdn ! drag coefficient INTEGER niter_bulk, iter C == external Functions c _RL exf_BulkCdn c external exf_BulkCdn c _RL exf_BulkqSat c external exf_BulkqSat c _RL exf_BulkRhn c external exf_BulkRhn DATA ssq0, ssq1, ssq2 & / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 / DATA p0 / 1013. _d 0 / C--- Compute turbulent surface fluxes ht = 2. _d 0 c hq = 2. _d 0 hu = 10. _d 0 c zref = 10. _d 0 zice = 0.0005 _d 0 aln = log(ht/zref) niter_bulk = 5 c cdalton = 0.0346000 _d 0 czol = zref*xkar*gravity c psim_fac=5. _d 0 c umin=1. _d 0 lath=Lvap if (iceornot.gt.0) lath=Lvap+Lfresh Tsf=Tsf_in+Tf0kel C- Wind speed if (us.eq.0. _d 0) then us = sqrt(uw*uw + vw*vw) endif usm = max(us,umin) c t0 = Ta*(1. _d 0 + humid_fac*Qa) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| cQQ ssq = 0.622*6.11*exp(22.47*(1.d0-Tf0kel/tsf))/p0 c ssq = 3.797915 _d 0*exp( c & lath*(7.93252 _d -6 - 2.166847 _d -3/Tsf) c & )/p0 ssq = ssq0*exp( lath*(ssq1-ssq2/Tsf) ) / p0 deltap = ta - tsf + gamma_blk*ht delq = Qa - ssq C-- initialize estimate exchange coefficients rdn=xkar/(log(zref/zice)) rhn=rdn ren=rdn C-- calculate turbulent scales ustar=rdn*usm tstar=rhn*deltap qstar=ren*delq C-- iteration with psi-functions to find transfer coefficients do iter=1,niter_bulk huol = czol/ustar**2 *(tstar/t0 + & qstar/(1. _d 0/humid_fac+Qa)) huol = sign( min(abs(huol),10. _d 0), huol) stable = 5. _d -1 + sign(5. _d -1 , huol) xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0) x = sqrt(xsq) psimh = -5. _d 0*huol*stable + (1. _d 0-stable)* & (2. _d 0*log(5. _d -1*(1. _d 0+x)) + & 2. _d 0*log(5. _d -1*(1. _d 0+xsq)) - & 2. _d 0*atan(x) + pi*.5 _d 0) psixh = -5. _d 0*huol*stable + (1. _d 0-stable)* & (2. _d 0*log(5. _d -1*(1. _d 0+xsq))) C-- Update the transfer coefficients rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar) rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar) re = rh C-- Update ustar, tstar, qstar using updated, shifted coefficients. ustar = rd*usm qstar = re*delq tstar = rh*deltap enddo tau = rhoa*ustar**2 tau = tau*us/usm csha = rhoa*cpair*us*rh*rd clha = rhoa*lath*us*re*rd fsha = csha*deltap flha = clha*delq evp = -flha/lath C-- Upward long wave radiation cQQ mixratio=Qa/(1-Qa) cQQ ea=p0*mixratio/(0.62197+mixratio) cQQ flwupa=-0.985*stefan*tsf**4 cQQ & *(0.39-0.05*sqrt(ea)) cQQ & *(1-0.6*nc**2) if (iceornot.eq.0) then flwupa=ocean_emissivity*stefan*tsf**4 dflwupdT=4. _d 0*ocean_emissivity*stefan*tsf**3 elseif (iceornot.eq.2) then flwupa=snow_emissivity*stefan*tsf**4 dflwupdT=4. _d 0*snow_emissivity*stefan*tsf**3 else flwupa=ice_emissivity*stefan*tsf**4 dflwupdT=4. _d 0*ice_emissivity*stefan*tsf**3 endif cQQ dflhdT = -clha*Tf0kel*ssq*22.47/(tsf**2) c dflhdT = -clha*Lath*ssq/(Rvap*tsf**2) c dflhdT = -clha*ssq*Lath*2.166847 _d -3/(Tsf**2) dEvdT = clha*ssq*ssq2/(Tsf*Tsf) dflhdT = -lath*dEvdT dfshdT = -csha cQQ dflwupdT= 4.*0.985*stefan*tsf**3 cQQ & *(0.39-0.05*sqrt(ea)) cQQ & *(1-0.6*nc**2) c total derivative with respect to surface temperature df0dT=-dflwupdT+dfshdT+dflhdT C-- Wind stress at center points C- in-lining of function: exf_BulkCdn(umps) = cdrag_1/umps + cdrag_2 + cdrag_3*umps bulkf_Cdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm ust = rhoa*bulkf_Cdn*us*uw vst = rhoa*bulkf_Cdn*us*vw #endif /*ALLOW_BULK_FORCE*/ RETURN END