Winterwx

physics/GFS_MP_generic_post.F90

@@ -20,10 +20,11 @@ module GFS_MP_generic_post
 !> @{
       subroutine GFS_MP_generic_post_run(                                                                                 &
         im, levs, kdt, nrcm, nncl, ntcw, ntrac, imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_nssl,    &
-        imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, cpllnd, progsigma, con_g, rainmin, dtf, frain, rainc, &
+        imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, cpllnd, progsigma, con_g, rainmin, dtf, frain, rainc,    &
         rain1, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, snow, graupel, save_t, save_q, rain0, ice0, snow0,&
         graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp, totprcp, totice,   &
         totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, rain_cpl, rainc_cpl, snow_cpl, pwat,               &
+        frzr, frzrb, frozr, frozrb, tsnowp, tsnowpb, rhonewsn1, vrbliceden_noah, iopt_snf,                                & 
         drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv,                             &
         graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, dfi_radar_max_intervals,                  &
         dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d,dqdt_qmicro, lssav, num_dfi_radar,       &
@@ -37,7 +38,7 @@ subroutine GFS_MP_generic_post_run(
       integer, intent(in) :: im, levs, kdt, nrcm, nncl, ntcw, ntrac, num_dfi_radar, index_of_process_dfi_radar
       integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_mg, imp_physics_fer_hires
       integer, intent(in) :: imp_physics_nssl
-      logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, cpllnd, progsigma
+      logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, cpllnd, progsigma, vrbliceden_noah
       integer, intent(in) :: index_of_temperature,index_of_process_mp
 
       integer                                                :: dfi_radar_max_intervals
@@ -70,7 +71,7 @@ subroutine GFS_MP_generic_post_run(
       real(kind=kind_phys), dimension(:),      intent(inout) :: drain_cpl, dsnow_cpl
 
       ! Rainfall variables previous time step
-      integer, intent(in) :: lsm, lsm_ruc, lsm_noahmp
+      integer, intent(in) :: lsm, lsm_ruc, lsm_noahmp, iopt_snf
       real(kind=kind_phys), dimension(:),      intent(inout) :: raincprv
       real(kind=kind_phys), dimension(:),      intent(inout) :: rainncprv
       real(kind=kind_phys), dimension(:),      intent(inout) :: iceprv
@@ -81,6 +82,13 @@ subroutine GFS_MP_generic_post_run(
       real(kind=kind_phys), dimension(:),      intent(inout) :: diceprv
       real(kind=kind_phys), dimension(:),      intent(inout) :: dsnowprv
       real(kind=kind_phys), dimension(:),      intent(inout) :: dgraupelprv
+      real(kind=kind_phys), dimension(:),      intent(inout) :: frzr
+      real(kind=kind_phys), dimension(:),      intent(inout) :: frzrb
+      real(kind=kind_phys), dimension(:),      intent(inout) :: frozr
+      real(kind=kind_phys), dimension(:),      intent(inout) :: frozrb
+      real(kind=kind_phys), dimension(:),      intent(inout) :: tsnowp
+      real(kind=kind_phys), dimension(:),      intent(inout) :: tsnowpb
+      real(kind=kind_phys), dimension(:),      intent(inout) :: rhonewsn1
       real(kind=kind_phys), dimension(:,:),    intent(inout) :: dqdt_qmicro
       real(kind=kind_phys), dimension(:,:),    intent(inout) :: prevsq
       real(kind=kind_phys),                    intent(in)    :: dtp
@@ -101,6 +109,10 @@ subroutine GFS_MP_generic_post_run(
       real(kind=kind_phys) :: crain, csnow, onebg, tem, total_precip, tem1, tem2, ttend
       real(kind=kind_phys), dimension(im) :: domr, domzr, domip, doms, t850, work1
 
+      real :: snowrat,grauprat,icerat,curat,prcpncfr,prcpcufr
+      real :: rhonewsnow,rhoprcpice,rhonewgr,rhonewice
+      real(kind=kind_phys), parameter :: rhowater = 1000.0_kind_phys
+
       ! Initialize CCPP error handling variables
       errmsg = ''
       errflg = 0
@@ -111,6 +123,68 @@ subroutine GFS_MP_generic_post_run(
         rain(i) = rainc(i) + frain * rain1(i) ! time-step convective plus explicit
       enddo
 
+! compute surface snowfall, graupel/sleet, freezing rain and precip ice density
+      if (imp_physics == imp_physics_gfdl .or. imp_physics == imp_physics_thompson .or. imp_physics == imp_physics_nssl ) then
+         do i = 1, im
+            if (gt0(i,1) .le. 273) then
+               frzr(i) = frzr(i) + rain0(i)
+               frzrb(i) = frzrb(i) + rain0(i)
+            endif
+            tsnowp(i)  = tsnowp(i)  + snow0(i)
+            tsnowpb(i)  = tsnowpb(i)  + snow0(i)
+            frozr(i) = frozr(i) + graupel0(i)
+            frozrb(i) = frozrb(i) + graupel0(i)
+         enddo
+!Compute the variable precip ice density for specific LSM schemes and options
+         if ( lsm == lsm_ruc .or. lsm == lsm_noahmp .and. iopt_snf == 5 .or. vrbliceden_noah == .true.) then
+            snowrat = 0.
+            grauprat = 0.
+            icerat = 0.
+            prcpncfr = 0.
+            prcpcufr = 0.
+            curat = 0.
+            prcpncfr = 0.
+            prcpcufr = 0.
+            rhonewsnow = 0.
+            rhonewgr = 0.
+            rhonewice = 0.
+            rhoprcpice = 0.
+            do i = 1, im
+               rhonewsn1(i)= 200.
+               prcpncfr = rain1(i)*sr(i)
+               if(sr(i) > 0..and. gt0(i,1) < 273.) then
+                  prcpcufr = max(0.,rainc(i)*sr(i))
+               else
+                  if(gt0(i,1) < 273.) then
+                     prcpcufr = max(0.,rainc(i))
+                  else
+                     prcpcufr = 0.
+                  endif  ! gt0(i,1) < 273.
+               endif ! frzfrac > 0.
+!
+               if((prcpncfr + prcpcufr) > 0.) then
+! -- calculate snow, graupel and ice fractions in falling frozen precip
+                  snowrat=min(1.,max(0.,snow0(i)/(prcpncfr + prcpcufr)))
+                  grauprat=min(1.,max(0.,graupel0(i)/(prcpncfr + prcpcufr)))
+                  icerat=min(1.,max(0.,(prcpncfr-snow0(i)-graupel0(i)) &
+                     /(prcpncfr + prcpcufr)))
+                  curat=min(1.,max(0.,(prcpcufr/(prcpncfr + prcpcufr))))
+
+                  rhonewsnow=min(125.,1000.0/max(8.,(17.*tanh((276.65-gt0(i,1))*0.15))))
+                  rhonewgr=min(500.,rhowater/max(2.,(3.5*tanh((274.15-gt0(i,1))*0.3333))))
+                  rhonewice=rhonewsnow
+!--- compute density of "precip ice" from weighted contribution
+!                 of snow, graupel and ice fractions
+
+                  rhoprcpice = min(500.,max(58.8,(rhonewsnow*snowrat +  &
+                      rhonewgr*grauprat + rhonewice*icerat + rhonewgr*curat)))
+! from now on rhonewsn1 is the density of falling frozen precipitation
+                  rhonewsn1(i)=rhoprcpice
+               endif
+            enddo
+         endif
+      endif
+
 !> - If requested (e.g. Zhao-Carr MP scheme), call calpreciptype() to calculate dominant 
 !! precipitation type.
       ! DH* TODO - Fix wrong code in non-CCPP build (GFS_physics_driver)

physics/GFS_MP_generic_post.meta

@@ -278,6 +278,76 @@
   type = real
   kind = kind_phys
   intent = inout
+[frzr]
+  standard_name = lwe_thickness_of_sfc_freezing_rain_amount
+  long_name = accumulated surface freezing rain
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[frzrb]
+  standard_name = lwe_thickness_of_sfc_freezing_rain_amount_in_bucket
+  long_name = accumulated surface freezing rain in bucket
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[frozr]
+  standard_name = lwe_thickness_of_sfc_graupel_amount
+  long_name = accumulated surface graupel
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[frozrb]
+  standard_name = lwe_thickness_of_sfc_graupel_amount_in_bucket
+  long_name = accumulated surface graupel in bucket
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[tsnowp]
+  standard_name = lwe_thickness_of_sfc_snow_amount
+  long_name = accumulated surface snow
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[tsnowpb]
+  standard_name = lwe_thickness_of_sfc_snow_amount_in_bucket
+  long_name = accumulated surface snow in bucket
+  units = m
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[rhonewsn1]
+  standard_name = lwe_density_of_precip_ice
+  long_name = density of precipitation ice
+  units = kg m-3
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[vrbliceden_noah]
+  standard_name = flag_for_vrbl_prcp_ice_den
+  long_name = flag for variable precip ice density
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[iopt_snf]
+  standard_name = control_for_land_surface_scheme_precipitation_type_partition
+  long_name = choice for precipitation partition option (see noahmp module for definition)
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [save_t]
   standard_name = air_temperature_save
   long_name = air temperature before entering a physics scheme

physics/GFS_debug.F90

@@ -1285,7 +1285,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%qss_ice             ', Interstitial%qss_ice                 )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%qss_land            ', Interstitial%qss_land                )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%qss_water           ', Interstitial%qss_water               )
-                     call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%radar_reset         ', Interstitial%radar_reset             )
+                     call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%fullradar_diag      ', Interstitial%fullradar_diag          )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%raddt               ', Interstitial%raddt                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%raincd              ', Interstitial%raincd                  )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%raincs              ', Interstitial%raincs                  )

physics/lsm_noah.f

@@ -220,7 +220,8 @@ subroutine lsm_noah_run                                           &
      &       lheatstrg, isot, ivegsrc,                                  &
      &       bexppert, xlaipert, vegfpert,pertvegf,                     &  ! sfc perts, mgehne
      &       albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd,        &  
-     &       adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd,                &  
+     &       adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd, rhonewsn1,     &  
+     &       vrbliceden_noah,                                           &
 !  ---  in/outs:
      &       weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc,        &
      &       canopy, trans, tsurf, zorl,                                &
@@ -252,7 +253,7 @@ subroutine lsm_noah_run                                           &
      &                  -1.0_kind_phys, -2.0_kind_phys /
 
 !  ---  input:
-      integer, intent(in) :: im, km, isot, ivegsrc
+      integer, intent(in) :: im, km, isot, ivegsrc 
       real (kind=kind_phys), intent(in) :: grav, cp, hvap, rd, eps,     &
      &       epsm1, rvrdm1
       real (kind=kind_phys), intent(in) :: pertvegf
@@ -265,13 +266,13 @@ subroutine lsm_noah_run                                           &
      &       snoalb, sfalb, zf,                                         &
      &       bexppert, xlaipert, vegfpert,                              &
      &       albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd,        &
-     &       adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd
+     &       adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd, rhonewsn1
 
       real (kind=kind_phys),  intent(in) :: delt
 
       logical, dimension(:), intent(in) :: flag_iter, flag_guess, land
 
-      logical, intent(in) :: lheatstrg
+      logical, intent(in) :: lheatstrg, vrbliceden_noah
 
 !  ---  in/out:
       real (kind=kind_phys), dimension(:), intent(inout) :: weasd,      &
@@ -292,7 +293,7 @@ subroutine lsm_noah_run                                           &
 !  ---  locals:
       real (kind=kind_phys), dimension(im) :: rch, rho,                 &
      &       q0, qs1, theta1,       weasd_old, snwdph_old,              &
-     &       tprcp_old, srflag_old, tskin_old, canopy_old
+     &       tprcp_old, srflag_old, tskin_old, canopy_old              
 
       real (kind=kind_phys), dimension(km) :: et, sldpth, stsoil,       &
      &       smsoil, slsoil
@@ -309,8 +310,8 @@ subroutine lsm_noah_run                                           &
      &       smcdry, smcref, smcmax, sneqv, snoalb1d, snowh,            &
      &       snomlt, sncovr, soilw, soilm, ssoil, tsea, th2, tbot,      &
      &       xlai, zlvl, swdn, tem, z0, bexpp, xlaip, vegfp,            &
-     &       mv, sv, alphav, betav, vegftmp, cpinv, hvapi, elocp
-
+     &       mv, sv, alphav, betav, vegftmp, cpinv, hvapi, elocp,       &
+     &       rhonewsn 
       integer :: couple, ice, nsoil, nroot, slope, stype, vtype
       integer :: i, k, iflag
 !
@@ -326,7 +327,6 @@ subroutine lsm_noah_run                                           &
       errflg = 0
 
 !> - Save land-related prognostic fields for guess run.
-
       do i = 1, im
         if (land(i) .and. flag_guess(i)) then
           weasd_old(i)  = weasd(i)
@@ -335,7 +335,6 @@ subroutine lsm_noah_run                                           &
           canopy_old(i) = canopy(i)
           tprcp_old(i)  = tprcp(i)
           srflag_old(i) = srflag(i)
-
           do k = 1, km
             smc_old(i,k) = smc(i,k)
             stc_old(i,k) = stc(i,k)
@@ -361,7 +360,6 @@ subroutine lsm_noah_run                                           &
           sbsno(i)  = zero
           snowc(i)  = zero
           snohf(i)  = zero
-
 !> - Initialize variables wind, q, and rh at level 1.
 
           q0(i)   = max(q1(i), qmin)   !* q1=specific humidity at level 1 (kg/kg)
@@ -371,7 +369,6 @@ subroutine lsm_noah_run                                           &
           qs1(i) = fpvs( t1(i) )        !* qs1=sat. humidity at level 1 (kg/kg)
           qs1(i) = max(eps*qs1(i) / (prsl1(i)+epsm1*qs1(i)), qmin)
           q0 (i) = min(qs1(i), q0(i))
-
           do k = 1, km
             zsoil(i,k) = zsoil_noah(k)
           enddo
@@ -524,7 +521,8 @@ subroutine lsm_noah_run                                           &
 !> - Apply perturbation of soil type b parameter and leave area index.
           bexpp  = bexppert(i)                   ! sfc perts, mgehne
           xlaip  = xlaipert(i)                   ! sfc perts, mgehne
-
+!> - New snow depth variables 
+          rhonewsn = rhonewsn1(i)
 !> - Call Noah LSM gfssflx().
 
           call gfssflx                                                  & ! ccppdox: these is sflx in mpbl
@@ -533,6 +531,7 @@ subroutine lsm_noah_run                                           &
      &       swdn, solnet, lwdn, sfcems, sfcprs, sfctmp,                &
      &       sfcspd, prcp, q2, q2sat, dqsdt2, th2, ivegsrc,             &
      &       vtype, stype, slope, shdmin1d, alb, snoalb1d,              &
+     &       rhonewsn, vrbliceden_noah,                                 &
      &       bexpp, xlaip,                                              & ! sfc-perts, mgehne
      &       lheatstrg,                                                 &
 !  ---  input/outputs:

physics/lsm_noah.meta

@@ -486,6 +486,21 @@
   type = real
   kind = kind_phys
   intent = in
+[rhonewsn1]
+  standard_name = lwe_density_of_precip_ice
+  long_name = density of precipitation ice
+  units = kg m-3
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[vrbliceden_noah]
+  standard_name = flag_for_vrbl_prcp_ice_den
+  long_name = flag for variable precip ice density
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [weasd]
   standard_name = water_equivalent_accumulated_snow_depth_over_land
   long_name = water equiv of acc snow depth over land

physics/lsm_ruc.F90

@@ -331,7 +331,7 @@ subroutine lsm_ruc_run                                            & ! inputs
      &       prsl1, zf, wind, shdmin, shdmax,                           &
      &       srflag, sfalb_lnd_bck, snoalb,                             &
      &       isot, ivegsrc, fice, smcwlt2, smcref2,                     &
-     &       min_lakeice, min_seaice, oceanfrac,                        &
+     &       min_lakeice, min_seaice, oceanfrac, rhonewsn1,             &
      ! --- constants
      &       con_cp, con_rd, con_rv, con_g, con_pi, con_hvap,           &
      &       con_fvirt,                                                 &
@@ -418,7 +418,7 @@ subroutine lsm_ruc_run                                            & ! inputs
 
 !  ---  in
       real (kind=kind_phys), dimension(:), intent(in) ::                 &
-     &       rainnc, rainc, ice, snow, graupel
+     &       rainnc, rainc, ice, snow, graupel, rhonewsn1
 !  ---  in/out:
 !  --- on RUC levels
       real (kind=kind_phys), dimension(:,:), intent(inout) ::            &
@@ -494,7 +494,8 @@ subroutine lsm_ruc_run                                            & ! inputs
      &     soilt_lnd, tbot,                                             &
      &     xlai, swdn, z0_lnd, znt_lnd, rhosnfr, infiltr,               &
      &     precipfr, snfallac_lnd, acsn,                                &
-     &     qsfc_lnd, qsg_lnd, qvg_lnd, qcg_lnd, soilt1_lnd, chklowq
+     &     qsfc_lnd, qsg_lnd, qvg_lnd, qcg_lnd, soilt1_lnd, chklowq,    &
+     &     rhonewsn
      ! ice
       real (kind=kind_phys),dimension (im,1)        ::                  &
      &     albbck_ice, alb_ice, chs_ice, flhc_ice, flqc_ice,            &
@@ -845,6 +846,7 @@ subroutine lsm_ruc_run                                            & ! inputs
         rainncv(i,j)    = rhoh2o * rainnc(i)                   ! total time-step explicit precip 
         graupelncv(i,j) = rhoh2o * graupel(i)
         snowncv(i,j)    = rhoh2o * snow(i)
+        rhonewsn(i,j) = rhonewsn1(i)
         if (debug_print) then
         !-- diagnostics for a test point with known lat/lon
         if (abs(xlat_d(i)-testptlat).lt.2.5 .and.   &
@@ -1125,7 +1127,7 @@ subroutine lsm_ruc_run                                            & ! inputs
      &          conflx2(i,1,j), sfcprs(i,1,j), sfctmp(i,1,j), q2(i,1,j),     &
      &          qcatm(i,1,j), rho2(i,1,j), semis_bck(i,j), lwdn(i,j),        &
      &          swdn(i,j), solnet_lnd(i,j), sfcems_lnd(i,j), chklowq(i,j),   &
-     &          chs_lnd(i,j), flqc_lnd(i,j), flhc_lnd(i,j),                  &
+     &          chs_lnd(i,j), flqc_lnd(i,j), flhc_lnd(i,j), rhonewsn(i,j),   &
 !  ---  input/outputs:
      &          wet(i,j), cmc(i,j), shdfac(i,j), alb_lnd(i,j), znt_lnd(i,j), &
      &          z0_lnd(i,j), snoalb1d_lnd(i,j), albbck_lnd(i,j),             &
@@ -1399,7 +1401,7 @@ subroutine lsm_ruc_run                                            & ! inputs
      &          conflx2(i,1,j), sfcprs(i,1,j), sfctmp(i,1,j), q2(i,1,j),     &
      &          qcatm(i,1,j), rho2(i,1,j), semis_bck(i,j), lwdn(i,j),        &
      &          swdn(i,j), solnet_ice(i,j), sfcems_ice(i,j), chklowq(i,j),   &
-     &          chs_ice(i,j), flqc_ice(i,j), flhc_ice(i,j),                  &
+     &          chs_ice(i,j), flqc_ice(i,j), flhc_ice(i,j), rhonewsn(i,j),   &
 !  ---  input/outputs:
      &          wet_ice(i,j), cmc(i,j), shdfac(i,j), alb_ice(i,j),           &
      &          znt_ice(i,j), z0_ice(i,j), snoalb1d_ice(i,j),                &