Update to Microphysics for HRRRv5

src/core_atmosphere/Registry.xml

@@ -631,6 +631,7 @@
 			<var name="tday_accum"/>
 			<var name="tyear_mean"/>
 			<var name="tyear_accum"/>
+			<var name="refl10cm"/>
 			<var name="refl10cm_max"/>
 			<var name="i_rainnc"/>
 			<var name="rainncv"/>
@@ -1086,6 +1087,7 @@
 			<var name="olrtoa"/>
 			<var name="rainc"/>
 			<var name="rainnc"/>
+			<var name="refl10cm"/>
 			<var name="refl10cm_max"/>
 			<var name="refl10cm_1km"/>
 			<var name="refl10cm_1km_max"/>
@@ -2237,6 +2239,16 @@
                      description="=0:activatevmixing of moisture species separately;= 1:activate mixing of total water in the mynn PBL scheme"
                      possible_values="0,1"/>
 
+		<nml_option name="config_thompson_aerosol_aware" type="logical" default_value="true" in_defaults="false"
+                     units="-"
+                     description="Logical flag to turn on/off prognostic cloud droplet and aerosol number concentrations"
+                     possible_values=".true. or .false."/>
+
+		<nml_option name="config_thompson_hail_aware" type="logical" default_value="false" in_defaults="false"
+                     units="-"
+                     description="Logical flag to turn on/off prognostic graupel number concentration and rime density"
+                     possible_values=".true. or .false."/>
+
                 <nml_option name="config_bucket_radt" type="real" default_value="1.0e9" in_defaults="false"
                      units="-"
                      description="threshold above which accumulated radiation diagnostics are reset"
@@ -2303,6 +2315,10 @@
                 <!-- ... PARAMETERIZATION OF CLOUD MICROPHYSICS:                                                        -->
                 <!-- ================================================================================================== -->
 
+		<var name="refl10cm" type="real" dimensions="nVertLevels nCells Time" units="dBZ"
+                     description="10 cm radar reflectivity"
+                     packages="mp_thompson_in;mp_wsm6_in"/>
+
                 <var name="refl10cm_max" type="real" dimensions="nCells Time" units="dBZ"
                      description="10 cm maximum radar reflectivity"
                      packages="mp_thompson_in;mp_wsm6_in"/>

src/core_atmosphere/physics/Makefile

@@ -113,7 +113,8 @@ mpas_atmphys_driver_lsm.o: \
 	mpas_atmphys_landuse.o \
 	mpas_atmphys_lsm_noahinit.o \
 	mpas_atmphys_lsm_rucinit.o \
-	mpas_atmphys_vars.o
+	mpas_atmphys_lsm_shared.o \
+	mpas_atmphys_vars.o 
 
 mpas_atmphys_driver_microphysics.o: \
 	mpas_atmphys_constants.o \

src/core_atmosphere/physics/mpas_atmphys_driver.F

@@ -280,15 +280,18 @@ subroutine physics_driver(domain,itimestep,xtime_s)
 !$OMP END PARALLEL DO
        call deallocate_lsm(config_frac_seaice,config_lsm_scheme)
 
-       call allocate_seaice
+       if(config_lsm_scheme .ne. 'sf_ruc') then
+          call allocate_seaice
 !$OMP PARALLEL DO
-       do thread=1,nThreads
-          call driver_seaice(block%configs,diag_physics,sfc_input, &
-                             cellSolveThreadStart(thread),cellSolveThreadEnd(thread))
-       enddo
+          do thread=1,nThreads
+             call driver_seaice(block%configs,diag_physics,sfc_input, &
+                                cellSolveThreadStart(thread),cellSolveThreadEnd(thread))
+          enddo
 !$OMP END PARALLEL DO
-       call deallocate_seaice,config_lsm_scheme
-    endif
+          call deallocate_seaice
+       endif ! sf_ruc
+    endif ! lsm off
+
 
     !call to pbl schemes:
     if(config_pbl_scheme .ne. 'off' .and. config_sfclayer_scheme .ne. 'off') then

src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F

@@ -12,6 +12,7 @@ module mpas_atmphys_driver_lsm
  use mpas_timer, only : mpas_timer_start, mpas_timer_stop
 
  use mpas_atmphys_constants
+ use mpas_atmphys_lsm_shared,only: correct_tsk_over_seaice
  use mpas_atmphys_landuse, only: isurban
  use mpas_atmphys_lsm_noahinit
  use mpas_atmphys_lsm_rucinit
