Add MYNN surface layer scheme

.gitattributes

@@ -0,0 +1 @@
+*.pf linguist-language=Fortran-Free-Form

.github/workflows/unit-tests.yaml

@@ -35,7 +35,7 @@ jobs:
       - name: Build atmospheric_physics
         run: |
           cmake \
-            -DCMAKE_PREFIX_PATH=/home/runner/work/atmospheric_physics/atmospheric_physics/pFUnit/build/installed \
+            -DCMAKE_PREFIX_PATH=$GITHUB_WORKSPACE/pFUnit/build/installed \
             -DATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE=ON \
             -B./build \
             -S./test/unit-test

.gitmodules

@@ -1,7 +1,7 @@
 [submodule "mmm-physics"]
 	path = schemes/mmm/mmm_physics
 	url = https://github.com/NCAR/MMM-physics.git
-	fxtag = 20250616-MPASv8.3
+	fxtag = 20251208-CAM-SIMA
 	fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/NCAR/MMM-physics.git
 [submodule "pumas"]

schemes/mmm/CMakeLists.txt

@@ -1,8 +1,8 @@
-cmake_minimum_required(VERSION 3.20)
+cmake_minimum_required(VERSION 3.21)
 
 # `mmm_physics_compat` has not been integrated into the CMake build of any top level projects yet,
 # and this CMakeLists.txt file is currently for unit testing purposes only.
-# Making a change to this CMakeLists.txt file will not impact the build of a parent project at this time.
+# Therefore, making a change to this CMakeLists.txt file will not impact the build of a parent project at this time.
 project(mmm_physics_compat
     VERSION
         0.1.0
@@ -19,11 +19,13 @@ target_sources(mmm_physics_compat
         ccpp_kind_types.F90
         mmm_physics_compat.F90
 )
-target_compile_options(mmm_physics_compat
-    PRIVATE
-        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:GNU>>:-fbacktrace -fcheck=all -std=f2018 -Wall -Wextra -Wpedantic>
-)
 target_include_directories(mmm_physics_compat
     INTERFACE
         ${CMAKE_CURRENT_BINARY_DIR}
 )
+target_compile_options(mmm_physics_compat
+    PRIVATE
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:GNU>>:-fbacktrace -fcheck=all -ffpe-trap=invalid,overflow,zero -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -std=f2018 -Wall -Wextra -Wpedantic>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:Intel>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:IntelLLVM>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+)

schemes/mmm/cu_ntiedtke_compat.F90

@@ -14,21 +14,19 @@ module cu_ntiedtke_compat
     !! \htmlinclude cu_ntiedtke_compat_pre_run.html
     subroutine cu_ntiedtke_compat_pre_run( &
             cflx, exner, landfrac, &
-            lhf, shf, &
             rthdynten, rthblten, rthratenlw, rthratensw, &
             rqvdynten, rqvblten, &
             lndj, &
-            ptf, pqvf, hfx, evap, &
+            ptf, pqvf, evap, &
             errmsg, errflg)
         use ccpp_kinds, only: kind_phys
         use ccpp_scheme_utils, only: ccpp_constituent_index
 
         real(kind_phys), intent(in) :: cflx(:, :), exner(:, :), landfrac(:), &
-                                       lhf(:), shf(:), &
                                        rthdynten(:, :), rthblten(:, :), rthratenlw(:, :), rthratensw(:, :), &
                                        rqvdynten(:, :), rqvblten(:, :)
         integer, intent(out) :: lndj(:)
-        real(kind_phys), intent(out) :: ptf(:, :), pqvf(:, :), hfx(:), evap(:)
+        real(kind_phys), intent(out) :: ptf(:, :), pqvf(:, :), evap(:)
         character(*), intent(out) :: errmsg
         integer, intent(out) :: errflg
 
@@ -42,15 +40,14 @@ subroutine cu_ntiedtke_compat_pre_run( &
 
         ptf(:, :) = (rthdynten(:, :) + rthblten(:, :) + rthratenlw(:, :) + rthratensw(:, :)) * exner(:, :)
         pqvf(:, :) = rqvdynten(:, :) + rqvblten(:, :)
-        hfx(:) = lhf(:) + shf(:)
         evap(:) = 0.0_kind_phys
 
         call ccpp_constituent_index( &
             'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', water_vapor_mixing_ratio_index, errflg, errmsg)
 
         if (errflg /= 0 .or. &
             water_vapor_mixing_ratio_index < lbound(cflx, 2) .or. water_vapor_mixing_ratio_index > ubound(cflx, 2)) then
-            errmsg = 'Failed to find desired constituent flux from cflx'
+            errmsg = 'cu_ntiedtke_compat_pre_run: Failed to find desired constituent flux from cflx'
             errflg = 1
 
             return
@@ -126,36 +123,54 @@ subroutine cu_ntiedtke_compat_run( &
         allocate(pu_local, source=pu, errmsg=errmsg, stat=errflg)
 
         if (errflg /= 0) then
+            errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pu_local"' // new_line('') // &
+                'Allocation returned with error: ' // trim(adjustl(errmsg))
+
             return
         end if
 
         allocate(pv_local, source=pv, errmsg=errmsg, stat=errflg)
 
         if (errflg /= 0) then
+            errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pv_local"' // new_line('') // &
+                'Allocation returned with error: ' // trim(adjustl(errmsg))
+
             return
         end if
 
         allocate(pt_local, source=pt, errmsg=errmsg, stat=errflg)
 
         if (errflg /= 0) then
+            errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pt_local"' // new_line('') // &
+                'Allocation returned with error: ' // trim(adjustl(errmsg))
+
             return
         end if
 
         allocate(pqv_local, source=pqv, errmsg=errmsg, stat=errflg)
 
         if (errflg /= 0) then
+            errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pqv_local"' // new_line('') // &
+                'Allocation returned with error: ' // trim(adjustl(errmsg))
+
             return
         end if
 
         allocate(pqc_local, source=pqc, errmsg=errmsg, stat=errflg)
 
         if (errflg /= 0) then
+            errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pqc_local"' // new_line('') // &
+                'Allocation returned with error: ' // trim(adjustl(errmsg))
+
             return
         end if
 
         allocate(pqi_local, source=pqi, errmsg=errmsg, stat=errflg)
 
         if (errflg /= 0) then
+            errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pqi_local"' // new_line('') // &
+                'Allocation returned with error: ' // trim(adjustl(errmsg))
+
             return
         end if
 

schemes/mmm/cu_ntiedtke_compat.meta

@@ -23,18 +23,6 @@
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent)
   intent = in
-[ lhf ]
-  standard_name = surface_upward_latent_heat_flux_from_coupler
-  units = W m-2
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent)
-  intent = in
-[ shf ]
-  standard_name = surface_upward_sensible_heat_flux_from_coupler
-  units = W m-2
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent)
-  intent = in
 [ rthdynten ]
   standard_name = tendency_of_air_potential_temperature_due_to_dynamics
   units = K s-1
