fix the frozen precip density issue

HelinWei-NOAA · Apr 18, 2022 · 8261064 · 8261064
1 parent ab3c5e6
commit 8261064
Show file tree

Hide file tree

Showing 4 changed files with 38 additions and 2 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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