Add constituent tendency capability (#584)

Adds constituent tendency array and enables the use of individual
constituent tendency variables.

In order to facilitate time-splitting, this PR adds/enables the
following:
- ccpp_constituent_tendencies: standard name for constituent tendencies
array
- standard names of the form "tendency_of_CONSTITUENT" now handled by
the framework
- must also include "constituent = True" in metadata for the tendency

Testing:
unit tests: Added doctests to new check_tendency_variables function in
host_cap.F90
system tests: Modified advection_test to use tendency variable and
tendency array

Author: peverwhee

scripts/constituents.py

@@ -276,17 +276,28 @@ def write_constituent_routines(self, outfile, indent, suite_name, err_vars):
         for evar in err_vars:
             evar.write_def(outfile, indent+1, self, dummy=True)
         # end for
+        # Figure out how many unique (non-tendency) constituent variables we have
+        const_num = 0
+        for std_name, _ in self.items():
+            if not std_name.startswith('tendency_of_'):
+                const_num += 1
+            # end if
+        # end for
         if self:
             outfile.write("! Local variables", indent+1)
             outfile.write("integer :: index", indent+1)
-            stmt = f"allocate({self.constituent_prop_array_name()}({len(self)}))"
+            stmt = f"allocate({self.constituent_prop_array_name()}({const_num}))"
             outfile.write(stmt, indent+1)
             outfile.write("index = 0", indent+1)
         # end if
         for evar in err_vars:
             self.__init_err_var(evar, outfile, indent+1)
         # end for
         for std_name, var in self.items():
+            if std_name.startswith('tendency_of_'):
+                # Skip tendency variables
+                continue
+            # end if
             outfile.write("index = index + 1", indent+1)
             long_name = var.get_prop_value('long_name')
             units = var.get_prop_value('units')
@@ -336,7 +347,7 @@ def write_constituent_routines(self, outfile, indent, suite_name, err_vars):
         outfile.write(stmt, indent+1)
         outfile.write(f"call {self.constituent_prop_init_consts()}({local_call})", indent+2)
         outfile.write("end if", indent+1)
-        outfile.write(f"{fname} = {len(self)}", indent+1)
+        outfile.write(f"{fname} = {const_num}", indent+1)
         outfile.write(f"end function {fname}", indent)
         outfile.write(f"\n! {border}\n", 1)
         # Return the name of a constituent given an index

scripts/host_cap.py

@@ -20,7 +20,7 @@
 from metavar import CCPP_LOOP_VAR_STDNAMES
 from fortran_tools import FortranWriter
 from parse_tools import CCPPError
-from parse_tools import ParseObject, ParseSource, ParseContext
+from parse_tools import ParseObject, ParseSource, ParseContext, ParseSyntaxError
 
 ###############################################################################
 _HEADER = "cap for {host_model} calls to CCPP API"
@@ -284,6 +284,7 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
     vert_layer_dim = "vertical_layer_dimension"
     vert_interface_dim = "vertical_interface_dimension"
     array_layer = "vars_layer"
+    tend_layer = "vars_layer_tend"
     # Table preamble (leave off ccpp-table-properties header)
     ddt_mdata = [
         #"[ccpp-table-properties]",
@@ -300,6 +301,9 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
         " type = real", " kind = kind_phys"]
     # Add entries for each constituent (once per standard name)
     const_stdnames = set()
+    tend_stdnames = set()
+    const_vars = set()
+    tend_vars = set()
     for suite in suite_list:
         if run_env.verbose:
             lmsg = "Adding constituents from {} to {}"
@@ -308,12 +312,13 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
         scdict = suite.constituent_dictionary()
         for cvar in scdict.variable_list():
             std_name = cvar.get_prop_value('standard_name')
-            if std_name not in const_stdnames:
+            if std_name not in const_stdnames and std_name not in tend_stdnames:
                 # Add a metadata entry for this constituent
                 # Check dimensions and figure vertical dimension
                 # Currently, we only support variables with first dimension,
                 #   horizontal_dimension, and second (optional) dimension,
                 #   vertical_layer_dimension or vertical_interface_dimension
+                is_tend_var = 'tendency_of' in std_name
                 dims = cvar.get_dimensions()
                 if (len(dims) < 1) or (len(dims) > 2):
                     emsg = "Unsupported constituent dimensions, '{}'"
@@ -329,7 +334,11 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
                 if len(dims) > 1:
                     vdim = dims[1].split(':')[-1]
                     if vdim == vert_layer_dim:
-                        cvar_array_name = array_layer
+                        if is_tend_var:
+                            cvar_array_name = tend_layer
+                        else:
+                            cvar_array_name = array_layer
+                        # end if
                     else:
                         emsg = "Unsupported vertical constituent dimension, "
                         emsg += "'{}', must be '{}' or '{}'"
@@ -340,8 +349,13 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
                     emsg = f"Unsupported 2-D variable, '{std_name}'"
                     raise CCPPError(emsg)
                 # end if
-                # First, create an index variable for <cvar>
-                ind_std_name = "index_of_{}".format(std_name)
+                # Create an index variable for <cvar>
+                if is_tend_var:
+                    const_std_name = std_name.split("tendency_of_")[1]
+                else:
+                    const_std_name = std_name
+                # end if
+                ind_std_name = f"index_of_{const_std_name}"
                 loc_name = f"{cvar_array_name}(:,:,{ind_std_name})"
                 ddt_mdata.append(f"[ {loc_name} ]")
                 ddt_mdata.append(f" standard_name = {std_name}")
@@ -352,10 +366,44 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
                 vtype = cvar.get_prop_value('type')
                 vkind = cvar.get_prop_value('kind')
                 ddt_mdata.append(f" type = {vtype} | kind = {vkind}")
-                const_stdnames.add(std_name)
+                if is_tend_var:
+                    tend_vars.add(cvar)
+                    tend_stdnames.add(std_name)
+                else:
+                    const_vars.add(cvar)
+                    const_stdnames.add(std_name)
+                # end if
+
             # end if
         # end for
     # end for
+    # Check that all tendency variables are valid
+    for tendency_variable in tend_vars:
+        tend_stdname = tendency_variable.get_prop_value('standard_name')
+        tend_const_name = tend_stdname.split('tendency_of_')[1]
+        found = False
+        # Find the corresponding constituent variable
+        for const_variable in const_vars:
+            const_stdname = const_variable.get_prop_value('standard_name')
+            if const_stdname == tend_const_name:
+                found = True
+                compat = tendency_variable.compatible(const_variable, run_env, is_tend=True)
+                if not compat:
+                    errstr = f"Tendency variable, '{tend_stdname}'"
+                    errstr += f", incompatible with associated state variable '{tend_const_name}'"
+                    errstr += f". Reason: '{compat.incompat_reason}'"
+                    raise ParseSyntaxError(errstr, token=tend_stdname,
+                                           context=tendency_variable.context)
+                # end if
+            # end if
+        # end for
+        if not found:
+            # error because we couldn't find the associated constituent
+            errstr = f"No associated state variable for tendency variable, '{tend_stdname}'"
+            raise ParseSyntaxError(errstr, token=tend_stdname,
+                                   context=tendency_variable.context)
+        # end if
+    # end for
     # Parse this table using a fake filename
     parse_obj = ParseObject(f"{host_model.name}_constituent_mod.meta",
                             ddt_mdata)

scripts/metavar.py

@@ -229,7 +229,9 @@ class Var:
                                             optional_in=True, default_in=False),
                            VariableProperty('molar_mass', float,
                                             optional_in=True, default_in=0.0,
-                                            check_fn_in=check_molar_mass)]
+                                            check_fn_in=check_molar_mass),
+                           VariableProperty('constituent', bool,
+                                            optional_in=True, default_in=False)]
 
     __constituent_prop_dict = {x.name : x for x in __constituent_props}
 