@@ -187,6 +188,9 @@ subroutine allocate_lsm(config_frac_seaice,config_lsm_scheme)
  if(.not.allocated(xland_p)     ) allocate(xland_p(ims:ime,jms:jme)     )
  if(.not.allocated(z0_p)        ) allocate(z0_p(ims:ime,jms:jme)        )
  if(.not.allocated(znt_p)       ) allocate(znt_p(ims:ime,jms:jme)       )
+ if(.not.allocated(t2m_p)       ) allocate(t2m_p(ims:ime,jms:jme)       )
+ if(.not.allocated(th2m_p)      ) allocate(th2m_p(ims:ime,jms:jme)      )
+ if(.not.allocated(q2_p)        ) allocate(q2_p(ims:ime,jms:jme)        )
  if(.not.allocated(flxsnow_p)   ) allocate(flxsnow_p(ims:ime,jms:jme)   )
  if(.not.allocated(fvbsnow_p)   ) allocate(fvbsnow_p(ims:ime,jms:jme)   )
  if(.not.allocated(fbursnow_p)  ) allocate(fbursnow_p(ims:ime,jms:jme)  )
@@ -301,6 +305,9 @@ subroutine deallocate_lsm(config_frac_seaice,config_lsm_scheme)
  if(allocated(xland_p)     ) deallocate(xland_p     )
  if(allocated(z0_p)        ) deallocate(z0_p        )
  if(allocated(znt_p)       ) deallocate(znt_p       )
+ if(allocated(t2m_p)       ) deallocate(t2m_p       )
+ if(allocated(th2m_p)      ) deallocate(th2m_p      )
+ if(allocated(q2_p)        ) deallocate(q2_p        )
  if(allocated(flxsnow_p)   ) deallocate(flxsnow_p   )
  if(allocated(fvbsnow_p)   ) deallocate(fvbsnow_p   )
  if(allocated(fbursnow_p)  ) deallocate(fbursnow_p  )
@@ -367,6 +374,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  integer,intent(in):: its,ite
 
 !local pointers:
+ logical, pointer              :: config_frac_seaice
  character(len=StrKIND),pointer:: config_microp_scheme,    &
                                   config_convection_scheme, &
                                   config_lsm_scheme
@@ -398,6 +406,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_config(configs,'config_convection_scheme',config_convection_scheme)
  call mpas_pool_get_config(configs,'config_microp_scheme'    ,config_microp_scheme    )
  call mpas_pool_get_config(configs,'config_lsm_scheme'       ,config_lsm_scheme       )
+ call mpas_pool_get_config(configs,'config_frac_seaice'      ,config_frac_seaice      )
 
  call mpas_pool_get_array(diag_physics,'acsnom'           ,acsnom           )
  call mpas_pool_get_array(diag_physics,'acsnow'           ,acsnow           )
@@ -559,6 +568,9 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
     udrunoff_p(i,j)   = udrunoff(i)
     z0_p(i,j)         = z0(i)
     znt_p(i,j)        = znt(i)
+    t2m_p(i,j)        = t2m(i)
+    th2m_p(i,j)       = th2m(i)
+    q2_p(i,j)         = q2(i)
 
     isltyp_p(i,j)     = isltyp(i)
     ivgtyp_p(i,j)     = ivgtyp(i)
@@ -728,6 +740,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
 
  call mpas_pool_get_config(configs,'config_microp_scheme',config_microp_scheme)
  call mpas_pool_get_config(configs,'config_lsm_scheme',config_lsm_scheme)
+ call mpas_pool_get_config(configs,'config_frac_seaice',config_frac_seaice)
 
  call mpas_pool_get_array(diag_physics,'acsnom'    ,acsnom    )
  call mpas_pool_get_array(diag_physics,'acsnow'    ,acsnow    )
@@ -856,6 +869,9 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
     udrunoff(i)   = udrunoff_p(i,j)
     z0(i)         = z0_p(i,j)
     znt(i)        = znt_p(i,j)
+    t2m(i)        = t2m_p(i,j)
+    th2m(i)       = th2m_p(i,j)
+    q2(i)         = q2_p(i,j)
 
     snoalb(i)     = snoalb_p(i,j)
     sfc_albbck(i) = sfc_albbck_p(i,j)

src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F

@@ -147,6 +147,9 @@ subroutine allocate_microphysics(configs)
        if(.not.allocated(reice_p)   ) allocate(reice_p(ims:ime,kms:kme,jms:jme)  )
        if(.not.allocated(resnow_p)  ) allocate(resnow_p(ims:ime,kms:kme,jms:jme) )
 
