From b47219764ebfedca15fe78f9faf9b6bef9d8d526 Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Wed, 1 May 2024 17:39:19 -0400 Subject: [PATCH 01/10] Fixed a bug where albedo over the ocean is computed incorrectly when ialb = 2 --- gsmphys/radiation_surface.f | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/gsmphys/radiation_surface.f b/gsmphys/radiation_surface.f index 11b6bd57..1656940c 100644 --- a/gsmphys/radiation_surface.f +++ b/gsmphys/radiation_surface.f @@ -664,7 +664,13 @@ subroutine setalb & !> - Calculate snow cover input directly for land model, no !! conversion needed. - fsno0 = f_zero + fsno0 = sncovr(i) + + if (nint(slmsk(i))==0 .and. + & (tsknf(i)>con_tice .or. + & ldisable_radiation_quasi_sea_ice)) then + fsno0 = f_zero + endif if (nint(slmsk(i)) == 2) then asnow = 0.02*snowf(i) @@ -674,7 +680,7 @@ subroutine setalb & endif fsno1 = f_one - fsno0 - flnd0 = f_one + flnd0 = min(f_one, facsf(i)+facwf(i)) fsea0 = max(f_zero, f_one-flnd0) fsno = fsno0 fsea = fsea0 * fsno1 @@ -728,6 +734,9 @@ subroutine setalb & asnvb = asnvd asnnb = asnnd endif + else + asnvb = snoalb(i) + asnnb = snoalb(i) endif !> - Calculate direct sea surface albedo, use fanglin's zenith angle From 08e918b7b5b2eda2bd51229bba50dcb8e12d29fd Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Fri, 3 May 2024 10:38:06 -0400 Subject: [PATCH 02/10] Revert changes on snow cover --- gsmphys/radiation_surface.f | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/gsmphys/radiation_surface.f b/gsmphys/radiation_surface.f index 1656940c..4fdf9f07 100644 --- a/gsmphys/radiation_surface.f +++ b/gsmphys/radiation_surface.f @@ -664,13 +664,7 @@ subroutine setalb & !> - Calculate snow cover input directly for land model, no !! conversion needed. - fsno0 = sncovr(i) - - if (nint(slmsk(i))==0 .and. - & (tsknf(i)>con_tice .or. - & ldisable_radiation_quasi_sea_ice)) then - fsno0 = f_zero - endif + fsno0 = f_zero if (nint(slmsk(i)) == 2) then asnow = 0.02*snowf(i) @@ -734,9 +728,6 @@ subroutine setalb & asnvb = asnvd asnnb = asnnd endif - else - asnvb = snoalb(i) - asnnb = snoalb(i) endif !> - Calculate direct sea surface albedo, use fanglin's zenith angle From f106601468552afe08e49802ebe1b23638661fa4 Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Wed, 8 May 2024 11:00:48 -0400 Subject: [PATCH 03/10] https://github.com/NCAR/noahmp/commit/e1b83ad4c55e5e186be1cf726a222c78d6abdcc2 --- GFS_layer/GFS_physics_driver.F90 | 2 +- GFS_layer/GFS_typedefs.F90 | 6 + gsmphys/module_sf_noahmp_glacier.f90 | 497 +++++++++++++++------------ gsmphys/sfc_noahmp_drv.f | 21 +- 4 files changed, 298 insertions(+), 228 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index f815c6ac..c2e669df 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1320,7 +1320,7 @@ subroutine GFS_physics_driver & Model%iopt_dveg, Model%iopt_crs, Model%iopt_btr, & Model%iopt_run, Model%iopt_sfc, Model%iopt_frz, & Model%iopt_inf, Model%iopt_rad, Model%iopt_alb, & - Model%iopt_snf, Model%iopt_tbot, Model%iopt_stc, & + Model%iopt_snf, Model%iopt_tbot, Model%iopt_stc, Model%iopt_gla, & grid%xlat, xcosz, Model%yearlen, Model%julian, Model%imn,& Sfcprop%drainncprv, Sfcprop%draincprv, Sfcprop%dsnowprv, & Sfcprop%dgraupelprv, Sfcprop%diceprv, & diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index d7f8110b..3d7b2e70 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -609,6 +609,7 @@ module GFS_typedefs integer :: iopt_snf !rainfall & snowfall (1-jordan91; 2->bats; 3->noah) integer :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->noah) integer :: iopt_stc !snow/soil temperature time scheme (only layer 1) + integer :: iopt_gla !glacier option (1->phase change; 2->simple) !--- tuning parameters for physical parameterizations logical :: ras !< flag for ras convection scheme @@ -2228,6 +2229,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: iopt_snf = 1 !rainfall & snowfall (1-jordan91; 2->bats; 3->noah) integer :: iopt_tbot = 2 !lower boundary of soil temperature (1->zero-flux; 2->noah) integer :: iopt_stc = 1 !snow/soil temperature time scheme (only layer 1) + integer :: iopt_gla = 1 !glacier option (1->phase change; 2->simple) !--- tuning parameters for physical parameterizations logical :: ras = .false. !< flag for ras convection scheme @@ -2479,6 +2481,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! Noah MP options iopt_dveg,iopt_crs,iopt_btr,iopt_run,iopt_sfc, iopt_frz, & iopt_inf, iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc, & + iopt_gla, & !--- physical parameterizations ras, trans_trac, old_monin, cnvgwd, mstrat, moist_adj, & cscnv, cal_pre, do_aw, do_shoc, shocaftcnv, shoc_cld, & @@ -2692,6 +2695,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%iopt_snf = iopt_snf Model%iopt_tbot = iopt_tbot Model%iopt_stc = iopt_stc + Model%iopt_gla = iopt_gla !--- tuning parameters for physical parameterizations @@ -3070,6 +3074,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & print *,'iopt_snf = ', Model%iopt_snf print *,'iopt_tbot = ',Model%iopt_tbot print *,'iopt_stc = ', Model%iopt_stc + print *,'iopt_gla = ', Model%iopt_gla @@ -3389,6 +3394,7 @@ subroutine control_print(Model) print *, ' iopt_snf : ', Model%iopt_snf print *, ' iopt_tbot : ', Model%iopt_tbot print *, ' iopt_stc : ', Model%iopt_stc + print *, ' iopt_gla : ', Model%iopt_gla endif diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index cdc251b4..6bb9eeec 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -27,57 +27,6 @@ module noahmp_glacier_globals real (kind=kind_phys), parameter :: denice = 917. !density of ice (kg/m3) ! =====================================options for different schemes================================ -! options for dynamic vegetation: -! 1 -> off (use table lai; use fveg = shdfac from input) -! 2 -> on (together with opt_crs = 1) -! 3 -> off (use table lai; calculate fveg) -! 4 -> off (use table lai; use maximum vegetation fraction) - - integer :: dveg != 2 ! - -! options for canopy stomatal resistance -! 1-> ball-berry; 2->jarvis - - integer :: opt_crs != 1 !(must 1 when dveg = 2) - -! options for soil moisture factor for stomatal resistance -! 1-> noah (soil moisture) -! 2-> clm (matric potential) -! 3-> ssib (matric potential) - - integer :: opt_btr != 1 !(suggested 1) - -! options for runoff and groundwater -! 1 -> topmodel with groundwater (niu et al. 2007 jgr) ; -! 2 -> topmodel with an equilibrium water table (niu et al. 2005 jgr) ; -! 3 -> original surface and subsurface runoff (free drainage) -! 4 -> bats surface and subsurface runoff (free drainage) - - integer :: opt_run != 1 !(suggested 1) - -! options for surface layer drag coeff (ch & cm) -! 1->m-o ; 2->original noah (chen97); 3->myj consistent; 4->ysu consistent. - - integer :: opt_sfc != 1 !(1 or 2 or 3 or 4) - -! options for supercooled liquid water (or ice fraction) -! 1-> no iteration (niu and yang, 2006 jhm); 2: koren's iteration - - integer :: opt_frz != 1 !(1 or 2) - -! options for frozen soil permeability -! 1 -> linear effects, more permeable (niu and yang, 2006, jhm) -! 2 -> nonlinear effects, less permeable (old) - - integer :: opt_inf != 1 !(suggested 1) - -! options for radiation transfer -! 1 -> modified two-stream (gap = f(solar angle, 3d structure ...)<1-fveg) -! 2 -> two-stream applied to grid-cell (gap = 0) -! 3 -> two-stream applied to vegetated fraction (gap=1-fveg) - - integer :: opt_rad != 1 !(suggested 1) - ! options for ground snow surface albedo ! 1-> bats; 2 -> class @@ -99,6 +48,10 @@ module noahmp_glacier_globals integer :: opt_stc != 1 !(suggested 1) +! options for glacier treatment +! 1 -> include phase change of ice; 2 -> ice treatment more like original noah + + integer :: opt_gla != 1 !(suggested 1) ! adjustable parameters for snow processes real (kind=kind_phys), parameter :: z0sno = 0.002 !snow surface roughness length (m) (0.002) @@ -159,7 +112,7 @@ subroutine noahmp_glacier (& fsa ,fsr ,fira ,fsh ,fgev ,ssoil , & ! out : trad ,edir ,runsrf ,runsub ,sag ,albedo , & ! out : qsnbot ,ponding ,ponding1,ponding2,t2m ,q2e , & ! out : - emissi, fpice ,ch2b , esnow ,albsnd ,albsni) + emissi, fpice ,ch2b ,qmelt ,esnow ,albsnd ,albsni) ! -------------------------------------------------------------------------------------------------- ! initial code: guo-yue niu, oct. 2007 @@ -222,7 +175,7 @@ subroutine noahmp_glacier (& real (kind=kind_phys) , intent(out) :: runsub !baseflow (saturation excess) [mm/s] real (kind=kind_phys) , intent(out) :: sag !solar rad absorbed by ground (w/m2) real (kind=kind_phys) , intent(out) :: albedo !surface albedo [-] - real (kind=kind_phys) , intent(out) :: qsnbot !snowmelt [mm/s] + real (kind=kind_phys) , intent(out) :: qsnbot !total liquid water (snowmelt + rain through pack)out of snowpack bottom [mm/s] real (kind=kind_phys) , intent(out) :: ponding!surface ponding [mm] real (kind=kind_phys) , intent(out) :: ponding1!surface ponding [mm] real (kind=kind_phys) , intent(out) :: ponding2!surface ponding [mm] @@ -234,7 +187,8 @@ subroutine noahmp_glacier (& real (kind=kind_phys) , intent(out) :: esnow real (kind=kind_phys), dimension(1:2) , intent(out) :: albsnd !snow albedo (direct) real (kind=kind_phys), dimension(1:2) , intent(out) :: albsni !snow albedo (diffuse) - + real (kind=kind_phys) , intent(out) :: qmelt !internal pack melt due to phase change [mm/s] + ! local integer :: iz !do-loop index integer, dimension(-nsnow+1:nsoil) :: imelt !phase change index [1-melt; 2-freeze] @@ -252,7 +206,6 @@ subroutine noahmp_glacier (& real (kind=kind_phys) :: qdew !ground surface dew rate [mm/s] real (kind=kind_phys) :: qvap !ground surface evap. rate [mm/s] real (kind=kind_phys) :: lathea !latent heat [j/kg] - real (kind=kind_phys) :: qmelt !internal pack melt real (kind=kind_phys) :: swdown !downward solar [w/m2] real (kind=kind_phys) :: beg_wb !beginning water for error check real (kind=kind_phys) :: zbot = -8.0 @@ -287,9 +240,9 @@ subroutine noahmp_glacier (& smc ,snice ,snliq ,albold ,cm ,ch , & !inout tauss ,qsfc , & !inout imelt ,snicev ,snliqv ,epore ,qmelt ,ponding, & !out - sag ,fsa ,fsr ,fira ,fsh ,fgev , & !out + sag ,fsa ,fsr ,fira ,fsh ,fgev , & !out trad ,t2m ,ssoil ,lathea ,q2e ,emissi , & - ch2b ,albsnd ,albsni ) !out + ch2b ,albsnd ,albsni ) !out sice = max(0.0, smc - sh2o) sneqvo = sneqv @@ -303,20 +256,25 @@ subroutine noahmp_glacier (& call water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in qvap ,qdew ,ficeold,zsoil , & !in isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout - dzsnso ,sh2o ,sice ,ponding,zsnso , & !inout + dzsnso ,sh2o ,sice ,ponding,zsnso ,fsh , & !inout runsrf ,runsub ,qsnow ,ponding1 ,ponding2,qsnbot,fpice,esnow & !out ) -! if(maxval(sice) < 0.0001) then -! write(message,*) "glacier has melted at:",iloc,jloc," are you sure this should be a glacier point?" -! call wrf_debug(10,trim(message)) -! end if + if(opt_gla == 2) then + edir = qvap - qdew + fgev = edir * lathea + end if + + if(maxval(sice) < 0.0001) then + write(message,*) "glacier has melted at:",iloc,jloc," are you sure this should be a glacier point?" + !call wrf_debug(10,trim(message)) + end if ! water and energy balance check call error_glacier (iloc ,jloc ,swdown ,fsa ,fsr ,fira , & fsh ,fgev ,ssoil ,sag ,prcp ,edir , & - runsrf ,runsub ,sneqv ,dt ,beg_wb ) + runsrf ,runsub ,sneqv ,dt ,beg_wb ) if(snowh <= 1.e-6 .or. sneqv <= 1.e-3) then snowh = 0.0 @@ -395,7 +353,7 @@ subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !i imelt ,snicev ,snliqv ,epore ,qmelt ,ponding, & !out sag ,fsa ,fsr ,fira ,fsh ,fgev , & !out trad ,t2m ,ssoil ,lathea ,q2e ,emissi , & - ch2b ,albsnd ,albsni ) !out + ch2b ,albsnd ,albsni ) !out ! -------------------------------------------------------------------------------------------------- ! -------------------------------------------------------------------------------------------------- @@ -448,7 +406,7 @@ subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !i real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(out) :: snicev !partial volume ice [m3/m3] real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(out) :: snliqv !partial volume liq. water [m3/m3] real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(out) :: epore !effective porosity [m3/m3] - real (kind=kind_phys) , intent(out) :: qmelt !snowmelt [mm/s] + real (kind=kind_phys) , intent(out) :: qmelt !snowmelt due to phase change [mm/s] real (kind=kind_phys) , intent(out) :: ponding!pounding at ground [mm] real (kind=kind_phys) , intent(out) :: sag !solar rad. absorbed by ground (w/m2) real (kind=kind_phys) , intent(out) :: fsa !tot. absorbed solar radiation (w/m2) @@ -526,11 +484,11 @@ subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !i call glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z0mg , & !in zlvl ,zpd ,qair ,sfctmp ,rhoair ,sfcprs , & !in - ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in - eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in - cm ,ch ,tg ,qsfc , & !inout - fira ,fsh ,fgev ,ssoil , & !out - t2m ,q2e ,ch2b) !out + ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in + eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in + cm ,ch ,tg ,qsfc , & !inout + fira ,fsh ,fgev ,ssoil , & !out + t2m ,q2e ,ch2b) !out !energy balance at surface: sag=(irb+shb+evb+ghb) @@ -551,7 +509,7 @@ subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !i call tsnosoi_glacier (nsoil ,nsnow ,isnow ,dt ,tbot , & !in ssoil ,snowh ,zbot ,zsnso ,df , & !in - hcpct , & !in + hcpct , & !in stc ) !inout ! adjusting snow surface temperature @@ -559,7 +517,7 @@ subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !i if (snowh > 0.05 .and. tg > tfrz) tg = tfrz end if -! energy released or consumed by snow & frozen soil +! energy released or consumed by snow & ice call phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & !in dzsnso , & !in @@ -747,18 +705,23 @@ subroutine radiation_glacier (dt ,tg ,sneqvo ,sneqv ,cosz , & !i albice(1) = 0.80 !albedo land ice: 1=vis, 2=nir albice(2) = 0.55 -! snow age +! compute snow albedo only if cosz >0 +! to be consistent with the main noahmp code, currently do not include snow aging when sun is not present +! this needs more future work + if (cosz > 0.0) then - call snow_age_glacier (dt,tg,sneqvo,sneqv,tauss,fage) + ! snow age -! snow albedos: age even when sun is not present + call snow_age_glacier (dt,tg,sneqvo,sneqv,tauss,fage) - if(opt_alb == 1) & - call snowalb_bats_glacier (nband,cosz,fage,albsnd,albsni) - if(opt_alb == 2) then - call snowalb_class_glacier(nband,qsnow,dt,alb,albold,albsnd,albsni) - albold = alb - end if + if(opt_alb == 1) & + call snowalb_bats_glacier (nband,cosz,fage,albsnd,albsni) + if(opt_alb == 2) then + call snowalb_class_glacier(nband,qsnow,dt,alb,albold,albsnd,albsni) + albold = alb + end if + + endif ! zero summed solar fluxes @@ -939,8 +902,8 @@ end subroutine snowalb_class_glacier ! ================================================================================================== subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z0m , & !in zlvl ,zpd ,qair ,sfctmp ,rhoair ,sfcprs , & !in - ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in - eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in + ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in + eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in cm ,ch ,tgb ,qsfc , & !inout irb ,shb ,evb ,ghb , & !out t2mb ,q2b ,ehb2) !out @@ -1081,7 +1044,11 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z end if csh = rhoair*cpair/rahb - cev = rhoair*cpair/gamma/(rsurf+rawb) + if(snowh > 0.0 .or. opt_gla == 1) then + cev = rhoair*cpair/gamma/(rsurf+rawb) + else + cev = 0.0 ! don't allow any sublimation of glacier in opt_gla=2 + end if ! surface fluxes and dtg @@ -1120,9 +1087,13 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z ! if snow on ground and tg > tfrz: reset tg = tfrz. reevaluate ground fluxes. sice = smc - sh2o - if(opt_stc == 1) then - if ((maxval(sice) > 0.0 .or. snowh > 0.0) .and. tgb > tfrz) then + if(opt_stc == 1 .or. opt_stc ==3) then + if ((maxval(sice) > 0.0 .or. snowh > 0.0) .and. tgb > tfrz .and. opt_gla == 1) then tgb = tfrz + t = tdc(tgb) ! mb: recalculate estg + call esat(t, esatw, esati, dsatw, dsati) + estg = esati + qsfc = 0.622*(estg*rhsur)/(sfcprs-0.378*(estg*rhsur)) irb = cir * tgb**4 - emg*lwdn shb = csh * (tgb - sfctmp) evb = cev * (estg*rhsur - eair ) !estg reevaluate ? @@ -1227,9 +1198,10 @@ subroutine sfcdif1_glacier(iter ,zlvl ,zpd ,z0h ,z0m , & !in real (kind=kind_phys), intent(inout) :: fh !sen heat stability correction, weighted by prior iters real (kind=kind_phys), intent(inout) :: fm2 !sen heat stability correction, weighted by prior iters real (kind=kind_phys), intent(inout) :: fh2 !sen heat stability correction, weighted by prior iters + real (kind=kind_phys), intent(inout) :: fv !friction velocity (m/s) ! outputs - real (kind=kind_phys), intent(out) :: fv !friction velocity (m/s) + real (kind=kind_phys), intent(out) :: cm !drag coefficient for momentum real (kind=kind_phys), intent(out) :: ch !drag coefficient for heat real (kind=kind_phys), intent(out) :: ch2 !drag coefficient for heat @@ -1269,7 +1241,7 @@ subroutine sfcdif1_glacier(iter ,zlvl ,zpd ,z0h ,z0m , & !in tmpch2 = log((2.0 + z0h) / z0h) if(iter == 1) then - fv = 0.1 + fv = 0.0 moz = 0.0 mol = 0.0 moz2 = 0.0 @@ -1278,7 +1250,6 @@ subroutine sfcdif1_glacier(iter ,zlvl ,zpd ,z0h ,z0m , & !in tmp1 = vkc * (grav/tvir) * h/(rhoair*cpair) if (abs(tmp1) .le. mpe) tmp1 = mpe mol = -1. * fv**3 / tmp1 - if (abs(mol) .le. mpe) mol = mpe moz = min( (zlvl-zpd)/mol, 1.) moz2 = min( (2.0 + z0h)/mol, 1.) endif @@ -1359,7 +1330,7 @@ end subroutine sfcdif1_glacier ! ================================================================================================== subroutine tsnosoi_glacier (nsoil ,nsnow ,isnow ,dt ,tbot , & !in ssoil ,snowh ,zbot ,zsnso ,df , & !in - hcpct , & !in + hcpct , & !in stc ) !inout ! -------------------------------------------------------------------------------------------------- ! compute snow (up to 3l) and soil (4l) temperature. note that snow temperatures @@ -1497,7 +1468,7 @@ subroutine hrt_glacier (nsnow ,nsoil ,isnow ,zsnso , & !in if (k == isnow+1) then ai(k) = 0.0 ci(k) = - df(k) * ddz(k) / denom(k) - if (opt_stc == 1) then + if (opt_stc == 1 .or. opt_stc == 3) then bi(k) = - ci(k) end if if (opt_stc == 2) then @@ -1693,12 +1664,106 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & mliq(j) = snliq(j) end do + do j = isnow+1,0 ! all snow layers; do ice later + imelt(j) = 0 + hm(j) = 0. + xm(j) = 0. + wice0(j) = mice(j) + wliq0(j) = mliq(j) + wmass0(j) = mice(j) + mliq(j) + enddo + + do j = isnow+1,0 + if (mice(j) > 0. .and. stc(j) >= tfrz) then ! melting + imelt(j) = 1 + endif + if (mliq(j) > 0. .and. stc(j) < tfrz) then ! freezing + imelt(j) = 2 + endif + + enddo + +! calculate the energy surplus and loss for melting and freezing + + do j = isnow+1,0 + if (imelt(j) > 0) then + hm(j) = (stc(j)-tfrz)/fact(j) + stc(j) = tfrz + endif + + if (imelt(j) == 1 .and. hm(j) < 0.) then + hm(j) = 0. + imelt(j) = 0 + endif + if (imelt(j) == 2 .and. hm(j) > 0.) then + hm(j) = 0. + imelt(j) = 0 + endif + xm(j) = hm(j)*dt/hfus + enddo + +! the rate of melting and freezing for snow without a layer, opt_gla==1 treated below + +if (opt_gla == 2) then + + if (isnow == 0 .and. sneqv > 0. .and. stc(1) >= tfrz) then + hm(1) = (stc(1)-tfrz)/fact(1) ! available heat + stc(1) = tfrz ! set t to freezing + xm(1) = hm(1)*dt/hfus ! total snow melt possible + + temp1 = sneqv + sneqv = max(0.,temp1-xm(1)) ! snow remaining + propor = sneqv/temp1 ! fraction melted + snowh = max(0.,propor * snowh) ! new snow height + heatr(1) = hm(1) - hfus*(temp1-sneqv)/dt ! excess heat + if (heatr(1) > 0.) then + xm(1) = heatr(1)*dt/hfus + stc(1) = stc(1) + fact(1)*heatr(1) ! re-heat ice + else + xm(1) = 0. ! heat used up + hm(1) = 0. + endif + qmelt = max(0.,(temp1-sneqv))/dt ! melted snow rate + xmf = hfus*qmelt ! melted snow energy + ponding = temp1-sneqv ! melt water + endif + +end if ! opt_gla == 2 + +! the rate of melting and freezing for snow + + do j = isnow+1,0 + if (imelt(j) > 0 .and. abs(hm(j)) > 0.) then + + heatr(j) = 0. + if (xm(j) > 0.) then + mice(j) = max(0., wice0(j)-xm(j)) + heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt + else if (xm(j) < 0.) then + mice(j) = min(wmass0(j), wice0(j)-xm(j)) + heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt + endif + + mliq(j) = max(0.,wmass0(j)-mice(j)) + + if (abs(heatr(j)) > 0.) then + stc(j) = stc(j) + fact(j)*heatr(j) + if (mliq(j)*mice(j)>0.) stc(j) = tfrz + endif + + qmelt = qmelt + max(0.,(wice0(j)-mice(j)))/dt + + endif + enddo + +if (opt_gla == 1) then ! operate on the ice layers + do j = 1, nsoil ! all soil layers mliq(j) = sh2o(j) * dzsnso(j) * 1000. mice(j) = (smc(j) - sh2o(j)) * dzsnso(j) * 1000. end do - do j = isnow+1,nsoil ! all layers + do j = 1,nsoil ! all layers imelt(j) = 0 hm(j) = 0. xm(j) = 0. @@ -1707,7 +1772,7 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & wmass0(j) = mice(j) + mliq(j) enddo - do j = isnow+1,nsoil + do j = 1,nsoil if (mice(j) > 0. .and. stc(j) >= tfrz) then ! melting imelt(j) = 1 endif @@ -1725,7 +1790,7 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & ! calculate the energy surplus and loss for melting and freezing - do j = isnow+1,nsoil + do j = 1,nsoil if (imelt(j) > 0) then hm(j) = (stc(j)-tfrz)/fact(j) stc(j) = tfrz @@ -1753,20 +1818,20 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (heatr(1) > 0.) then xm(1) = heatr(1)*dt/hfus hm(1) = heatr(1) - imelt(1) = 1 + imelt(1) = 1 else xm(1) = 0. hm(1) = 0. - imelt(1) = 0 + imelt(1) = 0 endif qmelt = max(0.,(temp1-sneqv))/dt xmf = hfus*qmelt ponding = temp1-sneqv endif -! the rate of melting and freezing for snow and soil +! the rate of melting and freezing for soil - do j = isnow+1,nsoil + do j = 1,nsoil if (imelt(j) > 0 .and. abs(hm(j)) > 0.) then heatr(j) = 0. @@ -1804,21 +1869,21 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz)) then do j = 1,nsoil if ( stc(j) > tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) do k = 1,nsoil - if (j .ne. k .and. stc(k) < tfrz .and. heatr(j) > 0.1) then - heatr(k) = (stc(k)-tfrz)/fact(k) - if (abs(heatr(k)) > heatr(j)) then ! layer absorbs all - heatr(k) = heatr(k) + heatr(j) - stc(k) = tfrz + heatr(k)*fact(k) - heatr(j) = 0.0 + if (j .ne. k .and. stc(k) < tfrz .and. heatr(j) > 0.1) then + heatr(k) = (stc(k)-tfrz)/fact(k) + if (abs(heatr(k)) > heatr(j)) then ! layer absorbs all + heatr(k) = heatr(k) + heatr(j) + stc(k) = tfrz + heatr(k)*fact(k) + heatr(j) = 0.0 else - heatr(j) = heatr(j) + heatr(k) - heatr(k) = 0.0 - stc(k) = tfrz + heatr(j) = heatr(j) + heatr(k) + heatr(k) = 0.0 + stc(k) = tfrz end if - end if - end do + end if + end do stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -1829,21 +1894,21 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz)) then do j = 1,nsoil if ( stc(j) < tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) do k = 1,nsoil - if (j .ne. k .and. stc(k) > tfrz .and. heatr(j) < -0.1) then - heatr(k) = (stc(k)-tfrz)/fact(k) - if (heatr(k) > abs(heatr(j))) then ! layer absorbs all - heatr(k) = heatr(k) + heatr(j) - stc(k) = tfrz + heatr(k)*fact(k) - heatr(j) = 0.0 + if (j .ne. k .and. stc(k) > tfrz .and. heatr(j) < -0.1) then + heatr(k) = (stc(k)-tfrz)/fact(k) + if (heatr(k) > abs(heatr(j))) then ! layer absorbs all + heatr(k) = heatr(k) + heatr(j) + stc(k) = tfrz + heatr(k)*fact(k) + heatr(j) = 0.0 else - heatr(j) = heatr(j) + heatr(k) - heatr(k) = 0.0 - stc(k) = tfrz + heatr(j) = heatr(j) + heatr(k) + heatr(k) = 0.0 + stc(k) = tfrz end if - end if - end do + end if + end do stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -1854,25 +1919,25 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) > tfrz) .and. any(mice(1:4) > 0.)) then do j = 1,nsoil if ( stc(j) > tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) xm(j) = heatr(j)*dt/hfus do k = 1,nsoil - if (j .ne. k .and. mice(k) > 0. .and. xm(j) > 0.1) then - if (mice(k) > xm(j)) then ! layer absorbs all - mice(k) = mice(k) - xm(j) - xmf = xmf + hfus * xm(j)/dt - stc(k) = tfrz - xm(j) = 0.0 + if (j .ne. k .and. mice(k) > 0. .and. xm(j) > 0.1) then + if (mice(k) > xm(j)) then ! layer absorbs all + mice(k) = mice(k) - xm(j) + xmf = xmf + hfus * xm(j)/dt + stc(k) = tfrz + xm(j) = 0.0 else - xm(j) = xm(j) - mice(k) - xmf = xmf + hfus * mice(k)/dt - mice(k) = 0.0 - stc(k) = tfrz + xm(j) = xm(j) - mice(k) + xmf = xmf + hfus * mice(k)/dt + mice(k) = 0.0 + stc(k) = tfrz end if mliq(k) = max(0.,wmass0(k)-mice(k)) - end if - end do - heatr(j) = xm(j)*hfus/dt + end if + end do + heatr(j) = xm(j)*hfus/dt stc(j) = tfrz + heatr(j)*fact(j) end if end do @@ -1883,29 +1948,31 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & if (any(stc(1:4) < tfrz) .and. any(mliq(1:4) > 0.)) then do j = 1,nsoil if ( stc(j) < tfrz ) then - heatr(j) = (stc(j)-tfrz)/fact(j) + heatr(j) = (stc(j)-tfrz)/fact(j) xm(j) = heatr(j)*dt/hfus do k = 1,nsoil - if (j .ne. k .and. mliq(k) > 0. .and. xm(j) < -0.1) then - if (mliq(k) > abs(xm(j))) then ! layer absorbs all - mice(k) = mice(k) - xm(j) - xmf = xmf + hfus * xm(j)/dt - stc(k) = tfrz - xm(j) = 0.0 + if (j .ne. k .and. mliq(k) > 0. .and. xm(j) < -0.1) then + if (mliq(k) > abs(xm(j))) then ! layer absorbs all + mice(k) = mice(k) - xm(j) + xmf = xmf + hfus * xm(j)/dt + stc(k) = tfrz + xm(j) = 0.0 else - xm(j) = xm(j) + mliq(k) - xmf = xmf - hfus * mliq(k)/dt - mice(k) = wmass0(k) - stc(k) = tfrz + xm(j) = xm(j) + mliq(k) + xmf = xmf - hfus * mliq(k)/dt + mice(k) = wmass0(k) + stc(k) = tfrz end if mliq(k) = max(0.,wmass0(k)-mice(k)) - end if - end do - heatr(j) = xm(j)*hfus/dt + end if + end do + heatr(j) = xm(j)*hfus/dt stc(j) = tfrz + heatr(j)*fact(j) end if end do end if + +end if ! opt_gla == 1 do j = isnow+1,0 ! snow snliq(j) = mliq(j) @@ -1913,10 +1980,14 @@ subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & end do do j = 1, nsoil ! soil + if(opt_gla == 1) then sh2o(j) = mliq(j) / (1000. * dzsnso(j)) sh2o(j) = max(0.0,min(1.0,sh2o(j))) ! smc(j) = (mliq(j) + mice(j)) / (1000. * dzsnso(j)) - smc(j) = 1.0 + elseif(opt_gla == 2) then + sh2o(j) = 0.0 ! ice, assume all frozen...forever + end if + smc(j) = 1.0 end do end subroutine phasechange_glacier @@ -1924,7 +1995,7 @@ end subroutine phasechange_glacier subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in qvap ,qdew ,ficeold,zsoil , & !in isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout - dzsnso ,sh2o ,sice ,ponding,zsnso , & !inout + dzsnso ,sh2o ,sice ,ponding,zsnso ,fsh , & !inout runsrf ,runsub ,qsnow ,ponding1 ,ponding2,qsnbot,fpice,esnow & !out ) !out ! ---------------------------------------------------------------------- @@ -1940,8 +2011,8 @@ subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in real (kind=kind_phys), intent(in) :: dt !main time step (s) real (kind=kind_phys), intent(in) :: prcp !precipitation (mm/s) real (kind=kind_phys), intent(in) :: sfctmp !surface air temperature [k] - real (kind=kind_phys), intent(in) :: qvap !soil surface evaporation rate[mm/s] - real (kind=kind_phys), intent(in) :: qdew !soil surface dew rate[mm/s] + real (kind=kind_phys), intent(inout) :: qvap !soil surface evaporation rate[mm/s] + real (kind=kind_phys), intent(inout) :: qdew !soil surface dew rate[mm/s] real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: ficeold !ice fraction at last timestep real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: zsoil !layer-bottom depth from soil surf (m) @@ -1957,6 +2028,7 @@ subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sice !soil ice content [m3/m3] real (kind=kind_phys) , intent(inout) :: ponding ![mm] real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: zsnso !layer-bottom depth from snow surf [m] + real (kind=kind_phys), intent(inout) :: fsh !total sensible heat (w/m2) [+ to atm] ! output real (kind=kind_phys), intent(out) :: runsrf !surface runoff [mm/s] @@ -1964,7 +2036,7 @@ subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in real (kind=kind_phys), intent(out) :: qsnow !snow at ground srf (mm/s) [+] real (kind=kind_phys), intent(out) :: ponding1 real (kind=kind_phys), intent(out) :: ponding2 - real (kind=kind_phys), intent(out) :: qsnbot !melting water out of snow bottom [mm/s] + real (kind=kind_phys), intent(out) :: qsnbot !total liquid water (snowmelt + rain through pack)out of snowpack bottom [mm/s] real (kind=kind_phys), intent(out) :: fpice !precipitation frozen fraction real (kind=kind_phys), intent(out) :: esnow ! @@ -2040,35 +2112,17 @@ subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in ! sublimation, frost, evaporation, and dew -! qsnsub = 0. -! if (sneqv > 0.) then -! qsnsub = min(qvap, sneqv/dt) -! endif -! qseva = qvap-qsnsub - -! qsnfro = 0. -! if (sneqv > 0.) then -! qsnfro = qdew -! endif -! qsdew = qdew - qsnfro - qsnsub = qvap ! send total sublimation/frost to snowwater and deal with it there qsnfro = qdew esnow = qsnsub*2.83e+6 -! print *, 'qvap',qvap,qvap*dt -! print *, 'qsnsub',qsnsub,qsnsub*dt -! print *, 'qseva',qseva,qseva*dt -! print *, 'qsnfro',qsnfro,qsnfro*dt -! print *, 'qdew',qdew,qdew*dt -! print *, 'qsdew',qsdew,qsdew*dt -!print *, 'before snowwater', sneqv,snowh,snice,snliq,sh2o,sice call snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in snowhin,qsnow ,qsnfro ,qsnsub ,qrain , & !in ficeold,zsoil , & !in isnow ,snowh ,sneqv ,snice ,snliq , & !inout sh2o ,sice ,stc ,dzsnso ,zsnso , & !inout + fsh , & !inout qsnbot ,snoflow,ponding1 ,ponding2) !out !print *, 'after snowwater', sneqv,snowh,snice,snliq,sh2o,sice !print *, 'ponding', ponding,ponding1,ponding2 @@ -2084,20 +2138,30 @@ subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in endif - replace = 0.0 - do ilev = 1,nsoil + if(opt_gla == 1) then + replace = 0.0 + do ilev = 1,nsoil replace = replace + dzsnso(ilev)*(sice(ilev) - sice_save(ilev) + sh2o(ilev) - sh2o_save(ilev)) - end do - replace = replace * 1000.0 / dt ! convert to [mm/s] + end do + replace = replace * 1000.0 / dt ! convert to [mm/s] - sice = min(1.0,sice_save) + sice = min(1.0,sice_save) + elseif(opt_gla == 2) then + sice = 1.0 + end if sh2o = 1.0 - sice -!print *, 'replace', replace + ! use runsub as a water balancer, snoflow is snow that disappears, replace is ! water from below that replaces glacier loss - runsub = snoflow + replace + if(opt_gla == 1) then + runsub = snoflow + replace + elseif(opt_gla == 2) then + runsub = snoflow + qvap = qsnsub + qdew = qsnfro + end if end subroutine water_glacier ! ================================================================================================== @@ -2107,6 +2171,7 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in ficeold,zsoil , & !in isnow ,snowh ,sneqv ,snice ,snliq , & !inout sh2o ,sice ,stc ,dzsnso ,zsnso , & !inout + fsh , & !inout qsnbot ,snoflow,ponding1 ,ponding2) !out ! ---------------------------------------------------------------------- implicit none @@ -2119,8 +2184,8 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in real (kind=kind_phys), intent(in) :: sfctmp !surface air temperature [k] real (kind=kind_phys), intent(in) :: snowhin!snow depth increasing rate (m/s) real (kind=kind_phys), intent(in) :: qsnow !snow at ground srf (mm/s) [+] - real (kind=kind_phys), intent(in) :: qsnfro !snow surface frost rate[mm/s] - real (kind=kind_phys), intent(in) :: qsnsub !snow surface sublimation rate[mm/s] + real (kind=kind_phys), intent(inout) :: qsnfro !snow surface frost rate[mm/s] + real (kind=kind_phys), intent(inout) :: qsnsub !snow surface sublimation rate[mm/s] real (kind=kind_phys), intent(in) :: qrain !snow surface rain rate[mm/s] real (kind=kind_phys), dimension(-nsnow+1:0) , intent(in) :: ficeold!ice fraction at last timestep real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: zsoil !layer-bottom depth from soil surf (m) @@ -2136,9 +2201,10 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc !snow layer temperature [k] real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso !snow/soil layer thickness [m] real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: zsnso !layer-bottom depth from snow surf [m] + real (kind=kind_phys), intent(inout) :: fsh !total sensible heat (w/m2) [+ to atm] ! output - real (kind=kind_phys), intent(out) :: qsnbot !melting water out of snow bottom [mm/s] + real (kind=kind_phys), intent(out) :: qsnbot !total liquid water (snowmelt + rain through pack)out of snowpack bottom [mm/s] real (kind=kind_phys), intent(out) :: snoflow!glacier flow [mm] real (kind=kind_phys), intent(out) :: ponding1 real (kind=kind_phys), intent(out) :: ponding2 @@ -2184,14 +2250,14 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in qrain , & !in isnow ,dzsnso ,snowh ,sneqv ,snice , & !inout snliq ,sh2o ,sice ,stc , & !inout - ponding1 ,ponding2 , & !inout + ponding1 ,ponding2 ,fsh , & !inout qsnbot ) !out !to obtain equilibrium state of snow in glacier region - if(sneqv > 2000.) then ! 2000 mm -> maximum water depth + if(sneqv > 5000.0) then ! 5000 mm -> maximum water depth bdsnow = snice(0) / dzsnso(0) - snoflow = (sneqv - 2000.) + snoflow = (sneqv - 5000.0) snice(0) = snice(0) - snoflow dzsnso(0) = dzsnso(0) - snoflow/bdsnow snoflow = snoflow / dt @@ -2744,7 +2810,7 @@ subroutine snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in qrain , & !in isnow ,dzsnso ,snowh ,sneqv ,snice , & !inout snliq ,sh2o ,sice ,stc , & !inout - ponding1 ,ponding2 , & !inout + ponding1 ,ponding2 ,fsh , & !inout qsnbot ) !out ! ---------------------------------------------------------------------- ! renew the mass of ice lens (snice) and liquid (snliq) of the @@ -2757,13 +2823,13 @@ subroutine snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in integer, intent(in) :: nsnow !maximum no. of snow layers[=3] integer, intent(in) :: nsoil !no. of soil layers[=4] real (kind=kind_phys), intent(in) :: dt !time step - real (kind=kind_phys), intent(in) :: qsnfro !snow surface frost rate[mm/s] - real (kind=kind_phys), intent(in) :: qsnsub !snow surface sublimation rate[mm/s] + real (kind=kind_phys), intent(inout) :: qsnfro !snow surface frost rate[mm/s] + real (kind=kind_phys), intent(inout) :: qsnsub !snow surface sublimation rate[mm/s] real (kind=kind_phys), intent(in) :: qrain !snow surface rain rate[mm/s] ! output - real (kind=kind_phys), intent(out) :: qsnbot !melting water out of snow bottom [mm/s] + real (kind=kind_phys), intent(out) :: qsnbot !total liquid water (snowmelt + rain through pack)out of snowpack bottom [mm/s] ! input and output @@ -2778,6 +2844,7 @@ subroutine snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc !snow layer temperature [k] real (kind=kind_phys), intent(inout) :: ponding1 real (kind=kind_phys), intent(inout) :: ponding2 + real (kind=kind_phys), intent(inout) :: fsh !total sensible heat (w/m2) [+ to atm] ! local variables: @@ -2794,7 +2861,13 @@ subroutine snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in !for the case when sneqv becomes '0' after 'combine' if(sneqv == 0.) then - sice(1) = sice(1) + (qsnfro-qsnsub)*dt/(dzsnso(1)*1000.) + if(opt_gla == 1) then + sice(1) = sice(1) + (qsnfro-qsnsub)*dt/(dzsnso(1)*1000.) + elseif(opt_gla == 2) then + fsh = fsh - (qsnfro-qsnsub)*hsub + qsnfro = 0.0 + qsnsub = 0.0 + end if end if ! for shallow snow without a layer @@ -2803,10 +2876,16 @@ subroutine snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in ! to aviod this problem. if(isnow == 0 .and. sneqv > 0.) then - temp = sneqv - sneqv = sneqv - qsnsub*dt + qsnfro*dt - propor = sneqv/temp - snowh = max(0.,propor * snowh) + if(opt_gla == 1) then + temp = sneqv + sneqv = sneqv - qsnsub*dt + qsnfro*dt + propor = sneqv/temp + snowh = max(0.,propor * snowh) + elseif(opt_gla == 2) then + fsh = fsh - (qsnfro-qsnsub)*hsub + qsnfro = 0.0 + qsnsub = 0.0 + end if if(sneqv < 0.) then sice(1) = sice(1) + sneqv/(dzsnso(1)*1000.) @@ -2889,7 +2968,7 @@ end subroutine snowh2o_glacier ! ================================================================================================== subroutine error_glacier (iloc ,jloc ,swdown ,fsa ,fsr ,fira , & fsh ,fgev ,ssoil ,sag ,prcp ,edir , & - runsrf ,runsub ,sneqv ,dt ,beg_wb ) + runsrf ,runsub ,sneqv ,dt ,beg_wb ) ! -------------------------------------------------------------------------------------------------- ! check surface energy balance and water balance ! -------------------------------------------------------------------------------------------------- @@ -2947,41 +3026,24 @@ subroutine error_glacier (iloc ,jloc ,swdown ,fsa ,fsr ,fira , & end subroutine error_glacier ! ================================================================================================== - subroutine noahmp_options_glacier(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & - iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc ) + subroutine noahmp_options_glacier(iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iopt_gla ) implicit none - integer, intent(in) :: idveg !dynamic vegetation (1 -> off ; 2 -> on) with opt_crs = 1 - integer, intent(in) :: iopt_crs !canopy stomatal resistance (1-> ball-berry; 2->jarvis) - integer, intent(in) :: iopt_btr !soil moisture factor for stomatal resistance (1-> noah; 2-> clm; 3-> ssib) - integer, intent(in) :: iopt_run !runoff and groundwater (1->simgm; 2->simtop; 3->schaake96; 4->bats) - integer, intent(in) :: iopt_sfc !surface layer drag coeff (ch & cm) (1->m-o; 2->chen97) - integer, intent(in) :: iopt_frz !supercooled liquid water (1-> ny06; 2->koren99) - integer, intent(in) :: iopt_inf !frozen soil permeability (1-> ny06; 2->koren99) - integer, intent(in) :: iopt_rad !radiation transfer (1->gap=f(3d,cosz); 2->gap=0; 3->gap=1-fveg) integer, intent(in) :: iopt_alb !snow surface albedo (1->bats; 2->class) integer, intent(in) :: iopt_snf !rainfall & snowfall (1-jordan91; 2->bats; 3->noah) integer, intent(in) :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->noah) - integer, intent(in) :: iopt_stc !snow/soil temperature time scheme (only layer 1) ! 1 -> semi-implicit; 2 -> full implicit (original noah) + integer, intent(in) :: iopt_gla ! glacier option (1->phase change; 2->simple) ! ------------------------------------------------------------------------------------------------- - dveg = idveg - - opt_crs = iopt_crs - opt_btr = iopt_btr - opt_run = iopt_run - opt_sfc = iopt_sfc - opt_frz = iopt_frz - opt_inf = iopt_inf - opt_rad = iopt_rad opt_alb = iopt_alb opt_snf = iopt_snf opt_tbot = iopt_tbot opt_stc = iopt_stc + opt_gla = iopt_gla end subroutine noahmp_options_glacier @@ -2995,3 +3057,4 @@ module module_sf_noahmp_glacier end module module_sf_noahmp_glacier + diff --git a/gsmphys/sfc_noahmp_drv.f b/gsmphys/sfc_noahmp_drv.f index 1359edab..32cd1ae2 100644 --- a/gsmphys/sfc_noahmp_drv.f +++ b/gsmphys/sfc_noahmp_drv.f @@ -9,6 +9,7 @@ subroutine noahmpdrv & & shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, & & idveg,iopt_crs, iopt_btr, iopt_run, iopt_sfc, iopt_frz, & & iopt_inf,iopt_rad, iopt_alb, iopt_snf,iopt_tbot,iopt_stc, & + & iopt_gla, & & xlatin,xcoszin, iyrlen, julian,imon, & & rainn_mp,rainc_mp,snow_mp,graupel_mp,ice_mp, & @@ -97,7 +98,8 @@ subroutine noahmpdrv & integer, intent(in) :: idveg, iopt_crs,iopt_btr,iopt_run, & & iopt_sfc,iopt_frz,iopt_inf,iopt_rad, & - & iopt_alb,iopt_snf,iopt_tbot,iopt_stc + & iopt_alb,iopt_snf,iopt_tbot,iopt_stc, & + & iopt_gla real (kind=kind_phys), intent(in) :: julian integer, intent(in) :: iyrlen @@ -173,7 +175,7 @@ subroutine noahmpdrv & & sfcems, sheat, shdfac, shdmin1d, shdmax1d, smcwlt, & & smcdry, smcref, smcmax, sneqv, snoalb1d, snowh, & & snomlt, sncovr, soilw, soilm, ssoil, tsea, th2, & - & xlai, zlvl, swdn, tem, psfc,fdown,t2v,tbot + & xlai, zlvl, swdn, tem, psfc,fdown,t2v,tbot, qmelt real (kind=kind_phys) :: pconv,pnonc,pshcv,psnow,pgrpl,phail real (kind=kind_phys) :: lat,cosz,uu,vv,swe @@ -555,7 +557,6 @@ subroutine noahmpdrv & !-- old ! do k = 1, km -! stsoil(k) = stc(i,k) smsoil(k) = smc(i,k) slsoil(k) = slc(i,k) enddo @@ -585,8 +586,7 @@ subroutine noahmpdrv & tbot = min(tbot,263.15) call noahmp_options_glacier & - & (idveg ,iopt_crs ,iopt_btr, iopt_run ,iopt_sfc ,iopt_frz, & - & iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc ) + & (iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iopt_gla ) call noahmp_glacier ( & & i ,1 ,cosz ,nsnow ,nsoil ,delt , & ! in : time/space/model-related @@ -598,7 +598,8 @@ subroutine noahmpdrv & & fsa ,fsr ,fira ,fsh ,fgev ,ssoil , & ! out : & trad ,edir ,runsrf ,runsub ,sag ,albedo , & ! out : albedo is surface albedo & qsnbot ,ponding ,ponding1,ponding2,t2mb ,q2b , & ! out : - & emissi ,fpice ,ch2b ,esnow ,albd, albi) + & emissi ,fpice ,ch2b ,qmelt ,esnow ,albd , & + & albi) ! ! in/out and outs @@ -799,10 +800,10 @@ subroutine noahmpdrv & sfcemis(i) = emissi if(albedo .gt. 0.0) then sfalb(i) = albedo - albdvis(i) = albd(1) - albdnir(i) = albd(2) - albivis(i) = albi(1) - albinir(i) = albi(2) + albdvis(i) = albd(1) + albdnir(i) = albd(2) + albivis(i) = albi(1) + albinir(i) = albi(2) end if stm(i) = (0.1*smsoil(1)+0.3*smsoil(2)+0.6*smsoil(3)+ & From 03766b5ab749f4983453eedb7c88e51fa5047c45 Mon Sep 17 00:00:00 2001 From: Spencer Clark Date: Mon, 20 May 2024 12:27:34 -0400 Subject: [PATCH 04/10] https://github.com/NCAR/ccpp-physics/commit/3407bf1bc11c4c37e6bf28307f3ce1ae178d9a55 --- gsmphys/module_sf_noahmp_glacier.f90 | 1 + gsmphys/module_sf_noahmplsm.f90 | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index 6bb9eeec..88c16b9e 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -2448,6 +2448,7 @@ subroutine compact_glacier (nsnow ,nsoil ,dt ,stc ,snice , & !in ! the change in dz due to compaction dzsnso(j) = dzsnso(j)*(1.+pdzdtc) + dzsnso(j) = min(max(dzsnso(j),(snliq(j)+snice(j))/500.0),(snliq(j)+snice(j))/50.0) ! limit adjustment to a reasonable density end if ! pressure of overlying snow diff --git a/gsmphys/module_sf_noahmplsm.f90 b/gsmphys/module_sf_noahmplsm.f90 index ccc34ef1..47201b29 100644 --- a/gsmphys/module_sf_noahmplsm.f90 +++ b/gsmphys/module_sf_noahmplsm.f90 @@ -6472,6 +6472,7 @@ subroutine compact (parameters,nsnow ,nsoil ,dt ,stc ,snice , & !in ! the change in dz due to compaction dzsnso(j) = dzsnso(j)*(1.+pdzdtc) + dzsnso(j) = min(max(dzsnso(j),(snliq(j)+snice(j))/500.0),(snliq(j)+snice(j))/50.0) ! limit adjustment to a reasonable density end if ! pressure of overlying snow @@ -6629,6 +6630,10 @@ subroutine snowh2o (parameters,nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in end if end do + do j = isnow+1, 0 + dzsnso(j) = min(max(dzsnso(j),(snliq(j)+snice(j))/500.0),(snliq(j)+snice(j))/50.0) ! limit adjustment to a reasonable density + end do + ! liquid water from snow bottom to soil qsnbot = qout / dt ! mm/s From 47553147fc180efba3b146942885a22ff6eb0dca Mon Sep 17 00:00:00 2001 From: Spencer Clark Date: Mon, 20 May 2024 12:19:17 -0400 Subject: [PATCH 05/10] https://github.com/NCAR/ccpp-physics/commit/d344faa1c6b42a8faeadd4b359b90b817393afd9 --- gsmphys/module_sf_noahmp_glacier.f90 | 2 ++ gsmphys/module_sf_noahmplsm.f90 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index 88c16b9e..ffcfb80f 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -2267,8 +2267,10 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in if(isnow /= 0) then sneqv = 0. + snowh = 0. do iz = isnow+1,0 sneqv = sneqv + snice(iz) + snliq(iz) + snowh = snowh + dzsnso(iz) enddo end if diff --git a/gsmphys/module_sf_noahmplsm.f90 b/gsmphys/module_sf_noahmplsm.f90 index 47201b29..90e015df 100644 --- a/gsmphys/module_sf_noahmplsm.f90 +++ b/gsmphys/module_sf_noahmplsm.f90 @@ -5912,8 +5912,10 @@ subroutine snowwater (parameters,nsnow ,nsoil ,imelt ,dt ,zsoil , & !in if(isnow < 0) then ! mb: only do for multi-layer sneqv = 0. + snowh = 0. do iz = isnow+1,0 sneqv = sneqv + snice(iz) + snliq(iz) + snowh = snowh + dzsnso(iz) enddo end if From 8d9e9a66ab808f30af16ed9b61837a1d59f7477b Mon Sep 17 00:00:00 2001 From: Spencer Clark Date: Mon, 20 May 2024 14:52:04 -0400 Subject: [PATCH 06/10] https://github.com/NCAR/ccpp-physics/commit/363c82cf4a4759387909dafec6a14d0fff622e1c https://github.com/NCAR/ccpp-physics/commit/73d8f0c750d24df9dc62e6c362a5c7bc62acb216 --- gsmphys/module_sf_noahmp_glacier.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index ffcfb80f..45393d4c 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -2212,6 +2212,7 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in ! local integer :: iz real (kind=kind_phys) :: bdsnow !bulk density of snow (kg/m3) + real (kind=kind_phys),parameter :: mwd = 100. !< maximum water depth (mm) ! ---------------------------------------------------------------------- snoflow = 0.0 ponding1 = 0.0 @@ -2255,9 +2256,9 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in !to obtain equilibrium state of snow in glacier region - if(sneqv > 5000.0) then ! 5000 mm -> maximum water depth + if(sneqv > mwd .and. isnow /= 0) then ! 100 mm -> maximum water depth bdsnow = snice(0) / dzsnso(0) - snoflow = (sneqv - 5000.0) + snoflow = (sneqv - mwd) snice(0) = snice(0) - snoflow dzsnso(0) = dzsnso(0) - snoflow/bdsnow snoflow = snoflow / dt From be1b33a757c651ed478a3204c8c7f9c36ea585ee Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Fri, 24 May 2024 02:02:20 +0000 Subject: [PATCH 07/10] revert maximum water depth to 5 meters --- gsmphys/module_sf_noahmp_glacier.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index 45393d4c..2d486439 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -2212,7 +2212,7 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in ! local integer :: iz real (kind=kind_phys) :: bdsnow !bulk density of snow (kg/m3) - real (kind=kind_phys),parameter :: mwd = 100. !< maximum water depth (mm) + real (kind=kind_phys),parameter :: mwd = 5000. !< maximum water depth (mm) ! ---------------------------------------------------------------------- snoflow = 0.0 ponding1 = 0.0 From 7a61cf9c9fd29b9c245e9348fa1751177ebe39ee Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Fri, 24 May 2024 18:22:45 +0000 Subject: [PATCH 08/10] Update file GFS_typedefs.F90 --- GFS_layer/GFS_typedefs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index 3d7b2e70..71dc0262 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -2229,7 +2229,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: iopt_snf = 1 !rainfall & snowfall (1-jordan91; 2->bats; 3->noah) integer :: iopt_tbot = 2 !lower boundary of soil temperature (1->zero-flux; 2->noah) integer :: iopt_stc = 1 !snow/soil temperature time scheme (only layer 1) - integer :: iopt_gla = 1 !glacier option (1->phase change; 2->simple) + integer :: iopt_gla = 2 !glacier option (1->phase change; 2->simple) !--- tuning parameters for physical parameterizations logical :: ras = .false. !< flag for ras convection scheme From a4ea3ac0a29f4086bdc8d5cf73e7d3c9a27e83b4 Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Fri, 24 May 2024 18:25:51 +0000 Subject: [PATCH 09/10] Update file module_sf_noahmp_glacier.f90 --- gsmphys/module_sf_noahmp_glacier.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index 2d486439..748ac4d1 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -2256,7 +2256,7 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in !to obtain equilibrium state of snow in glacier region - if(sneqv > mwd .and. isnow /= 0) then ! 100 mm -> maximum water depth + if(sneqv > mwd .and. isnow /= 0) then bdsnow = snice(0) / dzsnso(0) snoflow = (sneqv - mwd) snice(0) = snice(0) - snoflow From 3c7a76a02751feab30b7b012a0e8abbdb8d4412d Mon Sep 17 00:00:00 2001 From: Kai-Yuan Cheng Date: Fri, 24 May 2024 18:32:00 +0000 Subject: [PATCH 10/10] Update file module_sf_noahmp_glacier.f90 --- gsmphys/module_sf_noahmp_glacier.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsmphys/module_sf_noahmp_glacier.f90 b/gsmphys/module_sf_noahmp_glacier.f90 index 748ac4d1..41a09d59 100644 --- a/gsmphys/module_sf_noahmp_glacier.f90 +++ b/gsmphys/module_sf_noahmp_glacier.f90 @@ -51,7 +51,7 @@ module noahmp_glacier_globals ! options for glacier treatment ! 1 -> include phase change of ice; 2 -> ice treatment more like original noah - integer :: opt_gla != 1 !(suggested 1) + integer :: opt_gla != 2 !(suggested 2) ! adjustable parameters for snow processes real (kind=kind_phys), parameter :: z0sno = 0.002 !snow surface roughness length (m) (0.002)