@@ -372,7 +374,7 @@ def __init__(self, prop_dict, source, run_env, context=None,
                                    context=self.context) from cperr
         # end try
 
-    def compatible(self, other, run_env):
+    def compatible(self, other, run_env, is_tend=False):
         """Return a VarCompatObj object which describes the equivalence,
         compatibility, or incompatibility between <self> and <other>.
         """
@@ -395,7 +397,7 @@ def compatible(self, other, run_env):
         compat = VarCompatObj(sstd_name, stype, skind, sunits, sdims, sloc_name, stopp,
                               ostd_name, otype, okind, ounits, odims, oloc_name, otopp,
                               run_env,
-                              v1_context=self.context, v2_context=other.context)
+                              v1_context=self.context, v2_context=other.context, is_tend=is_tend)
         if (not compat) and (run_env.logger is not None):
             incompat_str = compat.incompat_reason
             if incompat_str is not None:

scripts/suite_objects.py

@@ -1454,6 +1454,7 @@ def write_var_debug_check(self, var, internal_var, cldicts, outfile, errcode, er
         intent = svar.get_prop_value('intent')
         if intent == 'out' and allocatable:
             return
+        # end if
 
         # Get the condition on which the variable is active
         (conditional, _) = var.conditional(cldicts)

scripts/var_props.py

@@ -857,7 +857,7 @@ class VarCompatObj:
     def __init__(self, var1_stdname, var1_type, var1_kind, var1_units,
                  var1_dims, var1_lname, var1_top, var2_stdname, var2_type, var2_kind,
                  var2_units, var2_dims, var2_lname, var2_top, run_env, v1_context=None,
-                 v2_context=None):
+                 v2_context=None, is_tend=False):
         """Initialize this object with information on the equivalence and/or
            conformability of two variables.
         variable 1 is described by <var1_stdname>, <var1_type>, <var1_kind>,
@@ -866,6 +866,8 @@ def __init__(self, var1_stdname, var1_type, var1_kind, var1_units,
            <var2_units>, <var2_dims>, <var2_lname>, <var2_top>, and <v2_context>.
         <run_env> is the CCPPFrameworkEnv object used here to verify kind
            equivalence or to produce kind transformations.
+        <is_tend> is a flag where, if true, we are validating a tendency variable (var1)
+           against it's equivalent state variable (var2)
         """
         self.__equiv = True   # No transformation required
         self.__compat = True  # Callable with transformation