@@ -89,12 +77,6 @@
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent, vertical_layer_dimension)
   intent = out
-[ hfx ]
-  standard_name = surface_upward_heat_flux
-  units = W m-2
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent)
-  intent = out
 [ evap ]
   standard_name = surface_upward_water_vapor_flux
   units = kg m-2 s-1
@@ -271,7 +253,7 @@
   dimensions = (horizontal_loop_extent)
   intent = in
 [ hfx ]
-  standard_name = surface_upward_heat_flux
+  standard_name = surface_upward_sensible_heat_flux_from_coupler
   units = W m-2
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent)

schemes/mmm/mmm_physics_compat.F90

@@ -3,6 +3,7 @@ module mmm_physics_compat
     implicit none
 
     private
+    public :: mmm_physics_compat_init
     public :: mmm_physics_compat_run
     public :: mmm_physics_accumulate_tendencies_timestep_init
     public :: mmm_physics_accumulate_tendencies_run
@@ -12,22 +13,53 @@ module mmm_physics_compat
     public :: geopotential_height_wrt_sfc_at_if_to_msl_run
     public :: geopotential_height_wrt_sfc_to_msl_run
 contains
+    !> \section arg_table_mmm_physics_compat_init Argument Table
+    !! \htmlinclude mmm_physics_compat_init.html
+    pure subroutine mmm_physics_compat_init( &
+            isfflx, isftcflx, iz0tlnd, &
+            spp_pbl, &
+            xice_threshold, &
+            errmsg, errflg)
+        use ccpp_kinds, only: kind_phys
+
+        integer, intent(out) :: isfflx, isftcflx, iz0tlnd
+        logical, intent(out) :: spp_pbl
+        real(kind_phys), intent(out) :: xice_threshold
+        character(*), intent(out) :: errmsg
+        integer, intent(out) :: errflg
+
+        errmsg = ''
+        errflg = 0
+
+        ! These options are hardcoded to the same values as in the MPAS physics driver.
+        ! There are other possible values for them in WRF, but the following combination is the only supported one in MPAS.
+        isfflx = 1
+        isftcflx = 0
+        iz0tlnd = 0
+        spp_pbl = .false.
+        xice_threshold = 0.02_kind_phys
+    end subroutine mmm_physics_compat_init
+
     !> \section arg_table_mmm_physics_compat_run Argument Table
     !! \htmlinclude mmm_physics_compat_run.html
     pure subroutine mmm_physics_compat_run( &
             nstep, &
             dt, &
             theta_curr, theta_prev, qv_curr, qv_prev, &
+            icefrac, xice_threshold, landfrac, &
             scheme_name, &
             rthdynten, rqvdynten, &
+            xland, &
             errmsg, errflg)
         use ccpp_kinds, only: kind_phys
 
         integer, intent(in) :: nstep
         real(kind_phys), intent(in) :: dt, &
