Support for NSSL 2-Moment Microphysics

.gitignore

@@ -2,3 +2,4 @@
 *.mod
 *.a
 *.pyc
+*.i90

ccpp/config/ccpp_prebuild_config.py

@@ -89,6 +89,7 @@
         },
     }
 
+
 # Add all physics scheme files relative to basedir
 SCHEME_FILES = [
     # Relative path to source (from where ccpp_prebuild.py is called) : [ list of physics sets in which scheme may be called ];
@@ -159,6 +160,7 @@
     'ccpp/physics/physics/mp_thompson_pre.F90',
     'ccpp/physics/physics/mp_thompson.F90',
     'ccpp/physics/physics/mp_thompson_post.F90',
+    'ccpp/physics/physics/mp_nssl_2mom.F90',
     'ccpp/physics/physics/ozphys.f',
     'ccpp/physics/physics/ozphys_2015.f',
     'ccpp/physics/physics/precpd.f',

ccpp/data/GFS_typedefs.F90

@@ -129,6 +129,7 @@ module GFS_typedefs
     real(kind=kind_phys) :: dt_phys              !< physics  time step in seconds
 !--- restart information
     logical :: restart                           !< flag whether this is a coldstart (.false.) or a warmstart/restart (.true.)
+    logical :: cycling                           !< flag whether this is a coldstart (.false.) or a cycled run (.true.)
 !--- hydrostatic/non-hydrostatic flag
     logical :: hydrostatic                       !< flag whether this is a hydrostatic or non-hydrostatic run
 !--- blocking data
@@ -695,6 +696,7 @@ module GFS_typedefs
     integer              :: imp_physics                    !< choice of microphysics scheme
     integer              :: imp_physics_gfdl = 11          !< choice of GFDL     microphysics scheme
     integer              :: imp_physics_thompson = 8       !< choice of Thompson microphysics scheme
+    integer              :: imp_physics_nssl = 17          !< choice of NSSL     microphysics scheme
     integer              :: imp_physics_wsm6 = 6           !< choice of WSMG     microphysics scheme
     integer              :: imp_physics_zhao_carr = 99     !< choice of Zhao-Carr microphysics scheme
     integer              :: imp_physics_zhao_carr_pdf = 98 !< choice of Zhao-Carr microphysics scheme with PDF clouds
@@ -1091,12 +1093,16 @@ module GFS_typedefs
     integer              :: ntrw            !< tracer index for rain water
     integer              :: ntsw            !< tracer index for snow water
     integer              :: ntgl            !< tracer index for graupel
+    integer              :: nthl            !< tracer index for hail
     integer              :: ntclamt         !< tracer index for cloud amount
     integer              :: ntlnc           !< tracer index for liquid number concentration
     integer              :: ntinc           !< tracer index for ice    number concentration
     integer              :: ntrnc           !< tracer index for rain   number concentration
     integer              :: ntsnc           !< tracer index for snow   number concentration
     integer              :: ntgnc           !< tracer index for graupel number concentration
+    integer              :: nthnc           !< tracer index for hail   number concentration
+    integer              :: nqvolg          !< tracer index for graupel particle volume m3 kg-1
+    integer              :: nqvolh          !< tracer index for hail   particle volume m3 kg-1
     integer              :: ntke            !< tracer index for kinetic energy
     integer              :: nto             !< tracer index for oxygen ion
     integer              :: nto2            !< tracer index for oxygen
@@ -1158,6 +1164,7 @@ module GFS_typedefs
     integer              :: kdt             !< current forecast iteration
     logical              :: first_time_step !< flag signaling first time step for time integration routine
     logical              :: restart         !< flag whether this is a coldstart (.false.) or a warmstart/restart (.true.)
+    logical              :: cycling         !< flag whether this is a coldstart (.false.) or a cycled run (.true.)
     logical              :: hydrostatic     !< flag whether this is a hydrostatic or non-hydrostatic run
     integer              :: jdat(1:8)       !< current forecast date and time
                                             !< (yr, mon, day, t-zone, hr, min, sec, mil-sec)
@@ -2952,6 +2959,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     logical              :: aux2d_time_avg(1:naux2dmax) = .false. !< flags for time averaging of auxiliary 2d arrays
     logical              :: aux3d_time_avg(1:naux3dmax) = .false. !< flags for time averaging of auxiliary 3d arrays
 
+    logical              :: cycling        = .false.         !< flag to activate extra cycling procedures
     real(kind=kind_phys) :: fhcyc          = 0.              !< frequency for surface data cycling (hours)
     integer              :: thermodyn_id   =  1              !< valid for GFS only for get_prs/phi
     integer              :: sfcpress_id    =  1              !< valid for GFS only for get_prs/phi
@@ -4203,12 +4211,16 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%ntrw             = get_tracer_index(Model%tracer_names, 'rainwat',    Model%me, Model%master, Model%debug)
     Model%ntsw             = get_tracer_index(Model%tracer_names, 'snowwat',    Model%me, Model%master, Model%debug)
     Model%ntgl             = get_tracer_index(Model%tracer_names, 'graupel',    Model%me, Model%master, Model%debug)
+    Model%nthl             = get_tracer_index(Model%tracer_names, 'hailwat',    Model%me, Model%master, Model%debug)
     Model%ntclamt          = get_tracer_index(Model%tracer_names, 'cld_amt',    Model%me, Model%master, Model%debug)
     Model%ntlnc            = get_tracer_index(Model%tracer_names, 'water_nc',   Model%me, Model%master, Model%debug)
     Model%ntinc            = get_tracer_index(Model%tracer_names, 'ice_nc',     Model%me, Model%master, Model%debug)
     Model%ntrnc            = get_tracer_index(Model%tracer_names, 'rain_nc',    Model%me, Model%master, Model%debug)
     Model%ntsnc            = get_tracer_index(Model%tracer_names, 'snow_nc',    Model%me, Model%master, Model%debug)
     Model%ntgnc            = get_tracer_index(Model%tracer_names, 'graupel_nc', Model%me, Model%master, Model%debug)
+    Model%nthnc            = get_tracer_index(Model%tracer_names, 'hail_nc',    Model%me, Model%master, Model%debug)
+    Model%nqvolg           = get_tracer_index(Model%tracer_names, 'graupel_pv', Model%me, Model%master, Model%debug)
+    Model%nqvolh           = get_tracer_index(Model%tracer_names, 'hail_pv',    Model%me, Model%master, Model%debug)
     Model%ntke             = get_tracer_index(Model%tracer_names, 'sgs_tke',    Model%me, Model%master, Model%debug)
     Model%nqrimef          = get_tracer_index(Model%tracer_names, 'q_rimef',    Model%me, Model%master, Model%debug)
     Model%ntwa             = get_tracer_index(Model%tracer_names, 'liq_aero',   Model%me, Model%master, Model%debug)
@@ -4315,6 +4327,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%kdt              = 0
     Model%first_time_step  = .true.
     Model%restart          = restart
+    Model%cycling          = cycling
     Model%hydrostatic      = hydrostatic
     Model%jdat(1:8)        = jdat(1:8)
     allocate(Model%si(Model%levr+1))
@@ -4665,6 +4678,20 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                                           ' num_p3d =',Model%num_p3d, &
                                           ' num_p2d =',Model%num_p2d
 
+    elseif (Model%imp_physics == Model%imp_physics_nssl) then ! NSSL microphysics
+      Model%npdf3d  = 0
+      Model%num_p3d = 3
+      Model%num_p2d = 1
+      Model%pdfcld  = .false.
+      Model%shcnvcw = .false.
+      Model%ncnd    = 5
+      Model%nleffr = 1
+      Model%nieffr = 2
+      Model%nseffr = 3
+      if (Model%me == Model%master) print *,' Using NSSL double moment', &
+                                          ' microphysics', ' num_p3d =',Model%num_p3d, &
+                                          ' num_p2d =',Model%num_p2d
+
     else if (Model%imp_physics == Model%imp_physics_mg) then        ! Morrison-Gettelman Microphysics
       Model%npdf3d  = 0
       Model%num_p3d = 5
@@ -4962,6 +4989,9 @@ subroutine control_print(Model)
         print *, ' ttendlim          : ', Model%ttendlim
         print *, ' '
       endif
+      if (Model%imp_physics == Model%imp_physics_nssl) then
+        print *, ' NSSL microphysical parameters'
+      endif
       if (Model%imp_physics == Model%imp_physics_mg) then
         print *, ' M-G microphysical parameters'
         print *, ' fprcp             : ', Model%fprcp
@@ -5201,12 +5231,16 @@ subroutine control_print(Model)
       print *, ' ntrw              : ', Model%ntrw
       print *, ' ntsw              : ', Model%ntsw
       print *, ' ntgl              : ', Model%ntgl
+      print *, ' nthl              : ', Model%nthl
       print *, ' ntclamt           : ', Model%ntclamt
       print *, ' ntlnc             : ', Model%ntlnc
       print *, ' ntinc             : ', Model%ntinc
       print *, ' ntrnc             : ', Model%ntrnc
       print *, ' ntsnc             : ', Model%ntsnc
       print *, ' ntgnc             : ', Model%ntgnc
+      print *, ' nthnc             : ', Model%nthnc
+      print *, ' nqvolg            : ', Model%nqvolg
+      print *, ' nqvolh            : ', Model%nqvolh
       print *, ' ntke              : ', Model%ntke
       print *, ' nto               : ', Model%nto
       print *, ' nto2              : ', Model%nto2
@@ -5254,6 +5288,7 @@ subroutine control_print(Model)
       print *, ' sec               : ', Model%sec
       print *, ' first_time_step   : ', Model%first_time_step
       print *, ' restart           : ', Model%restart
+      print *, ' cycling           : ', Model%cycling
       print *, ' hydrostatic       : ', Model%hydrostatic
     endif
 
@@ -5850,6 +5885,7 @@ end subroutine diag_create
   subroutine diag_rad_zero(Diag, Model)
     class(GFS_diag_type)               :: Diag
     type(GFS_control_type), intent(in) :: Model
+    integer :: i
 
     Diag%fluxr        = zero
     Diag%topfsw%upfxc = zero
@@ -5869,6 +5905,7 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
     logical,optional, intent(in)       :: linit, iauwindow_center
 
     logical set_totprcp
+    integer :: i
 
     !--- In/Out
     Diag%srunoff    = zero
@@ -6524,7 +6561,7 @@ subroutine interstitial_create (Interstitial, IM, Model)
     end if
 !
     ! Allocate arrays that are conditional on physics choices
-    if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson) then
+    if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson .or. Model%imp_physics == Model%imp_physics_nssl) then
        allocate (Interstitial%graupelmp  (IM))
        allocate (Interstitial%icemp      (IM))
        allocate (Interstitial%rainmp     (IM))
