Merge branch 'dev/gfdl' into dyn_split_RK2b

.github/actions/macos-setup/action.yml

@@ -16,3 +16,18 @@ runs:
         brew install netcdf-fortran
         brew install mpich
         echo "::endgroup::"
+
+    # NOTE: Floating point exceptions are currently disabled due to an error in
+    # HDF5 1.4.3.  They will be re-enabled when the default brew version has
+    # been updated to a working version.
+
+    - name: Set compiler flags
+      shell: bash
+      run: |
+        cd .testing
+        echo "FCFLAGS_DEBUG = -g -O0 -Wextra -Wno-compare-reals -fbacktrace -fcheck=bounds" >> config.mk
+        echo "FCFLAGS_REPRO = -g -O2 -fbacktrace" >> config.mk
+        echo "FCFLAGS_INIT = -finit-real=snan -finit-integer=2147483647 -finit-derived" >> config.mk
+        echo "FCFLAGS_FMS = -g -fbacktrace -O0" >> config.mk
+        cat config.mk
+        echo "::endgroup::"

.github/actions/testing-setup/action.yml

@@ -31,17 +31,6 @@ runs:
         REPORT_ERROR_LOGS=true make deps/lib/libFMS.a -s -j
         echo "::endgroup::"
 
-    - name: Store compiler flags used in Makefile
-      shell: bash
-      run: |
-        echo "::group::config.mk"
-        cd .testing
-        echo "FCFLAGS_DEBUG=-g -O0 -Wextra -Wno-compare-reals -fbacktrace -ffpe-trap=invalid,zero,overflow -fcheck=bounds" >> config.mk
-        echo "FCFLAGS_REPRO=-g -O2 -fbacktrace" >> config.mk
-        echo "FCFLAGS_INIT=-finit-real=snan -finit-integer=2147483647 -finit-derived" >> config.mk
-        cat config.mk
-        echo "::endgroup::"
-
     - name: Compile MOM6 in symmetric memory mode
       shell: bash
       run: |

.github/actions/ubuntu-setup/action.yml

@@ -17,3 +17,15 @@ runs:
         sudo apt-get install libopenmpi-dev
         sudo apt-get install linux-tools-common
         echo "::endgroup::"
+
+    - name: Store compiler flags used in Makefile
+      shell: bash
+      run: |
+        echo "::group::config.mk"
+        cd .testing
+        echo "FCFLAGS_DEBUG = -g -O0 -Wextra -Wno-compare-reals -fbacktrace -ffpe-trap=invalid,zero,overflow -fcheck=bounds" >> config.mk
+        echo "FCFLAGS_REPRO = -g -O2 -fbacktrace" >> config.mk
+        echo "FCFLAGS_INIT = -finit-real=snan -finit-integer=2147483647 -finit-derived" >> config.mk
+        echo "FCFLAGS_FMS = -g -fbacktrace -O0" >> config.mk
+        cat config.mk
+        echo "::endgroup::"

ac/Makefile.in

@@ -19,7 +19,6 @@ SRC_DIRS = @SRC_DIRS@
 
 -include Makefile.dep
 
-
 # Generate Makefile from template
 Makefile: @srcdir@/ac/Makefile.in config.status
 	./config.status
