Updates to Thompson-Eidhammer microphysics and MPAS to support future…

… microphysics features (#2) * Add microphysics submodule and updated graupel lookup table support (#1) * microphysics Makefiles for submodule * Add submodule support for microphysics * Fix bug in Makefile directory path * Thompson table updates for two-moment graupel * Add aerosol aware microphysics capability (#2) * Add microphysics aerosol information to initialization * Link microphysics aerosol tables for runtime * Update submodule * Fix link to microphysics aerosol tables * Add water- and ice-friendly aerosols to model for transport (#3) * Update pointer for microphysics * Update pointer to microphysics submodule * Add microphysics aerosol information to LBCs * Add aerosol-aware microphysics to LBCs (#4) * Add microphysics aerosol information to LBCs * Support for water- and ice-friendly aerosols in LBCs * Add documentation to subroutine init_microphysics_aerosols_lbc * Time interpolation for microphysics aerosols ICs * Add time interpolation for microphysics aerosols to LBCs * Aerosol time interpolation (#5) * Add microphysics aerosol information to LBCs * Time interpolation for microphysics aerosols ICs * Add time interpolation for microphysics aerosols to LBCs * Nc microphysics (#6) * Updates for aerosol aware * Add support for nc microphysics * Add nc to interface * point to develop branch of Thompson submodule * Import nwfa and nifa into microphysics (#7) * Add microphysics submodule makefile (#8) * Import nwfa and nifa into microphysics * Add makefile support for microphysics submodule * Update location for microphysics submodule * remove GSL_cloud_physics directory * Point to microphysics submodule * Point to microphysics submodule * Remove old submodule * Bug fix in Makefile to deal with old version of microphysics
NOAA-GSL · Mar 29, 2024 · f95490e · f95490e
1 parent cbba5a4
commit f95490e
Show file tree

Hide file tree

Showing 14 changed files with 552 additions and 20 deletions.
diff --git a/.gitmodules b/.gitmodules
@@ -0,0 +1,4 @@
+[submodule "src/core_atmosphere/physics/GSL_cloud_physics"]
+	path = src/core_atmosphere/physics/GSL_cloud_physics
+	url = https://github.com/AndersJensen-NOAA/GSL_cloud_physics.git
+	branch = develop
diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
@@ -1480,9 +1480,21 @@
                              description="Cloud ice number concentration"
                              packages="bl_mynn_in;mp_thompson_in"/>
 
+                        <var name="nc" array_group="number" units="nb kg^{-1}"
+                             description="Cloud number concentration"
+                             packages="mp_thompson_in"/>
+
                         <var name="nr" array_group="number" units="nb kg^{-1}"
                              description="Rain number concentration"
                              packages="mp_thompson_in"/>
+
+			<var name="nwfa" array_group="passive" units="nb kg^{-1}"
+                             description="Number of water-friendly aerosols for microphysics"
+                             packages="mp_thompson_in"/>
+
+			<var name="nifa" array_group="passive" units="nb kg^{-1}"
+                             description="Number of ice-friendly aerosols for microphysics"
+                             packages="mp_thompson_in"/>
                 </var_array>
 #endif
 
@@ -1802,9 +1814,21 @@
                              description="Tendency of cloud ice number concentration multiplied by dry air density divided by d(zeta)/dz"
                              packages="bl_mynn_in;mp_thompson_in"/>
 
+			<var name="tend_nc" name_in_code="nc" array_group="number" units="nb m^{-3} s^{-1}" 
+                             description="Tendency of cloud number concentration multiplied by dry air density divided by d(zeta)/dz"
+                             packages="mp_thompson_in"/>
+
                         <var name="tend_nr" name_in_code="nr" array_group="number" units="nb m^{-3} s^{-1}" 
                              description="Tendency of rain number concentration multiplied by dry air density divided by d(zeta)/dz"
                              packages="mp_thompson_in"/>
+
+			<var name="tend_nwfa" name_in_code="nwfa" array_group="passive" units="nb m^{-3} s^{-1}" 
+                             description="Tendency of water-friendly aerosols multiplied by dry air density divided by d(zeta)/dz"
+                             packages="mp_thompson_in"/>
+
+			<var name="tend_nifa" name_in_code="nifa" array_group="passive" units="nb m^{-3} s^{-1}" 
+                             description="Tendency of ice-friendly aerosols multiplied by dry air density divided by d(zeta)/dz"
+                             packages="mp_thompson_in"/>
                 </var_array>
 #endif
         </var_struct>
@@ -1858,11 +1882,20 @@
                         <var name="lbc_qg" name_in_code="qg" array_group="moist"  packages="mp_thompson_in;mp_wsm6_in" units="kg kg^{-1} s^{-1}"
                              description="Lateral boundary tendency of graupel mixing ratio"/>
 
+                        <var name="lbc_nc" name_in_code="nc" array_group="number" packages="mp_thompson_in" units="m^{-3} s^{-1}"
+                             description="Lateral boundary tendency of cloud number concentration"/>
+
                         <var name="lbc_nr" name_in_code="nr" array_group="number" packages="mp_thompson_in" units="m^{-3} s^{-1}"
                              description="Lateral boundary tendency of rain number concentration"/>
 
                         <var name="lbc_ni" name_in_code="ni" array_group="number" packages="bl_mynn_in;mp_thompson_in" units="m^{-3} s^{-1}"
                              description="Lateral boundary tendency of ice number concentration"/>
+
+                        <var name="lbc_nwfa" name_in_code="nwfa" array_group="passive" packages="mp_thompson_in" units="m^{-3} s^{-1}"
+                             description="Lateral boundary tendency of wfa"/>
+
+                        <var name="lbc_nifa" name_in_code="nifa" array_group="passive" packages="mp_thompson_in" units="m^{-3} s^{-1}"
+                             description="Lateral boundary tendency of ifa"/>
                 </var_array>
         </var_struct>
 

diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
@@ -367,7 +367,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
       character (len=StrKIND), pointer :: config_convection_scheme
 
       integer, pointer :: num_scalars, index_qv, nCells, nCellsSolve, nEdges, nEdgesSolve, nVertices, nVerticesSolve, nVertLevels
-      integer, pointer :: index_qc, index_qr, index_qi, index_qs, index_qg, index_nr, index_ni
+      integer, pointer :: index_qc, index_qr, index_qi, index_qs, index_qg, index_nc, index_nr, index_ni, index_nwfa, index_nifa
 
       character(len=StrKIND), pointer :: config_IAU_option
 
@@ -468,9 +468,12 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
          call mpas_pool_get_dimension(state, 'index_qr', index_qr)
          call mpas_pool_get_dimension(state, 'index_qi', index_qi)
          call mpas_pool_get_dimension(state, 'index_qs', index_qs)
-         call mpas_pool_get_dimension(state, 'index_qg', index_qg)  
+         call mpas_pool_get_dimension(state, 'index_qg', index_qg)
+         call mpas_pool_get_dimension(state, 'index_nc', index_nc)
          call mpas_pool_get_dimension(state, 'index_nr', index_nr)
          call mpas_pool_get_dimension(state, 'index_ni', index_ni)
+         call mpas_pool_get_dimension(state, 'index_nwfa', index_nwfa)
+         call mpas_pool_get_dimension(state, 'index_nifa', index_nifa)
       endif
 
       !
@@ -921,12 +924,21 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
                   if (index_qg > 0) then
                      scalars_driving(index_qg,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'qg', rk_timestep(rk_step) )
                   end if
+                  if (index_nc > 0) then
+                     scalars_driving(index_nc,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nc', rk_timestep(rk_step) )
+                  end if
                   if (index_nr > 0) then
                      scalars_driving(index_nr,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nr', rk_timestep(rk_step) )
                   end if
                   if (index_ni > 0) then
                      scalars_driving(index_ni,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'ni', rk_timestep(rk_step) )
                   end if
+                  if (index_nwfa > 0) then
+                     scalars_driving(index_nwfa,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nwfa', rk_timestep(rk_step) )
+                  end if
+                  if (index_nifa > 0) then
+                     scalars_driving(index_nifa,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nifa', rk_timestep(rk_step) )
+                  end if
 
                   !$OMP PARALLEL DO
                   do thread=1,nThreads
@@ -1088,12 +1100,21 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
                if (index_qg > 0) then
                   scalars_driving(index_qg,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'qg', rk_timestep(rk_step) )
                end if
+               if (index_nc > 0) then
+                  scalars_driving(index_nc,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nc', rk_timestep(rk_step) )
+               end if
                if (index_nr > 0) then
                   scalars_driving(index_nr,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nr', rk_timestep(rk_step) )
                end if
                if (index_ni > 0) then
                   scalars_driving(index_ni,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'ni', rk_timestep(rk_step) )
                end if
+               if (index_nwfa > 0) then
+                  scalars_driving(index_nwfa,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nwfa', rk_timestep(rk_step) )
+               end if
+               if (index_nifa > 0) then
+                  scalars_driving(index_nifa,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nifa', rk_timestep(rk_step) )
+               end if
 
 !$OMP PARALLEL DO
                do thread=1,nThreads
@@ -1248,12 +1269,21 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
          if (index_qg > 0) then
             scalars_driving(index_qg,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'qg', dt )
          end if
+         if (index_nc > 0) then
+            scalars_driving(index_nc,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nc', dt )
+         end if
          if (index_nr > 0) then
             scalars_driving(index_nr,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nr', dt )
          end if
          if (index_ni > 0) then
             scalars_driving(index_ni,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'ni', dt )
          end if
+         if (index_nwfa > 0) then
+            scalars_driving(index_nwfa,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nwfa', dt )
+         end if
+         if (index_nifa > 0) then
+            scalars_driving(index_nifa,1:nVertLevels,1:nCells+1) = mpas_atm_get_bdy_state( clock, block, nVertLevels, nCells, 'nifa', dt )
+         end if
 
 !$OMP PARALLEL DO
          do thread=1,nThreads

diff --git a/src/core_atmosphere/physics/GSL_cloud_physics b/src/core_atmosphere/physics/GSL_cloud_physics
diff --git a/src/core_atmosphere/physics/Makefile b/src/core_atmosphere/physics/Makefile
@@ -1,10 +1,10 @@
-.SUFFIXES: .F .o
+.SUFFIXES: .F .o .F90
 
 ifeq ($(CORE),atmosphere)
 COREDEF = -Dmpas
 endif
 
-all: lookup_tables core_physics_init core_physics_mmm core_physics_wrf core_physics
+all: lookup_tables core_physics_init core_physics_mmm core_microphysics core_physics_wrf core_physics
 
 dummy:
 	echo "****** compiling physics ******"
@@ -52,13 +52,17 @@ lookup_tables:
 core_physics_mmm: core_physics_init
 	(cd physics_mmm; $(MAKE) all)
 
+core_microphysics: core_physics_init core_physics_mmm
+	(if [ -f physics_wrf/module_mp_thompson.F ]; then rm physics_wrf/module_mp_thompson.F; else echo "No old module_mp_thompson.F file found"; fi)
+	(cd GSL_cloud_physics; cp ./drivers/mpas/module_mp_thompson.F90 .; $(MAKE) all COREDEF="$(COREDEF)")
+
 core_physics_wrf: core_physics_init core_physics_mmm
 	(cd physics_wrf; $(MAKE) all COREDEF="$(COREDEF)")
 
 core_physics_init: $(OBJS_init)
 	ar -ru libphys.a $(OBJS_init)
 
-core_physics: core_physics_wrf
+core_physics: core_microphysics core_physics_wrf
 	($(MAKE) phys_interface COREDEF="$(COREDEF)")
 	ar -ru libphys.a $(OBJS)
 
@@ -199,6 +203,7 @@ clean:
 	$(RM) *.o *.mod *.f90 libphys.a
 	( cd physics_wrf; $(MAKE) clean )
 	( cd physics_mmm; $(MAKE) clean )
+	( cd GSL_cloud_physics; $(MAKE) clean )
 	@# Certain systems with intel compilers generate *.i files
 	@# This removes them during the clean process
 	$(RM) *.i
@@ -207,7 +212,16 @@ clean:
 	$(RM) $@ $*.mod
 ifeq "$(GEN_F90)" "true"
 	$(CPP) $(CPPFLAGS) $(COREDEF) $(HYDROSTATIC) $(CPPINCLUDES) $< > $*.f90
-	$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I.. -I../../framework -I../../external/esmf_time_f90
+	$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./GSL_cloud_physics -I.. -I../../framework -I../../external/esmf_time_f90
+else
+	$(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./GSL_cloud_physics -I.. -I../../framework -I../../external/esmf_time_f90
+endif
+
+.F90.o:
+	$(RM) $@ $*.mod
+ifeq "$(GEN_F90)" "true"
+	$(CPP) $(CPPFLAGS) $(COREDEF) $(HYDROSTATIC) $(CPPINCLUDES) $< > $*.f90
+	$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./GSL_cloud_physics -I.. -I../../framework -I../../external/esmf_time_f90
 else
-	$(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I.. -I../../framework -I../../external/esmf_time_f90
+	$(FC) $(CPPFLAGS) $(COREDEF) $(HYDROSATIC) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I./physics_mmm -I./physics_wrf -I./GSL_cloud_physics -I.. -I../../framework -I../../external/esmf_time_f90
 endif
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F
@@ -19,6 +19,7 @@ module mpas_atmphys_driver_microphysics
 !wrf physics:
  use module_mp_kessler
  use module_mp_thompson
+ use module_mp_thompson_utils
  use module_mp_wsm6,only: wsm6
  use mp_wsm6,only: mp_wsm6_init,refl10cm_wsm6
 
@@ -153,7 +154,10 @@ subroutine allocate_microphysics(configs)
           if(.not.allocated(ntc_p)) allocate(ntc_p(ims:ime,jms:jme))
           if(.not.allocated(muc_p)) allocate(muc_p(ims:ime,jms:jme))
           if(.not.allocated(ni_p) ) allocate(ni_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(nc_p) ) allocate(nc_p(ims:ime,kms:kme,jms:jme))
           if(.not.allocated(nr_p) ) allocate(nr_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(nwfa_p) ) allocate(nwfa_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(nifa_p) ) allocate(nifa_p(ims:ime,kms:kme,jms:jme))
 
           if(.not.allocated(rainprod_p)) allocate(rainprod_p(ims:ime,kms:kme,jms:jme))
           if(.not.allocated(evapprod_p)) allocate(evapprod_p(ims:ime,kms:kme,jms:jme))
@@ -227,7 +231,10 @@ subroutine deallocate_microphysics(configs)
           if(allocated(ntc_p)) deallocate(ntc_p)
           if(allocated(muc_p)) deallocate(muc_p)
           if(allocated(ni_p) ) deallocate(ni_p )
+          if(allocated(nc_p) ) deallocate(nc_p )
           if(allocated(nr_p) ) deallocate(nr_p )
+          if(allocated(nwfa_p) ) deallocate(nwfa_p )
+          if(allocated(nifa_p) ) deallocate(nifa_p )
 
           if(allocated(rainprod_p)) deallocate(rainprod_p)
           if(allocated(evapprod_p)) deallocate(evapprod_p)
@@ -365,6 +372,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
                   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     , &

diff --git a/src/core_atmosphere/physics/mpas_atmphys_finalize.F b/src/core_atmosphere/physics/mpas_atmphys_finalize.F
@@ -57,7 +57,7 @@ subroutine mp_thompson_deallocate
  if(allocated(tcg_racg) ) deallocate(tcg_racg ) 
  if(allocated(tmr_racg) ) deallocate(tmr_racg )
  if(allocated(tcr_gacr) ) deallocate(tcr_gacr )
- if(allocated(tmg_gacr) ) deallocate(tmg_gacr )
+! if(allocated(tmg_gacr) ) deallocate(tmg_gacr )
  if(allocated(tnr_racg) ) deallocate(tnr_racg )
  if(allocated(tnr_gacr) ) deallocate(tnr_gacr )
  if(allocated(tcs_racs1)) deallocate(tcs_racs1)