@@ -6678,6 +6715,11 @@ subroutine interstitial_setup_tracers(Interstitial, Model)
       Interstitial%nncl = 5
     endif
 
+    if (Model%imp_physics == Model%imp_physics_nssl) then
+       Interstitial%nvdiff = Model%ntrac
+       Interstitial%nncl = 5
+    end if
+
     if (Model%imp_physics == Model%imp_physics_mg) then
       if (abs(Model%fprcp) == 1) then
         Interstitial%nncl = 4                          ! MG2 with rain and snow
@@ -6761,13 +6803,17 @@ subroutine interstitial_setup_tracers(Interstitial, Model)
       do n=2,Model%ntrac
         if ( n /= Model%ntcw  .and. n /= Model%ntiw  .and. n /= Model%ntclamt .and. &
              n /= Model%ntrw  .and. n /= Model%ntsw  .and. n /= Model%ntrnc   .and. &
-             n /= Model%ntsnc .and. n /= Model%ntgl  .and. n /= Model%ntgnc) then
+!            n /= Model%ntlnc .and. n /= Model%ntinc .and.                          &
+             n /= Model%ntsnc .and. n /= Model%ntgl  .and. n /= Model%ntgnc   .and. &
+             n /= Model%nthl  .and. n /= Model%nthnc .and.                          &
+             n /= Model%nqvolg .and. n /= Model%nqvolh ) then
           tracers = tracers + 1
           if (Model%ntke  == n ) then
             Interstitial%otspt(tracers+1,1) = .false.
             Interstitial%ntk = tracers
           endif
