PHYSICS.050 Documentation for Source Terms for Runs A087 1999/09/10 Energy Reference Level X is used: static energy of dry air is measured from 0øK, but energy of liquid water is measured from 0øC. Constants --------- ELHE (J/kg) = latent heat of vaporization at 0øC = 2500000 ELHM (J/kg) = latent heat of melting at 0øC = 334000 SHCD (J/kgøC) = specific heat capacity of dry air = 1003.5 SHCV (J/kgøC) = specific heat capacity of water vapor = 1911 SHCW (J/kgøC) = specific heat capacity of liquid water = 4185 SHCI (J/kgøC) = specific heat capacity of ice = 2060 RDRY (J/kgøC) = gas constant for dry air = 287 RVAP (J/kgøC) = gas constant for water vapor = 461.5 RHOI (kg/m^3) = density of ice = 916.6 RHOS (kg/m^3) = density of snow = 100. TKF (K) = freexing point of water = 273.16 ALMAI (W/møC) = thermal conductivity of ice = 2.1762 ALMAS (W/møC) = thermal conductivity of snow = .35 ACE1I (kg/mý) = ice thickness of mass layer 1 = 182. PPVSF (Pa) = saturated partial pressure of water vapor at 0øC = 610.571 S1byG1 = ratio of surface level from first layer to ground = RRT3 WSBYW1 = ratio of surface wind speed to first layer speed = .75 CIAX = surface wind factor in cross isobar angle = .2 WSADD (m/s) = addition to wind speed used in surface fluxes = 2 DXYP (mý) = horizontal area of grid box dt (s) = source term time step XSI1,XSI2 = fractions of mass layer 1 assigned to thermal layers 1,2 RKAP = exponent of exner function = RDRY/SHCD RRT3 = 1/û3 Prognostic Variables -------------------- MA1 (kg/mý) = mass of lowest atmospheric layer U1 (m/s) = zonal velocity of lowest atmospheric layer V1 (m/s) = meridional velocity of lowest atmospheric layer H0M1 (J) = mean potential enthalpy of lowest atmospheric layer HZM1 (J) = vertical gradient of potential enthalpy of lowest layer Q0M1 (kg) = mean water vapor of lowest atmospheric layer QZM1 (kg) = vertical gradient of water vapor of lowest atmospheric layer MSI1 (kg/mý) = mass of ice and snow of layers 1 and 2 HSI1 (J/mý) = enthalpy minus latent heat of layer 1 HSI2 (J/mý) = enthalpy minus latent heat of layer 2 Derived variables ----------------- MSUM (kg/mý) = total column mass including MSTRAT P1 (Pa) = GRAV*(MSUM-MA1/2) = pressure at center of layer 1 PS (Pa) = GRAV*[MSUM-MA1*(1-S1byG1)/2] = pressure of surface level PG (Pa) = GRAV*MSUM = pressure of ground PPV (Pa) = RVAP*MV*PR / (RDRY*MD+RVAP*MV) = partial pressure of water vapor PPVG (Pa) = PPVSF exp [ELHV*(1/TKF - 1/TKG)/RVAP] = PPVSAT(TKG,ELHV) PK1 = P1^RKAP PKUP = PK1 - .25*RKAP*GRAV*MA1*P1^(RKAP-1) PKDN = PK1 + .25*RKAP*GRAV*MA1*P1^(RKAP-1) PGK = PG^RKAP EMUP (J) = .5*H0M*PKUP + .25*HZM*PK = sensible heat of upper half air layer EMDN (J) = .5*H0M*PKDN - .25*HZM*PK = sensible heat of lower half air layer QG = MV / (MD+MV) = RDRY*PPVG / [RVAP*PG+(RDRY-RVAP)*PPVG] QS = (Q0M1-QZM1*S1byG1) / MA1*DXYP = specific humidity of surf. HS (J/kg) = (H0M1-HZM1*S1byG1)*PGK / MA1*DXYP = pot. spec. enthalpy of surf. TKS (øK) = (H0M1-HZM1*S1byG1)*PGK / MA1*DXYP*SHCS = surface air temperature TCG1 (øC) = (HSI1/.5*MSI1 + ELHM) / SHCI TCG2 (øC) = (HSI2/.5*MSI1 + ELHM) / SHCI TKG (øK) = TCG1 + TKF HCG1 (J/øCmý) = SHCI*XSI1*MSI1 = heat capacity of sea ice layer 1 HCG2 (J/øCmý) = SHCI*XSI2*MSI1 = heat capacity of sea ice layer 2 SNOW (kg/mý) = MSI1 - ACE1I = mass of snow SHCS (J/kgøC) = SHCD + QS *(SHCV-SHCD) = specific heat capacity of surf. air ELHV (J/kg) = ELHE + TCG*(SHCV-SHCW) = latent heat of vaporization ELHV (J/kg) = ELHS + TCG*(SHCV-SHCI) = minus latent heat of sublimation ELHQ (J/kg) = ELHE + TCG* SHCV = specific energy of atmo. water vapor ELHG (J/kg) = ELHE + TCG* SHCV = specific enthalpy due to vapor ELHG (J/kg) = ELHS + TCG* SHCV = specific enthalpy due to vapor W1SQ (mý/sý) = U1*U1 + V1*V1 WSSQ (mý/sý) = W1SQ*WSBYW1*WSBYW1 WS (m/s) = WSSQ^.5 BETAS (m^3/kg) = [RDRY+(RVAP-RDRY)*QS]*TKS / PG = pot. spec. vol. of surf. BETAG (m^3/kg) = [RDRY+(RVAP-RDRY)*QG*WET]*TKG / PG = pot. spec. vol. of grnd. RIGS = (PG-PS)*(BETAS-BETAG) / WSSQ CDN = .0009 + .0000025*WSSQ (over open water) CDN = .0024 (over sea ice) CDN = .008 - .035 (over land) CDM = CDH = CDN*2/(1+4*RIGSý) CIA = cross isobar angle between first layer wind and surface wind = = TWOPI*[.09375-.03125/(1+4*RIGSý)] / (1+WS*CIAX) COSC = cos(CIA) SINC = sin(CIA)*HEMI US (m/s) = (U1*COSC-V1*SINC)*WSBYW1 VS (m/s) = (V1*COSC+U1*SINC)*WSBYW1 Energy ------ SE (J/mý) = static energy GE (J/mý) = geopotential energy SE-GE = [SHCD*MD + SHCV*MV]*TK + [SHCW*ML + SHCI*MI]*TC + + [ELHE - SHCV*TKF]*MV - ELHM*MI Division of SE-GE into components --------------------------------- ELVA (J/mý) = atmos. latent energy of vapor = [ELHE - SHCV*TKF]*MV ELIA (J/mý) = atmos. latent energy of ice = - ELHM*MI EKA (J/mý) = atmos. gaseous enthalpy = [SHCD*MD + SHCV*MV]*TK ECA (J/mý) = atmos. condensate enthalpy = [SHCW*ML + SHCI*MI]*TC H0M (J) = atmos. potential enthalpy = EKA*DXYP / PK ELIG (J/mý) = ground latent energy of ice = - ELHM*MI ECG (J/mý) = ground condensate enthalpy = [SHCW*ML + SHCI*MI]*TC Recreation of H0M and HZM from EMUP and EMDN -------------------------------------------- EMUP + EMDN = .5*H0M*PKUP + .25*HZM*PK + .5*H0M*PKDN - .25*HZM*PK = H0M*PK H0M = (EMUP+EMDN)/PK EKUP*PKDN - EKDN*PKUP = = (.5*H0M*PKUP+.25*HZM*PK)*PKDN - (.5*H0M*PKDN-.25*HZM*PK)*PKUP = = .25*HZM*PK*(PKDN+PKUP) = .5*HZM*PKý HZM = 2*(EKUP*PKDN - EKDN*PKUP) / PKý Change in H0M and HZM due to radiation -------------------------------------- RHR (W/mý) = radiation heating rate which is constant in the vertical = = SRHR*COSZ + TRHR EMUPnew = EMUPold + .5*dT*RHR*DXYP EMDNnew = EMDNold + .5*dT*RHR*DXYP H0Mnew = (EMUPnew+EMDNnew)/PK = H0Mold + dt*RHR*DXYP/PK HZMnew = 2*(EKUPnew*PKDN - EKDNnew*PKUP)/PKý = = 2*(EKUPold*PKDN - EKDNold*PKUP)/PKý + dt*RHR*DXYP*(PKDN-PKUP)/PKý = = HZMold + dt*RHR*DXYP*(.5*RKAP*GRAV*MA*PKm1)/PKý Surface Fluxes -------------- WVF (kg/mýs) = (QS-QG*WET)*CDH*(WS+WSADD)/BETAS = water vapor or dew flux SRH (W/mý) = solar radiation absorbed by ground TRH (W/mý) = TRDN - STBO*TKG^4 = thermal radiation absorbed by ground SNH (W/mý) = (TKS-TKG)*SHCS*CDH*(WS+WSADD)/BETAS = sensible heat flux EVH (W/mý) = ELHQ*WVF = water vapor heat flux GVH (W/mý) = ELHG*WVF = ground enthalpy change due to dew F1 (W/mý) = (TCG1-TCG2) / RRT3*(ACE1I/RHOI*ALAMI+SNOW/RHOS*ALAMS) TAUX (N/mý) = US*CDM*(WS+WSADD)/BETAS = momentum stress TAUY (N/mý) = VS*CDM*(WS+WSADD)/BETAS = momentum stress Derivatives ----------- dQG/dTG ÷ [QG*ELHV / RVAP*TKGý] * [1 - (RDRY-RVAP)*PPVG/RVAP*PG] ÷ ÷ QG*ELHV / RVAP*TKGý (see calculation of dQG/dTG later) dWV/dQS = CDH*(WS+WSADD)/BETAS dWV/dQG = - WET*CDH*(WS+WSADD)/BETAS dWV/dTG = dWV/dQG * dQG/dTG dTR/dTG = - 4*STBO*TKG^3 dSN/dTG = - SHCS*CDH*(WS+WSADD)/BETAS dSN/dTS = SHCS*CDH*(WS+WSADD)/BETAS = - dSN/dTG dF1/dTG = 1 / .5*(ACE1I/RHOI*ALAMI+SNOW/RHOS*ALAMS) dF1/dT2 = -1 / .5*(ACE1I/RHOI*ALAMI+SNOW/RHOS*ALAMS) = - dF1/dTG Explicit Time Integration ------------------------- EKUP = EKUP EKDN = EKDN - dt*(SNH+SHCV*TKG*WVF)*DXYP H0M1 = H0M1 - dt*(SNH+SHCV*TKG*WVF)*DXYP / PK HZM1 = HZM1 + dt*(SNH+SHCV*TKG*WVF)*DXYP*2? / PK HCG1new = HCG1old + dt*SHCI*WVF TCG1new*HCG1new = TCG1old*HCG1old + dt*(SRH+TRH+SNH+GVH-F1) TCG1new*HCG1old = TCG1old*HCG1old + dt*(SRH+TRH+SNH+ELHG*WVF-F1) - - TGC1new*dt*SHCI*WVF ÷ ÷ TCG1old*HCG1old + dt*[SRH+TRH+SNH+(ELHS+SHCV*TCG)*WVF-F1] - - TGC1old*dt*SHCI*WVF = = TCG1*HCG1 + dt*(SRH+TRH+SNH+ELHV*WVF-F1) TCG1 = TCG1 + dt*(SRH+TRH+SNH+ELHV*WVF-F1) / HCG1 TCG2 = TCG2 + dt*F1 / HCG2 TKS ÷ TKS - dt*SNH*(1+2*S1byG1)*PGK / MA1*PK*SHCS Q0M1 = Q0M1 - dt*WVF*DXYP QZM1 = QZM1 + dt*WVF*DXYP*2? QS = QS - dt*WVF*(1+2*S1byG1) / MA1 U1 = U1 - dt*TAUX/MA1 V1 = V1 - dt*TAUY/MA1 Implicit Time Integration (a=0: explicit; a=1: fully implicit) ------------------------- dTS*MA1*SHCS*PK = -(1+2*S1byG1)*dt*PGK*[SNH + a*dTG*dSN/dTG + a*dTS*dSN/dTS] dQS*MA1 = -(1+2*S1byG1)*dt* [WVF + a*dTG*dWV/dTG + a*dQS*dWV/dQS] dT2*HCG2 = dt*[F1 + a*dTG*dF1/dTG + a*dT2*dF1/dT2] dTS*[MA1*SHCS*PK + (1+2*S1byG1)*dt*PGK*a*dSN/dTS] = -(1+2*S1byG1)*dt*PGK*[SNH + a*dTG*dSN/dTG] dQS*[MA1 + (1+2*S1byG1)*dt*a*dWV/dQS] = -(1+2*S1byG1)*dt*[WVF + a*dTG*dWV/dTG] dT2*[HCG2 - dt*a*dF1/dT2] = dt*[F1 + a*dTG*dF1/dTG] TSDEN = (1+2*S1byG1)*dt*a*PGK*dSN/dTS + MA1*SHCS*PK = = -(1+2*S1byG1)*dt*a*PGK*dSN/dTG + MA1*SHCS*PK QSDEN = (1+2*S1byG1)*dt*a*dWV/dQS + MA1 TSCON = -(1+2*S1byG1)*dt*SNH*PGK / TSDEN QSCON = -(1+2*S1byG1)*dt*WVF / QSDEN TSMUL = -(1+2*S1byG1)*dt*a*PGK*dSN/dTG / TSDEN QSMUL = -(1+2*S1byG1)*dt*a*dWV/dTG] / QSDEN T2DEN = HCG2 - dt*a*dF1/dT2 = HCG2 + dt*a*dF1/dTG T2CON = dt*F1 / T2DEN T2MUL = dt*a*dF1/dTG / T2DEN dTS = TSCON + dTG*TSMUL dQS = QSCON + dTG*QSMUL dT2 = T2CON + dTG*T2MUL dTG*HCG1 = dt*[SRH+TRH+SNH+ELHV*WVF-F1 + + a*dTG*(dTR/dTG+dSN/dTG+ELHV*dWV/dTG-dF1/dTG) + + a*dTS*dSN/dTS + a*dQS*ELHV*dWV/dQS - a*dT2*dF1/dT2] dTG*[HCG1 - dt*a*(dTR/dTG+dSN/dTG+ELHV*dWV/dTG-dF1/dTG)] = = dt*[SRH+TRH+SNH+ELHV*WVF-F1 + a*dSN/dTS*(TSCON+dTG*TSMUL) + + a*ELHV*dWV/dQS*(QSCON+dTG*QSMUL) - a*dF1/dT2*(T2CON+dTG*T2MUL)] dTG*[HCG1 - dt*a*(dTR/dTG+dSN/dTG+ELHV*dWV/dTG-dF1/dTG + + dSN/dTS*TSMUL + ELHV*dWV/dQS*QSMUL - dF1/dT2*T2MUL)] = = dt*[SRH+TRH+SNH+ELHV*WVF-F1 + + a*dSN/dTS*TSCON + a*ELHV*dWV/dQS*QSCON - a*dF1/dT2*T2CON] dTG*[HCG1 - dt*a*(dTR/dTG+dSN/dTG+ELHV*dWV/dTG-dF1/dTG - - TSMUL*dSN/dTG + QSMUL*ELHV*dWV/dQS + T2MUL*dF1/dTG)] = dt*[SRH+TRH+SNH+ELHV*WVF-F1 - - a*TSCON*dSN/dTG + a*QSCON*ELHV*dWV/dQS + a*T2CON*dF1/dTG] TRH*dT = dt*[TRH + a* dTG*dTR/dTG] SNH*dt = dt*[SNH + a*(dTG*dSN/dTG + dTS*dSN/dTS)] = = dt*[SNH + a*(dTG - dTS)*dSN/dTG] WVF*dt = dt*[WVF + a*(dTG*dWV/dTG + dQS*dWV/dQS)] = DEW = - EVAP F1 *dt = dt*[F1 + a*(dTG*dF1/dTG + dT2*dF1/dT2)] = = dt*[F1 + a*(dTG - dT2)*dF1/dTG] Approximation of dQG/dTG ------------------------ dPPVG/dTG = PPVG*ELHV / RVAP*TKGý dQG/dPPVG = {[RVAP*PG + (RDRY-RVAP)*PPVG]*RDRY - RDRY*PPVG*(RDRY-RVAP)] / [RVAP*PG + (RDRY-RVAP)*PPVG]ý = = RDRY*RVAP*PG / [RVAP*PG + (RDRY-RVAP)*PPVG]ý ÷ ÷ RDRY*RVAP*PG / [RVAPý*PGý + 2*RVAP*PG*(RDRY-RVAP)*PPVG] = = RDRY / [RVAP*PG + 2*(RDRY-RVAP)*PPVG] = = RDRY/RVAP*PG / [1 + 2*(RDRY-RVAP)*PPVG/RVAP*PG] ÷ ÷ RDRY/RVAP*PG * [1 - 2*(RDRY-RVAP)*PPVG/RVAP*PG] = ÷ RDRY/RVAP*PG * [1 - (RDRY-RVAP)*PPVG/RVAP*PG] - - RDRY*(RDRY-RVAP)*PPVG/RVAPý*PGý ÷ ÷ RDRY / RVAP*PG*[1 - (RDRY-RVAP)*PPVG/RVAP*PG] - - RDRY*(RDRY-RVAP)*PPVG/RVAPý*PGý = ÷ RDRY / [RVAP*PG* - (RDRY-RVAP)*PPVG] - - RDRY*(RDRY-RVAP)*PPVG/RVAPý*PGý = = QG/PPVG - RDRY*(RDRY-RVAP)*PPVG/RVAPý*PGý ÷ ÷ QG/PPVG - QG*(RDRY-RVAP)/RVAP*PG dQG/dTG = dQG/dPPVG * dPPVG/dTG ÷ ÷ [QG/PPVG - QG*(RDRY-RVAP)/RVAP*PG] * [PPVG*ELHV / RVAP*TKGý] = = [ELHV*QG / RVAP*TKGý] * [1 - (RDRY-RVAP)*PPVG/RVAP*PG] Components of downward surface fluxes from atmosphere to ground --------------------------------------------------------------- DEW (kg/mý) = water vapor flux from atmosphere to ground SNHDT (J/mý) = dry heating from atmosphere to ground = dECG = - dEKA EVHDT (J/mý) = energy to condense water vapor to liquid and drop to ground = = ELHV*DEW + SHCW*TCG*DEW = (ELHE + SHCV*TCG)*DEW Over liquid: dECG = SNHDT + (ELHE + SHCV*TCG)*DEW = - dELVA - dEKA dELVA = - (ELHE - SHCV*TKF)*DEW dEKA = - SNHDT - SHCV*TKG*DEW Over ice: dELIG + dECG = SNHDT + (ELHE + SHCV*TCG)*DEW = - dELVA - dEKA dELIG = - ELHM*DEW dECG = SNHDT + (ELHS + SHCV*TCG)*DEW dELVA = - (ELHE - SHCV*TKF)*DEW dEKA = - SNHDT - SHCV*TKG*DEW Update of atmospheric variables due to surface fluxes ----------------------------------------------------- MA1new = MA1old - DEW Q0M1new = Q0M1old - DEW*DXYP QZM1new = QZM1old + DEW*DXYP*2 H0M1new = (H0Mold*PKold/DXYP - SNHDT - SHCV*TKG*DEW)*DXYP / PKnew HZM1new = HZM1old + SNHDT*DXYP*2 / PKnew