subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut) 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: is_ok_z0h_rn, . is_ok_density_kotlyakov, . prescribed_z0m_snow, . iflag_z0m_snow, . iflag_tsurf_inlandsis, . iflag_temp_inlandsis, . discret_xf, buf_sph_pol,buf_siz_pol IMPLICIT NONE logical SnoMod logical BloMod logical debut 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,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 real rstar,alors ! 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 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 real e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations real zm1, zm2, coefslope ! variables for surface temperature extrapolation ! for Aeolian erosion and blowing snow integer nit ,iit real Fac ! Correc. factor for drift ratio real dusuth,signus real sss__F,sss__N real sss__K,sss__G real us_127,us_227,us_327,us_427,us_527 real VVa_OK, usuth0 real ssstar real SblPom real rCd10n ! Square root of drag coefficient real DendOK ! Dendricity Switch real SaltOK ! Saltation Switch real MeltOK ! Saltation Switch (Melting Snow) real SnowOK ! Pack Top Switch real SaltM1,SaltM2,SaltMo,SaltMx ! Saltation Parameters real ShearX, ShearS ! Arg. Max Shear Stress real Por_BS ! Snow Porosity real Salt_us ! New thresh.friction velocity u*t real Fac_Mo,ArguSi,FacRho ! Numerical factors for u*t real SaltSI(klonv,0:nsno) ! Snow Drift Index ! real MIN_Mo ! Minimum Mobility Fresh Fallen * character*3 qsalt_param ! Switch for saltation flux param. character*3 usth_param ! Switch for u*t param ! Internal DATA ! ============= data T__Min / 200.00/ ! Minimum realistic Temperature data TaPole / 268.15/ ! Maximum Polar Temperature (value from C. Agosta) data roSMin / 300. / ! 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) ! for aerolian erosion data SblPom/ 1.27/ ! Lower Boundary Height Parameter C + ! for Suspension C + ! Pommeroy, Gray and Landine 1993, C + ! J. Hydrology, 144(8) p.169 data nit / 5 / ! us(is0,uth) recursivity: Nb Iterations cc#AE data qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p data qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray cc#AE data usth_param/"lis"/ ! u*t from Liston et al. 2007 data usth_param/"gal"/ ! u*t from Gallee et al. 2001 data SaltMx/-5.83e-2/ 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 IF (SnoMod) THEN C +--Aeolian erosion and Blowing Snow C +================================== DO ikl=1,knonv usthSV(ikl) = 1.0e+2 END DO IF (BloMod) THEN if (klonv==1) then if(isnoSV(1)>=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)))300.) then Por_BS = 1.000 - ro__SV(ikl,isn) /ro_Ice else Por_BS = 1.000 - 300. /ro_Ice endif ShearX = Por_BS/max(epsi,1.-Por_BS) Fac_Mo = exp(-ShearX+ShearS) C + Gallee et al., 2001 eq 5, p5 if (usth_param == "gal") then Salt_us = (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085 Salt_us = Salt_us * Fac_Mo C +... Salt_us : Extension of Guyomarc'h & Merindol 1998 with C +... de Montmollin (1978). Gallee et al. 2001 endif if (usth_param == "lis") then !Liston et al. 2007 if(ro__SV(ikl,isn)>300.) then Salt_us = 0.005*exp(0.013*ro__SV(ikl,isn)) else Salt_us = 0.01*exp(0.003*ro__SV(ikl,isn)) endif endif SnowOK = 1 -min(1,iabs(isn-isnoSV(ikl))) !Switch new vs old snow usthSV(ikl) = SnowOK * (Salt_us) . + (1.-SnowOK)* usthSV(ikl) END DO ! Feeback between blowing snow turbulent Scale u* (commented here ! since ustar is an input variable (not in/out) of inlandsis) ! ----------------------------------------------------------------- ! VVa_OK = max(0.000001, VVaSBL(ikl)) ! sss__N = vonkar * VVa_OK ! sss__F = (sqrCm0(ikl) - psim_z + psim_0) ! usuth0 = sss__N /sss__F ! u* if NO Blow. Snow ! sss__G = 0.27417 * gravit ! ! ______________ _____ ! ! Newton-Raphson (! Iteration, BEGIN) ! ! ~~~~~~~~~~~~~~ ~~~~~ ! DO iit=1,nit ! sss__K = gravit * r_Turb * A_Turb *za__SV(ikl) ! . *rCDmSV(ikl)*rCDmSV(ikl) ! . /(1.+0.608*QaT_SV(ikl)-qsnoSV(ikl)) ! us_127 = exp( SblPom *log(us__SV(ikl))) ! us_227 = us_127 * us__SV(ikl) ! us_327 = us_227 * us__SV(ikl) ! us_427 = us_327 * us__SV(ikl) ! us_527 = us_427 * us__SV(ikl) ! us__SV(ikl) = us__SV(ikl) ! . - ( us_527 *sss__F /sss__N ! . - us_427 ! . - us_227 *qsnoSV(ikl)*sss__K ! . + (us__SV(ikl)*us__SV(ikl)-usthSV(ikl)*usthSV(ikl))/sss__G) ! . /( us_427*5.27*sss__F /sss__N ! . - us_327*4.27 ! . - us_127*2.27*qsnoSV(ikl)*sss__K ! . + us__SV(ikl)*2.0 /sss__G) ! us__SV(ikl)= min(us__SV(ikl),usuth0) ! us__SV(ikl)= max(us__SV(ikl),epsi ) ! rCDmSV(ikl)= us__SV(ikl)/VVa_OK ! ! #AE sss__F = vonkar /rCDmSV(ikl) ! ENDDO ! ! ______________ ___ ! ! Newton-Raphson (! Iteration, END ) ! ! ~~~~~~~~~~~~~~ ~~~ ! us_127 = exp( SblPom *log(us__SV(ikl))) ! us_227 = us_127 * us__SV(ikl) ! ! Momentum Turbulent Scale u*: 0-Limit in case of no Blow. Snow ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! dusuth = us__SV(ikl) - usthSV(ikl) ! u* - uth* ! signus = max(sign(unun,dusuth),zero) ! 1 <=> u* - uth* > 0 ! us__SV(ikl) = ! ! . us__SV(ikl) *signus + ! u* (_BS) ! . usuth0 ! u* (nBS) ! . *(1.-signus) ! ! Blowing Snow Turbulent Scale ss* ! --------------------------------------- hSalSV(ikl) = 8.436e-2 * us__SV(ikl)**SblPom if (qsalt_param == "pom") then qSalSV(ikl) = (us__SV(ikl)**2 - usthSV(ikl)**2) *signus . / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25) endif if (qsalt_param == "bin") then qSalSV(ikl) = (us__SV(ikl) * us__SV(ikl) . -usthSV(ikl) * usthSV(ikl))*signus . * 0.535 / (hSalSV(ikl) * gravit) endif qSalSV(ikl) = qSalSV(ikl)/rht_SV(ikl) ! conversion kg/m3 to kg/kg ssstar = rCDmSV(ikl) * (qsnoSV(ikl) - qSalSV(ikl)) . * r_Turb !Bintanja 2000, BLM !r_Turb compensates for an overestim. of the blown snow part. fall velocity uss_SV(ikl) = min(zero , us__SV(ikl) *ssstar) uss_SV(ikl) = max(-0.0001 , uss_SV(ikl)) ENDIF ! BloMod C + ------------------------------------------------------ C +--Buffer Layer C + ----------------------------------------------------- DO ikl=1,knonv c BufsSV(ikl) [mm w.e.] i.e, i.e., [kg/m2] d_Bufs = max(dsn_SV(ikl) *dt__SV,0.) ! dsn_SV(ikl) = 0. ! Bufs_N = BufsSV(ikl) +d_Bufs ! C +--Snow Density C + ^^^^^^^^^^^^ Polair = zero c #NP Polair = max(zero, ! c #NP. sign(unun,TaPole ! c #NP. -TaT_SV(ikl))) ! Polair = max(zero, ! . sign(unun,TaPole ! . -TaT_SV(ikl))) ! Buf_ro = max( rosMin, ! Fallen Snow Density . roSn_1+roSn_2* (TaT_SV(ikl)-TfSnow) ! [kg/m3] . +roSn_3*sqrt( VV__SV(ikl))) ! Pahaut (CEN), Etienne: use wind speed at first model level instead of 10m wind c #NP BufPro = max( rosMin, ! Fallen Snow Density c #NP. 104. *sqrt( max( VV10SV(ikl)-6.0,0.0))) ! Kotlyakov (1961) ! C.Agosta option for snow density, same as for BS i.e. ! is_ok_density_kotlyakov=.false. c #BS density_kotlyakov = .false. !C.Amory BS 2018 C + ... Fallen Snow Density, Adapted for Antarctica if (is_ok_density_kotlyakov) then tt_tmp = TaT_SV(ikl)-TfSnow !vv_tmp = VV10SV(ikl) vv_tmp=VV__SV(ikl) ! Etienne: use wind speed at first model level instead of 10m wind C + ... [ A compromise between C + ... Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ] if (tt_tmp>=-10) then BufPro = max( rosMin, . 104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961) else vv_virt = (tt_c*vv_tmp+vv_c*(tt_tmp+10)) . /(tt_c+tt_tmp+10) BufPro = 104. *sqrt( max( vv_virt-6.0,0.0)) endif else C + ... [ density derived from observations of the first 50cm of C + ... snow - cf. Rajashree Datta - and multiplied by 0.8 ] C + ... C. Agosta, 2016-09 cc #SD BufPro = 149.2 + 6.84*VV10SV(ikl) + 0.48*Tsrfsv(ikl) cc #SD BufPro = 125 + 14*VV10SV(ikl) + 0.6*Tsrfsv(ikl) !MAJ CK and CAm ! BufPro = 200 + 21 * VV10SV(ikl)!CK 29/07/19 BufPro = 200 + 21 * VV__SV(ikl)!Etienne: use wind speed at first model level instead of 10m wind endif Bros_N = (1. - Polair) * Buf_ro ! Temperate Snow . + Polair * BufPro ! Polar Snow Bros_N = max( 20.,max(rosMin, Bros_N)) Bros_N = min(400.,min(rosMax-1,Bros_N)) ! for dz_min in SISVAT_zSn ! Density of deposited blown snow ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (BloMod) then Bros_N = frsno ro_new = ro__SV(ikl,max(1,isnoSV(ikl))) ro_new = max(Bros_N,min(roBdSV,ro_new)) Fac = 1-((ro__SV(ikl,max(1,isnoSV(ikl))) . -roBdSV)/(500.-roBdSV)) Fac = max(0.,min(1.,Fac)) dsnbSV(ikl) = Fac*dsnbSV(ikl) Bros_N = Bros_N * (1.0-dsnbSV(ikl)) . + ro_new * dsnbSV(ikl) endif ! Time averaged Density of deposited blown Snow ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ BrosSV(ikl) =(Bros_N * d_Bufs ! . +BrosSV(ikl)* BufsSV(ikl))! . / max(epsi,Bufs_N) ! C +-- S.Falling Snow Properties (computed as in SISVAT_zAg) C + ^^^^^^^^^^^^^^^^^^^^^^^ Buf_G1 = max(-G1_dSV, ! Temperate Snow . min(Dendr1*VV__SV(ikl)-Dendr2, ! Dendricity . Dendr3 )) ! Buf_G2 = min( Spher4, ! Temperate Snow . max(Spher1*VV__SV(ikl)+Spher2, ! Sphericity . Spher3 )) ! ! EV: now control buf_sph_pol and bug_siz_pol in physiq.def Buf_G1 = (1. - Polair) * Buf_G1 ! Temperate Snow . + Polair * buf_sph_pol ! Polar Snow Buf_G2 = (1. - Polair) * Buf_G2 ! Temperate Snow . + Polair * buf_siz_pol ! Polar Snow G1 = Buf_G1 ! NO Blown Snow G2 = Buf_G2 ! NO Blown Snow IF (BloMod) THEN ! S.1. Meme Type de Neige / same Grain Type ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ SameOK = max(zero, . sign(unun, Buf_G1 *G1_dSV . - eps_21 )) G1same = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV) G2same = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV) ! Blowing Snow Properties: G1_dSV, ADSdSV ! S.2. Types differents / differents Types ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic zroNEW = typ__1 *(1.0-dsnbSV(ikl)) ! fract.Dendr.Lay. . + (1.-typ__1) * dsnbSV(ikl) ! G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay. . + (1.-typ__1) *G1_dSV ! G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay. . + (1.-typ__1) *ADSdSV ! zroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) ! fract.Spher.Lay. . + typ__1 * dsnbSV(ikl) ! G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay. . + typ__1 *G1_dSV ! G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay. . + typ__1 *ADSdSV ! 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 ENDIF ! 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(option XF in MAR) ! ========================================== if (discret_xf.AND.klonv==1) then if(isnoSV(1)>=1.or.NLaysv(1)>=1) then C + ********** call SISVAT_zSn C + ********** endif else C + ********** call SISVAT_zSn C + ********** endif C + ********** ! #ve call SISVAT_wEq('_zSn ',0) C + ********** ! 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 ! SnoMod ! 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 (calculated from drags given by LMDZ) ! Commented because already calculated in surf_inlandsis_mod ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! 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_temp_inlandsis == 0) then call SISVAT_TSo else DO ikl=1,knonv Tsf_SV(ikl)=Tsrfsv(ikl) END DO 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 ! ^^^^^^^^^^^^^^^^^^^^ IF (iflag_tsurf_inlandsis == 0) THEN Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl)) ELSE IF (iflag_tsurf_inlandsis > 0) THEN ! 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) == 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) ELSE !(default) Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl)) END IF END DO ! Snow Pack Properties (sphericity, dendricity, size) ! =================================================== IF (SnoMod) THEN if (discret_xf .AND. klonv==1) then if(isnoSV(1)>=1) then C + ********** call SISVAT_GSn C + ********** endif else C + ********** call SISVAT_GSn C + ********** endif 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 + ----------------------------- ! ETIENNE WARNING: changes have been made wrt original SISVAT C +--Land+Sea-Ice / Ice-free Sea Mask C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv IcIndx(ikl) = 0 ENDDO DO isn=1,nsno DO ikl=1,knonv IcIndx(ikl) = max(IcIndx(ikl), . isn*max(0, . sign(1, . int(ro__SV(ikl,isn)-900.)))) ENDDO ENDDO DO ikl=1,knonv LISmsk = 1. ! in inlandsis, land only IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) ) SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0) 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 ! Calculation of snow roughness length !===================================== IF (iflag_z0m_snow == 0) THEN Z0m_Sn=prescribed_z0m_snow ELSE IF (iflag_z0m_snow == 1) THEN Z0m_Sn=Z0enSV(ikl) ELSE IF (iflag_z0m_snow == 2) THEN 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) < ta1) then Z0_obs = z01 else if (TaT_SV(ikl)>=ta1 .and. TaT_SV(ikl)=ta2 .and. TaT_SV(ikl) 0, melting induce smooth surface Z0_obs = exp(coefc*TaT_SV(ikl) + coefd) else Z0_obs = z03 endif Z0m_Sn=Z0_obs ELSE Z0m_Sn=0.500e-3 ! default=0.500e-3m (tuning of MAR) ENDIF ! 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) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0e_SV(ikl) = Z0enSV(ikl) C +--Momentum Roughness Length (Etienne: changes wrt original SISVAT) C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0mnSV(ikl) = Z0m_nu *(1-SnoMsk) ! Ice z0 . + (Z0m_Sn)*SnoMsk ! Snow Sastrugi Form and 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 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant. C +--Roughness Length for Scalars C + ---------------------------- Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4 IF (is_ok_z0h_rn) THEN rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol rstar = max(epsi,min(rstar,R_1000)) alors = log(rstar) rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar)) . +(1. - max(zero,sign(unun,0.135e0 - rstar))) . *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar)) . + 0.317e0 . *(1. - max(zero,sign(unun,2.500e0 - rstar)))) rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar)) . +(1. - max(zero,sign(unun,0.135e0 - rstar))) . *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar)) . - 0.565 . *(1. - max(zero,sign(unun,2.500e0 - rstar)))) rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar)) . +(1. - max(zero,sign(unun,0.135e0 - rstar))) . *(0. * max(zero,sign(unun,2.500e0 - rstar)) . - 0.183 . *(unun - max(zero,sign(unun,2.500e0 - rstar)))) !XF #RN (is_ok_z0h_rn) does not work well over bare ice !XF MAR is then too warm and not enough melt if(ro__SV(ikl,isnoSV(ikl))>50 . .and.ro__SV(ikl,isnoSV(ikl))