-          if (Model%ntlnc == n .or. Model%ntinc == n .or. Model%ntrnc == n .or. Model%ntsnc == n .or. Model%ntgnc == n)    &
+          if (Model%ntlnc == n .or. Model%ntinc == n .or. Model%ntrnc == n .or. Model%ntsnc == n .or. &
+               Model%ntgnc == n .or. Model%nthnc == n .or. Model%nqvolg == n .or. Model%nqvolh == n)    &
 !           if (ntlnc == n .or. ntinc == n .or. ntrnc == n .or. ntsnc == n .or.&
 !               ntrw  == n .or. ntsw  == n .or. ntgl  == n)                    &
                   Interstitial%otspt(tracers+1,1) = .false.
@@ -7138,7 +7184,7 @@ subroutine interstitial_phys_reset (Interstitial, Model)
     end if
 !
     ! Reset fields that are conditional on physics choices
-    if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson) then
+    if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson .or. Model%imp_physics == Model%imp_physics_nssl) then
        Interstitial%graupelmp = clear_val
        Interstitial%icemp     = clear_val
        Interstitial%rainmp    = clear_val
@@ -7525,8 +7571,8 @@ subroutine interstitial_print(Interstitial, Model, mpirank, omprank, blkno)
     end if
 !
     ! Print arrays that are conditional on physics choices
-    if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson) then
-       write (0,*) 'Interstitial_print: values specific to GFDL/Thompson microphysics'
+    if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson .or. Model%imp_physics == Model%imp_physics_nssl) then
+       write (0,*) 'Interstitial_print: values specific to GFDL/Thompson/NSSL microphysics'
        write (0,*) 'sum(Interstitial%graupelmp    ) = ', sum(Interstitial%graupelmp       )
        write (0,*) 'sum(Interstitial%icemp        ) = ', sum(Interstitial%icemp           )
        write (0,*) 'sum(Interstitial%rainmp       ) = ', sum(Interstitial%rainmp          )

