Add constituent tendency capability

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, var in self.items():
+            if 'tendency_of' not in std_name:
+                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 'tendency_of' in std_name:
+                # 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

@@ -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 = "index_of_{}".format(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,23 @@ 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
+    errmsg, errflg = check_tendency_variables(tend_vars, const_vars)
+    if errflg != 0:
+        # There is at least one invalid tendency variable - this is a fatal error
+        raise CCPPError(errmsg)
+    # end if
     # Parse this table using a fake filename
     parse_obj = ParseObject(f"{host_model.name}_constituent_mod.meta",
                             ddt_mdata)
@@ -440,6 +467,117 @@ def add_constituent_vars(cap, host_model, suite_list, run_env):
 
     return const_dict
 
+###############################################################################
+def check_tendency_variables(tend_vars, const_vars):
+###############################################################################
+    """Return an error flag and relevant error message if there are mismatches
+    between the tendency variable and the constituent variable (standard name,
+    units)
+    >>> tendency_vars = [Var({'local_name':'qv_tend', 'standard_name':'tendency_of_water_vapor', 'units':'kg kg-1 s-1', 'dimensions':'(horizontal_dimension, vertical_layer_dimension)', 'type':'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV)]
+    >>> const_vars = [Var({'local_name':'qv', 'standard_name':'water_vapor', 'units': 'kg kg-1', 'dimensions': '(horizontal_dimension, vertical_layer_dimension)', 'type': 'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV)]
+    >>> check_tendency_variables(tendency_vars, const_vars)
+    ('', 0)
+    >>> tendency_vars.append(Var({'local_name':'qc_tend', 'standard_name':'tendency_of_cloud_liquid_water_mixing_ratio', 'units':'kg kg-1 s-1', 'dimensions':'(horizontal_dimension, vertical_layer_dimension)', 'type':'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV))
+    >>> const_vars.append(Var({'local_name':'qc', 'standard_name':'cloud_liquid_water_mixing_ratio', 'units': 'kg kg-1', 'dimensions': '(horizontal_dimension, vertical_layer_dimension)', 'type': 'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV))
+    >>> check_tendency_variables(tendency_vars, const_vars)
+    ('', 0)
+    >>> tendency_vars = [Var({'local_name':'qv_tend', 'standard_name':'tendency_of_water_vapor', 'units':'mol m-3 s-1', 'dimensions':'(horizontal_dimension, vertical_layer_dimension)', 'type':'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV)]
+    >>> check_tendency_variables(tendency_vars, const_vars)
+    ("\\nMismatch tendency variable units 'mol m-3 s-1' for variable 'tendency_of_water_vapor'. No variable transforms supported for constituent tendencies. Tendency units should be 'kg kg-1 s-1' to match constituent.", 1)
+    >>> tendency_vars.append(Var({'local_name':'qc_tend', 'standard_name':'tendency_of_cloud_liquid_water_mixing_ratio', 'units':'kg kg-1 s-1', 'dimensions':'(horizontal_dimension, vertical_interface_dimension)', 'type':'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV))
+    >>> check_tendency_variables(tendency_vars, const_vars)
+    ("\\nMismatch tendency variable units 'mol m-3 s-1' for variable 'tendency_of_water_vapor'. No variable transforms supported for constituent tendencies. Tendency units should be 'kg kg-1 s-1' to match constituent.\\nMismatch tendency variable dimension for dim 2: 'vertical_interface_dimension' for variable 'tendency_of_cloud_liquid_water_mixing_ratio'. Dimension should match constituent dimension of 'vertical_layer_dimension'", 1)
+    >>> tendency_vars = [Var({'local_name':'qv_tend', 'standard_name':'tendency_of_water_vapor', 'units':'kg kg-1 s-1', 'dimensions':'(horizontal_dimension, vertical_layer_dimension)', 'type':'character', 'kind':'len=32'}, _API_SOURCE, _API_DUMMY_RUN_ENV)]
+    >>> check_tendency_variables(tendency_vars, const_vars)
+    ("\\nMismatch tendency variable type 'character' for variable 'tendency_of_water_vapor'. Type should match constituent type of 'real'\\nMismatch tendency variable kind 'len=32' for variable 'tendency_of_water_vapor'. Kind should match constituent kind of 'kind_phys'", 1)
+    >>> tendency_vars = [Var({'local_name':'qv_tend', 'standard_name':'tendency_of_water_vapor_mixing_ratio', 'units':'kg kg-1 s-1', 'dimensions':'(horizontal_dimension, vertical_layer_dimension)', 'type':'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV)]
+    >>> tendency_vars.append(Var({'local_name':'qc_tend', 'standard_name':'tendency_of_cloud_liquid', 'units':'kg kg-1 s-1', 'dimensions':'(horizontal_dimension, vertical_layer_dimension)', 'type':'real', 'kind':'kind_phys'}, _API_SOURCE, _API_DUMMY_RUN_ENV))
+    >>> check_tendency_variables(tendency_vars, const_vars)
+    ("\\nInvalid tendency variable 'tendency_of_water_vapor_mixing_ratio'. All constituent tendency variable standard names must be of the form 'tendency_of_<constituent_stdname>'\\nInvalid tendency variable 'tendency_of_cloud_liquid'. All constituent tendency variable standard names must be of the form 'tendency_of_<constituent_stdname>'", 1)
+    """
+    errflg = 0
+    errmsg = ''
+
+    for tend_var in tend_vars:
+        prop_error = False
+        valid_var = False
+        tend_stdname = tend_var.get_prop_value('standard_name')
+        tend_units = tend_var.get_prop_value('units')
+        tend_stdname_split = tend_stdname.split('tendency_of_')[1]
+        tend_dimensions = tend_var.get_dimensions()
+        tend_type = tend_var.get_prop_value('type')
+        tend_kind = tend_var.get_prop_value('kind')
+        for const_var in const_vars:
+            const_stdname = const_var.get_prop_value('standard_name')
+            const_units = const_var.get_prop_value('units')
+            const_dimensions = const_var.get_dimensions()
+            const_type = const_var.get_prop_value('type')
+            const_kind = const_var.get_prop_value('kind')
+            if tend_stdname_split == const_stdname:
+                # We found a match! Check units
+                if tend_units.split('s-1')[0].strip() != const_units:
+                    errmsg += f"\nMismatch tendency variable units '{tend_units}'"
+                    errmsg += f" for variable '{tend_stdname}'."
+                    errmsg += f" No variable transforms supported for constituent tendencies."
+                    errmsg += f" Tendency units should be '{const_units} s-1' to match constituent."
+                    prop_error = True
+                # end if
+                # Check the dimensions
+                for index, dim in enumerate(tend_dimensions):
+                    # dimension replacement for horizontal_loop_extent
+                    if dim == 'horizontal_loop_extent':
+                        tend_dim = 'horizontal_dimension'
+                    else:
+                        tend_dim = dim
+                    # end if
+                    if const_dimensions[index] == 'horizontal_loop_extent':
+                        const_dim = 'horizontal_dimension'
+                    else:
+                        const_dim = const_dimensions[index]
+                    # end if
+                    if tend_dim != const_dim:
+                        errmsg += f"\nMismatch tendency variable dimension for dim {index + 1}: '{tend_dim}'"
+                        errmsg += f" for variable '{tend_stdname}'."
+                        errmsg += f" Dimension should match constituent dimension of '{const_dim}'"
+                        prop_error = True
+                    # end if
+                # end for
+                # Check the type
+                if const_type != tend_type:
+                    errmsg += f"\nMismatch tendency variable type '{tend_type}'"
+                    errmsg += f" for variable '{tend_stdname}'."
+                    errmsg += f" Type should match constituent type of '{const_type}'"
+                    prop_error = True
+                # end if
+                # Check the kind
+                if const_kind != tend_kind:
+                    errmsg += f"\nMismatch tendency variable kind '{tend_kind}'"
+                    errmsg += f" for variable '{tend_stdname}'."
+                    errmsg += f" Kind should match constituent kind of '{const_kind}'"
+                    prop_error = True
+                # end if
+                if prop_error == 0:
+                    valid_var = True
+                    exit
+                # end if
+            # end if
+        # end for
+        if not valid_var and not prop_error:
+            # Tendency standard name doesn't match an existing constituent
+            errmsg += f"\nInvalid tendency variable '{tend_stdname}'."
+            errmsg += " All constituent tendency variable standard names must "
+            errmsg += "be of the form 'tendency_of_<constituent_stdname>'"
+            errflg = 1
+        elif prop_error:
+            errflg = 1
+        # end if
+    # end for
+    if errflg != 0:
+        return errmsg, errflg
+    # end if
+
+    return errmsg, errflg
+
 ###############################################################################
 def suite_part_call_list(host_model, const_dict, suite_part, subst_loop_vars):
 ###############################################################################

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}
 

scripts/suite_objects.py

@@ -1433,6 +1433,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)

src/ccpp_constituent_prop_mod.F90

@@ -148,6 +148,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
@@ -1465,12 +1466,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()
@@ -1521,6 +1531,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

@@ -34,6 +34,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