@@ -33,7 +32,7 @@ rwildcard=$(foreach d,$(wildcard $(1:=/*)),$(call rwildcard,$d,$2) $(filter $(su
 .PHONY: depend
 depend: Makefile.dep
 Makefile.dep: $(MAKEDEP) $(call rwildcard,$(SRC_DIRS),*.h *.c *.inc *.F90)
-	$(PYTHON) $(MAKEDEP) -o Makefile.dep -e $(SRC_DIRS)
+	$(PYTHON) $(MAKEDEP) $(DEFS) -o Makefile.dep -e $(SRC_DIRS)
 
 
 # Delete any files associated with configuration (including the Makefile).

ac/makedep

@@ -13,9 +13,16 @@ import sys
 # Pre-compile re searches
 re_module = re.compile(r"^ *module +([a-z_0-9]+)")
 re_use = re.compile(r"^ *use +([a-z_0-9]+)")
+re_cpp_define = re.compile(r"^ *# *define +[_a-zA-Z][_a-zA-Z0-9]")
+re_cpp_undef = re.compile(r"^ *# *undef +[_a-zA-Z][_a-zA-Z0-9]")
+re_cpp_ifdef = re.compile(r"^ *# *ifdef +[_a-zA-Z][_a-zA-Z0-9]")
+re_cpp_ifndef = re.compile(r"^ *# *ifndef +[_a-zA-Z][_a-zA-Z0-9]")
+re_cpp_if = re.compile(r"^ *# *if +")
+re_cpp_else = re.compile(r"^ *# *else")
+re_cpp_endif = re.compile(r"^ *# *endif")
 re_cpp_include = re.compile(r"^ *# *include *[<\"']([a-zA-Z_0-9\.]+)[>\"']")
 re_f90_include = re.compile(r"^ *include +[\"']([a-zA-Z_0-9\.]+)[\"']")
-re_program = re.compile(r"^ *[pP][rR][oO][gG][rR][aA][mM] +([a-zA-Z_0-9]+)")
+re_program = re.compile(r"^ *program +([a-z_0-9]+)", re.IGNORECASE)
 re_end = re.compile(r"^ *end *(module|procedure) ", re.IGNORECASE)
 # NOTE: This excludes comments and tokens with substrings containing `function`
 # or `subroutine`, but will fail if the keywords appear in other contexts.
@@ -26,7 +33,7 @@ re_procedure = re.compile(
 
 
 def create_deps(src_dirs, skip_dirs, makefile, debug, exec_target, fc_rule,
-                link_externals, script_path):
+                link_externals, defines):
     """Create "makefile" after scanning "src_dis"."""
 
     # Scan everything Fortran related
@@ -66,7 +73,7 @@ def create_deps(src_dirs, skip_dirs, makefile, debug, exec_target, fc_rule,
     o2mods, o2uses, o2h, o2inc, o2prg, prg2o, mod2o = {}, {}, {}, {}, {}, {}, {}
     externals, all_modules = [], []
     for f in F90_files:
-        mods, used, cpp, inc, prg, has_externals = scan_fortran_file(f)
+        mods, used, cpp, inc, prg, has_externals = scan_fortran_file(f, defines)
         # maps object file to modules produced
         o2mods[object_file(f)] = mods
         # maps module produced to object file
@@ -272,10 +279,16 @@ def nested_inc(inc_files, f2F):
     return inc_files + sorted(set(hlst)), used_mods
 
 
-def scan_fortran_file(src_file):
+def scan_fortran_file(src_file, defines=None):
     """Scan the Fortran file "src_file" and return lists of module defined,
     module used, and files included."""
     module_decl, used_modules, cpp_includes, f90_includes, programs = [], [], [], [], []
+
+    cpp_defines = defines if defines is not None else []
+
+    cpp_macros = [define.split('=')[0] for define in cpp_defines]
+    cpp_group_stack = []
+
     with io.open(src_file, 'r', errors='replace') as file:
         lines = file.readlines()
 
@@ -285,7 +298,72 @@ def scan_fortran_file(src_file):
         file_has_externals = False
             # True if the file contains any external objects
 
+        cpp_exclude = False
+            # True if the parser excludes the subsequent lines
+
+        cpp_group_stack = []
+            # Stack of condition group exclusion states
+
         for line in lines:
+            # Start of #ifdef condition group
+            match = re_cpp_ifdef.match(line)
+            if match:
+                cpp_group_stack.append(cpp_exclude)
+
+                # If outer group is excluding or macro is missing, then exclude
+                macro = line.lstrip()[1:].split()[1]
+                cpp_exclude = cpp_exclude or macro not in cpp_macros
+
+            # Start of #ifndef condition group
+            match = re_cpp_ifndef.match(line)
+            if match:
+                cpp_group_stack.append(cpp_exclude)
+
+                # If outer group is excluding or macro is present, then exclude
+                macro = line.lstrip()[1:].split()[1]
+                cpp_exclude = cpp_exclude or macro in cpp_macros
+
+            # Start of #if condition group
+            match = re_cpp_if.match(line)
+            if match:
+                cpp_group_stack.append(cpp_exclude)
+
+                # XXX: Don't attempt to parse #if statements, but store the state.
+                # if/endif stack.  For now, assume that these always fail.
+                cpp_exclude = False
+
+            # Complement #else condition group
+            match = re_cpp_else.match(line)
+            if match:
+                # Reverse the exclude state, if there is no outer exclude state
+                outer_grp_exclude = cpp_group_stack and cpp_group_stack[-1]
+                cpp_exclude = not cpp_exclude or outer_grp_exclude
+
+            # Restore exclude state when exiting conditional block
+            match = re_cpp_endif.match(line)
+            if match:
+                cpp_exclude = cpp_group_stack.pop()
+
+            # Skip lines inside of false condition blocks
+            if cpp_exclude:
+                continue
+
+            # Activate a new macro (ignoring the value)
+            match = re_cpp_define.match(line)
+            if match:
+                new_macro = line.lstrip()[1:].split()[1]
+                cpp_macros.append(new_macro)
+
+            # Deactivate a macro
+            match = re_cpp_undef.match(line)
+            if match:
+                new_macro = line.lstrip()[1:].split()[1]
+                try:
+                    cpp_macros.remove(new_macro)
+                except:
+                    # Ignore missing macros (for now?)
+                    continue
+
             match = re_module.match(line.lower())
             if match:
                 if match.group(1) not in 'procedure':   # avoid "module procedure" statements
@@ -404,8 +482,13 @@ parser.add_argument(
     action='append',
     help="Skip directory in source code search."
 )
+parser.add_argument(
+    '-D', '--define',
+    action='append',
+    help="Apply preprocessor define macros (of the form -DMACRO[=value])",
+)
 args = parser.parse_args()
 
 # Do the thing
 create_deps(args.path, args.skip, args.makefile, args.debug, args.exec_target,
-            args.fc_rule, args.link_externals, sys.argv[0])
+            args.fc_rule, args.link_externals, args.define)

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -1522,7 +1522,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
   endif
   call get_param(param_file, mdl, "SPEAR_DTFREEZE_DS", CS%SPEAR_dTf_dS, &
                  "The derivative of the freezing temperature with salinity.", &
-                 units="deg C PSU-1", default=-0.054, scale=US%degC_to_C*US%S_to_ppt, &
+                 units="degC ppt-1", default=-0.054, scale=US%degC_to_C*US%S_to_ppt, &
                  do_not_log=.not.CS%trestore_SPEAR_ECDA)
   call get_param(param_file, mdl, "RESTORE_FLUX_RHO", CS%rho_restore, &
                  "The density that is used to convert piston velocities into salt or heat "//&

config_src/drivers/solo_driver/MOM_surface_forcing.F90

@@ -1907,19 +1907,19 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
       call get_param(param_file, mdl, "SST_NORTH", CS%T_north, &
                  "With buoy_config linear, the sea surface temperature "//&
                  "at the northern end of the domain toward which to "//&
-                 "to restore.", units="deg C", default=0.0, scale=US%degC_to_C)
+                 "to restore.", units="degC", default=0.0, scale=US%degC_to_C)
       call get_param(param_file, mdl, "SST_SOUTH", CS%T_south, &
                  "With buoy_config linear, the sea surface temperature "//&
                  "at the southern end of the domain toward which to "//&
-                 "to restore.", units="deg C", default=0.0, scale=US%degC_to_C)
+                 "to restore.", units="degC", default=0.0, scale=US%degC_to_C)
       call get_param(param_file, mdl, "SSS_NORTH", CS%S_north, &
                  "With buoy_config linear, the sea surface salinity "//&
                  "at the northern end of the domain toward which to "//&
-                 "to restore.", units="PSU", default=35.0, scale=US%ppt_to_S)
+                 "to restore.", units="ppt", default=35.0, scale=US%ppt_to_S)
       call get_param(param_file, mdl, "SSS_SOUTH", CS%S_south, &
                  "With buoy_config linear, the sea surface salinity "//&
                  "at the southern end of the domain toward which to "//&
-                 "to restore.", units="PSU", default=35.0, scale=US%ppt_to_S)
+                 "to restore.", units="ppt", default=35.0, scale=US%ppt_to_S)
     endif
     call get_param(param_file, mdl, "RESTORE_FLUX_RHO", CS%rho_restore, &
                  "The density that is used to convert piston velocities into salt or heat "//&

src/core/MOM_open_boundary.F90

@@ -3919,7 +3919,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
           segment%Htot(I,j) = segment%Htot(I,j) + segment%h(I,j,k)
           segment%dZtot(I,j) = segment%dZtot(I,j) + dz(i+ishift,j,k)
         enddo
-        segment%Cg(I,j) = sqrt(GV%g_prime(1) * segment%dZtot(I,j))
+        segment%Cg(I,j) = sqrt(GV%g_prime(1) * max(0.0, segment%dZtot(I,j)))
       enddo
     else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S)
       allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB), source=0.0)
@@ -3933,7 +3933,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
           segment%Htot(i,J) = segment%Htot(i,J) + segment%h(i,J,k)
           segment%dZtot(i,J) = segment%dZtot(i,J) + dz(i,j+jshift,k)
         enddo
-        segment%Cg(i,J) = sqrt(GV%g_prime(1) * segment%dZtot(i,J))
+        segment%Cg(i,J) = sqrt(GV%g_prime(1) * max(0.0, segment%dZtot(i,J)))
       enddo
     endif
 

src/equation_of_state/MOM_EOS.F90

@@ -1517,8 +1517,8 @@ subroutine EOS_init(param_file, EOS, US)
                  "temperature.", units="kg m-3 K-1", default=-0.2)
     call get_param(param_file, mdl, "DRHO_DS", EOS%dRho_dS, &
                  "When EQN_OF_STATE="//trim(EOS_LINEAR_STRING)//", "//&
-                 "this is the partial derivative of density with "//&
-                 "salinity.", units="kg m-3 PSU-1", default=0.8)
+                 "this is the partial derivative of density with salinity.", &
+                 units="kg m-3 ppt-1", default=0.8)
     call EOS_manual_init(EOS, form_of_EOS=EOS_LINEAR, Rho_T0_S0=EOS%Rho_T0_S0, dRho_dT=EOS%dRho_dT, dRho_dS=EOS%dRho_dS)
   endif
   if (EOS%form_of_EOS == EOS_WRIGHT) then
@@ -1563,17 +1563,17 @@ subroutine EOS_init(param_file, EOS, US)
     call get_param(param_file, mdl, "TFREEZE_S0_P0",EOS%TFr_S0_P0, &
                  "When TFREEZE_FORM="//trim(TFREEZE_LINEAR_STRING)//", "//&
                  "this is the freezing potential temperature at "//&
-                 "S=0, P=0.", units="deg C", default=0.0)
+                 "S=0, P=0.", units="degC", default=0.0)
     call get_param(param_file, mdl, "DTFREEZE_DS",EOS%dTFr_dS, &
                  "When TFREEZE_FORM="//trim(TFREEZE_LINEAR_STRING)//", "//&
                  "this is the derivative of the freezing potential "//&
                  "temperature with salinity.", &
-                 units="deg C PSU-1", default=-0.054)
+                 units="degC ppt-1", default=-0.054)
     call get_param(param_file, mdl, "DTFREEZE_DP",EOS%dTFr_dP, &
                  "When TFREEZE_FORM="//trim(TFREEZE_LINEAR_STRING)//", "//&
                  "this is the derivative of the freezing potential "//&
                  "temperature with pressure.", &
-                 units="deg C Pa-1", default=0.0)
+                 units="degC Pa-1", default=0.0)
   endif
 
   if ((EOS%form_of_EOS == EOS_TEOS10 .or. EOS%form_of_EOS == EOS_ROQUET_RHO .or. &
@@ -1694,7 +1694,7 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS)
   type(EOS_type),        intent(in)    :: EOS !< Equation of state structure
 
   real, parameter :: Sref_Sprac = (35.16504/35.0) ! The TEOS 10 conversion factor to go from
-                                    ! practical salinity to reference salinity [nondim]
+                                    ! practical salinity to reference salinity [PSU ppt-1]
   integer :: i, j, k
 
   if ((EOS%form_of_EOS /= EOS_TEOS10) .and. (EOS%form_of_EOS /= EOS_ROQUET_RHO) .and. &
@@ -1808,20 +1808,20 @@ end subroutine pot_temp_to_cons_temp
 !! temperature uses this same scaling, but this can be replaced by the factor given by scale.
 subroutine abs_saln_to_prac_saln(S, prSaln, EOS, dom, scale)
   real, dimension(:), intent(in)    :: S        !< Absolute salinity [S ~> ppt]
-  real, dimension(:), intent(inout) :: prSaln   !< Practical salinity [S ~> ppt]
+  real, dimension(:), intent(inout) :: prSaln   !< Practical salinity [S ~> PSU]
   type(EOS_type),     intent(in)    :: EOS      !< Equation of state structure
   integer, dimension(2), optional, intent(in) :: dom  !< The domain of indices to work on, taking
                                                 !! into account that arrays start at 1.
   real,     optional, intent(in)    :: scale    !< A multiplicative factor by which to scale the output
-                                                !! practical in place of with scaling stored
+                                                !! practical salinities in place of with scaling stored
                                                 !! in EOS.  A value of 1.0 returns salinities in [PSU],
                                                 !! while the default is equivalent to EOS%ppt_to_S.
 
   ! Local variables
   real, dimension(size(S)) :: Sa    ! Salinity converted to [ppt]
-  real :: S_scale ! A factor to convert practical salinity from ppt to the desired units [S ppt-1 ~> 1]
+  real :: S_scale ! A factor to convert practical salinity from ppt to the desired units [S PSU-1 ~> 1]
   real, parameter :: Sprac_Sref = (35.0/35.16504) ! The TEOS 10 conversion factor to go from
-                                    ! reference salinity to practical salinity [nondim]
+                                    ! reference salinity to practical salinity [PSU ppt-1]
   integer :: i, is, ie
 
   if (present(dom)) then
@@ -1848,21 +1848,21 @@ end subroutine abs_saln_to_prac_saln
 !! use the dimensionally rescaling as specified within the EOS type.  The output potential
 !! temperature uses this same scaling, but this can be replaced by the factor given by scale.
 subroutine prac_saln_to_abs_saln(S, absSaln, EOS, dom, scale)
-  real, dimension(:), intent(in)    :: S        !< Practical salinity [S ~> ppt]
+  real, dimension(:), intent(in)    :: S        !< Practical salinity [S ~> PSU]
   real, dimension(:), intent(inout) :: absSaln  !< Absolute salinity [S ~> ppt]
   type(EOS_type),     intent(in)    :: EOS      !< Equation of state structure
   integer, dimension(2), optional, intent(in) :: dom  !< The domain of indices to work on, taking
                                                 !! into account that arrays start at 1.
   real,     optional, intent(in)    :: scale    !< A multiplicative factor by which to scale the output
-                                                !! practical in place of with scaling stored
-                                                !! in EOS.  A value of 1.0 returns salinities in [PSU],
+                                                !! absolute salnities in place of with scaling stored
+                                                !! in EOS.  A value of 1.0 returns salinities in [ppt],
                                                 !! while the default is equivalent to EOS%ppt_to_S.
 
   ! Local variables
   real, dimension(size(S)) :: Sp    ! Salinity converted to [ppt]
-  real :: S_scale ! A factor to convert practical salinity from ppt to the desired units [S ppt-1 ~> 1]
+  real :: S_scale ! A factor to convert absolute salinity from ppt to the desired units [S ppt-1 ~> 1]
   real, parameter :: Sref_Sprac = (35.16504/35.0) ! The TEOS 10 conversion factor to go from
-                                    ! practical salinity to reference salinity [nondim]
+                                    ! practical salinity to reference salinity [PSU ppt-1]
   integer :: i, is, ie
 
   if (present(dom)) then