ccpp/data/GFS_typedefs.meta

@@ -195,6 +195,13 @@
   dimensions = (horizontal_loop_extent,vertical_dimension)
   type = real
   kind = kind_phys
+[qgrs(:,:,index_for_hail)]
+  standard_name = hail_mixing_ratio
+  long_name = ratio of mass of hail to mass of dry air plus vapor (without condensates)
+  units = kg kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
 [qgrs(:,:,index_for_ozone)]
   standard_name = ozone_mixing_ratio
   long_name = ozone mixing ratio
@@ -254,6 +261,27 @@
   dimensions = (horizontal_loop_extent,vertical_dimension)
   type = real
   kind = kind_phys
+[qgrs(:,:,index_for_hail_number_concentration)]
+  standard_name = hail_number_concentration
+  long_name = number concentration of hail
+  units = kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+[qgrs(:,:,index_for_graupel_particle_volume)]
+  standard_name = graupel_particle_volume
+  long_name = volume of graupel particles
+  units = kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+[qgrs(:,:,index_for_hail_particle_volume)]
+  standard_name = hail_particle_volume
+  long_name = volume of hail particles
+  units = kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
 [qgrs(:,:,index_for_turbulent_kinetic_energy)]
   standard_name = turbulent_kinetic_energy
   long_name = turbulent kinetic energy
@@ -383,6 +411,13 @@
   dimensions = (horizontal_loop_extent,vertical_dimension)
   type = real
   kind = kind_phys
+[gq0(:,:,index_for_hail)]
+  standard_name = hail_mixing_ratio_updated_by_physics
+  long_name = ratio of mass of hail to mass of dry air plus vapor (without condensates) updated by physics
+  units = kg kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
 [gq0(:,:,index_for_mass_weighted_rime_factor)]
   standard_name = mass_weighted_rime_factor_updated_by_physics
   long_name = mass weighted rime factor updated by physics
@@ -442,6 +477,27 @@
   dimensions = (horizontal_loop_extent,vertical_dimension)
   type = real
   kind = kind_phys
+[gq0(:,:,index_for_hail_number_concentration)]
+  standard_name = hail_number_concentration_updated_by_physics
+  long_name = number concentration of hail updated by physics
+  units = kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+[gq0(:,:,index_for_graupel_particle_volume)]
+  standard_name = graupel_particle_volume_updated_by_physics
+  long_name = volume of graupel particles updated by physics
+  units = kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+[gq0(:,:,index_for_hail_particle_volume)]
+  standard_name = hail_particle_volume_updated_by_physics
+  long_name = hail particle volume updated by physics
+  units = kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
 [gq0(:,:,index_for_cloud_amount)]
   standard_name = cloud_fraction_updated_by_physics
   long_name = cloud fraction updated by physics
@@ -2870,6 +2926,12 @@
   units = flag
   dimensions = ()
   type = integer
+[imp_physics_nssl]
+  standard_name = flag_for_nssl_microphysics_scheme
+  long_name = choice of NSSL microphysics scheme
+  units = flag
+  dimensions = ()
+  type = integer
 [imp_physics_wsm6]
   standard_name = flag_for_wsm6_microphysics_scheme
   long_name = choice of WSM6 microphysics scheme
@@ -4365,6 +4427,12 @@
   units = index
   dimensions = ()
   type = integer
+[nthl]
+  standard_name = index_for_hail
+  long_name = tracer index for hail
+  units = index
+  dimensions = ()
+  type = integer
 [ntclamt]
   standard_name = index_for_cloud_amount
   long_name = tracer index for cloud amount integer
@@ -4401,6 +4469,24 @@
   units = index
   dimensions = ()
   type = integer
+[nthnc]
+  standard_name = index_for_hail_number_concentration
+  long_name = tracer index for hail number concentration
+  units = index
+  dimensions = ()
+  type = integer
+[nqvolg]
+  standard_name = index_for_graupel_particle_volume
+  long_name = tracer index for graupel particle volume
+  units = index
+  dimensions = ()
+  type = integer
+[nqvolh]
+  standard_name = index_for_hail_particle_volume
+  long_name = tracer index for hail particle volume
+  units = index
+  dimensions = ()
+  type = integer
 [ntke]
   standard_name = index_for_turbulent_kinetic_energy
   long_name = tracer index for turbulent kinetic energy
@@ -4711,6 +4797,12 @@
   units = flag
   dimensions = ()
   type = logical
+[cycling]
+  standard_name = flag_for_cycling
+  long_name = flag for cycling or coldstart
+  units = flag
+  dimensions = ()
+  type = logical
 [hydrostatic]
   standard_name = flag_for_hydrostatic_solver
   long_name = flag for hydrostatic solver from dynamics