@@ -881,7 +883,13 @@ def __init__(self, var1_stdname, var1_type, var1_kind, var1_units,
         self.has_vert_transforms = False
         incompat_reason = list()
         # First, check for fatal incompatibilities
-        if var1_stdname != var2_stdname:
+        # If it's a tendency variable, the standard name should be of the
+        #  form "tendency_of_var2_stdname"
+        if is_tend and not var1_stdname.startswith('tendency_of'):
+            self.__equiv = False
+            self.__compat = False
+            incompat_reason.append('not a tendency variable')
+        if not is_tend and  var1_stdname != var2_stdname:
             self.__equiv = False
             self.__compat = False
             incompat_reason.append("standard names")
@@ -944,7 +952,20 @@ def __init__(self, var1_stdname, var1_type, var1_kind, var1_units,
                 var2_units = 'none'
             # end if
             # Check units argument
-            if var1_units != var2_units:
+            if is_tend:
+                # A tendency variable's units should be "<var2_units> s-1"
+                tendency_split_units = var1_units.split('s-1')[0].strip()
+                if tendency_split_units != var2_units:
+                    # We don't currently support unit conversions for tendency variables
+                    emsg = f"\nMismatch tendency variable units '{var1_units}'"
+                    emsg += f" for variable '{var1_stdname}'."
+                    emsg += " No variable transforms supported for tendencies."
+                    emsg += f" Tendency units should be '{var2_units} s-1' to match state variable."
+                    self.__equiv = False
+                    self.__compat = False
+                    incompat_reason.append(emsg)
+                # end if
+            elif var1_units != var2_units:
                 self.__equiv = False
                 # Try to find a set of unit conversions
                 self.__unit_transforms = self._get_unit_convstrs(var1_units,

src/ccpp_constituent_prop_mod.F90

@@ -151,6 +151,7 @@ module ccpp_constituent_prop_mod
       ! These fields are public to allow for efficient (i.e., no copying)
       !   usage even though it breaks object independence
       real(kind_phys), allocatable     :: vars_layer(:,:,:)
+      real(kind_phys), allocatable     :: vars_layer_tend(:,:,:)
       real(kind_phys), allocatable     :: vars_minvalue(:)
       ! An array containing all the constituent metadata
       ! Each element contains a pointer to a constituent from the hash table
@@ -1493,12 +1494,21 @@ subroutine ccp_model_const_data_lock(this, ncols, num_layers, errcode, errmsg)
          call handle_allocate_error(astat, 'vars_layer',                   &
               subname, errcode=errcode, errmsg=errmsg)
          errcode_local = astat
+         if (astat == 0) then
+            allocate(this%vars_layer_tend(ncols, num_layers, this%hash_table%num_values()), &
+                 stat=astat)
+            call handle_allocate_error(astat, 'vars_layer_tend',           &
+                 subname, errcode=errcode, errmsg=errmsg)
+            errcode_local = astat
+         end if
          if (astat == 0) then
             allocate(this%vars_minvalue(this%hash_table%num_values()), stat=astat)
             call handle_allocate_error(astat, 'vars_minvalue',             &
                  subname, errcode=errcode, errmsg=errmsg)
             errcode_local = astat
          end if
+         ! Initialize tendencies to 0
+         this%vars_layer_tend(:,:,:) = 0._kind_phys
          if (errcode_local == 0) then
             this%num_layers = num_layers
             do index = 1, this%hash_table%num_values()
@@ -1549,6 +1559,9 @@ subroutine ccp_model_const_reset(this, clear_hash_table)
       if (allocated(this%vars_minvalue)) then
          deallocate(this%vars_minvalue)
       end if
+      if (allocated(this%vars_layer_tend)) then
+         deallocate(this%vars_layer_tend)
+      end if
       if (allocated(this%const_metadata)) then
          if (clear_table) then
             do index = 1, size(this%const_metadata, 1)

src/ccpp_constituent_prop_mod.meta

@@ -45,6 +45,12 @@
   state_variable = true
   dimensions = (horizontal_dimension, vertical_layer_dimension, number_of_ccpp_constituents)
   type = real | kind = kind_phys
+[ vars_layer_tend ]
+  standard_name = ccpp_constituent_tendencies
+  long_name = Array of constituent tendencies managed by CCPP Framework
+  units = none
+  dimensions = (horizontal_dimension, vertical_layer_dimension, number_of_ccpp_constituents)
+  type = real | kind = kind_phys
 [ const_metadata ]
   standard_name = ccpp_constituent_properties
   units = None

test/advection_test/apply_constituent_tendencies.F90

@@ -0,0 +1,39 @@
+module apply_constituent_tendencies
+
+  use ccpp_kinds, only: kind_phys
+
+  implicit none
+  private
+
+  public :: apply_constituent_tendencies_run
+
+CONTAINS
+
+   !> \section arg_table_apply_constituent_tendencies_run Argument Table
+   !!! \htmlinclude apply_constituent_tendencies_run.html
+   subroutine apply_constituent_tendencies_run(const_tend, const, errcode, errmsg)
+      ! Dummy arguments
+      real(kind_phys),    intent(inout) :: const_tend(:,:,:)  ! constituent tendency array
+      real(kind_phys),    intent(inout) :: const(:,:,:)       ! constituent state array
+      integer,            intent(out)   :: errcode
+      character(len=512), intent(out)   :: errmsg
+
+      ! Local variables
+      integer :: klev, jcnst, icol
+
+      errcode = 0
+      errmsg = ''
+
+      do icol = 1, size(const_tend, 1)
+         do klev = 1, size(const_tend, 2)
+            do jcnst = 1, size(const_tend, 3)
+               const(icol, klev, jcnst) = const(icol, klev, jcnst) + const_tend(icol, klev, jcnst)
+            end do
+         end do
+      end do
+
+      const_tend = 0._kind_phys
+
+   end subroutine apply_constituent_tendencies_run
+
+end module apply_constituent_tendencies

test/advection_test/apply_constituent_tendencies.meta

@@ -0,0 +1,36 @@
+#####################################################################
+[ccpp-table-properties]
+  name = apply_constituent_tendencies
+  type = scheme
+[ccpp-arg-table]
+  name = apply_constituent_tendencies_run
+  type = scheme
+[ const_tend ]
+  standard_name = ccpp_constituent_tendencies
+  long_name = ccpp constituent tendencies
+  units = none
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+  intent = inout
+[ const ]
+  standard_name = ccpp_constituents
+  long_name = ccpp constituents
+  units = none
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+  intent = inout
+[ errcode ]
+  standard_name = ccpp_error_code
+  long_name = Error flag for error handling in CCPP
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+[ errmsg ]
+  standard_name = ccpp_error_message
+  long_name = Error message for error handling in CCPP
+  units = none
+  type = character | kind = len=512
+  dimensions = ()
+  intent = out
+#########################################################