+       ! 3D reflectivity
+       if(.not.allocated(refl10cm_p)  ) allocate(refl10cm_p(ims:ime,kms:kme,jms:jme) )
+
     microp2_select: select case(microp_scheme)
 
        case("mp_thompson")
@@ -224,6 +227,8 @@ subroutine deallocate_microphysics(configs)
        if(allocated(reice_p)   ) deallocate(reice_p   )
        if(allocated(resnow_p)  ) deallocate(resnow_p  )
 
+       if(allocated(refl10cm_p)  ) deallocate(refl10cm_p  )
+
     microp2_select: select case(microp_scheme)
 
        case("mp_thompson")
@@ -250,7 +255,7 @@ subroutine deallocate_microphysics(configs)
  end subroutine deallocate_microphysics
 
 !=================================================================================================================
- subroutine microphysics_init(dminfo,configs,mesh,sfc_input,diag_physics,state)
+ subroutine microphysics_init(dminfo,configs,mesh,sfc_input,diag,diag_physics,state)
 !=================================================================================================================
 
 !input arguments:
@@ -262,10 +267,12 @@ subroutine microphysics_init(dminfo,configs,mesh,sfc_input,diag_physics,state)
 !inout arguments:
  type(mpas_pool_type),intent(inout):: diag_physics
  type(mpas_pool_type),intent(inout):: state
+ type(mpas_pool_type),intent(in):: diag
 
 !local pointer:
  character(len=StrKIND),pointer:: microp_scheme
-
+ logical,pointer:: thompson_aerosol_aware,thompson_hail_aware
+
 !CCPP-compliant flags:
  character(len=StrKIND):: errmsg
  integer:: errflg
@@ -277,12 +284,14 @@ subroutine microphysics_init(dminfo,configs,mesh,sfc_input,diag_physics,state)
  errflg = 0
 
  call mpas_pool_get_config(configs,'config_microp_scheme',microp_scheme)
-
+ call mpas_pool_get_config(configs,'config_thompson_aerosol_aware',thompson_aerosol_aware)
+ call mpas_pool_get_config(configs,'config_thompson_hail_aware',thompson_hail_aware)
+
   microp_select: select case(microp_scheme)
 
      case("mp_thompson")
-        call thompson_init(l_mp_tables)
-        call init_thompson_clouddroplets_forMPAS(mesh,sfc_input,diag_physics,state)
+        call thompson_init(l_mp_tables, thompson_hail_aware, thompson_aerosol_aware)
+        call init_thompson_clouddroplets_forMPAS(mesh,sfc_input,diag,diag_physics,state,configs)
 
      case("mp_wsm6")
         call mp_wsm6_init(den0=rho_a,denr=rho_r,dens=rho_s,cl=cliq,cpv=cpv, &
@@ -316,7 +325,8 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
 
 !local pointers:
  character(len=StrKIND),pointer:: microp_scheme
-
+ logical,pointer:: thompson_aerosol_aware,thompson_hail_aware
+
 !local variables and arrays:
  integer:: istep
 
@@ -333,7 +343,9 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
  errflg = 0
 
  call mpas_pool_get_config(configs,'config_microp_scheme',microp_scheme)
-
+ call mpas_pool_get_config(configs,'config_thompson_aerosol_aware',thompson_aerosol_aware)
+ call mpas_pool_get_config(configs,'config_thompson_hail_aware',thompson_hail_aware)
+
 !... allocation of microphysics arrays:
 !$OMP MASTER
  call allocate_microphysics(configs)
@@ -369,7 +381,9 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
        istep = 1
        call mpas_timer_start('Thompson')
        do while (istep .le. n_microp)
-       call mp_gt_driver( &
+
+       if (thompson_aerosol_aware) then
+          call thompson_3d_to_1d_driver( &
                   th        = th_p        , qv         = qv_p         , qc         = qc_p         , &
                   qr        = qr_p        , qi         = qi_p         , qs         = qs_p         , &
                   qg        = qg_p        , ni         = ni_p         , nr         = nr_p         , &
@@ -381,11 +395,30 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
                   sr        = sr_p        , rainprod   = rainprod_p   , evapprod   = evapprod_p   , &
                   re_cloud  = recloud_p   , re_ice     = reice_p      , re_snow    = resnow_p     , &
                   has_reqc  = has_reqc    , has_reqi   = has_reqi     , has_reqs   = has_reqs     , &
-                  ntc       = ntc_p       , muc        = muc_p        ,                             &
+                  ntc       = ntc_p       , muc        = muc_p        , refl_10cm  = refl10cm_p   , &
                   ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde           , &
                   ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme           , &
                   its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte             &
-                           )
+                  )
+       else
+          call thompson_3d_to_1d_driver( &
+                  th        = th_p        , qv         = qv_p         , qc         = qc_p         , &
+                  qr        = qr_p        , qi         = qi_p         , qs         = qs_p         , &
+                  qg        = qg_p        , ni         = ni_p         , nr         = nr_p         , &
+!!!                  nc        = nc_p        , nwfa       = nwfa_p       , nifa       = nifa_p       , &
+                  pii       = pi_p        , p          = pres_p       , dz         = dz_p         , &
+                  w         = w_p         , dt_in      = dt_microp    , itimestep  = itimestep    , &
+                  rainnc    = rainnc_p    , rainncv    = rainncv_p    , snownc     = snownc_p     , &
+                  snowncv   = snowncv_p   , graupelnc  = graupelnc_p  , graupelncv = graupelncv_p , &
+                  sr        = sr_p        , rainprod   = rainprod_p   , evapprod   = evapprod_p   , &
+                  re_cloud  = recloud_p   , re_ice     = reice_p      , re_snow    = resnow_p     , &
+                  has_reqc  = has_reqc    , has_reqi   = has_reqi     , has_reqs   = has_reqs     , &
+                  ntc       = ntc_p       , muc        = muc_p        , refl_10cm  = refl10cm_p   , &
+                  ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde           , &
+                  ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme           , &
+                  its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte             &
+                  )
+       endif
        istep = istep + 1
        enddo
        call mpas_timer_stop('Thompson')
@@ -627,7 +660,7 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite)
 
 !local pointers:
  character(len=StrKIND),pointer:: microp_scheme
