From 826106487303a1fea723094f370226a6ead253c8 Mon Sep 17 00:00:00 2001 From: helin wei Date: Mon, 18 Apr 2022 13:00:32 +0000 Subject: [PATCH] fix the frozen precip density issue --- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 5 +++-- gfsphysics/GFS_layer/GFS_typedefs.F90 | 4 ++++ gfsphysics/physics/sfc_drv.f | 11 +++++++++++ gfsphysics/physics/sflx.f | 20 ++++++++++++++++++++ 4 files changed, 38 insertions(+), 2 deletions(-) diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 2bacf0edf..aafee84ea 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1822,6 +1822,7 @@ subroutine GFS_physics_driver & Radtend%sfalb, flag_iter, flag_guess, Model%lheatstrg, & Model%isot, Model%ivegsrc, & bexp1d, xlai1d, vegf1d, Model%pertvegf, & + Sfcprop%dgraupelprv, Sfcprop%diceprv, & ! --- input/output: weasd3(:,1), snowd3(:,1), tsfc3(:,1), tprcp3(:,1), & Sfcprop%srflag, smsoil, stsoil, slsoil, Sfcprop%canopy, & @@ -5246,7 +5247,7 @@ subroutine GFS_physics_driver & ! --- get the amount of different precip type for Noah MP ! --- convert from m/dtp to mm/s - if (Model%lsm==Model%lsm_noahmp) then + if (Model%lsm==Model%lsm_noahmp .or. Model%lsm==Model%lsm_noah) then if (Model%imp_physics == Model%imp_physics_mg .or. & Model%imp_physics == Model%imp_physics_gfdl) then !GJF: Should all precipitation rates have the same denominator below? @@ -5266,7 +5267,7 @@ subroutine GFS_physics_driver & Sfcprop%dgraupelprv(:) = 0.0 Sfcprop%diceprv(:) = 0.0 endif - end if ! if (Model%lsm == Model%lsm_noahmp) + end if ! if (Model%lsm == Model%lsm_noahmp .or. Model%lsm==Model%lsm_noah) if (Model%cal_pre) then ! hchuang: add dominant precipitation type algorithm ! diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 499be39b0..efbf64dcc 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -2308,6 +2308,10 @@ subroutine sfcprop_create (Sfcprop, IM, Model) Sfcprop%smoiseq = clear_val Sfcprop%zsnsoxy = clear_val + endif + + if (Model%lsm == Model%lsm_noahmp .or. Model%lsm == Model%lsm_noah) then + allocate(Sfcprop%draincprv (IM)) allocate(Sfcprop%drainncprv (IM)) allocate(Sfcprop%diceprv (IM)) diff --git a/gfsphysics/physics/sfc_drv.f b/gfsphysics/physics/sfc_drv.f index e5626362f..86143cb27 100644 --- a/gfsphysics/physics/sfc_drv.f +++ b/gfsphysics/physics/sfc_drv.f @@ -145,6 +145,7 @@ subroutine sfc_drv & & shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, & & lheatstrg, isot, ivegsrc, & & bexppert, xlaipert, vegfpert,pertvegf, & ! sfc perts, mgehne + & graupel_mp, ice_mp, & ! --- in/outs: & weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & & canopy, trans, tsurf, zorl, & @@ -190,6 +191,9 @@ subroutine sfc_drv & & snoalb, sfalb, zf, & bexppert, xlaipert, vegfpert + real (kind=kind_phys), dimension(im), intent(in) :: & + & graupel_mp, ice_mp & + real (kind=kind_phys), intent(in) :: delt logical, dimension(im), intent(in) :: flag_iter, flag_guess, land @@ -231,6 +235,8 @@ subroutine sfc_drv & & xlai, zlvl, swdn, tem, z0, bexpp, xlaip, vegfp, & & mv,sv,alphav,betav,vegftmp + real (kind=kind_phys) :: graupel_prcp, ice_prcp + integer :: couple, ice, nsoil, nroot, slope, stype, vtype integer :: i, k, iflag @@ -320,6 +326,10 @@ subroutine sfc_drv & ! ffrozp = 0.0 ! endif ffrozp = srflag(i) + + graupel_prcp = graupel_mp(i) ! grpl part of MP precip [mm/s] + ice_prcp = ice_mp(i) ! ice part of MP precip [mm/s] + ice = 0 zlvl = zf(i) @@ -451,6 +461,7 @@ subroutine sfc_drv & & vtype, stype, slope, shdmin1d, alb, snoalb1d, & & bexpp, xlaip, & ! sfc-perts, mgehne & lheatstrg, & + & graupel_prcp, ice_prcp, & ! --- input/outputs: & tbot, cmc, tsea, stsoil, smsoil, slsoil, sneqv, chx, cmx, & & z0, & diff --git a/gfsphysics/physics/sflx.f b/gfsphysics/physics/sflx.f index 24cd8209e..9798b526c 100644 --- a/gfsphysics/physics/sflx.f +++ b/gfsphysics/physics/sflx.f @@ -8,6 +8,7 @@ subroutine sflx & & vegtyp, soiltyp, slopetyp, shdmin, alb, snoalb, & & bexpp, xlaip, & ! sfc-perts, mgehne & lheatstrg, & + & graupel_prcp, ice_prcp, & ! --- input/outputs: & tbot, cmc, t1, stc, smc, sh2o, sneqv, ch, cm,z0, & ! --- outputs: @@ -206,6 +207,8 @@ subroutine sflx & logical, intent(in) :: lheatstrg + real (kind=kind_phys), intent(in) :: ice_prcp, graupel_prcp + ! --- input/outputs: real (kind=kind_phys), intent(inout) :: tbot, cmc, t1, sneqv, & & stc(nsoil), smc(nsoil), sh2o(nsoil), ch, cm @@ -2707,6 +2710,7 @@ subroutine snow_new ! --- locals: real(kind=kind_phys) :: dsnew, snowhc, hnewc, newsnc, tempc + real(kind=kind_phys) :: dgnew, dinew, dfnew ! !===> ... begin here @@ -2728,6 +2732,22 @@ subroutine snow_new dsnew = 0.05 + 0.0017*(tempc + 15.0)**1.5 endif + dgnew = 1000.0 / max(2.,(3.5*tanh((274.15-sfctmp)*0.3333))) + dgnew = min(500.0, dgnew) ! kg/m3 (from RUC model) + dgnew = dgnew / 1000.0 ! convert units to [m-liq/m-snow] +! dinew = 250.0 / 1000.0 ! assume cloud ice is 250 kg/m3 + dfnew = 500.0 / 1000.0 ! assume freezing rain is 500 kg/m3 + + if (snowng) then + dsnew = (dsnew*(prcp * ffrozp - graupel_prcp ) + + & dgnew*graupel_prcp)/(prcp * ffrozp) +! dsnew = (dsnew*(prcp * ffrozp - graupel_prcp - ice_prcp) + +! & dgnew*graupel_prcp + dinew*ice_prcp)/(prcp * ffrozp) + endif + if (frzgra) then + dsnew = dfnew + endif + ! --- ... adjustment of snow density depending on new snowfall hnewc = newsnc / dsnew