subroutine INLANDSIS(SnoMod,BloMod,jjtime) USE dimphy !--------------------------------------------------------------------------+ ! INLANDSIS module | ! Simplified SISVAT module, containing ice and snow processes for | ! ice-covered surfaces | ! version MARv3, november 2020 | ! SubRoutine INLANDSIS contains the fortran 77 code of the | ! Soil/Ice Snow Vegetation Atmosphere Transfer Scheme | ! | !--------------------------------------------------------------------------+ ! PARAMETERS: klonv: Total Number of columns = | ! ^^^^^^^^^^ = Total Number of continental grid boxes | ! X Number of Mosaic Cell per grid box | ! | ! INPUT: daHost : Date Host Model | ! ^^^^^ | ! | ! INPUT: LSmask : 1: Land MASK | ! ^^^^^ 0: Sea MASK | ! isotSV = 0,...,12: Soil Type | ! 0: Water, Liquid (Sea, Lake) | ! 12: Water, Solid (Ice) | ! | ! INPUT: coszSV : Cosine of the Sun Zenithal Distance [-] | ! ^^^^^ sol_SV : Surface Downward Solar Radiation [W/m2] | ! IRd_SV : Surface Downward Longwave Radiation [W/m2] | ! drr_SV : Rain Intensity [kg/m2/s] | ! dsn_SV : Snow Intensity [mm w.e./s] | ! dsnbSV : Snow Intensity, Drift Fraction [-] | ! dbs_SV : Drift Amount [mm w.e.] | ! za__SV : Surface Boundary Layer (SBL) Height [m] | ! VV__SV :(SBL Top) Wind Velocity [m/s] | ! TaT_SV : SBL Top Temperature [K] | ! rhT_SV : SBL Top Air Density [kg/m3] | ! QaT_SV : SBL Top Specific Humidity [kg/kg] | ! qsnoSV : SBL Mean Snow Content [kg/kg] | ! alb0SV : Soil Basic Albedo [-] | ! slopSV : Surface Slope [-] | ! dt__SV : Time Step [s] | ! | ! INPUT / isnoSV = total Nb of Ice/Snow Layers | ! OUTPUT: ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer | ! ^^^^^^ iiceSV = total Nb of Ice Layers | ! istoSV = 0,...,5 : Snow History (see istdSV data) | ! | ! INPUT / alb_SV : Surface Albedo [-] | ! OUTPUT: emi_SV : Surface Emissivity [-] | ! ^^^^^^ IRs_SV : Soil IR Flux (negative) [W/m2] | ! LMO_SV : Monin-Obukhov Scale [m] | ! us__SV : Friction Velocity [m/s] | ! uts_SV : Temperature Turbulent Scale [m/s] | ! uqs_SV : Specific Humidity Velocity [m/s] | ! uss_SV : Blowing Snow Turbulent Scale [m/s] | ! usthSV : Blowing Snow Erosion Threshold [m/s] | ! Z0m_SV : Momentum Roughness Length [m] | ! Z0mmSV : Momentum Roughness Length (time mean) [m] | ! Z0mnSV : Momentum Roughness Length (instantaneous)[m] | ! Z0SaSV : Sastrugi Roughness Length [m] | ! Z0e_SV : Erosion Snow Roughness Length [m] | ! Z0emSV : Erosion Snow Roughness Length (time mean) [m] | ! Z0enSV : Erosion Snow Roughness Length (instantaneous)[m] | ! Z0roSV : Subgrid Topo Roughness Length [m] | ! Z0h_SV : Heat Roughness Length [m] | ! TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)| ! & Snow Temperatures (layers 1,2,...,nsno) [K] | ! ro__SV : Soil/Snow Volumic Mass [kg/m3] | ! eta_SV : Soil/Snow Water Content [m3/m3] | ! G1snSV : snow dendricity/sphericity | ! G2snSV : snow sphericity/grain size | ! dzsnSV : Snow Layer Thickness [m] | ! agsnSV : Snow Age [day] | ! BufsSV : Snow Buffer Layer [kg/m2] .OR. [mm] | ! BrosSV : Snow Buffer Layer Density [kg/m3] | ! BG1sSV : Snow Buffer Layer Dendricity / Sphericity [-] | ! BG2sSV : Snow Buffer Layer Sphericity / Size [-] [0.1 mm] | ! rusnSV : Surficial Water [kg/m2] .OR. [mm] | ! | ! OUTPUT: no__SV : OUTPUT file Unit Number [-] | ! ^^^^^^ i___SV : OUTPUT point i Coordinate [-] | ! j___SV : OUTPUT point j Coordinate [-] | ! n___SV : OUTPUT point n Coordinate [-] | ! lwriSV : OUTPUT point vec Index [-] | ! | ! OUTPUT: IRu_SV : Upward IR Flux (+, upw., effective) [K] | ! ^^^^^^ hSalSV : Saltating Layer Height [m] | ! qSalSV : Saltating Snow Concentration [kg/kg] | ! RnofSV : RunOFF Intensity [kg/m2/s] | ! | ! Internal Variables: | ! ^^^^^^^^^^^^^^^^^^ | ! NLaysv = New Snow Layer Switch [-] | ! albisv : Snow/Ice/Water/Soil Integrated Albedo [-] | ! SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] | ! TBr_sv : Brightness Temperature [K] | ! IRupsv : Upward IR Flux (-, upw.) [W/m2] | ! ram_sv : Aerodynamic Resistance for Momentum [s/m] | ! rah_sv : Aerodynamic Resistance for Heat [s/m] | ! Evp_sv : Evaporation [kg/m2] | ! EvT_sv : Evapotranspiration [kg/m2] | ! HSs_sv : Surface Sensible Heat Flux + => absorb.[W/m2] | ! HLs_sv : Surface Latent Heat Flux + => absorb.[W/m2] | ! Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] | ! Tsrfsv : Surface Temperature [K] | ! sEX_sv : Verticaly Integrated Extinction Coefficient [-] | ! LSdzsv : Vertical Discretization Factor [-] | ! = 1. Soil | ! = 1000. Ocean | ! z_snsv : Snow Pack Thickness [m] | ! zzsnsv : Snow Pack Thickness [m] | ! albssv : Soil Albedo [-] | ! Eso_sv : Soil+Snow Emissivity [-] | ! Khydsv : Soil Hydraulic Conductivity [m/s] | ! | ! ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] | ! ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] | ! ETSo_d : Snow/Soil Energy Power Forcing [W/m2] | ! EqSn_0 : Snow Energy, before Phase Change [J/m2] | ! EqSn_1 : Snow Energy, after Phase Change [J/m2] | ! EqSn_d : Snow Energy, net Forcing [J/m2] | ! Enrsvd : SVAT Energy Power Forcing [W/m2] | ! Enrbal : SVAT Energy Balance [W/m2] | ! Wats_0 : Soil Water, before Forcing [mm] | ! Wats_1 : Soil Water, after Forcing [mm] | ! Wats_d : Soil Water Forcing [mm] | ! SIWm_0 : Snow initial Mass [mm w.e.] | ! SIWm_1 : Snow final Mass [mm w.e.] | ! SIWa_i : Snow Atmos. initial Forcing [mm w.e.] | ! SIWa_f : Snow Atmos. final Forcing(noConsumed)[mm w.e.] | ! SIWe_i : SnowErosion initial Forcing [mm w.e.] | ! SIWe_f : SnowErosion final Forcing(noConsumed)[mm w.e.] | ! SIsubl : Snow sublimed/deposed Mass [mm w.e.] | ! SImelt : Snow Melted Mass [mm w.e.] | ! SIrnof : Surficial Water + Run OFF Change [mm w.e.] | ! SIvAcr : Sea-Ice vertical Acretion [mm w.e.] | ! Watsvd : SVAT Water Forcing [mm] | ! Watbal : SVAT Water Balance [W/m2] | ! | ! vk2 : Square of Von Karman Constant [-] | ! sqrCm0 : Factor of Neutral Drag Coeffic.Momentum [s/m] | ! sqrCh0 : Factor of Neutral Drag Coeffic.Heat [s/m] | ! EmiSol : Soil Emissivity [-] | ! EmiSno : Snow Emissivity [-] | ! EmiWat : Water Emissivity [-] | ! Z0mLnd : Land Roughness Length [m] | ! sqrrZ0 : u*t/u* | ! f_eff : Marticorena & B. 1995 JGR (20) | ! A_Fact : Fundamental * Roughness | ! Z0mBSn : BSnow Roughness Length [m] | ! Z0mBS0 : Mimimum BSnow Roughness Length (blown* ) [m] | ! Z0m_Sn : Snow Roughness Length (surface) [m] | ! Z0m_S0 : Mimimum Snow Roughness Length [m] | ! Z0m_S1 : Maximum Snow Roughness Length [m] | ! Z0_GIM : Minimum GIMEX Roughness Length [m] | ! Z0_ICE : Sea Ice ISW Roughness Length [m] | ! | ! | !--------------------------------------------------------------------------+ ! Global Variables ! ================ USE VARphy USE VAR_SV USE VARdSV USE VAR0SV USE VARxSV USE VARySV USE VARtSV USE surface_data, only: iflag_tsurf_inlandsis IMPLICIT NONE logical SnoMod logical BloMod integer jjtime ! Internal Variables ! ================== ! Non Local ! --------- real TBr_sv(klonv) ! Brightness Temperature real IRdwsv(klonv) ! DOWNward IR Flux real IRupsv(klonv) ! UPward IR Flux real d_Bufs,Bufs_N ! Buffer Snow Layer Increment real Buf_ro,Bros_N ! Buffer Snow Layer Density real BufPro ! Buffer Snow Layer Density real Buf_G1,BG1__N ! Buffer Snow Layer Dendr/Sphe[-] real Buf_G2,BG2__N ! Buffer Snow Layer Spher/Size[-] real Bdzssv(klonv) ! Buffer Snow Layer Thickness real z_snsv(klonv) ! Snow-Ice, current Thickness ! Local ! ----- integer iwr integer ikl ,isn ,isl ,ist ! integer ist__s,ist__w ! Soil/Water Body Identifier integer growth ! Seasonal Mask integer LISmsk ! Land+Ice / Open Sea Mask integer LSnMsk ! Snow-Ice / No Snow-Ice Mask integer IceMsk,IcIndx(klonv) ! Ice / No Ice Mask integer SnoMsk ! Snow / No Snow Mask real roSMin,roSMax,roSn_1,roSn_2,roSn_3 ! Fallen Snow Density (PAHAUT) real Dendr1,Dendr2,Dendr3 ! Fallen Snow Dendric.(GIRAUD) real Spher1,Spher2,Spher3,Spher4 ! Fallen Snow Spheric.(GIRAUD) real Polair ! Polar Snow Switch real PorSno,Por_BS,Salt_f,PorRef ! c #sw real PorVol,rWater ! c #sw real rusNEW,rdzNEW,etaNEW ! real ro_new ! real TaPole ! Maximum Polar Temperature real T__Min ! Minimum realistic Temperature real EmiSol ! Emissivity of Soil real EmiSno ! Emissivity of Snow real EmiWat ! Emissivity of a Water Area real vk2 ! Square of Von Karman Constant real u2star !(u*)**2 real Z0mLnd ! Land Roughness Length c #ZN real sqrrZ0 ! u*t/u* real f_eff ! Marticorena & B. 1995 JGR (20) real A_Fact ! Fundamental * Roughness real Z0m_nu ! Smooth R Snow Roughness Length real Z0mBSn ! BSnow Roughness Length real Z0mBS0 ! Mimimum BSnow Roughness Length real Z0m_S0 ! Mimimum Snow Roughness Length real Z0m_S1 ! Maximum Snow Roughness Length c #SZ real Z0Sa_N ! Regime Snow Roughness Length c #SZ real Z0SaSi ! 1.IF Rgm Snow Roughness Length c #GL real Z0_GIM ! Mimimum GIMEX Roughness Length real Z0_ICE ! Ice ISW Roughness Length real Z0m_Sn,Z0m_90 ! Snow Surface Roughness Length real SnoWat ! Snow Layer Switch c #RN real rstar,alors ! c #RN real rstar0,rstar1,rstar2 ! real SameOK ! 1. => Same Type of Grains real G1same ! Averaged G1, same Grains real G2same ! Averaged G2, same Grains real typ__1 ! 1. => Lay1 Type: Dendritic real zroNEW ! dz X ro, if fresh Snow real G1_NEW ! G1, if fresh Snow real G2_NEW ! G2, if fresh Snow real zroOLD ! dz X ro, if old Snow real G1_OLD ! G1, if old Snow real G2_OLD ! G2, if old Snow real SizNEW ! Size, if fresh Snow real SphNEW ! Spheric.,if fresh Snow real SizOLD ! Size, if old Snow real SphOLD ! Spheric.,if old Snow real Siz_av ! Averaged Grain Size real Sph_av ! Averaged Grain Spher. real Den_av ! Averaged Grain Dendr. real DendOK ! 1. => Average is Dendr. real G1diff ! Averaged G1, diff. Grains real G2diff ! Averaged G2, diff. Grains real G1 ! Averaged G1 real G2 ! Averaged G2 real param ! Polynomial fit z0=f(T) real Z0_obs ! Fit Z0_obs=f(T) (m) real tamin ! min T of linear fit (K) real tamax ! max T of linear fit (K) real coefa,coefb,coefc,coefd ! Coefs for z0=f(T) real ta1,ta2,ta3 ! Air temperature thresholds real z01,z02,z03 ! z0 thresholds real tt_c,vv_c ! Critical param. real tt_tmp,vv_tmp,vv_virt ! Temporary variables logical density_kotlyakov ! .true. if Kotlyakov 1961 real e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations real zm1, zm2, coefslope ! variables for surface temperature extrapolation ! Internal DATA ! ============= data T__Min / 200.00/ ! Minimum realistic Temperature data TaPole / 263.15/ ! Maximum Polar Temperature data roSMin / 30. / ! Minimum Snow Density data roSMax / 400. / ! Max Fresh Snow Density data tt_c / -2.0 / ! Critical Temp. (degC) data vv_c / 14.3 / ! Critical Wind speed (m/s) data roSn_1 / 109. / ! Fall.Sno.Density, Indep. Param. data roSn_2 / 6. / ! Fall.Sno.Density, Temper.Param. data roSn_3 / 26. / ! Fall.Sno.Density, Wind Param. data Dendr1 / 17.12/ ! Fall.Sno.Dendric.,Wind 1/Param. data Dendr2 / 128. / ! Fall.Sno.Dendric.,Wind 2/Param. data Dendr3 / -20. / ! Fall.Sno.Dendric.,Indep. Param. data Spher1 / 7.87/ ! Fall.Sno.Spheric.,Wind 1/Param. data Spher2 / 38. / ! Fall.Sno.Spheric.,Wind 2/Param. data Spher3 / 50. / ! Fall.Sno.Spheric.,Wind 3/Param. data Spher4 / 90. / ! Fall.Sno.Spheric.,Indep. Param. data EmiSol / 0.99999999/ ! 0.94Emissivity of Soil data EmiWat / 0.99999999/ ! Emissivity of a Water Area data EmiSno / 0.99999999/ ! Emissivity of Snow ! DATA Emissivities ! Pielke, 1984, pp. 383,409 data Z0mBS0 / 0.5e-6/ ! MINimum Snow Roughness Length ! for Momentum if Blowing Snow ! Gallee et al. 2001 BLM 99 (19) data Z0m_S0/ 0.00005/ ! MINimum Snow Roughness Length ! MegaDunes included data Z0m_S1/ 0.030 / ! MAXimum Snow Roughness Length ! (Sastrugis) c #GL data Z0_GIM/ 0.0013/ ! Ice Min Z0 = 0.0013 m (Broeke) ! ! Old Ice Z0 = 0.0500 m (Bruce) ! ! 0.0500 m (Smeets) ! ! 0.1200 m (Broeke) data Z0_ICE/ 0.0010/ ! Sea-Ice Z0 = 0.0010 m (Andreas) ! ! (Ice Station Weddel -- ISW) vk2 = vonKrm * vonKrm ! Square of Von Karman Constant ! BEGIN.main. ! =========================== ! "Soil" Humidity of Water Bodies ! =============================== DO ikl=1,knonv ist = isotSV(ikl) ! Soil Type ist__s = min(ist, 1) ! 1 => Soil ist__w = 1 - ist__s ! 1 => Water Body DO isl=-nsol,0 eta_SV(ikl,isl) = eta_SV(ikl,isl) * ist__s ! Soil . + etadSV(ist) * ist__w ! Water Body END DO ! Vertical Discretization Factor ! ============================== LSdzsv(ikl) = ist__s ! Soil . + OcndSV * ist__w ! Water Body END DO ! Blowing Particles Threshold Friction velocity ! ============================================= c #AE usthSV(ikl) = 1.0e+2 ! END DO !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! Contribution of Snow to the Surface Snow Pack ! ============================================= IF (SnoMod) THEN C +--Blowing Snow C + ------------ IF (BloMod) then if (klonv.eq.1) then if(isnoSV(1).ge.2 .and. . TsisSV(1,max(1,isnoSV(1)))<273. .and. . ro__SV(1,max(1,isnoSV(1)))<500. .and. . eta_SV(1,max(1,isnoSV(1))) Dendritic c #BS zroNEW = typ__1 *(1.0-dsnbSV(ikl)) ! fract.Dendr.Lay. c #BS. + (1.-typ__1) * dsnbSV(ikl) ! c #BS G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay. c #BS. + (1.-typ__1) *G1_dSV ! c #BS G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay. c #BS. + (1.-typ__1) *ADSdSV ! c #BS zroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) ! fract.Spher.Lay. c #BS. + typ__1 * dsnbSV(ikl) ! c #BS G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay. c #BS. + typ__1 *G1_dSV ! c #BS G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay. c #BS. + typ__1 *ADSdSV ! c #BS SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay. c #BS. +(1.+G1_NEW /G1_dSV) ! c #BS. *(G2_NEW *DScdSV/G1_dSV ! c #BS. +(1.-G2_NEW /G1_dSV)*DFcdSV) ! c #BS SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay. c #BS SizOLD = G2_OLD ! Size Spher.Lay. c #BS SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay. c #BS Siz_av = (zroNEW*SizNEW+zroOLD*SizOLD) ! Averaged Size c #BS Sph_av = min( zroNEW*SphNEW+zroOLD*SphOLD ! c #BS. , unun) ! Averaged Sphericity c #BS Den_av = min((Siz_av -( Sph_av *DScdSV ! c #BS. +(1.-Sph_av)*DFcdSV)) ! c #BS. / (DDcdSV -( Sph_av *DScdSV ! c #BS. +(1.-Sph_av)*DFcdSV)) ! c #BS. , unun) ! c #BS DendOK = max(zero, ! c #BS. sign(unun, Sph_av *DScdSV ! Small Grains c #BS. +(1.-Sph_av)*DFcdSV ! Faceted Grains c #BS. - Siz_av )) ! C +... REMARQUE: le type moyen (dendritique ou non) depend C + ^^^^^^^^ de la comparaison avec le diametre optique C + d'une neige recente de dendricite nulle C +... REMARK: the mean type (dendritic or not) depends C + ^^^^^^ on the comparaison with the optical diameter C + of a recent snow having zero dendricity c #BS G1diff =( -DendOK *Den_av c #BS. +(1.-DendOK)*Sph_av) *G1_dSV c #BS G2diff = DendOK *Sph_av *G1_dSV c #BS. +(1.-DendOK)*Siz_av c #BS G1 = SameOK *G1same c #BS. +(1.-SameOK)*G1diff c #BS G2 = SameOK *G2same c #BS. +(1.-SameOK)*G2diff ! S.1. Meme Type de Neige / same Grain Type ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ SameOK = max(zero, . sign(unun, Buf_G1 *BG1sSV(ikl) . - eps_21 )) G1same = (d_Bufs*Buf_G1+BufsSV(ikl)*BG1sSV(ikl)) . /max(epsi,Bufs_N) G2same = (d_Bufs*Buf_G2+BufsSV(ikl)*BG2sSV(ikl)) . /max(epsi,Bufs_N) ! S.2. Types differents / differents Types ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic zroNEW =( typ__1 *d_Bufs ! fract.Dendr.Lay. . + (1.-typ__1) *BufsSV(ikl)) ! . /max(epsi,Bufs_N) ! G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay. . + (1.-typ__1) *BG1sSV(ikl) ! G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay. . + (1.-typ__1) *BG2sSV(ikl) ! zroOLD =((1.-typ__1) *d_Bufs ! fract.Spher.Lay. . + typ__1 *BufsSV(ikl)) ! . /max(epsi,Bufs_N) ! G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay. . + typ__1 *BG1sSV(ikl) ! G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay. . + typ__1 *BG2sSV(ikl) ! SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay. . +(1.+G1_NEW /G1_dSV) ! . *(G2_NEW *DScdSV/G1_dSV ! . +(1.-G2_NEW /G1_dSV)*DFcdSV) ! SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay. SizOLD = G2_OLD ! Size Spher.Lay. SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay. Siz_av = ( zroNEW *SizNEW+zroOLD*SizOLD) ! Averaged Size Sph_av = min( zroNEW *SphNEW+zroOLD*SphOLD ! . , unun ) ! Averaged Sphericity Den_av = min((Siz_av - ( Sph_av *DScdSV ! . +(1.-Sph_av)*DFcdSV)) ! . / (DDcdSV - ( Sph_av *DScdSV ! . +(1.-Sph_av)*DFcdSV)) ! . , unun )! DendOK = max(zero, ! . sign(unun, Sph_av *DScdSV ! Small Grains . +(1.-Sph_av)*DFcdSV ! Faceted Grains . - Siz_av )) ! C +... REMARQUE: le type moyen (dendritique ou non) depend C + ^^^^^^^^ de la comparaison avec le diametre optique C + d'une neige recente de dendricite nulle C +... REMARK: the mean type (dendritic or not) depends C + ^^^^^^ on the comparaison with the optical diameter C + of a recent snow having zero dendricity G1diff =( -DendOK *Den_av . +(1.-DendOK)*Sph_av) *G1_dSV G2diff = DendOK *Sph_av *G1_dSV . +(1.-DendOK)*Siz_av G1 = SameOK *G1same . +(1.-SameOK)*G1diff G2 = SameOK *G2same . +(1.-SameOK)*G2diff BG1sSV(ikl) = G1 ! . * Bufs_N/max(epsi,Bufs_N) ! BG2sSV(ikl) = G2 ! . * Bufs_N/max(epsi,Bufs_N) ! C +--Update of Buffer Layer Content & Decision about creating a new snow layer C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ BufsSV(ikl) = Bufs_N ! [mm w.e.] NLaysv(ikl) = min(unun, ! . max(zero, ! Allows to create . sign(unun,BufsSV(ikl) ! a new snow Layer . -SMndSV )) ! if Buffer > SMndSV . *max(zero, ! Except if * Erosion . sign(unun,0.50 ! dominates . -dsnbSV(ikl))) ! . +max(zero, ! Allows to create . sign(unun,BufsSV(ikl) ! a new snow Layer . -SMndSV*3.00))) ! is Buffer > SMndSV*3 Bdzssv(ikl) = 1.e-3*BufsSV(ikl)*ro_Wat ! [mm w.e.] -> [m w.e.] . /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m] END DO ! Snow Pack Discretization ! ======================== ! ********** call SISVAT_zSn ! ********** ! ********** ! #ve call SISVAT_wEq('_zSn ',0) ! ********** ! Add a new Snow Layer ! ==================== DO ikl=1,knonv isnoSV(ikl) = isnoSV(ikl) +NLaysv(ikl) isn = isnoSV(ikl) dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl)) . + Bdzssv(ikl) * NLaysv(ikl) TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl)) . + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl) ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl)) . + Brossv(ikl) * NLaysv(ikl) eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl)) ! + 0. agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl)) ! + 0. G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl)) . + BG1ssv(ikl) * NLaysv(ikl) G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl)) . + BG2ssv(ikl) * NLaysv(ikl) istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl)) . + max(zer0,sign(un_1,TaT_SV(ikl) . -Tf_Sno-eps_21)) * istdSV(2) . * NLaysv(ikl) BufsSV(ikl) = BufsSV(ikl) * (1-NLaysv(ikl)) NLaysv(ikl) = 0 END DO ! Snow Pack Thickness ! ------------------- DO ikl=1,knonv z_snsv(ikl) = 0.0 END DO DO isn=1,nsno DO ikl=1,knonv z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn) zzsnsv(ikl,isn) = z_snsv(ikl) END DO END DO END IF ! Soil Albedo: Soil Humidity Correction ! ========================================== ! REFERENCE: McCumber and Pielke (1981), Pielke (1984) ! ^^^^^^^^^ DO ikl=1,knonv albssv(ikl) = . alb0SV(ikl) *(1.0-min(half,eta_SV( ikl,0) . /etadSV(isotSV(ikl)))) ! REMARK: Albedo of Water Surfaces (isotSV=0): ! ^^^^^^ alb0SV := 2 X effective value, while ! eta_SV := etadSV END DO ! Snow Pack Optical Properties ! ============================ IF (SnoMod) THEN ! ****** call SnOptP(jjtime) ! ****** ELSE DO ikl=1,knonv sEX_sv(ikl,1) = 1.0 sEX_sv(ikl,0) = 0.0 albisv(ikl) = albssv(ikl) END DO END IF ! Soil optical properties ! ============================= !Etienne: as in inlandis we do not call vgopt, we need to define !the albedo alb_SV and to calculate the !absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv DO ikl=1,klonv e_pRad = 2.5 * coszSV(ikl) ! exponential argument, ! V/nIR radiation partitioning, ! DR97, 2, eqn (2.53) & (2.54) e1pRad = 1.-exp(-e_pRad) ! exponential, V/nIR Rad. Part. exdRad= 1. ! Visible Part of the Solar Radiation Spectrum (V, 0.4--0.7mi.m) A_Rad0 = 0.25 + 0.697 * e1pRad ! Absorbed Vis. Radiation absg_V = (1.-albisv(ikl))*(A_Rad0*exdRad) ! ! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m) A_Rad0 = 0.80 + 0.185 * e1pRad ! Absorbed nIR. Radiation absgnI = (1.-albisv(ikl))*(A_Rad0*exdRad) ! SoSosv(ikl) = (absg_V+absgnI)*0.5d0 alb_SV(ikl) = albisv(ikl) END DO ! ********** ! #ve call SISVAT_wEq('SnOptP',0) ! ********** ! Surface Emissivity (Etienne: simplified calculation for landice) ! ============================================================= ! DO ikl=1,knonv LSnMsk = min( 1,isnoSV(ikl)) Eso_sv(ikl)= EmiSol*(1-LSnMsk)+EmiSno*LSnMsk ! Sol+Sno Emissivity emi_SV(ikl)= EmiSol*(1-LSnMsk) + EmiSno*LSnMsk END DO ! Upward IR (INPUT, from previous time step) ! =================================================================== DO ikl=1,knonv ! #e1 Enrsvd(ikl) = - IRs_SV(ikl) IRupsv(ikl) = IRs_SV(ikl) END DO ! Turbulence ! ========== ! Latent Heat of Vaporization/Sublimation ! --------------------------------------- DO ikl=1,knonv SnoWat = min(isnoSV(ikl),0) Lx_H2O(ikl) = . (1.-SnoWat) * LhvH2O . + SnoWat *(LhsH2O * (1.-eta_SV(ikl,isnoSV(ikl))) . +LhvH2O * eta_SV(ikl,isnoSV(ikl)) ) END DO ! Aerodynamic Resistance ! ^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6)) rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6)) END DO ! Soil Energy Balance ! ===================== if (iflag_tsurf_inlandsis .eq. 0) then call SISVAT_TSo else call SISVAT_TS2 end if ! ********** ! #ve call SISVAT_wEq('_TSo ',0) ! ********** ! Soil Water Potential ! ------------------------ DO isl=-nsol,0 DO ikl=1,knonv ist = isotSV(ikl) ! Soil Type psi_sv(ikl,isl) = psidSV(ist) ! DR97, Eqn.(3.34) . *(etadSV(ist) /max(eps6,eta_SV(ikl,isl))) ! . **bCHdSV(ist) ! ! Soil Hydraulic Conductivity ! --------------------------- Khydsv(ikl,isl) = s2__SV(ist) ! DR97, Eqn.(3.35) . *(eta_SV(ikl,isl)**(2.*bCHdSV(ist)+3.)) ! END DO END DO ! Melting / Refreezing in the Snow Pack ! ===================================== IF (SnoMod) THEN ! ********** call SISVAT_qSn ! ********** ! ********** ! #ve call SISVAT_wEq('_qSn ',0) ! ********** ! Snow Pack Thickness ! ------------------- DO ikl=1,knonv z_snsv(ikl) = 0.0 END DO DO isn=1,nsno DO ikl=1,knonv z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn) zzsnsv(ikl,isn) = z_snsv(ikl) END DO END DO ! Energy in Excess is added to the first Soil Layer ! ------------------------------------------------- DO ikl=1,knonv z_snsv(ikl) = max(zer0, . sign(un_1,eps6-z_snsv(ikl))) TsisSV(ikl,0) = TsisSV(ikl,0) + EExcsv(ikl) . /(rocsSV(isotSV(ikl)) . +rcwdSV*eta_SV(ikl,0)) EExcsv(ikl) = 0. END DO END IF ! Soil Water Balance ! ===================== ! ********** call SISVAT_qSo ! #m0. (Wats_0,Wats_1,Wats_d) ! ********** ! Surface Fluxes ! ===================== DO ikl=1,knonv IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR ! IRdwsv(ikl)=tau_sv(ikl) *IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR ! . +(1.0-tau_sv(ikl))*IRd_SV(ikl)*Evg_sv(ikl) ! ! Etienne, remove vegetation component IRupsv(ikl) = IRupsv(ikl) ! Upward IR IRu_SV(ikl) = -IRupsv(ikl) ! Upward IR . +IRd_SV(ikl) ! (effective) . -IRdwsv(ikl) ! (positive) TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo)) ! Brightness ! ! Temperature uts_SV(ikl) = (HSv_sv(ikl) +HSs_sv(ikl)) ! u*T* . /(rhT_SV(ikl) *cp) ! uqs_SV(ikl) = (HLv_sv(ikl) +HLs_sv(ikl)) ! u*q* . /(rhT_SV(ikl) *LhvH2O) ! LMO_SV(ikl) = TaT_SV(ikl)*(us__SV(ikl)**3) . /gravit/uts_SV(ikl)/vonKrm ! MO length ! Surface Temperature ! ^^^^^^^^^^^^^^^^^^^^ ! Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl)) ! Etienne: extrapolation from the two uppermost levels: if (isnoSV(ikl) >=2) then zm1=-dzsnSV(ikl,isnoSV(ikl))/2. zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.) else if (isnoSV(ikl) .EQ. 1) then zm1=-dzsnSV(ikl,isnoSV(ikl))/2. zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.) else zm1=-dz_dSV(0)/2. zm2=-(dz_dSV(0)+dz_dSV(-1)/2.) end if coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1)) . /(zm1-zm2) Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1) END DO ! Snow Pack Properties (sphericity, dendricity, size) ! =================================================== IF (SnoMod) THEN ! ********** call SISVAT_GSn ! ********** ! ********** ! #ve call SISVAT_wEq('_GSn ',0) ! ********** END IF ! Roughness Length for next time step !==================================== ! Note that in INLANDSIS, we treat only ice covered surfaces so calculation ! of z0 is much simpler (no subgrid fraction of ocean or land) ! old calculations are commented below C +--Roughness Length for Momentum C + ----------------------------- C +--Land+Sea-Ice / Ice-free Sea Mask C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,klonv IcIndx(ikl) = 0 ENDDO DO isn=1,nsno DO ikl=1,klonv IcIndx(ikl) = max(IcIndx(ikl), . isn*max(0, . sign(1, . int(ro__SV(ikl,isn)-900.)))) ENDDO ENDDO DO ikl=1,klonv LISmsk = 1. ! in inlandsis, land only IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) ) SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0) Z0mLnd =max( Z0_ICE , 5.e-5 ) ! Min set := Z0 on * C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03)) C +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ u2star = us__SV(ikl) *us__SV(ikl) Z0mBSn = u2star *0.536e-3 - 61.8e-6 Z0mBSn = max(Z0mBS0 ,Z0mBSn) C +--Z0 Smooth + Saltat. Regime C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0enSV(ikl) = Z0m_nu . + Z0mBSn C +--Rough Snow Surface Roughness Length (Typical Value) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c #tz Z0m_Sn = 0.250e-3 ! Andreas 1995, CRREL Report 95-16, fig.1&p.2 ! z0r~(10-d)*exp(-vonkar/sqrt(1.5e-03))-5.e-5 Z0m_Sn = 2.000e-3 ! Calibration of MAR c #TZ Z0m_Sn = 1.000e-3 ! Exemple Tuning in RACMO c #TZ Z0m_Sn = 0.500e-3 ! Exemple Tuning in MAR C +--Rough Snow Surface Roughness Length (Variable Sastrugi Height) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A_Fact = 1.0000 ! Andreas et al., 2004, p.4 ! ams.confex.com/ams/pdfpapers/68601.pdf ! Parameterization of z0 dependance on Temperature (C. Amory, 2017) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013 coefa = 0.1658 !0.1862 !Ant coefb = -50.3869 !-55.7718 !Ant ta1 = 253.15 !255. Ant ta2 = 273.15 ta3 = 273.15+3 z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm z02 = exp(coefa*ta2 + coefb) !~6 !~7 mm z03 = z01 coefc = log(z03/z02)/(ta3-ta2) coefd = log(z03)-coefc*ta3 if (TaT_SV(ikl) .lt. ta1) then Z0_obs = z01 else if (TaT_SV(ikl).ge.ta1 .and. TaT_SV(ikl).lt.ta2) then Z0_obs = exp(coefa*TaT_SV(ikl) + coefb) else if (TaT_SV(ikl).ge.ta2 .and. TaT_SV(ikl).lt.ta3) then ! if st > 0, melting induce smooth surface Z0_obs = exp(coefc*TaT_SV(ikl) + coefd) else Z0_obs = z03 endif ! pour le moment, on choisit une valeur fixe Z0_obs = 1.000e-3 cCA Snow roughness lenght deduced from observations cCA (parametrization if no Blowing Snow) cCA ----------------------------------- C. Agosta 09-2016 ----- cCA Substract Z0enSV(ikl) because re-added later in Z0mnSV(ikl) Z0m_Sn = Z0_obs - Z0enSV(ikl) cCA ----------------------------------------------------------- param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING c #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING c #SZ. * max(zero,sign(unun,TfSnow-eps9 c #SZ. -TsisSV(ikl , isnoSV(ikl)))) !!#SZ Z0SaSi = max(zero,sign(unun,Z0Sa_N ))! 1 if erosion c #SZ Z0SaSi = max(zero,sign(unun,zero -eps9 -uss_SV(ikl)))! c #SZ Z0Sa_N = max(zero, Z0Sa_N) c #SZ Z0SaSV(ikl) = c #SZ. max(Z0SaSV(ikl) ,Z0SaSV(ikl) c #SZ. + Z0SaSi*(Z0Sa_N-Z0SaSV(ikl))*exp(-dt__SV/43200.)) c #SZ. - min(dz0_SV(ikl) , Z0SaSV(ikl)) c #SZ A_Fact = Z0SaSV(ikl) * 5.0/0.15 ! A=5 if h~10cm C +... CAUTION: The influence of the sastrugi direction is not yet included c #SZ Z0m_Sn = Z0SaSV(ikl) ! c #SZ. - Z0m_nu ! C +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46) p.222) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ c #ZN sqrrZ0 = usthSV(ikl)/max( us__SV(ikl),0.001) c #ZN sqrrZ0 = min( sqrrZ0 ,0.999) c #ZN Z0mBSn = 0.55 *0.55 *exp(-sqrrZ0 *sqrrZ0) c #ZN. *us__SV(ikl)* us__SV(ikl)*grvinv*0.5 C +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46) p.222) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ c #ZN Z0enSV(ikl) = (Z0m_nu ** sqrrZ0 ) c #ZN. * (Z0mBSn **(1.-sqrrZ0)) c #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu) C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf) c #ZA Z0m_nu = 0.135*akmol / max(us__SV(ikl) , epsi) C +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf) c #ZA Z0mBSn = 0.035*u2star *grvinv C +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004 ! ( used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf) ! ^^^^^^^^^^^^^^^^^^^^^^^^^^ c #ZA Z0enSV(ikl) = Z0m_nu c #ZA. + Z0mBSn C +--Z0 Rough Regime over Snow (Andreas etAl., 2004 C + (.NOT. used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf) ! ^^^^^^^^^^^^^^^^^^^^^^^^^^ !!#ZA u2star = (us__SV(ikl) -0.1800) / 0.1 !!#ZA Z0m_Sn =A_Fact*Z0mBSn *exp(-u2star*u2star) c #ZA Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.) c #ZA. *exp(-0.4/sqrt(.00275+.00001*max(0.,VVs_SV(ikl)-5.))) c #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45. c #ZA. - DDs_SV(ikl)*DDs_SV(ikl)* Z0m_90 /(90.*90.) C +--Z0 (Erosion) over Snow (instantaneous or time average) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0e_SV(ikl) = Z0enSV(ikl) Z0e_SV(ikl) = Z0emSV(ikl) C +--Momentum Roughness Length C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Contribution of Z0mnSV(ikl) = Z0mLnd ! land Form . + (Z0m_Sn ! Sastrugi Form . + Z0enSV(ikl)) *SnoMsk ! Snow Erosion C +--GIS Roughness Length C + ^^^^^^^^^^^^^^^^^^^^^ c #GL Z0mnSV(ikl) = c #GL. (1-LSmask(ikl)) * Z0mnSV(ikl) c #GL. + LSmask(ikl) * max(Z0mnSV(ikl),max(Z0_GIM, c #GL. Z0_GIM+ c #GL. (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.) ! c #GL. /(920.00 -600.))) ! C +--Mom. Roughness Length, Instantaneous OR Box Moving Average in Time C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant. ! Z0m_SV(ikl) = Z0mmSV(ikl) ! Z0mnSV Average C +--Corrected Threshold Friction Velocity before Erosion ! Marticorena and C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Bergametti 1995 ! not used anymore since Marticorena and Bergametti disabled !CK 18/07/2018 cc #BS Z0e_SV(ikl) = min(Z0m_SV(ikl),Z0e_SV(ikl)) ! cc #MB f_eff= log(0.35*(0.1 /Z0e_SV(ikl))**0.8) ! JGR 100 cc #MB f_eff=1.-(log( Z0m_SV(ikl)/Z0e_SV(ikl) ))! (20) p. 16420 cc #MB. /(max( f_eff ,epsi ))! p.16426 2nd ? cc #MB f_eff= max( f_eff ,epsi )! CONTROL cc #MB f_eff=1.0 -(1.0 - f_eff) /5.00 ! TUNING cc #MB f_eff= min( f_eff ,1.00 )! cc #MB usthSV(ikl) = usthSV(ikl)/f_eff ! C +--Roughness Length for Scalars C + ---------------------------- Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4 c #SH Z0hnSV(ikl) = Z0mnSV(ikl)/100.0 C + Z0h = Z0m /100.0 over the Sahel C + (Taylor & Clark, QJRMS 127,p864) c #RN rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol c #RN rstar = max(epsi,min(rstar,thous)) c #RN alors = log(rstar) c #RN rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar)) c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar))) c #RN. *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar)) c #RN. + 0.317e0 c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar)))) c #RN rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar)) c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar))) c #RN. *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar)) c #RN. - 0.565 c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar)))) c #RN rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar)) c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar))) c #RN. *(0. * max(zero,sign(unun,2.500e0 - rstar)) c #RN. - 0.183 c #RN. *(unun - max(zero,sign(unun,2.500e0 - rstar)))) cXF #RN does not work over bare ice cXF MAR is then too warm and not enough melt c #RN if(ro__SV(ikl,isnoSV(ikl))>50 c #RN. .and.ro__SV(ikl,isnoSV(ikl))