MYNN-EDMF update

CMakeLists.txt

@@ -6,7 +6,7 @@ project(ccpp_physics
 
 #------------------------------------------------------------------------------
 set(PACKAGE "ccpp-physics")
-set(AUTHORS "Grant Firl" "Dom Heinzeller" "Man Zhang" "Mike Kavulich" "Chunxi Zhang")
+set(AUTHORS "Grant Firl" "Dustin Swales" "Man Zhang" "Mike Kavulich" )
 
 #------------------------------------------------------------------------------
 # Set OpenMP flags for C/C++/Fortran

physics/mynnedmf_wrapper.meta

@@ -125,13 +125,6 @@
   dimensions = ()
   type = logical
   intent = in
-[lheatstrg]
-  standard_name = flag_for_canopy_heat_storage_in_land_surface_scheme
-  long_name = flag for canopy heat storage parameterization
-  units = flag
-  dimensions = ()
-  type = logical
-  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP
@@ -303,14 +296,22 @@
   type = real
   kind = kind_phys
   intent = inout
-[qgrs_ice_cloud]
+[qgrs_ice]
   standard_name = cloud_ice_mixing_ratio
   long_name = ratio of mass of ice water to mass of dry air plus vapor (without condensates)
   units = kg kg-1
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
   intent = inout
+[qgrs_snow]
+  standard_name = snow_mixing_ratio
+  long_name = ratio of mass of snow water to mass of dry air plus vapor (without condensates)
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
 [qgrs_cloud_droplet_num_conc]
   standard_name = mass_number_concentration_of_cloud_liquid_water_particles_in_air
   long_name = number concentration of cloud droplets (liquid)
@@ -367,6 +368,14 @@
   type = real
   kind = kind_phys
   intent = in
+[prsi]
+  standard_name = air_pressure_at_interface
+  long_name = air pressure at model layer interfaces
+  units = Pa
+  dimensions = (horizontal_loop_extent,vertical_interface_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
 [exner]
   standard_name = dimensionless_exner_function
   long_name = Exner function at layers
@@ -1017,14 +1026,22 @@
   type = real
   kind = kind_phys
   intent = inout
-[dqdt_ice_cloud]
+[dqdt_ice]
   standard_name = process_split_cumulative_tendency_of_cloud_ice_mixing_ratio
   long_name = cloud condensed water mixing ratio tendency due to model physics
   units = kg kg-1 s-1
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
   intent = inout
+[dqdt_snow]
+  standard_name = process_split_cumulative_tendency_of_snow_mixing_ratio
+  long_name = ratio of mass of snow water tendency to mass of dry air plus vapor (without condensates) due to model physics
+  units = kg kg-1 s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
 [dqdt_ozone]
   standard_name = process_split_cumulative_tendency_of_ozone_mixing_ratio
   long_name = ozone mixing ratio tendency due to model physics
@@ -1151,6 +1168,13 @@
   dimensions = ()
   type = integer
   intent = in
+[ntsw]
+  standard_name = index_of_snow_mixing_ratio_in_tracer_concentration_array
+  long_name = tracer index for snow water
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [ntlnc]
   standard_name = index_of_mass_number_concentration_of_cloud_droplets_in_tracer_concentration_array
   long_name = tracer index for liquid number concentration
@@ -1210,12 +1234,12 @@
   type = real
   kind = kind_phys
   intent = in
-[bl_mynn_tkebudget]
+[tke_budget]
   standard_name = control_for_tke_budget_output
   long_name = flag for activating TKE budget
   units = flag
   dimensions = ()
-  type = logical
+  type = integer
   intent = in
 [bl_mynn_tkeadvect]
   standard_name = flag_for_tke_advection
@@ -1329,6 +1353,13 @@
   dimensions = ()
   type = integer
   intent = in
+[imp_physics_fa]
+  standard_name = identifier_for_fer_hires_microphysics_scheme
+  long_name = choice of Ferrier-Aligo microphysics scheme
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = in
 [imp_physics_nssl]
   standard_name = identifier_for_nssl_microphysics_scheme
   long_name = choice of NSSL 2-moment microphysics scheme
@@ -1359,7 +1390,7 @@
   type = real
   kind = kind_phys
   intent = inout
-[rrfs_smoke]
+[rrfs_sd]
   standard_name = do_smoke_coupling
   long_name = flag controlling rrfs_smoke collection (default off)
   units = flag
@@ -1373,7 +1404,7 @@
   dimensions = ()
   type = logical
   intent = in
-[fire_turb]
+[enh_mix]
   standard_name = do_planetary_boundary_layer_fire_enhancement
   long_name = flag for rrfs smoke mynn enh vermix
   units = flag

physics/noahmpdrv.F90

@@ -77,13 +77,6 @@ subroutine noahmpdrv_init(lsm, lsm_noahmp, me, isot, ivegsrc, &
           return
         end if
 
-        if (.not. do_mynnsfclay .and. do_mynnedmf) then
-          errmsg = 'Problem : do_mynnsfclay = .false.' // &
-                   'but mynnpbl is .true.. Exiting ...'
-          errflg = 1
-          return
-        end if
-
         if ( do_mynnsfclay .and. .not. do_mynnedmf) then
           errmsg = 'Problem : do_mynnsfclay = .true.' // &
                    'but mynnpbl is .false.. Exiting ...'

physics/sgscloud_radpre.F90

@@ -35,16 +35,17 @@ module sgscloud_radpre
 !!
 !>\section sgscloud_radpre_mod  SGS Cloud Scheme Pre General Algorithm
       subroutine sgscloud_radpre_run(    &
-           im,dt,levs,                   &
+           im,dt,fhswr,levs,             &
            flag_init,flag_restart,       &
            con_g, con_pi, eps, epsm1,    &
            r_v, cpv, rcp,                &
            xlv, xlf, cp,                 &
            do_mynnedmf,                  &
            qc, qi, qv, T3D, P3D, exner,  &
            qr, qs, qg,                   &
-           qci_conv,ud_mf,               &
+           qci_conv,qlc,qli,ud_mf,       &
            imfdeepcnv, imfdeepcnv_gf,    &
+           imfdeepcnv_sas,               &
            qc_save, qi_save, qs_save,    &
            qc_bl,qi_bl,cldfra_bl,        &
            delp,clouds1,clouds2,clouds3, &
@@ -53,6 +54,7 @@ subroutine sgscloud_radpre_run(    &
            nlay, plyr, xlat, dz,de_lgth, &
            cldsa,mtopa,mbota,            &
            imp_physics, imp_physics_gfdl,&
+           imp_physics_fa,               &
            iovr,                         &
            errmsg, errflg                )
 
@@ -67,17 +69,18 @@ subroutine sgscloud_radpre_run(    &
       real(kind=kind_phys), intent(in) :: con_g, con_pi, eps, epsm1
       real(kind=kind_phys), intent(in) :: r_v, cpv, rcp
       real(kind=kind_phys), intent(in) :: xlv, xlf, cp
-      real(kind=kind_phys), intent(in) :: dt
+      real(kind=kind_phys), intent(in) :: dt,fhswr
       real :: xls, xlvcp, xlscp !derived below
       real(kind=kind_phys)             :: gfac
       integer,             intent(in)  :: im, levs, imfdeepcnv, imfdeepcnv_gf, &
-           &               nlay, imp_physics, imp_physics_gfdl
+           &  nlay, imfdeepcnv_sas, imp_physics, imp_physics_gfdl, imp_physics_fa
       logical,             intent(in)  :: flag_init, flag_restart, do_mynnedmf
 
       real(kind=kind_phys), dimension(:,:), intent(inout) :: qc, qi
       real(kind=kind_phys), dimension(:,:), intent(inout) :: qr, qs, qg
-      ! qci_conv only allocated if GF is used
+      ! note: qci_conv only allocated if GF is used
       real(kind=kind_phys), dimension(:,:), intent(inout) :: qci_conv
+      real(kind=kind_phys), dimension(:,:), intent(inout) :: qlc, qli !for SAS
       real(kind=kind_phys), dimension(:,:), intent(in)    :: ud_mf
       real(kind=kind_phys), dimension(:,:), intent(in)    :: T3D,delp
       real(kind=kind_phys), dimension(:,:), intent(in)    :: qv,P3D,exner
@@ -112,7 +115,8 @@ subroutine sgscloud_radpre_run(    &
       real :: rhgrid,h2oliq,qsat,tem1,tem2,clwt,es,onemrh,value
 
       !Chaboureau and Bechtold (2002 and 2005)
-      real :: a, f, sigq, qmq, qt, xl, tlk, th, thl, rsl, cpm, cb_cf
+      real :: a, f, sigq, qmq, qt, xl, th, thl, rsl, cpm, cb_cf
+      real(kind=kind_phys) :: tlk
 
       !Option to convective cloud fraction
       integer, parameter :: conv_cf_opt = 0  !0: C-B, 1: X-R
@@ -188,7 +192,7 @@ subroutine sgscloud_radpre_run(    &
               !endif
 
               if (qc(i,k) < 1.e-6 .and. cldfra_bl(i,k)>0.001) then
-                qc(i,k) = qc_bl(i,k)*cldfra_bl(i,k)
+                qc(i,k) = qc_bl(i,k)
 
                 !eff radius cloud water (microns) from Miles et al. (2007)
                 if (nint(slmsk(i)) == 1) then !land
@@ -206,8 +210,8 @@ subroutine sgscloud_radpre_run(    &
               !~700 mb and decrease snow to zero by ~300 mb 
               snow_frac = min(0.5, max((p3d(i,k)-30000.0),0.0)/140000.0)
               ice_frac  = 1.0 - snow_frac
-              if (qi(i,k) < 1.e-8 .and. cldfra_bl(i,k)>0.001) then
-                qi(i,k) = ice_frac*qi_bl(i,k)*cldfra_bl(i,k)
+              if (qi(i,k) < 1.e-9 .and. cldfra_bl(i,k)>0.001) then
+                qi(i,k) = ice_frac*qi_bl(i,k)
 
                 !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b)
                 if(qi(i,k)>1.E-8)clouds5(i,k)=max(173.45 + 2.14*Tc, 20.)
@@ -219,8 +223,8 @@ subroutine sgscloud_radpre_run(    &
                 clouds4(i,k) = max(0.0, qi(i,k) * gfac * delp(i,k))
               endif
 
-              if (qs(i,k) < 1.e-8 .and. cldfra_bl(i,k)>0.001) then
-                qs(i,k) = snow_frac*qi_bl(i,k)*cldfra_bl(i,k)
+              if (qs(i,k) < 1.e-9 .and. cldfra_bl(i,k)>0.001) then
+                qs(i,k) = snow_frac*qi_bl(i,k)
 
                 !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b)
                 if(qs(i,k)>1.E-8)clouds9(i,k)=max(2.*(173.45 + 2.14*Tc), 50.)
@@ -270,7 +274,6 @@ subroutine sgscloud_radpre_run(    &
         if (imfdeepcnv == imfdeepcnv_gf) then
           do k = 1, levs
             do i = 1, im
-              !if ( qci_conv(i,k) > 0. .AND. (qi(i,k) < 1E-7 .AND. qc(i,k) < 1E-7 ) ) then
               if ( qci_conv(i,k) > 0. ) then
                 Tk = T3D(i,k)
                 Tc = Tk - 273.15
@@ -321,10 +324,15 @@ subroutine sgscloud_radpre_run(    &
                    sigq = SQRT(sigq**2 + 1e-10)   ! combined conv + background components
                    qmq  = a * (qt - qsat)         ! saturation deficit/excess;
                                                   !   the numerator of Q1
-                   cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.01),0.99)
+                   cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.0),0.99)
+                   if (qci_conv(i,k) .lt. 1e-9) cb_cf = 0.0
                    if (do_mynnedmf .and. qmq .ge. 0.0) then
                       ! leverage C-B stratus clouds from MYNN in saturated conditions
-                      clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
+                      if (cb_cf .gt. 0.0) then
+                         clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
+                      else
+                         !default to MYNN clouds - already specified
+                      endif
                    else                           ! unsaturated
                       clouds1(i,k) = cb_cf
                    endif
@@ -354,7 +362,101 @@ subroutine sgscloud_radpre_run(    &
               endif ! qci_conv
             enddo
           enddo
-        endif ! imfdeepcnv_gf
+
+        elseif (imfdeepcnv == imfdeepcnv_sas) then
+
+          do k = 1, levs
+            do i = 1, im
+              h2oliq  = qlc(i,k)+qli(i,k)
+              if ( h2oliq > 0. ) then
+                Tk = T3D(i,k)
+                Tc = Tk - 273.15
+
+                !Partition the convective clouds into water & frozen species
+                liqfrac = min(1., max(0., (Tk-244.)/29.))
+
+                qc(i,k) = qc(i,k)+qlc(i,k)
+                !split ice & snow 50-50%
+                qi(i,k) = qi(i,k)+0.5*qli(i,k)
+                qs(i,k) = qs(i,k)+0.5*qli(i,k)
+
+                !eff radius cloud water (microns)
+                if (nint(slmsk(i)) == 1) then !land
+                  if(qc(i,k)>1.E-8)clouds3(i,k)=5.4
+                else
+                  !from Miles et al.
+                  if(qc(i,k)>1.E-8)clouds3(i,k)=9.6
+                endif
+                !from Mishra et al. (2014, JGR Atmos), assume R_sno = 2*R_ice
+                if(qi(i,k)>1.e-8)clouds5(i,k)=max(     173.45 + 2.14*Tc , 20.)
+                if(qs(i,k)>1.e-8)clouds9(i,k)=max(2.0*(173.45 + 2.14*Tc), 50.)
+
+                if ( conv_cf_opt .eq. 0 ) then
+                   !print *,'Chab-Bechtold cloud fraction used'
+                   !Alternatively, use Chaboureau-Bechtold (CB) convective component
+                   !Based on both CB2002 and CB2005.
+                   xl  = xlv*liqfrac + xls*(1.-liqfrac)  ! blended heat capacity
+                   tlk = t3d(i,k) - xlvcp/exner(i,k)*qc(i,k) &
+                       &          - xlscp/exner(i,k)*qi(i,k)! liquid temp
+                   ! get saturation water vapor mixing ratio at tl and p
+                   es  = min( p3d(i,k), fpvs( tlk ) )   ! fpvs and prsl in pa
+                   qsat= max( QMIN, eps*es / (p3d(i,k) + epsm1*es) )
+                   rsl = xl*qsat / (r_v*tlk**2)   ! slope of C-C curve at t = tl
+                                                  ! CB02, Eqn. 4
+                   qt  = qc(i,k) + qi(i,k) + qv(i,k) !total water
+                   cpm = cp + qt*cpv              ! CB02, sec. 2, para. 1
+                   a   = 1./(1. + xl*rsl/cpm)     ! CB02 variable "a"
+                   !Now calculate convective component of the cloud fraction:
+                   if (a > 0.0) then
+                      f = min(1.0/a, 4.0)         ! f is the vertical profile
+                   else                           ! scaling function (CB2005)
+                      f = 1.0
+                   endif
+                   sigq = 1.5E-3 * ud_mf(i,k)/dt * f
+                   !sigq = 3.E-3 * ud_mf(i,k)/dt * f
+                   sigq = SQRT(sigq**2 + 1e-10)   ! combined conv + background components
+                   qmq  = a * (qt - qsat)         ! saturation deficit/excess;
+                                                  !   the numerator of Q1
+                   cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.0),0.99)
+                   if (h2oliq .lt. 1e-9) cb_cf = 0.0
+                   if (do_mynnedmf .and. qmq .ge. 0.0) then
+                      ! leverage C-B stratus clouds from MYNN in saturated conditions
+                      if (cb_cf .gt. 0.0) then
+                         clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
+                      else
+                         !default to MYNN clouds - already specified
+                      endif
+                   else                           ! unsaturated
+                      clouds1(i,k) = cb_cf
+                   endif
+                else
+                   !print *,'SAS with Xu-Randall cloud fraction'
+                   ! Xu-Randall (1996) cloud fraction
+                   es     = min( p3d(i,k),  fpvs( t3d(i,k) ) )  ! fpvs and prsl in pa
+                   qsat   = max( QMIN, eps*es / (p3d(i,k) + epsm1*es) )
+                   rhgrid = max( 0., min( 1.00, qv(i,k)/qsat ) )
+                   h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
+                   clwt   = 1.0e-6 * (p3d(i,k)*0.00001)
+
+                   if (h2oliq > clwt) then
+                      onemrh= max( 1.e-10, 1.0-rhgrid )
+                      tem1  = min(max((onemrh*qsat)**0.49,0.0001),1.0)  !jhan
+                      tem1  = 100.0 / tem1
+                      value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
+                      tem2  = sqrt( sqrt(rhgrid) )
+
+                      clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
+                   else
+                      clouds1(i,k) = 0.0
+                   endif
+                   !print*,"XuRandla- cf:",clouds1(i,k)," rh:",rhgrid," qt:",h2oliq
+                   !print*,"XuRandlb- clwt:",clwt," qsat:",qsat," p:",p3d(i,k)
+                endif ! end convective cf choice
+              endif ! qlc/qli check
+            enddo
+          enddo
+
+        endif ! convection scheme check
 
       endif ! timestep > 1