diff --git a/cime_config/buildlib b/cime_config/buildlib index 214e6182..33f1d3ec 100755 --- a/cime_config/buildlib +++ b/cime_config/buildlib @@ -189,9 +189,9 @@ def _setup_mpas(case: Case) -> None: shutil.copytree(mpas_dycore_src_root, mpas_dycore_bld_root, copy_function=_copy2_as_needed, dirs_exist_ok=True) shutil.move(os.path.join(mpas_dycore_bld_root, "Makefile"), os.path.join(mpas_dycore_bld_root, "Makefile.CESM")) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "driver", "dyn_mpas_subdriver.F90")), os.path.join(mpas_dycore_bld_root, "driver")) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile")), mpas_dycore_bld_root) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile.in.CESM")), mpas_dycore_bld_root) + + shutil.copytree(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "assets")), mpas_dycore_bld_root, copy_function=_copy2_as_needed, dirs_exist_ok=True) + shutil.copytree(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "driver")), os.path.join(mpas_dycore_bld_root, "driver"), copy_function=_copy2_as_needed, dirs_exist_ok=True) def _copy2_as_needed(src: str, dst: str) -> None: """ diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index fce70728..3054ebb6 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -149,8 +149,16 @@ ${DIN_LOC_ROOT}/atm/waccm/ic/waccm5_1850_ne30np4_L70_0001-01-11-00000_c151217.nc ${DIN_LOC_ROOT}/atm/waccm/ic/fw2000_ne30np4_L70_c181221.nc ${DIN_LOC_ROOT}/atm/cam/inic/gaus/cami_0000-09-01_64x128_L30_c031210.nc - ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L32_notopo_coords_c201216.nc - ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L32_notopo_coords_c201125.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa30_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L58_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L58_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60_L58_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L93_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L93_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60_L93_notopo_coords_c240814.nc diff --git a/src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch b/src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch new file mode 100644 index 00000000..c0d4e1e5 --- /dev/null +++ b/src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch @@ -0,0 +1,190 @@ +From 0000000000000000000000000000000000000000 Mon Sep 17 00:00:00 2001 +From: Kuan-Chih Wang +Date: Mon, 1 Apr 2024 16:46:50 -0600 +Subject: [PATCH] Prefix all MPAS namelist group and option names + +The following transformations are performed for each MPAS namelist group and +option name: + +1. Leading `config_` is removed recursively from the name. Case-insensitive. +2. Leading `mpas_` is removed recursively from the name. Case-insensitive. +3. Prepend `mpas_` to the name. + +As a result, it is easier to distinguish MPAS namelist groups and options from +CAM-SIMA ones. Note that only namelist I/O is affected. Internally, MPAS still +refers to its namelist options by their original names due to compatibility reasons. + +This downstream patch is maintained by CAM-SIMA for its particular use case of MPAS. +--- + src/tools/registry/gen_inc.c | 80 +++++++++++++++++++++++++++++++++++- + src/tools/registry/gen_inc.h | 4 ++ + 2 files changed, 82 insertions(+), 2 deletions(-) + +diff --git a/src/tools/registry/gen_inc.c b/src/tools/registry/gen_inc.c +index 94f5f714..95b2502c 100644 +--- a/src/tools/registry/gen_inc.c ++++ b/src/tools/registry/gen_inc.c +@@ -15,6 +15,10 @@ + #include "fortprintf.h" + #include "utility.h" + ++#ifdef MPAS_CAM_DYCORE ++#include ++#endif ++ + void process_core_macro(const char *macro, const char *val, va_list ap); + void process_domain_macro(const char *macro, const char *val, va_list ap); + +@@ -696,8 +700,18 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + ezxml_t nmlrecs_xml, nmlopt_xml; + + const char *const_core; +- const char *nmlrecname, *nmlrecindef, *nmlrecinsub; +- const char *nmloptname, *nmlopttype, *nmloptval, *nmloptunits, *nmloptdesc, *nmloptposvals, *nmloptindef; ++#ifdef MPAS_CAM_DYCORE ++ const char *orinmlrecname; ++ const char *orinmloptname; ++ // Fortran variable names have a length limit of 63 characters. +1 for the terminating null character. ++ char nmlrecname[64]; ++ char nmloptname[64]; ++#else ++ const char *nmlrecname; ++ const char *nmloptname; ++#endif ++ const char *nmlrecindef, *nmlrecinsub; ++ const char *nmlopttype, *nmloptval, *nmloptunits, *nmloptdesc, *nmloptposvals, *nmloptindef; + + char pool_name[1024]; + char core_string[1024]; +@@ -743,7 +757,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + + // Parse Namelist Records + for (nmlrecs_xml = ezxml_child(registry, "nml_record"); nmlrecs_xml; nmlrecs_xml = nmlrecs_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmlrecname = ezxml_attr(nmlrecs_xml, "name"); ++ transform_name(nmlrecname, sizeof(nmlrecname), orinmlrecname); ++#else + nmlrecname = ezxml_attr(nmlrecs_xml, "name"); ++#endif + nmlrecindef = ezxml_attr(nmlrecs_xml, "in_defaults"); + nmlrecinsub = ezxml_attr(nmlrecs_xml, "in_subpool"); + +@@ -777,7 +796,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + + // Define variable definitions prior to reading the namelist in. + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + nmlopttype = ezxml_attr(nmlopt_xml, "type"); + nmloptval = ezxml_attr(nmlopt_xml, "default_value"); + nmloptunits = ezxml_attr(nmlopt_xml, "units"); +@@ -809,7 +833,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + // Define the namelist block, to read the namelist record in. + fortprintf(fd, " namelist /%s/ &\n", nmlrecname); + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + if(nmlopt_xml->next){ + fortprintf(fd, " %s, &\n", nmloptname); + } else { +@@ -840,7 +869,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + // Define broadcast calls for namelist values. + fortprintf(fd, " if (ierr <= 0) then\n"); + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + nmlopttype = ezxml_attr(nmlopt_xml, "type"); + + if(strncmp(nmlopttype, "real", 1024) == 0){ +@@ -858,7 +892,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + fortprintf(fd, " call mpas_log_write(' The following values will be used for variables in this record:')\n"); + fortprintf(fd, " call mpas_log_write(' ')\n"); + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + nmlopttype = ezxml_attr(nmlopt_xml, "type"); + + if (strncmp(nmlopttype, "character", 1024) == 0) { +@@ -885,10 +924,18 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + fortprintf(fd, "\n"); + + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++ ++ fortprintf(fd, " call mpas_pool_add_config(%s, '%s', %s)\n", pool_name, orinmloptname, nmloptname); ++ fortprintf(fcg, " call mpas_pool_get_config(configPool, '%s', %s)\n", orinmloptname, nmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); + + fortprintf(fd, " call mpas_pool_add_config(%s, '%s', %s)\n", pool_name, nmloptname, nmloptname); + fortprintf(fcg, " call mpas_pool_get_config(configPool, '%s', %s)\n", nmloptname, nmloptname); ++#endif + } + fortprintf(fd, "\n"); + fortprintf(fcg, "\n"); +@@ -2532,3 +2579,32 @@ int parse_structs_from_registry(ezxml_t registry)/*{{{*/ + + return 0; + }/*}}}*/ ++ ++ ++#ifdef MPAS_CAM_DYCORE ++// Perform transformations for namelist group and option names. ++void transform_name(char *new_name, const size_t new_name_size, const char *old_name) { ++ const char *const new_prefix = "mpas_"; ++ const char *const old_prefix = "config_"; ++ size_t size = 0; ++ ++ if (!new_name || !old_name || new_name_size == 0) return; ++ ++ // Remove all leading whitespaces by moving pointer forward. ++ while (*old_name != '\0' && isspace((unsigned char) *old_name)) old_name++; ++ ++ // Remove all leading `config_` by moving pointer forward. ++ while (strncasecmp(old_name, old_prefix, strlen(old_prefix)) == 0) old_name += strlen(old_prefix); ++ ++ // Remove all leading `mpas_` by moving pointer forward. ++ while (strncasecmp(old_name, new_prefix, strlen(new_prefix)) == 0) old_name += strlen(new_prefix); ++ ++ *new_name = '\0'; ++ size = snprintf(NULL, 0, "%s%s", new_prefix, old_name) + 1; ++ snprintf(new_name, size > new_name_size ? new_name_size : size, "%s%s", new_prefix, old_name); ++ ++ // Remove all trailing whitespaces by zeroing (nulling) out. ++ new_name += strlen(new_name) - 1; ++ while (*new_name != '\0' && isspace((unsigned char) *new_name)) *new_name-- = '\0'; ++} ++#endif +diff --git a/src/tools/registry/gen_inc.h b/src/tools/registry/gen_inc.h +index fc94e78b..69a76ace 100644 +--- a/src/tools/registry/gen_inc.h ++++ b/src/tools/registry/gen_inc.h +@@ -38,3 +38,7 @@ int push_attributes(ezxml_t currentPosition); + int merge_structs_and_var_arrays(ezxml_t currentPosition); + int merge_streams(ezxml_t registry); + int parse_structs_from_registry(ezxml_t registry); ++ ++#ifdef MPAS_CAM_DYCORE ++void transform_name(char *new_name, const size_t new_name_size, const char *old_name); ++#endif +-- +2.43.0 + diff --git a/src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch b/src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch new file mode 100644 index 00000000..4f664d88 --- /dev/null +++ b/src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch @@ -0,0 +1,42 @@ +From 0000000000000000000000000000000000000000 Mon Sep 17 00:00:00 2001 +From: Kuan-Chih Wang +Date: Thu, 1 Aug 2024 17:09:58 -0600 +Subject: [PATCH] Disable physics for MPAS dycore-only build + +When building MPAS as a dycore, all physics-related components are disabled. +This build configuration is usually used by CAM/CAM-SIMA, and can be achieved by +defining the `MPAS_CAM_DYCORE` macro in `CPPFLAGS`. + +The `PHYSICS` variable controls whether physics are enabled in MPAS, but its logic +is decoupled from the `MPAS_CAM_DYCORE` macro. Disabling physics in MPAS currently +requires manual interventions. + +Therefore, automatically disable physics when the `MPAS_CAM_DYCORE` macro is found +in `CPPFLAGS`. + +This downstream patch is maintained by CAM-SIMA for its particular use case of MPAS. +--- + src/core_atmosphere/Makefile | 7 +++++-- + 1 file changed, 5 insertions(+), 2 deletions(-) + +diff --git a/src/core_atmosphere/Makefile b/src/core_atmosphere/Makefile +index 8d9f4f1a..cac8255e 100644 +--- a/src/core_atmosphere/Makefile ++++ b/src/core_atmosphere/Makefile +@@ -4,8 +4,11 @@ + # To build a dycore-only MPAS-Atmosphere model, comment-out or delete + # the definition of PHYSICS, below + # +-PHYSICS=-DDO_PHYSICS +- ++# If MPAS_CAM_DYCORE is found in CPPFLAGS, PHYSICS will become undefined automatically ++# ++ifeq ($(findstring MPAS_CAM_DYCORE,$(CPPFLAGS)),) ++ PHYSICS = -DDO_PHYSICS ++endif + + ifdef PHYSICS + PHYSCORE = physcore +-- +2.43.0 + diff --git a/src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch b/src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch new file mode 100644 index 00000000..d88d02ac --- /dev/null +++ b/src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch @@ -0,0 +1,94 @@ +From 0000000000000000000000000000000000000000 Mon Sep 17 00:00:00 2001 +From: Kuan-Chih Wang +Date: Wed, 13 Nov 2024 17:19:10 -0700 +Subject: [PATCH] Declare constants at the native precision of MPAS + +When building MPAS as a dycore, certain constants in the `mpas_constants` module are +imported from the `physconst` module, which is a part of CAM/CAM-SIMA. However, +multiple issues arise if the precision of those constants differs from MPAS. + +For example, building MPAS in single precision mode with CAM-SIMA fails due to +multiple occurrences of type mismatch between actual and dummy arguments. + +mpas_geometry_utils.F:885:157: + + 885 | call mpas_log_write('$r', MPAS_LOG_ERR, realArgs=(/mpas_triangle_signed_area_sphere(a,b,c,sphereRadius) - pii/2.0_RKIND*sphereRadius*sphereRadius/)) + | 1 +Error: Type mismatch in argument 'realargs' at (1); passed REAL(8) to REAL(4) + +Here, `pii` is declared by CAM-SIMA to be double precision, and it causes unintended +floating-point promotion in the expression. + +The solution is to ensure that constants in the `mpas_constants` module are declared +at the native precision of MPAS. + +This downstream patch is maintained by CAM-SIMA for its particular use case of MPAS. +--- + src/framework/mpas_constants.F | 44 ++++++++++++++++++++++++++-------- + 1 file changed, 34 insertions(+), 10 deletions(-) + +diff --git a/src/framework/mpas_constants.F b/src/framework/mpas_constants.F +index c98cb810..68164bab 100644 +--- a/src/framework/mpas_constants.F ++++ b/src/framework/mpas_constants.F +@@ -23,16 +23,30 @@ module mpas_constants + use mpas_kind_types + + #ifdef MPAS_CAM_DYCORE +- use physconst, only : pii => pi +- use physconst, only : gravity => gravit +- use physconst, only : omega +- use physconst, only : a => rearth +- use physconst, only : cp => cpair +- use physconst, only : rgas => rair +- use physconst, only : rv => rh2o +- real (kind=RKIND) :: rvord = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived +- real (kind=RKIND) :: cv = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived +- real (kind=RKIND) :: cvpm = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ use physconst, only: external_pii => pi ++ use physconst, only: external_a => rearth ++ use physconst, only: external_omega => omega ++ use physconst, only: external_gravity => gravit ++ use physconst, only: external_rgas => rair ++ use physconst, only: external_rv => rh2o ++ use physconst, only: external_cp => cpair ++ private :: external_pii ++ private :: external_a ++ private :: external_omega ++ private :: external_gravity ++ private :: external_rgas ++ private :: external_rv ++ private :: external_cp ++ real (kind=RKIND), protected :: pii = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: a = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: omega = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: gravity = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: rgas = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: rv = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: cp = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: rvord = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: cv = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: cvpm = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived + #else + real (kind=RKIND), parameter :: pii = 3.141592653589793_RKIND !< Constant: Pi + real (kind=RKIND), parameter :: a = 6371229.0_RKIND !< Constant: Spherical Earth radius [m] +@@ -77,6 +91,16 @@ module mpas_constants + implicit none + + #ifdef MPAS_CAM_DYCORE ++ ! Convert external constants to the native precision of MPAS (i.e., `RKIND`). ++ ++ pii = real(external_pii, RKIND) ++ a = real(external_a, RKIND) ++ omega = real(external_omega, RKIND) ++ gravity = real(external_gravity, RKIND) ++ rgas = real(external_rgas, RKIND) ++ rv = real(external_rv, RKIND) ++ cp = real(external_cp, RKIND) ++ + ! + ! In the case of CAM-MPAS, rgas may depend on a CAM namelist option, + ! so physical constants that depend on rgas must be computed here after +-- +2.43.0 + diff --git a/src/dynamics/mpas/Makefile b/src/dynamics/mpas/assets/Makefile similarity index 100% rename from src/dynamics/mpas/Makefile rename to src/dynamics/mpas/assets/Makefile diff --git a/src/dynamics/mpas/Makefile.in.CESM b/src/dynamics/mpas/assets/Makefile.in.CESM similarity index 86% rename from src/dynamics/mpas/Makefile.in.CESM rename to src/dynamics/mpas/assets/Makefile.in.CESM index 8e4ba049..4b1c5708 100644 --- a/src/dynamics/mpas/Makefile.in.CESM +++ b/src/dynamics/mpas/assets/Makefile.in.CESM @@ -49,6 +49,10 @@ export CPPFLAGS := -D_MPI \ ifneq ($(strip $(PIODEF)),) export CPPFLAGS += $(strip $(PIODEF)) endif +# Uncomment below for MPI Fortran 2008 interface support. There is currently no corresponding option in CIME to enable this. +# export CPPFLAGS += -DMPAS_USE_MPI_F08 +# Uncomment below for single precision mode support. There is currently no corresponding option in CIME to enable this. +# export CPPFLAGS += -DSINGLE_PRECISION export LINKER := $(strip $(FC)) export SCC := $(strip $(CC)) export SCXX := $(strip $(CXX)) @@ -71,7 +75,17 @@ all: @echo ' `make libmpas-clean ESM="CESM" LIBROOT="..." SRCROOT="..."`' .PHONY: libmpas-prepare -libmpas-prepare: libmpas-archiver-script.txt libmpas-no-physics libmpas-prefix-namelist-groups libmpas-preview +libmpas-prepare: libmpas-apply-patch libmpas-archiver-script.txt libmpas-preview + +# Apply patches. +.PHONY: libmpas-apply-patch +libmpas-apply-patch: + @for file in *.patch; do \ + if git apply --check -p2 "$${file}"; then \ + echo "Applying $${file}"; \ + git apply -p2 "$${file}"; \ + fi; \ + done # Combine multiple static libraries into `libmpas.a` via archiver/MRI script. This requires GNU or GNU-like archiver (`ar`) program. libmpas-archiver-script.txt: @@ -83,19 +97,8 @@ libmpas-archiver-script.txt: @echo "save" >> $(@) @echo "end" >> $(@) -# Do not use built-in MPAS/WRF physics. -.PHONY: libmpas-no-physics -libmpas-no-physics: - @sed -E -i -e "s/^ *PHYSICS=.+$$/PHYSICS=/g" core_atmosphere/Makefile - -# Prefix `mpas_` to MPAS namelist groups to avoid naming collisions. -.PHONY: libmpas-prefix-namelist-groups -libmpas-prefix-namelist-groups: - @sed -E -i -e "s/(^ *< *nml_record.+name=)\"(mpas_)?(\w+)\"/\1\"mpas_\3\"/g" core_atmosphere/Registry.xml - @sed -E -i -e "s/(^ *< *nml_record.+name=)\"(mpas_)?(\w+)\"/\1\"mpas_\3\"/g" core_atmosphere/diagnostics/Registry_*.xml - .PHONY: libmpas-preview -libmpas-preview: +libmpas-preview: libmpas-apply-patch libmpas-archiver-script.txt @echo "Previewing build options for $(LIBROOT)/libmpas.a:" @echo "AR = $(AR)" @echo "ARFLAGS = $(ARFLAGS)" diff --git a/src/dynamics/mpas/assets/generate_namelist_definition.py b/src/dynamics/mpas/assets/generate_namelist_definition.py new file mode 100755 index 00000000..12dd0c54 --- /dev/null +++ b/src/dynamics/mpas/assets/generate_namelist_definition.py @@ -0,0 +1,301 @@ +#!/usr/bin/env python3 + +''' +Generate XML namelist definition file for MPAS dynamical core in CAM-SIMA. +''' + +import argparse +import textwrap +import xml.etree.ElementTree as ET + +# These namelist groups are irrelevant to CAM-SIMA for its particular use case of MPAS. +# Hide them to prevent users from inadvertently modifying them. +EXCLUDED_NAMELIST_GROUP = [ + 'assimilation', + 'iau', + 'limited_area', + 'physics', + 'restart' +] +# These namelist options are forcefully controlled by CAM-SIMA at run-time. +# Hide them to prevent users from inadvertently modifying them. +EXCLUDED_NAMELIST_OPTION = [ + 'config_calendar_type', + 'config_do_restart', + 'config_run_duration', + 'config_start_time', + 'config_stop_time' +] +# List all overridden namelist options below. +# `OVERRIDDEN_NAMELIST_OPTION` is a dictionary with `str` as keys and `list[str]` as values. +OVERRIDDEN_NAMELIST_OPTION = { + 'mpas_block_decomp_file_prefix': [ + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa30.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa12.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15-3.graph.info.part.' + ], + 'mpas_coef_3rd_order': [ + '1.0' + ], + 'mpas_dt': [ + '1800.0', + '900.0', + '450.0', + '225.0' + ], + 'mpas_len_disp': [ + '480000.0', + '120000.0', + '60000.0', + '30000.0' + ], + 'mpas_number_cam_damping_levels': [ + '0' + ], + 'mpas_number_rayleigh_damp_u_levels': [ + '5' + ], + 'mpas_print_detailed_minmax_vel': [ + '.true.' + ], + 'mpas_rayleigh_damp_u': [ + '.true.' + ] +} + +INDENT_PER_LEVEL = ' ' * 4 +NEW_PREFIX = 'mpas_' +OLD_PREFIX = 'config_' + +def parse_argument() -> argparse.Namespace: + ''' + Parse command line arguments. + ''' + + parser = argparse.ArgumentParser( + description='Generate XML namelist definition file for MPAS dynamical core in CAM-SIMA.' + ) + + parser.add_argument( + '-r', '--registry', + default='Registry.xml', + type=str, + required=False, + help='XML MPAS registry file.', + dest='reg_xml' + ) + parser.add_argument( + '-n', '--namelist', + default='Namelist.xml', + type=str, + required=False, + help='XML CAM-SIMA namelist definition file.', + dest='nml_xml' + ) + parser.add_argument( + '-s', '--schema', + default=None, + type=str, + required=False, + help='XML schema for CAM-SIMA namelist definition file.', + dest='nml_xsd' + ) + + argument = parser.parse_args() + + return argument + +def parse_xml(xml_file: str) -> ET.ElementTree: + ''' + Parse XML file into element tree. + ''' + + xml_et = ET.parse(xml_file) + + return xml_et + +def validate_xml(xml_file: str, xsd_file: str) -> bool: + ''' + Validate XML file against XSD file. + ''' + + # Only import `xmlschema` if XML validation is requested. Shush pylint about it. + import xmlschema # pylint: disable=import-outside-toplevel + + xml_schema = xmlschema.XMLSchema(xsd_file) + + return xml_schema.is_valid(xml_file) + +def transform_name(name: str) -> str: + ''' + Change prefix of namelist option/group name. + ''' + + while name.startswith(OLD_PREFIX): + name = name[len(OLD_PREFIX):] + + while name.startswith(NEW_PREFIX): + name = name[len(NEW_PREFIX):] + + name = NEW_PREFIX + name + + return name + +def translate_element_tree(reg_xml_et: ET.ElementTree) -> ET.ElementTree: + ''' + Translate MPAS registry into namelist definition. + ''' + + # `entry_id_pg` is the root element in namelist definition. + entry_id_pg_element = ET.Element('entry_id_pg', {'version': '0.1'}) + + comment_element = ET.Comment( + '\n' + + INDENT_PER_LEVEL * 2 + 'MPAS dycore' + '\n' + + '\n' + + INDENT_PER_LEVEL * 2 + 'Note to developers/maintainers:' + '\n' + + INDENT_PER_LEVEL * 2 + 'This file is auto-generated from the MPAS registry. Do not edit directly.' + '\n' + + INDENT_PER_LEVEL * 2 + 'Instead, use the Python script at `src/dynamics/mpas/assets/generate_namelist_definition.py`.' + '\n' + + INDENT_PER_LEVEL + ) + entry_id_pg_element.append(comment_element) + + for namelist_group in reg_xml_et.findall('nml_record'): + if namelist_group.attrib['name'].strip().lower() in EXCLUDED_NAMELIST_GROUP: + continue + + for namelist_option in namelist_group.findall('nml_option'): + if namelist_option.attrib['name'].strip().lower() in EXCLUDED_NAMELIST_OPTION: + continue + + # The `entry_id_pg` root element contains many `entry` elements. + # Each `entry` element describes a namelist option, indicated by its `id` attribute. + entry_element = ET.SubElement(entry_id_pg_element, 'entry', {'id': transform_name(namelist_option.attrib['name'].strip().lower())}) + + # The `category` element. + category_element = ET.SubElement(entry_element, 'category') + category_element.text = 'mpas' + + # The `desc` element. + desc_text = ' '.join(namelist_option.attrib['description'].strip(' .').split()) + desc_text = desc_text[0].upper() + desc_text[1:] + desc_text = '\n' + textwrap.fill(desc_text, 80, initial_indent=INDENT_PER_LEVEL * 3, subsequent_indent=INDENT_PER_LEVEL * 3) + '\n' + INDENT_PER_LEVEL * 2 + + desc_element = ET.SubElement(entry_element, 'desc') + desc_element.text = desc_text + + # The `group` element. + group_element = ET.SubElement(entry_element, 'group') + group_element.text = transform_name(namelist_group.attrib['name'].strip().lower()) + + # The `type` element. + # The `values` element and its containing `value` element. + type_text = namelist_option.attrib['type'].strip().lower() + value_text = namelist_option.attrib['default_value'] + + # Do some sanitization. + if type_text.startswith('ch'): + type_text = 'char*256' + value_text = value_text.strip() + + if not value_text.isascii(): + raise ValueError('"' + value_text + '" is not ASCII') + elif type_text.startswith('in'): + type_text = 'integer' + value_text = canonicalize_int(value_text) + elif type_text.startswith('lo'): + type_text = 'logical' + value_text = value_text.strip().lower() + + if value_text.startswith(('t', '.t')): + value_text = '.true.' + elif value_text.startswith(('f', '.f')): + value_text = '.false.' + else: + raise ValueError('"' + value_text + '" does not represent a logical') + elif type_text.startswith('re'): + type_text = 'real' + value_text = canonicalize_real(value_text) + else: + raise ValueError('"' + type_text + '" is not a supported type') + + type_element = ET.SubElement(entry_element, 'type') + type_element.text = type_text + + values_element = ET.SubElement(entry_element, 'values') + + if entry_element.attrib['id'] in OVERRIDDEN_NAMELIST_OPTION: + for value_string in OVERRIDDEN_NAMELIST_OPTION[entry_element.attrib['id']]: + value_element = ET.fromstring(value_string) + values_element.append(value_element) + else: + value_element = ET.SubElement(values_element, 'value') + value_element.text = value_text + + # Sort the `entry` elements for result stability except for the comment element at index 0. + entry_id_pg_element[1:] = sorted(entry_id_pg_element.iterfind('entry'), key=lambda entry: entry.get('id')) + nml_xml_et = ET.ElementTree(entry_id_pg_element) + + return nml_xml_et + +def canonicalize_int(s: str) -> str: + ''' + Canonicalize a string that represents an integer. + ''' + + try: + i = int(s) + except Exception as e: + raise ValueError('"' + s + '" does not represent an integer') from e + + return str(i) + +def canonicalize_real(s: str) -> str: + ''' + Canonicalize a string that represents a real. + ''' + + try: + f = float(s) + except Exception as e: + raise ValueError('"' + s + '" does not represent a real') from e + + return str(f) + +def write_element_tree(xml_file: str, xml_et: ET.ElementTree) -> None: + ''' + Write element tree into XML file. + ''' + + ET.indent(xml_et, space=INDENT_PER_LEVEL) + + xml_et.write( + xml_file, + encoding='UTF-8', + xml_declaration=True, + default_namespace=None, + method='xml', + short_empty_elements=False + ) + + # The file written by `ElementTree.write()` contains no newline at end of file. Add it manually. + with open(xml_file, 'a', encoding='utf-8') as nml_xml_file: + nml_xml_file.write('\n') + +if __name__ == '__main__': + arg = parse_argument() + reg_xml_element_tree = parse_xml(arg.reg_xml) + nml_xml_element_tree = translate_element_tree(reg_xml_element_tree) + + write_element_tree(arg.nml_xml, nml_xml_element_tree) + print('Generated ' + arg.nml_xml) + + if arg.nml_xsd is not None: + if validate_xml(arg.nml_xml, arg.nml_xsd): + print('Successfully validated ' + arg.nml_xml + ' against ' + arg.nml_xsd) + else: + print('Failed to validate ' + arg.nml_xml + ' against ' + arg.nml_xsd) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index ca8e4159..b4ec0921 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -8,62 +8,21 @@ module dyn_mpas_subdriver !------------------------------------------------------------------------------- use, intrinsic :: iso_fortran_env, only: output_unit - - ! Modules from external libraries. + ! Module(s) from external libraries. #ifdef MPAS_USE_MPI_F08 use mpi_f08, only: mpi_comm_null, mpi_comm_rank, mpi_success, & mpi_comm_type => mpi_comm, operator(==) #else use mpi, only: mpi_comm_null, mpi_comm_rank, mpi_success #endif - use pio, only: pio_char, pio_int, pio_real, pio_double, & - file_desc_t, iosystem_desc_t, pio_file_is_open, pio_iosystem_is_active, & - pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr - - ! Modules from MPAS. - use atm_core, only: atm_compute_output_diagnostics, atm_do_timestep, atm_mpas_init_block - use atm_core_interface, only: atm_setup_core, atm_setup_domain - use atm_time_integration, only: mpas_atm_dynamics_init - use mpas_atm_dimensions, only: mpas_atm_set_dims - use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group - use mpas_atm_threading, only: mpas_atm_threading_init - use mpas_attlist, only: mpas_modify_att - use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2 - use mpas_constants, only: mpas_constants_compute_derived - use mpas_derived_types, only: core_type, domain_type, & - mpas_pool_type, mpas_pool_field_info_type, & - mpas_pool_character, mpas_pool_real, mpas_pool_integer, & - mpas_stream_type, mpas_stream_noerr, & - mpas_time_type, mpas_timeinterval_type, mpas_now, mpas_start_time, & - mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, & - field0dchar, field1dchar, & - field0dinteger, field1dinteger, field2dinteger, field3dinteger, & - field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal - use mpas_dmpar, only: mpas_dmpar_exch_halo_field, & - mpas_dmpar_max_int, mpas_dmpar_sum_int - use mpas_domain_routines, only: mpas_allocate_domain - use mpas_field_routines, only: mpas_allocate_scratch_field - use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2 - use mpas_io_streams, only: mpas_createstream, mpas_closestream, mpas_streamaddfield, & - mpas_readstream, mpas_writestream, mpas_writestreamatt - use mpas_kind_types, only: rkind, r4kind, r8kind, strkind - use mpas_pool_routines, only: mpas_pool_create_pool, mpas_pool_destroy_pool, mpas_pool_get_subpool, & - mpas_pool_add_config, mpas_pool_get_config, & - mpas_pool_get_array, & - mpas_pool_add_dimension, mpas_pool_get_dimension, & - mpas_pool_get_field, mpas_pool_get_field_info, & - mpas_pool_initialize_time_levels, mpas_pool_shift_time_levels - use mpas_stream_inquiry, only: mpas_stream_inquiry_new_streaminfo - use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex - use mpas_string_utils, only: mpas_string_replace - use mpas_timekeeping, only: mpas_advance_clock, mpas_get_clock_time, mpas_get_time, & - mpas_set_timeinterval, & - operator(+), operator(<) - use mpas_vector_operations, only: mpas_initialize_vectors + ! Module(s) from MPAS. + use mpas_derived_types, only: core_type, domain_type + use mpas_kind_types, only: rkind, strkind implicit none private + public :: mpas_dynamical_core_real_kind public :: mpas_dynamical_core_type abstract interface @@ -75,6 +34,9 @@ subroutine model_error_if(message, file, line) end subroutine model_error_if end interface + !> The native floating-point precision of MPAS dynamical core. + integer, parameter :: mpas_dynamical_core_real_kind = rkind + !> The "class" of MPAS dynamical core. !> Important data structures like states of MPAS dynamical core are encapsulated inside this derived type to prevent misuse. !> Type-bound procedures provide well-defined APIs for CAM-SIMA to interact with MPAS dynamical core. @@ -128,6 +90,7 @@ end subroutine model_error_if procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4 procedure, pass, public :: run => dyn_mpas_run + procedure, pass, public :: final => dyn_mpas_final ! Accessor subroutines for users to access internal states of MPAS dynamical core. @@ -189,6 +152,8 @@ end subroutine model_error_if !> !> Here, it is described as: !> var_info_type(name="xCell", type="real", rank=1) + !> However, note that MPAS treats the "Time" dimension specially. It is implemented as 1-d pointer arrays of + !> custom derived types. For a variable with the "Time" dimension, its rank needs to be subtracted by one. type :: var_info_type private @@ -277,40 +242,54 @@ end subroutine model_error_if ! the `atm_init_coupled_diagnostics` subroutine in MPAS. ! If a variable first appears on the LHS of an equation, it should be in restart. ! If a variable first appears on the RHS of an equation, it should be in input. + ! The remaining ones of interest should be in output. !> This list corresponds to the "input" stream in MPAS registry. !> It consists of variables that are members of the "diag" and "state" struct. !> Only variables that are specific to the "input" stream are included. type(var_info_type), parameter :: input_var_info_list(*) = [ & - var_info_type('Time' , 'real' , 1), & + var_info_type('Time' , 'real' , 0), & var_info_type('initial_time' , 'character' , 0), & - var_info_type('rho' , 'real' , 3), & - var_info_type('rho_base' , 'real' , 3), & + var_info_type('rho' , 'real' , 2), & + var_info_type('rho_base' , 'real' , 2), & var_info_type('scalars' , 'real' , 3), & - var_info_type('theta' , 'real' , 3), & - var_info_type('theta_base' , 'real' , 3), & - var_info_type('u' , 'real' , 3), & - var_info_type('w' , 'real' , 3), & - var_info_type('xtime' , 'character' , 1) & + var_info_type('theta' , 'real' , 2), & + var_info_type('theta_base' , 'real' , 2), & + var_info_type('u' , 'real' , 2), & + var_info_type('w' , 'real' , 2), & + var_info_type('xtime' , 'character' , 0) & ] !> This list corresponds to the "restart" stream in MPAS registry. !> It consists of variables that are members of the "diag" and "state" struct. !> Only variables that are specific to the "restart" stream are included. type(var_info_type), parameter :: restart_var_info_list(*) = [ & - var_info_type('exner' , 'real' , 1), & - var_info_type('exner_base' , 'real' , 1), & - var_info_type('pressure_base' , 'real' , 1), & - var_info_type('pressure_p' , 'real' , 1), & - var_info_type('rho_p' , 'real' , 1), & - var_info_type('rho_zz' , 'real' , 1), & - var_info_type('rtheta_base' , 'real' , 1), & - var_info_type('rtheta_p' , 'real' , 1), & - var_info_type('ru' , 'real' , 1), & - var_info_type('ru_p' , 'real' , 1), & - var_info_type('rw' , 'real' , 1), & - var_info_type('rw_p' , 'real' , 1), & - var_info_type('theta_m' , 'real' , 1) & + var_info_type('exner' , 'real' , 2), & + var_info_type('exner_base' , 'real' , 2), & + var_info_type('pressure_base' , 'real' , 2), & + var_info_type('pressure_p' , 'real' , 2), & + var_info_type('rho_p' , 'real' , 2), & + var_info_type('rho_zz' , 'real' , 2), & + var_info_type('rtheta_base' , 'real' , 2), & + var_info_type('rtheta_p' , 'real' , 2), & + var_info_type('ru' , 'real' , 2), & + var_info_type('ru_p' , 'real' , 2), & + var_info_type('rw' , 'real' , 2), & + var_info_type('rw_p' , 'real' , 2), & + var_info_type('theta_m' , 'real' , 2) & + ] + + !> This list corresponds to the "output" stream in MPAS registry. + !> It consists of variables that are members of the "diag" struct. + !> Only variables that are specific to the "output" stream are included. + type(var_info_type), parameter :: output_var_info_list(*) = [ & + var_info_type('divergence' , 'real' , 2), & + var_info_type('pressure' , 'real' , 2), & + var_info_type('relhum' , 'real' , 2), & + var_info_type('surface_pressure' , 'real' , 1), & + var_info_type('uReconstructMeridional' , 'real' , 2), & + var_info_type('uReconstructZonal' , 'real' , 2), & + var_info_type('vorticity' , 'real' , 2) & ] contains !> Print a debug message with optionally the value(s) of a variable. @@ -461,6 +440,11 @@ end function stringify ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_unit, mpas_log_unit) + ! Module(s) from MPAS. + use atm_core_interface, only: atm_setup_core, atm_setup_domain + use mpas_domain_routines, only: mpas_allocate_domain + use mpas_framework, only: mpas_framework_init_phase1 + class(mpas_dynamical_core_type), intent(inout) :: self #ifdef MPAS_USE_MPI_F08 type(mpi_comm_type), intent(in) :: mpi_comm @@ -619,11 +603,7 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, & subname, __LINE__) end select - call mpas_pool_get_config(self % domain_ptr % configs, 'config_calendar_type', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_calendar_type"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_calendar_type') config_pointer_c = trim(adjustl(mpas_calendar)) call self % debug_print('config_calendar_type = ', [config_pointer_c]) @@ -632,43 +612,27 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, & ! MPAS represents date and time in ISO 8601 format. However, the separator between date and time is `_` ! instead of standard `T`. ! Format in `YYYY-MM-DD_hh:mm:ss` is acceptable. - call mpas_pool_get_config(self % domain_ptr % configs, 'config_start_time', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_start_time"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_start_time') config_pointer_c = stringify(start_date_time(1:3), '-') // '_' // stringify(start_date_time(4:6), ':') call self % debug_print('config_start_time = ', [config_pointer_c]) nullify(config_pointer_c) - call mpas_pool_get_config(self % domain_ptr % configs, 'config_stop_time', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_stop_time"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_stop_time') config_pointer_c = stringify(stop_date_time(1:3), '-') // '_' // stringify(stop_date_time(4:6), ':') call self % debug_print('config_stop_time = ', [config_pointer_c]) nullify(config_pointer_c) ! Format in `DD_hh:mm:ss` is acceptable. - call mpas_pool_get_config(self % domain_ptr % configs, 'config_run_duration', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_run_duration"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_run_duration') config_pointer_c = stringify([run_duration(1)]) // '_' // stringify(run_duration(2:4), ':') call self % debug_print('config_run_duration = ', [config_pointer_c]) nullify(config_pointer_c) ! Reflect current run type to MPAS. - call mpas_pool_get_config(self % domain_ptr % configs, 'config_do_restart', config_pointer_l) - - if (.not. associated(config_pointer_l)) then - call self % model_error('Failed to find config "config_do_restart"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_l, 'cfg', 'config_do_restart') if (initial_run) then ! Run type is initial run. @@ -699,6 +663,12 @@ end subroutine dyn_mpas_read_namelist ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase2(self, pio_iosystem) + ! Module(s) from external libraries. + use pio, only: iosystem_desc_t, pio_iosystem_is_active + ! Module(s) from MPAS. + use mpas_framework, only: mpas_framework_init_phase2 + use mpas_stream_inquiry, only: mpas_stream_inquiry_new_streaminfo + class(mpas_dynamical_core_type), intent(in) :: self type(iosystem_desc_t), pointer, intent(in) :: pio_iosystem @@ -788,6 +758,13 @@ end subroutine dyn_mpas_init_phase2 ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file) + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open + ! Module(s) from MPAS. + use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2 + use mpas_derived_types, only: mpas_io_pnetcdf, mpas_pool_type + use mpas_pool_routines, only: mpas_pool_add_config, mpas_pool_add_dimension, mpas_pool_get_dimension + class(mpas_dynamical_core_type), intent(inout) :: self integer, intent(in) :: number_of_constituents type(file_desc_t), pointer, intent(in) :: pio_file @@ -884,6 +861,10 @@ end subroutine dyn_mpas_init_phase3 ! !------------------------------------------------------------------------------- subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species) + ! Module(s) from MPAS. + use mpas_derived_types, only: field3dreal, mpas_pool_type + use mpas_pool_routines, only: mpas_pool_add_dimension, mpas_pool_get_field + class(mpas_dynamical_core_type), intent(inout) :: self character(*), intent(in) :: constituent_name(:) logical, intent(in) :: is_water_species(:) @@ -1142,6 +1123,14 @@ end subroutine dyn_mpas_define_scalar ! !------------------------------------------------------------------------------- subroutine dyn_mpas_read_write_stream(self, pio_file, stream_mode, stream_name) + ! Module(s) from external libraries. + use pio, only: file_desc_t + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type, mpas_stream_noerr, mpas_stream_type + use mpas_io_streams, only: mpas_closestream, mpas_readstream, mpas_writestream + use mpas_pool_routines, only: mpas_pool_destroy_pool + use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex + class(mpas_dynamical_core_type), intent(in) :: self type(file_desc_t), pointer, intent(in) :: pio_file character(*), intent(in) :: stream_mode @@ -1248,6 +1237,17 @@ end subroutine dyn_mpas_read_write_stream ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file, stream_mode, stream_name) + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open + ! Module(s) from MPAS. + use mpas_derived_types, only: field0dchar, field1dchar, & + field0dinteger, field1dinteger, field2dinteger, field3dinteger, & + field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal, & + mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, & + mpas_pool_type, mpas_stream_noerr, mpas_stream_type + use mpas_io_streams, only: mpas_createstream, mpas_streamaddfield + use mpas_pool_routines, only: mpas_pool_add_config, mpas_pool_create_pool, mpas_pool_get_field + class(mpas_dynamical_core_type), intent(in) :: self type(mpas_pool_type), pointer, intent(out) :: mpas_pool type(mpas_stream_type), pointer, intent(out) :: mpas_stream @@ -1583,6 +1583,9 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file !> Helper subroutine for adding a 0-d stream attribute by calling `mpas_writestreamatt` with error checking. !> (KCW, 2024-03-14) subroutine add_stream_attribute_0d(attribute_name, attribute_value) + ! Module(s) from MPAS. + use mpas_io_streams, only: mpas_writestreamatt + character(*), intent(in) :: attribute_name class(*), intent(in) :: attribute_value @@ -1620,6 +1623,9 @@ end subroutine add_stream_attribute_0d !> Helper subroutine for adding a 1-d stream attribute by calling `mpas_writestreamatt` with error checking. !> (KCW, 2024-03-14) subroutine add_stream_attribute_1d(attribute_name, attribute_value) + ! Module(s) from MPAS. + use mpas_io_streams, only: mpas_writestreamatt + character(*), intent(in) :: attribute_name class(*), intent(in) :: attribute_value(:) @@ -1752,6 +1758,8 @@ pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_l allocate(var_info_list, source=input_var_info_list) case ('restart') allocate(var_info_list, source=restart_var_info_list) + case ('output') + allocate(var_info_list, source=output_var_info_list) case default allocate(var_info_list(0)) @@ -1775,6 +1783,13 @@ pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_l var_info_list_buffer = pack(restart_var_info_list, var_name_list == trim(adjustl(stream_name_fragment))) var_info_list = [var_info_list, var_info_list_buffer] end if + + var_name_list = output_var_info_list % name + + if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then + var_info_list_buffer = pack(output_var_info_list, var_name_list == trim(adjustl(stream_name_fragment))) + var_info_list = [var_info_list, var_info_list_buffer] + end if end select end function parse_stream_name_fragment @@ -1874,6 +1889,17 @@ end function index_unique ! !------------------------------------------------------------------------------- subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compatible, pio_file, var_info) + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open, & + pio_char, pio_int, pio_real, pio_double, & + pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr + ! Module(s) from MPAS. + use mpas_derived_types, only: field0dchar, field1dchar, & + field0dinteger, field1dinteger, field2dinteger, field3dinteger, & + field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal + use mpas_kind_types, only: r4kind, r8kind + use mpas_pool_routines, only: mpas_pool_get_field + class(mpas_dynamical_core_type), intent(in) :: self logical, allocatable, intent(out) :: var_is_present(:) logical, allocatable, intent(out) :: var_is_tkr_compatible(:) @@ -2239,11 +2265,15 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa cycle end if case ('real') + ! When MPAS dynamical core is compiled at single precision, pairing it with double precision input data + ! is not allowed to prevent loss of precision. if (rkind == r4kind .and. vartype /= pio_real) then cycle end if - if (rkind == r8kind .and. vartype /= pio_double) then + ! When MPAS dynamical core is compiled at double precision, pairing it with single and double precision + ! input data is allowed. + if (rkind == r8kind .and. vartype /= pio_real .and. vartype /= pio_double) then cycle end if case default @@ -2285,6 +2315,13 @@ end subroutine dyn_mpas_check_variable_status ! !------------------------------------------------------------------------------- subroutine dyn_mpas_exchange_halo(self, field_name) + ! Module(s) from MPAS. + use mpas_derived_types, only: field1dinteger, field2dinteger, field3dinteger, & + field1dreal, field2dreal, field3dreal, field4dreal, field5dreal, & + mpas_pool_field_info_type, mpas_pool_integer, mpas_pool_real + use mpas_dmpar, only: mpas_dmpar_exch_halo_field + use mpas_pool_routines, only: mpas_pool_get_field, mpas_pool_get_field_info + class(mpas_dynamical_core_type), intent(in) :: self character(*), intent(in) :: field_name @@ -2468,6 +2505,10 @@ end subroutine dyn_mpas_exchange_halo ! !------------------------------------------------------------------------------- subroutine dyn_mpas_compute_unit_vector(self) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_vector_operations, only: mpas_initialize_vectors + class(mpas_dynamical_core_type), intent(in) :: self character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_unit_vector' @@ -2652,6 +2693,21 @@ end subroutine dyn_mpas_compute_edge_wind ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase4(self, coupling_time_interval) + ! Module(s) from MPAS. + use atm_core, only: atm_mpas_init_block + use atm_time_integration, only: mpas_atm_dynamics_init + use mpas_atm_dimensions, only: mpas_atm_set_dims + use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group + use mpas_atm_threading, only: mpas_atm_threading_init + use mpas_attlist, only: mpas_modify_att + use mpas_constants, only: mpas_constants_compute_derived + use mpas_derived_types, only: field0dreal, field2dreal, mpas_pool_type, mpas_time_type, & + mpas_start_time + use mpas_field_routines, only: mpas_allocate_scratch_field + use mpas_pool_routines, only: mpas_pool_get_field, mpas_pool_initialize_time_levels + use mpas_string_utils, only: mpas_string_replace + use mpas_timekeeping, only: mpas_get_clock_time, mpas_get_time + class(mpas_dynamical_core_type), intent(inout) :: self integer, intent(in) :: coupling_time_interval ! Set the time interval, in seconds, over which MPAS dynamical core ! should integrate each time it is called to run. @@ -2922,6 +2978,15 @@ end subroutine dyn_mpas_init_phase4 ! !------------------------------------------------------------------------------- subroutine dyn_mpas_run(self) + ! Module(s) from MPAS. + use atm_core, only: atm_compute_output_diagnostics, atm_do_timestep + use mpas_derived_types, only: mpas_pool_type, mpas_time_type, mpas_timeinterval_type, & + mpas_now + use mpas_pool_routines, only: mpas_pool_shift_time_levels + use mpas_timekeeping, only: mpas_advance_clock, mpas_get_clock_time, mpas_get_time, & + mpas_set_timeinterval, & + operator(+), operator(<) + class(mpas_dynamical_core_type), intent(inout) :: self character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_run' @@ -3012,6 +3077,143 @@ subroutine dyn_mpas_run(self) call self % debug_print(subname // ' completed') end subroutine dyn_mpas_run + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_final + ! + !> \brief Finalizes MPAS dynamical core as well as its framework + !> \author Michael Duda + !> \date 29 February 2020 + !> \details + !> This subroutine finalizes and cleans up MPAS dynamical core as well as its + !> framework that was set up during initialization. Finalization happens in + !> reverse chronological order. + !> Essentially, it closely follows what is done in `atm_core_finalize` and + !> `mpas_finalize`, except that here, there is no need to call MPAS diagnostics + !> manager or MPAS stream manager. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-10-10) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_final(self) + ! Module(s) from MPAS. + use atm_time_integration, only: mpas_atm_dynamics_finalize + use mpas_atm_halos, only: atm_destroy_halo_groups, exchange_halo_group + use mpas_atm_threading, only: mpas_atm_threading_finalize + use mpas_decomp, only: mpas_decomp_destroy_decomp_list + use mpas_derived_types, only: field2dreal + use mpas_field_routines, only: mpas_deallocate_scratch_field + use mpas_framework, only: mpas_framework_finalize + use mpas_log, only: mpas_log_finalize + use mpas_pool_routines, only: mpas_pool_get_field + use mpas_timekeeping, only: mpas_destroy_clock + use mpas_timer, only: mpas_timer_write_header, mpas_timer_write, mpas_timer_finalize + + class(mpas_dynamical_core_type), intent(inout) :: self + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_final' + integer :: ierr + type(field2dreal), pointer :: field_2d_real + + call self % debug_print(subname // ' entered') + + nullify(field_2d_real) + + ! First, wind down MPAS dynamical core by calling its own finalization procedures. + + ! Some additional "scratch" fields are needed for interoperability with CAM-SIMA, but they are not finalized by + ! `mpas_atm_dynamics_finalize`. Finalize them below. + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_uzonal', field_2d_real, timelevel=1) + call mpas_deallocate_scratch_field(field_2d_real) + nullify(field_2d_real) + + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_umerid', field_2d_real, timelevel=1) + call mpas_deallocate_scratch_field(field_2d_real) + nullify(field_2d_real) + + ! Opposite to `mpas_atm_dynamics_init`. + call mpas_atm_dynamics_finalize(self % domain_ptr) + + ! Opposite to `atm_build_halo_groups`. + call atm_destroy_halo_groups(self % domain_ptr, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to destroy halo exchange groups', subname, __LINE__) + end if + + nullify(exchange_halo_group) + + ! Opposite to `mpas_atm_threading_init`. + call mpas_atm_threading_finalize(self % domain_ptr % blocklist, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to clean up OpenMP threading', subname, __LINE__) + end if + + ! Opposite to `mpas_create_clock`, which was called by `atm_simulation_clock_init`, then `atm_setup_clock`. + call mpas_destroy_clock(self % domain_ptr % clock, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to clean up clock', subname, __LINE__) + end if + + ! Opposite to `mpas_decomp_create_decomp_list`, which was called by `atm_setup_decompositions`. + call mpas_decomp_destroy_decomp_list(self % domain_ptr % decompositions) + + deallocate(self % domain_ptr % streaminfo) + + ! Write timing information to log. + call mpas_timer_write_header() + call mpas_timer_write() + + ! Opposite to `mpas_timer_init`, which was called by `mpas_framework_init_phase2`. + call mpas_timer_finalize(self % domain_ptr) + + ! Opposite to `mpas_log_init`, which was called by `atm_setup_log`. + call mpas_log_finalize(ierr) + + if (ierr /= 0) then + call self % model_error('Failed to clean up log', subname, __LINE__) + end if + + ! Opposite to `mpas_framework_init_phase1` and `mpas_framework_init_phase2`. + call mpas_framework_finalize(self % domain_ptr % dminfo, self % domain_ptr) + + call self % debug_print(subname // ' completed') + + call self % debug_print('Successful finalization of MPAS dynamical core') + + ! Second, clean up this MPAS dynamical core instance. + + ! Initialized by `dyn_mpas_init_phase4`. + self % coupling_time_interval = 0 + self % number_of_time_steps = 0 + + ! Initialized by `dyn_mpas_define_scalar`. + deallocate(self % constituent_name) + deallocate(self % index_constituent_to_mpas_scalar) + deallocate(self % index_mpas_scalar_to_constituent) + deallocate(self % is_water_species) + + ! Initialized by `dyn_mpas_init_phase3`. + self % number_of_constituents = 0 + + ! Initialized by `dyn_mpas_init_phase1`. + self % log_unit = output_unit + self % mpi_comm = mpi_comm_null + self % mpi_rank = 0 + self % mpi_rank_root = .false. + + nullify(self % model_error) + + deallocate(self % corelist % domainlist) + deallocate(self % corelist) + + nullify(self % corelist) + nullify(self % domain_ptr) + + self % debug_output = .false. + end subroutine dyn_mpas_final + !------------------------------------------------------------------------------- ! function dyn_mpas_get_constituent_name ! @@ -3231,6 +3433,9 @@ end subroutine dyn_mpas_get_local_mesh_dimension subroutine dyn_mpas_get_global_mesh_dimension(self, & ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, & sphere_radius) + ! Module(s) from MPAS. + use mpas_dmpar, only: mpas_dmpar_max_int, mpas_dmpar_sum_int + class(mpas_dynamical_core_type), intent(in) :: self integer, intent(out) :: ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max real(rkind), intent(out) :: sphere_radius @@ -3278,6 +3483,10 @@ end subroutine dyn_mpas_get_global_mesh_dimension !> It is used by the `dyn_mpas_get_variable_{pointer,value}_*` subroutines to draw a variable from a pool. !> (KCW, 2024-03-21) subroutine dyn_mpas_get_pool_pointer(self, pool_pointer, pool_name) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_subpool + class(mpas_dynamical_core_type), intent(in) :: self type(mpas_pool_type), pointer, intent(out) :: pool_pointer character(*), intent(in) :: pool_name @@ -3323,6 +3532,10 @@ end subroutine dyn_mpas_get_pool_pointer ! !------------------------------------------------------------------------------- subroutine dyn_mpas_get_variable_pointer_c0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_config + class(mpas_dynamical_core_type), intent(in) :: self character(strkind), pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3351,6 +3564,10 @@ subroutine dyn_mpas_get_variable_pointer_c0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_c0 subroutine dyn_mpas_get_variable_pointer_c1(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self character(strkind), pointer, intent(out) :: variable_pointer(:) character(*), intent(in) :: pool_name @@ -3373,6 +3590,10 @@ subroutine dyn_mpas_get_variable_pointer_c1(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_c1 subroutine dyn_mpas_get_variable_pointer_i0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_config, mpas_pool_get_dimension + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3404,6 +3625,10 @@ subroutine dyn_mpas_get_variable_pointer_i0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i0 subroutine dyn_mpas_get_variable_pointer_i1(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_dimension + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer(:) character(*), intent(in) :: pool_name @@ -3432,6 +3657,10 @@ subroutine dyn_mpas_get_variable_pointer_i1(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i1 subroutine dyn_mpas_get_variable_pointer_i2(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer(:, :) character(*), intent(in) :: pool_name @@ -3454,6 +3683,10 @@ subroutine dyn_mpas_get_variable_pointer_i2(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i2 subroutine dyn_mpas_get_variable_pointer_i3(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer(:, :, :) character(*), intent(in) :: pool_name @@ -3476,6 +3709,10 @@ subroutine dyn_mpas_get_variable_pointer_i3(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i3 subroutine dyn_mpas_get_variable_pointer_l0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_config + class(mpas_dynamical_core_type), intent(in) :: self logical, pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3502,6 +3739,10 @@ subroutine dyn_mpas_get_variable_pointer_l0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_l0 subroutine dyn_mpas_get_variable_pointer_r0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_config + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3530,6 +3771,10 @@ subroutine dyn_mpas_get_variable_pointer_r0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r0 subroutine dyn_mpas_get_variable_pointer_r1(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:) character(*), intent(in) :: pool_name @@ -3552,6 +3797,10 @@ subroutine dyn_mpas_get_variable_pointer_r1(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r1 subroutine dyn_mpas_get_variable_pointer_r2(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :) character(*), intent(in) :: pool_name @@ -3574,6 +3823,10 @@ subroutine dyn_mpas_get_variable_pointer_r2(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r2 subroutine dyn_mpas_get_variable_pointer_r3(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :, :) character(*), intent(in) :: pool_name @@ -3596,6 +3849,10 @@ subroutine dyn_mpas_get_variable_pointer_r3(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r3 subroutine dyn_mpas_get_variable_pointer_r4(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :, :, :) character(*), intent(in) :: pool_name @@ -3618,6 +3875,10 @@ subroutine dyn_mpas_get_variable_pointer_r4(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r4 subroutine dyn_mpas_get_variable_pointer_r5(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :, :, :, :) character(*), intent(in) :: pool_name diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 3911af4f..f6e3a10a 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -1,46 +1,6 @@ module dyn_comp - use dyn_mpas_subdriver, only: mpas_dynamical_core_type - - ! Modules from CAM-SIMA. - use air_composition, only: thermodynamic_active_species_num, & - thermodynamic_active_species_liq_num, & - thermodynamic_active_species_ice_num, & - thermodynamic_active_species_idx, thermodynamic_active_species_idx_dycore, & - thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, & - thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore - use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: const_name, const_is_dry, const_is_water_species, num_advected, readtrace - use cam_control_mod, only: initial_run - use cam_field_read, only: cam_read_field - use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id - use cam_initfiles, only: initial_file_get_id, topo_file_get_id - use cam_instance, only: atm_id - use cam_logfile, only: debug_output, debugout_none, iulog - use cam_pio_utils, only: clean_iodesc_list - use dyn_tests_utils, only: vc_height - use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, constant_pi => pi, & - constant_rd => rair, constant_rv => rh2o, & - deg_to_rad - use inic_analytic, only: analytic_ic_active, dyn_set_inic_col - use runtime_obj, only: runtime_options - use spmd_utils, only: iam, masterproc, mpicom - use string_utils, only: stringify - use time_manager, only: get_start_date, get_stop_date, get_step_size, get_run_duration, timemgr_get_calendar_cf - use vert_coord, only: pver, pverp - - ! Modules from CCPP. - use cam_ccpp_cap, only: cam_constituents_array - use ccpp_kinds, only: kind_phys - use phys_vars_init_check, only: mark_as_initialized, std_name_len - use physics_types, only: phys_state - - ! Modules from CESM Share. - use shr_file_mod, only: shr_file_getunit - use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8 - use shr_pio_mod, only: shr_pio_getiosys - - ! Modules from external libraries. - use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind, mpas_dynamical_core_type implicit none @@ -51,10 +11,11 @@ module dyn_comp public :: dyn_readnl public :: dyn_init public :: dyn_run - ! public :: dyn_final + public :: dyn_final public :: dyn_debug_print - public :: dyn_exchange_constituent_state + public :: dyn_exchange_constituent_states + public :: dyn_inquire_mesh_dimensions public :: reverse public :: mpas_dynamical_core public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels @@ -79,14 +40,19 @@ module dyn_comp type(mpas_dynamical_core_type) :: mpas_dynamical_core ! Local and global mesh dimensions of MPAS dynamical core. - integer :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels - integer :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max - real(kind_r8) :: sphere_radius + integer, protected :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels + integer, protected :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max + real(kind_dyn_mpas), protected :: sphere_radius contains !> Print a debug message with optionally the value(s) of a variable. !> If `printer` is not supplied, the MPI root rank will print. Otherwise, the designated MPI rank will print instead. !> (KCW, 2024-02-03) subroutine dyn_debug_print(message, variable, printer) + ! Module(s) from CAM-SIMA. + use cam_logfile, only: debug_output, debugout_none, iulog + use spmd_utils, only: iam, masterproc + use string_utils, only: stringify + character(*), intent(in) :: message class(*), optional, intent(in) :: variable(:) integer, optional, intent(in) :: printer @@ -121,10 +87,24 @@ end subroutine dyn_debug_print ! ! Called by `read_namelist` in `src/control/runtime_opts.F90`. subroutine dyn_readnl(namelist_path) + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: endrun + use cam_control_mod, only: initial_run + use cam_instance, only: atm_id + use cam_logfile, only: debug_output, debugout_none, iulog + use spmd_utils, only: mpicom + use time_manager, only: get_start_date, get_stop_date, get_run_duration, timemgr_get_calendar_cf + ! Module(s) from CESM Share. + use shr_file_mod, only: shr_file_getunit + use shr_kind_mod, only: len_cs => shr_kind_cs + use shr_pio_mod, only: shr_pio_getiosys + ! Module(s) from external libraries. + use pio, only: iosystem_desc_t + character(*), intent(in) :: namelist_path character(*), parameter :: subname = 'dyn_comp::dyn_readnl' - character(kind_cs) :: cam_calendar + character(len_cs) :: cam_calendar integer :: log_unit(2) integer :: start_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. stop_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. @@ -187,9 +167,27 @@ end subroutine dyn_readnl ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) - use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS - use phys_vars_init_check, only: mark_as_initialized - use physics_types, only: dycore_energy_consistency_adjust + ! Module(s) from CAM-SIMA. + use air_composition, only: thermodynamic_active_species_num, & + thermodynamic_active_species_liq_num, & + thermodynamic_active_species_ice_num, & + thermodynamic_active_species_idx, thermodynamic_active_species_idx_dycore, & + thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, & + thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore + use cam_abortutils, only: check_allocate + use cam_constituents, only: const_name, const_is_water_species, num_advected, readtrace + use cam_control_mod, only: initial_run + use cam_initfiles, only: initial_file_get_id, topo_file_get_id + use cam_pio_utils, only: clean_iodesc_list + use cam_thermo_formula, only: energy_formula_dycore, energy_formula_dycore_mpas + use inic_analytic, only: analytic_ic_active + use physics_types, only: dycore_energy_consistency_adjust + use runtime_obj, only: runtime_options + use time_manager, only: get_step_size + ! Module(s) from CCPP. + use phys_vars_init_check, only: std_name_len + ! Module(s) from external libraries. + use pio, only: file_desc_t type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in @@ -207,14 +205,13 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) nullify(pio_init_file) nullify(pio_topo_file) - ! Set dynamical core energy formula for use in cam_thermo. - energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS - call mark_as_initialized('total_energy_formula_for_dycore') + ! Set the energy formula of dynamical core to MPAS for use in `cam_thermo`. + energy_formula_dycore = energy_formula_dycore_mpas - ! Dynamical core energy is not consistent with CAM physics and requires - ! temperature and temperature tendency adjustment at end of physics. + ! The total energy of dynamical core, which uses "MPAS formula" as set above, is not consistent with + ! that of CAM physics, which uses "FV formula". Therefore, temperature and temperature tendency adjustments + ! are needed at the end of each physics time step. dycore_energy_consistency_adjust = .true. - call mark_as_initialized('flag_for_dycore_energy_consistency_adjustment') allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) @@ -267,9 +264,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! Perform default initialization for all constituents. ! Subsequently, they can be overridden depending on the namelist option (below) and ! the actual availability (checked and handled by MPAS). - call dyn_debug_print('Calling dyn_exchange_constituent_state') + call dyn_debug_print('Calling dyn_exchange_constituent_states') - call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.false.) + call dyn_exchange_constituent_states(direction='e', exchange=.true., conversion=.false.) ! Namelist option that controls if constituents are to be read from the file. if (readtrace) then @@ -288,7 +285,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) end if call clean_iodesc_list() - call mark_variable_as_initialized() + call mark_variables_as_initialized() nullify(pio_init_file) nullify(pio_topo_file) @@ -308,6 +305,16 @@ end subroutine dyn_init !> Otherwise, if topography file is not used, check that the surface geometric height in MPAS is zero. !> (KCW, 2024-05-10) subroutine check_topography_data(pio_file) + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_field_read, only: cam_read_field + use cam_logfile, only: debug_output, debugout_none + use dynconst, only: constant_g => gravit + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open + type(file_desc_t), pointer, intent(in) :: pio_file character(*), parameter :: subname = 'dyn_comp::check_topography_data' @@ -316,7 +323,7 @@ subroutine check_topography_data(pio_file) real(kind_r8), parameter :: error_tolerance = 1.0E-3_kind_r8 ! Error tolerance for consistency check. real(kind_r8), allocatable :: surface_geometric_height(:) ! Computed from topography file. real(kind_r8), allocatable :: surface_geopotential(:) ! Read from topography file. - real(kind_r8), pointer :: zgrid(:, :) ! From MPAS. Geometric height (meters) at layer interfaces. + real(kind_dyn_mpas), pointer :: zgrid(:, :) ! From MPAS. Geometric height (meters) at layer interfaces. nullify(zgrid) @@ -348,7 +355,7 @@ subroutine check_topography_data(pio_file) surface_geometric_height(:) = surface_geopotential(:) / constant_g ! Surface geometric height in MPAS should match the values in topography file. - if (any(abs(zgrid(1, 1:ncells_solve) - surface_geometric_height) > error_tolerance)) then + if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8) - surface_geometric_height(:)) > error_tolerance)) then call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__) end if @@ -358,7 +365,7 @@ subroutine check_topography_data(pio_file) call dyn_debug_print('Topography file is not used') ! Surface geometric height in MPAS should be zero. - if (any(abs(zgrid(1, 1:ncells_solve)) > error_tolerance)) then + if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8)) > error_tolerance)) then call endrun('Surface geometric height in MPAS is not zero', subname, __LINE__) end if end if @@ -369,14 +376,17 @@ end subroutine check_topography_data !> Set analytic initial condition for MPAS. !> (KCW, 2024-05-22) subroutine set_analytic_initial_condition() + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition' integer, allocatable :: global_grid_index(:) real(kind_r8), allocatable :: buffer_2d_real(:, :), buffer_3d_real(:, :, :) real(kind_r8), allocatable :: lat_rad(:), lon_rad(:) - real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces. - ! Dimension and vertical index orders follow CAM-SIMA convention. - real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. - ! Dimension and vertical index orders follow MPAS convention. + real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow CAM-SIMA convention. + real(kind_dyn_mpas), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow MPAS convention. call init_shared_variables() @@ -391,9 +401,15 @@ subroutine set_analytic_initial_condition() !> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. !> (KCW, 2024-05-13) subroutine init_shared_variables() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id + use dynconst, only: deg_to_rad + use vert_coord, only: pverp + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variables' - integer :: ierr integer :: i + integer :: ierr integer, pointer :: indextocellid(:) real(kind_r8), pointer :: lat_deg(:), lon_deg(:) @@ -442,7 +458,7 @@ subroutine init_shared_variables() ! Vertical index order is reversed between CAM-SIMA and MPAS. do i = 1, ncells_solve - z_int(i, :) = reverse(zgrid(:, i)) + z_int(i, :) = reverse(real(zgrid(:, i), kind_r8)) end do end subroutine init_shared_variables @@ -461,10 +477,16 @@ end subroutine final_shared_variables !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces). !> (KCW, 2024-05-13) subroutine set_mpas_state_u() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_tests_utils, only: vc_height + use inic_analytic, only: dyn_set_inic_col + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u' - integer :: ierr integer :: i - real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :) + integer :: ierr + real(kind_dyn_mpas), pointer :: ucellzonal(:, :), ucellmeridional(:, :) call dyn_debug_print('Setting MPAS state "u"') @@ -482,7 +504,7 @@ subroutine set_mpas_state_u() ! Vertical index order is reversed between CAM-SIMA and MPAS. do i = 1, ncells_solve - ucellzonal(:, i) = reverse(buffer_2d_real(i, :)) + ucellzonal(:, i) = real(reverse(buffer_2d_real(i, :)), kind_dyn_mpas) end do buffer_2d_real(:, :) = 0.0_kind_r8 @@ -491,21 +513,21 @@ subroutine set_mpas_state_u() ! Vertical index order is reversed between CAM-SIMA and MPAS. do i = 1, ncells_solve - ucellmeridional(:, i) = reverse(buffer_2d_real(i, :)) + ucellmeridional(:, i) = real(reverse(buffer_2d_real(i, :)), kind_dyn_mpas) end do deallocate(buffer_2d_real) nullify(ucellzonal, ucellmeridional) - call mpas_dynamical_core % compute_edge_wind(.false.) + call mpas_dynamical_core % compute_edge_wind(wind_tendency=.false.) end subroutine set_mpas_state_u !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces). !> (KCW, 2024-05-13) subroutine set_mpas_state_w() character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_w' - real(kind_r8), pointer :: w(:, :) + real(kind_dyn_mpas), pointer :: w(:, :) call dyn_debug_print('Setting MPAS state "w"') @@ -513,7 +535,7 @@ subroutine set_mpas_state_w() call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1) - w(:, 1:ncells_solve) = 0.0_kind_r8 + w(:, 1:ncells_solve) = 0.0_kind_dyn_mpas nullify(w) @@ -524,6 +546,13 @@ end subroutine set_mpas_state_w !> Set MPAS state `scalars` (i.e., constituents). !> (KCW, 2024-05-17) subroutine set_mpas_state_scalars() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_constituents, only: num_advected + use dyn_tests_utils, only: vc_height + use inic_analytic, only: dyn_set_inic_col + use vert_coord, only: pver + ! CCPP standard name of `qv`, which denotes water vapor mixing ratio. character(*), parameter :: constituent_qv_standard_name = & 'water_vapor_mixing_ratio_wrt_dry_air' @@ -533,7 +562,7 @@ subroutine set_mpas_state_scalars() integer :: ierr integer, allocatable :: constituent_index(:) integer, pointer :: index_qv - real(kind_r8), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) call dyn_debug_print('Setting MPAS state "scalars"') @@ -560,7 +589,7 @@ subroutine set_mpas_state_scalars() do j = 1, num_advected ! Vertical index order is reversed between CAM-SIMA and MPAS. scalars(j, :, i) = & - reverse(buffer_3d_real(i, :, mpas_dynamical_core % map_constituent_index(j))) + real(reverse(buffer_3d_real(i, :, mpas_dynamical_core % map_constituent_index(j))), kind_dyn_mpas) end do end do @@ -574,7 +603,7 @@ subroutine set_mpas_state_scalars() ! Convert specific humidity to water vapor mixing ratio. scalars(index_qv, :, 1:ncells_solve) = & - scalars(index_qv, :, 1:ncells_solve) / (1.0_kind_r8 - scalars(index_qv, :, 1:ncells_solve)) + scalars(index_qv, :, 1:ncells_solve) / (1.0_kind_dyn_mpas - scalars(index_qv, :, 1:ncells_solve)) end if deallocate(buffer_3d_real) @@ -590,6 +619,13 @@ end subroutine set_mpas_state_scalars !> Set MPAS state `rho` (i.e., dry air density) and `theta` (i.e., potential temperature). !> (KCW, 2024-05-19) subroutine set_mpas_state_rho_theta() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_tests_utils, only: vc_height + use dynconst, only: constant_p0 => pref, constant_rd => rair, constant_rv => rh2o + use inic_analytic, only: dyn_set_inic_col + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_theta' integer :: i, k integer :: ierr @@ -603,9 +639,9 @@ subroutine set_mpas_state_rho_theta() real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column. ! Be advised that it is not virtual temperature. ! See doi:10.5065/1DFH-6P97 and doi:10.1175/MWR-D-11-00215.1 for details. - real(kind_r8), pointer :: rho(:, :) - real(kind_r8), pointer :: theta(:, :) - real(kind_r8), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: rho(:, :) + real(kind_dyn_mpas), pointer :: theta(:, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) call dyn_debug_print('Setting MPAS state "rho" and "theta"') @@ -654,16 +690,16 @@ subroutine set_mpas_state_rho_theta() ! Set `rho` and `theta` column by column. This way, peak memory usage can be reduced. do i = 1, ncells_solve - qv_mid_col(:) = scalars(index_qv, :, i) + qv_mid_col(:) = real(scalars(index_qv, :, i), kind_r8) tm_mid_col(:) = t_mid(:, i) * (1.0_kind_r8 + constant_rv / constant_rd * qv_mid_col(:)) ! Piecewise integrate hypsometric equation to derive `p_mid_col(1)`. ! The formulation used here is exact. p_mid_col(1) = p_by_hypsometric_equation( & p_sfc(i), & - zgrid(1, i), & + real(zgrid(1, i), kind_r8), & tm_mid_col(1) / (1.0_kind_r8 + qv_mid_col(1)), & - 0.5_kind_r8 * (zgrid(2, i) + zgrid(1, i))) + 0.5_kind_r8 * real(zgrid(2, i) + zgrid(1, i), kind_r8)) ! Piecewise integrate hypsometric equation to derive subsequent `p_mid_col(k)`. ! The formulation used here is exact. @@ -671,16 +707,16 @@ subroutine set_mpas_state_rho_theta() p_mid_col(k) = p_by_hypsometric_equation( & p_by_hypsometric_equation( & p_mid_col(k - 1), & - 0.5_kind_r8 * (zgrid(k, i) + zgrid(k - 1, i)), & + 0.5_kind_r8 * real(zgrid(k, i) + zgrid(k - 1, i), kind_r8), & tm_mid_col(k - 1) / (1.0_kind_r8 + qv_mid_col(k - 1)), & - zgrid(k, i)), & - zgrid(k, i), & + real(zgrid(k, i), kind_r8)), & + real(zgrid(k, i), kind_r8), & tm_mid_col(k) / (1.0_kind_r8 + qv_mid_col(k)), & - 0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i))) + 0.5_kind_r8 * real(zgrid(k + 1, i) + zgrid(k, i), kind_r8)) end do - rho(:, i) = p_mid_col(:) / (constant_rd * tm_mid_col(:)) - theta(:, i) = theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0) + rho(:, i) = real(p_mid_col(:) / (constant_rd * tm_mid_col(:)), kind_dyn_mpas) + theta(:, i) = real(theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0), kind_dyn_mpas) end do deallocate(p_mid_col) @@ -702,15 +738,20 @@ end subroutine set_mpas_state_rho_theta !> Set MPAS state `rho_base` (i.e., base state dry air density) and `theta_base` (i.e., base state potential temperature). !> (KCW, 2024-05-21) subroutine set_mpas_state_rho_base_theta_base() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dynconst, only: constant_p0 => pref, constant_rd => rair + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_base_theta_base' integer :: i, k integer :: ierr real(kind_r8), parameter :: t_base = 250.0_kind_r8 ! Base state temperature (K) of dry isothermal atmosphere. ! The value used here is identical to MPAS. real(kind_r8), allocatable :: p_base(:) ! Base state pressure (Pa) at layer midpoints of each column. - real(kind_r8), pointer :: rho_base(:, :) - real(kind_r8), pointer :: theta_base(:, :) - real(kind_r8), pointer :: zz(:, :) + real(kind_dyn_mpas), pointer :: rho_base(:, :) + real(kind_dyn_mpas), pointer :: theta_base(:, :) + real(kind_dyn_mpas), pointer :: zz(:, :) call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"') @@ -734,11 +775,11 @@ subroutine set_mpas_state_rho_base_theta_base() constant_p0, & 0.0_kind_r8, & t_base, & - 0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i))) + 0.5_kind_r8 * real(zgrid(k + 1, i) + zgrid(k, i), kind_r8)) end do - rho_base(:, i) = p_base(:) / (constant_rd * t_base * zz(:, i)) - theta_base(:, i) = theta_by_poisson_equation(p_base, t_base, constant_p0) + rho_base(:, i) = real(p_base(:) / (constant_rd * t_base * real(zz(:, i), kind_r8)), kind_dyn_mpas) + theta_base(:, i) = real(theta_by_poisson_equation(p_base, t_base, constant_p0), kind_dyn_mpas) end do deallocate(p_base) @@ -760,6 +801,9 @@ end subroutine set_mpas_state_rho_base_theta_base !> `t_v` is the mean virtual temperature between `z_1` and `z_2`. !> (KCW, 2024-07-02) pure elemental function p_by_hypsometric_equation(p_1, z_1, t_v, z_2) result(p_2) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_g => gravit, constant_rd => rair + real(kind_r8), intent(in) :: p_1, z_1, t_v, z_2 real(kind_r8) :: p_2 @@ -774,6 +818,9 @@ end function p_by_hypsometric_equation !> Poisson equation. !> (KCW, 2024-07-02) pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_cpd => cpair, constant_rd => rair + real(kind_r8), intent(in) :: p_1, t_1, p_0 real(kind_r8) :: t_0 @@ -786,18 +833,29 @@ end subroutine set_analytic_initial_condition !> If `exchange` is `.true.` and `direction` is "i" or "import", set physics state `constituents` from MPAS state `scalars`. !> Think of it as "exporting/importing constituent states in CAM-SIMA to/from MPAS". !> Otherwise, if `exchange` is `.false.`, no exchange is performed at all. - !> If `conversion` is `.true.`, appropriate conversion is performed for constituent mixing ratio that has different + !> If `conversion` is `.true.`, appropriate conversion is performed for constituent mixing ratios that have different !> definitions between CAM-SIMA and MPAS (i.e., dry/moist). !> Otherwise, if `conversion` is `.false.`, no conversion is performed at all. !> This subroutine is intentionally designed to have these elaborate controls due to complications in CAM-SIMA. !> Some procedures in CAM-SIMA expect constituent states to be dry, while the others expect them to be moist. !> (KCW, 2024-09-26) - subroutine dyn_exchange_constituent_state(direction, exchange, conversion) + subroutine dyn_exchange_constituent_states(direction, exchange, conversion) + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_constituents, only: const_is_dry, const_is_water_species, num_advected + use physics_types, only: phys_state + use vert_coord, only: pver + ! Module(s) from CCPP. + use cam_ccpp_cap, only: cam_constituents_array + use ccpp_kinds, only: kind_phys + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + character(*), intent(in) :: direction logical, intent(in) :: exchange logical, intent(in) :: conversion - character(*), parameter :: subname = 'dyn_comp::dyn_exchange_constituent_state' + character(*), parameter :: subname = 'dyn_comp::dyn_exchange_constituent_states' integer :: i, j integer :: ierr integer, allocatable :: is_water_species_index(:) @@ -805,7 +863,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) logical, allocatable :: is_water_species(:) real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water species mixing ratios. - real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. + real(kind_dyn_mpas), pointer :: scalars(:, :, :) ! This points to MPAS memory. select case (trim(adjustl(direction))) case ('e', 'export') @@ -880,13 +938,13 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) if (exchange) then ! Vertical index order is reversed between CAM-SIMA and MPAS. scalars(j, :, i) = & - reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))) + real(reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))), kind_dyn_mpas) end if if (conversion .and. is_conversion_needed(mpas_dynamical_core % map_constituent_index(j))) then ! Equation 8 in doi:10.1029/2017MS001257. scalars(j, :, i) = & - scalars(j, :, i) * sigma_all_q(:) + real(real(scalars(j, :, i), kind_r8) * sigma_all_q(:), kind_dyn_mpas) end if end do end do @@ -897,7 +955,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) do i = 1, ncells_solve if (conversion .and. any(is_conversion_needed)) then ! The summation term of equation 8 in doi:10.1029/2017MS001257. - sigma_all_q(:) = reverse(1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1)) + sigma_all_q(:) = reverse(1.0_kind_r8 + sum(real(scalars(is_water_species_index, :, i), kind_r8), 1)) end if ! `j` is indexing into `constituents`, so it is regarded as constituent index. @@ -905,7 +963,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) if (exchange) then ! Vertical index order is reversed between CAM-SIMA and MPAS. constituents(i, :, j) = & - reverse(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i)) + reverse(real(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i), kind_r8)) end if if (conversion .and. is_conversion_needed(j)) then @@ -929,17 +987,54 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. call mpas_dynamical_core % exchange_halo('scalars') end if - end subroutine dyn_exchange_constituent_state + end subroutine dyn_exchange_constituent_states + + !> Inquire local and global mesh dimensions. Save them as protected module variables. + !> (KCW, 2024-11-21) + subroutine dyn_inquire_mesh_dimensions() + ! Module(s) from CAM-SIMA. + use string_utils, only: stringify + + character(*), parameter :: subname = 'dyn_comp::dyn_inquire_mesh_dimensions' - !> Mark everything in the `physics_{state,tend}` derived types along with constituents as initialized - !> to prevent physics from attempting to read them from a file. These variables are to be exchanged later - !> during dynamics-physics coupling. + call dyn_debug_print('Inquiring local and global mesh dimensions') + + call mpas_dynamical_core % get_local_mesh_dimension( & + ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels) + + call mpas_dynamical_core % get_global_mesh_dimension( & + ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, & + sphere_radius) + + call dyn_debug_print('ncells_global = ' // stringify([ncells_global])) + call dyn_debug_print('nedges_global = ' // stringify([nedges_global])) + call dyn_debug_print('nvertices_global = ' // stringify([nvertices_global])) + call dyn_debug_print('nvertlevels = ' // stringify([nvertlevels])) + call dyn_debug_print('ncells_max = ' // stringify([ncells_max])) + call dyn_debug_print('nedges_max = ' // stringify([nedges_max])) + call dyn_debug_print('sphere_radius = ' // stringify([sphere_radius])) + end subroutine dyn_inquire_mesh_dimensions + + !> Mark everything in the `physics_types` module along with constituents as initialized + !> to prevent physics from attempting to read them from a file. !> (KCW, 2024-05-23) - subroutine mark_variable_as_initialized() - character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized' + subroutine mark_variables_as_initialized() + ! Module(s) from CAM-SIMA. + use cam_constituents, only: const_name, num_advected + ! Module(s) from CCPP. + use phys_vars_init_check, only: mark_as_initialized + + character(*), parameter :: subname = 'dyn_comp::mark_variables_as_initialized' integer :: i - ! CCPP standard names of physical quantities in the `physics_{state,tend}` derived types. + ! The variables below are managed by dynamics interface. + ! We are responsible for initializing and updating them. + + ! These variables are to be set during dynamics initialization. + call mark_as_initialized('flag_for_dycore_energy_consistency_adjustment') + call mark_as_initialized('total_energy_formula_for_dycore') + + ! These variables are to be set during dynamics-physics coupling. call mark_as_initialized('air_pressure') call mark_as_initialized('air_pressure_at_interface') call mark_as_initialized('air_pressure_of_dry_air') @@ -960,6 +1055,7 @@ subroutine mark_variable_as_initialized() call mark_as_initialized('reciprocal_of_air_pressure_thickness') call mark_as_initialized('reciprocal_of_air_pressure_thickness_of_dry_air') call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure') + call mark_as_initialized('specific_heat_of_air_used_in_dycore') call mark_as_initialized('surface_air_pressure') call mark_as_initialized('surface_geopotential') call mark_as_initialized('surface_pressure_of_dry_air') @@ -972,10 +1068,10 @@ subroutine mark_variable_as_initialized() call mark_as_initialized(trim(adjustl(const_name(i)))) end do - call mark_as_initialized('specific_heat_of_air_used_in_dycore') + ! The variables below are not managed by dynamics interface. They are used by external CCPP physics schemes. + ! While we are not responsible for initializing or updating them, we still need to help mark them as initialized. - ! These energy variables are calculated by check_energy_timestep_init - ! but need to be marked here + ! These variables are to be set externally by the `check_energy_chng` CCPP physics scheme. call mark_as_initialized('vertically_integrated_total_energy_at_end_of_physics_timestep') call mark_as_initialized('vertically_integrated_total_energy_using_dycore_energy_formula') call mark_as_initialized('vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep') @@ -983,7 +1079,7 @@ subroutine mark_variable_as_initialized() call mark_as_initialized('vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep') call mark_as_initialized('vertically_integrated_total_water') call mark_as_initialized('vertically_integrated_total_water_at_start_of_physics_timestep') - end subroutine mark_variable_as_initialized + end subroutine mark_variables_as_initialized !> Run MPAS dynamical core to integrate the dynamical states with time. !> (KCW, 2024-07-11) @@ -994,13 +1090,80 @@ subroutine dyn_run() call mpas_dynamical_core % run() end subroutine dyn_run - ! Not used for now. Intended to be called by `stepon_final` in `src/dynamics/mpas/stepon.F90`. - ! subroutine dyn_final() - ! end subroutine dyn_final + !> Finalize MPAS dynamical core as well as its framework. + !> (KCW, 2024-10-04) + subroutine dyn_final() + character(*), parameter :: subname = 'dyn_comp::dyn_final' + + ! Quick hack for dumping variables from MPAS dynamical core. + ! Remove it once history and restart are wired up in CAM-SIMA. + call dyn_variable_dump() + + ! After this point, do not access anything under MPAS dynamical core or runtime errors will ensue. + call mpas_dynamical_core % final() + end subroutine dyn_final + + subroutine dyn_variable_dump() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_instance, only: atm_id + use physics_types, only: phys_state + ! Module(s) from CESM Share. + use shr_pio_mod, only: shr_pio_getioformat, shr_pio_getiosys, shr_pio_getiotype + ! Module(s) from external libraries. + use pio, only: file_desc_t, iosystem_desc_t, pio_createfile, pio_closefile, pio_clobber, pio_noerr + + character(*), parameter :: subname = 'dyn_comp::dyn_variable_dump' + integer :: ierr + integer :: pio_ioformat, pio_iotype + real(kind_dyn_mpas), pointer :: surface_pressure(:) + type(file_desc_t), pointer :: pio_file + type(iosystem_desc_t), pointer :: pio_iosystem + + nullify(pio_file) + nullify(pio_iosystem) + nullify(surface_pressure) + + call mpas_dynamical_core % get_variable_pointer(surface_pressure, 'diag', 'surface_pressure') + + surface_pressure(1:ncells_solve) = real(phys_state % ps(:), kind_dyn_mpas) + + nullify(surface_pressure) + + call mpas_dynamical_core % exchange_halo('surface_pressure') + + allocate(pio_file, stat=ierr) + call check_allocate(ierr, subname, 'pio_file', 'dyn_comp', __LINE__) + + pio_iosystem => shr_pio_getiosys(atm_id) + + pio_ioformat = shr_pio_getioformat(atm_id) + pio_ioformat = ior(pio_ioformat, pio_clobber) + + pio_iotype = shr_pio_getiotype(atm_id) + + ierr = pio_createfile(pio_iosystem, pio_file, pio_iotype, 'dyn_variable_dump.nc', pio_ioformat) + + if (ierr /= pio_noerr) then + call endrun('Failed to create file for variable dumping', subname, __LINE__) + end if + + call mpas_dynamical_core % read_write_stream(pio_file, 'w', 'invariant+input+restart+output') + + call pio_closefile(pio_file) + + deallocate(pio_file) + + nullify(pio_file) + nullify(pio_iosystem) + end subroutine dyn_variable_dump !> Helper function for reversing the order of elements in `array`. !> (KCW, 2024-07-17) pure function reverse(array) + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + real(kind_r8), intent(in) :: array(:) real(kind_r8) :: reverse(size(array)) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index a2066e16..aa176a73 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -1,32 +1,4 @@ module dyn_coupling - ! Modules from CAM-SIMA. - use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: const_is_water_species, const_qmin, num_advected - use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update - use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_MPAS - use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, reverse, mpas_dynamical_core, & - ncells_solve - use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, & - constant_rd => rair, constant_rv => rh2o - use runtime_obj, only: cam_runtime_opts - use vert_coord, only: pver, pverp - - ! Modules from CCPP. - use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties - use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t - use ccpp_kinds, only: kind_phys - use geopotential_temp, only: geopotential_temp_run - use physics_types, only: cappav, cpairv, rairv, zvirv, & - dtime_phys, lagrangian_vertical, & - phys_state, phys_tend - use physics_types, only: cp_or_cv_dycore - use qneg, only: qneg_run - use static_energy, only: update_dry_static_energy_run - use string_utils, only: stringify - - ! Modules from CESM Share. - use shr_kind_mod, only: kind_cx => shr_kind_cx, kind_r8 => shr_kind_r8 - implicit none private @@ -38,6 +10,13 @@ module dyn_coupling !> The other coupling direction is implemented by its counterpart, `physics_to_dynamics_coupling`. !> (KCW, 2024-07-31) subroutine dynamics_to_physics_coupling() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_states, ncells_solve + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling' integer :: column_index integer, allocatable :: is_water_species_index(:) @@ -64,17 +43,17 @@ subroutine dynamics_to_physics_coupling() real(kind_r8), allocatable :: u_mid_col(:), & ! Eastward wind velocity (m s-1). v_mid_col(:), & ! Northward wind velocity (m s-1). omega_mid_col(:) ! Vertical wind velocity (Pa s-1). - real(kind_r8), pointer :: exner(:, :) - real(kind_r8), pointer :: rho_zz(:, :) - real(kind_r8), pointer :: scalars(:, :, :) - real(kind_r8), pointer :: theta_m(:, :) - real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :) - real(kind_r8), pointer :: zgrid(:, :) - real(kind_r8), pointer :: zz(:, :) + real(kind_dyn_mpas), pointer :: exner(:, :) + real(kind_dyn_mpas), pointer :: rho_zz(:, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: theta_m(:, :) + real(kind_dyn_mpas), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :) + real(kind_dyn_mpas), pointer :: zgrid(:, :) + real(kind_dyn_mpas), pointer :: zz(:, :) call init_shared_variables() - call dyn_exchange_constituent_state(direction='i', exchange=.true., conversion=.false.) + call dyn_exchange_constituent_states(direction='i', exchange=.true., conversion=.false.) call dyn_debug_print('Setting physics state variables column by column') @@ -93,6 +72,12 @@ subroutine dynamics_to_physics_coupling() !> `set_physics_state_column` internal subroutines. !> (KCW, 2024-07-20) subroutine init_shared_variables() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_constituents, only: const_is_water_species, num_advected + use dyn_comp, only: mpas_dynamical_core + use vert_coord, only: pver, pverp + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variables' integer :: i integer :: ierr @@ -198,6 +183,10 @@ end subroutine final_shared_variables !> should be called in pairs. !> (KCW, 2024-07-30) subroutine update_shared_variables(i) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_g => gravit, constant_rd => rair, constant_rv => rh2o + use vert_coord, only: pver, pverp + integer, intent(in) :: i character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variables' @@ -209,10 +198,10 @@ subroutine update_shared_variables(i) ! Compute thermodynamic variables. ! By definition. - z_int_col(:) = zgrid(:, i) + z_int_col(:) = real(zgrid(:, i), kind_r8) dz_col(:) = z_int_col(2:pverp) - z_int_col(1:pver) - qv_mid_col(:) = scalars(index_qv, :, i) - rhod_mid_col(:) = rho_zz(:, i) * zz(:, i) + qv_mid_col(:) = real(scalars(index_qv, :, i), kind_r8) + rhod_mid_col(:) = real(rho_zz(:, i) * zz(:, i), kind_r8) ! Equation 5 in doi:10.1029/2017MS001257. rho_mid_col(:) = rhod_mid_col(:) * sigma_all_q_mid_col(:) @@ -222,7 +211,7 @@ subroutine update_shared_variables(i) dp_col(:) = -rho_mid_col(:) * constant_g * dz_col(:) ! By definition of Exner function. Also see below. - tm_mid_col(:) = theta_m(:, i) * exner(:, i) + tm_mid_col(:) = real(theta_m(:, i) * exner(:, i), kind_r8) ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. @@ -250,15 +239,20 @@ subroutine update_shared_variables(i) ! Compute momentum variables. ! By definition. - u_mid_col(:) = ucellzonal(:, i) - v_mid_col(:) = ucellmeridional(:, i) - omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * (w(1:pver, i) + w(2:pverp, i)) + u_mid_col(:) = real(ucellzonal(:, i), kind_r8) + v_mid_col(:) = real(ucellmeridional(:, i), kind_r8) + omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * real(w(1:pver, i) + w(2:pverp, i), kind_r8) end subroutine update_shared_variables !> Set variables for the specific column, indicated by `i`, in the `physics_state` derived type. !> This subroutine and `update_shared_variables` should be called in pairs. !> (KCW, 2024-07-30) subroutine set_physics_state_column(i) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: reverse + use dynconst, only: constant_g => gravit + use physics_types, only: phys_state + integer, intent(in) :: i character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_column' @@ -294,8 +288,29 @@ end subroutine set_physics_state_column !> Set variables in the `physics_state` derived type by calling external procedures. !> (KCW, 2024-07-30) subroutine set_physics_state_external() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_constituents, only: const_qmin, num_advected + use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update + use cam_thermo_formula, only: energy_formula_dycore_mpas + use dyn_comp, only: mpas_dynamical_core + use dynconst, only: constant_g => gravit + use physics_types, only: cappav, cp_or_cv_dycore, cpairv, lagrangian_vertical, phys_state, rairv, zvirv + use runtime_obj, only: cam_runtime_opts + use string_utils, only: stringify + use vert_coord, only: pver, pverp + ! Module(s) from CCPP. + use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_kinds, only: kind_phys + use geopotential_temp, only: geopotential_temp_run + use qneg, only: qneg_run + use static_energy, only: update_dry_static_energy_run + ! Module(s) from CESM Share. + use shr_kind_mod, only: len_cx => shr_kind_cx + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_external' - character(kind_cx) :: cerr + character(len_cx) :: cerr integer :: i integer :: ierr real(kind_phys), allocatable :: minimum_constituents(:) @@ -333,15 +348,10 @@ subroutine set_physics_state_external() call cam_thermo_dry_air_update( & constituents, phys_state % t, ncells_solve, pver, cam_runtime_opts % update_thermodynamic_variables()) - ! update cp_or_cv_dycore in SIMA state. - ! (note: at this point q is dry) + ! Update `cp_or_cv_dycore` by calling `cam_thermo_water_update`. + ! Note that this subroutine expects constituents to be dry. call cam_thermo_water_update( & - mmr = constituents, & ! dry MMR - ncol = ncells_solve, & - pver = pver, & - energy_formula = ENERGY_FORMULA_DYCORE_MPAS, & - cp_or_cv_dycore = cp_or_cv_dycore & - ) + constituents, ncells_solve, pver, energy_formula_dycore_mpas, cp_or_cv_dycore) ! This variable name is really misleading. It actually represents the reciprocal of Exner function ! with respect to surface pressure. This definition is sometimes used for boundary layer work. See @@ -352,7 +362,7 @@ subroutine set_physics_state_external() end do ! Note that constituents become moist after this. - call dyn_exchange_constituent_state(direction='i', exchange=.false., conversion=.true.) + call dyn_exchange_constituent_states(direction='i', exchange=.false., conversion=.true.) ! Impose minimum limits on constituents. call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr) @@ -364,7 +374,7 @@ subroutine set_physics_state_external() end if ! Set `zi` (i.e., geopotential height at layer interfaces) and `zm` (i.e., geopotential height at layer midpoints). - ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_update`. + ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_dry_air_update`. call geopotential_temp_run( & pver, lagrangian_vertical, pver, 1, pverp, 1, num_advected, & phys_state % lnpint, phys_state % pint, phys_state % pmid, phys_state % pdel, phys_state % rpdel, phys_state % t, & @@ -378,7 +388,7 @@ subroutine set_physics_state_external() end if ! Set `dse` (i.e., dry static energy). - ! Note that `cpairv` is updated externally by `cam_thermo_update`. + ! Note that `cpairv` is updated externally by `cam_thermo_dry_air_update`. call update_dry_static_energy_run( & pver, constant_g, phys_state % t, phys_state % zm, phys_state % phis, phys_state % dse, cpairv, ierr, cerr) @@ -399,17 +409,24 @@ end subroutine dynamics_to_physics_coupling !> The other coupling direction is implemented by its counterpart, `dynamics_to_physics_coupling`. !> (KCW, 2024-09-20) subroutine physics_to_dynamics_coupling() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_exchange_constituent_states + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling' integer, pointer :: index_qv real(kind_r8), allocatable :: qv_prev(:, :) ! Water vapor mixing ratio (kg kg-1) ! before being updated by physics. - real(kind_r8), pointer :: rho_zz(:, :) - real(kind_r8), pointer :: scalars(:, :, :) - real(kind_r8), pointer :: zz(:, :) + real(kind_dyn_mpas), pointer :: rho_zz(:, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: zz(:, :) call init_shared_variables() - call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.true.) + call dyn_exchange_constituent_states(direction='e', exchange=.true., conversion=.true.) call set_mpas_physics_tendency_ru() call set_mpas_physics_tendency_rho() @@ -420,6 +437,11 @@ subroutine physics_to_dynamics_coupling() !> Initialize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. !> (KCW, 2024-09-13) subroutine init_shared_variables() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, ncells_solve + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variables' integer :: ierr @@ -441,8 +463,8 @@ subroutine init_shared_variables() 'dyn_coupling', __LINE__) ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta` - ! needs it. This must be done before calling `dyn_exchange_constituent_state`. - qv_prev(:, :) = scalars(index_qv, :, 1:ncells_solve) + ! needs it. This must be done before calling `dyn_exchange_constituent_states`. + qv_prev(:, :) = real(scalars(index_qv, :, 1:ncells_solve), kind_r8) end subroutine init_shared_variables !> Finalize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. @@ -462,9 +484,13 @@ end subroutine final_shared_variables !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. !> (KCW, 2024-09-11) subroutine set_mpas_physics_tendency_ru() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_debug_print, reverse, mpas_dynamical_core, ncells_solve + use physics_types, only: phys_tend + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_ru' integer :: i - real(kind_r8), pointer :: u_tendency(:, :), v_tendency(:, :) + real(kind_dyn_mpas), pointer :: u_tendency(:, :), v_tendency(:, :) call dyn_debug_print('Setting MPAS physics tendency "tend_ru_physics"') @@ -476,21 +502,24 @@ subroutine set_mpas_physics_tendency_ru() ! Vertical index order is reversed between CAM-SIMA and MPAS. ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type. do i = 1, ncells_solve - u_tendency(:, i) = reverse(phys_tend % dudt_total(i, :)) * rho_zz(:, i) - v_tendency(:, i) = reverse(phys_tend % dvdt_total(i, :)) * rho_zz(:, i) + u_tendency(:, i) = real(reverse(phys_tend % dudt_total(i, :)) * real(rho_zz(:, i), kind_r8), kind_dyn_mpas) + v_tendency(:, i) = real(reverse(phys_tend % dvdt_total(i, :)) * real(rho_zz(:, i), kind_r8), kind_dyn_mpas) end do nullify(u_tendency, v_tendency) - call mpas_dynamical_core % compute_edge_wind(.true.) + call mpas_dynamical_core % compute_edge_wind(wind_tendency=.true.) end subroutine set_mpas_physics_tendency_ru !> Set MPAS physics tendency `tend_rho_physics` (i.e., "coupled" tendency of dry air density due to physics). !> In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. !> (KCW, 2024-09-11) subroutine set_mpas_physics_tendency_rho() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, ncells_solve + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rho' - real(kind_r8), pointer :: rho_tendency(:, :) + real(kind_dyn_mpas), pointer :: rho_tendency(:, :) call dyn_debug_print('Setting MPAS physics tendency "tend_rho_physics"') @@ -499,7 +528,7 @@ subroutine set_mpas_physics_tendency_rho() call mpas_dynamical_core % get_variable_pointer(rho_tendency, 'tend_physics', 'tend_rho_physics') ! The material derivative of `rho` (i.e., dry air density) is zero for incompressible fluid. - rho_tendency(:, 1:ncells_solve) = 0.0_kind_r8 + rho_tendency(:, 1:ncells_solve) = 0.0_kind_dyn_mpas nullify(rho_tendency) @@ -511,6 +540,13 @@ end subroutine set_mpas_physics_tendency_rho !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. !> (KCW, 2024-09-19) subroutine set_mpas_physics_tendency_rtheta() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_comp, only: dyn_debug_print, reverse, mpas_dynamical_core, ncells_solve + use dynconst, only: constant_rd => rair, constant_rv => rh2o + use physics_types, only: dtime_phys, phys_tend + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rtheta' integer :: i integer :: ierr @@ -523,8 +559,8 @@ subroutine set_mpas_physics_tendency_rtheta() real(kind_r8), allocatable :: t_col_prev(:), t_col_curr(:) ! Temperature (K). real(kind_r8), allocatable :: theta_col_prev(:), theta_col_curr(:) ! Potential temperature (K). real(kind_r8), allocatable :: thetam_col_prev(:), thetam_col_curr(:) ! Modified "moist" potential temperature (K). - real(kind_r8), pointer :: theta_m(:, :) - real(kind_r8), pointer :: theta_m_tendency(:, :) + real(kind_dyn_mpas), pointer :: theta_m(:, :) + real(kind_dyn_mpas), pointer :: theta_m_tendency(:, :) call dyn_debug_print('Setting MPAS physics tendency "tend_rtheta_physics"') @@ -561,11 +597,11 @@ subroutine set_mpas_physics_tendency_rtheta() ! Set `theta_m_tendency` column by column. This way, peak memory usage can be reduced. do i = 1, ncells_solve - qv_col_curr(:) = scalars(index_qv, :, i) + qv_col_curr(:) = real(scalars(index_qv, :, i), kind_r8) qv_col_prev(:) = qv_prev(:, i) - rhod_col(:) = rho_zz(:, i) * zz(:, i) + rhod_col(:) = real(rho_zz(:, i) * zz(:, i), kind_r8) - thetam_col_prev(:) = theta_m(:, i) + thetam_col_prev(:) = real(theta_m(:, i), kind_r8) theta_col_prev(:) = thetam_col_prev(:) / (1.0_kind_r8 + constant_rv / constant_rd * qv_col_prev(:)) t_col_prev(:) = t_of_theta_rhod_qv(theta_col_prev, rhod_col, qv_col_prev) @@ -575,7 +611,8 @@ subroutine set_mpas_physics_tendency_rtheta() theta_col_curr(:) = theta_of_t_rhod_qv(t_col_curr, rhod_col, qv_col_curr) thetam_col_curr(:) = theta_col_curr(:) * (1.0_kind_r8 + constant_rv / constant_rd * qv_col_curr(:)) - theta_m_tendency(:, i) = (thetam_col_curr(:) - thetam_col_prev(:)) * rho_zz(:, i) / dtime_phys + theta_m_tendency(:, i) = & + real((thetam_col_curr(:) - thetam_col_prev(:)) * real(rho_zz(:, i), kind_r8) / dtime_phys, kind_dyn_mpas) end do deallocate(qv_col_prev, qv_col_curr) @@ -597,6 +634,10 @@ end subroutine set_mpas_physics_tendency_rtheta !> `t == t_of_theta_rhod_qv(theta_of_t_rhod_qv(t, rhod, qv), rhod, qv)`. !> (KCW, 2024-09-13) pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_cpd => cpair, constant_p0 => pref, & + constant_rd => rair, constant_rv => rh2o + real(kind_r8), intent(in) :: theta, rhod, qv real(kind_r8) :: t @@ -631,6 +672,10 @@ end function t_of_theta_rhod_qv !> `theta == theta_of_t_rhod_qv(t_of_theta_rhod_qv(theta, rhod, qv), rhod, qv)`. !> (KCW, 2024-09-13) pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_cpd => cpair, constant_p0 => pref, & + constant_rd => rair, constant_rv => rh2o + real(kind_r8), intent(in) :: t, rhod, qv real(kind_r8) :: theta diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 07277b1c..874b751d 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -1,31 +1,6 @@ module dyn_grid - ! Modules from CAM-SIMA. - use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: num_advected - use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register, & - horiz_coord_t, horiz_coord_create, & - max_hcoordname_len - use cam_history_support, only: add_vert_coord - use cam_initfiles, only: initial_file_get_id - use cam_map_utils, only: kind_imap => imap - use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & - ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, & - ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, & - sphere_radius - use dynconst, only: constant_p0 => pref, constant_pi => pi, rad_to_deg, dynconst_init - use physics_column_type, only: kind_pcol, physics_column_t - use physics_grid, only: phys_decomp, phys_grid_init - use ref_pres, only: ref_pres_init - use spmd_utils, only: iam - use std_atm_profile, only: std_atm_pres - use string_utils, only: stringify - use vert_coord, only: pver, pverp, vert_coord_init - - ! Modules from CESM Share. - use shr_kind_mod, only: kind_r8 => shr_kind_r8 - - ! Modules from external libraries. - use pio, only: file_desc_t + ! Module(s) from CAM-SIMA. + use cam_grid_support, only: max_hcoordname_len implicit none @@ -52,6 +27,17 @@ module dyn_grid ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine model_grid_init() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: endrun + use cam_constituents, only: num_advected + use cam_initfiles, only: initial_file_get_id + use dyn_comp, only: dyn_debug_print, dyn_inquire_mesh_dimensions, mpas_dynamical_core, nvertlevels + use dynconst, only: dynconst_init + use string_utils, only: stringify + use vert_coord, only: pver, vert_coord_init + ! Module(s) from external libraries. + use pio, only: file_desc_t + character(*), parameter :: subname = 'dyn_grid::model_grid_init' type(file_desc_t), pointer :: pio_file @@ -79,26 +65,13 @@ subroutine model_grid_init() ! Read time-invariant (e.g., grid/mesh) variables. call mpas_dynamical_core % read_write_stream(pio_file, 'r', 'invariant') + nullify(pio_file) + ! Compute local east, north and edge-normal unit vectors whenever time-invariant (e.g., grid/mesh) variables are read. call mpas_dynamical_core % compute_unit_vector() - ! Inquire local and global mesh dimensions and save them as module variables. - call dyn_debug_print('Inquiring local and global mesh dimensions') - - call mpas_dynamical_core % get_local_mesh_dimension( & - ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels) - - call mpas_dynamical_core % get_global_mesh_dimension( & - ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, & - sphere_radius) - - call dyn_debug_print('ncells_global = ' // stringify([ncells_global])) - call dyn_debug_print('nedges_global = ' // stringify([nedges_global])) - call dyn_debug_print('nvertices_global = ' // stringify([nvertices_global])) - call dyn_debug_print('nvertlevels = ' // stringify([nvertlevels])) - call dyn_debug_print('ncells_max = ' // stringify([ncells_max])) - call dyn_debug_print('nedges_max = ' // stringify([nedges_max])) - call dyn_debug_print('sphere_radius = ' // stringify([sphere_radius])) + ! Inquire local and global mesh dimensions. + call dyn_inquire_mesh_dimensions() ! Check for consistency in numbers of vertical layers. if (nvertlevels /= pver) then @@ -127,6 +100,20 @@ end subroutine model_grid_init !> Initialize reference pressure for use by physics. !> (KCW, 2024-03-25) subroutine init_reference_pressure() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_history_support, only: add_vert_coord + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core + use dynconst, only: constant_p0 => pref + use ref_pres, only: ref_pres_init + use std_atm_profile, only: std_atm_pres + use string_utils, only: stringify + use vert_coord, only: pver, pverp + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind + character(*), parameter :: subname = 'dyn_grid::init_reference_pressure' ! Number of pure pressure levels at model top. integer, parameter :: num_pure_p_lev = 0 @@ -141,7 +128,7 @@ subroutine init_reference_pressure() ! `dzw` denotes the delta/difference between `zw`. ! `rdzw` denotes the reciprocal of `dzw`. real(kind_r8), allocatable :: dzw(:) - real(kind_r8), pointer :: rdzw(:) + real(kind_dyn_mpas), pointer :: rdzw(:) real(kind_r8), pointer :: zu(:) ! CANNOT be safely deallocated because `add_vert_coord` ! just uses pointers to point at it internally. real(kind_r8), pointer :: zw(:) ! CANNOT be safely deallocated because `add_vert_coord` @@ -157,7 +144,7 @@ subroutine init_reference_pressure() allocate(dzw(pver), stat=ierr) call check_allocate(ierr, subname, 'dzw(pver)', 'dyn_grid', __LINE__) - dzw(:) = 1.0_kind_r8 / rdzw + dzw(:) = 1.0_kind_r8 / real(rdzw(:), kind_r8) nullify(rdzw) @@ -228,15 +215,28 @@ end subroutine init_reference_pressure !> Provide grid and mapping information between global and local indexes to physics by calling `phys_grid_init`. !> (KCW, 2024-03-27) subroutine init_physics_grid() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_comp, only: mpas_dynamical_core, ncells_global, ncells_solve, sphere_radius + use dynconst, only: constant_pi => pi, rad_to_deg + use physics_column_type, only: kind_pcol, physics_column_t + use physics_grid, only: phys_grid_init + use spmd_utils, only: iam + use string_utils, only: stringify + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind + character(*), parameter :: subname = 'dyn_grid::init_physics_grid' character(max_hcoordname_len), allocatable :: dyn_attribute_name(:) integer :: hdim1_d, hdim2_d ! First and second horizontal dimensions of physics grid. integer :: i integer :: ierr - integer, pointer :: indextocellid(:) ! Global indexes of cell centers. - real(kind_r8), pointer :: areacell(:) ! Cell areas (square meters). - real(kind_r8), pointer :: latcell(:) ! Cell center latitudes (radians). - real(kind_r8), pointer :: loncell(:) ! Cell center longitudes (radians). + integer, pointer :: indextocellid(:) ! Global indexes of cell centers. + real(kind_dyn_mpas), pointer :: areacell(:) ! Cell areas (square meters). + real(kind_dyn_mpas), pointer :: latcell(:) ! Cell center latitudes (radians). + real(kind_dyn_mpas), pointer :: loncell(:) ! Cell center longitudes (radians). type(physics_column_t), allocatable :: dyn_column(:) ! Grid and mapping information between global and local indexes. nullify(areacell) @@ -261,12 +261,13 @@ subroutine init_physics_grid() ! Column information. dyn_column(i) % lat_rad = real(latcell(i), kind_pcol) dyn_column(i) % lon_rad = real(loncell(i), kind_pcol) - dyn_column(i) % lat_deg = real(latcell(i) * rad_to_deg, kind_pcol) - dyn_column(i) % lon_deg = real(loncell(i) * rad_to_deg, kind_pcol) + dyn_column(i) % lat_deg = real(real(latcell(i), kind_r8) * rad_to_deg, kind_pcol) + dyn_column(i) % lon_deg = real(real(loncell(i), kind_r8) * rad_to_deg, kind_pcol) ! Cell areas normalized to unit sphere. dyn_column(i) % area = real(areacell(i) / (sphere_radius ** 2), kind_pcol) ! Cell weights normalized to unity. - dyn_column(i) % weight = real(areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2), kind_pcol) + dyn_column(i) % weight = & + real(real(areacell(i), kind_r8) / (4.0_kind_r8 * constant_pi * real(sphere_radius, kind_r8) ** 2), kind_pcol) ! File decomposition. ! For unstructured grid, `coord_indices` is not used by `phys_grid_init`. @@ -308,19 +309,35 @@ end subroutine init_physics_grid !> * "mpas_vertex": Grid that is centered at MPAS "vertex" points. !> (KCW, 2024-03-28) subroutine define_cam_grid() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_grid_support, only: cam_grid_attribute_register, cam_grid_register, & + horiz_coord_create, horiz_coord_t + use cam_map_utils, only: kind_imap => imap + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & + ncells_global, nedges_global, nvertices_global, & + ncells_solve, nedges_solve, nvertices_solve, & + sphere_radius + use dynconst, only: constant_pi => pi, rad_to_deg + use string_utils, only: stringify + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind + character(*), parameter :: subname = 'dyn_grid::define_cam_grid' integer :: i integer :: ierr integer, pointer :: indextocellid(:) ! Global indexes of cell centers. integer, pointer :: indextoedgeid(:) ! Global indexes of edge nodes. integer, pointer :: indextovertexid(:) ! Global indexes of vertex nodes. - real(kind_r8), pointer :: areacell(:) ! Cell areas (square meters). - real(kind_r8), pointer :: latcell(:) ! Cell center latitudes (radians). - real(kind_r8), pointer :: latedge(:) ! Edge node latitudes (radians). - real(kind_r8), pointer :: latvertex(:) ! Vertex node latitudes (radians). - real(kind_r8), pointer :: loncell(:) ! Cell center longitudes (radians). - real(kind_r8), pointer :: lonedge(:) ! Edge node longitudes (radians). - real(kind_r8), pointer :: lonvertex(:) ! Vertex node longitudes (radians). + real(kind_dyn_mpas), pointer :: areacell(:) ! Cell areas (square meters). + real(kind_dyn_mpas), pointer :: latcell(:) ! Cell center latitudes (radians). + real(kind_dyn_mpas), pointer :: latedge(:) ! Edge node latitudes (radians). + real(kind_dyn_mpas), pointer :: latvertex(:) ! Vertex node latitudes (radians). + real(kind_dyn_mpas), pointer :: loncell(:) ! Cell center longitudes (radians). + real(kind_dyn_mpas), pointer :: lonedge(:) ! Edge node longitudes (radians). + real(kind_dyn_mpas), pointer :: lonvertex(:) ! Vertex node longitudes (radians). ! Global grid indexes. CAN be safely deallocated because its values are copied internally by ! `cam_grid_attribute_register` and `horiz_coord_create`. @@ -367,9 +384,9 @@ subroutine define_cam_grid() global_grid_index(:) = int(indextocellid(1:ncells_solve), kind_imap) lat_coord => horiz_coord_create('latCell', 'nCells', ncells_global, 'latitude', 'degrees_north', & - 1, ncells_solve, latcell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(latcell, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lonCell', 'nCells', ncells_global, 'longitude', 'degrees_east', & - 1, ncells_solve, loncell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(loncell, kind_r8) * rad_to_deg, map=global_grid_index) allocate(cell_area(ncells_solve), stat=ierr) call check_allocate(ierr, subname, 'cell_area(ncells_solve)', 'dyn_grid', __LINE__) @@ -379,8 +396,8 @@ subroutine define_cam_grid() call check_allocate(ierr, subname, 'global_grid_map(3, ncells_solve)', 'dyn_grid', __LINE__) do i = 1, ncells_solve - cell_area(i) = areacell(i) - cell_weight(i) = areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2) + cell_area(i) = real(areacell(i), kind_r8) + cell_weight(i) = real(areacell(i), kind_r8) / (4.0_kind_r8 * constant_pi * real(sphere_radius, kind_r8) ** 2) global_grid_map(1, i) = int(i, kind_imap) global_grid_map(2, i) = int(1, kind_imap) @@ -396,6 +413,7 @@ subroutine define_cam_grid() call cam_grid_attribute_register('mpas_cell', 'cell_weight', 'MPAS cell weight', 'nCells', cell_weight, & map=global_grid_index) + nullify(areacell) nullify(cell_area, cell_weight) nullify(lat_coord, lon_coord) @@ -403,16 +421,15 @@ subroutine define_cam_grid() ! Standard CAM-SIMA coordinate and dimension names are used. lat_coord => horiz_coord_create('lat', 'ncol', ncells_global, 'latitude', 'degrees_north', & - 1, ncells_solve, latcell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(latcell, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lon', 'ncol', ncells_global, 'longitude', 'degrees_east', & - 1, ncells_solve, loncell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(loncell, kind_r8) * rad_to_deg, map=global_grid_index) call dyn_debug_print('Registering grid "cam_cell" with id ' // stringify([dyn_grid_id('cam_cell')])) call cam_grid_register('cam_cell', dyn_grid_id('cam_cell'), lat_coord, lon_coord, global_grid_map, & unstruct=.true., block_indexed=.false.) - nullify(areacell) nullify(indextocellid) nullify(latcell, loncell) @@ -433,9 +450,9 @@ subroutine define_cam_grid() global_grid_index(:) = int(indextoedgeid(1:nedges_solve), kind_imap) lat_coord => horiz_coord_create('latEdge', 'nEdges', nedges_global, 'latitude', 'degrees_north', & - 1, nedges_solve, latedge * rad_to_deg, map=global_grid_index) + 1, nedges_solve, real(latedge, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lonEdge', 'nEdges', nedges_global, 'longitude', 'degrees_east', & - 1, nedges_solve, lonedge * rad_to_deg, map=global_grid_index) + 1, nedges_solve, real(lonedge, kind_r8) * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nedges_solve), stat=ierr) call check_allocate(ierr, subname, 'global_grid_map(3, nedges_solve)', 'dyn_grid', __LINE__) @@ -471,9 +488,9 @@ subroutine define_cam_grid() global_grid_index(:) = int(indextovertexid(1:nvertices_solve), kind_imap) lat_coord => horiz_coord_create('latVertex', 'nVertices', nvertices_global, 'latitude', 'degrees_north', & - 1, nvertices_solve, latvertex * rad_to_deg, map=global_grid_index) + 1, nvertices_solve, real(latvertex, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lonVertex', 'nVertices', nvertices_global, 'longitude', 'degrees_east', & - 1, nvertices_solve, lonvertex * rad_to_deg, map=global_grid_index) + 1, nvertices_solve, real(lonvertex, kind_r8) * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nvertices_solve), stat=ierr) call check_allocate(ierr, subname, 'global_grid_map(3, nvertices_solve)', 'dyn_grid', __LINE__) @@ -500,6 +517,9 @@ end subroutine define_cam_grid !> Helper function for returning grid id given its name. !> (KCW, 2024-03-27) pure function dyn_grid_id(name) + ! Module(s) from CAM-SIMA. + use physics_grid, only: phys_decomp + character(*), intent(in) :: name integer :: dyn_grid_id diff --git a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml index 874b6d86..c59570bf 100644 --- a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml +++ b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml @@ -1,710 +1,628 @@ - + - - - - char*256 - mpas - mpas_nhyd_model - - Time integration scheme in MPAS. - - Possible values: `SRK3' - Default: SRK3 - - - SRK3 - - - - integer - mpas - mpas_nhyd_model - - Order for RK time integration in MPAS. - - Possible values: 2 or 3 - Default: 2 - - - 2 - - - - real - mpas - mpas_nhyd_model - - Model time step, seconds in MPAS. - - Possible values: Positive real values - Default: 720.0 - - - 720.0 - - - - logical - mpas - mpas_nhyd_model - - Whether to super-cycle scalar transport in MPAS. - - Possible values: Logical values - Default: true - - - true - - - - integer - mpas - mpas_nhyd_model - - Number of acoustic steps per full RK step in MPAS. - - Possible values: Positive, even integer values, typically 2 or 6 depending - on transport splitting - Default: 2 - - - 2 - - - - integer - mpas - mpas_nhyd_model - - When config_split_dynamics_transport = T, the number of RK steps per - transport step in MPAS. - - Possible values: Positive integer values - Default: 3 - - - 3 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for horizontal diffusion of momentum in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of momentum in - MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for vertical diffusion of momentum in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for horizontal diffusion of theta in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of theta in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for vertical diffusion of theta in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - char*256 - mpas - mpas_nhyd_model - - Formulation of horizontal mixing in MPAS. - - Possible values: `2d_fixed' or `2d_smagorinsky' - Default: 2d_smagorinsky - - - 2d_smagorinsky - - - - real - mpas - mpas_nhyd_model - - Horizontal length scale, used by the Smagorinsky formulation of horizontal - diffusion and by 3-d divergence damping in MPAS. - - Possible values: Positive real values. A zero value implies that the - length scale is prescribed by the nominalMinDc value in the input file. - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - Scaling coefficient of $\delta x^3$ to obtain $\nabla^4$ diffusion - coefficient in MPAS. - - Possible values: Non-negative real values - Default: 0.05 - - - 0.05 - - - - real - mpas - mpas_nhyd_model - - Scaling factor for the divergent component of $\nabla^4 u$ calculation in - MPAS. - - Possible values: Positive real values - Default: 10.0 - - - 10.0 - - - - integer - mpas - mpas_nhyd_model - - Horizontal advection order for w in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Horizontal advection order for theta in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Horizontal advection order for scalars in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for normal velocities (u) in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for w in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for theta in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for scalars in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - logical - mpas - mpas_nhyd_model - - Whether to advect scalar fields in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - logical - mpas - mpas_nhyd_model - - Whether to enable positive-definite advection of scalars in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_nhyd_model - - Whether to enable monotonic limiter in scalar advection in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - real - mpas - mpas_nhyd_model - - Upwinding coefficient in the 3rd order advection scheme in MPAS. - - Possible values: 0 $\leq$ config_coef_3rd_order $\leq$ 1 - Default: 0.25 - - - 0.25 - - - - real - mpas - mpas_nhyd_model - - Dimensionless empirical parameter relating the strain tensor to the eddy - viscosity in the Smagorinsky turbulence model in MPAS. - - Possible values: Real values typically in the range 0.1 to 0.4 - Default: 0.125 - - - 0.125 - - - - logical - mpas - mpas_nhyd_model - - Mix full $\theta$ and $u$ fields, or mix perturbation from intitial state - in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - real - mpas - mpas_nhyd_model - - Off-centering parameter for the vertically implicit acoustic integration - in MPAS. - - Possible values: Positive real values - Default: 0.1 - - - 0.1 - - - - real - mpas - mpas_nhyd_model - - 3-d divergence damping coefficient in MPAS. - - Possible values: Positive real values - Default: 0.1 - - - 0.1 - - - - real - mpas - mpas_nhyd_model - - Amount of upwinding in APVM in MPAS. - - Possible values: 0 $\leq$ config_apvm_upwinding $\leq$ 1 - Default: 0.5 - - - 0.5 - - - - logical - mpas - mpas_nhyd_model - - Scale eddy viscosities with mesh-density function for horizontal diffusion - in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - real - mpas - mpas_nhyd_model - - Coefficient for the divergent component of the Laplacian filter of - momentum in the relaxation zone in MPAS. - - Possible values: Positive real values - Default: 6.0 - - - 6.0 - - - - real - mpas - mpas_damping - - Height MSL to begin w-damping profile in MPAS. - - Possible values: Positive real values - Default: 22000.0 - - - 22000.0 - - - - real - mpas - mpas_damping - - Maximum w-damping coefficient at model top in MPAS. - - Possible values: 0 $\leq$ config_xnutr $\leq$ 1 - Default: 0.2 - - - 0.2 - - - - real - mpas - mpas_damping - - Coefficient for scaling the 2nd-order horizontal mixing in the mpas_cam - absorbing layer in MPAS. - - Possible values: 0 $\leq$ config_mpas_cam_coef $\leq$ 1, standard value is - 0.2 - Default: 0.0 - - - 0.0 - - - - integer - mpas - mpas_damping - - Number of layers in which to apply cam 2nd-order horizontal filter top of - model; viscosity linearly ramps to zero by layer number from the top in - MPAS. - - Possible values: Positive integer values - Default: 4 - - - 4 - - - - logical - mpas - mpas_damping - - Whether to apply Rayleigh damping on horizontal velocity in the top-most - model levels. The number of levels is specified by the - config_number_rayleigh_damp_u_levels option, and the damping timescale is - specified by the config_rayleigh_damp_u_timescale_days option. in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - real - mpas - mpas_damping - - Timescale, in days (86400 s), for the Rayleigh damping on horizontal - velocity in the top-most model levels. in MPAS. - - Possible values: Positive real values - Default: 5.0 - - - 5.0 - - - - integer - mpas - mpas_damping - - Number of layers in which to apply Rayleigh damping on horizontal velocity - at top of model; damping linearly ramps to zero by layer number from the - top in MPAS. - - Possible values: Positive integer values - Default: 6 - - - 6 - - - - logical - mpas - mpas_limited_area - - Whether to apply lateral boundary conditions in MPAS. - - Possible values: true or false; this option must be set to true for - limited-area simulations and false for global simulations - Default: false - - - false - - - - char*256 - mpas - mpas_decomposition - - Prefix of graph decomposition file, to be suffixed with the MPI task count - in MPAS. - - Possible values: Any valid filename - Default: x1.40962.graph.info.part. - - - x1.40962.graph.info.part. - - - - logical - mpas - mpas_restart - - Whether this run of the model is to restart from a previous restart file - or not in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_printout - - Whether to print the global min/max of horizontal normal velocity and - vertical velocity each timestep in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - logical - mpas - mpas_printout - - Whether to print the global min/max of horizontal normal velocity and - vertical velocity each timestep, along with the location in the domain - where those extrema occurred in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_printout - - Whether to print the global min/max of scalar fields each timestep in - MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_assimilation - - Whether this run is within the JEDI data assimilation framework; used to - add temperature and specific humidity as diagnostics in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - + + + mpas + + Amount of upwinding in APVM + + mpas_nhyd_model + real + + 0.5 + + + + mpas + + Prefix of graph decomposition file, to be suffixed with the MPI task + count + + mpas_decomposition + char*256 + + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa30.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa12.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15-3.graph.info.part. + + + + mpas + + Coefficient for scaling the 2nd-order horizontal mixing in the + mpas_cam absorbing layer + + mpas_damping + real + + 0.0 + + + + mpas + + Upwinding coefficient in the 3rd order advection scheme + + mpas_nhyd_model + real + + 1.0 + + + + mpas + + Scaling factor for the divergent component of $\nabla^4 u$ + calculation + + mpas_nhyd_model + real + + 10.0 + + + + mpas + + Model time step, seconds + + mpas_nhyd_model + real + + 1800.0 + 900.0 + 450.0 + 225.0 + + + + mpas + + When config_split_dynamics_transport = T, the number of RK steps per + transport step + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Off-centering parameter for the vertically implicit acoustic + integration + + mpas_nhyd_model + real + + 0.1 + + + + mpas + + Whether to use an explicit mapping of blocks to MPI tasks + + mpas_decomposition + logical + + .false. + + + + mpas + + $\nabla^2$ eddy viscosity for horizontal diffusion of momentum + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of momentum + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + Scale eddy viscosities with mesh-density function for horizontal + diffusion + + mpas_nhyd_model + logical + + .true. + + + + mpas + + $\nabla^2$ eddy viscosity for horizontal diffusion of theta + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of theta + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + Method to use for exchanging halos + + mpas_development + char*256 + + mpas_halo + + + + mpas + + Formulation of horizontal mixing + + mpas_nhyd_model + char*256 + + 2d_smagorinsky + + + + mpas + + Horizontal length scale, used by the Smagorinsky formulation of + horizontal diffusion and by 3-d divergence damping + + mpas_nhyd_model + real + + 480000.0 + 120000.0 + 60000.0 + 30000.0 + + + + mpas + + Mix full $\theta$ and $u$ fields, or mix perturbation from intitial + state + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Whether to enable monotonic limiter in scalar advection + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Number of halo layers for fields + + mpas_nhyd_model + integer + + 2 + + + + mpas + + Number of layers in which to apply cam 2nd-order horizontal filter + top of model; viscosity linearly ramps to zero by layer number from + the top + + mpas_damping + integer + + 0 + + + + mpas + + Number of blocks to assign to each MPI task + + mpas_decomposition + integer + + 0 + + + + mpas + + Number of acoustic steps per full RK step + + mpas_nhyd_model + integer + + 2 + + + + mpas + + Number of layers in which to apply Rayleigh damping on horizontal + velocity at top of model; damping linearly ramps to zero by layer + number from the top + + mpas_damping + integer + + 5 + + + + mpas + + Number of tasks to perform file I/O + + mpas_io + integer + + 0 + + + + mpas + + Stride between file I/O tasks + + mpas_io + integer + + 1 + + + + mpas + + Whether to enable positive-definite advection of scalars + + mpas_nhyd_model + logical + + .false. + + + + mpas + + Whether to print the global min/max of horizontal normal velocity + and vertical velocity each timestep, along with the location in the + domain where those extrema occurred + + mpas_printout + logical + + .true. + + + + mpas + + Whether to print the global min/max of scalar fields each timestep + + mpas_printout + logical + + .false. + + + + mpas + + Whether to print the global min/max of horizontal normal velocity + and vertical velocity each timestep + + mpas_printout + logical + + .true. + + + + mpas + + Prefix of block mapping file + + mpas_decomposition + char*256 + + graph.info.part. + + + + mpas + + Whether to apply Rayleigh damping on horizontal velocity in the top- + most model levels. The number of levels is specified by the + config_number_rayleigh_damp_u_levels option, and the damping + timescale is specified by the config_rayleigh_damp_u_timescale_days + option + + mpas_damping + logical + + .true. + + + + mpas + + Timescale, in days (86400 s), for the Rayleigh damping on horizontal + velocity in the top-most model levels + + mpas_damping + real + + 5.0 + + + + mpas + + Coefficient for the divergent component of the Laplacian filter of + momentum in the relaxation zone + + mpas_nhyd_model + real + + 6.0 + + + + mpas + + Filename used to store most recent restart time stamp + + mpas_io + char*256 + + restart_timestamp + + + + mpas + + Horizontal advection order for scalars + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Whether to advect scalar fields + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Vertical advection order for scalars + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Dimensionless empirical parameter relating the strain tensor to the + eddy viscosity in the Smagorinsky turbulence model + + mpas_nhyd_model + real + + 0.125 + + + + mpas + + 3-d divergence damping coefficient + + mpas_nhyd_model + real + + 0.1 + + + + mpas + + Whether to super-cycle scalar transport + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Horizontal advection order for theta + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Vertical advection order for theta + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Time integration scheme + + mpas_nhyd_model + char*256 + + SRK3 + + + + mpas + + Order for RK time integration + + mpas_nhyd_model + integer + + 2 + + + + mpas + + Vertical advection order for normal velocities (u) + + mpas_nhyd_model + integer + + 3 + + + + mpas + + $\nabla^2$ eddy viscosity for vertical diffusion of momentum + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + $\nabla^2$ eddy viscosity for vertical diffusion of theta + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + Scaling coefficient of $\delta x^3$ to obtain $\nabla^4$ diffusion + coefficient + + mpas_nhyd_model + real + + 0.05 + + + + mpas + + Horizontal advection order for w + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Vertical advection order for w + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Maximum w-damping coefficient at model top + + mpas_damping + real + + 0.2 + + + + mpas + + Height MSL to begin w-damping profile + + mpas_damping + real + + 22000.0 + + diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index d089636e..a89c7c75 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -1,15 +1,4 @@ module stepon - ! Modules from CAM-SIMA. - use camsrfexch, only: cam_out_t - use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run - use dyn_coupling, only: dynamics_to_physics_coupling, physics_to_dynamics_coupling - use physics_types, only: physics_state, physics_tend - use runtime_obj, only: runtime_options - use time_manager, only: get_step_size - - ! Modules from CCPP. - use ccpp_kinds, only: kind_phys - implicit none private @@ -22,6 +11,10 @@ module stepon contains ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine stepon_init(cam_runtime_opts, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t + use runtime_obj, only: runtime_options + type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out @@ -29,6 +22,15 @@ end subroutine stepon_init ! Called by `cam_timestep_init` in `src/control/cam_comp.F90`. subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t + use dyn_coupling, only: dynamics_to_physics_coupling + use physics_types, only: physics_state, physics_tend + use runtime_obj, only: runtime_options + use time_manager, only: get_step_size + ! Module(s) from CCPP. + use ccpp_kinds, only: kind_phys + real(kind_phys), intent(out) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts type(physics_state), intent(in) :: phys_state @@ -44,6 +46,12 @@ end subroutine stepon_timestep_init ! Called by `cam_run2` in `src/control/cam_comp.F90`. subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t + use dyn_coupling, only: physics_to_dynamics_coupling + use physics_types, only: physics_state, physics_tend + use runtime_obj, only: runtime_options + type(runtime_options), intent(in) :: cam_runtime_opts type(physics_state), intent(in) :: phys_state type(physics_tend), intent(in) :: phys_tend @@ -55,6 +63,14 @@ end subroutine stepon_run2 ! Called by `cam_run3` in `src/control/cam_comp.F90`. subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use camsrfexch, only: cam_out_t + use dyn_comp, only: dyn_export_t, dyn_import_t, dyn_run + use physics_types, only: physics_state + use runtime_obj, only: runtime_options + ! Module(s) from CCPP. + use ccpp_kinds, only: kind_phys + real(kind_phys), intent(in) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts type(cam_out_t), intent(in) :: cam_out @@ -67,8 +83,14 @@ end subroutine stepon_run3 ! Called by `cam_final` in `src/control/cam_comp.F90`. subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t, dyn_final + use runtime_obj, only: runtime_options + type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out + + call dyn_final() end subroutine stepon_final end module stepon diff --git a/test/existing-test-failures.txt b/test/existing-test-failures.txt index 975b8681..2328b59d 100644 --- a/test/existing-test-failures.txt +++ b/test/existing-test-failures.txt @@ -1,3 +1 @@ -SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_intel.cam-outfrq_kessler_mpas_derecho (Overall: FAIL) -SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_gnu.cam-outfrq_kessler_mpas_derecho (Overall: FAIL) - - will fail until MPAS is fully integrated +All clear!