- real(kind=RKIND),dimension(:),pointer:: refl10cm_max,refl10cm_1km,refl10cm_1km_max
+ real(kind=RKIND),dimension(:),pointer:: refl10cm,refl10cm_max,refl10cm_1km,refl10cm_1km_max
 
 !local variables and arrays:
  integer:: i,j,k,kp
@@ -637,7 +670,8 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite)
 !-----------------------------------------------------------------------------------------------------------------
 
  call mpas_pool_get_config(configs,'config_microp_scheme',microp_scheme)
-
+
+ call mpas_pool_get_array(diag_physics,'refl10cm',refl10cm)
  call mpas_pool_get_array(diag_physics,'refl10cm_max',refl10cm_max)
  call mpas_pool_get_array(diag_physics,'refl10cm_1km',refl10cm_1km)
  call mpas_pool_get_array(diag_physics,'refl10cm_1km_max',refl10cm_1km_max)
@@ -718,11 +752,12 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite)
              qs1d(k)  = qs_p(i,k,j)
              qg1d(k)  = qg_p(i,k,j)
              nr1d(k)  = nr_p(i,k,j)
-             dBZ1d(k) = -35._RKIND
+             dBZ1d(k) = refl10cm_p(i,k,j)
              zp(k) = z_p(i,k,j) - z_p(i,1,j)+0.5*dz_p(i,1,j) ! height AGL
           enddo
-
-          call calc_refl10cm(qv1d,qc1d,qr1d,nr1d,qs1d,qg1d,t1d,p1d,dBZ1d,kts,kte,i,j)
+
+          ! AAJ This function is called directly from the microphysics
+          ! call calc_refl10cm(qv1d,qc1d,qr1d,nr1d,qs1d,qg1d,t1d,p1d,dBZ1d,kts,kte,i,j)
 
           kp = 1
           do k = kts,kte

src/core_atmosphere/physics/mpas_atmphys_init.F

@@ -363,7 +363,7 @@ subroutine physics_init(dminfo,clock,configs,mesh,diag,tend,state,time_lev,diag_
 
 !initialization of cloud microphysics processes:
  if(config_microp_scheme .ne. 'off') &
-    call microphysics_init(dminfo,configs,mesh,sfc_input,diag_physics,state)
+    call microphysics_init(dminfo,configs,mesh,sfc_input,diag,diag_physics,state)
 
 !initialization of PBL processes:
  if(config_pbl_scheme .ne. 'off') call init_pbl(configs)