Merge branch 'dev/gfdl' into NonBous_Leith_QG_visc

.testing/tools/parse_perf.py

@@ -58,9 +58,16 @@ def parse_perf_report(perf_data_path):
 
             # get per-symbol count
             else:
-                tokens = line.split()
-                symbol = tokens[2]
-                period = int(tokens[3])
+                try:
+                    tokens = line.split()
+                    symbol = tokens[2]
+                    period = int(tokens[3])
+                except ValueError:
+                    print("parse_perf.py: Error extracting symbol count",
+                            file=sys.stderr)
+                    print("line:", repr(line), file=sys.stderr)
+                    print("tokens:", tokens, file=sys.stderr)
+                    raise
 
                 profile[event_name]['symbol'][symbol] = period
 

ac/makedep

@@ -34,6 +34,8 @@ re_procedure = re.compile(
 
 
 # Preprocessor expression tokenization
+# NOTE: Labels and attributes could be assigned here, but for now we just use
+#   the token string as the label.
 cpp_scanner = re.Scanner([
   (r'defined', lambda scanner, token: token),
   (r'[_A-Za-z][_0-9a-zA-Z]*', lambda scanner, token: token),
@@ -56,13 +58,15 @@ cpp_scanner = re.Scanner([
   (r'&', lambda scanner, token: token),
   (r'\|\|', lambda scanner, token: token),
   (r'\|', lambda scanner, token: token),
-  (r'^\#if', None),
+  (r'^ *\# *if', None),
   (r'\s+', None),
 ])
 
 
 cpp_operate = {
+    '(': lambda x: x,
     '!': lambda x: not x,
+    'defined': lambda x, y: x in y,
     '*': lambda x, y: x * y,
     '/': lambda x, y: x // y,
     '+': lambda x, y: x + y,
@@ -85,6 +89,7 @@ cpp_operate = {
 cpp_op_rank = {
     '(': 13,
     '!': 12,
+    'defined': 12,
     '*': 11,
     '/': 11,
     '+': 10,
@@ -527,7 +532,7 @@ def cpp_expr_eval(expr, macros=None):
     if macros is None:
         macros = {}
 
-    results, remainder = cpp_scanner.scan(expr)
+    results, remainder = cpp_scanner.scan(expr.strip())
 
     # Abort if any characters are not tokenized
     if remainder:
@@ -545,72 +550,59 @@ def cpp_expr_eval(expr, macros=None):
 
     tokens = iter(results)
     for tok in tokens:
-        # Evaluate "defined()" statements
-        if tok == 'defined':
-            tok = next(tokens)
-
-            parens = tok == '('
-            if parens:
-                tok = next(tokens)
+        if tok in cpp_op_rank.keys():
+            while cpp_op_rank[tok] <= cpp_op_rank[prior_op]:
 
-            # NOTE: Any key in `macros` is considered to be set, even if the
-            # value is None.
-            value = tok in macros
+                # Unary operators are "look ahead" so we always skip them.
+                # (However, `op` below could be a unary operator.)
+                if tok in ('!', 'defined', '('):
+                    break
 
-            # Negation
-            while prior_op == '!':
+                second = stack.pop()
                 op = stack.pop()
-                assert op == '!'
-                value = cpp_operate[op](value)
-                prior_op = stack[-1] if stack else None
-
-            stack.append(value)
 
-            if parens:
-                tok = next(tokens)
-                assert tok == ')'
+                if op == '(':
+                    value = second
 
-        elif tok.isdigit():
-            value = int(tok)
-            stack.append(value)
+                elif op == '!':
+                    if isinstance(second, str):
+                        if second.isidentifier():
+                            second = macros.get(second, '0')
+                        if second.isdigit():
+                            second = int(second)
 
-        elif tok.isidentifier():
-            # "Identifiers that are not macros, which are all considered to be
-            # the number zero." (CPP manual, 4.2.2)
-            value = macros.get(tok, '0')
-            if value.isdigit():
-                value = int(value)
-            stack.append(value)
+                    value = cpp_operate[op](second)
 
-        elif tok in cpp_op_rank.keys():
-            while cpp_op_rank[tok] <= cpp_op_rank[prior_op]:
+                elif op == 'defined':
+                    value = cpp_operate[op](second, macros)
 
-                # Skip unary prefix operators (only '!' at the moment)
-                if tok == '!':
-                    break
+                else:
+                    first = stack.pop()
 
-                second = stack.pop()
-                op = stack.pop()
-                first = stack.pop()
+                    if isinstance(first, str):
+                        if first.isidentifier():
+                            first = macros.get(first, '0')
+                        if first.isdigit():
+                            first = int(first)
 
-                value = cpp_operate[op](first, second)
-                prior_op = stack[-1] if stack else None
+                    if isinstance(second, str):
+                        if second.isidentifier():
+                            second = macros.get(second, '0')
+                        if second.isdigit():
+                            second = int(second)
 
-                if prior_op == '(':
-                    prior_op = None
-                    if tok == ')':
-                        stack.pop()
+                    value = cpp_operate[op](first, second)
 
+                prior_op = stack[-1] if stack else None
                 stack.append(value)
 
-            if tok == ')':
-                prior_op = stack[-2] if stack and len(stack) > 1 else None
-            else:
+            # The ) "operator" has already been applied, so it can be dropped.
+            if tok != ')':
                 stack.append(tok)
                 prior_op = tok
 
-                if prior_op in ('(',):
-                    prior_op = None
+        elif tok.isdigit() or tok.isidentifier():
+            stack.append(tok)
 
         else:
             print("Unsupported token:", tok)

src/ALE/MOM_ALE.F90

@@ -40,7 +40,7 @@ module MOM_ALE
 use MOM_remapping,        only : remapping_core_h, remapping_core_w
 use MOM_remapping,        only : remappingSchemesDoc, remappingDefaultScheme
 use MOM_remapping,        only : interpolate_column, reintegrate_column
-use MOM_remapping,        only : remapping_CS, dzFromH1H2
+use MOM_remapping,        only : remapping_CS, dzFromH1H2, remapping_set_param
 use MOM_string_functions, only : uppercase, extractWord, extract_integer
 use MOM_tracer_registry,  only : tracer_registry_type, tracer_type, MOM_tracer_chkinv
 use MOM_unit_scaling,     only : unit_scale_type
@@ -147,6 +147,7 @@ module MOM_ALE
 public pre_ALE_adjustments
 public ALE_remap_init_conds
 public ALE_register_diags
+public ALE_set_extrap_boundaries
 
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
@@ -176,6 +177,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   logical           :: force_bounds_in_subcell
   logical           :: local_logical
   logical           :: remap_boundary_extrap
+  logical           :: init_boundary_extrap
   type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
                                                          ! for sharing parameters.
 
@@ -225,6 +227,10 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
                  "If true, values at the interfaces of boundary cells are "//&
                  "extrapolated instead of piecewise constant", default=.false.)
+  call get_param(param_file, mdl, "INIT_BOUNDARY_EXTRAP", init_boundary_extrap, &
+                 "If true, values at the interfaces of boundary cells are "//&
+                 "extrapolated instead of piecewise constant during initialization."//&
+                 "Defaults to REMAP_BOUNDARY_EXTRAP.", default=remap_boundary_extrap)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)
@@ -237,13 +243,13 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
 
   call initialize_remapping( CS%remapCS, string, &
-                             boundary_extrapolation=remap_boundary_extrap, &
+                             boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
                              answer_date=CS%answer_date)
   call initialize_remapping( CS%vel_remapCS, vel_string, &
-                             boundary_extrapolation=remap_boundary_extrap, &
+                             boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
@@ -308,6 +314,18 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   if (CS%show_call_tree) call callTree_leave("ALE_init()")
 end subroutine ALE_init
 
+!> Sets the boundary extrapolation set for the remapping type.
+subroutine ALE_set_extrap_boundaries( param_file, CS)
+  type(param_file_type),   intent(in) :: param_file !< Parameter file
+  type(ALE_CS),            pointer    :: CS         !< Module control structure
+
+  logical :: remap_boundary_extrap
+  call get_param(param_file, "MOM_ALE", "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
+                 "If true, values at the interfaces of boundary cells are "//&
+                 "extrapolated instead of piecewise constant", default=.false.)
+  call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap)
+end subroutine ALE_set_extrap_boundaries
+
 !> Initialize diagnostics for the ALE module.
 subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   type(time_type),target,     intent(in)  :: Time  !< Time structure

src/core/MOM.F90

@@ -56,6 +56,7 @@ module MOM
 use MOM_ALE,                   only : ALE_remap_tracers, ALE_remap_velocities
 use MOM_ALE,                   only : ALE_remap_set_h_vel, ALE_remap_set_h_vel_via_dz
 use MOM_ALE,                   only : ALE_update_regrid_weights, pre_ALE_diagnostics, ALE_register_diags
+use MOM_ALE,                   only : ALE_set_extrap_boundaries
 use MOM_ALE_sponge,            only : rotate_ALE_sponge, update_ALE_sponge_field
 use MOM_barotropic,            only : Barotropic_CS
 use MOM_boundary_update,       only : call_OBC_register, OBC_register_end, update_OBC_CS
@@ -3120,8 +3121,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
       endif
     endif
   endif
-  if ( CS%use_ALE_algorithm ) call ALE_updateVerticalGridType( CS%ALE_CSp, GV )
-
+  if ( CS%use_ALE_algorithm ) then
+   call ALE_set_extrap_boundaries (param_file, CS%ALE_CSp)
+   call callTree_waypoint("returned from ALE_init() (initialize_MOM)")
+   call ALE_updateVerticalGridType( CS%ALE_CSp, GV )
+  endif
   ! The basic state variables have now been fully initialized, so update their halos and
   ! calculate any derived thermodynmics quantities.
 

src/core/MOM_CoriolisAdv.F90

@@ -902,20 +902,20 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv, Wav
         if (associated(AD%rv_x_u)) then
           do J=Jsq,Jeq ; do i=is,ie
             AD%rv_x_u(i,J,k) = -G%IdyCv(i,J) * C1_12 * &
-              ((q2(I,J) + q2(I-1,J) + q2(I-1,J-1)) * uh(I-1,j,k) + &
-               (q2(I-1,J) + q2(I,J) + q2(I,J-1)) * uh(I,j,k) + &
-               (q2(I-1,J) + q2(I,J+1) + q2(I,J)) * uh(I,j+1,k) + &
-               (q2(I,J) + q2(I-1,J+1) + q2(I-1,J)) * uh(I-1,j+1,k))
+              (((((q2(I,J) + q2(I-1,J-1)) + q2(I-1,J)) * uh(I-1,j,k)) + &
+                (((q2(I-1,J) + q2(I,J+1)) + q2(I,J)) * uh(I,j+1,k))) + &
+               ((((q2(I-1,J) + q2(I,J-1)) + q2(I,J)) * uh(I,j,k))+ &
+                (((q2(I,J) + q2(I-1,J+1)) + q2(I-1,J)) * uh(I-1,j+1,k))))
           enddo ; enddo
         endif
 
         if (associated(AD%rv_x_v)) then
           do j=js,je ; do I=Isq,Ieq
             AD%rv_x_v(I,j,k) = G%IdxCu(I,j) * C1_12 * &
-              ((q2(I+1,J) + q2(I,J) + q2(I,J-1)) * vh(i+1,J,k) + &
-               (q2(I-1,J) + q2(I,J) + q2(I,J-1)) * vh(i,J,k) + &
-               (q2(I-1,J-1) + q2(I,J) + q2(I,J-1)) * vh(i,J-1,k) + &
-               (q2(I+1,J-1) + q2(I,J) + q2(I,J-1)) * vh(i+1,J-1,k))
+              (((((q2(I+1,J) + q2(I,J-1)) + q2(I,J)) * vh(i+1,J,k)) + &
+                (((q2(I-1,J-1) + q2(I,J)) + q2(I,J-1)) * vh(i,J-1,k))) + &
+               ((((q2(I-1,J) + q2(I,J-1)) + q2(I,J)) * vh(i,J,k)) + &
+                (((q2(I+1,J-1) + q2(I,J)) + q2(I,J-1)) * vh(i+1,J-1,k))))
           enddo ; enddo
         endif
       endif