diff --git a/.gitmodules b/.gitmodules
index 7ace1182e349..bf9030ea873d 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -70,6 +70,9 @@
path = components/ww3/src/WW3
url = git@github.com:E3SM-Project/WW3.git
branch = e3sm
+[submodule "components/eam/src/physics/crm/pam/external"]
+ path = components/eam/src/physics/crm/pam/external
+ url = git@github.com:E3SM-Project/PAM.git
[submodule "externals/haero"]
path = externals/haero
url = git@github.com:eagles-project/haero.git
diff --git a/cime_config/machines/cmake_macros/gnu.cmake b/cime_config/machines/cmake_macros/gnu.cmake
index eae59e3e4b9f..3325f9966344 100644
--- a/cime_config/machines/cmake_macros/gnu.cmake
+++ b/cime_config/machines/cmake_macros/gnu.cmake
@@ -1,6 +1,9 @@
if (CMAKE_SYSTEM_PROCESSOR MATCHES "^aarch64")
string(APPEND CFLAGS "-mcmodel=small")
string(APPEND FFLAGS "-mcmodel=small")
+elseif (CMAKE_SYSTEM_PROCESSOR MATCHES "arm64")
+ string(APPEND CFLAGS " -mcmodel=large")
+ string(APPEND FFLAGS " -mcmodel=large")
else()
string(APPEND CFLAGS "-mcmodel=medium")
string(APPEND FFLAGS "-mcmodel=medium")
diff --git a/cime_config/machines/cmake_macros/gnugpu_frontier.cmake b/cime_config/machines/cmake_macros/gnugpu_frontier.cmake
index 8cb67348c46e..31bc97afd0e3 100644
--- a/cime_config/machines/cmake_macros/gnugpu_frontier.cmake
+++ b/cime_config/machines/cmake_macros/gnugpu_frontier.cmake
@@ -21,7 +21,7 @@ string(APPEND CMAKE_OPTS " -DPIO_ENABLE_TOOLS:BOOL=OFF")
string(APPEND CXXFLAGS " -I$ENV{MPICH_DIR}/include --offload-arch=gfx90a")
string(APPEND SLIBS " -Wl,--copy-dt-needed-entries -L/opt/cray/pe/gcc-libs -lgfortran -L$ENV{MPICH_DIR}/lib -lmpi -L$ENV{CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa ")
-string(APPEND KOKKOS_OPTIONS " -DKokkos_ENABLE_HIP=On -DKokkos_ARCH_ZEN3=On -DKokkos_ARCH_VEGA90A=On")
+string(APPEND KOKKOS_OPTIONS " -DKokkos_ENABLE_HIP=On -DKokkos_ARCH_ZEN3=On -DKokkos_ARCH_VEGA90A=On -DKokkos_ENABLE_OPENMP=Off")
set(USE_HIP "TRUE")
string(APPEND HIP_FLAGS "${CXXFLAGS} -munsafe-fp-atomics -x hip")
diff --git a/cime_config/machines/cmake_macros/gnugpu_pm-gpu.cmake b/cime_config/machines/cmake_macros/gnugpu_pm-gpu.cmake
index b362b984082e..b696962159cf 100644
--- a/cime_config/machines/cmake_macros/gnugpu_pm-gpu.cmake
+++ b/cime_config/machines/cmake_macros/gnugpu_pm-gpu.cmake
@@ -7,6 +7,7 @@ endif()
string(APPEND CPPDEFS " -DTHRUST_IGNORE_CUB_VERSION_CHECK")
string(APPEND SLIBS " -lblas -llapack")
string(APPEND CUDA_FLAGS " -ccbin CC -O2 -arch sm_80 --use_fast_math")
+string(APPEND KOKKOS_OPTIONS " -DKokkos_ARCH_AMPERE80=On -DKokkos_ENABLE_CUDA=On -DKokkos_ENABLE_CUDA_LAMBDA=On -DKokkos_ENABLE_SERIAL=ON -DKokkos_ENABLE_OPENMP=Off")
if (NOT DEBUG)
string(APPEND CFLAGS " -O2")
string(APPEND FFLAGS " -O2")
diff --git a/cime_config/machines/cmake_macros/gnugpu_summit.cmake b/cime_config/machines/cmake_macros/gnugpu_summit.cmake
index 52faa0d641ad..bcc77264128b 100644
--- a/cime_config/machines/cmake_macros/gnugpu_summit.cmake
+++ b/cime_config/machines/cmake_macros/gnugpu_summit.cmake
@@ -14,3 +14,4 @@ string(APPEND SLIBS " -L$ENV{ESSL_PATH}/lib64 -lessl -L$ENV{OLCF_NETLIB_LAPACK_R
set(MPICXX "mpiCC")
set(PIO_FILESYSTEM_HINTS "gpfs")
set(USE_CUDA "TRUE")
+string(APPEND KOKKOS_OPTIONS " -DKokkos_ARCH_POWER9=On -DKokkos_ARCH_VOLTA70=On -DKokkos_ENABLE_CUDA=On -DKokkos_ENABLE_CUDA_LAMBDA=On -DKokkos_ENABLE_OPENMP=Off")
diff --git a/cime_config/tests.py b/cime_config/tests.py
index 3ed0e4d59687..cccc5a56f534 100644
--- a/cime_config/tests.py
+++ b/cime_config/tests.py
@@ -339,6 +339,7 @@
"ERP_Ln9.ne4pg2_oQU480.WCYCL20TRNS-MMF1.allactive-mmf_fixed_subcycle",
"ERS_Ln9.ne4pg2_ne4pg2.FRCE-MMF1.eam-cosp_nhtfrq9",
"SMS_Ln5.ne4_ne4.FSCM-ARM97-MMF1",
+ "SMS_Ln3.ne4pg2_ne4pg2.F2010-MMF2",
)
},
diff --git a/components/cmake/build_model.cmake b/components/cmake/build_model.cmake
index 1615fba9bb03..0049dad3f6d9 100644
--- a/components/cmake/build_model.cmake
+++ b/components/cmake/build_model.cmake
@@ -116,6 +116,22 @@ function(build_model COMP_CLASS COMP_NAME)
cmake/atm/../../eam/src/physics/crm/crm_ecpp_output_module.F90 )
endif()
+ if (USE_PAM)
+ message(STATUS "Building PAM")
+ # PAM_HOME is where the samxx source code lives
+ set(PAM_HOME ${CMAKE_CURRENT_SOURCE_DIR}/../../eam/src/physics/crm/pam)
+ # PAM_BIN is where the samxx library will live
+ set(PAM_BIN ${CMAKE_CURRENT_BINARY_DIR}/pam)
+ # Build the static samxx library
+ add_subdirectory(${PAM_HOME} ${PAM_BIN})
+ # Add samxx F90 files to the main E3SM build
+ set(SOURCES ${SOURCES} cmake/atm/../../eam/src/physics/crm/pam/params.F90
+ cmake/atm/../../eam/src/physics/crm/crm_ecpp_output_module.F90
+ cmake/atm/../../eam/src/physics/crm/pam/pam_driver.F90)
+ # Pam interface need to include modules from pam
+ include_directories(${PAM_BIN}/external/pam_core)
+ endif()
+
# Add rrtmgp++ source code if asked for
if (USE_RRTMGPXX)
message(STATUS "Building RRTMGPXX")
@@ -245,6 +261,9 @@ function(build_model COMP_CLASS COMP_NAME)
if (USE_SAMXX)
target_link_libraries(${TARGET_NAME} PRIVATE samxx)
endif()
+ if (USE_PAM)
+ target_link_libraries(${TARGET_NAME} PRIVATE pam_driver)
+ endif()
if (USE_RRTMGPXX)
target_link_libraries(${TARGET_NAME} PRIVATE rrtmgp rrtmgp_interface)
endif()
diff --git a/components/cmake/common_setup.cmake b/components/cmake/common_setup.cmake
index 2290ec8c715d..67c583f94888 100644
--- a/components/cmake/common_setup.cmake
+++ b/components/cmake/common_setup.cmake
@@ -35,6 +35,14 @@ if (NOT HAS_SAMXX EQUAL -1)
set(USE_SAMXX TRUE)
endif()
+# Look for -crm pam in the CAM_CONFIG_OPTS CIME variable
+# If it's found, then enable USE_PAM
+string(FIND "${CAM_CONFIG_OPTS}" "-crm pam" HAS_PAM)
+if (NOT HAS_PAM EQUAL -1)
+ # The following is for the PAM code:
+ set(USE_PAM TRUE)
+endif()
+
string(FIND "${CAM_CONFIG_OPTS}" "-rrtmgpxx" HAS_RRTMGPXX)
if (NOT HAS_RRTMGPXX EQUAL -1)
# The following is for the RRTMGPXX code:
@@ -42,7 +50,7 @@ if (NOT HAS_RRTMGPXX EQUAL -1)
endif()
# If samxx or rrtmgpxx is being used, then YAKL must be used as well
-if (USE_SAMXX OR USE_RRTMGPXX)
+if (USE_SAMXX OR USE_RRTMGPXX OR USE_PAM)
set(USE_YAKL TRUE)
else()
set(USE_YAKL FALSE)
diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist
index 90772b100042..fccfb9f4011b 100755
--- a/components/eam/bld/build-namelist
+++ b/components/eam/bld/build-namelist
@@ -1705,13 +1705,9 @@ if ( !($phys eq 'ideal' or $phys eq 'adiabatic') ) {
# Cloud optics; for MMF be specific about which optics package to use
if ($rad_pkg eq 'rrtmg' or $rad_pkg eq 'rrtmgp') {
if ($cfg->get('use_MMF')) {
- if ($cfg->get('MMF_microphysics_scheme') eq 'sam1mom') {
- add_default($nl, 'liqcldoptics', 'val' => 'slingo');
- add_default($nl, 'icecldoptics', 'val' => 'ebertcurry');
- } else {
- add_default($nl, 'liqcldoptics', 'val' => 'gammadist');
- add_default($nl, 'icecldoptics', 'val' => 'mitchell');
- }
+ # NOTE: eventually we should switch PAM optics to gammadist+mitchell
+ add_default($nl, 'liqcldoptics', 'val' => 'slingo');
+ add_default($nl, 'icecldoptics', 'val' => 'ebertcurry');
} else {
add_default($nl, 'liqcldoptics');
add_default($nl, 'icecldoptics');
diff --git a/components/eam/bld/config_files/definition.xml b/components/eam/bld/config_files/definition.xml
index dcaf803a647c..2cbd3d33ebe7 100644
--- a/components/eam/bld/config_files/definition.xml
+++ b/components/eam/bld/config_files/definition.xml
@@ -345,15 +345,15 @@ Switch to enable explicit scalar momentum transport for 2D CRM in E3SM-MMF: 0=of
Switch to enable or disable ECPP in CAM: 0=off, 1=on.
-
-Switch to enable or disable fractional cloudiness in m2005 microphysics in CRM: 0=off, 1=on.
-
-
+
MMF CRM cloud microphysics scheme
-
+
MMF CRM model
+
+MMF PAM dycor option
+
MMF CRM advection scheme
diff --git a/components/eam/bld/configure b/components/eam/bld/configure
index 15ec67275278..ecd7561c72a0 100755
--- a/components/eam/bld/configure
+++ b/components/eam/bld/configure
@@ -272,20 +272,24 @@ OPTIONS
Options for MMF:
- -use_MMF Build the MMF configuration (super-parameterized EAM)
- -use_MMF_VT enable CRM variance transport (VT)
- -use_MMF_ESMT enable CRM explicit scalar momentum transport (ESMT)
- -crm_nx CRM's x-grid.
- -crm_ny CRM's y-grid.
- -crm_nz CRM's z-grid.
- -crm_dx CRM's horizontal grid spacing.
- -crm_dt CRM's timestep.
- -crm_nx_rad number of averaged columns for CRM radiation calculation in x direction
- -crm_ny_rad number of averaged columns for CRM radiation calculation in y direction
- -MMF_microphysics_scheme CRM microphysics package name [sam1mom | m2005 ].
- -use_ECPP use CRM clouds for vertical transport, aqueous chemistry and wet removable of aerosols
- -crm_adv CRM advection scheme [MPDATA | UM5]
- -crm CRM model [sam | samomp | samxx]
+ -use_MMF Build the MMF configuration (super-parameterized EAM)
+ -use_MMF_VT enable CRM variance transport (VT)
+ -use_MMF_ESMT enable CRM explicit scalar momentum transport (ESMT)
+ -use_ECPP enable explicit-clouds-parameterized-pollutants scheme
+ -crm_nx CRM's x-grid.
+ -crm_ny CRM's y-grid.
+ -crm_nz CRM's z-grid.
+ -crm_dx CRM's horizontal grid spacing.
+ -crm_dt CRM's timestep.
+ -crm_nx_rad # of columns for CRM radiation calculation in x direction
+ -crm_ny_rad # of columns for CRM radiation calculation in y direction
+ -MMF_microphysics_scheme CRM micro scheme - depends on CRM model [ sam1mom | p3 ].
+ -crm MMF CRM model [ sam | samxx | pam]
+ -crm_adv SAM CRM advection scheme [ MPDATA | UM5]
+ -pam_dycor PAM CRM dycor option [ awfl | spam ]
+
+ Options for radiation:
+
-rrtmgpxx Use RRTMGP++ code
EOF
}
@@ -358,9 +362,10 @@ GetOptions(
"crm_dt=s" => \$opts{'crm_dt'},
"crm_nx_rad=s" => \$opts{'crm_nx_rad'},
"crm_ny_rad=s" => \$opts{'crm_ny_rad'},
- "MMF_microphysics_scheme=s" => \$opts{'MMF_microphysics_scheme'},
+ "MMF_microphysics_scheme=s" => \$opts{'MMF_microphysics_scheme'},
"crm_adv=s" => \$opts{'crm_adv'},
"crm=s" => \$opts{'crm'},
+ "pam_dycor=s" => \$opts{'pam_dycor'},
"rrtmgpxx" => \$opts{'rrtmgpxx'},
"debug" => \$opts{'debug'},
"rain_evap_to_coarse_aero" => \$opts{'rain_evap_to_coarse_aero'},
@@ -1293,6 +1298,7 @@ if (defined $opts{'use_MMF'}) {
if (defined $opts{'crm_nx_rad'}) { $cfg_ref->set('crm_nx_rad', $opts{'crm_nx_rad'}); }
if (defined $opts{'crm_ny_rad'}) { $cfg_ref->set('crm_ny_rad', $opts{'crm_ny_rad'}); }
if (defined $opts{'crm_adv'}) { $cfg_ref->set('crm_adv', $opts{'crm_adv'}); }
+ if (defined $opts{'pam_dycor'}) { $cfg_ref->set('pam_dycor', $opts{'pam_dycor'}); }
$cfg_ref->set('MMF_microphysics_scheme', $opts{'MMF_microphysics_scheme'});
$cfg_ref->set('crm', $opts{'crm'});
if (defined $opts{'use_MMF_VT'}) { $cfg_ref->set('use_MMF_VT', 1); }
@@ -1812,12 +1818,16 @@ if (defined $opts{'use_MMF'}) {
my $crm_nx_rad = $cfg_ref->get('crm_nx_rad');
my $crm_ny_rad = $cfg_ref->get('crm_ny_rad');
my $crm = $cfg_ref->get('crm');
+ my $pam_dycor = $cfg_ref->get('pam_dycor');
my $MMF_microphysics_scheme = $cfg_ref->get('MMF_microphysics_scheme');
my $crm_adv = $cfg_ref->get('crm_adv');
my $yes3Dval = 1;
if ($crm_ny eq 1) {$yes3Dval = 0;}
if (defined $opts{'crm_adv'}) { $cfg_cppdefs .= " -D_$crm_adv "; }
- $cfg_cppdefs .= " -D$MMF_microphysics_scheme -DYES3DVAL=$yes3Dval";
+ $cfg_cppdefs .= " -DYES3DVAL=$yes3Dval";
+ if ($MMF_microphysics_scheme eq 'sam1mom') {
+ $cfg_cppdefs .= " -D$MMF_microphysics_scheme ";
+ }
$cfg_cppdefs .= " -DCRM_NX=$crm_nx -DCRM_NY=$crm_ny -DCRM_NZ=$crm_nz ";
$cfg_cppdefs .= " -DCRM_DX=$crm_dx -DCRM_DT=$crm_dt ";
if (defined $opts{'use_ECPP'}) { $cfg_cppdefs .= " -DECPP " }
@@ -1832,13 +1842,14 @@ if (defined $opts{'use_MMF'}) {
$cfg_cppdefs .= " -DCRM_NY_RAD=$crm_ny "
}
if ($crm_ny eq 1) {$yes3Dval = 0;}
- my $crm = $cfg_ref->get('crm');
if ($crm eq 'sam') {
$cfg_cppdefs .= " -DMMF_SAM "
- } elsif ($crm eq 'samomp') {
- $cfg_cppdefs .= " -DMMF_SAMOMP"
} elsif ($crm eq 'samxx') {
$cfg_cppdefs .= " -DMMF_SAMXX "
+ } elsif ($crm eq 'pam') {
+ $cfg_cppdefs .= " -DMMF_PAM ";
+ if ($pam_dycor eq 'awfl') { $cfg_cppdefs .= " -DMMF_PAM_DYCOR_AWFL " }
+ if ($pam_dycor eq 'spam') { $cfg_cppdefs .= " -DMMF_PAM_DYCOR_SPAM " }
}
}
@@ -2555,6 +2566,7 @@ sub write_filepath
my $ice = $cfg_ref->get('ice');
my $rof = $cfg_ref->get('rof');
my $caseroot = $cfg_ref->get('caseroot');
+ my $crm = $cfg_ref->get('crm');
# Root directory
my $camsrcdir = "$cam_root/components";
@@ -2589,6 +2601,7 @@ sub write_filepath
# interface and base model code
print $fh "$camsrcdir/eam/src/physics/crm\n";
+ print $fh "$camsrcdir/eam/src/physics/crm/dummy_modules\n";
# MMF-specific radiation drivers (overrides default drivers that exist
# in each individual radiation package; done this way to avoid
@@ -2603,7 +2616,6 @@ sub write_filepath
print $fh "$camsrcdir/eam/src/physics/crm/ecpp\n";
}
- my $crm = $cfg_ref->get('crm');
if ($crm eq 'sam') {
#####################################################
## -crm sam will GLOB directories for CMake
@@ -2614,41 +2626,9 @@ sub write_filepath
my $MMF_microphysics_scheme = $cfg_ref->get('MMF_microphysics_scheme');
if ($MMF_microphysics_scheme eq 'sam1mom') {
print $fh "$camsrcdir/eam/src/physics/crm/sam/MICRO_SAM1MOM\n";
- } elsif ($MMF_microphysics_scheme eq 'm2005') {
- print $fh "$camsrcdir/eam/src/physics/crm/sam/MICRO_M2005\n";
}
-
# turbulence closure
print $fh "$camsrcdir/eam/src/physics/crm/sam/SGS_TKE\n";
-
- # advection scheme
- my $crm_adv = 'MPDATA';
- if (defined $opts{'crm_adv'}) { $crm_adv = $cfg_ref->get('crm_adv'); }
- if ($crm_adv eq 'MPDATA') {
- print $fh "$camsrcdir/eam/src/physics/crm/sam/ADV_MPDATA\n";
- } elsif ($crm_adv eq 'UM5') {
- print $fh "$camsrcdir/eam/src/physics/crm/sam/ADV_UM5\n";
- } else {
- die "configure ERROR: invalid option - crm_adv: $crm_adv.\n";
- }
-
- } elsif ($crm eq 'samomp') {
- #####################################################
- ## -crm samomp will GLOB directories for CMake
- ######################################################
- print $fh "$camsrcdir/eam/src/physics/crm/samomp\n";
- #
- # microphysics
- my $MMF_microphysics_scheme = $cfg_ref->get('MMF_microphysics_scheme');
- if ($MMF_microphysics_scheme eq 'sam1mom') {
- print $fh "$camsrcdir/eam/src/physics/crm/samomp/MICRO_SAM1MOM\n";
- } elsif ($MMF_microphysics_scheme eq 'm2005') {
- print $fh "$camsrcdir/eam/src/physics/crm/samomp/MICRO_M2005\n";
- }
-
- # turbulence closure
- print $fh "$camsrcdir/eam/src/physics/crm/samomp/SGS_TKE\n";
-
# advection scheme
my $crm_adv = 'MPDATA';
if (defined $opts{'crm_adv'}) { $crm_adv = $cfg_ref->get('crm_adv'); }
@@ -2659,12 +2639,6 @@ sub write_filepath
} else {
die "configure ERROR: invalid option - crm_adv: $crm_adv.\n";
}
-
- } elsif ($crm eq 'samxx') {
- #####################################################
- ## -crm samxx will add_subdirectory in CMake rather
- ## than GLOBing.
- #####################################################
}
}
@@ -2725,11 +2699,12 @@ sub write_filepath
print $fh "$camsrcdir/eam/src/physics/silhs\n";
}
- # set Fortran P3 path, needed even if microphys is not P3
- if ($is_scream_config eq '1') {
- print $fh "$camsrcdir/eam/src/physics/p3/scream\n";
- } else {
- print $fh "$camsrcdir/eam/src/physics/p3/eam\n";
+ if (not(defined $opts{'use_MMF'})) {
+ if ($is_scream_config eq '1') {
+ print $fh "$camsrcdir/eam/src/physics/p3/scream\n";
+ } else {
+ print $fh "$camsrcdir/eam/src/physics/p3/eam\n";
+ }
}
print $fh "$camsrcdir/eam/src/dynamics/$dyn\n";
@@ -2867,6 +2842,7 @@ sub write_cppdefs
my $fh = new IO::File;
my $dyn = $cfg_ref->get('dyn');
+ my $crm = $cfg_ref->get('crm');
print "Writing CPPDEFS for EAM in $file, with $dyn and $opts{'dyn_target'}\n";
diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml
index 4a1bcb2ed2a5..68d315aa0cf0 100755
--- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml
+++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml
@@ -869,6 +869,7 @@
sam1mom
m2005
+p3
.false.
.true.
.false.
@@ -1143,15 +1144,21 @@
.false.
-.false.
-.false.
-0
-.true.
-.true.
-
-.false.
-.false.
-2
+.false.
+.false.
+0
+2
+.true.
+.true.
+2
+.true.
+.true.
+2
+.true.
+.true.
+0
+.false.
+.false.
.true.
@@ -1511,7 +1518,6 @@
4
4
-
- -use_MMF -crm samxx -nlev 60 -crm_nz 50
- -crm_dt 10 -crm_dx 2000 -crm_nx 64 -crm_ny 1
- -crm_nx_rad 4 -crm_ny_rad 1 -rad rrtmgp -rrtmgpxx
+ -crm samxx -crm_dt 10
+ -crm pam -pam_dycor spam -crm_dt 8
+ -crm pam -pam_dycor awfl -crm_dt 8
+ -use_MMF -nlev 60 -crm_nz 50
+ -crm_dx 2000 -crm_nx 64 -crm_ny 1
-MMF_microphysics_scheme sam1mom -chem none
- -use_MMF_VT
- -use_MMF_ESMT
+ -MMF_microphysics_scheme p3 -chem none
+ -crm_nx_rad 4 -crm_ny_rad 1 -rad rrtmgp -rrtmgpxx
+ -use_MMF_VT
+ -use_MMF_ESMT
-aquaplanet
-aquaplanet -rce
-scam
@@ -141,10 +145,10 @@
aquaplanet_SCREAM-LR
aquaplanet_SCREAM-HR
- 20TR_MMF-1mom_CMIP6
- 1950_MMF-1mom_CMIP6
- 2010_mmf_1mom
- 2010_mmf_2mom
+
+ 20TR_MMF-1mom_CMIP6
+ 1950_MMF-1mom_CMIP6
+ 2010_mmf_1mom
aquaplanet_MMF-1mom
RCEMIP_EAMv1
scam_arm97_MMF-1mom
diff --git a/components/eam/cime_config/config_compsets.xml b/components/eam/cime_config/config_compsets.xml
index 442e545ab530..314f9a147945 100644
--- a/components/eam/cime_config/config_compsets.xml
+++ b/components/eam/cime_config/config_compsets.xml
@@ -136,11 +136,26 @@
ARM97_EAM%SCM-MMF1_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+ FSCM-ARM97-MMF2
+ ARM97_EAM%SCM-MMF2_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+
+
+ FSCM-ARM97-MMF2-AWFL
+ ARM97_EAM%SCM-MMF2-AWFL_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+
FSCM-RICO-MMF1
RICO_EAM%SCM-MMF1_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+ FSCM-RICO-MMF2
+ RICO_EAM%SCM-MMF2_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+
@@ -169,6 +184,11 @@
FAQP-MMF1
2000_EAM%AQP-MMF1_SLND_SICE_DOCN%AQP1_SROF_SGLC_SWAV
+
+
+ FAQP-MMF2
+ 2000_EAM%AQP-MMF2_SLND_SICE_DOCN%AQP1_SROF_SGLC_SWAV
+
@@ -254,6 +274,11 @@
2000_EAM%RCE-MMF1_SLND_SICE_DOCN%AQPCONST_SROF_SGLC_SWAV
+
+ FRCE-MMF2
+ 2000_EAM%RCE-MMF2_SLND_SICE_DOCN%AQPCONST_SROF_SGLC_SWAV
+
+
@@ -284,12 +309,27 @@
2010_EAM%MMF1_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+ F2010-MMF2
+ 2010_EAM%MMF2_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+
+
+ F2010-MMF2-AWFL
+ 2010_EAM%MMF2-AWFL_ELM%SP_CICE%PRES_DOCN%DOM_SROF_SGLC_SWAV
+
+
F20TR-MMF1
20TR_EAM%CMIP6-MMF1_ELM%SP_MPASSI%PRES_DOCN%DOM_MOSART_SGLC_SWAV
+
+ F20TR-MMF1
+ 20TR_EAM%CMIP6-MMF2_ELM%SP_MPASSI%PRES_DOCN%DOM_MOSART_SGLC_SWAV
+
+
diff --git a/components/eam/src/control/cam_comp.F90 b/components/eam/src/control/cam_comp.F90
index 940e9fdea51d..ddb453e36073 100644
--- a/components/eam/src/control/cam_comp.F90
+++ b/components/eam/src/control/cam_comp.F90
@@ -262,7 +262,7 @@ subroutine cam_run1(cam_in, cam_out)
!
call t_barrierf ('sync_phys_run1', mpicom)
call t_startf ('phys_run1')
-#if defined(MMF_SAMXX)
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
call phys_run1(phys_state, dtime, phys_tend, pbuf2d, cam_in, cam_out)
#else
call phys_run1(phys_state, dtime, phys_tend, pbuf2d, cam_in, cam_out, phys_diag)
@@ -301,7 +301,7 @@ subroutine cam_run2( cam_out, cam_in )
!
call t_barrierf ('sync_phys_run2', mpicom)
call t_startf ('phys_run2')
-#if defined(MMF_SAMXX)
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
call phys_run2(phys_state, dtime, phys_tend, pbuf2d, cam_out, cam_in)
#else
call phys_run2(phys_state, dtime, phys_tend, pbuf2d, cam_out, cam_in, phys_diag)
@@ -484,7 +484,7 @@ subroutine cam_final( cam_out, cam_in )
#endif
call t_startf ('phys_final')
-#if defined(MMF_SAMXX)
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
call phys_final( phys_state, phys_tend , pbuf2d )
#else
call phys_final( phys_state, phys_tend , pbuf2d, phys_diag )
diff --git a/components/eam/src/physics/cam/phys_control.F90 b/components/eam/src/physics/cam/phys_control.F90
index c933305f8208..c4d4b7fba7fa 100644
--- a/components/eam/src/physics/cam/phys_control.F90
+++ b/components/eam/src/physics/cam/phys_control.F90
@@ -491,7 +491,7 @@ subroutine phys_ctl_readnl(nlfile)
! check MMF parameters
if (use_MMF) then
! Check settings for MMF_microphysics_scheme
- if ( .not.(MMF_microphysics_scheme .eq. 'm2005' .or. &
+ if ( .not.(MMF_microphysics_scheme .eq. 'p3' .or. &
MMF_microphysics_scheme .eq. 'sam1mom' )) then
write(iulog,*)'phys_setopts: illegal value of MMF_microphysics_scheme:', MMF_microphysics_scheme
call endrun('phys_setopts: illegal value of MMF_microphysics_scheme')
diff --git a/components/eam/src/physics/crm/crm_history.F90 b/components/eam/src/physics/crm/crm_history.F90
index e7534c2d3adc..755f71d74569 100644
--- a/components/eam/src/physics/crm/crm_history.F90
+++ b/components/eam/src/physics/crm/crm_history.F90
@@ -85,6 +85,7 @@ subroutine crm_history_init(species_class)
character(8) :: unit
#endif
+ character(len=10),dimension(3) :: dims_rad_3D = (/'crm_nx_rad','crm_ny_rad','crm_nz '/)
character(len=6), dimension(3) :: dims_crm_3D = (/'crm_nx','crm_ny','crm_nz'/)
character(len=6), dimension(2) :: dims_crm_2D = (/'crm_nx','crm_ny'/)
@@ -137,12 +138,26 @@ subroutine crm_history_init(species_class)
call addfld('MMF_NC ',(/'lev'/), 'A', '/kg', 'Cloud water dropet number from CRM')
call addfld('MMF_NI ',(/'lev'/), 'A', '/kg', 'Cloud ice crystal number from CRM')
call addfld('MMF_NR ',(/'lev'/), 'A', '/kg', 'Rain particle number from CRM')
+ call addfld('MMF_RHOD',(/'lev'/), 'A', '/kg', 'Dry density from CRM')
+ call addfld('MMF_RHOV',(/'lev'/), 'A', '/kg', 'Vapor density from CRM')
call addfld('CRM_QR ',dims_crm_3D, 'A', 'kg/kg','Rain mixing ratio from CRM' )
call addfld('CRM_NC ',dims_crm_3D, 'A', '/kg', 'Cloud water dropet number from CRM' )
call addfld('CRM_NI ',dims_crm_3D, 'A', '/kg', 'Cloud ice crystal number from CRM' )
call addfld('CRM_NR ',dims_crm_3D, 'A', '/kg', 'Rain particle number from CRM' )
+
+ call addfld('MMF_LIQ_ICE',(/'lev'/),'A', '/kg', 'liq-ice phase change tendency from P3')
+ call addfld('MMF_VAP_LIQ',(/'lev'/),'A', '/kg', 'vap-liq phase change tendency from P3')
+ call addfld('MMF_VAP_ICE',(/'lev'/),'A', '/kg', 'vap-ice phase change tendency from P3')
endif
+ !----------------------------------------------------------------------------
+ ! CRM radiation columns
+ call addfld('CRM_RAD_T' ,dims_rad_3D, 'A', 'K', 'CRM rad column temperature' )
+ call addfld('CRM_RAD_QV' ,dims_rad_3D, 'A', 'kg/kg', 'CRM rad column water vapor' )
+ call addfld('CRM_RAD_QC' ,dims_rad_3D, 'A', 'kg/kg', 'CRM rad column cloud liquid' )
+ call addfld('CRM_RAD_QI' ,dims_rad_3D, 'A', 'kg/kg', 'CRM rad column cloud ice' )
+ call addfld('CRM_RAD_CLD',dims_rad_3D, 'A', 'kg/kg', 'CRM rad cloud fracton' )
+
!----------------------------------------------------------------------------
! ECPP output variables
#ifdef ECPP
@@ -178,6 +193,7 @@ subroutine crm_history_init(species_class)
call addfld('MMF_DQ', (/'lev'/), 'A','kg/kg/s','Q tendency due to CRM' )
call addfld('MMF_DQC', (/'lev'/), 'A','kg/kg/s','QC tendency due to CRM' )
call addfld('MMF_DQI', (/'lev'/), 'A','kg/kg/s','QI tendency due to CRM' )
+ call addfld('MMF_DQR', (/'lev'/), 'A','kg/kg/s','QR tendency due to CRM' )
call addfld('MMF_MC', (/'lev'/), 'A','kg/m2/s','Total mass flux from CRM' )
call addfld('MMF_MCUP', (/'lev'/), 'A','kg/m2/s','Updraft mass flux from CRM' )
call addfld('MMF_MCDN', (/'lev'/), 'A','kg/m2/s','Downdraft mass flux from CRM' )
@@ -188,6 +204,7 @@ subroutine crm_history_init(species_class)
call addfld('DU_CRM', (/'lev'/), 'A','/s', 'detrainment from updraft from CRM')
call addfld('EU_CRM', (/'lev'/), 'A','/s', 'entraiment rate from updraft')
call addfld('ED_CRM', (/'lev'/), 'A','/s', 'entraiment rate from downdraft')
+ call addfld('MMF_QV', (/'lev'/), 'A','kg/kg', 'Water vapor from CRM' )
call addfld('MMF_QC', (/'lev'/), 'A','kg/kg', 'Cloud water from CRM' )
call addfld('MMF_QI', (/'lev'/), 'A','kg/kg', 'Cloud ice from CRM' )
call addfld('MMF_QR', (/'lev'/), 'A','kg/kg', 'Rain from CRM' )
@@ -207,7 +224,36 @@ subroutine crm_history_init(species_class)
call addfld('MMF_QPFALL', (/'lev'/), 'A','kg/kg/s','Prec. water fall-out from CRM' )
call addfld('MMF_QPSRC', (/'lev'/), 'A','kg/kg/s','Prec. water source from CRM' )
call addfld('MMF_QTLS', (/'lev'/), 'A','kg/kg/s','L.S. Total Water Forcing in CRM' )
- call addfld('MMF_TLS', (/'lev'/), 'A','kg/kg/s','L.S. LIWSE Forcing in CRM' )
+ call addfld('MMF_TLS', (/'lev'/), 'A','K/s', 'L.S. Temperature Forcing in CRM' )
+ call addfld('MMF_RHODLS', (/'lev'/), 'A','kg/kg/s','L.S. dry density Forcing in CRM' )
+ call addfld('MMF_RHOVLS', (/'lev'/), 'A','kg/kg/s','L.S. vapor density Forcing in CRM' )
+ call addfld('MMF_RHOLLS', (/'lev'/), 'A','kg/kg/s','L.S. cloud liquid density Forcing in CRM' )
+ call addfld('MMF_RHOILS', (/'lev'/), 'A','kg/kg/s','L.S. cloud ice density Forcing in CRM' )
+
+ call addfld('MMF_DT_SGS ',(/'lev'/),'A','K/s' ,'CRM temperature tendency from SGS')
+ call addfld('MMF_DQV_SGS ',(/'lev'/),'A','kg/kg/s','CRM water vapor tendency from SGS')
+ call addfld('MMF_DQC_SGS ',(/'lev'/),'A','kg/kg/s','CRM cloud water tendency from SGS')
+ call addfld('MMF_DQI_SGS ',(/'lev'/),'A','kg/kg/s','CRM cloud ice tendency from SGS')
+ call addfld('MMF_DQR_SGS ',(/'lev'/),'A','kg/kg/s','CRM liquid rain tendency from SGS')
+
+ call addfld('MMF_DT_MICRO ',(/'lev'/),'A','K/s' ,'CRM temperature tendency from micro')
+ call addfld('MMF_DQV_MICRO',(/'lev'/),'A','kg/kg/s','CRM water vapor tendency from micro')
+ call addfld('MMF_DQC_MICRO',(/'lev'/),'A','kg/kg/s','CRM cloud water tendency from micro')
+ call addfld('MMF_DQI_MICRO',(/'lev'/),'A','kg/kg/s','CRM cloud ice tendency from micro')
+ call addfld('MMF_DQR_MICRO',(/'lev'/),'A','kg/kg/s','CRM liquid rain tendency from micro')
+
+ call addfld('MMF_DT_DYCOR ',(/'lev'/),'A','K/s' ,'CRM temperature tendency from dycor')
+ call addfld('MMF_DQV_DYCOR ',(/'lev'/),'A','kg/kg/s','CRM water vapor tendency from dycor')
+ call addfld('MMF_DQC_DYCOR ',(/'lev'/),'A','kg/kg/s','CRM cloud water tendency from dycor')
+ call addfld('MMF_DQI_DYCOR ',(/'lev'/),'A','kg/kg/s','CRM cloud ice tendency from dycor')
+ call addfld('MMF_DQR_DYCOR ',(/'lev'/),'A','kg/kg/s','CRM liquid rain tendency from dycor')
+
+ call addfld('MMF_DT_SPONGE ',(/'lev'/),'A','K/s' ,'CRM temperature tendency from sponge')
+ call addfld('MMF_DQV_SPONGE',(/'lev'/),'A','kg/kg/s','CRM water vapor tendency from sponge')
+ call addfld('MMF_DQC_SPONGE',(/'lev'/),'A','kg/kg/s','CRM cloud water tendency from sponge')
+ call addfld('MMF_DQI_SPONGE',(/'lev'/),'A','kg/kg/s','CRM cloud ice tendency from sponge')
+ call addfld('MMF_DQR_SPONGE',(/'lev'/),'A','kg/kg/s','CRM liquid rain tendency from sponge')
+
call addfld('MMF_TVFLUX', (/'lev'/), 'A', 'W/m2', 'Buoyancy Flux from CRM' )
call addfld('MMF_BUOY', (/'lev'/), 'A', 'W/m3', 'Buoyancy Term from CRM' )
call addfld('MMF_BUOYSD', (/'lev'/), 'A', 'W/m3', 'Std Dev of Buoyancy Term from CRM' )
@@ -358,7 +404,7 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
real(r8) :: MMF_DT_out(pcols,pver) ! CRM heating tendency
integer :: lchnk ! chunk identifier
integer :: ncol ! number of atmospheric columns
- integer :: ixcldliq, ixcldice ! liquid and ice constituent indices
+ integer :: ixcldliq,ixcldice,ixrain ! cloud liquid, ice, rain constituent indices
integer :: i, k, icol ! loop iterators
logical :: use_ECPP ! explicit cloud parameterized pollutants
logical :: use_MMF_VT ! CRM variance transport
@@ -377,6 +423,7 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
call cnst_get_ind('CLDLIQ', ixcldliq)
call cnst_get_ind('CLDICE', ixcldice)
+ call cnst_get_ind('RAINQM', ixrain)
call phys_getopts(use_ECPP_out = use_ECPP)
call phys_getopts(use_MMF_VT_out = use_MMF_VT)
@@ -391,6 +438,7 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
call outfld('MMF_DQ ',ptend%q(1:ncol,1,1), ncol, lchnk )
call outfld('MMF_DQC ',ptend%q(1:ncol,1,ixcldliq), ncol, lchnk )
call outfld('MMF_DQI ',ptend%q(1:ncol,1,ixcldice), ncol, lchnk )
+ call outfld('MMF_DQR ',ptend%q(1:ncol,1,ixrain ), ncol, lchnk )
! CRM radiative heating rate
! NOTE: We output the radiative heating rates (MMF_QRS and MMF_QRL) here
@@ -414,7 +462,7 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
call outfld('CRM_V ',crm_state%v_wind (icol_beg:icol_end,:,:,:), ncol, lchnk )
call outfld('CRM_W ',crm_state%w_wind (icol_beg:icol_end,:,:,:), ncol, lchnk )
call outfld('CRM_T ',crm_state%temperature(icol_beg:icol_end,:,:,:), ncol, lchnk )
- call outfld('CRM_QV ',crm_state%qv(icol_beg:icol_end,:,:,:) , ncol, lchnk )
+ call outfld('CRM_QV ',crm_state%qv (icol_beg:icol_end,:,:,:), ncol, lchnk )
!----------------------------------------------------------------------------
! Turbulence parameter on CRM grid
@@ -423,22 +471,33 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
!----------------------------------------------------------------------------
! CRM condensate and precipitation on CRM grid
- call outfld('CRM_QC ',crm_output%qcl (icol_beg:icol_end,:,:,:),ncol, lchnk )
- call outfld('CRM_QI ',crm_output%qci (icol_beg:icol_end,:,:,:),ncol, lchnk )
call outfld('CRM_PREC',crm_output%prec_crm(icol_beg:icol_end,:,:), ncol, lchnk )
if (MMF_microphysics_scheme .eq. 'sam1mom') then
- call outfld('CRM_QPC ',crm_output%qpl (icol_beg:icol_end,:,:,:),ncol, lchnk )
- call outfld('CRM_QPI ',crm_output%qpi (icol_beg:icol_end,:,:,:),ncol, lchnk )
+ call outfld('CRM_QC ',crm_output%qcl(icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_QI ',crm_output%qci(icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_QPC',crm_output%qpl(icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_QPI',crm_output%qpi(icol_beg:icol_end,:,:,:), ncol, lchnk )
end if
if (MMF_microphysics_scheme .eq. 'p3') then
+ call outfld('CRM_QC ',crm_state%qc(icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_QI ',crm_state%qi(icol_beg:icol_end,:,:,:), ncol, lchnk )
call outfld('CRM_NC ',crm_state%nc(icol_beg:icol_end,:,:,:), ncol, lchnk )
call outfld('CRM_NI ',crm_state%ni(icol_beg:icol_end,:,:,:), ncol, lchnk )
call outfld('CRM_NR ',crm_state%nr(icol_beg:icol_end,:,:,:), ncol, lchnk )
call outfld('CRM_QR ',crm_state%qr(icol_beg:icol_end,:,:,:), ncol, lchnk )
endif
+ !----------------------------------------------------------------------------
+ ! CRM radiation columns
+ call outfld('CRM_RAD_T' ,crm_rad%temperature(icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_RAD_QV' ,crm_rad%qv (icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_RAD_QC' ,crm_rad%qc (icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_RAD_QI' ,crm_rad%qi (icol_beg:icol_end,:,:,:), ncol, lchnk )
+ call outfld('CRM_RAD_CLD',crm_rad%cld (icol_beg:icol_end,:,:,:), ncol, lchnk )
+
!----------------------------------------------------------------------------
! CRM domain average condensate and precipitation
+ call outfld('MMF_QV ',crm_output%qv_mean(icol_beg:icol_end,:), ncol ,lchnk )
call outfld('MMF_QC ',crm_output%qc_mean(icol_beg:icol_end,:), ncol ,lchnk )
call outfld('MMF_QI ',crm_output%qi_mean(icol_beg:icol_end,:), ncol ,lchnk )
call outfld('MMF_QR ',crm_output%qr_mean(icol_beg:icol_end,:), ncol ,lchnk )
@@ -446,6 +505,8 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
call outfld('MMF_NC ',crm_output%nc_mean(icol_beg:icol_end,:), ncol, lchnk )
call outfld('MMF_NI ',crm_output%ni_mean(icol_beg:icol_end,:), ncol, lchnk )
call outfld('MMF_NR ',crm_output%nr_mean(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_RHOD ',crm_output%rho_d_mean(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_RHOV ',crm_output%rho_v_mean(icol_beg:icol_end,:), ncol, lchnk )
endif
!----------------------------------------------------------------------------
@@ -467,6 +528,40 @@ subroutine crm_history_out(state, ptend, crm_state, crm_rad, crm_output, &
call outfld('MMF_QPFALL',crm_output%qp_fall (icol_beg:icol_end,:), ncol, lchnk )
call outfld('MMF_QPSRC ',crm_output%qp_src (icol_beg:icol_end,:), ncol, lchnk )
call outfld('MMF_TLS ',crm_output%t_ls (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_RHODLS',crm_output%rho_d_ls (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_RHOVLS',crm_output%rho_v_ls (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_RHOLLS',crm_output%rho_l_ls (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_RHOILS',crm_output%rho_i_ls (icol_beg:icol_end,:), ncol, lchnk )
+
+ call outfld('MMF_DT_SGS ',crm_output%dt_sgs (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQV_SGS ',crm_output%dqv_sgs (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQC_SGS ',crm_output%dqc_sgs (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQI_SGS ',crm_output%dqi_sgs (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQR_SGS ',crm_output%dqr_sgs (icol_beg:icol_end,:), ncol, lchnk )
+
+ call outfld('MMF_DT_MICRO ',crm_output%dt_micro (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQV_MICRO',crm_output%dqv_micro(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQC_MICRO',crm_output%dqc_micro(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQI_MICRO',crm_output%dqi_micro(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQR_MICRO',crm_output%dqr_micro(icol_beg:icol_end,:), ncol, lchnk )
+
+ call outfld('MMF_DT_DYCOR', crm_output%dt_dycor (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQV_DYCOR', crm_output%dqv_dycor (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQC_DYCOR', crm_output%dqc_dycor (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQI_DYCOR', crm_output%dqi_dycor (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQR_DYCOR', crm_output%dqr_dycor (icol_beg:icol_end,:), ncol, lchnk )
+
+ call outfld('MMF_DT_SPONGE', crm_output%dt_sponge (icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQV_SPONGE',crm_output%dqv_sponge(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQC_SPONGE',crm_output%dqc_sponge(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQI_SPONGE',crm_output%dqi_sponge(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_DQR_SPONGE',crm_output%dqr_sponge(icol_beg:icol_end,:), ncol, lchnk )
+
+ if (MMF_microphysics_scheme .eq. 'p3') then
+ call outfld('MMF_LIQ_ICE',crm_output%liq_ice_exchange(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_VAP_LIQ',crm_output%vap_liq_exchange(icol_beg:icol_end,:), ncol, lchnk )
+ call outfld('MMF_VAP_ICE',crm_output%vap_ice_exchange(icol_beg:icol_end,:), ncol, lchnk )
+ endif
! NOTE: these should overwrite cloud outputs from non-MMF routines
call outfld('CLOUD ',crm_output%cld (icol_beg:icol_end,:), ncol, lchnk )
diff --git a/components/eam/src/physics/crm/crm_input_module.F90 b/components/eam/src/physics/crm/crm_input_module.F90
index 122be63d6d17..ddb75156e1ba 100644
--- a/components/eam/src/physics/crm/crm_input_module.F90
+++ b/components/eam/src/physics/crm/crm_input_module.F90
@@ -33,10 +33,10 @@ module crm_input_module
real(crm_rknd), allocatable :: tau00 (:) ! large-scale surface stress (N/m2)
real(crm_rknd), allocatable :: wndls (:) ! large-scale surface wind (m/s)
real(crm_rknd), allocatable :: bflxls (:) ! large-scale surface buoyancy flux (K m/s)
- real(crm_rknd), allocatable :: fluxu00(:) ! surface momenent fluxes [N/m2]
- real(crm_rknd), allocatable :: fluxv00(:) ! surface momenent fluxes [N/m2]
- real(crm_rknd), allocatable :: fluxt00(:) ! surface sensible heat fluxes [K Kg/ (m2 s)]
- real(crm_rknd), allocatable :: fluxq00(:) ! surface latent heat fluxes [ kg/(m2 s)]
+ real(crm_rknd), allocatable :: fluxu00(:) ! surface momenent fluxes [N/m2]
+ real(crm_rknd), allocatable :: fluxv00(:) ! surface momenent fluxes [N/m2]
+ real(crm_rknd), allocatable :: fluxt00(:) ! surface sensible heat fluxes [W/m2]
+ real(crm_rknd), allocatable :: fluxq00(:) ! surface latent heat fluxes [W/m2]
real(crm_rknd), allocatable :: ul_esmt(:,:) ! input u for ESMT
real(crm_rknd), allocatable :: vl_esmt(:,:) ! input v for ESMT
@@ -46,9 +46,9 @@ module crm_input_module
real(crm_rknd), allocatable :: u_vt(:,:) ! CRM input of variance used for forcing tendency
! inputs for P3
- real(crm_rknd), allocatable :: nccn(:,:) ! CCN number concentration [kg-1]
- real(crm_rknd), allocatable :: nc_nuceat_tend(:,:) ! activated CCN number tendency [kg-1 s-1]
- real(crm_rknd), allocatable :: ni_activated(:,:) ! activated ice nuclei concentration [kg-1]
+ real(crm_rknd), allocatable :: nccn_prescribed(:,:)! CCN number concentration [#/kg]
+ real(crm_rknd), allocatable :: nc_nuceat_tend(:,:) ! activated CCN number tendency [#/kg/s]
+ real(crm_rknd), allocatable :: ni_activated(:,:) ! activated ice nuclei concentration [#/kg]
end type crm_input_type
!------------------------------------------------------------------------------------------------
@@ -116,10 +116,10 @@ subroutine crm_input_initialize(input, ncrms, nlev, MMF_microphysics_scheme)
call prefetch(input%u_vt)
if (trim(MMF_microphysics_scheme).eq.'p3') then
- if (.not. allocated(input%nccn )) allocate(input%nccn(ncrms,nlev))
- if (.not. allocated(input%nc_nuceat_tend)) allocate(input%nc_nuceat_tend(ncrms,nlev))
- if (.not. allocated(input%ni_activated )) allocate(input%ni_activated(ncrms,nlev))
- call prefetch(input%nccn)
+ if (.not. allocated(input%nccn_prescribed)) allocate(input%nccn_prescribed(ncrms,nlev))
+ if (.not. allocated(input%nc_nuceat_tend )) allocate(input%nc_nuceat_tend(ncrms,nlev))
+ if (.not. allocated(input%ni_activated )) allocate(input%ni_activated(ncrms,nlev))
+ call prefetch(input%nccn_prescribed)
call prefetch(input%nc_nuceat_tend)
call prefetch(input%ni_activated)
end if
@@ -156,9 +156,9 @@ subroutine crm_input_initialize(input, ncrms, nlev, MMF_microphysics_scheme)
input%u_vt = 0
if (trim(MMF_microphysics_scheme).eq.'p3') then
- input%nccn = 0
- input%nc_nuceat_tend = 0
- input%ni_activated = 0
+ input%nccn_prescribed = 0
+ input%nc_nuceat_tend = 0
+ input%ni_activated = 0
end if
end subroutine crm_input_initialize
@@ -197,9 +197,9 @@ subroutine crm_input_finalize(input, MMF_microphysics_scheme)
if (allocated(input%q_vt)) deallocate(input%q_vt)
if (allocated(input%u_vt)) deallocate(input%u_vt)
- if (allocated(input%nccn )) deallocate(input%nccn)
- if (allocated(input%nc_nuceat_tend)) deallocate(input%nc_nuceat_tend)
- if (allocated(input%ni_activated )) deallocate(input%ni_activated)
+ if (allocated(input%nccn_prescribed)) deallocate(input%nccn_prescribed)
+ if (allocated(input%nc_nuceat_tend )) deallocate(input%nc_nuceat_tend)
+ if (allocated(input%ni_activated )) deallocate(input%ni_activated)
end subroutine crm_input_finalize
diff --git a/components/eam/src/physics/crm/crm_output_module.F90 b/components/eam/src/physics/crm/crm_output_module.F90
index 9cc39940da10..07af5cd0eec8 100644
--- a/components/eam/src/physics/crm/crm_output_module.F90
+++ b/components/eam/src/physics/crm/crm_output_module.F90
@@ -49,11 +49,16 @@ module crm_output_module
! do not need to be calculated here, and might make more sense to calculate at the
! crm_physics_tend level. For example, I think tendencies should be calculated in
! crm_physics_tend, from, for example, something like crm_output%uwind - crm_input%uwind.
+ real(crm_rknd), allocatable :: qv_mean(:,:) ! mean cloud water
real(crm_rknd), allocatable :: qc_mean(:,:) ! mean cloud water
real(crm_rknd), allocatable :: qi_mean(:,:) ! mean cloud ice
real(crm_rknd), allocatable :: qr_mean(:,:) ! mean rain
real(crm_rknd), allocatable :: qs_mean(:,:) ! mean snow
real(crm_rknd), allocatable :: qg_mean(:,:) ! mean graupel
+ real(crm_rknd), allocatable :: qm_mean(:,:) ! mean ice rime mass
+ real(crm_rknd), allocatable :: bm_mean(:,:) ! mean ice rime volume
+ real(crm_rknd), allocatable :: rho_d_mean(:,:) ! mean dry density
+ real(crm_rknd), allocatable :: rho_v_mean(:,:) ! mean vapor density
real(crm_rknd), allocatable :: nc_mean(:,:) ! mean cloud water (#/kg)
real(crm_rknd), allocatable :: ni_mean(:,:) ! mean cloud ice (#/kg)
@@ -77,6 +82,11 @@ module crm_output_module
real(crm_rknd), allocatable :: cld (:,:) ! cloud fraction
real(crm_rknd), allocatable :: gicewp(:,:) ! ice water path
real(crm_rknd), allocatable :: gliqwp(:,:) ! ice water path
+
+ real(crm_rknd), allocatable :: liq_ice_exchange(:,:) ! P3 liq-ice phase change tendency
+ real(crm_rknd), allocatable :: vap_liq_exchange(:,:) ! P3 vap-liq phase change tendency
+ real(crm_rknd), allocatable :: vap_ice_exchange(:,:) ! P3 vap-ice phase change tendency
+
real(crm_rknd), allocatable :: mctot (:,:) ! cloud mass flux
real(crm_rknd), allocatable :: mcup (:,:) ! updraft cloud mass flux
real(crm_rknd), allocatable :: mcdn (:,:) ! downdraft cloud mass flux
@@ -112,11 +122,40 @@ module crm_output_module
real(crm_rknd), allocatable :: t_ls (:,:) ! tend of lwse due to large-scale [kg/kg/s] ???
real(crm_rknd), allocatable :: prectend (:) ! column integrated tend in precip water+ice [kg/m2/s]
real(crm_rknd), allocatable :: precstend (:) ! column integrated tend in precip ice [kg/m2/s]
- real(crm_rknd), allocatable :: taux (:) ! zonal CRM surface stress perturbation [N/m2]
- real(crm_rknd), allocatable :: tauy (:) ! merid CRM surface stress perturbation [N/m2]
+ real(crm_rknd), allocatable :: taux (:) ! zonal CRM surface stress perturbation [N/m2]
+ real(crm_rknd), allocatable :: tauy (:) ! merid CRM surface stress perturbation [N/m2]
real(crm_rknd), allocatable :: z0m (:) ! surface stress [N/m2]
real(crm_rknd), allocatable :: subcycle_factor(:) ! crm cpu efficiency
+ real(crm_rknd), allocatable :: dt_sgs (:,:) ! CRM temperature tendency from SGS [K/s]
+ real(crm_rknd), allocatable :: dqv_sgs (:,:) ! CRM water vapor tendency from SGS [kg/kg/s]
+ real(crm_rknd), allocatable :: dqc_sgs (:,:) ! CRM cloud water tendency from SGS [kg/kg/s]
+ real(crm_rknd), allocatable :: dqi_sgs (:,:) ! CRM cloud ice tendency from SGS [kg/kg/s]
+ real(crm_rknd), allocatable :: dqr_sgs (:,:) ! CRM liquid rain tendency from SGS [kg/kg/s]
+
+ real(crm_rknd), allocatable :: dt_micro (:,:) ! CRM temperature tendency from micro [K/s]
+ real(crm_rknd), allocatable :: dqv_micro (:,:) ! CRM water vapor tendency from micro [kg/kg/s]
+ real(crm_rknd), allocatable :: dqc_micro (:,:) ! CRM cloud water tendency from micro [kg/kg/s]
+ real(crm_rknd), allocatable :: dqi_micro (:,:) ! CRM cloud ice tendency from micro [kg/kg/s]
+ real(crm_rknd), allocatable :: dqr_micro (:,:) ! CRM liquid rain tendency from micro [kg/kg/s]
+
+ real(crm_rknd), allocatable :: dt_dycor (:,:)
+ real(crm_rknd), allocatable :: dqv_dycor (:,:)
+ real(crm_rknd), allocatable :: dqc_dycor (:,:)
+ real(crm_rknd), allocatable :: dqi_dycor (:,:)
+ real(crm_rknd), allocatable :: dqr_dycor (:,:)
+
+ real(crm_rknd), allocatable :: dt_sponge (:,:)
+ real(crm_rknd), allocatable :: dqv_sponge(:,:)
+ real(crm_rknd), allocatable :: dqc_sponge(:,:)
+ real(crm_rknd), allocatable :: dqi_sponge(:,:)
+ real(crm_rknd), allocatable :: dqr_sponge(:,:)
+
+ real(crm_rknd), allocatable :: rho_d_ls (:,:) ! large-scale forcing of dry density [kg/m3/s]
+ real(crm_rknd), allocatable :: rho_v_ls (:,:) ! large-scale forcing of vapor density [kg/m3/s]
+ real(crm_rknd), allocatable :: rho_l_ls (:,:) ! large-scale forcing of vapor density [kg/m3/s]
+ real(crm_rknd), allocatable :: rho_i_ls (:,:) ! large-scale forcing of vapor density [kg/m3/s]
+
end type crm_output_type
contains
@@ -162,11 +201,16 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
! NOTE: this output had a bug in the previous implementation
if (.not. allocated(output%cldtop)) allocate(output%cldtop(ncol,nlev))
+ if (.not. allocated(output%qv_mean)) allocate(output%qv_mean(ncol,nlev))
if (.not. allocated(output%qc_mean)) allocate(output%qc_mean(ncol,nlev))
if (.not. allocated(output%qi_mean)) allocate(output%qi_mean(ncol,nlev))
if (.not. allocated(output%qr_mean)) allocate(output%qr_mean(ncol,nlev))
if (.not. allocated(output%qs_mean)) allocate(output%qs_mean(ncol,nlev))
if (.not. allocated(output%qg_mean)) allocate(output%qg_mean(ncol,nlev))
+ if (.not. allocated(output%qm_mean)) allocate(output%qm_mean(ncol,nlev))
+ if (.not. allocated(output%bm_mean)) allocate(output%bm_mean(ncol,nlev))
+ if (.not. allocated(output%rho_d_mean)) allocate(output%rho_d_mean(ncol,nlev))
+ if (.not. allocated(output%rho_v_mean)) allocate(output%rho_v_mean(ncol,nlev))
call prefetch(output%qcl)
call prefetch(output%qci)
@@ -193,11 +237,16 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
call prefetch(output%precsc)
call prefetch(output%precsl)
call prefetch(output%cldtop)
+ call prefetch(output%qv_mean)
call prefetch(output%qc_mean)
call prefetch(output%qi_mean)
call prefetch(output%qr_mean)
call prefetch(output%qs_mean)
call prefetch(output%qg_mean)
+ call prefetch(output%qm_mean)
+ call prefetch(output%bm_mean)
+ call prefetch(output%rho_d_mean)
+ call prefetch(output%rho_v_mean)
if (.not. allocated(output%nc_mean)) allocate(output%nc_mean(ncol,nlev))
if (.not. allocated(output%ni_mean)) allocate(output%ni_mean(ncol,nlev))
@@ -220,6 +269,11 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
if (.not. allocated(output%cld )) allocate(output%cld (ncol,nlev)) ! cloud fraction
if (.not. allocated(output%gicewp)) allocate(output%gicewp(ncol,nlev)) ! ice water path
if (.not. allocated(output%gliqwp)) allocate(output%gliqwp(ncol,nlev)) ! ice water path
+
+ if (.not. allocated(output%liq_ice_exchange)) allocate(output%liq_ice_exchange(ncol,nlev)) ! P3 liq-ice phase change tendency
+ if (.not. allocated(output%vap_liq_exchange)) allocate(output%vap_liq_exchange(ncol,nlev)) ! P3 vap-liq phase change tendency
+ if (.not. allocated(output%vap_ice_exchange)) allocate(output%vap_ice_exchange(ncol,nlev)) ! P3 vap-ice phase change tendency
+
if (.not. allocated(output%mctot )) allocate(output%mctot (ncol,nlev)) ! cloud mass flux
if (.not. allocated(output%mcup )) allocate(output%mcup (ncol,nlev)) ! updraft cloud mass flux
if (.not. allocated(output%mcdn )) allocate(output%mcdn (ncol,nlev)) ! downdraft cloud mass flux
@@ -258,6 +312,35 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
if (.not. allocated(output%z0m )) allocate(output%z0m (ncol))
if (.not. allocated(output%subcycle_factor)) allocate(output%subcycle_factor(ncol))
+ if (.not. allocated(output%dt_sgs )) allocate(output%dt_sgs (ncol,nlev))
+ if (.not. allocated(output%dqv_sgs )) allocate(output%dqv_sgs (ncol,nlev))
+ if (.not. allocated(output%dqc_sgs )) allocate(output%dqc_sgs (ncol,nlev))
+ if (.not. allocated(output%dqi_sgs )) allocate(output%dqi_sgs (ncol,nlev))
+ if (.not. allocated(output%dqr_sgs )) allocate(output%dqr_sgs (ncol,nlev))
+
+ if (.not. allocated(output%dt_micro )) allocate(output%dt_micro (ncol,nlev))
+ if (.not. allocated(output%dqv_micro )) allocate(output%dqv_micro (ncol,nlev))
+ if (.not. allocated(output%dqc_micro )) allocate(output%dqc_micro (ncol,nlev))
+ if (.not. allocated(output%dqi_micro )) allocate(output%dqi_micro (ncol,nlev))
+ if (.not. allocated(output%dqr_micro )) allocate(output%dqr_micro (ncol,nlev))
+
+ if (.not. allocated(output%dt_dycor )) allocate(output%dt_dycor (ncol,nlev))
+ if (.not. allocated(output%dqv_dycor )) allocate(output%dqv_dycor (ncol,nlev))
+ if (.not. allocated(output%dqc_dycor )) allocate(output%dqc_dycor (ncol,nlev))
+ if (.not. allocated(output%dqi_dycor )) allocate(output%dqi_dycor (ncol,nlev))
+ if (.not. allocated(output%dqr_dycor )) allocate(output%dqr_dycor (ncol,nlev))
+
+ if (.not. allocated(output%dt_sponge )) allocate(output%dt_sponge (ncol,nlev))
+ if (.not. allocated(output%dqv_sponge )) allocate(output%dqv_sponge (ncol,nlev))
+ if (.not. allocated(output%dqc_sponge )) allocate(output%dqc_sponge (ncol,nlev))
+ if (.not. allocated(output%dqi_sponge )) allocate(output%dqi_sponge (ncol,nlev))
+ if (.not. allocated(output%dqr_sponge )) allocate(output%dqr_sponge (ncol,nlev))
+
+ if (.not. allocated(output%rho_d_ls )) allocate(output%rho_d_ls (ncol,nlev))
+ if (.not. allocated(output%rho_v_ls )) allocate(output%rho_v_ls (ncol,nlev))
+ if (.not. allocated(output%rho_l_ls )) allocate(output%rho_l_ls (ncol,nlev))
+ if (.not. allocated(output%rho_i_ls )) allocate(output%rho_i_ls (ncol,nlev))
+
call prefetch(output%sltend )
call prefetch(output%qltend )
call prefetch(output%qcltend )
@@ -273,11 +356,17 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
call prefetch(output%cld )
call prefetch(output%gicewp )
call prefetch(output%gliqwp )
+
+ call prefetch(output%liq_ice_exchange)
+ call prefetch(output%vap_liq_exchange)
+ call prefetch(output%vap_ice_exchange)
+
call prefetch(output%mctot )
call prefetch(output%mcup )
call prefetch(output%mcdn )
call prefetch(output%mcuup )
call prefetch(output%mcudn )
+
call prefetch(output%mu_crm )
call prefetch(output%md_crm )
call prefetch(output%du_crm )
@@ -309,6 +398,29 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
call prefetch(output%z0m )
call prefetch(output%subcycle_factor )
+ call prefetch(output%dt_sgs)
+ call prefetch(output%dqv_sgs)
+ call prefetch(output%dqc_sgs)
+ call prefetch(output%dqi_sgs)
+ call prefetch(output%dt_micro)
+ call prefetch(output%dqv_micro)
+ call prefetch(output%dqc_micro)
+ call prefetch(output%dqi_micro)
+
+ call prefetch(output%dt_dycor )
+ call prefetch(output%dqv_dycor )
+ call prefetch(output%dqc_dycor )
+ call prefetch(output%dqi_dycor )
+ call prefetch(output%dt_sponge )
+ call prefetch(output%dqv_sponge)
+ call prefetch(output%dqc_sponge)
+ call prefetch(output%dqi_sponge)
+
+ call prefetch(output%rho_d_ls)
+ call prefetch(output%rho_v_ls)
+ call prefetch(output%rho_l_ls)
+ call prefetch(output%rho_i_ls)
+
! Initialize
output%qcl = 0
output%qci = 0
@@ -341,11 +453,16 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
output%precsc = 0
output%precsl = 0
+ output%qv_mean = 0
output%qc_mean = 0
output%qi_mean = 0
output%qr_mean = 0
output%qs_mean = 0
output%qg_mean = 0
+ output%qm_mean = 0
+ output%bm_mean = 0
+ output%rho_d_mean = 0
+ output%rho_v_mean = 0
output%nc_mean = 0
output%ni_mean = 0
@@ -368,6 +485,11 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
output%cld = 0
output%gicewp = 0
output%gliqwp = 0
+
+ output%liq_ice_exchange = 0
+ output%vap_liq_exchange = 0
+ output%vap_ice_exchange = 0
+
output%mctot = 0
output%mcup = 0
output%mcdn = 0
@@ -403,11 +525,34 @@ subroutine crm_output_initialize(output, ncol, nlev, crm_nx, crm_ny, crm_nz, MMF
output%t_ls = 0
output%prectend = 0
output%precstend = 0
- output%taux = 0
- output%tauy = 0
+ output%taux = 0
+ output%tauy = 0
output%z0m = 0
output%subcycle_factor = 0
+ output%dt_sgs = 0
+ output%dqv_sgs = 0
+ output%dqc_sgs = 0
+ output%dqi_sgs = 0
+ output%dt_micro = 0
+ output%dqv_micro = 0
+ output%dqc_micro = 0
+ output%dqi_micro = 0
+
+ output%dt_dycor = 0
+ output%dqv_dycor = 0
+ output%dqc_dycor = 0
+ output%dqi_dycor = 0
+ output%dt_sponge = 0
+ output%dqv_sponge = 0
+ output%dqc_sponge = 0
+ output%dqi_sponge = 0
+
+ output%rho_d_ls = 0
+ output%rho_v_ls = 0
+ output%rho_l_ls = 0
+ output%rho_i_ls = 0
+
end subroutine crm_output_initialize
!------------------------------------------------------------------------------------------------
subroutine crm_output_finalize(output, MMF_microphysics_scheme)
@@ -442,11 +587,16 @@ subroutine crm_output_finalize(output, MMF_microphysics_scheme)
if (allocated(output%precsc)) deallocate(output%precsc)
if (allocated(output%precsl)) deallocate(output%precsl)
+ if (allocated(output%qv_mean)) deallocate(output%qv_mean)
if (allocated(output%qc_mean)) deallocate(output%qc_mean)
if (allocated(output%qi_mean)) deallocate(output%qi_mean)
if (allocated(output%qr_mean)) deallocate(output%qr_mean)
if (allocated(output%qs_mean)) deallocate(output%qs_mean)
if (allocated(output%qg_mean)) deallocate(output%qg_mean)
+ if (allocated(output%qm_mean)) deallocate(output%qm_mean)
+ if (allocated(output%bm_mean)) deallocate(output%bm_mean)
+ if (allocated(output%rho_d_mean)) deallocate(output%rho_d_mean)
+ if (allocated(output%rho_v_mean)) deallocate(output%rho_v_mean)
if (allocated(output%nc_mean)) deallocate(output%nc_mean)
if (allocated(output%ni_mean)) deallocate(output%ni_mean)
@@ -469,6 +619,11 @@ subroutine crm_output_finalize(output, MMF_microphysics_scheme)
if (allocated(output%cld)) deallocate(output%cld)
if (allocated(output%gicewp)) deallocate(output%gicewp)
if (allocated(output%gliqwp)) deallocate(output%gliqwp)
+
+ if (allocated(output%liq_ice_exchange)) deallocate(output%liq_ice_exchange)
+ if (allocated(output%vap_liq_exchange)) deallocate(output%vap_liq_exchange)
+ if (allocated(output%vap_ice_exchange)) deallocate(output%vap_ice_exchange)
+
if (allocated(output%mctot)) deallocate(output%mctot)
if (allocated(output%mcup)) deallocate(output%mcup)
if (allocated(output%mcdn)) deallocate(output%mcdn)
@@ -507,6 +662,29 @@ subroutine crm_output_finalize(output, MMF_microphysics_scheme)
if (allocated(output%z0m)) deallocate(output%z0m)
if (allocated(output%subcycle_factor)) deallocate(output%subcycle_factor)
+ if (allocated(output%dt_sgs )) deallocate(output%dt_sgs)
+ if (allocated(output%dqv_sgs )) deallocate(output%dqv_sgs)
+ if (allocated(output%dqc_sgs )) deallocate(output%dqc_sgs)
+ if (allocated(output%dqi_sgs )) deallocate(output%dqi_sgs)
+ if (allocated(output%dt_micro )) deallocate(output%dt_micro)
+ if (allocated(output%dqv_micro)) deallocate(output%dqv_micro)
+ if (allocated(output%dqc_micro)) deallocate(output%dqc_micro)
+ if (allocated(output%dqi_micro)) deallocate(output%dqi_micro)
+
+ if (allocated(output%dt_dycor )) deallocate(output%dt_dycor )
+ if (allocated(output%dqv_dycor )) deallocate(output%dqv_dycor )
+ if (allocated(output%dqc_dycor )) deallocate(output%dqc_dycor )
+ if (allocated(output%dqi_dycor )) deallocate(output%dqi_dycor )
+ if (allocated(output%dt_sponge )) deallocate(output%dt_sponge )
+ if (allocated(output%dqv_sponge)) deallocate(output%dqv_sponge)
+ if (allocated(output%dqc_sponge)) deallocate(output%dqc_sponge)
+ if (allocated(output%dqi_sponge)) deallocate(output%dqi_sponge)
+
+ if (allocated(output%rho_d_ls)) deallocate(output%rho_d_ls)
+ if (allocated(output%rho_v_ls)) deallocate(output%rho_v_ls)
+ if (allocated(output%rho_l_ls)) deallocate(output%rho_l_ls)
+ if (allocated(output%rho_i_ls)) deallocate(output%rho_i_ls)
+
end subroutine crm_output_finalize
!------------------------------------------------------------------------------------------------
diff --git a/components/eam/src/physics/crm/crm_physics.F90 b/components/eam/src/physics/crm/crm_physics.F90
index d62243fd1954..66e5b57e28eb 100644
--- a/components/eam/src/physics/crm/crm_physics.F90
+++ b/components/eam/src/physics/crm/crm_physics.F90
@@ -4,7 +4,10 @@ module crm_physics
!---------------------------------------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use shr_sys_mod, only: shr_sys_flush
+ use shr_const_mod, only: SHR_CONST_RDAIR,SHR_CONST_RWV
+ use spmd_utils, only: masterproc
use cam_abortutils, only: endrun
+ use cam_control_mod, only: nsrest ! restart flag
use cam_logfile, only: iulog
use physics_types, only: physics_state, physics_tend
use ppgrid, only: begchunk, endchunk, pcols, pver, pverp
@@ -65,6 +68,12 @@ module crm_physics
integer :: crm_t_prev_idx = -1
integer :: crm_q_prev_idx = -1
+ integer :: crm_shoc_tk_idx = -1
+ integer :: crm_shoc_tkh_idx = -1
+ integer :: crm_shoc_wthv_idx = -1
+ integer :: crm_shoc_relvar_idx = -1
+ integer :: crm_shoc_cldfrac_idx = -1
+
contains
!===================================================================================================
!===================================================================================================
@@ -73,7 +82,6 @@ subroutine crm_physics_register()
!---------------------------------------------------------------------------------------------------
! Purpose: add necessary fields into physics buffer
!---------------------------------------------------------------------------------------------------
- use spmd_utils, only: masterproc
use physconst, only: cpair, mwh2o, mwdry
use physics_buffer, only: dyn_time_lvls, pbuf_add_field, dtype_r8, pbuf_get_index
use phys_control, only: phys_getopts, use_gw_convect
@@ -83,9 +91,11 @@ subroutine crm_physics_register()
crm_nx_rad, crm_ny_rad
use constituents, only: cnst_add
use crm_history, only: crm_history_register
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
+ use gator_mod, only: gator_init
+#endif
#if defined(MMF_SAMXX)
use cpp_interface_mod, only: setparm
- use gator_mod, only: gator_init
#elif defined(MMF_SAM) || defined(MMF_SAMOMP)
use setparm_mod , only: setparm
#endif
@@ -106,7 +116,7 @@ subroutine crm_physics_register()
integer, dimension(4) :: dims_crm_rad
integer :: cnst_ind ! dummy for adding new constituents for variance transport
!----------------------------------------------------------------------------
-#if defined(MMF_SAMXX)
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
call gator_init()
#endif
@@ -140,7 +150,9 @@ subroutine crm_physics_register()
!----------------------------------------------------------------------------
! Setup CRM internal parameters
!----------------------------------------------------------------------------
+#if defined(MMF_SAMXX) || defined(MMF_SAM) || defined(MMF_SAMOMP)
call setparm()
+#endif
if (use_MMF_VT) then
! add variance tracers
@@ -194,6 +206,7 @@ subroutine crm_physics_register()
call pbuf_add_field('CRM_V', 'global',dtype_r8,dims_crm_3D,idx)
call pbuf_add_field('CRM_W', 'global',dtype_r8,dims_crm_3D,idx)
call pbuf_add_field('CRM_T', 'global',dtype_r8,dims_crm_3D,idx)
+ call pbuf_add_field('CRM_RHO', 'global',dtype_r8,dims_crm_3D,idx)
! Radiation
call pbuf_add_field('CRM_T_RAD', 'physpkg',dtype_r8,dims_crm_rad,crm_t_rad_idx)
@@ -238,6 +251,11 @@ subroutine crm_physics_register()
call pbuf_add_field('CRM_T_PREV','global', dtype_r8,dims_crm_3D,crm_t_prev_idx)
call pbuf_add_field('CRM_Q_PREV','global', dtype_r8,dims_crm_3D,crm_q_prev_idx)
if (prog_modal_aero) call pbuf_add_field('RATE1_CW2PR_ST','physpkg',dtype_r8,dims_gcm_2D,idx)
+ call pbuf_add_field('CRM_SHOC_TK' ,'global', dtype_r8,dims_crm_3D,crm_shoc_tk_idx)
+ call pbuf_add_field('CRM_SHOC_THH' ,'global', dtype_r8,dims_crm_3D,crm_shoc_tkh_idx)
+ call pbuf_add_field('CRM_SHOC_WTHV' ,'global', dtype_r8,dims_crm_3D,crm_shoc_wthv_idx)
+ call pbuf_add_field('CRM_SHOC_RELVAR' ,'global', dtype_r8,dims_crm_3D,crm_shoc_relvar_idx)
+ call pbuf_add_field('CRM_SHOC_CLDFRAC','global', dtype_r8,dims_crm_3D,crm_shoc_cldfrac_idx)
end if
! CRM rad stuff specific to COSP; this does not strictly need to be in
@@ -323,12 +341,12 @@ subroutine crm_physics_init(state, pbuf2d, species_class)
use phys_control, only: phys_getopts, use_gw_convect
use phys_grid, only: get_ncols_p
use crm_history, only: crm_history_init
- use cam_history, only: addfld
- use constituents, only: apcnst, bpcnst, cnst_name, cnst_longname, cnst_get_ind
+ use cam_history, only: addfld, hist_fld_active
+ use constituents, only: apcnst, bpcnst, cnst_name, cnst_longname, cnst_get_ind
+ use constituents, only: pcnst, cnst_get_ind
#ifdef ECPP
use module_ecpp_ppdriver2, only: papampollu_init
#endif
- use constituents, only: pcnst, cnst_get_ind
!----------------------------------------------------------------------------
! interface variables
! NOTE - species_class is an input so it needs to be outside of ifdef MODAL_AERO for 1-mom micro
@@ -339,11 +357,11 @@ subroutine crm_physics_init(state, pbuf2d, species_class)
! local variables
integer :: m, mm, c
integer :: ncnst
- integer :: ierror ! Error code
logical :: use_ECPP
logical :: use_MMF_VT
character(len=16) :: MMF_microphysics_scheme
- integer :: lchnk, ncol
+ integer :: ncol
+ logical :: pam_stat_fields_active
!----------------------------------------------------------------------------
call phys_getopts(use_ECPP_out = use_ECPP)
call phys_getopts(use_MMF_VT_out = use_MMF_VT)
@@ -366,12 +384,15 @@ subroutine crm_physics_init(state, pbuf2d, species_class)
ncnst = size(cnst_names)
do m = 1, ncnst
call cnst_get_ind(cnst_names(m), mm)
- if ( any(mm == (/ ixcldliq, ixcldice, ixrain, ixsnow /)) ) then
+ if ( any(mm == (/ ixcldliq, ixcldice, ixrain, ixsnow, ixcldrim /)) ) then
! mass mixing ratios
call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm) )
else if ( any(mm == (/ ixnumliq, ixnumice, ixnumrain, ixnumsnow /)) ) then
! number concentrations
call addfld(cnst_name(mm), (/ 'lev' /), 'A', '1/kg', cnst_longname(mm) )
+ else if ( any(mm == (/ ixrimvol /)) ) then
+ ! rime volume
+ call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'm3/kg', cnst_longname(mm) )
else
call endrun( "crm_physics_init: Could not call addfld for constituent with unknown units.")
endif
@@ -392,18 +413,18 @@ subroutine crm_physics_init(state, pbuf2d, species_class)
call addfld(bpcnst(ixcldliq ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixcldliq ))//' before physics' )
call addfld(apcnst(ixcldice ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixcldice ))//' after physics' )
call addfld(bpcnst(ixcldice ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixcldice ))//' before physics' )
- call addfld(apcnst(ixnumliq ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixnumliq ))//' after physics' )
- call addfld(bpcnst(ixnumliq ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixnumliq ))//' before physics' )
- call addfld(apcnst(ixnumice ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixnumice ))//' after physics' )
- call addfld(bpcnst(ixnumice ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixnumice ))//' before physics' )
+ call addfld(apcnst(ixnumliq ),(/'lev'/),'A', '/kg',trim(cnst_name(ixnumliq ))//' after physics' )
+ call addfld(bpcnst(ixnumliq ),(/'lev'/),'A', '/kg',trim(cnst_name(ixnumliq ))//' before physics' )
+ call addfld(apcnst(ixnumice ),(/'lev'/),'A', '/kg',trim(cnst_name(ixnumice ))//' after physics' )
+ call addfld(bpcnst(ixnumice ),(/'lev'/),'A', '/kg',trim(cnst_name(ixnumice ))//' before physics' )
call addfld(apcnst(ixrain ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixrain ))//' after physics' )
call addfld(bpcnst(ixrain ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixrain ))//' before physics' )
- call addfld(apcnst(ixnumrain),(/'lev'/),'A','kg/kg',trim(cnst_name(ixnumrain))//' after physics' )
- call addfld(bpcnst(ixnumrain),(/'lev'/),'A','kg/kg',trim(cnst_name(ixnumrain))//' before physics' )
+ call addfld(apcnst(ixnumrain),(/'lev'/),'A', '/kg',trim(cnst_name(ixnumrain))//' after physics' )
+ call addfld(bpcnst(ixnumrain),(/'lev'/),'A', '/kg',trim(cnst_name(ixnumrain))//' before physics' )
call addfld(apcnst(ixcldrim ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixcldrim ))//' after physics' )
call addfld(bpcnst(ixcldrim ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixcldrim ))//' before physics' )
- call addfld(apcnst(ixrimvol ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixrimvol ))//' after physics' )
- call addfld(bpcnst(ixrimvol ),(/'lev'/),'A','kg/kg',trim(cnst_name(ixrimvol ))//' before physics' )
+ call addfld(apcnst(ixrimvol ),(/'lev'/),'A','m3/kg',trim(cnst_name(ixrimvol ))//' after physics' )
+ call addfld(bpcnst(ixrimvol ),(/'lev'/),'A','m3/kg',trim(cnst_name(ixrimvol ))//' before physics' )
end if
if (use_MMF_VT) then
@@ -473,6 +494,12 @@ end subroutine crm_physics_init
!===================================================================================================
subroutine crm_physics_final()
+#if defined(MMF_PAM)
+ use gator_mod, only: gator_finalize
+ use pam_driver_mod, only: pam_finalize
+ call gator_finalize()
+ call pam_finalize()
+#endif
#if defined(MMF_SAMXX)
use gator_mod, only: gator_finalize
call gator_finalize()
@@ -495,11 +522,14 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
use physics_types, only: physics_state, physics_tend, physics_ptend, physics_ptend_init
use camsrfexch, only: cam_in_t, cam_out_t
use time_manager, only: is_first_step, get_nstep
- use crmdims, only: crm_nx, crm_ny, crm_nz, crm_nx_rad, crm_ny_rad
+ use crmdims, only: crm_nx, crm_ny, crm_nz, crm_nx_rad, crm_ny_rad, crm_dx, crm_dy, crm_dt
use physconst, only: cpair, latvap, latice, gravit, cappa, pi
use constituents, only: pcnst, cnst_get_ind
#if defined(MMF_SAMXX)
use cpp_interface_mod, only: crm
+#elif defined(MMF_PAM)
+ use pam_fortran_interface
+ use pam_driver_mod, only: pam_driver
#elif defined(MMF_SAM) || defined(MMF_SAMOMP)
use crm_module, only: crm
#endif
@@ -523,7 +553,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
use iso_c_binding, only: c_bool
use phys_grid, only: get_rlon_p, get_rlat_p, get_gcol_p
- use spmd_utils, only: masterproc
use openacc_utils, only: prefetch
real(r8), intent(in ) :: ztodt ! global model time increment and CRM run length
@@ -616,6 +645,9 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
real(crm_rknd), dimension(pcols) :: wsx_tmp
real(crm_rknd), dimension(pcols) :: wsy_tmp
+ real(crm_rknd) :: tmp_dz
+ real(crm_rknd) :: tmp_dp
+
! CRM types
type(crm_state_type) :: crm_state
type(crm_rad_type) :: crm_rad
@@ -650,6 +682,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
real(crm_rknd), pointer :: crm_v (:,:,:,:) ! CRM v-wind component
real(crm_rknd), pointer :: crm_w (:,:,:,:) ! CRM w-wind component
real(crm_rknd), pointer :: crm_t (:,:,:,:) ! CRM temperature
+ real(crm_rknd), pointer :: crm_rho_d(:,:,:,:) ! CRM dry density
real(crm_rknd), pointer :: crm_qv(:,:,:,:) ! CRM water vapor mixing ratio (kg/kg of wet air)
real(crm_rknd), pointer :: crm_qp(:,:,:,:) ! 1-mom mass mixing ratio of precipitating condensate
@@ -665,6 +698,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
real(crm_rknd), pointer :: crm_bm(:,:,:,:) ! p3 rime volume
real(crm_rknd), pointer :: crm_q_prev(:,:,:,:) ! p3 previous qv
real(crm_rknd), pointer :: crm_t_prev(:,:,:,:) ! p3 previous t
+ real(crm_rknd), pointer :: crm_shoc_tk(:,:,:,:) ! SHOC Eddy coefficient for momentum [m2/s]
+ real(crm_rknd), pointer :: crm_shoc_tkh(:,:,:,:) ! SHOC Eddy coefficent for heat [m2/s]
+ real(crm_rknd), pointer :: crm_shoc_wthv(:,:,:,:) ! SHOC Buoyancy flux [K m/s]
+ real(crm_rknd), pointer :: crm_shoc_relvar(:,:,:,:) ! SHOC Relative cloud water variance
+ real(crm_rknd), pointer :: crm_shoc_cldfrac(:,:,:,:) ! SHOC Cloud fraction [-]
!------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------
@@ -694,13 +732,6 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
use_crm_accel = use_crm_accel_tmp
crm_accel_uv = crm_accel_uv_tmp
- if (masterproc) then
- if (use_crm_accel .and. trim(MMF_microphysics_scheme)/='sam1mom') then
- write(0,*) "CRM time step relaxation is only compatible with sam1mom microphysics"
- call endrun('crm main')
- endif
- endif
-
nstep = get_nstep()
itim = pbuf_old_tim_idx() ! "Old" pbuf time index (what does all this mean?)
@@ -807,6 +838,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_V'), crm_v)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_W'), crm_w)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_T'), crm_t)
+ call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_RHO'),crm_rho_d)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_QV'), crm_qv)
if (MMF_microphysics_scheme .eq. 'sam1mom') then
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_QP'), crm_qp)
@@ -822,6 +854,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_BM'), crm_bm)
call pbuf_get_field(pbuf_chunk, crm_t_prev_idx, crm_t_prev)
call pbuf_get_field(pbuf_chunk, crm_q_prev_idx, crm_q_prev)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_tk_idx , crm_shoc_tk)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_tkh_idx , crm_shoc_tkh)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_wthv_idx , crm_shoc_wthv)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_relvar_idx , crm_shoc_relvar)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_cldfrac_idx, crm_shoc_cldfrac)
end if
! initialize all water to zero (needed for ncol < i <= pcols)
@@ -849,7 +886,14 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
crm_w(i,:,:,k) = 0.
crm_t(i,:,:,k) = state(c)%t(i,m)
- ! Initialize microphysics arrays
+ ! diagnose dry density for PAM using:
+ ! pdel = pdel_dry + pdel*qv
+ ! pdel_dry = rho_dry*dz*g
+ tmp_dp = state(c)%pint(i,m+1) - state(c)%pint(i,m)
+ tmp_dz = state(c)%zi(i,m+1) - state(c)%zi(i,m)
+ crm_rho_d(i,:,:,k) = -1 * tmp_dp * ( 1 - state(c)%q(i,m,1) ) / ( tmp_dz * gravit )
+
+ ! Initialize microphysics arrays
if (MMF_microphysics_scheme .eq. 'sam1mom') then
crm_qv(i,:,:,k) = state(c)%q(i,m,1)!+state(c)%q(i,m,ixcldliq)+state(c)%q(i,m,ixcldice)
crm_qp(i,:,:,k) = 0.0_r8
@@ -866,6 +910,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
crm_bm(i,:,:,k) = 0.0_r8
crm_t_prev(i,:,:,k) = state(c)%t(i,m)
crm_q_prev(i,:,:,k) = state(c)%q(i,m,1)
+ crm_shoc_tk (i,:,:,k) = 0.0_r8
+ crm_shoc_tkh (i,:,:,k) = 0.0_r8
+ crm_shoc_wthv (i,:,:,k) = 0.0_r8
+ crm_shoc_relvar (i,:,:,k) = 0.0_r8
+ crm_shoc_cldfrac(i,:,:,k) = 0.0_r8
end if
end do
@@ -950,6 +999,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_V'), crm_v)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_W'), crm_w)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_T'), crm_t)
+ call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_RHO'),crm_rho_d)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_QV'), crm_qv)
if (MMF_microphysics_scheme .eq. 'sam1mom') then
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_QP'), crm_qp)
@@ -966,6 +1016,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_BM'), crm_bm)
call pbuf_get_field(pbuf_chunk, crm_t_prev_idx, crm_t_prev)
call pbuf_get_field(pbuf_chunk, crm_q_prev_idx, crm_q_prev)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_tk_idx , crm_shoc_tk)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_tkh_idx , crm_shoc_tkh)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_wthv_idx , crm_shoc_wthv)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_relvar_idx , crm_shoc_relvar)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_cldfrac_idx, crm_shoc_cldfrac)
end if
! copy pbuf data into crm_state
@@ -975,6 +1030,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
crm_state%v_wind (icrm,:,:,:) = crm_v (i,:,:,:)
crm_state%w_wind (icrm,:,:,:) = crm_w (i,:,:,:)
crm_state%temperature(icrm,:,:,:) = crm_t (i,:,:,:)
+ crm_state%rho_dry (icrm,:,:,:) = crm_rho_d(i,:,:,:)
crm_state%qv (icrm,:,:,:) = crm_qv(i,:,:,:)
if (MMF_microphysics_scheme .eq. 'sam1mom') then
crm_state%qp (icrm,:,:,:) = crm_qp(i,:,:,:)
@@ -991,16 +1047,23 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
crm_state%bm (icrm,:,:,:) = crm_bm(i,:,:,:)
crm_state%t_prev (icrm,:,:,:) = crm_t_prev(i,:,:,:)
crm_state%q_prev (icrm,:,:,:) = crm_q_prev(i,:,:,:)
+ crm_state%shoc_tk (icrm,:,:,:) = crm_shoc_tk (i,:,:,:)
+ crm_state%shoc_tkh (icrm,:,:,:) = crm_shoc_tkh (i,:,:,:)
+ crm_state%shoc_wthv (icrm,:,:,:) = crm_shoc_wthv (i,:,:,:)
+ crm_state%shoc_relvar (icrm,:,:,:) = crm_shoc_relvar (i,:,:,:)
+ crm_state%shoc_cldfrac(icrm,:,:,:) = crm_shoc_cldfrac(i,:,:,:)
end if
end do ! i=1,ncol
! Retrieve radiative heating tendency
+ ! Rad tendency was normalized for energy conservation
+ ! un-normalize before sending to CRM - this yields units of [K/s]
call pbuf_get_field(pbuf_chunk, crm_qrad_idx, crm_qrad)
do m = 1,crm_nz
k = pver-m+1
do i = 1,ncol
icrm = ncol_sum + i
- crm_rad%qrad(icrm,:,:,m) = crm_qrad(i,:,:,m) / state(c)%pdel(i,k) ! normalize for energy conservation
+ crm_rad%qrad(icrm,:,:,m) = crm_qrad(i,:,:,m) / state(c)%pdel(i,k)
end do
end do
@@ -1020,7 +1083,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
do ii = 1,crm_nx
if (MMF_microphysics_scheme .eq. 'p3') then
qli_hydro_before(c,i) = qli_hydro_before(c,i)+(crm_state%qr(icrm,ii,jj,m)) * dp_g
- ! TODO: how do we handle qi_hydro_before for P3?
+ ! P3 does not have a separate category for snow, so qi_hydro_before should be zero
end if
if (MMF_microphysics_scheme .eq. 'sam1mom') then
sfactor = max(0._r8,min(1._r8,(crm_state%temperature(icrm,ii,jj,m)-268.16)*1./(283.16-268.16)))
@@ -1075,15 +1138,18 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
! P3 input data
!------------------------------------------------------------------------------------------
if (MMF_microphysics_scheme .eq. 'p3') then
+ ! NOTE - nccn_prescribed is only used when do_prescribed_CCN=true
+ ! if do_prescribed_CCN=false then nc_nuceat_tend and ni_activated
+ ! must be provided by an aerosol activation scheme
+ ! TODO: build interface for SPA to get nccn_prescribed consistent with SCREAM
do i = 1,ncol
icrm = ncol_sum + i
do k = 1, pver
- crm_input%nccn (icrm,k) = 1e3
- crm_input%nc_nuceat_tend(icrm,k) = 1.0 ! npccn (i,l)
- crm_input%ni_activated (icrm,k) = 1.0 ! ni_activated(i,l)
+ crm_input%nccn_prescribed(icrm,k) = 50e6 ! 50 [#/cm3] 1e6 [cm3/m3] / (kg/m3) => 50e6 [#/kg]
+ crm_input%nc_nuceat_tend (icrm,k) = 0.0 ! [#/kg/s]
+ crm_input%ni_activated (icrm,k) = 0.0 ! [#/kg]
end do
end do
-
end if
!------------------------------------------------------------------------------------------
@@ -1109,10 +1175,10 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
icrm = ncol_sum + i
crm_input%tau00(icrm) = sqrt(wsx_tmp(i)**2 + wsy_tmp(i)**2)
crm_input%bflxls(icrm) = shf_tmp(i)/cpair + 0.61*state(c)%t(i,pver)*lhf_tmp(i)/latvap
- crm_input%fluxu00(icrm) = wsx_tmp(i) ! N/m2
- crm_input%fluxv00(icrm) = wsy_tmp(i) ! N/m2
- crm_input%fluxt00(icrm) = shf_tmp(i)/cpair ! K Kg/ (m2 s)
- crm_input%fluxq00(icrm) = lhf_tmp(i)/latvap ! Kg/(m2 s)
+ crm_input%fluxu00(icrm) = wsx_tmp(i) ! N/m2
+ crm_input%fluxv00(icrm) = wsy_tmp(i) ! N/m2
+ crm_input%fluxt00(icrm) = shf_tmp(i) ! W/m2 ( divide by cpair to get K kg/(m2 s) )
+ crm_input%fluxq00(icrm) = lhf_tmp(i) ! W/m2 ( divide by latvap to get kg/(m2 s) )
crm_input%wndls(icrm) = sqrt(state(c)%u(i,pver)**2 + state(c)%v(i,pver)**2)
end do
@@ -1126,11 +1192,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
air_density(i,1:pver) = state(c)%pmid(i,1:pver) / (287.15*state(c)%t(i,1:pver))
do k = 1, pver
do m = 1, ntot_amode
- call loadaer( state(c), pbuf_chunk, i, i, k, m, air_density, phase, &
- aerosol_num, aerosol_vol, aerosol_hygro)
- crm_input%naermod (icrm,k,m) = aerosol_num(i)
- crm_input%vaerosol(icrm,k,m) = aerosol_vol(i)
- crm_input%hygro (icrm,k,m) = aerosol_hygro(i)
+ call loadaer( state(c), pbuf_chunk, i, i, k, m, air_density, phase, &
+ aerosol_num, aerosol_vol, aerosol_hygro)
+ crm_input%naermod (icrm,k,m) = aerosol_num(i)
+ crm_input%vaerosol(icrm,k,m) = aerosol_vol(i)
+ crm_input%hygro (icrm,k,m) = aerosol_hygro(i)
end do
end do
end do
@@ -1164,7 +1230,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
end do ! c=begchunk, endchunk
#if defined(MMF_SAM) || defined(MMF_SAMOMP)
-
+
call t_startf ('crm_call')
call crm(ncrms, ztodt, pver, &
crm_input, crm_state, crm_rad, &
@@ -1173,12 +1239,12 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
use_MMF_VT_tmp, MMF_VT_wn_max, &
use_crm_accel_tmp, crm_accel_factor, crm_accel_uv_tmp)
call t_stopf('crm_call')
-
+
#elif defined(MMF_SAMXX)
- call t_startf ('crm_call')
! Fortran classes don't translate to C++ classes, we we have to separate
! this stuff out when calling the C++ routinte crm(...)
+ call t_startf ('crm_call')
call crm(ncrms, ncrms, ztodt, pver, crm_input%bflxls, crm_input%wndls, crm_input%zmid, crm_input%zint, &
crm_input%pmid, crm_input%pint, crm_input%pdel, crm_input%ul, crm_input%vl, &
crm_input%tl, crm_input%qccl, crm_input%qiil, crm_input%ql, crm_input%tau00, &
@@ -1209,6 +1275,198 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
use_crm_accel, crm_accel_factor, crm_accel_uv)
call t_stopf('crm_call')
+#elif defined(MMF_PAM)
+
+ call pam_mirror_array_readonly( 'latitude', latitude0 )
+ call pam_mirror_array_readonly( 'longitude', longitude0 )
+
+ call pam_mirror_array_readonly( 'input_bflxls', crm_input%bflxls )
+ call pam_mirror_array_readonly( 'input_wndls', crm_input%wndls )
+
+ call pam_mirror_array_readonly( 'input_shf', crm_input%fluxt00 )
+ call pam_mirror_array_readonly( 'input_lhf', crm_input%fluxq00 )
+
+ call pam_mirror_array_readonly( 'input_zmid', crm_input%zmid )
+ call pam_mirror_array_readonly( 'input_zint', crm_input%zint )
+ call pam_mirror_array_readonly( 'input_pmid', crm_input%pmid )
+ call pam_mirror_array_readonly( 'input_pint', crm_input%pint )
+ call pam_mirror_array_readonly( 'input_pdel', crm_input%pdel )
+ call pam_mirror_array_readonly( 'input_phis', crm_input%phis )
+
+ call pam_mirror_array_readonly( 'input_ul', crm_input%ul )
+ call pam_mirror_array_readonly( 'input_vl', crm_input%vl )
+ call pam_mirror_array_readonly( 'input_tl', crm_input%tl )
+ call pam_mirror_array_readonly( 'input_qccl', crm_input%qccl )
+ call pam_mirror_array_readonly( 'input_qiil', crm_input%qiil )
+ call pam_mirror_array_readonly( 'input_ql', crm_input%ql )
+ call pam_mirror_array_readonly( 'input_tau00', crm_input%tau00 )
+ ! call pam_mirror_array_readonly( 'input_ul_esmt', crm_input%ul_esmt )
+ ! call pam_mirror_array_readonly( 'input_vl_esmt', crm_input%vl_esmt )
+ ! call pam_mirror_array_readonly( 'input_t_vt', crm_input%t_vt )
+ ! call pam_mirror_array_readonly( 'input_q_vt', crm_input%q_vt )
+ ! call pam_mirror_array_readonly( 'input_u_vt', crm_input%u_vt )
+
+ call pam_mirror_array_readonly( 'input_nccn_prescribed',crm_input%nccn_prescribed )
+ call pam_mirror_array_readonly( 'input_nc_nuceat_tend', crm_input%nc_nuceat_tend )
+ call pam_mirror_array_readonly( 'input_ni_activated', crm_input%ni_activated )
+
+ call pam_mirror_array_readwrite( 'state_u_wind', crm_state%u_wind )
+ call pam_mirror_array_readwrite( 'state_v_wind', crm_state%v_wind )
+ call pam_mirror_array_readwrite( 'state_w_wind', crm_state%w_wind )
+ call pam_mirror_array_readwrite( 'state_temperature', crm_state%temperature )
+ call pam_mirror_array_readwrite( 'state_rho_dry', crm_state%rho_dry )
+ call pam_mirror_array_readwrite( 'state_qv', crm_state%qv )
+ call pam_mirror_array_readwrite( 'state_qc', crm_state%qc )
+ call pam_mirror_array_readwrite( 'state_nc', crm_state%nc )
+ call pam_mirror_array_readwrite( 'state_qr', crm_state%qr )
+ call pam_mirror_array_readwrite( 'state_nr', crm_state%nr )
+ call pam_mirror_array_readwrite( 'state_qi', crm_state%qi )
+ call pam_mirror_array_readwrite( 'state_ni', crm_state%ni )
+ call pam_mirror_array_readwrite( 'state_qm', crm_state%qm )
+ call pam_mirror_array_readwrite( 'state_bm', crm_state%bm )
+ call pam_mirror_array_readwrite( 'state_t_prev', crm_state%t_prev )
+ call pam_mirror_array_readwrite( 'state_q_prev', crm_state%q_prev )
+
+ ! SHOC variables
+ call pam_mirror_array_readwrite( 'state_shoc_tk', crm_state%shoc_tk )
+ call pam_mirror_array_readwrite( 'state_shoc_tkh', crm_state%shoc_tkh )
+ call pam_mirror_array_readwrite( 'state_shoc_wthv', crm_state%shoc_wthv )
+ call pam_mirror_array_readwrite( 'state_shoc_relvar', crm_state%shoc_relvar )
+ call pam_mirror_array_readwrite( 'state_shoc_cldfrac', crm_state%shoc_cldfrac )
+
+ ! Radiation tendency and output conditions
+ call pam_mirror_array_readwrite( 'rad_qrad', crm_rad%qrad )
+ call pam_mirror_array_readwrite( 'rad_temperature', crm_rad%temperature )
+ call pam_mirror_array_readwrite( 'rad_qv', crm_rad%qv )
+ call pam_mirror_array_readwrite( 'rad_qc', crm_rad%qc )
+ call pam_mirror_array_readwrite( 'rad_qi', crm_rad%qi )
+ call pam_mirror_array_readwrite( 'rad_nc', crm_rad%nc )
+ call pam_mirror_array_readwrite( 'rad_ni', crm_rad%ni )
+ call pam_mirror_array_readwrite( 'rad_cld', crm_rad%cld )
+
+ call pam_mirror_array_readwrite( 'output_clear_rh', crm_clear_rh )
+
+ ! call pam_mirror_array_readwrite( 'output_prectend', crm_output%prectend, '' )
+ ! call pam_mirror_array_readwrite( 'output_precstend', crm_output%precstend, '' )
+ call pam_mirror_array_readwrite( 'output_cld', crm_output%cld, '' )
+ ! call pam_mirror_array_readwrite( 'output_cldtop', crm_output%cldtop, '' )
+ call pam_mirror_array_readwrite( 'output_gicewp', crm_output%gicewp, '' )
+ call pam_mirror_array_readwrite( 'output_gliqwp', crm_output%gliqwp, '' )
+ call pam_mirror_array_readwrite( 'output_liq_ice_exchange',crm_output%liq_ice_exchange,'' )
+ call pam_mirror_array_readwrite( 'output_vap_liq_exchange',crm_output%vap_liq_exchange,'' )
+ call pam_mirror_array_readwrite( 'output_vap_ice_exchange',crm_output%vap_ice_exchange,'' )
+ ! call pam_mirror_array_readwrite( 'output_mctot', crm_output%mctot, '' )
+ ! call pam_mirror_array_readwrite( 'output_mcup', crm_output%mcup, '' )
+ ! call pam_mirror_array_readwrite( 'output_mcdn', crm_output%mcdn, '' )
+ ! call pam_mirror_array_readwrite( 'output_mcuup', crm_output%mcuup, '' )
+ ! call pam_mirror_array_readwrite( 'output_mcudn', crm_output%mcudn, '' )
+
+ call pam_mirror_array_readwrite( 'output_qv_mean', crm_output%qv_mean, '' )
+ call pam_mirror_array_readwrite( 'output_qc_mean', crm_output%qc_mean, '' )
+ call pam_mirror_array_readwrite( 'output_qr_mean', crm_output%qr_mean, '' )
+ call pam_mirror_array_readwrite( 'output_qi_mean', crm_output%qi_mean, '' )
+ call pam_mirror_array_readwrite( 'output_nc_mean', crm_output%nc_mean, '' )
+ call pam_mirror_array_readwrite( 'output_nr_mean', crm_output%nr_mean, '' )
+ call pam_mirror_array_readwrite( 'output_ni_mean', crm_output%ni_mean, '' )
+ call pam_mirror_array_readwrite( 'output_qm_mean', crm_output%qm_mean, '' )
+ call pam_mirror_array_readwrite( 'output_bm_mean', crm_output%bm_mean, '' )
+ call pam_mirror_array_readwrite( 'output_rho_d_mean', crm_output%rho_d_mean, '' )
+ call pam_mirror_array_readwrite( 'output_rho_v_mean', crm_output%rho_v_mean, '' )
+
+ call pam_mirror_array_readwrite( 'output_qt_ls', crm_output%qt_ls, '' )
+ call pam_mirror_array_readwrite( 'output_t_ls', crm_output%t_ls, '' )
+ call pam_mirror_array_readwrite( 'output_rho_v_ls', crm_output%rho_v_ls, '' )
+ call pam_mirror_array_readwrite( 'output_rho_d_ls', crm_output%rho_d_ls, '' )
+ call pam_mirror_array_readwrite( 'output_rho_l_ls', crm_output%rho_l_ls, '' )
+ call pam_mirror_array_readwrite( 'output_rho_i_ls', crm_output%rho_i_ls, '' )
+ ! call pam_mirror_array_readwrite( 'output_jt_crm', crm_output%jt_crm, '' )
+ ! call pam_mirror_array_readwrite( 'output_mx_crm', crm_output%mx_crm, '' )
+ ! call pam_mirror_array_readwrite( 'output_cltot', crm_output%cltot, '' )
+ ! call pam_mirror_array_readwrite( 'output_clhgh', crm_output%clhgh, '' )
+ ! call pam_mirror_array_readwrite( 'output_clmed', crm_output%clmed, '' )
+ ! call pam_mirror_array_readwrite( 'output_cllow', crm_output%cllow, '' )
+ call pam_mirror_array_readwrite( 'output_sltend', crm_output%sltend, '' )
+ call pam_mirror_array_readwrite( 'output_qltend', crm_output%qltend, '' )
+ call pam_mirror_array_readwrite( 'output_qcltend', crm_output%qcltend, '' )
+ call pam_mirror_array_readwrite( 'output_qiltend', crm_output%qiltend, '' )
+ ! call pam_mirror_array_readwrite( 'output_t_vt_tend', crm_output%t_vt_tend, '' )
+ ! call pam_mirror_array_readwrite( 'output_q_vt_tend', crm_output%q_vt_tend, '' )
+ ! call pam_mirror_array_readwrite( 'output_u_vt_tend', crm_output%u_vt_tend, '' )
+ ! call pam_mirror_array_readwrite( 'output_t_vt_ls', crm_output%t_vt_ls, '' )
+ ! call pam_mirror_array_readwrite( 'output_q_vt_ls', crm_output%q_vt_ls, '' )
+ ! call pam_mirror_array_readwrite( 'output_u_vt_ls', crm_output%u_vt_ls, '' )
+ call pam_mirror_array_readwrite( 'output_ultend', crm_output%ultend, '' )
+ call pam_mirror_array_readwrite( 'output_vltend', crm_output%vltend, '' )
+ ! call pam_mirror_array_readwrite( 'output_tk', crm_output%tk, '' )
+ ! call pam_mirror_array_readwrite( 'output_tkh', crm_output%tkh, '' )
+ ! call pam_mirror_array_readwrite( 'output_qcl', crm_output%qcl, '' )
+ ! call pam_mirror_array_readwrite( 'output_qci', crm_output%qci, '' )
+ ! call pam_mirror_array_readwrite( 'output_qpl', crm_output%qpl, '' )
+ ! call pam_mirror_array_readwrite( 'output_qpi', crm_output%qpi, '' )
+ call pam_mirror_array_readwrite( 'output_precc', crm_output%precc, '' )
+ call pam_mirror_array_readwrite( 'output_precl', crm_output%precl, '' )
+ call pam_mirror_array_readwrite( 'output_precsc', crm_output%precsc, '' )
+ call pam_mirror_array_readwrite( 'output_precsl', crm_output%precsl, '' )
+ call pam_mirror_array_readwrite( 'output_prec_crm', crm_output%prec_crm, '' )
+
+ call pam_mirror_array_readwrite( 'output_dt_sgs', crm_output%dt_sgs, '' )
+ call pam_mirror_array_readwrite( 'output_dqv_sgs', crm_output%dqv_sgs, '' )
+ call pam_mirror_array_readwrite( 'output_dqc_sgs', crm_output%dqc_sgs, '' )
+ call pam_mirror_array_readwrite( 'output_dqi_sgs', crm_output%dqi_sgs, '' )
+ call pam_mirror_array_readwrite( 'output_dqr_sgs', crm_output%dqr_sgs, '' )
+
+ call pam_mirror_array_readwrite( 'output_dt_micro', crm_output%dt_micro, '' )
+ call pam_mirror_array_readwrite( 'output_dqv_micro', crm_output%dqv_micro, '' )
+ call pam_mirror_array_readwrite( 'output_dqc_micro', crm_output%dqc_micro, '' )
+ call pam_mirror_array_readwrite( 'output_dqi_micro', crm_output%dqi_micro, '' )
+ call pam_mirror_array_readwrite( 'output_dqr_micro', crm_output%dqr_micro, '' )
+
+ call pam_mirror_array_readwrite( 'output_dt_dycor', crm_output%dt_dycor, '' )
+ call pam_mirror_array_readwrite( 'output_dqv_dycor', crm_output%dqv_dycor, '' )
+ call pam_mirror_array_readwrite( 'output_dqc_dycor', crm_output%dqc_dycor, '' )
+ call pam_mirror_array_readwrite( 'output_dqi_dycor', crm_output%dqi_dycor, '' )
+ call pam_mirror_array_readwrite( 'output_dqr_dycor', crm_output%dqr_dycor, '' )
+
+ call pam_mirror_array_readwrite( 'output_dt_sponge', crm_output%dt_sponge, '' )
+ call pam_mirror_array_readwrite( 'output_dqv_sponge', crm_output%dqv_sponge, '' )
+ call pam_mirror_array_readwrite( 'output_dqc_sponge', crm_output%dqc_sponge, '' )
+ call pam_mirror_array_readwrite( 'output_dqi_sponge', crm_output%dqi_sponge, '' )
+ call pam_mirror_array_readwrite( 'output_dqr_sponge', crm_output%dqr_sponge, '' )
+
+ call pam_mirror_array_readonly( 'global_column_id', gcolp )
+
+ call pam_set_option('ncrms', ncrms )
+ call pam_set_option('gcm_nlev', pver )
+ call pam_set_option('crm_nz',crm_nz )
+ call pam_set_option('crm_nx',crm_nx )
+ call pam_set_option('crm_ny',crm_ny )
+ call pam_set_option('rad_nx',crm_nx_rad)
+ call pam_set_option('rad_ny',crm_ny_rad)
+ call pam_set_option('crm_dx',crm_dx )
+ call pam_set_option('crm_dy',crm_dy )
+ call pam_set_option('gcm_dt',ztodt )
+ call pam_set_option('crm_dt',crm_dt )
+
+ call pam_register_dimension('gcm_lev',pver)
+
+ call pam_set_option('use_MMF_VT', use_MMF_VT_tmp )
+ call pam_set_option('MMF_VT_wn_max', MMF_VT_wn_max )
+ call pam_set_option('use_MMF_ESMT', use_MMF_ESMT_tmp )
+
+ call pam_set_option('use_crm_accel', use_crm_accel_tmp )
+ call pam_set_option('crm_accel_uv', crm_accel_uv_tmp)
+ call pam_set_option('crm_accel_factor', crm_accel_factor )
+
+ call pam_set_option('enable_physics_tend_stats', .false. )
+
+ call pam_set_option('is_first_step', (nstep<=1) )
+ call pam_set_option('is_restart', (nsrest>0) )
+ call pam_set_option('am_i_root', masterproc )
+
+ call t_startf ('crm_call')
+ call pam_driver()
+ call t_stopf('crm_call')
+
#endif
deallocate(longitude0)
@@ -1273,12 +1531,22 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
ptend(c)%lq(ixnumliq) = .TRUE.
ptend(c)%lq(ixnumice) = .TRUE.
- if (use_ECPP) then
- ptend(c)%lq(ixrain) = .TRUE.
- ptend(c)%lq(ixnumrain) = .TRUE.
- ptend(c)%lq(ixcldrim) = .TRUE.
- ptend(c)%lq(ixrimvol) = .TRUE.
- end if
+ ptend(c)%lq(ixrain) = .TRUE.
+ ptend(c)%lq(ixnumrain) = .TRUE.
+ ptend(c)%lq(ixcldrim) = .TRUE.
+ ptend(c)%lq(ixrimvol) = .TRUE.
+
+ do i = 1, ncol
+ do k = 1, crm_nz
+ m = pver-k+1
+ ptend(c)%q(i,m,ixnumliq) = 0
+ ptend(c)%q(i,m,ixnumice) = 0
+ ptend(c)%q(i,m,ixrain) = 0
+ ptend(c)%q(i,m,ixnumrain) = 0
+ ptend(c)%q(i,m,ixcldrim) = 0
+ ptend(c)%q(i,m,ixrimvol) = 0
+ end do
+ end do
do i = 1, ncol
icrm = ncol_sum + i
@@ -1286,24 +1554,20 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
m = pver-k+1
do ii = 1, crm_nx
do jj = 1, crm_ny
- ptend(c)%q(i,m,ixnumliq) = ptend(c)%q(i,m,ixnumliq) + crm_state%nc(icrm,ii,jj,k)
- ptend(c)%q(i,m,ixnumice) = ptend(c)%q(i,m,ixnumice) + crm_state%ni(icrm,ii,jj,k)
- if (use_ECPP) then
- ptend(c)%q(i,m,ixrain) = ptend(c)%q(i,m,ixrain) + crm_state%qr(icrm,ii,jj,k)
- ptend(c)%q(i,m,ixnumrain) = ptend(c)%q(i,m,ixnumrain) + crm_state%nr(icrm,ii,jj,k)
- ptend(c)%q(i,m,ixcldrim) = ptend(c)%q(i,m,ixcldrim) + crm_state%qm(icrm,ii,jj,k)
- ptend(c)%q(i,m,ixrimvol) = ptend(c)%q(i,m,ixrimvol) + crm_state%bm(icrm,ii,jj,k)
- end if
+ ptend(c)%q(i,m,ixnumliq) = ptend(c)%q(i,m,ixnumliq ) + crm_state%nc(icrm,ii,jj,k)
+ ptend(c)%q(i,m,ixnumice) = ptend(c)%q(i,m,ixnumice ) + crm_state%ni(icrm,ii,jj,k)
+ ptend(c)%q(i,m,ixrain) = ptend(c)%q(i,m,ixrain ) + crm_state%qr(icrm,ii,jj,k)
+ ptend(c)%q(i,m,ixnumrain) = ptend(c)%q(i,m,ixnumrain) + crm_state%nr(icrm,ii,jj,k)
+ ptend(c)%q(i,m,ixcldrim) = ptend(c)%q(i,m,ixcldrim ) + crm_state%qm(icrm,ii,jj,k)
+ ptend(c)%q(i,m,ixrimvol) = ptend(c)%q(i,m,ixrimvol ) + crm_state%bm(icrm,ii,jj,k)
end do
end do
- ptend(c)%q(i,m,ixnumliq) = (ptend(c)%q(i,m,ixnumliq) /(crm_nx*crm_ny) - state(c)%q(i,m,ixnumliq)) /ztodt
- ptend(c)%q(i,m,ixnumice) = (ptend(c)%q(i,m,ixnumice) /(crm_nx*crm_ny) - state(c)%q(i,m,ixnumice)) /ztodt
- if (use_ECPP) then
- ptend(c)%q(i,m,ixrain) = (ptend(c)%q(i,m,ixrain) /(crm_nx*crm_ny) - state(c)%q(i,m,ixrain)) /ztodt
- ptend(c)%q(i,m,ixnumrain) = (ptend(c)%q(i,m,ixnumrain)/(crm_nx*crm_ny) - state(c)%q(i,m,ixnumrain))/ztodt
- ptend(c)%q(i,m,ixcldrim) = (ptend(c)%q(i,m,ixcldrim) /(crm_nx*crm_ny) - state(c)%q(i,m,ixcldrim)) /ztodt
- ptend(c)%q(i,m,ixrimvol) = (ptend(c)%q(i,m,ixrimvol) /(crm_nx*crm_ny) - state(c)%q(i,m,ixrimvol)) /ztodt
- end if
+ ptend(c)%q(i,m,ixnumliq ) = ( ptend(c)%q(i,m,ixnumliq )/(crm_nx*crm_ny) - state(c)%q(i,m,ixnumliq ) )/ztodt
+ ptend(c)%q(i,m,ixnumice ) = ( ptend(c)%q(i,m,ixnumice )/(crm_nx*crm_ny) - state(c)%q(i,m,ixnumice ) )/ztodt
+ ptend(c)%q(i,m,ixrain ) = ( ptend(c)%q(i,m,ixrain )/(crm_nx*crm_ny) - state(c)%q(i,m,ixrain ) )/ztodt
+ ptend(c)%q(i,m,ixnumrain) = ( ptend(c)%q(i,m,ixnumrain)/(crm_nx*crm_ny) - state(c)%q(i,m,ixnumrain) )/ztodt
+ ptend(c)%q(i,m,ixcldrim ) = ( ptend(c)%q(i,m,ixcldrim )/(crm_nx*crm_ny) - state(c)%q(i,m,ixcldrim ) )/ztodt
+ ptend(c)%q(i,m,ixrimvol ) = ( ptend(c)%q(i,m,ixrimvol )/(crm_nx*crm_ny) - state(c)%q(i,m,ixrimvol ) )/ztodt
end do
end do
@@ -1330,9 +1594,14 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
end do
end do
+#if defined(MMF_PAM)
+ ! Currently PAM is coupled the CRM to all levels, so only add rad to levels above CRM
+ ptend(c)%s(1:ncol, 1:pver-crm_nz) = qrs(1:ncol,1:pver-crm_nz) + qrl(1:ncol,1:pver-crm_nz)
+#else
! The radiation tendencies in the GCM levels above the CRM and the top 2 CRM levels are set to
! be zero in the CRM, So add radiation tendencies to these levels
ptend(c)%s(1:ncol, 1:pver-crm_nz+2) = qrs(1:ncol,1:pver-crm_nz+2) + qrl(1:ncol,1:pver-crm_nz+2)
+#endif
! This will be used to check energy conservation
mmf_rad_flux(c,:ncol) = 0.0_r8
@@ -1462,6 +1731,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_V'), crm_v)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_W'), crm_w)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_T'), crm_t)
+ call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_RHO'), crm_rho_d)
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_QV'), crm_qv)
if (MMF_microphysics_scheme .eq. 'sam1mom') then
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_QP'), crm_qp)
@@ -1478,6 +1748,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
call pbuf_get_field(pbuf_chunk, pbuf_get_index('CRM_BM'), crm_bm)
call pbuf_get_field(pbuf_chunk, crm_t_prev_idx, crm_t_prev)
call pbuf_get_field(pbuf_chunk, crm_q_prev_idx, crm_q_prev)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_tk_idx , crm_shoc_tk)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_tkh_idx , crm_shoc_tkh)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_wthv_idx , crm_shoc_wthv)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_relvar_idx , crm_shoc_relvar)
+ call pbuf_get_field(pbuf_chunk, crm_shoc_cldfrac_idx, crm_shoc_cldfrac)
end if
do i = 1,ncol
@@ -1486,6 +1761,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
crm_v (i,:,:,:) = crm_state%v_wind (icrm,:,:,:)
crm_w (i,:,:,:) = crm_state%w_wind (icrm,:,:,:)
crm_t (i,:,:,:) = crm_state%temperature(icrm,:,:,:)
+ crm_rho_d(i,:,:,:) = crm_state%rho_dry (icrm,:,:,:)
crm_qv(i,:,:,:) = crm_state%qv (icrm,:,:,:)
if (MMF_microphysics_scheme .eq. 'sam1mom') then
crm_qp(i,:,:,:) = crm_state%qp(icrm,:,:,:)
@@ -1502,6 +1778,11 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
crm_bm(i,:,:,:) = crm_state%bm (icrm,:,:,:)
crm_t_prev(i,:,:,:) = crm_state%t_prev(icrm,:,:,:)
crm_q_prev(i,:,:,:) = crm_state%q_prev(icrm,:,:,:)
+ crm_shoc_tk (i,:,:,:) = crm_state%shoc_tk (icrm,:,:,:)
+ crm_shoc_tkh (i,:,:,:) = crm_state%shoc_tkh (icrm,:,:,:)
+ crm_shoc_wthv (i,:,:,:) = crm_state%shoc_wthv (icrm,:,:,:)
+ crm_shoc_relvar (i,:,:,:) = crm_state%shoc_relvar (icrm,:,:,:)
+ crm_shoc_cldfrac(i,:,:,:) = crm_state%shoc_cldfrac(icrm,:,:,:)
end if
end do
@@ -1521,7 +1802,7 @@ subroutine crm_physics_tend(ztodt, state, tend, ptend, pbuf2d, cam_in, cam_out,
do ii = 1,crm_nx
if(MMF_microphysics_scheme .eq. 'p3') then
qli_hydro_after(c,i) = qli_hydro_after(c,i)+(crm_state%qr(icrm,ii,jj,m)) * dp_g
- ! TODO: how do we handle qi_hydro_after for P3?
+ ! P3 does not have a separate category for snow, so qi_hydro_after should be zero
end if
if(MMF_microphysics_scheme .eq. 'sam1mom') then
sfactor = max(0._r8,min(1._r8,(crm_state%temperature(icrm,ii,jj,m)-268.16)*1./(283.16-268.16)))
@@ -1628,3 +1909,4 @@ end subroutine crm_surface_flux_bypass_tend
!==================================================================================================
end module crm_physics
+
\ No newline at end of file
diff --git a/components/eam/src/physics/crm/crm_state_module.F90 b/components/eam/src/physics/crm/crm_state_module.F90
index 82db8d682466..63aea4648039 100644
--- a/components/eam/src/physics/crm/crm_state_module.F90
+++ b/components/eam/src/physics/crm/crm_state_module.F90
@@ -25,6 +25,7 @@ module crm_state_module
real(crm_rknd), allocatable :: v_wind(:,:,:,:) ! CRM v-wind component
real(crm_rknd), allocatable :: w_wind(:,:,:,:) ! CRM w-wind component
real(crm_rknd), allocatable :: temperature(:,:,:,:) ! CRM temperature
+ real(crm_rknd), allocatable :: rho_dry(:,:,:,:) ! CRM dry density
real(crm_rknd), allocatable :: qv(:,:,:,:) ! CRM water vapor
! 1-moment microphsics variables
@@ -44,9 +45,16 @@ module crm_state_module
real(crm_rknd), allocatable :: bm(:,:,:,:) ! averaged riming volume
! "previous" state variables needed for P3
- real(crm_rknd), allocatable :: t_prev(:,:,:,:) ! previous CRM time step temperature
+ real(crm_rknd), allocatable :: t_prev(:,:,:,:) ! previous CRM time step temperature
real(crm_rknd), allocatable :: q_prev(:,:,:,:) ! previous CRM time step water vapor
+ ! SHOC quantities that need to persist between CRM calls
+ real(crm_rknd), allocatable :: shoc_tk (:,:,:,:) ! eddy coefficient for momentum [m2/s]
+ real(crm_rknd), allocatable :: shoc_tkh (:,:,:,:) ! eddy coefficient for heat [m2/s]
+ real(crm_rknd), allocatable :: shoc_wthv (:,:,:,:) ! buoyancy flux [K m/s]
+ real(crm_rknd), allocatable :: shoc_relvar (:,:,:,:) ! relative cloud water variance
+ real(crm_rknd), allocatable :: shoc_cldfrac(:,:,:,:) ! Cloud fraction
+
end type crm_state_type
!------------------------------------------------------------------------------------------------
contains
@@ -63,12 +71,14 @@ subroutine crm_state_initialize(state,ncrms,crm_nx,crm_ny,crm_nz,MMF_microphysic
if (.not. allocated(state%v_wind)) allocate(state%v_wind(ncrms,crm_nx,crm_ny,crm_nz))
if (.not. allocated(state%w_wind)) allocate(state%w_wind(ncrms,crm_nx,crm_ny,crm_nz))
if (.not. allocated(state%temperature)) allocate(state%temperature(ncrms,crm_nx,crm_ny,crm_nz))
+ if (.not. allocated(state%rho_dry)) allocate(state%rho_dry(ncrms,crm_nx,crm_ny,crm_nz))
if (.not. allocated(state%qv)) allocate(state%qv(ncrms,crm_nx,crm_ny,crm_nz))
call prefetch(state%u_wind)
call prefetch(state%v_wind)
call prefetch(state%w_wind)
call prefetch(state%temperature)
+ call prefetch(state%rho_dry)
call prefetch(state%qv)
if (trim(MMF_microphysics_scheme) .eq. 'sam1mom') then
@@ -99,6 +109,19 @@ subroutine crm_state_initialize(state,ncrms,crm_nx,crm_ny,crm_nz,MMF_microphysic
call prefetch(state%bm)
call prefetch(state%t_prev)
call prefetch(state%q_prev)
+
+ ! SHOC variables
+ if (.not. allocated(state%shoc_tk )) allocate(state%shoc_tk (ncrms,crm_nx,crm_ny,crm_nz))
+ if (.not. allocated(state%shoc_tkh )) allocate(state%shoc_tkh (ncrms,crm_nx,crm_ny,crm_nz))
+ if (.not. allocated(state%shoc_wthv )) allocate(state%shoc_wthv (ncrms,crm_nx,crm_ny,crm_nz))
+ if (.not. allocated(state%shoc_relvar )) allocate(state%shoc_relvar (ncrms,crm_nx,crm_ny,crm_nz))
+ if (.not. allocated(state%shoc_cldfrac )) allocate(state%shoc_cldfrac (ncrms,crm_nx,crm_ny,crm_nz))
+ call prefetch(state%shoc_tk )
+ call prefetch(state%shoc_tkh )
+ call prefetch(state%shoc_relvar )
+ call prefetch(state%shoc_wthv )
+ call prefetch(state%shoc_cldfrac )
+
end if
end subroutine crm_state_initialize
@@ -111,6 +134,7 @@ subroutine crm_state_finalize(state, MMF_microphysics_scheme)
if (allocated(state%v_wind)) deallocate(state%v_wind)
if (allocated(state%w_wind)) deallocate(state%w_wind)
if (allocated(state%temperature)) deallocate(state%temperature)
+ if (allocated(state%rho_dry)) deallocate(state%rho_dry)
if (allocated(state%qv)) deallocate(state%qv)
if (allocated(state%qp)) deallocate(state%qp)
@@ -127,6 +151,13 @@ subroutine crm_state_finalize(state, MMF_microphysics_scheme)
if (allocated(state%bm)) deallocate(state%bm)
if (allocated(state%t_prev)) deallocate(state%t_prev)
if (allocated(state%q_prev)) deallocate(state%q_prev)
+
+ ! SHOC variables
+ if (allocated(state%shoc_tk )) deallocate(state%shoc_tk )
+ if (allocated(state%shoc_tkh )) deallocate(state%shoc_tkh )
+ if (allocated(state%shoc_wthv )) deallocate(state%shoc_wthv )
+ if (allocated(state%shoc_relvar )) deallocate(state%shoc_relvar )
+ if (allocated(state%shoc_cldfrac )) deallocate(state%shoc_cldfrac )
end subroutine crm_state_finalize
end module crm_state_module
diff --git a/components/eam/src/physics/crm/dummy_modules/debug_info.F90 b/components/eam/src/physics/crm/dummy_modules/debug_info.F90
new file mode 100644
index 000000000000..340871911368
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/debug_info.F90
@@ -0,0 +1,5 @@
+module debug_info
+!-------------------------------------------------------------------------------
+! Dummy module to override src/physics/cam/debug_info.F90
+!-------------------------------------------------------------------------------
+end module debug_info
diff --git a/components/eam/src/physics/crm/dummy_modules/micro_p3_utils.F90 b/components/eam/src/physics/crm/dummy_modules/micro_p3_utils.F90
new file mode 100644
index 000000000000..c7c88a700ff3
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/micro_p3_utils.F90
@@ -0,0 +1,5 @@
+module micro_p3_utils
+!-------------------------------------------------------------------------------
+! Dummy module to override src/physics/cam/micro_p3_utils.F90
+!-------------------------------------------------------------------------------
+end module micro_p3_utils
diff --git a/components/eam/src/physics/crm/dummy_modules/microp_driver.F90 b/components/eam/src/physics/crm/dummy_modules/microp_driver.F90
new file mode 100644
index 000000000000..b5449469744d
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/microp_driver.F90
@@ -0,0 +1,41 @@
+module microp_driver
+!-------------------------------------------------------------------------------
+! Dummy module to avoid building the corresponding module in src/physics/cam
+!-------------------------------------------------------------------------------
+use shr_kind_mod, only: r8 => shr_kind_r8
+implicit none
+private
+save
+public :: microp_driver_readnl
+public :: microp_driver_init_cnst
+public :: microp_driver_implements_cnst
+!===============================================================================
+contains
+!===============================================================================
+subroutine microp_driver_readnl(nlfile)
+ character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
+ ! Read in namelist for microphysics scheme
+ !-----------------------------------------------------------------------
+ return
+end subroutine microp_driver_readnl
+!===============================================================================
+function microp_driver_implements_cnst(name)
+ ! Return true if specified constituent is implemented by the
+ ! microphysics package
+ character(len=*), intent(in) :: name ! constituent name
+ logical :: microp_driver_implements_cnst ! return value
+ !-----------------------------------------------------------------------
+ microp_driver_implements_cnst = .false.
+end function microp_driver_implements_cnst
+!===============================================================================
+subroutine microp_driver_init_cnst(name, q, gcid)
+ ! Initialize the microphysics constituents, if they are
+ ! not read from the initial file.
+ character(len=*), intent(in) :: name ! constituent name
+ real(r8), intent(out) :: q(:,:) ! mass mixing ratio (gcol, plev)
+ integer, intent(in) :: gcid(:) ! global column id
+ !-----------------------------------------------------------------------
+ return
+end subroutine microp_driver_init_cnst
+!===============================================================================
+end module microp_driver
diff --git a/components/eam/src/physics/crm/dummy_modules/scream_abortutils.F90 b/components/eam/src/physics/crm/dummy_modules/scream_abortutils.F90
new file mode 100644
index 000000000000..0408fe2e4c59
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/scream_abortutils.F90
@@ -0,0 +1,5 @@
+module scream_abortutils
+!-------------------------------------------------------------------------------
+! Dummy module to override src/physics/cam/scream_abortutils.F90
+!-------------------------------------------------------------------------------
+end module scream_abortutils
diff --git a/components/eam/src/physics/crm/dummy_modules/shoc.F90 b/components/eam/src/physics/crm/dummy_modules/shoc.F90
new file mode 100644
index 000000000000..4ff4cdf18ba2
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/shoc.F90
@@ -0,0 +1,5 @@
+module shoc
+!-------------------------------------------------------------------------------
+! Dummy module to override src/physics/cam/shoc.F90
+!-------------------------------------------------------------------------------
+end module shoc
diff --git a/components/eam/src/physics/crm/dummy_modules/shoc_intr.F90 b/components/eam/src/physics/crm/dummy_modules/shoc_intr.F90
new file mode 100644
index 000000000000..884e6eb1ee5c
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/shoc_intr.F90
@@ -0,0 +1,13 @@
+module shoc_intr
+!-------------------------------------------------------------------------------
+! Dummy module to override src/physics/cam/shoc_intr.F90
+!-------------------------------------------------------------------------------
+public :: mmf_micro_p3_utils_init
+contains
+!===============================================================================
+subroutine shoc_readnl(nlfile)
+ character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
+ return
+end subroutine shoc_readnl
+!===============================================================================
+end module shoc_intr
diff --git a/components/eam/src/physics/crm/dummy_modules/wv_sat_scream.F90 b/components/eam/src/physics/crm/dummy_modules/wv_sat_scream.F90
new file mode 100644
index 000000000000..eb0a97ce0ab6
--- /dev/null
+++ b/components/eam/src/physics/crm/dummy_modules/wv_sat_scream.F90
@@ -0,0 +1,5 @@
+module wv_sat_scream
+!-------------------------------------------------------------------------------
+! Dummy module to override src/physics/cam/wv_sat_scream.F90
+!-------------------------------------------------------------------------------
+end module wv_sat_scream
diff --git a/components/eam/src/physics/crm/pam/CMakeLists.txt b/components/eam/src/physics/crm/pam/CMakeLists.txt
new file mode 100644
index 000000000000..95075debd4fd
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/CMakeLists.txt
@@ -0,0 +1,132 @@
+
+set(PAM_DRIVER_SRC
+ pam_driver.F90
+ pam_driver.cpp
+ pam_feedback.h
+ pam_radiation.h
+ pam_statistics.h
+ pam_state.h
+ params.F90)
+
+add_library(pam_driver
+ ${PAM_DRIVER_SRC})
+
+if (${USE_CUDA})
+ # PAM will be CUDA-linked with device symbols resolved at library creation
+ set_target_properties(pam_driver
+ PROPERTIES
+ LINKER_LANGUAGE
+ CUDA
+ CUDA_SEPARABLE_COMPILATION OFF
+ CUDA_RESOLVE_DEVICE_SYMBOLS ON
+ )
+endif()
+
+string(FIND "${CAM_CONFIG_OPTS}" "-pam_dycor awfl" USE_AWFL)
+string(FIND "${CAM_CONFIG_OPTS}" "-pam_dycor spam" USE_SPAM)
+
+if (NOT USE_AWFL EQUAL -1)
+ set(PAM_DYCORE awfl)
+endif()
+
+if (NOT USE_SPAM EQUAL -1)
+ set(PAM_DYCORE spam)
+endif()
+
+set(PAM_MICRO p3)
+set(PAM_SGS shoc)
+set(PAM_RAD forced)
+set(PAM_SCREAM_USE_CXX True)
+set(SCREAM_HOME ${CMAKE_CURRENT_SOURCE_DIR}/../../../../../..)
+
+# removing this once the cime cmake_macros aren't using obsolete
+# Kokkos settings like KOKKOS_OPTIONS.
+if (KOKKOS_OPTIONS)
+ unset(KOKKOS_OPTIONS)
+endif()
+
+if (PAM_SCREAM_USE_CXX)
+ message(STATUS "*** building WITH C++ scream interface ***")
+ target_compile_definitions(pam_driver PUBLIC "-DP3_CXX -DSHOC_CXX")
+endif()
+
+if ("${PAM_DYCORE}" STREQUAL "spam")
+ set(PAMC_MODEL extrudedmodel)
+ set(PAMC_HAMIL man) # man = anelastic / mce_rho = moist compressible
+ set(PAMC_THERMO constkappavirpottemp)
+ set(PAMC_IO none)
+endif()
+
+if (NOT PAM_SCREAM_USE_CXX)
+ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+ set(PAM_Fortran_FLAGS -g )
+ else()
+ set(PAM_Fortran_FLAGS -g -ffree-line-length-none -ffixed-line-length-none )
+ endif()
+endif()
+
+
+if (PAM_SCREAM_USE_CXX)
+ set(EKAT_ENABLE_YAML_PARSER OFF CACHE BOOL "" FORCE)
+ add_subdirectory(external/physics/scream_cxx_p3_shoc)
+ add_subdirectory(external/physics/scream_cxx_interfaces)
+ # add these here to allow driver to intialize P3
+ include_directories(${SCREAM_HOME}/components/eamxx/src/physics/p3 )
+ include_directories(${SCREAM_HOME}/components/eamxx/src/physics/p3/impl )
+ include_directories(${SCREAM_HOME}/components/eamxx/src/physics/share )
+ include_directories(${SCREAM_HOME}/components/eamxx/src )
+endif()
+
+add_subdirectory(external/pam_core)
+add_subdirectory(external/physics)
+add_subdirectory(external/dynamics)
+
+target_include_directories(pam_driver PUBLIC
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_CURRENT_BINARY_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/pam_core
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/pam_core/modules
+ ${CMAKE_CURRENT_BINARY_DIR}/external/pam_core/modules
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/pam_core/pam_interface
+ ${CMAKE_CURRENT_BINARY_DIR}/external/pam_core/pam_interface
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/physics/micro/p3
+ ${CMAKE_CURRENT_BINARY_DIR}/external/physics/micro/p3
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/physics/sgs/shoc
+ ${CMAKE_CURRENT_BINARY_DIR}/external/physics/sgs/shoc
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/physics/radiation/forced
+ ${CMAKE_CURRENT_BINARY_DIR}/external/physics/radiation/forced
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/dynamics/${PAM_DYCORE}
+ ${CMAKE_CURRENT_BINARY_DIR}/external/dynamics/${PAM_DYCORE}
+ )
+
+if (PAM_SCREAM_USE_CXX)
+target_include_directories(pam_driver PUBLIC
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/physics/scream_cxx_p3_shoc
+ ${CMAKE_CURRENT_BINARY_DIR}/external/physics/scream_cxx_p3_shoc
+ ${CMAKE_CURRENT_SOURCE_DIR}/external/physics/scream_cxx_interfaces
+ ${CMAKE_CURRENT_BINARY_DIR}/external/physics/scream_cxx_interfaces
+ )
+endif()
+
+include(${YAKL_HOME}/yakl_utils.cmake)
+yakl_process_target(pam_driver)
+
+if (NOT PAM_SCREAM_USE_CXX)
+ set_target_properties(pam_driver PROPERTIES
+ Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules
+ )
+endif()
+
+add_compile_options($<$: -g >)
+
+target_compile_options(pam_driver PUBLIC
+ $<$:${PAM_Fortran_FLAGS}>
+)
+
+if (PAM_SCREAM_USE_CXX)
+ target_link_libraries(pam_driver pam_core physics dynamics pam_scream_cxx_interfaces ekat p3 shoc physics_share scream_share)
+else()
+ target_link_libraries(pam_driver pam_core physics dynamics )
+endif()
+
+
diff --git a/components/eam/src/physics/crm/pam/external b/components/eam/src/physics/crm/pam/external
new file mode 160000
index 000000000000..ce614fcd8d1b
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/external
@@ -0,0 +1 @@
+Subproject commit ce614fcd8d1b38e7e638e55e5fcb220531aca178
diff --git a/components/eam/src/physics/crm/pam/pam_accelerate.h b/components/eam/src/physics/crm/pam/pam_accelerate.h
new file mode 100644
index 000000000000..a48ca628f691
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_accelerate.h
@@ -0,0 +1,234 @@
+#pragma once
+
+#include "pam_coupler.h"
+
+void pam_accelerate_nstop( pam::PamCoupler &coupler, int &nstop) {
+ auto crm_accel_factor = coupler.get_option("crm_accel_factor");
+ if(nstop%static_cast((1+crm_accel_factor)) != 0) {
+ printf("pam_accelerate_nstop: Error: (1+crm_accel_factor) does not divide equally into nstop: %4.4d crm_accel_factor: %6.1f \n",nstop, crm_accel_factor);
+ exit(-1);
+ } else {
+ nstop = nstop / (1 + crm_accel_factor);
+ }
+}
+
+
+inline void pam_accelerate_init( pam::PamCoupler &coupler ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto nens = coupler.get_option("ncrms");
+ auto nz = coupler.get_option("crm_nz");
+ auto ny = coupler.get_option("crm_ny");
+ auto nx = coupler.get_option("crm_nx");
+ auto crm_accel_uv = coupler.get_option("crm_accel_uv");
+ //------------------------------------------------------------------------------------------------
+ dm_device.register_and_allocate("accel_save_t", "saved temperature for MSA", {nz,nens}, {"z","nens"} );
+ dm_device.register_and_allocate("accel_save_r", "saved dry density for MSA", {nz,nens}, {"z","nens"} );
+ dm_device.register_and_allocate("accel_save_q", "saved total water for MSA", {nz,nens}, {"z","nens"} );
+ dm_device.register_and_allocate("accel_save_u", "saved uvel for MSA", {nz,nens}, {"z","nens"} );
+ dm_device.register_and_allocate("accel_save_v", "saved vvel for MSA", {nz,nens}, {"z","nens"} );
+ //------------------------------------------------------------------------------------------------
+}
+
+
+inline void pam_accelerate_diagnose( pam::PamCoupler &coupler ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto nens = coupler.get_option("ncrms");
+ auto nz = coupler.get_option("crm_nz");
+ auto ny = coupler.get_option("crm_ny");
+ auto nx = coupler.get_option("crm_nx");
+ auto crm_accel_uv = coupler.get_option("crm_accel_uv");
+ //------------------------------------------------------------------------------------------------
+ auto temp = dm_device.get("temp" );
+ auto rhod = dm_device.get("density_dry");
+ auto rhov = dm_device.get("water_vapor");
+ auto rhol = dm_device.get("cloud_water");
+ auto rhoi = dm_device.get("ice" );
+ auto uvel = dm_device.get("uvel" );
+ auto vvel = dm_device.get("vvel" );
+ auto accel_save_t = dm_device.get("accel_save_t");
+ auto accel_save_r = dm_device.get("accel_save_r");
+ auto accel_save_q = dm_device.get("accel_save_q");
+ auto accel_save_u = dm_device.get("accel_save_u");
+ auto accel_save_v = dm_device.get("accel_save_v");
+ //------------------------------------------------------------------------------------------------
+ // compute horizontal means needed later for mean-state acceleration
+ parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
+ accel_save_t(k,n) = 0.0;
+ accel_save_r(k,n) = 0.0;
+ accel_save_q(k,n) = 0.0;
+ if (crm_accel_uv) {
+ accel_save_u(k,n) = 0.0;
+ accel_save_v(k,n) = 0.0;
+ }
+ });
+ real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+ yakl::atomicAdd( accel_save_t(k,n), temp(k,j,i,n) * r_nx_ny );
+ yakl::atomicAdd( accel_save_r(k,n), rhod(k,j,i,n) * r_nx_ny );
+ yakl::atomicAdd( accel_save_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny );
+ if (crm_accel_uv) {
+ yakl::atomicAdd( accel_save_u(k,n), uvel(k,j,i,n) * r_nx_ny );
+ yakl::atomicAdd( accel_save_v(k,n), vvel(k,j,i,n) * r_nx_ny );
+ }
+ });
+ //------------------------------------------------------------------------------------------------
+}
+
+
+inline void pam_accelerate( pam::PamCoupler &coupler, int nstep, int &nstop ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ using yakl::ScalarLiveOut;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto nens = coupler.get_option("ncrms");
+ auto nz = coupler.get_option("crm_nz");
+ auto ny = coupler.get_option("crm_ny");
+ auto nx = coupler.get_option("crm_nx");
+ //------------------------------------------------------------------------------------------------
+ auto temp = dm_device.get("temp" );
+ auto rhod = dm_device.get("density_dry");
+ auto rhov = dm_device.get("water_vapor");
+ auto rhol = dm_device.get("cloud_water");
+ auto rhoi = dm_device.get("ice" );
+ auto uvel = dm_device.get("uvel" );
+ auto vvel = dm_device.get("vvel" );
+ auto accel_save_t = dm_device.get("accel_save_t");
+ auto accel_save_r = dm_device.get("accel_save_r");
+ auto accel_save_q = dm_device.get("accel_save_q");
+ auto accel_save_u = dm_device.get("accel_save_u");
+ auto accel_save_v = dm_device.get("accel_save_v");
+ //------------------------------------------------------------------------------------------------
+ bool crm_accel_uv = coupler.get_option("crm_accel_uv");
+ real crm_accel_factor = coupler.get_option("crm_accel_factor");
+ //------------------------------------------------------------------------------------------------
+ real2d hmean_t ("hmean_t", nz,nens);
+ real2d hmean_r ("hmean_r", nz,nens);
+ real2d hmean_q ("hmean_q", nz,nens);
+ real2d hmean_u ("hmean_u", nz,nens);
+ real2d hmean_v ("hmean_v", nz,nens);
+ real2d ttend_acc("ttend_acc", nz,nens);
+ real2d rtend_acc("rtend_acc", nz,nens);
+ real2d qtend_acc("qtend_acc", nz,nens);
+ real2d utend_acc("utend_acc", nz,nens);
+ real2d vtend_acc("vtend_acc", nz,nens);
+ real2d qpoz ("qpoz", nz,nens);
+ real2d qneg ("qneg", nz,nens);
+ //------------------------------------------------------------------------------------------------
+ real constexpr dtemp_max = 5; // temperature tendency max threshold => 5 K following UP-CAM
+ real constexpr temp_min = 50; // temperature minimum minthreshold => 50 K following UP-CAM
+ //------------------------------------------------------------------------------------------------
+ // Compute the horizontal mean for each variable
+ parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
+ hmean_t(k,n) = 0.0;
+ hmean_r(k,n) = 0.0;
+ hmean_q(k,n) = 0.0;
+ if (crm_accel_uv) {
+ hmean_u(k,n) = 0.0;
+ hmean_v(k,n) = 0.0;
+ }
+ });
+ real r_nx_ny = 1._fp/(nx*ny); // precompute reciprocal to avoid costly divisions
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+ yakl::atomicAdd( hmean_t(k,n), temp(k,j,i,n) * r_nx_ny );
+ yakl::atomicAdd( hmean_r(k,n), rhod(k,j,i,n) * r_nx_ny );
+ yakl::atomicAdd( hmean_q(k,n), ( rhov(k,j,i,n) + rhol(k,j,i,n) + rhoi(k,j,i,n) ) * r_nx_ny );
+ if (crm_accel_uv) {
+ yakl::atomicAdd( hmean_u(k,n), uvel(k,j,i,n) * r_nx_ny );
+ yakl::atomicAdd( hmean_v(k,n), vvel(k,j,i,n) * r_nx_ny );
+ }
+ });
+ //------------------------------------------------------------------------------------------------
+ // Compute the accelerated tendencies
+ // NOTE - these are tendencies multiplied by the time step (i.e. d/dt * crm_dt)
+ ScalarLiveOut ceaseflag_liveout(false);
+ parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
+ ttend_acc(k,n) = hmean_t(k,n) - accel_save_t(k,n);
+ rtend_acc(k,n) = hmean_r(k,n) - accel_save_r(k,n);
+ qtend_acc(k,n) = hmean_q(k,n) - accel_save_q(k,n);
+ if (crm_accel_uv) {
+ utend_acc(k,n) = hmean_u(k,n) - accel_save_u(k,n);
+ vtend_acc(k,n) = hmean_v(k,n) - accel_save_v(k,n);
+ }
+ if (abs(ttend_acc(k,n)) > dtemp_max) {
+ ceaseflag_liveout = true;
+ }
+ });
+ bool ceaseflag = ceaseflag_liveout.hostRead();
+ //------------------------------------------------------------------------------------------------
+ // If acceleration tendencies are insane then just abort the acceleration
+ if (ceaseflag) {
+ // When temperature tendency threshold is triggered the acceleration will
+ // not be applied for the remainder of the CRM integration.
+ // The number of CRM steps for this integration (nstop) must be updated
+ // to ensure the CRM integration duration is unchanged.
+ //
+ // The effective MSA timestep is dt_a = crm_dt * (1 + crm_accel_factor). When
+ // ceaseflag is triggered at nstep, we've taken (nstep - 1) previous steps of
+ // size crm_dt * (1 + crm_accel_factor). The current step, and all future
+ // steps, will revert to size crm_dt. Therefore, the total crm integration
+ // time remaining after this step is
+ // time_remaining = crm_run_time - (nstep - 1)* dt_a + crm_dt
+ // nsteps_remaining = time_remaining / crm_dt
+ // updated nstop = nstep + nsteps_remaining
+ // Because we set nstop = crm_run_time / dt_a in crm_accel_nstop, subbing
+ // crm_run_time = nstop * dt_a and working through algebra yields
+ // updated nstop = nstop + (nstop - nstep + 1) * crm_accel_factor.
+ // This only can happen once!
+ coupler.set_option("crm_acceleration_ceaseflag",true);
+ int nstop_old = nstop;
+ int nstop_new = nstop_old + (nstop_old - nstep + 1)*crm_accel_factor;
+ printf("pam_accelerate: MSA tendencies too large - nstop increased from %d to %d \n",nstop_old, nstop_new);
+ nstop = nstop_new;
+ return;
+ }
+ //------------------------------------------------------------------------------------------------
+ // Apply the accelerated tendencies
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+
+ temp(k,j,i,n) = temp(k,j,i,n) + crm_accel_factor * ttend_acc(k,n);
+ rhod(k,j,i,n) = rhod(k,j,i,n) + crm_accel_factor * rtend_acc(k,n);
+ rhov(k,j,i,n) = rhov(k,j,i,n) + crm_accel_factor * qtend_acc(k,n);
+ if (crm_accel_uv) {
+ uvel(k,j,i,n) = uvel(k,j,i,n) + crm_accel_factor * utend_acc(k,n);
+ vvel(k,j,i,n) = vvel(k,j,i,n) + crm_accel_factor * vtend_acc(k,n);
+ }
+ // apply limiter on temperature
+ temp(k,j,i,n) = std::max( temp_min, temp(k,j,i,n) );
+ });
+ //------------------------------------------------------------------------------------------------
+ // Evaluate and address instances of negative water vapor
+ parallel_for( SimpleBounds<2>(nz,nens) , YAKL_LAMBDA (int k, int n) {
+ qpoz(k,n) = 0.0;
+ qneg(k,n) = 0.0;
+ });
+ // accumulate positive and negative water vapor density values in each layer
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+ if (rhov(k,j,i,n) < 0.0) {
+ yakl::atomicAdd( qneg(k,n) , rhov(k,j,i,n) );
+ } else {
+ yakl::atomicAdd( qpoz(k,n) , rhov(k,j,i,n) );
+ }
+ });
+ // clip or adjust points with negative water vapor density
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+ real factor;
+ // check if all moisture is depleted in the layer
+ if (qpoz(k,n) + qneg(k,n) <= 0.0) {
+ rhov(k,j,i,n) = 0.0;
+ } else {
+ // Clip vapor density values at 0 and remove the negative excess in each layer
+ // proportionally from the positive qt fields in the layer
+ factor = 1.0 + qneg(k,n) / qpoz(k,n);
+ rhov(k,j,i,n) = std::max(0.0, rhov(k,j,i,n) * factor);
+ }
+ });
+ //------------------------------------------------------------------------------------------------
+}
+
diff --git a/components/eam/src/physics/crm/pam/pam_debug.h b/components/eam/src/physics/crm/pam/pam_debug.h
new file mode 100644
index 000000000000..5186d51c6a9d
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_debug.h
@@ -0,0 +1,211 @@
+#pragma once
+
+#include "pam_coupler.h"
+
+// These routines were helpful for debugging the coupling
+// between PAM and E3SM, so we kept them here for future use
+
+// register and initialize various quantities for statistical calculations
+inline void pam_debug_init( pam::PamCoupler &coupler ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto nx = coupler.get_option("crm_nx");
+ auto ny = coupler.get_option("crm_ny");
+ auto nz = coupler.get_option("crm_nz");
+ auto nens = coupler.get_option("ncrms");
+ //------------------------------------------------------------------------------------------------
+ dm_device.register_and_allocate("debug_save_temp", "saved temp for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} );
+ dm_device.register_and_allocate("debug_save_rhod", "saved rhod for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} );
+ dm_device.register_and_allocate("debug_save_rhov", "saved rhov for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} );
+ dm_device.register_and_allocate("debug_save_rhoc", "saved rhoc for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} );
+ dm_device.register_and_allocate("debug_save_rhoi", "saved rhoi for debug", {nz,ny,nx,nens}, {"z","y","x","nens"} );
+ auto debug_save_temp = dm_device.get("debug_save_temp");
+ auto debug_save_rhod = dm_device.get("debug_save_rhod");
+ auto debug_save_rhov = dm_device.get("debug_save_rhov");
+ auto debug_save_rhoc = dm_device.get("debug_save_rhoc");
+ auto debug_save_rhoi = dm_device.get("debug_save_rhoi");
+ //------------------------------------------------------------------------------------------------
+ auto temp = dm_device.get("temp");
+ auto rhod = dm_device.get("density_dry");
+ auto rhov = dm_device.get("water_vapor");
+ auto rhoc = dm_device.get("cloud_water");
+ auto rhoi = dm_device.get("ice");
+ //------------------------------------------------------------------------------------------------
+ parallel_for("copy data to saved debug variables", SimpleBounds<4>(nz,ny,nx,nens),
+ YAKL_LAMBDA (int k, int j, int i, int iens) {
+ debug_save_temp(k,j,i,iens) = temp(k,j,i,iens);
+ debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens);
+ debug_save_rhov(k,j,i,iens) = rhov(k,j,i,iens);
+ debug_save_rhoc(k,j,i,iens) = rhoc(k,j,i,iens);
+ debug_save_rhoi(k,j,i,iens) = rhoi(k,j,i,iens);
+ });
+ //------------------------------------------------------------------------------------------------
+}
+
+// main routine to encapsulate checking for bad values across many variables
+void pam_debug_check_state( pam::PamCoupler &coupler, int id, int nstep ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto &dm_host = coupler.get_data_manager_host_readwrite();
+ auto nx = coupler.get_option("crm_nx");
+ auto ny = coupler.get_option("crm_ny");
+ auto nz = coupler.get_option("crm_nz");
+ auto nens = coupler.get_option("ncrms");
+ auto temp = dm_device.get("temp");
+ auto rhod = dm_device.get("density_dry");
+ auto rhov = dm_device.get("water_vapor");
+ auto rhoc = dm_device.get("cloud_water");
+ auto rhoi = dm_device.get("ice");
+ auto debug_save_temp = dm_device.get("debug_save_temp");
+ auto debug_save_rhod = dm_device.get("debug_save_rhod");
+ auto debug_save_rhov = dm_device.get("debug_save_rhov");
+ auto debug_save_rhoc = dm_device.get("debug_save_rhoc");
+ auto debug_save_rhoi = dm_device.get("debug_save_rhoi");
+ real grav = 9.80616;
+ auto input_phis = dm_host.get("input_phis").createDeviceCopy();
+ auto lat = dm_host.get("latitude" ).createDeviceCopy();
+ auto lon = dm_host.get("longitude" ).createDeviceCopy();
+ //------------------------------------------------------------------------------------------------
+ // Check for invalid values
+ parallel_for("", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) {
+ // Check for NaNs
+ const auto is_nan_t_atm = isnan( temp(k,j,i,iens) );
+ const auto is_nan_r_atm = isnan( rhod(k,j,i,iens) );
+ const auto is_nan_q_atm = isnan( rhov(k,j,i,iens) );
+ if ( is_nan_t_atm || is_nan_r_atm || is_nan_q_atm ) {
+ auto phis = input_phis(iens)/grav;
+ printf("PAM-DEBUG nan-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g \n",
+ nstep,id,k,i,iens,lat(iens),lon(iens),phis,
+ temp(k,j,i,iens),
+ rhod(k,j,i,iens),
+ rhov(k,j,i,iens),
+ rhoc(k,j,i,iens),
+ rhoi(k,j,i,iens),
+ debug_save_temp(k,j,i,iens),
+ debug_save_rhod(k,j,i,iens),
+ debug_save_rhov(k,j,i,iens),
+ debug_save_rhoc(k,j,i,iens),
+ debug_save_rhoi(k,j,i,iens)
+ );
+ }
+ // Check for negative values
+ const auto is_neg_t_atm = temp(k,j,i,iens)<0;
+ const auto is_neg_r_atm = rhod(k,j,i,iens)<0;
+ const auto is_neg_q_atm = rhov(k,j,i,iens)<0;
+ if ( is_neg_t_atm || is_neg_r_atm || is_neg_q_atm ) {
+ auto phis = input_phis(iens)/grav;
+ printf("PAM-DEBUG neg-found - st:%3.3d id:%2.2d k:%3.3d i:%3.3d n:%3.3d y:%5.1f x:%5.1f ph:%6.1f -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g -- t:%8.2g rd:%8.2g rv:%8.2g rc:%8.2g ri:%8.2g \n",
+ nstep,id,k,i,iens,lat(iens),lon(iens),phis,
+ temp(k,j,i,iens),
+ rhod(k,j,i,iens),
+ rhov(k,j,i,iens),
+ rhoc(k,j,i,iens),
+ rhoi(k,j,i,iens),
+ debug_save_temp(k,j,i,iens),
+ debug_save_rhod(k,j,i,iens),
+ debug_save_rhov(k,j,i,iens),
+ debug_save_rhoc(k,j,i,iens),
+ debug_save_rhoi(k,j,i,iens)
+ );
+ }
+ // update saved previous values
+ debug_save_temp(k,j,i,iens) = temp(k,j,i,iens);
+ debug_save_rhod(k,j,i,iens) = rhod(k,j,i,iens);
+ debug_save_rhov(k,j,i,iens) = rhov(k,j,i,iens);
+ debug_save_rhoc(k,j,i,iens) = rhoc(k,j,i,iens);
+ debug_save_rhoi(k,j,i,iens) = rhoi(k,j,i,iens);
+ });
+ //------------------------------------------------------------------------------------------------
+}
+
+// print the min and max of PAM state variables to help look for problems
+// void print_state_min_max( pam::PamCoupler &coupler, std::string id ) {
+ // auto &dm_device = coupler.get_data_manager_device_readwrite();
+ // auto nz = coupler.get_option("crm_nz");
+ // auto nx = coupler.get_option("crm_nx");
+ // auto ny = coupler.get_option("crm_ny");
+ // auto nens = coupler.get_option("ncrms");
+ // // auto gcm_nlev = coupler.get_option("gcm_nlev");
+ // auto pmid = coupler.compute_pressure_array();
+ // auto zmid = dm_device.get("vertical_midpoint_height" );
+ // // auto gcm_rho_d = dm_device.get("gcm_density_dry");
+ // auto temp = dm_device.get("temp");
+ // auto crm_rho_v = dm_device.get("water_vapor");
+ // auto crm_rho_c = dm_device.get("cloud_water");
+ // auto crm_rho_d = dm_device.get("density_dry");
+ // for (int k=0; k("crm_nz");
+ auto nx = coupler.get_option("crm_nx");
+ auto ny = coupler.get_option("crm_ny");
+ auto nens = coupler.get_option("ncrms");
+ // auto pmid = coupler.compute_pressure_array();
+ // auto zmid = dm_device.get("vertical_midpoint_height" );
+ auto temp = dm_device.get("temp");
+ auto rho_v = dm_device.get("water_vapor");
+ auto rho_c = dm_device.get("cloud_water");
+ auto rho_d = dm_device.get("density_dry");
+ auto rho_i = dm_device.get("ice");
+ // parallel_for("", SimpleBounds<4>(nz,ny,nx,nens), YAKL_LAMBDA (int k, int j, int i, int iens) {
+ parallel_for("pam_debug_print_state", SimpleBounds<2>(nz,nx), YAKL_LAMBDA (int k, int i) {
+ int j = 0;
+ int n = 0;
+ printf("PAM-DEBUG %d - k:%d i:%d temp : %g rv: %g rc: %g ri: %g \n",
+ id,k,i,
+ temp(k,j,i,n),
+ rho_v(k,j,i,n),
+ rho_c(k,j,i,n),
+ rho_i(k,j,i,n)
+ );
+ });
+}
diff --git a/components/eam/src/physics/crm/pam/pam_driver.F90 b/components/eam/src/physics/crm/pam/pam_driver.F90
new file mode 100644
index 000000000000..ca6ffb49d1d6
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_driver.F90
@@ -0,0 +1,10 @@
+module pam_driver_mod
+ use iso_c_binding
+ implicit none
+ interface
+ subroutine pam_driver() bind(C,name="pam_driver")
+ end subroutine
+ subroutine pam_finalize() bind(C,name="pam_finalize")
+ end subroutine
+ end interface
+end module pam_driver_mod
diff --git a/components/eam/src/physics/crm/pam/pam_driver.cpp b/components/eam/src/physics/crm/pam/pam_driver.cpp
new file mode 100644
index 000000000000..f131bbd44f09
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_driver.cpp
@@ -0,0 +1,256 @@
+#include "pam_coupler.h"
+// #include "params.h"
+#include "Dycore.h"
+#include "Microphysics.h"
+#include "SGS.h"
+#include "radiation.h"
+#include "pam_interface.h"
+#include "perturb_temperature.h"
+#include "gcm_forcing.h"
+#include "pam_feedback.h"
+#include "pam_state.h"
+#include "pam_radiation.h"
+#include "pam_statistics.h"
+#include "pam_output.h"
+#include "pam_accelerate.h"
+#include "sponge_layer.h"
+#include "surface_friction.h"
+#include "scream_cxx_interface_finalize.h"
+
+#include "pam_hyperdiffusion.h"
+
+// Needed for p3_init
+#include "p3_functions.hpp"
+#include "p3_f90.hpp"
+
+#include "pam_debug.h"
+bool constexpr enable_check_state = false;
+
+extern "C" void pam_driver() {
+ //------------------------------------------------------------------------------------------------
+ using yakl::intrinsics::abs;
+ using yakl::intrinsics::maxval;
+ using yakl::atomicAdd;
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ auto &coupler = pam_interface::get_coupler();
+ //------------------------------------------------------------------------------------------------
+ // retreive coupler options
+ auto nens = coupler.get_option("ncrms");
+ auto gcm_nlev = coupler.get_option("gcm_nlev");
+ auto crm_nz = coupler.get_option("crm_nz");
+ auto crm_nx = coupler.get_option("crm_nx");
+ auto crm_ny = coupler.get_option("crm_ny");
+ auto gcm_dt = coupler.get_option("gcm_dt");
+ auto crm_dt = coupler.get_option("crm_dt");
+ auto is_first_step = coupler.get_option("is_first_step");
+ auto is_restart = coupler.get_option("is_restart");
+ bool use_crm_accel = coupler.get_option("use_crm_accel");
+ bool enable_physics_tend_stats = coupler.get_option("enable_physics_tend_stats");
+ //------------------------------------------------------------------------------------------------
+ // set various coupler options
+ coupler.set_option("gcm_physics_dt",gcm_dt);
+ #ifdef MMF_PAM_DPP
+ // this is leftover from debugging, but it might still be useful for testing values of crm_per_phys
+ coupler.set_option("crm_per_phys",MMF_PAM_DPP);
+ #else
+ coupler.set_option("crm_per_phys",2); // # of PAM-C dynamics steps per physics
+ #endif
+ coupler.set_option("sponge_num_layers",crm_nz*0.3); // depth of sponge layer
+ coupler.set_option("sponge_time_scale",60); // minimum damping timescale at top
+ coupler.set_option("crm_acceleration_ceaseflag",false);
+ //------------------------------------------------------------------------------------------------
+ // Allocate the coupler state and retrieve host/device data managers
+ coupler.allocate_coupler_state( crm_nz , crm_ny , crm_nx , nens );
+
+ // set up the grid - this needs to happen before initializing coupler objects
+ pam_state_set_grid(coupler);
+ //------------------------------------------------------------------------------------------------
+ // get seperate data manager objects for host and device
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto &dm_host = coupler.get_data_manager_host_readwrite();
+ //------------------------------------------------------------------------------------------------
+ // Create objects for dycor, microphysics, and turbulence and initialize them
+ bool verbose = is_first_step || is_restart;
+ Microphysics micro;
+ SGS sgs;
+ Dycore dycore;
+ Radiation rad;
+ micro .init(coupler);
+ sgs .init(coupler);
+ dycore.init(coupler,verbose); // pass is_first_step to control verbosity in PAM-C
+ rad .init(coupler);
+ //------------------------------------------------------------------------------------------------
+ // update coupler GCM state with input GCM state
+ pam_state_update_gcm_state(coupler);
+
+ // Copy input CRM state (saved by the GCM) to coupler
+ pam_state_copy_input_to_coupler(coupler);
+
+ // // update CRM dry density to match GCM and disable dry density forcing
+ // pam_state_update_dry_density(coupler);
+
+ // if debugging - initialize saved state variables and check initial CRM state
+ if (enable_check_state) {
+ pam_debug_init(coupler);
+ pam_debug_check_state(coupler, 0, 0);
+ }
+
+ // Compute CRM forcing tendencies
+ modules::compute_gcm_forcing_tendencies(coupler);
+
+ // Copy input radiation tendencies to coupler
+ pam_radiation_copy_input_to_coupler(coupler);
+
+ // initialize aggregated variables needed for radiation
+ pam_radiation_init(coupler);
+
+ // initialize aggregated variables for output statistics
+ pam_statistics_init(coupler);
+
+ // initialize variables for CRM mean-state acceleration
+ if (use_crm_accel) { pam_accelerate_init(coupler); }
+
+ // initilize surface "psuedo-friction" (psuedo => doesn't match "real" GCM friction)
+ auto input_tau = dm_host.get("input_tau00").createDeviceCopy();
+ auto input_bflx = dm_host.get("input_bflxls").createDeviceCopy();
+ modules::surface_friction_init(coupler, input_tau, input_bflx);
+
+ // Perturb the CRM at the only on first CRM call
+ if (is_first_step) {
+ auto global_column_id = dm_host.get("global_column_id").createDeviceCopy();
+ modules::perturb_temperature( coupler , global_column_id );
+ }
+
+ // Microphysics initialization - load lookup tables
+ #if defined(P3_CXX)
+ if (is_first_step || is_restart) {
+ auto am_i_root = coupler.get_option("am_i_root");
+ scream::p3::p3_init(/*write_tables=*/false, am_i_root);
+ pam::p3_init_lookup_tables(); // Load P3 lookup table data - avoid re-loading every CRM call
+ }
+ #endif
+
+ // dycor initialization
+ bool do_density_save_recall = false;
+ #if defined(MMF_PAM_DYCOR_SPAM)
+ // The PAM-C anelastic dycor can dramatically change the dry density due to the assumption that
+ // the total density does not change, so we will save it here and recall after the dycor
+ do_density_save_recall = true;
+ pam_state_set_reference_state(coupler);
+ dycore.pre_time_loop(coupler);
+ #elif defined(MMF_PAM_DYCOR_AWFL)
+ dycore.declare_current_profile_as_hydrostatic(coupler,/*use_gcm_data=*/true);
+ #endif
+
+ //------------------------------------------------------------------------------------------------
+ //------------------------------------------------------------------------------------------------
+ //------------------------------------------------------------------------------------------------
+
+ // set number of CRM steps
+ int nstop = int(gcm_dt/crm_dt);
+
+ // for mean-state acceleration adjust nstop and diagnose horizontal means
+ if (use_crm_accel) {
+ pam_accelerate_nstop(coupler,nstop);
+ pam_accelerate_diagnose(coupler);
+ };
+
+ // Run the CRM
+ real etime_crm = 0;
+ int nstep = 0;
+ // while (etime_crm < gcm_dt) {
+ while (nstep < nstop) {
+ if (crm_dt == 0.) { crm_dt = dycore.compute_time_step(coupler); }
+ if (etime_crm + crm_dt > gcm_dt) { crm_dt = gcm_dt - etime_crm; }
+
+ if (enable_check_state) { pam_debug_check_state(coupler, 1, nstep); }
+
+ // run a PAM time step
+ coupler.run_module( "apply_gcm_forcing_tendencies" , modules::apply_gcm_forcing_tendencies );
+ coupler.run_module( "radiation" , [&] (pam::PamCoupler &coupler) {rad .timeStep(coupler);} );
+ if (enable_check_state) { pam_debug_check_state(coupler, 2, nstep); }
+
+ // Dynamics
+ if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); }
+ if (do_density_save_recall) { pam_state_save_dry_density(coupler); }
+ coupler.run_module( "dycore", [&] (pam::PamCoupler &coupler) {dycore.timeStep(coupler);} );
+ if (do_density_save_recall) { pam_state_recall_dry_density(coupler); }
+ if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"dycor"); }
+ if (enable_check_state) { pam_debug_check_state(coupler, 3, nstep); }
+
+ // Sponge layer damping
+ if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); }
+ coupler.run_module( "sponge_layer", modules::sponge_layer );
+ if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sponge"); }
+ if (enable_check_state) { pam_debug_check_state(coupler, 4, nstep); }
+
+ // Apply hyperdiffusion to account for lack of horizontal mixing in SHOC
+ pam_hyperdiffusion(coupler);
+
+ // Turbulence - SHOC
+ coupler.run_module( "compute_surface_friction", modules::compute_surface_friction );
+ if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); }
+ coupler.run_module( "sgs", [&] (pam::PamCoupler &coupler) {sgs .timeStep(coupler);} );
+ if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"sgs"); }
+ if (enable_check_state) { pam_debug_check_state(coupler, 5, nstep); }
+
+ // Microphysics - P3
+ if (enable_physics_tend_stats) { pam_statistics_save_state(coupler); }
+ coupler.run_module( "micro", [&] (pam::PamCoupler &coupler) {micro .timeStep(coupler);} );
+ if (enable_physics_tend_stats) { pam_statistics_aggregate_tendency(coupler,"micro"); }
+ if (enable_check_state) { pam_debug_check_state(coupler, 6, nstep); }
+
+ // CRM mean state acceleration
+ if (use_crm_accel && !coupler.get_option("crm_acceleration_ceaseflag")) {
+ pam_accelerate(coupler, nstep, nstop);
+ pam_accelerate_diagnose(coupler);
+ }
+
+ // Diagnostic aggregation
+ pam_radiation_timestep_aggregation(coupler);
+ pam_statistics_timestep_aggregation(coupler);
+
+ etime_crm += crm_dt;
+ nstep += 1;
+ }
+ //------------------------------------------------------------------------------------------------
+ //------------------------------------------------------------------------------------------------
+ //------------------------------------------------------------------------------------------------
+
+ // Compute CRM feedback tendencies and copy to host
+ pam_feedback_compute_tendencies(coupler,gcm_dt);
+ pam_feedback_copy_to_host(coupler);
+
+ // Copy the final CRM state to the host to be saved for next time step
+ pam_state_copy_to_host(coupler);
+
+ // Compute horizontal means of CRM state variables and copy to host
+ pam_output_compute_means(coupler);
+ pam_output_copy_to_host(coupler);
+
+ // convert aggregated radiation quantities to means and copy to host
+ pam_radiation_compute_means(coupler);
+ pam_radiation_copy_to_host(coupler);
+
+ // convert aggregated diagnostic quantities to means and copy to host
+ pam_statistics_compute_means(coupler);
+ pam_statistics_copy_to_host(coupler);
+
+ //------------------------------------------------------------------------------------------------
+ // Finalize and clean up
+ micro .finalize(coupler);
+ sgs .finalize(coupler);
+ dycore.finalize(coupler);
+ rad .finalize(coupler);
+ pam_interface::finalize();
+ //------------------------------------------------------------------------------------------------
+}
+
+extern "C" void pam_finalize() {
+ #if defined(P3_CXX) || defined(SHOC_CXX)
+ pam::deallocate_scream_cxx_globals();
+ // if using SL tracer advection then COMPOSE will call Kokkos::finalize(), otherwise, call it here
+ // pam::call_kokkos_finalize();
+ #endif
+}
diff --git a/components/eam/src/physics/crm/pam/pam_feedback.h b/components/eam/src/physics/crm/pam/pam_feedback.h
new file mode 100644
index 000000000000..3c45153f79fd
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_feedback.h
@@ -0,0 +1,152 @@
+#pragma once
+
+#include "pam_coupler.h"
+
+// These routines are only called once at the end of the CRM call
+// to provide the tendencies and fields to couple the CRM and GCM
+
+// Terminology:
+// - GCM state at the beginning of the CRM call: state_gcm
+// - instantaneous horizontally-averaged CRM state: state_crm
+// The GCM forces the CRM in the MMF by computing
+// (state_gcm - state_crm) / gcm_dt
+// The forcing is applied to the CRM state at each CRM physics time step.
+// Similarly, the CRM provides a feedback tendency to the GCM by computing
+// (state_crm - state_gcm) / gcm_dt
+// using the final state of the CRM at the end of the integration
+
+
+// Compute feedback tendencies for the GCM
+inline void pam_feedback_compute_tendencies( pam::PamCoupler &coupler , real gcm_dt ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto &dm_host = coupler.get_data_manager_host_readwrite();
+ int crm_nz = dm_device.get_dimension_size("z" );
+ int crm_ny = dm_device.get_dimension_size("y" );
+ int crm_nx = dm_device.get_dimension_size("x" );
+ int nens = dm_device.get_dimension_size("nens");
+ int gcm_nlev = coupler.get_option("gcm_nlev");
+ //------------------------------------------------------------------------------------------------
+ // Get current CRM state
+ auto rho_d = dm_device.get("density_dry");
+ auto uvel = dm_device.get("uvel" );
+ auto vvel = dm_device.get("vvel" );
+ auto temp = dm_device.get("temp" );
+ auto rho_v = dm_device.get("water_vapor");
+ auto rho_c = dm_device.get("cloud_water");
+ auto rho_i = dm_device.get("ice" );
+ //------------------------------------------------------------------------------------------------
+ // Get input GCM state
+ auto gcm_ul = dm_host.get("input_ul").createDeviceCopy();
+ auto gcm_vl = dm_host.get("input_vl").createDeviceCopy();
+ auto gcm_tl = dm_host.get("input_tl").createDeviceCopy();
+ auto gcm_qv = dm_host.get("input_ql").createDeviceCopy();
+ auto gcm_qc = dm_host.get("input_qccl").createDeviceCopy();
+ auto gcm_qi = dm_host.get("input_qiil").createDeviceCopy();
+ //------------------------------------------------------------------------------------------------
+ // Create arrays to hold the current column average of the CRM internal columns
+ real2d crm_hmean_uvel ("crm_hmean_uvel" ,crm_nz,nens);
+ real2d crm_hmean_vvel ("crm_hmean_vvel" ,crm_nz,nens);
+ real2d crm_hmean_temp ("crm_hmean_temp" ,crm_nz,nens);
+ real2d crm_hmean_qv ("crm_hmean_qv" ,crm_nz,nens);
+ real2d crm_hmean_qc ("crm_hmean_qc" ,crm_nz,nens);
+ real2d crm_hmean_qi ("crm_hmean_qi" ,crm_nz,nens);
+ //------------------------------------------------------------------------------------------------
+ // We will be essentially reducing a summation to these variables, so initialize them to zero
+ parallel_for("Initialize horzontal means", SimpleBounds<2>(crm_nz,nens), YAKL_LAMBDA (int k_crm, int iens) {
+ crm_hmean_uvel (k_crm,iens) = 0;
+ crm_hmean_vvel (k_crm,iens) = 0;
+ crm_hmean_temp (k_crm,iens) = 0;
+ crm_hmean_qv (k_crm,iens) = 0;
+ crm_hmean_qc (k_crm,iens) = 0;
+ crm_hmean_qi (k_crm,iens) = 0;
+ });
+ //------------------------------------------------------------------------------------------------
+ // Compute horizontal means
+ real r_nx_ny = 1._fp/(crm_nx*crm_ny); // precompute reciprocal to avoid costly divisions
+ parallel_for("Horz mean of CRM state", SimpleBounds<4>(crm_nz,crm_ny,crm_nx,nens), YAKL_LAMBDA (int k_crm, int j, int i, int iens) {
+ real rho_total = rho_d(k_crm,j,i,iens) + rho_v(k_crm,j,i,iens);
+ atomicAdd( crm_hmean_uvel(k_crm,iens), uvel (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( crm_hmean_vvel(k_crm,iens), vvel (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( crm_hmean_temp(k_crm,iens), temp (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( crm_hmean_qv (k_crm,iens),(rho_v(k_crm,j,i,iens)/rho_total) * r_nx_ny );
+ atomicAdd( crm_hmean_qc (k_crm,iens),(rho_c(k_crm,j,i,iens)/rho_total) * r_nx_ny );
+ atomicAdd( crm_hmean_qi (k_crm,iens),(rho_i(k_crm,j,i,iens)/rho_total) * r_nx_ny );
+ });
+ //------------------------------------------------------------------------------------------------
+ // Create arrays to hold the feedback tendencies
+ dm_device.register_and_allocate("crm_feedback_tend_uvel", "feedback tendency of uvel", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("crm_feedback_tend_vvel", "feedback tendency of vvel", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("crm_feedback_tend_dse" , "feedback tendency of dse", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("crm_feedback_tend_qv" , "feedback tendency of qv", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("crm_feedback_tend_qc" , "feedback tendency of qc", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("crm_feedback_tend_qi" , "feedback tendency of qi", {gcm_nlev,nens},{"gcm_lev","nens"});
+ auto crm_feedback_tend_uvel = dm_device.get("crm_feedback_tend_uvel");
+ auto crm_feedback_tend_vvel = dm_device.get("crm_feedback_tend_vvel");
+ auto crm_feedback_tend_dse = dm_device.get("crm_feedback_tend_dse");
+ auto crm_feedback_tend_qv = dm_device.get("crm_feedback_tend_qv");
+ auto crm_feedback_tend_qc = dm_device.get("crm_feedback_tend_qc");
+ auto crm_feedback_tend_qi = dm_device.get("crm_feedback_tend_qi");
+ //------------------------------------------------------------------------------------------------
+ // Compute feedback tendencies
+ real cp_d = coupler.get_option("cp_d");
+ real r_gcm_dt = 1._fp / gcm_dt; // precompute reciprocal to avoid costly divisions
+ parallel_for( "Compute CRM feedback tendencies", SimpleBounds<2>(gcm_nlev,nens), YAKL_LAMBDA (int k_gcm, int iens) {
+ int k_crm = gcm_nlev-1-k_gcm;
+ // if (k_crm("crm_feedback_tend_uvel");
+ auto crm_feedback_tend_vvel = dm_device.get("crm_feedback_tend_vvel");
+ auto crm_feedback_tend_dse = dm_device.get("crm_feedback_tend_dse");
+ auto crm_feedback_tend_qv = dm_device.get("crm_feedback_tend_qv");
+ auto crm_feedback_tend_qc = dm_device.get("crm_feedback_tend_qc");
+ auto crm_feedback_tend_qi = dm_device.get("crm_feedback_tend_qi");
+ //------------------------------------------------------------------------------------------------
+ auto output_ultend_host = dm_host.get("output_ultend");
+ auto output_vltend_host = dm_host.get("output_vltend");
+ auto output_sltend_host = dm_host.get("output_sltend");
+ auto output_qvltend_host = dm_host.get("output_qltend");
+ auto output_qcltend_host = dm_host.get("output_qcltend");
+ auto output_qiltend_host = dm_host.get("output_qiltend");
+ //------------------------------------------------------------------------------------------------
+ // Copy the data to host
+ crm_feedback_tend_uvel.deep_copy_to(output_ultend_host);
+ crm_feedback_tend_vvel.deep_copy_to(output_vltend_host);
+ crm_feedback_tend_dse .deep_copy_to(output_sltend_host);
+ crm_feedback_tend_qv .deep_copy_to(output_qvltend_host);
+ crm_feedback_tend_qc .deep_copy_to(output_qcltend_host);
+ crm_feedback_tend_qi .deep_copy_to(output_qiltend_host);
+ yakl::fence();
+ //------------------------------------------------------------------------------------------------
+}
+
+
diff --git a/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h
new file mode 100644
index 000000000000..83dc83252c68
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_hyperdiffusion.h
@@ -0,0 +1,78 @@
+#pragma once
+#include "pam_coupler.h"
+
+inline void pam_hyperdiffusion( pam::PamCoupler &coupler ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ using yakl::ScalarLiveOut;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto nens = coupler.get_option("ncrms");
+ auto nz = coupler.get_option("crm_nz");
+ auto ny = coupler.get_option("crm_ny");
+ auto nx = coupler.get_option("crm_nx");
+ auto dx = coupler.get_option("crm_dx");
+ auto dy = coupler.get_option("crm_dy");
+ auto dt = coupler.get_option("crm_dt");
+ //------------------------------------------------------------------------------------------------
+ auto temp = dm_device.get("temp" );
+ auto rhod = dm_device.get("density_dry" );
+ auto rhov = dm_device.get("water_vapor" );
+ auto rhol = dm_device.get("cloud_water" );
+ auto rhoi = dm_device.get("ice" );
+ auto nliq = dm_device.get("cloud_water_num");
+ auto nice = dm_device.get("ice_num" );
+ auto uvel = dm_device.get("uvel" );
+ auto vvel = dm_device.get("vvel" );
+ //------------------------------------------------------------------------------------------------
+ real constexpr hd_timescale = 10.0; // damping time scale [sec]
+ //------------------------------------------------------------------------------------------------
+ real4d hd_temp("hd_temp",nz,ny,nx,nens);
+ real4d hd_rhod("hd_rhod",nz,ny,nx,nens);
+ real4d hd_rhov("hd_rhov",nz,ny,nx,nens);
+ real4d hd_rhol("hd_rhol",nz,ny,nx,nens);
+ real4d hd_rhoi("hd_rhoi",nz,ny,nx,nens);
+ real4d hd_nliq("hd_nliq",nz,ny,nx,nens);
+ real4d hd_nice("hd_nice",nz,ny,nx,nens);
+ real4d hd_uvel("hd_uvel",nz,ny,nx,nens);
+ real4d hd_vvel("hd_vvel",nz,ny,nx,nens);
+ //------------------------------------------------------------------------------------------------
+ if (ny==1) {
+ // Calculate hyperdiffusion to all fields in 2D domain configuration
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+ int im2 = i-2;
+ int im1 = i-1;
+ int ip1 = i+1;
+ int ip2 = i+2;
+ if (im2<0 ) { im2 += nx; }
+ if (im1<0 ) { im1 += nx; }
+ if (ip1>nx-1) { ip1 -= nx; }
+ if (ip2>nx-1) { ip2 -= nx; }
+ hd_temp(k,j,i,n) = temp(k,j,im2,n) - 4*temp(k,j,im1,n) + 6*temp(k,j,i,n) - 4*temp(k,j,ip1,n) + temp(k,j,ip2,n);
+ hd_rhod(k,j,i,n) = rhod(k,j,im2,n) - 4*rhod(k,j,im1,n) + 6*rhod(k,j,i,n) - 4*rhod(k,j,ip1,n) + rhod(k,j,ip2,n);
+ hd_rhov(k,j,i,n) = rhov(k,j,im2,n) - 4*rhov(k,j,im1,n) + 6*rhov(k,j,i,n) - 4*rhov(k,j,ip1,n) + rhov(k,j,ip2,n);
+ hd_rhol(k,j,i,n) = rhol(k,j,im2,n) - 4*rhol(k,j,im1,n) + 6*rhol(k,j,i,n) - 4*rhol(k,j,ip1,n) + rhol(k,j,ip2,n);
+ hd_rhoi(k,j,i,n) = rhoi(k,j,im2,n) - 4*rhoi(k,j,im1,n) + 6*rhoi(k,j,i,n) - 4*rhoi(k,j,ip1,n) + rhoi(k,j,ip2,n);
+ hd_nliq(k,j,i,n) = nliq(k,j,im2,n) - 4*nliq(k,j,im1,n) + 6*nliq(k,j,i,n) - 4*nliq(k,j,ip1,n) + nliq(k,j,ip2,n);
+ hd_nice(k,j,i,n) = nice(k,j,im2,n) - 4*nice(k,j,im1,n) + 6*nice(k,j,i,n) - 4*nice(k,j,ip1,n) + nice(k,j,ip2,n);
+ hd_uvel(k,j,i,n) = uvel(k,j,im2,n) - 4*uvel(k,j,im1,n) + 6*uvel(k,j,i,n) - 4*uvel(k,j,ip1,n) + uvel(k,j,ip2,n);
+ hd_vvel(k,j,i,n) = vvel(k,j,im2,n) - 4*vvel(k,j,im1,n) + 6*vvel(k,j,i,n) - 4*vvel(k,j,ip1,n) + vvel(k,j,ip2,n);
+ });
+ // } else {
+ // Implement support for 3D domain here
+ }
+ //------------------------------------------------------------------------------------------------
+ // Apply hyperdiffusion
+ parallel_for( SimpleBounds<4>(nz,ny,nx,nens) , YAKL_LAMBDA (int k, int j, int i, int n) {
+ temp(k,j,i,n) += -1 * dt * hd_temp(k,j,i,n) / ( 16 * hd_timescale );
+ rhod(k,j,i,n) += -1 * dt * hd_rhod(k,j,i,n) / ( 16 * hd_timescale );
+ rhov(k,j,i,n) += -1 * dt * hd_rhov(k,j,i,n) / ( 16 * hd_timescale );
+ rhol(k,j,i,n) += -1 * dt * hd_rhol(k,j,i,n) / ( 16 * hd_timescale );
+ rhoi(k,j,i,n) += -1 * dt * hd_rhoi(k,j,i,n) / ( 16 * hd_timescale );
+ nliq(k,j,i,n) += -1 * dt * hd_nliq(k,j,i,n) / ( 16 * hd_timescale );
+ nice(k,j,i,n) += -1 * dt * hd_nice(k,j,i,n) / ( 16 * hd_timescale );
+ uvel(k,j,i,n) += -1 * dt * hd_uvel(k,j,i,n) / ( 16 * hd_timescale );
+ vvel(k,j,i,n) += -1 * dt * hd_vvel(k,j,i,n) / ( 16 * hd_timescale );
+ });
+ //------------------------------------------------------------------------------------------------
+}
diff --git a/components/eam/src/physics/crm/pam/pam_output.h b/components/eam/src/physics/crm/pam/pam_output.h
new file mode 100644
index 000000000000..74b724b8ea70
--- /dev/null
+++ b/components/eam/src/physics/crm/pam/pam_output.h
@@ -0,0 +1,188 @@
+#pragma once
+
+#include "pam_coupler.h"
+
+// Compute horizontal means for feedback tendencies of variables that are not forced
+inline void pam_output_compute_means( pam::PamCoupler &coupler ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto &dm_host = coupler.get_data_manager_host_readwrite();
+ int crm_nz = dm_device.get_dimension_size("z" );
+ int crm_ny = dm_device.get_dimension_size("y" );
+ int crm_nx = dm_device.get_dimension_size("x" );
+ int nens = dm_device.get_dimension_size("nens");
+ int gcm_nlev = coupler.get_option("gcm_nlev");
+ //------------------------------------------------------------------------------------------------
+ // Get current CRM state
+ auto crm_rho_d = dm_device.get("density_dry");
+ auto crm_rho_v = dm_device.get("water_vapor");
+ auto crm_rho_c = dm_device.get("cloud_water");
+ auto crm_rho_r = dm_device.get("rain");
+ auto crm_rho_i = dm_device.get("ice");
+ auto crm_nc = dm_device.get("cloud_water_num");
+ auto crm_ni = dm_device.get("ice_num");
+ auto crm_nr = dm_device.get("rain_num");
+ auto crm_qm = dm_device.get("ice_rime");
+ auto crm_bm = dm_device.get("ice_rime_vol");
+ //------------------------------------------------------------------------------------------------
+ // Create arrays to hold the current column average of the CRM internal columns
+ dm_device.register_and_allocate("qv_mean", "domain mean qv", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("qc_mean", "domain mean qc", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("qi_mean", "domain mean qi", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("qr_mean", "domain mean qr", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("nc_mean", "domain mean nc", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("ni_mean", "domain mean ni", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("nr_mean", "domain mean nr", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("qm_mean", "domain mean qm", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("bm_mean", "domain mean bm", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("rho_d_mean", "domain mean rho_d", {gcm_nlev,nens},{"gcm_lev","nens"});
+ dm_device.register_and_allocate("rho_v_mean", "domain mean rho_v", {gcm_nlev,nens},{"gcm_lev","nens"});
+ auto qv_mean = dm_device.get("qv_mean");
+ auto qc_mean = dm_device.get("qc_mean");
+ auto qi_mean = dm_device.get("qi_mean");
+ auto qr_mean = dm_device.get("qr_mean");
+ auto nc_mean = dm_device.get("nc_mean");
+ auto ni_mean = dm_device.get("ni_mean");
+ auto nr_mean = dm_device.get("nr_mean");
+ auto qm_mean = dm_device.get("qm_mean");
+ auto bm_mean = dm_device.get("bm_mean");
+ auto rho_d_mean = dm_device.get("rho_d_mean");
+ auto rho_v_mean = dm_device.get("rho_v_mean");
+ //------------------------------------------------------------------------------------------------
+ // We will be essentially reducing a summation to these variables, so initialize them to zero
+ parallel_for("Initialize horzontal means", SimpleBounds<2>(gcm_nlev,nens), YAKL_LAMBDA (int k_gcm, int iens) {
+ qv_mean(k_gcm,iens) = 0;
+ qc_mean(k_gcm,iens) = 0;
+ qi_mean(k_gcm,iens) = 0;
+ qr_mean(k_gcm,iens) = 0;
+ nc_mean(k_gcm,iens) = 0;
+ ni_mean(k_gcm,iens) = 0;
+ nr_mean(k_gcm,iens) = 0;
+ qm_mean(k_gcm,iens) = 0;
+ bm_mean(k_gcm,iens) = 0;
+ rho_d_mean(k_gcm,iens) = 0;
+ rho_v_mean(k_gcm,iens) = 0;
+ });
+ //------------------------------------------------------------------------------------------------
+ // Compute horizontal means
+ real r_nx_ny = 1._fp/(crm_nx*crm_ny); // precompute reciprocal to avoid costly divisions
+ parallel_for("Horz mean of CRM state", SimpleBounds<4>(crm_nz,crm_ny,crm_nx,nens), YAKL_LAMBDA (int k_crm, int j, int i, int iens) {
+ int k_gcm = gcm_nlev-1-k_crm;
+ real rho_total = crm_rho_d(k_crm,j,i,iens) + crm_rho_v(k_crm,j,i,iens);
+ atomicAdd( qv_mean(k_gcm,iens), (crm_rho_v(k_crm,j,i,iens) / rho_total) * r_nx_ny );
+ atomicAdd( qc_mean(k_gcm,iens), (crm_rho_c(k_crm,j,i,iens) / rho_total) * r_nx_ny );
+ atomicAdd( qi_mean(k_gcm,iens), (crm_rho_i(k_crm,j,i,iens) / rho_total) * r_nx_ny );
+ atomicAdd( qr_mean(k_gcm,iens), (crm_rho_r(k_crm,j,i,iens) / rho_total) * r_nx_ny );
+ atomicAdd( nc_mean(k_gcm,iens), crm_nc (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( ni_mean(k_gcm,iens), crm_ni (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( nr_mean(k_gcm,iens), crm_nr (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( qm_mean(k_gcm,iens), crm_qm (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( bm_mean(k_gcm,iens), crm_bm (k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( rho_d_mean(k_gcm,iens), crm_rho_d(k_crm,j,i,iens) * r_nx_ny );
+ atomicAdd( rho_v_mean(k_gcm,iens), crm_rho_v(k_crm,j,i,iens) * r_nx_ny );
+ });
+ //------------------------------------------------------------------------------------------------
+}
+
+
+inline void pam_output_copy_to_host( pam::PamCoupler &coupler ) {
+ using yakl::c::parallel_for;
+ using yakl::c::SimpleBounds;
+ using yakl::atomicAdd;
+ auto &dm_device = coupler.get_data_manager_device_readwrite();
+ auto &dm_host = coupler.get_data_manager_host_readwrite();
+ auto nens = coupler.get_option("ncrms");
+ auto crm_nz = coupler.get_option("crm_nz");
+ auto crm_nx = coupler.get_option("crm_nx");
+ auto crm_ny = coupler.get_option("crm_ny");
+ auto gcm_nlev = coupler.get_option("gcm_nlev");
+ //------------------------------------------------------------------------------------------------
+ auto qv_mean = dm_device.get("qv_mean");
+ auto qc_mean = dm_device.get("qc_mean");
+ auto qi_mean = dm_device.get("qi_mean");
+ auto qr_mean = dm_device.get("qr_mean");
+ auto nc_mean = dm_device.get("nc_mean");
+ auto ni_mean = dm_device.get("ni_mean");
+ auto nr_mean = dm_device.get("nr_mean");
+ auto qm_mean = dm_device.get("qm_mean");
+ auto bm_mean = dm_device.get("bm_mean");
+ auto rho_d_mean = dm_device.get("rho_d_mean");
+ auto rho_v_mean = dm_device.get("rho_v_mean");
+ auto gcm_forcing_tend_temp = dm_device.get("gcm_forcing_tend_temp" );
+ auto gcm_forcing_tend_rho_d = dm_device.get