MITgcm |
|
|
|
|||
File indexing completed on 2018-03-02 18:38:13 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTCf664a6d8bb Jean*0001 #include "BULK_FORCE_OPTIONS.h" 0002 0003 CBOP 0004 C !ROUTINE: BULKF_FORMULA_LAY 0005 C !INTERFACE: 0006 SUBROUTINE BULKF_FORMULA_LAY( 0007 I uw, vw, ws, Ta, Qa, tsfCel, 0008 O flwupa, flha, fsha, df0dT, 0009 O ust, vst, evp, ssq, dEvdT, 0010 I iceornot, i,j,bi,bj,myThid ) 0011 0012 C !DESCRIPTION: \bv 0013 C *==========================================================* 0014 C | SUBROUTINE BULKF_FORMULA_LAY 0015 C | o Calculate bulk formula fluxes over open ocean or seaice 0016 C | Large and Yeager, 2004, NCAR/TN-460+STR. 0017 C *==========================================================* 0018 C \ev 0019 C 0020 C === Turbulent Fluxes : 0021 C * use the approach "B": shift coeff to height & stability of the 0022 C atmosphere state (instead of "C": shift temp & humid to the height 0023 C of wind, then shift the coeff to this height & stability of the atmos). 0024 C * similar to EXF (except over sea-ice) ; default parameter values 0025 C taken from Large & Yeager. 0026 C * assume that Qair & Tair inputs are from the same height (zq=zt) 0027 C * formulae in short: 0028 C wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v) 0029 C Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir 0030 C Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap 0031 C = -Evap * Lvap 0032 C with Ws = wind speed = sqrt(del.u^2 +del.v^2) ; 0033 C del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ; 0034 C Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number 0035 C respectively [no-units], function of height & stability 0036 0037 C !USES: 0038 IMPLICIT NONE 0039 C === Global variables === 0040 #include "EEPARAMS.h" 0041 #include "SIZE.h" 0042 #include "PARAMS.h" 0043 #include "BULKF_PARAMS.h" 0044 0045 C !INPUT/OUTPUT PARAMETERS: 0046 C input: 0047 _RL uw ! zonal wind speed (at grid center) [m/s] 0048 _RL vw ! meridional wind speed (at grid center) [m/s] 0049 _RL ws ! wind speed [m/s] at height zwd 0050 _RL Ta ! air temperature [K] at height zth 0051 _RL Qa ! specific humidity [kg/kg] at heigth zth 0052 _RL tsfCel ! sea-ice or sea surface temperature [oC] 0053 INTEGER iceornot ! 0=open water, 1=sea-ice, 2=sea-ice with snow 0054 INTEGER i,j, bi,bj !current grid point indices 0055 INTEGER myThid ! my Thread Id number 0056 C output: 0057 _RL flwupa ! upward long wave radiation (>0 upward) [W/m2] 0058 _RL flha ! latent heat flux (>0 downward) [W/m2] 0059 _RL fsha ! sensible heat flux (>0 downward) [W/m2] 0060 _RL df0dT ! derivative of heat flux with respect to Tsf [W/m2/K] 0061 _RL ust ! zonal wind stress (at grid center) [N/m2] 0062 _RL vst ! meridional wind stress (at grid center)[N/m2] 0063 _RL evp ! evaporation rate (over open water) [kg/m2/s] 0064 _RL ssq ! surface specific humidity [kg/kg] 0065 _RL dEvdT ! derivative of evap. with respect to Tsf [kg/m2/s/K] 0066 CEOP 0067 0068 #ifdef ALLOW_BULK_FORCE 0069 0070 C == Local variables == 0071 _RL dflhdT ! derivative of latent heat with respect to T 0072 _RL dfshdT ! derivative of sensible heat with respect to T 0073 _RL dflwupdT ! derivative of long wave with respect to T 0074 0075 _RL Tsf ! surface temperature [K] 0076 _RL Ts2 ! surface temperature square [K^2] 0077 c _RL ht ! height for air temperature [m] 0078 c _RL hq ! height for humidity [m] 0079 c _RL hu ! height for wind speed [m] 0080 c _RL zref ! reference height [m] 0081 _RL wsm ! limited wind speed [m/s] (> umin) 0082 _RL usn ! neutral, zref (=10m) wind speed [m/s] 0083 _RL usm ! usn but limited [m/s] (> umin) 0084 c _RL umin ! minimum wind speed used for drag-coeff [m/s] 0085 _RL lath ! latent heat of vaporization or sublimation [J/kg] 0086 _RL t0 ! virtual temperature [K] 0087 _RL delth ! potential temperature diff [K] 0088 _RL delq ! specific humidity difference [kg/kg] 0089 _RL ustar ! friction velocity [m/s] 0090 _RL tstar ! temperature scale [K] 0091 _RL qstar ! humidity scale [kg/kg] 0092 _RL rd ! = sqrt(Cd) [-] 0093 _RL re ! = Ce / sqrt(Cd) [-] 0094 _RL rh ! = Ch / sqrt(Cd) [-] 0095 _RL rdn, ren, rhn ! neutral, zref (=10m) values of rd, re, rh 0096 _RL stable ! = 1 if stable ; = 0 if unstable 0097 _RL huol ! stability parameter at zwd [-] (=z/Monin-Obuklov length) 0098 _RL htol ! stability parameter at zth [-] 0099 _RL x ! stability function [-] 0100 _RL xsq ! = x^2 [-] 0101 _RL psimh ! momentum stability function 0102 _RL psixh ! latent & sensib. stability function 0103 _RL czol ! = zref*Karman_cst*gravity 0104 _RL zwln ! = log(zwd/zref) 0105 _RL ztln ! = log(zth/zref) 0106 c _RL cdalton ! coeff to evaluate Dalton Number 0107 c _RL mixratio 0108 c _RL ea 0109 c _RL psim_fac 0110 _RL tau ! surface stress coef = rhoA * Ws * Cd 0111 _RL csha ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir 0112 _RL clha ! latent heat flx coef = rhoA * Ws * Ce * Lvap 0113 c _RL zice 0114 c _RL ssq0, ssq1, ssq2 ! constant used in saturated specific humidity 0115 c _RL p0 ! reference sea-level atmospheric pressure [mb] 0116 _RL qs1w, qs2w ! above freezing saturated specific humidity 0117 _RL qs1i, qs2i ! below freezing saturated specific humidity 0118 _RL tmpBlk 0119 _RL half, one, two 0120 INTEGER iter 0121 0122 C == external Functions 0123 0124 C-- Constant 0125 DATA half, one, two 0126 & / 0.5 _d 0 , 1. _d 0 , 2. _d 0 / 0127 c DATA ssq0, ssq1, ssq2 0128 c & / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 / 0129 c DATA p0 / 1013. _d 0 / 0130 DATA qs1w, qs2w 0131 & / 640.38 _d 3 , 5107.0 _d -0 / 0132 DATA qs1i, qs2i 0133 & / 11637.80 _d 3 , 5897.8 _d -0 / 0134 0135 C-- Set surface parameters : 0136 c zice = 0.0005 _d 0 0137 zwln = LOG(zwd/zref) 0138 ztln = LOG(zth/zref) 0139 czol = zref*xkar*gravity 0140 0141 C- Surface Temp. 0142 Tsf = tsfCel+Tf0kel 0143 Ts2 = Tsf*Tsf 0144 C- Wind speed 0145 IF (ws.EQ.0. _d 0) THEN 0146 ws = SQRT(uw*uw + vw*vw) 0147 ENDIF 0148 wsm = MAX(ws,umin) 0149 0150 C--- Compute turbulent surface fluxes 0151 C- Pot. Temp and saturated specific humidity 0152 t0 = Ta*(one + humid_fac*Qa) 0153 IF ( iceornot.EQ.0 ) THEN 0154 lath=Lvap 0155 ssq = saltQsFac*qs1w*EXP( -qs2w/Tsf ) / rhoA 0156 dEvdT = qs2w 0157 ELSE 0158 lath = Lvap+Lfresh 0159 ssq = qs1i*EXP( -qs2i/Tsf ) / rhoA 0160 dEvdT = qs2i 0161 ENDIF 0162 c ssq = ssq0*EXP( lath*(ssq1-ssq2/Tsf) ) / p0 0163 c dEvdT = lath*ssq2 0164 0165 delth = Ta + gamma_blk*zth - Tsf 0166 delq = Qa - ssq 0167 0168 C-- initial guess for exchange coefficients: 0169 C take U_N = del.U ; stability from del.Theta ; 0170 stable = half + SIGN(half, delth) 0171 tmpBlk = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm 0172 rdn = SQRT(tmpBlk) 0173 rhn = stable*cStantonS + (one-stable)*cStantonU 0174 ren = cDalton 0175 c rdn=xkar/(LOG(zref/zice)) 0176 c rhn=rdn 0177 c ren=rdn 0178 C-- calculate turbulent scales 0179 ustar=rdn*wsm 0180 tstar=rhn*delth 0181 qstar=ren*delq 0182 0183 C--- iterate with psi-functions to find transfer coefficients 0184 DO iter=1,blk_nIter 0185 0186 huol = ( tstar/t0 0187 & +qstar/(Qa + one/humid_fac) 0188 & )*czol/(ustar*ustar) 0189 huol = SIGN( MIN(abs(huol),10. _d 0), huol) 0190 stable = half + SIGN(half, huol) 0191 xsq = SQRT( ABS(one - huol*16. _d 0) ) 0192 x = SQRT(xsq) 0193 psimh = -5. _d 0*huol*stable 0194 & + (one-stable)* 0195 & ( LOG( (one + two*x + xsq)*(one+xsq)*.125 ) 0196 & -two*ATAN(x) + half*pi ) 0197 htol = huol*zth/zwd 0198 xsq = SQRT( ABS(one - htol*16. _d 0) ) 0199 psixh = -5. _d 0*htol*stable 0200 & + (one-stable)*( two*LOG(half*(one+xsq)) ) 0201 0202 C- Shift wind speed using old coefficient 0203 usn = ws/(one + rdn/xkar*(zwln-psimh) ) 0204 usm = MAX(usn, umin) 0205 0206 C- Update the 10m, neutral stability transfer coefficients 0207 tmpBlk = cdrag_1/usm + cdrag_2 + cdrag_3*usm 0208 rdn = SQRT(tmpBlk) 0209 rhn = stable*cStantonS + (one-stable)*cStantonU 0210 ren = cDalton 0211 0212 C- Shift all coefficients to the measurement height and stability. 0213 rd = rdn/(1. _d 0 + rdn*(zwln-psimh)/xkar) 0214 rh = rhn/(1. _d 0 + rhn*(ztln-psixh)/xkar) 0215 re = ren/(1. _d 0 + ren*(ztln-psixh)/xkar) 0216 0217 C-- Update ustar, tstar, qstar using updated, shifted coefficients. 0218 ustar = rd*wsm 0219 qstar = re*delq 0220 tstar = rh*delth 0221 0222 ENDDO 0223 0224 C- Coeff: 0225 tau = rhoA*rd*ws 0226 csha = cpAir*tau*rh 0227 clha = lath*tau*re 0228 0229 C- Turbulent Fluxes 0230 fsha = csha*delth 0231 flha = clha*delq 0232 evp = -flha/lath 0233 ust = tau*rd*uw ec561b83af Mart*0234 vst = tau*rd*vw f664a6d8bb Jean*0235 0236 C- surf.Temp derivative of turbulent Fluxes 0237 dEvdT = tau*re*ssq*dEvdT/Ts2 0238 dflhdT = -lath*dEvdT 0239 dfshdT = -csha 0240 0241 C--- Upward long wave radiation 0242 IF ( iceornot.EQ.0 ) THEN 0243 flwupa = ocean_emissivity*stefan*Ts2*Ts2 0244 dflwupdT= ocean_emissivity*stefan*Ts2*Tsf*4. _d 0 0245 ELSEIF (iceornot.EQ.2) THEN 0246 flwupa = snow_emissivity*stefan*Ts2*Ts2 0247 dflwupdT = snow_emissivity*stefan*Ts2*Tsf*4. _d 0 0248 ELSE 0249 flwupa = ice_emissivity*stefan*Ts2*Ts2 0250 dflwupdT = ice_emissivity*stefan*Ts2*Tsf*4. _d 0 0251 ENDIF 0252 0253 C- Total derivative with respect to surface temperature 0254 df0dT = -dflwupdT+dfshdT+dflhdT 0255 0256 #endif /*ALLOW_BULK_FORCE*/ 0257 0258 RETURN 0259 END
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine. The LXR team |
|