-                                       theta_curr(:, :), theta_prev(:, :), qv_curr(:, :), qv_prev(:, :)
+                                       theta_curr(:, :), theta_prev(:, :), qv_curr(:, :), qv_prev(:, :), &
+                                       icefrac(:), xice_threshold, landfrac(:)
         character(256), intent(out) :: scheme_name
-        real(kind_phys), intent(out) :: rthdynten(:, :), rqvdynten(:, :)
+        real(kind_phys), intent(out) :: rthdynten(:, :), rqvdynten(:, :), &
+                                        xland(:)
         character(*), intent(out) :: errmsg
         integer, intent(out) :: errflg
 
@@ -43,6 +75,16 @@ pure subroutine mmm_physics_compat_run( &
             rthdynten(:, :) = (theta_curr(:, :) - theta_prev(:, :)) / dt
             rqvdynten(:, :) = (qv_curr(:, :) - qv_prev(:, :)) / dt
         end if
+
+        ! For MMM physics, land mask (`xland`) is defined as
+        ! * xland = 1.0 for land cells, including sea ice cells.
+        ! * xland = 2.0 for water cells.
+        where (landfrac >= 0.5_kind_phys .or. &
+               icefrac >= xice_threshold)
+            xland = 1.0_kind_phys
+        elsewhere
+            xland = 2.0_kind_phys
+        end where
     end subroutine mmm_physics_compat_run
 
     !> \section arg_table_mmm_physics_accumulate_tendencies_timestep_init Argument Table
@@ -196,7 +238,8 @@ end subroutine mmm_physics_persist_states_timestep_final
     !> \section arg_table_compute_characteristic_grid_length_scale_init Argument Table
     !! \htmlinclude compute_characteristic_grid_length_scale_init.html
     pure subroutine compute_characteristic_grid_length_scale_init( &
-            omega, rearth, dx, &
+            omega, rearth, &
+            dx, &
             errmsg, errflg)
         use ccpp_kinds, only: kind_phys
 

schemes/mmm/mmm_physics_compat.meta

@@ -2,6 +2,54 @@
   name = mmm_physics_compat
   type = scheme
 
+[ccpp-arg-table]
+  name = mmm_physics_compat_init
+  type = scheme
+[ isfflx ]
+  standard_name = control_for_surface_flux_for_mmm_scheme
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+[ isftcflx ]
+  standard_name = control_for_surface_roughness_length_over_water_for_mmm_scheme
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+[ iz0tlnd ]
+  standard_name = control_for_surface_roughness_length_over_land_for_mmm_scheme
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+[ spp_pbl ]
+  standard_name = flag_for_stochastically_perturbed_parameterization_for_mynn_scheme
+  units = flag
+  type = logical
+  dimensions = ()
+  intent = out
+[ xice_threshold ]
+  standard_name = sea_ice_area_fraction_threshold_for_mmm_scheme
+  units = fraction
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = out
+[ errmsg ]
+  standard_name = ccpp_error_message
+  long_name = Error message for error handling in CCPP
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  long_name = Error flag for error handling in CCPP
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+
 [ccpp-arg-table]
   name = mmm_physics_compat_run
   type = scheme
@@ -41,6 +89,24 @@
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent, vertical_layer_dimension)
   intent = in
+[ icefrac ]
+  standard_name = sea_ice_area_fraction_from_coupler
+  units = fraction
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ xice_threshold ]
+  standard_name = sea_ice_area_fraction_threshold_for_mmm_scheme
+  units = fraction
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = in
+[ landfrac ]
+  standard_name = land_area_fraction_from_coupler
+  units = fraction
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = in
 [ scheme_name ]
   standard_name = scheme_name
   units = none
@@ -59,6 +125,12 @@
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent, vertical_layer_dimension)
   intent = out
+[ xland ]
+  standard_name = land_binary_mask_for_mmm_scheme
+  units = 1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = out
 [ errmsg ]
   standard_name = ccpp_error_message
   long_name = Error message for error handling in CCPP