Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix barotropic thickness on edge in split explicit subcycling #5195

Conversation

mark-petersen
Copy link
Contributor

@mark-petersen mark-petersen commented Sep 21, 2022

Barotropic thickness at the edge must match the sum of baroclinic thicknesses at that edge. This PR corrects a long-standing problem of noise in the vertVelocityTop variable, which is the Eulerian (fixed-frame) velocity for momentum. The vertTransportVelocityTop variable for tracer and thickness transport was smooth, both before and after this change.

Currently in the code, the thickness at an edge used for volume (SSH) transport in the barotropic subcycling is

thicknessSum = sshEdge + min(bottomDepth(cell1), bottomDepth(cell2))

That is, the thickness at an edge is taken as the shallower of the two neighboring cells. This was inconsistent with the baroclinic mode, which takes the average thickness of the two neighboring cells.

This PR is climate changing.

[CC]

@mark-petersen
Copy link
Contributor Author

The change on each edge for the thickness used for barotropic transport is up to half a layer, so 250/2 = 125m in deep cells:
bottomDepthSketch
Note this problem occurs when using the combination of partial bottom cells and split time stepping, which is the case for all our standard runs.

@mark-petersen
Copy link
Contributor Author

before and after:
image
image

@mark-petersen
Copy link
Contributor Author

Thanks to @dengwirda for discussions on this topic, and for questioning that line that specifies the barotropic thickness.

@mark-petersen
Copy link
Contributor Author

For documentation, I'll add here that noise in the vertical velocity does still exist, but to a much lower extent. Here is a simulation with the most plain options I could imagine:

  • RK4
  • z-level
  • no partial bottom cells
  • GM off
  • Redi off
  • config_pressure_gradient_type = 'pressure_and_zmid'

This is 57 days into the simulation:
image
note colorbar is tighter than the images above. Images with split-explicit look similar.

@mark-petersen mark-petersen force-pushed the mark-petersen/ocean/barotropic-edge-thickness branch from 3202079 to ef076d4 Compare September 24, 2022 03:10
@dengwirda
Copy link
Contributor

@mark-petersen and all, a QU30km MPAS-O standalone spin-up (4 years) confirms the improvement here, with vertVelocityTop now matching vertAleTransportTop (for a z-star coord.). Prior to this PR, vertVelocityTop contained significant grid-scale noise, due to the 2d vs 3d thickness incompatibility discussed above.

This change does not appear to degrade the energetics of the flow, with both the max and rms kinetic energy in the same range before and after this change. (Spurious damping of K had been a problem with other attempts to fix this issue).

Note that while we're focusing on improvements to vertVelocityTop here, this PR actually changes the way the (split) horizontal dynamics evolve, so changes are not isolated to vertical diagnostics. vertVelocityTop is simply a place where the change is most visible.

New |u|, vertVelocityTop and vertAleTransportTop:

qu30_speed_new
qu30_vert_vel_top_new
qu30_vert_ale_top_new

Old vertVelocityTop and vertAleTransportTop:

qu30_vert_vel_top_old
qu30_vert_ale_top_old

@dengwirda
Copy link
Contributor

Sticking with the question of barotropic vs baroclinic compatibility in the split integration schemes, we should also be careful with the thickness forcing terms S_h:

dh/dt + div(u*h) + dw/dz = S_h.

The rhs S_h term contains any thickness fluxes in the coupled system (runoff, evap., precip., etc).

When we integrate vertically to get the barotropic subsystem (for the split subcycles), I don't believe the forcing drops out:

dH/dt + div(U*H) = S_H, with S_H = SUM_layers(h*S_h) / SUM_layers(h).

Currently, we do take similar care for the forcing term for momentum (barotropicForcing), but not for thickness. Can we do:

barotropicForcingU = SUM_layers(h*S_u) / SUM_layers(h),
barotropicForcingH = SUM_layers(h*S_h) / SUM_layers(h),

to be consistent. The 2d thickness forcing is probably fairly small in most cases, but this seems like another opportunity for a 2d vs 3d thickness mismatch as per the PBC issue above. Would require coupled E3SM runs to test, not MPAS-O standalone where S_h = 0 usually.

@jonbob
Copy link
Contributor

jonbob commented Oct 3, 2022

I ran a 30-year case with this PR test-merged into recent next (and including PR #5172). MPAS-Analysis output comparing this run with a previous one which also included PR #5172 is available at:

20220929.PR5195.chrysalis

If anyone would prefer a direct comparison with obs as well, please let me know. @xylar and @dengwirda, including you in this review as well - thanks

@maltrud
Copy link
Contributor

maltrud commented Oct 3, 2022

@jonbob it looks to me like the AMOC might be a bit stronger now. I didn't see anything that raised a red flag for me.

@vanroekel
Copy link
Contributor

The changes here look mostly minor or to within the noise of different ensemble members. As @maltrud notes AMOC may be a bit stronger. This is definitely a bug and the fix is obvious. I'll approve this now.

Copy link
Contributor

@vanroekel vanroekel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Approving based on visual inspection and testing from @jonbob and @mark-petersen

@jonbob
Copy link
Contributor

jonbob commented Oct 3, 2022

Thanks @maltrud

@jonbob
Copy link
Contributor

jonbob commented Oct 3, 2022

also analysis vs obs:

20220929.PR5195.chrysalis_vs_obs1-30

@mark-petersen
Copy link
Contributor Author

Thanks @jonbob for the testing. With this PR, the AMOC is a bit stronger, but longer runs may turn out similar.
image
Other statistics in the 30-year run also look similar between master and this PR, for example ocean heat content and transports. At a glance:
image

jonbob added a commit that referenced this pull request Oct 4, 2022
…t (PR #5195)

Fix barotropic thickness on edge in split explicit subcycling

Barotropic thickness at the edge must match the sum of baroclinic
thicknesses at that edge. This PR corrects a long-standing problem of
noise in the vertVelocityTop variable, which is the Eulerian
(fixed-frame) velocity for momentum. The vertTransportVelocityTop
variable for tracer and thickness transport was smooth, both before and
after this change.

Currently in the code, the thickness at an edge used for volume (SSH)
transport in the barotropic subcycling is
   thicknessSum = sshEdge + min(bottomDepth(cell1), bottomDepth(cell2))
That is, the thickness at an edge is taken as the shallower of the two
neighboring cells. This was inconsistent with the baroclinic mode, which
takes the average thickness of the two neighboring cells.

[CC]
@jonbob
Copy link
Contributor

jonbob commented Oct 4, 2022

passes sanity testing, merged to next

@jonbob jonbob merged commit 4dcc8db into E3SM-Project:master Oct 5, 2022
@jonbob
Copy link
Contributor

jonbob commented Oct 5, 2022

merged to master and expected DIFFs blessed

@xylar xylar deleted the mark-petersen/ocean/barotropic-edge-thickness branch October 5, 2022 16:47
jonbob added a commit that referenced this pull request Oct 7, 2022
jonbob added a commit that referenced this pull request Oct 10, 2022
Add missing variables to omp pragmas from PR #5195

PR #5195 added new code, but had some variables missing from omp pragmas
that cause SMS_D_Ld3.T62_oQU120.CMPASO-IAF.chrysalis_intel failures.
This fixes that issue.

[non-BFB] for threaded tests
jonbob added a commit that referenced this pull request Oct 11, 2022
Add missing variables to omp pragmas from PR #5195

PR #5195 added new code, but had some variables missing from omp pragmas
that cause SMS_D_Ld3.T62_oQU120.CMPASO-IAF.chrysalis_intel failures.
This fixes that issue.

[non-BFB] for threaded tests
mark-petersen added a commit to mark-petersen/MPAS-Model that referenced this pull request Oct 26, 2022
mark-petersen added a commit to mark-petersen/MPAS-Model that referenced this pull request Oct 26, 2022
jonbob added a commit that referenced this pull request Dec 13, 2022
…5356)

Move bottomDepthEdge calculation to single loop over all edges

After #5195 was merged, the MPAS-Ocean standalone test
ocean/baroclinic_channel/10km/decomp_test failed to match between 4 and
8 partitions, but only for intel optimized. All compass nightly suite
tests passed for gnu debug, gnu optimized, intel debug.

This PR solves the problem by merging the computation of bottomDepthEdge
into a single edge loop. Previously it was split into two loops,
1:nEdgesOwned (with many other calculations) and another from
nEdgesOwned+1:nEdgesArray(4). The intel optimized compiler must have
changed order-of-operations in these two loops for different partitions.

Fixes #5219
[BFB]
jonbob added a commit that referenced this pull request Dec 14, 2022
Move bottomDepthEdge calculation to single loop over all edges

After #5195 was merged, the MPAS-Ocean standalone test
ocean/baroclinic_channel/10km/decomp_test failed to match between 4 and
8 partitions, but only for intel optimized. All compass nightly suite
tests passed for gnu debug, gnu optimized, intel debug.

This PR solves the problem by merging the computation of bottomDepthEdge
into a single edge loop. Previously it was split into two loops,
1:nEdgesOwned (with many other calculations) and another from
nEdgesOwned+1:nEdgesArray(4). The intel optimized compiler must have
changed order-of-operations in these two loops for different partitions.

Fixes #5219
[BFB]
singhbalwinder pushed a commit that referenced this pull request Jan 15, 2025
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736920654 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736920445 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736908121 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736906357 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736899098 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736898892 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736898759 -0800

parent 0075116ccf04bbcafb5129a0a9cdf5e733d16802
author Mark Petersen <mpetersen@lanl.gov> 1663706865 -0600
committer Balwinder Singh <balwindersingh@gmail.com> 1736893569 -0800

Fix barotropic edge thickness, assumes z-type coordinate

add bottomDepthEdge calculation

Change from z-level to general vertical coordinate

Add missing variables to omp pragmas

Corrects issue #5206, sea ice snicar_ad error

Fixes interpolation of shortwave parameters for the 5-band scheme.
This error was introduced in the Snicar_ad feature implementation.

nonBFB

Add Time coordinate

This variable has CF-compliant units and calendar.  To accommodate
the CF-compliant units, a new config option has been added to
supply the reference time, 0001-01-01 by default.

Add Time_bnds variable

Fix typos and white-space in time series stats

Clean up input and output field names

Always add Time and Time_bnds to time series stats

Compute Time and Time_bnds in time series stats

Add standard name

Change defunct Fortran labels to end do

change cycle to exit

Change goto to backwards do loop

Hommexx/SL: Implement HOMMEXX_MPI_ON_DEVICE=Off in SL transport.

Previously, the unstructured SL-internal communication rounds were done using
device pointers regardless of the value of HOMMEXX_MPI_ON_DEVICE. Now they are
done on host if HOMMEXX_MPI_ON_DEVICE is true.

tweaks for spock

cache file

fix build of unit tests without pio

added define for hip build

more fixes, intermediate

intermediate

revert nonsgnificant changes

reverting

clean up

clean up

reverting

fix ifdef, add prnts for the last pair

mods for cmake file, comment out printing pairs

set the last pair for hip manuall

add flags

newer environment

update kokkos build made with a new env

minor, 1 comment on LB, commenting out pow call that wasnt implemented

add crusher cmake

fix bug from merging

clean comments, remove if-statement in preqx push logic

switch to c++17

remove unused cmake var enable_cuda

change if HIP to ifdef HIP

fix cprnc path when not using cprnc_dir

summit bfb file for tests

minor fix in summit cmake

theta no-compose build

remove opt flag for 1 remap file

fix to the team size in diagnostics, not verified on crusher yet

diagn for repro_sum, not verified on crusher

fixing in vs inout issue for nstep_c

remove noopt flag from 1 file

use inout for some vars

start on organizing gpu build vars

remove duplicated code and fix ifdef

replace Hommexx_Cuda with HommexxGPU

clean up and fixes, to be tested

fixing bfb; confirmed for fhs1 test

Hommexx: Fix some CMake issues.

Using link_to_kokkos from Kokkos.cmake provides the correct include and lib(64)
directories. Use this macro for Hommexx targets.

As a result, crusher-gpumpi.cmake can be trimmed down substantially. Remove
lines in that file depending on E3SM_KOKKOS_PATH. Instead, set E3SM_KOKKOS_PATH
in your configuration script.

Hommexx/SL: Clean up use of HOMME_ENABLE_COMPOSE.

Include Config.hpp to see the symbol. Move some uses around.

Hommexx/SL: Generalize Cuda to GPU spaces.

Handle Cuda, HIP, and perhaps SYCL.

Hommexx: Update DIRK solver's tridiag solver to handle HIP.

Hommexx: Fix a warning in VR.

Hommexx: Make more link lines use link_to_kokkos.

link_to_kokkos can handle the lib layouts of various versions of
Kokkos, wherease linking directly to libkokkos is an error with
some Kokkos versions.

Homme: Slightly relax a tolerance in a unit test.

Homme/SL: Slightly relax a tolerance in a unit test.

Hommexx/SL: Fix a variable for lambda capture in two spots.

Worked on various Cuda devices, but Crusher needs this to be explicitly not a
reference to capture it correctly.

Hommexx: Get standalone Homme to build its own Kokkos.

This capability seemed broken because algorithms/src wasn't being included. Fix
this.

Also update summit-gpumpi.cmake to always build COMPOSE.

Hommexx: CUDA_ARCH equivalent in GllFvRemap.

fix preqx build

remove string args from prim_state calls

disable code around energy vars that are not used in cxx version

remove string args from glob_norms calls

a comment

cleaning team size logic

clean 1 line for HOMMEXX_ENABLE_GPU

address comments

fix homme standalone build not via cime

Hommexx/SL: (Re)fix a variable for lambda capture in two spots.

HIP can't handle the reference, but we still need it for the F90 dycore with
HORIZ_OPENMP enabled, since in that code path we enter all of our routines
within a threaded region and so cannot create new Kokkos::Serial views during
time stepping. (Recall Kokkos::Serial is necessary b/c Homme's threading is not
compatible with Kokkos::OpenMP.) Thus, if COMPOSE_PORT is not defined, take a
reference.

Updates the settings for JRA data

Removes the hard-coded start and end year for JRA from 1958/2016 to be
the start/end values specified by DATM_CLMNCEP_YR_START/DATM_CLMNCEP_YR_END.

[BFB]

Adds start/end years

Updates the two-way river-ocean test

Shorten test-job-id with --jenkins-id Jenkins script option

removed device_resident and fixed OpenACC in ocean eff density

   replaced OpenACC device_resident with usual enter/exit data
   added missing OpenACC data directives for pool data
   removed mesh pool and replaced with ocn mesh module
   other general cleanup

replaced OpenACC device_resident statements in ocn eq state routines

  replaced device_resident with the more generic create/delete data
   statements for portability
  some minor cleanup

replace device_resident and other OpenACC fixes for ocn vel diagnostics

   replaced OpenACC device_resident with more common create/delete
   fixed missing private directives in some cases
   other cleanup for long lines and formatting

replaced device_resident in ocn diag compute KPP routine

   also significant cleanup/formatting

replace device_resident and other fixes for ocn vert transport

   replaced OpenACC device_resident with more common create/delete
   removed meshPool and use mesh module variables
   initial (but disabled) OpenACC conversion of ALE routine
   a lot of cleanup of long lines, etc

back out some cleanup changes to preserve bit-for-bit

  some previous changes were not bit-for-bit due to aggressive
  optimization in Intel compiler, so backing out two minor changes

re-inserted some omp parallel directives

Add mpas_ocn_thick_ale.f90 for pgigpu runs

Add missing variables to omp pragmas, from PR #5195

Double PEs to prevent walltime timeouts on Ascent

Increase number of nodes for GPMPAS-JRA compset from 2 to 4

Increase nnodes for GPMPAS-JRA compset from 2 to 4 on Crusher

update to github.io in ocean and sea ice

update to github.io in mpas-albany-landice

update to github.io in mpas-framework

update to github.io in c and new additions

Remove unneeded mesh declaration

add/change input files for CBGCv2 land-atm 1850 compset for v2 CBGC spin-up runs

add/change input files for CBGCv2 land-atm 20tr compsets for CBGCv2 simulations

update paramfile file

fix errors in PR #5180

fix errors in PR #5180

fix errors in PR #5180

Changes to for ELM's SP mode

Old surface dataset file, which does not include the cft dimension, is used for SP mode.
The new surface dataset, which includes the cdf dimension, is used for BGC turned on.

Updates the test

Modify ELMBuildNamelist.pm to check 'use_cn' when setting 'fsurdat'

Modified the setting of create_crop_landunit to be true for any run with hgrid="r05" use_cn = ".true."

New flanduse_timeseries file on set if use_cn = .true.

Updated namelist build script to check value of "use_cn" when setting elm paramfile.
New clm_params_c211124.nc is only set if "use_cn = .true."

Update config_compsets.xml

Fix runoff configuration in the BGC land-atm compsets

Fixed using wrong file for ECA, crop, and FATES runs.
Removed "nu_com = 'RD'" check since it isn't use in lnd_in but is the default setting
in ELM source code.

move define to hommexx_config

revert changes in config files

Change name of cime directory

Add ww3 component directory from cesm

Consolidate ww3 cpl directory and delete .svn directories

 -

Update ww3 buildlib to agree with consolidated cpl directory structure

Change wave model component name from ww to ww3

Change order of servers in config_inputdata.xml

Setup test wave compsets

Bug fixes to get case.build to work with ww3

 - Fix syntax mistake in ww3 buildlib

 - Change name of user_nl_ww to user_nl_ww3

 - Update ww3 buildnml with correct component name and user
   namelist file name

Add test model grid to config_grids.xml for wave+datm case

Change parameters in namelist_definition_ww3.xml

 - Switch restart file so model starts from calm conditions

 - Set output frequency to 3 hours

Add "only" option to ww3 module use statements in wav_comp_mct.F90

 - This makes it more clear which variables are needed in the
   coupling interface and helps inform what changes are necessary
   for incorporating the latest version of ww3

Add cesm_splice_maps.py for creating ocean/wave remapping files

 - Takes in both a  bilinear and nearest neighbor remapping file
   and uses bilinear for the interior ocean and nearest neighbor
   in coastal regions.

   This avoids coupling issues between the waves and ocean caused by
   using bilinear interpolation at the coast.

 - Provided by Qing Li

Modify ww3 NetCDF output in w3iogomd.f90

 - Make output consistent with newer version

Remove unused subroutine from use statement in wav_comp_mct.F90

Reduce output interval to 1 hour in ww3 namelist_definition_ww3.xml

Switch output variables in wav_comp_mct.F90

 - Turn on peak frequency and direction

 - Turn off mean period and direction

Add WW3 v5 source code

 - Only add files in the bin, ftn and inp directories for now

Add WW3 v5 compliation to buildlib

 - Pre-processes wavewatch *.ftn files

 - Also, complies wavewarch executables (not necessary, this step
   should be taken out later)

 - Calls E3SM Makefile to create libwav.a

 - Libaray creation fails right now because of version incompatibilites
   with exsiting mct interface (wav_comp_mct.F90)

Comment out coupling variables from CESM

 - Will add these back in later

Update wav_comp_mct.F90 for WW3 v5

 - Comment out Langmuir enhancement (lamult) and mixed layer depth
   (HML) variables, since they are not in the new wavewatch source
   code.

 - Modify the field output specification. The outputs are now grouped
   into 9 different groups. The flag array to turn them on/off is
   now 2D

 - Add field output flag variables for coupling variables (not used)

 - Update call to w3init

 - Increase size of odat variable to account for coupling variables

Use the correct mod_def.ww3 file in buildnml

Add WW3 F90 source code files

1

Get WW3 gridded output working

 - The NOTYPE variable wasn't being assigned and defaulted to 0.
   The loop that goes over output types to call the write routines
   was being skipped over.

 - Setting the NOTYPE variable to 7 in wav_comp_mct.F90 fixes the
   issue.

Update WW3 bin directory after reinstall

Update MPI complier name in WW3 comp and link files

Don't remove *.f90 (lowercase f) files when running w3_clean

Minor updates to WW3  buildlib

 - Add w3_clean command

 - Get rid of redundant environment variable specification

 - Switch NETCDF_CONFIG variable to intel

Update pre-processed WW3 .F90 files for new switch configuration

Move w3cesmmd.f90 to WW3 cpl directory

 - This is the cleanest way to include this file in the build
   process

Revert "Don't remove *.f90 (lowercase f) files when running w3_clean"

This reverts commit 3a4455f081848e31f193b8adfd8e5c19396a34f1.

Change WW3 link and comp files to gnu

Update WW3 buildlib

 - Use WW3's w3_source script to pre-process the source code
   for the ww3_shel program only, instead of using w3_make to
   pre-process and compile all programs

 - w3_source creates a .tar.gz file that must be untared before
   calling the e3sm makefile

Change WW3 restart file output interval to 5 days

Update WW3 *.F90 files

 - Only include ww3_shel files since those are now the only ones
   generated by buildlib

 - ww3_grid files are also included, since that program was built
   using w3_make to create the mod_def.ww3 file

Add the wwatch3.env file

Update build scripts and coupler interface for WW3 v6.07

 - Add argument to call to w3init in wav_comp_mct.F90

 - Update paths to ww3 repo in buildlib. Also handle setting up the
   wwatch3.env file (if it doesn't already exist) by calling w3_setup

 - Point to correct mod_def.ww3 and restart.ww3 files in buildnml and
   namelist_definition_ww3.xml

Comment WW3 buildlib and remove debug output

Add WW3 github repo as submodule

Add switch_E3SM to control WW3 model options

 - This file will get copied to WW3 repo and when  w3_setup is run,
   it will cause this to be used as the switch file for pre-processing
   the .ftn source code

Add some file removal cleanup to WW3 buildlib

 - This makes sure that the WW3 git working directory is clean (no
   untracked files) after building

Remove ww3 v5 directory

Move old CESM WW3 .f90 files to ww3_cesm directory

Partially add wave radiation stress terms (Sxx, Sxy, Syy)

 - Add fields in seq_flds_mod.F90

 - Set w2x indicies in ww3_cpl_indicies.F90

 - Set x2o indicies in mpaso)cpl_indices.F90

 - Pack variables in wav_comp_mod.F90

 - No modifications to MPAS-O source code yet

Add wave radiation stress terms to ocn_comp_mct.F

 - Unpack variables from attribute vectors

Include ssh as a coupled variable between wav and ocn

Remove all .svn directories

Updates to WW3 buildlib

 - Pass in extra variables when invoking make

Fix compliation bug in wave_comp_mct.F90

 - Add additional variables in use statements

Changes to config_compilers.xml and config_machines.xml

Fix typo in ocn_comp_mct.F

Add compsets and grids for datm/wav testing

 - Use unstructured WaveWatchIII meshes

Add CFSv2 wind config

Update WW3 source to global unstrctured branch

Update ww3 buildlib to include compling ww3_grid

Update ww3 switch_E3SM

Update buildnml for ww3

Configure buildnml for ww3_grid

Build additional namelist for ww3_grid

 - Also run ww3_grid in buildlib

Begin to add new WW3 compset

Updates to MPAS-O coupling fields

 - Remove radiaiton stress
 - Add stokes drift

Update WW3 coupling fields

 - remove radiation stress

Finish adding wQU100 and bug fixes

 - read from ww3_in instead of wav_in in wav_comp_mct.F90

 - remove UOST from switch file for now

Update WW3 submodule

Fix w3_setup input file in ww3 buildlib

Move ww3_grid build after main make call in ww3 buildlib

Add waves to cime_config files

Remove cime_add-ww3

Update cime submodule hash and url/branch in .gitmodules

Update description in ww3 config_component.xml

Bug fix for UWAV-SANDY compset long name

Add cmake buildlib file, remove legacy buildlib

Progress getting ww3 to build

Comment to indicate where mesh.msh file is needed

Add more comments

Comment out ww3_grid build

Added compset F20TRC5-CMIP6-WW3 for Wavewatch with F case

Update WW3 submodule and related changes

 - E3SM WW3 branch has been rebased on most recent NOAA-EMC develop branch
   (fc3576f413acfe44f4c81e95eefd8ff7b5b21cdd)

   On top of this there are two important modifications:

     1) Fixes for rotated global unstructured meshes
     2) Re-organization of ww3_grid code into new module and subroutine

 - wav_comp_mct.F90 has been modified to update the interfaces to w3init
   and w3wave from the develop branch.

 - buildlib_cmake has be updated for the develop branch, w3_setup no
   longer uses the -c argument

 - buildlib_cmake now runs w3_source on the ww3_grid target to generate
   source code for the  w3gridmd module used to call w3grid from within
   wav_comp_init

 -  The wave_comp_init subroutine in wave_comp_mct.F was modified to
    call w3grid to generate the mod_def.ww3 and other related files
    prior to running the wave model. This is done serially.

Update compset for wQU225EC60to30 grid

Remove ww3_shel and ww3_grid .F90 files before building

Update build-namelist

 - Add path to .msh file for specific grid in ww3_grid.nml

 - Add path to station file

 - Turn on k propagation (for current coupling)

 - Reduce minimum depth to .25

 - Add ww3.input_data_list output

Add station output capability in wav_comp_mct.F90

Rotate winds into rotated wave coordinate system

Add namelist options for RTD grid

Update E3SM WW3 switches

Add CFSR_wQU225EC60to30 resolution

Add CFSR datm compset

Add wavewatch output and coupling flags to namelist

Rotate stokes drift vectors to standard pole

Clean up wav_comp_mct.F90 screen output (restrict to iaproc)

Adjust w2x variables for rotated grid

 - localize AnglD variable

 - rotate Stokes drift vector back to standard pole

 - switch to dedicated output processor to prevent overwriting
   local mean wave properties on ouptur processor

Add Stokes drift output to ww3_in namelist

Fix typo in stokes drift variable name in ocn_comp_mct.F

Add support for T62_oEC60to30v3_wQU225EC60to30 grid

 - Remove wDE4

Add hooks to enable active wave in mpas-o

Add pe layout for active wave component in current low-res G case

Clean up wav_comp_mct.F90: remove outfreq and histwr

Improve WW3 I/O

 - Make point and gridded output intervals namelist parameters

 - Write out restart files on driver interval

 - Read in restart files for continue and branch runs

Make ww3 restart equilvant to manual linking restart.ww3

Implement BFB restarts with ocean couplings

Add ocean to wave ssh field in E3SM's seq_flds_mod.F90

Add CFSv2 and CFSR datm configurations to E3SM datm

Subtree merge from MPAS-Model hash in OceanNGD/add-ww3

Update bld script to match Registry changes

Clean-up from subtree merge

Update cime branch to jonbob/e3sm/support-add-ww3 in .gitmodules

Subtree merge of Stokes drift MPAS-Ocean changes

Add partitioned Stokes drift information to ww3 namelists

Add partitioned stokes drift variable coupling

Fix compile bugs in wav_comp_mct.F90 and ww3_cpl_indices.f90

Add significant wave height as coupled field

Change wave coupling variable names

 - get rid of wave prefix

Remove original surface stokes drift coupling fields

Remove outdated comments

Add mpaso namelist and streams options for active wave langmuir

Add significant wave height coupling variable to ww3 namelist

Broadcast partitioned Stokes drift information

 - Prevents error with setting gridded output and coupling variable
   flags

Zero out significant wave height at non-sea points

Update WW3 submodule

 - Rebase RTD fixes

Remove USP as an output variable for ww3

Update cime submodule

 - Fix MCT build error that occurs with newer version of cime

Updates for G case + waves with Stokes drift coupling

 - Change PE layout

 - Add Langmuir number to ocean streams

 - Turn off wave model output to avoid MPI errors

 - Uncomment KPP enhancement factor calculations

 - Modify compset to include active sea ice

Update WW3 submodule

 - Extend the length of mesh filename variable

Fixes log file redirection error

 - Caused by the hardcoded unit number in w3grid

 - Also updates WW3 submodule to include w3grid calling argument
   change (unit numbers added)

Add ice damping limits and relative wind namelist parameters

Add UOST inputs to namelist and switch files

Add e3sm_wav_developer test group

Add unit number to read_stations_file subroutine

Use IC2/IS2 switches in WW3

Reorganize Langmuir number calculations

Add boundaryLayerDepth to restart stream

Add peak direction and frequency to coupling variables

Add B-case with waves compset

Add stokes drift omp private variables in mpas_ocn_vmix_cvmix.F

Fix typo for peak frequency/direction import

Update atm -> wav mapping file between n30np4 and wQU225EC60to30

Add WAV_WCYCL1850 compset

Set inflags2 = inflags1 in wav_comp_mct.F90

 - Fixes bug that prevented reduction of input source terms based on
   ice concentration

Add peak freqency and direction to ww3 coupling variable output

Update active wave coupling flags

 - Add config_use_active_wave which is separate from config_cvmix_kpp_use_active_wave

 - config_cvmix_kpp_n_wavenumber_partitions -> config_n_stokes_drift_wavenumber_partitions

 - Add surface stokes drift routine

Add waveOutput stream

Set icedisp default to .false.

Add resolution for v2 LR + waves

Update v2 LR + waves configuration

Update cime submodule in .gitmodules

Remove legacy WW3 code from CESM

Remove legacy ChangeLog from CESM

Fix whitespace issue

Add TL319_EC30to60E2r2_wQU225EC30to60E2r2 grid

Remove old ww3 build code

Move WW3 pe layouts

Add JRA-forced G-case + waves

Rename wave only compsets

Change location of ww3 submodule to E3SM-Project

Add array for Stokes drift wave numbers

Move E3SM switch file to WW3 repo

Move WW3 submodule

Fix .gitmodules whitespace

Update comments in stokes drift subroutines

Remove unsed file: namelist_definition_ww3.xml

Remove old ww3.buildlib and ww3.buildnml files

Remove unsed compset in eam and mpas-ocean config_compsets.xml

Remove ww3a grids

Add wave namelist options to build-namelist-*

Remove QU100 wave grid

Rename wave compsets

Add SWAV to pe layouts for non wave fully coupled cases

Update tests

Fix duplicate declaration for ssh after rebase

Fix merge unresolved merge conflict in mpas-o makefile

Improve description and error checking for config_use_active_wave

Detect wave-only cases in builiding namelist

Add tbot field to datm for wave only cases

 - Avoid errors associated with not specifying a field from the
   stifld list in datm_com_mod.F90

Remove test file names for ww3 in config_archive.xml

Update WW3 submodule to incorporate seg fault fix

Clean up logging and related issues

 - Update WW3 submodule

Update mct ggrid

Clean up unit numbers and w3init call

Update WW3 submodule

Update extended coverage test for waves

Add langmuirNumber output to MPAS-O stream for waves runs

Fix minor ww3 buildnml issues

 - Add uost files to input file list

 - Add cesm_namelist as -infile to build-namelist call so user_nl_ww3
   modifications are considered

Fix CF compliance issues in Registry

Change JRA compset name

Fix ssh coupling fields

Check for single layer in Mannings implicit vertical mix

Add to single layer to vmix_tend_implicit_variable

add one layer option to Rayleigh drag vmix

add one layer option to standard vmix

Build our own versions of mo_load_coefficients

Build our own versions of mo_load_coefficients so that we can use
updated RRTMGP absorption coefficient input files with new solar source
variables (maintains backward compatibility with old files).

Use reduced g-points for radiation

add zeroing of velocity to init_mode

Fixes error that occurs when ocn_rof_two_way = .false. and an element of an un-allocated o2r_rx is passed as an argument.

Removing whitespace changes

Hommexx: Print element state data in various model failure situations.

In the log file, one will see, e.g.,

 46: Bad dphi, dp3d, or vtheta_dp; label: 'Vertical remap: Negative (or nan) layer thickness detected, aborting!'; see hommexx.errlog.320.46

Then, for example,

$ head hommexx.errlog.320.46
label: Vertical remap: Negative (or nan) layer thickness detected, aborting!
time-level 1
lat  5.182512623200946e-01 lon  1.413716694115407e+00
ie 3 igll 3 jgll 3 lev 122: bad vtheta_dp bad dp3d bad dphi
level                   dphi                   dp3d              vtheta_dp
    0 -1.894494653842581e+04  6.511250007060127e+01  9.449394802849884e+04
...
  122  1.061523656197460e+02 -3.435052907124674e+01 -8.703189416949290e+03
...

This feature is for the theta-l_kokkos dycore. preqx_kokkos has an empty
implementation of the check routine, so its behavior on failure is unchanged.

Hommexx/SL: Fix hipcc build errors in new MPI-host-ptr code.

This fixes some errors hipcc detects in the work done in PR #5223.

Update google cloud machine to gnu gcc-12.2

Allows config_use_snow_liquid_ponds = false with snow tracers

With the flag false, the routine to drain snow liquid was omitted.
However, this routine is needed to properly keep track of snow liquid.
he flag now only controls with the drained snow-melt contributes to
pond formation.

BFB

e3sm_compile_wrap.py: Move to E3SM repo

[BFB]

fixes for kokkos 3.7

fixes for kokkos 3.7

Fix precision for writing field attributes

It is necessary that the precision of floating-point attributes
match the precision of the field they are associated with.  In
particular, when the `_FillValue` has a different precision from
the field it is associated with, PIO complains and will not add
the attribute.

public get_pottemp needed for scream

Refactor E3SM build system handling of KOKKOS_OPTIONS

KOKKOS_OPTIONS is now set in universal.cmake and the intent is for
other macros to string(append KOKKOS_OPTIONS ...) when they want
to customize.

Also KOKKOS_OPTION is not set up to use Kokkos' CMake build system
directly instead of going through generate_makefile.

Fix EAM config options for CMIP6 family

The current repo on setting eam default physics and chemistry
fails to recognize compset with physics like EAM%CMIP6-1pctCO2
and EAM%CMIP6-4xCO2. It is now change to match all EAM%CMIP6.

[BFB] except for compsets with additional atm physics aside from CMIP6

Added WCYCL1850-1pctCO2 and WCYCL1850-4xCO2 tests to wcprod test suite

Change to use separate entries for eam config opts with co2 modifiers

Adding fsurdat and finidat for 1950 conditions for arcticx4v1 mesh

Adding mask specification for previous arcticx4v1 finidat

Updates mpas-seaice graph partition files used with oARRM60to10 mesh

Fixup minLevel* redi with conditional

Fixup minLevel* ocn_diagnostic_solve_ssh

Fixup minLevel* ocn_surface_bulk_forcing salt terms

Do not compute Haney number when layerThicknessEdge = 0

Bugfix inactive top cell surface restoring

Add PE layouts for cryo configurations on chrysalis

Changing G-cases to stacked layouts

update livermore computing machine entries

Fixes for quartz

whitespace cleanup

For pm-cpu/pm-gpu and alvarez at NERSC, update a few versions of modules
to be the default.  mpich, netcdf-related modules, and nvidia compiler.

For pm-cpu and pm-gpu, add env variable MPICH_COLL_SYNC=MPI_Bcast
to avoid hanging during writing of certain files at various node counts
and resolutions
[bfb]

Add narwhal, warhawk, onyx machines

Fixed column abort bug

Fixed bug where an abort in step_therm1 would not abort the code
Cleaned up abort code in column
Added abort code for step_snow

Changes to column abort for threading

Removed do loop exit for column for do loops that have omp parallel pragmas
Added anyAbort flag to determine if any cell in column loops aborts
Only write error messages for first abort in loop

Added explicit bounds to array assginments for bgc interface arrays.

add explicit bounds to fix FP invalid

Removing Switch from mpas component build-namelist

Changing cryo ocn/ice grids to current grids

Add missing time level index to ssh

Add tracers to dam break initialization

Add pm-cpu/pm-gpu to list of machine that save additional system logs to provenance
[bfb]

Adding support for JRA55 v1.5 datm/drof forcing through 2020

Adding IAF_JRA_1p5 as case for rhoa calculation

change const to double prec

case.run.sh: Some fixes and upgrades

Needs to use python3.
Now supports get_timing.

[BFB]

Use shell variables instead of passing things to python through env

Support for single process single cell standalone test cases

Added several options for controlling which thermodynamic processes are active
Added option for conservation AM to include ocean heat, mass, energy content in analysis with mixed layer model
Added another mixed layer ocean model that mimics ocean model in E3SM for standalone testing
Added another forcing option to support idealized single cell test case

Changes to how active process flags implemented

Active column processes now toggled with multiplication by real rather than conditional for performance
Renamed ncar_free_surface forcing type to som_varying_thickness

Added new namelist options to build-namelist

Removed trailing whitespace from Registry namelist options

Changing default GM for cryo v2 ocn meshes

Hommexx: Add missing team_barriers to omega computation.

share: Protect bool typedef from stdbool.h inclusion.

Some compiler wrappers add headers that include stdbool.h. In a C build,
including stdbool.h conflicts with the typedef. Guard the typedef with a check
for whether 'true' has already been defined.

On crusher, disable archiving of old test data

Apparently our shared disk quota for /gpfs/alpine/cli133/proj-shared
is pretty limited and some of our nighltly runs are started to fail
with "Disk quota exceeded". With this limitation, we cannot afford
to archive old test data.

[BFB]

Fix bfb fail in elem_ops_ut_test

- Add PhysicalConstant for tref_lapse_rate (Real type)
- Use (along with tref) in ElementOps.hpp instead of hardcoding,
changing from float type

fix uninitialized values in MMF cosp interface

change arg dims from pcols to npoints populate_cosp_subcol()

Pass cloud optics by argument instead of pbuf

add pgrad_correction to theta-l_kokkos

Hommexx: Limit HIP team size in debug builds.

Also clean up the team_num_threads_vectors routine to consolidate
the rules.

Fixes SCREAM issue 2058.

[BFB]

Fix logic for wave output in mpaso streams

CMake: Set E3SM CXX standard globally

Not in macros.

[BFB]

For cori-knl, add PE layout for a new test in e3sm_hi_res suite.
SMS.T62_SOwISC12to60E2r4.GMPAS-IAF.cori-knl_intel

Add OneAPI compiler to Chrysalis

Update to ucp-1.12.0 modules

Update cmake to id icpx as IntelLLVM 2022.1.0 for kokkos

Update to OneAPI v.2022.2.0.20220730

Limit omp functions to threaded builds for OneAPI

Update OneAPI modules

Add an update to e3sm-master's scorpio for implicit function decl

Update JLSE to currently default oneapi module

Prior oneapi/eng-compiler/2022.01.30.008 is now broken/deprecated

Add comments to mpas shallow water code

Added endwb to restart file if do_budgets = .true.
Required for runs that restart at the beginning of a month.

Don't think begwb is necessary as it is computed at
beginning of timestep and only used if nstep <= 1

Update YAKL and RRTMGP

update FFT for vriance transport

update YAKL submodule

misc YAKL FFT updates

bug fix for ESMT fft

Update YAKL to fix ArrayIR bug

fix FFT initialization

new fix for FFT init problem

move fft objects to vars.h

another FFT fix

add YAKL fence before FFT cleanup

update YAKL submodule to fix FFT cleanup issue

Add missing_value_mask to state vars

Add missing_value_mask to vars in E3SM timeSeriesStatsMontlyOutput stream

Use missing_value_mask for E3SM timeSeriesStatsMontlyOutput stream

Add mask to all tracer state and tend variables

Add mask to gm vars

Add mask to few more vars

Add mask to tend vars

Add mask to layerThick vars

Add mask to diag vars

Add mask to inSitu* vars

Add masking to all non-restart time series stats streams

Turn masking on for ocn standalone stats

Add missing code to set threadNum before its use

ELM: Use intent(inout) to preserve values in early exit.

In PhotosynthesisMod.F90, early exits from brent and brent_PHS do
not assign gs_mol* outputs values. The desired behavior in those
cases is for the input values to be retained. Thus, this commit
changes intent(out) to intent(inout) for these variables.

The original code usually works because Fortran normally doesn't
overwrite these inputs. But if an option like crayftn's -hzero is
provided, intent(out) stack variables are 0-init'ed, causing a
bug.

For pm-cpu (and alvarez), adjust/add PE layouts.
Add PE layout for RRM test in e3sm_prod.
Minor adjustment to PE layout for compset A_WCYCL* with res ne30pg2_oECv3 on 7 nodes.
Add a large (L) layout for compset A_WCYCL* res ne30pg2_oECv3 on 58 nodes.
Add alvarez to same to match pm-cpu.
[bfb]

Updating to SCORPIO 1.4.1

Updating the version of SCORPIO submodule to 1.4.1

This release includes the following changes,

* Fix for a performance bug where the NetCDF header was
  reallocated multiple times, slowing down runs that
  create large output files
* Support for ADIOS BP5 file format
* Preliminary support (disabled by default) for HDF5
* Improved performance by reducing some collectives
* Support for Perlmutter/Crusher
* Minor bug fixes

On GCP only, add larger allowance for MEMLEAK detection in tests
and add a slightly smaller allowance for TPUT.

This fixes a fail for one test in e3sm_extra_coverage
ERP_Lm3.ne4_oQU240.F2010.gcp_gnu

For alvarez at NERSC, update modules to defaults.
Including using GNU 12.1

Add instance string to mpaso output / rst files

Add instance string to mpas-si output / rst files

Move setting of branch/hybrid init files inside of inst loop

Add ismip6_retreat calving option

Uses the retreat parameterization from Slater et al. (2019) that
sets retreat length as a function of subglacial discharge and ocean
thermal forcing. Compiles but does not yet run.

Use ismip6GroundedFaceMelt package for retreat parameterization

We need to activate the ismip6GroundedFaceMelt package when
config_calving = 'ismip6_retreat' in order to use
ismip6_2dThermalForcing and ismip6Runoff in the ISMIP6
retreat parameterization.

Add working ismip6_retreat parameterization

Add a functional ismip6 retreat parameterization, active when
config_calving = 'ismip6_retreat'. This uses ismip6Runoff and
ismip6_2dThermalForcing fields from the two most recent time-levels
to calculate a retreat rate based on the parameterization of Slater
et al. (2019) [https://doi.org/10.5194/tc-13-2489-2019]. Currently,
ismip6Runoff and ismip6_2dThermalForcing must have the same xtime
values and be read in a stream called 'ismip-gis'.

Remove check for ringing alarm in mpas_li_calving

Remove check for ringing alarm in mpas_li_calving.F when looking for
ismip6-gis stream. This was not necessary and made for more confusing
code structure.

Throw error if forcing times are the same

Add a check that ismip6_retreat is not using two forcings with the same
timestamp. This could result either from an error in creating the input
files or from faulty logic in the routine itself.

Don't read ismip6-gis stream every timestep

Only read ismip6-gis stream when mostRecentAccessTime advances, i.e.,
when clock passes a new input time. Still force a read of two
time-levels on init. This still uses saveActualWhen = .true. to get
forcing times, which will be removed in the near future.

Do not set mostRecentAccessTime after init

Remove unnecessary reset of mostRecentAccessTime after forcing read of
first two input times from forcing files.

Remove redundant read of ismip6-gis stream on init

Do not force redundant read of ismip6-gis to get forcings and
mostRecentAccessTime at init time.

Add TODO comments where updates will be needed in coupled E3SM runs

This code currently assumes units from external ISMIP6 forcings.
These conversions may need to change when we use fields passed from
the E3SM coupler.

Apply suggestions from code review

Fix variable descriptions and amend a log message.

Co-authored-by: Matt Hoffman <wyeast@gmail.com>

Check configs and streams for ismip6_retreat

Check that ismip6-gis stream is found, config_ismip6_retreat_k is
negative, and that config_front_mass_bal_grounded is none, or else throw
an error. Also add error checks after calls to mpas timekeeping
routines, update comments, and expand on routine description.

Remove unnecessary variables from input streams

ismip6_2dThermalForcingPrevious, ismip6_2dThermalForcingCurrent,
ismip6RunoffPrevious, and ismip6RunoffCurrent are not needed in
input streams in Registry.xml. Remove them.

Make ismip6-gis an immutable stream

This ensures that careless setup does not lead to clobbering other
input variables like sfcMassBal that might inadvertently be in the same
input file. Also fix typo that looked for 'ismip-gis' instead of
'ismip6-gis' stream.

Fix error in ismip6_retreat on restart

Add forcingTimeStamp to Registry and add to restart stream to
prevent ismip6_retreat from throwing an error on restart.

Add halo update on calvingVelocity

Testing indicates that a halo update is needed either on submergedArea
or on calvingVelocity before calling li_apply_front_ablation_velocity.
Those two halo updates give the same result, so updating calvingVelocity
because submergedArea is a local variable. Alternatively, submergedArea
could be added to Registry, but this is simpler. Behavior is still
different across different numbers of processors, even without velocity
solver or advection. This is likely due to an issue within li_apply_front_ablation_velocity,
and will be addressed in another PR.

Add reference_time to ismip6-gis stream

Use 1950-07-01_00:00:00 as reference time in immutable stream
ismip6-gis. This is the first time in the ISMIP6 forcings files.
Also clean up several calls to mpas_log_write.

Fix issues with restarts and immutable stream

Use timestepNumber to prevent resetting forcingTime with
forcingTimeStamp every timestep on a restart. Also use (stream_cursor %
valid) to throw an error if immutable stream 'ismip6-gis' does not
exist in streams.landice.

Include a mask for the applied SMB

Include a config_smbMask option that zeros out sfcMassBalApplied in bare land

Change the condition statement for checking ice-free ocean cells

Address PR review comments

which include
1) revert the change made in the commit 647ea7ad34c50aabf343c06dea9bd4b956a3dca0
  and change the condition statement for checking bare-land cells

2  change the name of the config option more descriptive and not camel-cased

Add spatially variable von Mises threshold stress

Add the ability to have the von Mises threshold stress for grounded
and floating ice read in from input files. Add config options
`config_grounded_von_Mises_threshold_stress_source` and
`config_floating_von_Mises_threshold_stress_source` which are set  to
'data' if the spatially varying field is to be red from an input file,
and to 'scalar' if it is to be set by config_grounded_von_Mises_threshold_stress
and config_grounded_von_Mises_threshold_stress.

Require valid von Mises threshold stress values from input streams

Add check that groundedVonMisesThresholdStress and floatingVonMisesThresholdStress
are >0.0 and >=0.0, respectively if read from an input stream. Also fix minor bug
that was passing err instead of err_tmp to li_apply_front_ablation_velocity.

Add halo update for calvingVelocity in von_Mises_calving

Results from von Mises calving were not bit-for-bit across
different numbers of processors. This fixes that issue. It is not
apparent why this is necessary here, and why it is not necessary
for eigencalving, so perhaps the issue is rooted in the fields
used to calculate calvingVelocity in the von_Mises_calving routine.

Turn threshold damage calving on or off for von Mises

Adds flexibility when combining von Mises calving and damage calving.
If config_calving = 'von_Mises_stress', config_use_damage_calving_multiplier = .false.
and config_damage_calving_method = 'threshold',
ice with damage > config_damage_threshold will be removed after
von Mises calving rate is applied in each timestep.

Turn damage threshold calving on or off for eigencalving

As for previous commit for von Mises calving, you can now use
config_damage_calving_method to determine whether threshold calving will be
applied along with eigencalving. If config_damage_calving_method ==
'threshold', then damage threshold calving is applied after the calving
rate from eigencalving is applied.

Add 'none' option for config_damage_calving_method

Add 'none' option for config_damage_calving_method, to be used when using
config_calving = 'von_Mises_stress' or config_calving = 'eigencalving' without
damage threshold calving. Also add log errors if config_calving = 'von_Mises_stress'
or config_calving = 'eigencalving' and config_damage_calving_method = 'calving_rate'
or anything other than 'threshold' and 'none'. Update description of config_damage_calving_method
in Registry.xml to explain the usage options.

Update config_calving description

Update config_calving description to note that 'von_Mises_stress' and
'eigencalving' options can be used in combination with damage threshold
calving, and refer to description of config_damage_calving_method for details.

Add submodule for Sea Level Model

At this point this points to a private repo and will need to be moved to
a public repo before merging.

Makefile changes to compile with Sea Level Model

This commit modifies the Makefile to support compiling with or without
the sea level model.  To include it, add "SLM=true" on the make call.
If that is included, the build system expects the SLM submodule to have been
initialized/updated already.  It also activates a cpp directive that
includes calls to the sea level code.

Add missing call to bedtopo init routine

Add missing includes and dependency to Makefile for SLM

This was an oversight in the previous Makefile commit

Add example call to sea level model

SeaLevelModel initializtion is called from mpas_li_bedtopo.F

Working on the mpas_li_bedtopo.F code - not much difference from the previous commit

Sea level model is called from mpas_li_bedtopo.F

Include mpas alarm for MALI-SLM coupling time interval

	modified:   src/Registry.xml
                    added config_uplift_method = 'sealevelmodel'
                    added config_slm_coupling_interval

	modified:   src/mode_forward/mpas_li_bedtopo.F
                    added a check and reset alarm for coupling time interval

	modified:   src/mode_forward/mpas_li_core.F
                    added a setup and add mpas alarm for coupling interval

Temporary commit to mpas_li_core.F

Include a variable 'slmTimeStep' as a sea-level model timestep.

Since MALI timstep is finer than SLM timestep, the slmTimeStep gets updated
by a counter value 1 when MPAS alarm for slmcouplinginterval is rining.

Add config options for map files between MPAS mesh and SLM grid

- Config options added to Registry.xml
- Clean up unnecessary comments in  the other codes.

Implement interpolation between MALI mesh and SLM grid

Interpolation subroutine is copied from MPAS-O
(https://github.com/E3SM-Project/E3SM/pull/4472/commits/3332ea1e02ecb03352c9f2b7f7f92e14dc4b8392)

Fix the inclusion of a wrong variable

Include MPI functions.

Works fine on a single processor but not on multiple processor yet

Add a new variable 'UpliftDiff' in Registry and li_bedtopo module

'UpliftDiff' represents the amount of topography change
calculated by the sea-level model over the  sea-level model timestep.

Coupled simulation works fine for running on single processor only.

Organize the wrapper subroutine for slmodel_init

Include the subroutine 'interpoate_init' and calls to MPI subroutines
inside the wrappter 'slmodel_init'.

 Organize the SLM wrapper subroutine 'slmodel_solve'

MALI runs on multiple processors using MPI

Clean up the previous version. But somehow doesn't run on multiple processor anymore.

Add MALI mesh mask for SLM to use in merging MALI data with Global dataset (Works on multi-processors)

Additionally, change variable names 'globalArray#' and 'gatheredArray#' for
variables thickness, BedTopography and sea-level change to be more specific.

Clean up the code, fix indentation and apply minor edits such as variable name change

Organize the code in a more logical and efficient way

This commit resolves reveiew comments for PR #21 ('Add a regional
sea-level projection capability to MALI')

- Set up MPI variables outside of routine 'interpolation_init'.
- Move reformatting command of 1D SL array to 2D array inside of
  routine 'interpolation'
- Include #ifdef USE_SEALEVELMODEL inside of the check for config_uplift_method
- other small fixes

Include a wrapper subroutine that reads in the SLM namelist

Also fix white trailing spaces and variable names to use Camel-case

Clean up and reorginize the code

	modified:   src/Registry.xml
- Add a note on config_slm_coupling_interval
	modified:   src/mode_forward/mpas_li_bedtopo.F
- Move reformatting command of 1D SL array to 2D SL array outside of
 routine 'interpolation'
- use err_tmp instead of err in calling subroutines
	modified:   src/mode_forward/mpas_li_core.F
- delete statements that are not used in the routines in the module

Change 'err' to 'err_tmp' in calling routines

Update sea-level model submodule

Change the variable name 'UpliftDiff' to 'TopoChange'

Associated change is the change in the sign convention when
updating bedTopography.

Clean up

Update sea-level model submodule

Add an error statement for when compile flag and namelist config have a conflict

Log an error when MALI is built with 'SLM=true' but 'config_uplift_method'
namelist option is not set to 'sealevelmodel'.

Change a variable name to camelCase

deallocate MPI variables across routines

Also,
- fix add checks for err_tmp
- fix variable names to camelCase

Remove white trailing spaces

Fix white trailing spaces

Fix white trailing spaces

Update SLM submodule URL and hash

We recently transitioned the SLM repo to
https://github.com/MALI-Dev/1DSeaLevelModel_FWTW instead of using a
version in Holly's fork.  This commit updates the submodule URL and also
updates the commit used in the submodule.

To update an existing clone, you need to do the following from the E3SM
repo (not the submodule):
git submodule sync --recursive
git submodule update --init --recursive
Any new clones just need the second command, like usual.

Clean up and add an address to the SLM Github repo

Update the SeaLevelModel submodule

Change variable name `unit_num` in MALI to avoid conflict with the one in SLM

Update the SeaLevelModel submodule

Fix the file name for the 'gridToMpasFile'

Also declare the real value variable `rhs` with (kind=RKIND)

Fix the while loop in subroutine 'interpolate'

This fix avoids the bound-check error when the code is compiled in DEBUG mode

Modify the src/Makefile such that 'make clean' removes files in the SLM submodule

Add a check for consistency in coupling interval prescribed in MALI and SLM settings

In MALI, coupled interval is setup in namelist.landice under
'config_slm_coupling_interval'. In SLM, it is defined in
namelist.sealevel under the variable 'dtime'. These two variables
should have the same value when running coupled simulations.
This fix checks the consistency between the two values, and gives
an error if the check fails.

Change the variable name 'topoChange' to 'bedTopographyChange'

because 'topoChange' is not specific enough (e.g., ice sheet surface
change is also topography)

Initialize forcingTimeOld on restart

This prevents a fatal error in timekeeping on (now) line 2008 when doing
a restart that does not align with a time level in the forcing file.
While this prevents a fatal error, the restart is not BFB.

Force ISMIP6-retreat i/o on restart to initialize TF properly

This commit replace the fix of the previous commit with a solution that
also solves the restart BFB error for the ISMIP6 retreat param. when the
restart time does not match a time in the ISMIP6 forcing file.  It does
so by always forcing the previous forcing data to be read on the initial
timestep of a restart run.  This would be a redundant read if the
initial time of the restart happens to be a time in the forcing file,
but is necessary otherwise.  Rather than add logic to figure that out, I
just have that extra read always happen on the initial time of a restart
run.  That has the benefit of making the initial timestep of a restart
and a cold start use the same logic in the ISMIP6 retreat routine.

Before this commit, the new
landice/humboldt/mesh-3km_restart_test/velo-none_calving-ismip6_retreat
test failed the validation step.  After this commit, that test passes.

Fix indexing error causing decomp test to fail for ISMIP6-retreat param

Presumably this also led to inaccurate results.

The new test
landice_humboldt_mesh-3km_decomposition_test_velo-none_calving-ismip6_retreat
failed prior to this commit but passes after.

Move mask update prior to calling any calving routines

This commit has no effect on functionality but reduces code redundancy
by calling for mask and geometry update prior to calling any of the
calving routines.  It deletes the mask updates in individual calving
routines as they are no longer needed.  There are no changes to calving
test results with this change, indicating that the mask update is not
necessary, but it is safest to include it to keep the calving module
indpendent of potential changes made elsewhere in the code.

Replace iCell with indexToCellID(iCell) in log messages

Replace iCell with indexToCellID(iCell) in log messages, since iCell
is not useful for debugging.

Add applied basal mass balance fields

Add applied basal mass balance fields, in which the calculated and prescribed
basal mass balance is limited by the ice thickness.

Use applied basal mass balance in global and regional stats

Use applied basal mass balance fields in global and regional stats and remove
grounded, floating, and ice masks from calculations.

Move applied basal mass balance calculations to apply_mass_balance subroutine

Calculate basalMassBalApplied, groundedBasalMassBalApplied, and
floatingBasalMassBalApplied within the apply_mass_balance subroutine
in mpas_li_advection.F. This catches the case in which sfcMassBal
helps to remove all ice from a cell. Here sfcMassBal arbitrarily
removes ice first, but this is unlikely to be a problem.

Fix typo in comment in mpas_li_advection.F

Co-authored-by: Matt Hoffman <wyeast@gmail.com>

Distribute unablatedVolumeDynCell among neighboring cells

If a cell is left with unablated volume after step 2c in li_apply_front_ablation_velocity,
distribute that volume among its dynamic neighbors. Testing with von Mises calving on the MISMIP+
domain indicates that this greatly reduces inaccuracies in the calving routine, but that very large
calving events can still trigger the >10% calving inaccuracy error. This could potentially
be alleviated by adding a do-while loop over the new section 3 in li_apply_front_ablation_velocity,
but we may want to control for extremely large calving events anyway.

Require unablatedVolumeDynCell(iCell))>0.0_RKIND when applying ablation in step 3

Require unablatedVolumeDynCell(iCell))>0.0_RKIND when applying ablation in step 3
in li_apply_front_ablation_velocity to avoid doing extra calculations where they
do not apply.

Use small tolerance value for cellVolume

Use a tolerance consistent with round-off level ice thickness in
li_apply_front_ablation_velocity instead of requiring cellVolume > 0.0.

Add config_distribute_unablatedVolumeDynCell option

Add a namelist option to turn distribution of unablatedVolumeDynCell
to neighboring dynamic cells on or off. Default is to have this feature
enabled, but it may be desirable to have it disabled to ensure that results
are not strongly dependent on this ad hoc cleanup, and that the primary
control on the accuracy of calving results comes from the length of
the timestep (i.e., that the code is not violating a CFL-like
condition based on ablation velocity).

Create a calving CFL limit and apply to adaptive timestepper

There is reason to think it is incorrect to allow calving to remove more
than a whole grid cell in a given timestep - a CFL condition for
calving.  This commit calculates a calving CFL and adds it to the
adaptive timestepper.  It's important to note that the calving CFL is
calculated from the *previous* timestep, because determining the calving
CFL prior to starting the timestep is prohibitively difficult.  That
would require doing advection and other operations to generate the
geometry that the calvingVelocity will be based on just to calculate the
CFL limit, and then going back and actually progressing through the time
step for real once the dt has been set.  So this addition works under
the principle that the calving CFL limit will not change a lot from
timestep to timestep, thus a lagged CFL limit is a useful way to adapt
the timestep.  In initial tests this appears to work for the simple
MISMIP+ test case.

Note that this additional work is required beyond this commit:
* Right now it is only calculated if using the von Mises calving law
* It is not calculated on the initial time
* Restarts are not handled
* There is not yet an option to disable this
* Adjustments may be desired for how the option for what fraction of the
calving CFL to use interacts with the general CFL dt fraction option
* some general clean up may be desired

Add calving CFL calculation to other calving laws

This commit adds the calving CFL calculation to all the other calving
laws that use a calvingVelocity: specified_velocity, eigencalving,
ISMIP6-retreat, and damage-rate
This calculation should be low cost, so I've enabled it in these
routines, regardless of if it is being used.  I anticipate its result
will be a useful diagnostic even if it is not being applied in the
timestepper.

Improve logic for calving CFL application

The calving CFL application can be now controlled with the option
config_adaptive_timestep_include_calving.  If that option is true, the
calculation is only applied if a calving law that calculates
calvingVelocity is actually enabled, to avoid attempting to apply the
calvingCFLDt variable when it is not being calculated.

This commit also changes config_adaptive_timestep_CFL_fraction so that
it is only applied to the advective CFL.  That is typically how it has
been used.  I did not change the option name to avoid breaking backwards
compatibility of namelist files.  With this change, there is no longer a
master adjustment - there are only adjustments to the individual
candidate CFL conditions.  Note that the diffusive CFL does not have a
fraction, but it is rarely used, so that is fine for now.

Move config_min_adaptive_timestep check before force interval adjustment

This commit fixes a longstanding annoyance where an error would be
triggered if the adaptive timestepper’s force interval calculation
resulted in a timestep less than the limit specified by
config_min_adaptive_timestep

Add calvingCFLdt to restarts

This is needed for the calving CFL to behave correctly on a restart.
This variable is only needed in a restart file if the calving CFL is
actually being applied, but including it always because it is only a
scalar value, and adding logic for that would be complicated.

Calculate calving on initial time when not a restart

This is needed to get the calving CFL to work properly on the first
timestep because the calving CFL calculation is made on the timestep
prior to its application.  To ensure there is no actual calving for
rate-based calving laws, the model deltat is temporarily set to 0 before
calling the calving routines.  Note that non-rate-based calving laws
(thickness, damage thresholds, etc.) will now be applied on the
initial time whereas before they were not.  But I think that behavior is
actually preferable.  We don't have any tests or routine configurations
using those calving laws currently, so that change will not disrupt
anything.

Initialize layerThickness

Required for von Mises calving to work on initial time

Adjust log output to make timesteps clearer

Improve CFL log messages

* Convert limiting timestep values reported to days
* Add a message saying which process controlled the timestep
* Add a new output variable for which process controlled the timestep so
that behavior could be analyzed more easily in post-processing if
desired.

Add variable and calculation for calving CFL ratio

This assesses if using the calving CFL from the previous timestep is
sufficient.

Make calving error threshold a nl option

Treat sea level more carefully, add some cleanup

These changes have only small impacts on the behavior of the calving
routine, but add some minor cleanup.

Update default nl values

Restore application of config_max_adaptive_timestep

It accidentally got dropped in an earlier commit.

Add addl scalars to globalStats

These aren't always relevant, but there is no harm in including them,
and it is convenient not to have to add them manually when desired.

Updates to Registry

The change to principalStrainRateRatio is unrelated to this branch, but
I happened to notice the missing Time dimension.

Minor cleanup from review

Fix bug in regional GL migration calculation

The regional mask was not applied, so the global GL migration flux was
being calculated for every region.

Add regionCellMasks to restart file

This is needed to allow regionalStats to restart correctly.

MALI: writing additional fields to ascii files

Check for dynamic ice that is at the end of non-dynamic peninsula

Do the flood-fill in two steps in remove_icebergs. First, fill from
grounded ice into the dynamic ice shelf. Then, fill from the dynamic
shelf into the non-dynamic cells. This will eliminate the case in
which an iceberg of dynamic cells is connected to the rest of the
ice shelf by non-dynamic cells. This may be an overly expensive
solution to an issue that only happens once in a while, but it
should prevent the velocity solver from failing to converge.

Reduce big number in diffusivity calc to avoid timekeeping overflow

This only clearly appeared using Intel.

Reduce big number in mpas_li_advection as well

Avoid overflow by reducing bigNumber from 1e16 to 1e11.

Reduce big number in diffusivity calc to avoid timekeeping overflow

This only clearly appeared using Intel.

Allocate arrays on all procs before using MPI_GATHER(V)

Intel complained if destination arrays used on MPI_GATHER and
MPI_GATHERV were not allocated before use.  It was not obvious that it
was necessary to do so, because those arguments are unused on all but
the destination processor (where we already have them allocated to the
proper size).  To avoid unnecessary memory usage, they are simply
allocated to length 1 on all other processors.  The deallocate commands
are now performed on all processors.

Update SLM submodule to 26c156d

This brings in two critical bug fixes:
* initializing arrays that were previously used before being initialized
* fixing out of bounds error in the spharmt library that prevented usage
  of debug mode

Fix calving unablated volume error message

Add output vars for calved volume and uncalved volume

These were previously just written to the log file, but it is very
useful to have these as time series in a netCDF output file.

Note that totalRatebasedCalvedVolume differs from totalCalvingFlux,
because the former is just the result of the
li_apply_front_ablation_velocity routine, whereas the latter includes
calving from any calving process (damage, flow speed limit, iceberg
detection, etc.)

Note that the two variables are currently only being calculated for von
Mises calving and should be extended to the other rate-based calving
laws.

Add output variable for reporting if Albany converges or not

And add to the globalStats stream so it becomes a standard output
field.

Apply new calving quality metrics to the other rate-based calving laws

This commit applies the new variables totalRatebasedCalvedVolume and
totalRatebasedUncalvedVolume to the remaining rate-based calving laws.

Do not modify floatingBasalMassBal

Do not modify floatingBasalMassBal, which is an input variable. If
internal melting is added to groundedBasalMassBal in a floating cell,
add it to basalMassBal instead of floatingBasalMassBal.

Initialize thermalCellMask in mpas_li_advection.F

Initialize thermalCellMask in li_advection_thickness_tracers. The
previous commit neglected to initialize and deallocate this array, and
thus it was not being applied properly.

Add internalMeltRate directly to basalMassBal

Add internalMeltRate directly to basalMassBal instead of adding it
to groundedBasalMasBal and then repartitioning between grounded and floating.
Also define internalMeltRate in Registry.xml.

Rename internalMeltRate to drainedInternalMeltRate

The name internalMeltRate was misleading because there is a component
that is not accounted for in this variable, which is the the amount
needed to fill the maximum water fraction.

Add check of config_thermal_calculate_bmb to new drained melt BMB assignment

The assignment of drained melt to the BMB field previously was
conditional on a check to config_thermal_calculate_bmb, but that
dependence was lost with the reorganization in this branch.  This commit
re-applies that to the new assignment location.

Initialize variables that are used when MALI is coupled to the SLM

Fix check on vM stress values

The check needs to only cover 1:nCells, because the array values for the
nCells+1 garbage index are not valid.  Without this change, the code
would die when reading in a data field of vM values.

Fix halo bug in restore_calving

As originall written, the restore_calving routine only calculates over
owned cells but does not do a halo update, and a halo update never
happens over temperature.  This led to garbage values (potentially 0K)
being used by Albany along processor boundaries after restore calving
resets floating ice that has gone away.  This commit fixes this
bug by calculating restore_calving over all cells (owned+halo).  No
additional halo update is necessary.

Update grounding line mask definition for cells, edges, and vertices

Define grounding line for cells as grounded cells with at least one
dynamic floating neighbor, and for cells and vertices as those with
one grounded and one dynamic floating cell. Previously the grounding
line was defined as where a grounded cell was adjacent to a floating
(dynamic or non-dynamic) or open ocean cell. This definition is
problematic because it prevents us from closing a mass budget, since flux
at the grounding line in the previous definition was not necessarily
leading to loss of grounded ice.

Update logic for ablation velocity after redefining grounding line

Update logic for ablation velocity after redefining grounding line.
We will still use applyToGroundingLine to mean apply the ablation
velocity at the margin of grounded ice, whether there is an ice shelf
or not. However, the new definition of the grounding line required
a change from using li_mask_is_grounding_line to li_mask_is_grounded_ice
in one place.

Do not melt floating non-dynamic fringe

When zeroing melt on floating non-dynamic fringe, use
li_mask_is_grounded_ice instead of li_mask_is_grounding_line. Using
the grounding line mask here is not consistent with the new definition
of the grounding line, and leads to ice-shelf melt being applied to
the floating non-dynamic fringe, resulting in spuriously high mass loss.

Changes to hydro model to comply with new GL mask change

After various attempts, it was found that the best way to ensure BFB
behavior in the hydro model was to re-create the old GL mask (that
includes both grounded to floating and grounded open ocean) as a new
mask variable in the hydro model.  More work could be done to make a
more univerasl hydro mask and apply it uniformly, but this commit does
the minimum to keep the hydro model BFB using the existing logic in the
hydro model.

Change two instances of cycle to exit

One is new in this branch, the other was longstanding.  In both cases,
the old usage was not hurting anything, but was not actually doing
anything.  With this change, the code will stop looing over neighbors
when it finds what it's looking for, which will be slightly more
efficient.  COMPASS testing indicates no change in results, as should be
the case.

Update comments with month of grounding line mask update

Previously said the change occured in May, but it is being merged in Oct
2022.

Add new field fluxAcrossGroundingLineOnCells

This is the GL flux definition required for the ligroundf field required
by ISMIP6.  Calculating it within MALI allows us to calculate it fully
accurately and eliminates a slow post-processing calculation.

Update bld files with output from scripts to match Registry changes

Replace some where/elsewhere/elsewhere coding to make pgi happy

Clean up "*-nersc" build targets

Remove `pgi-nersc` and `cray-nersc`, which a…
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
CC PR is climate changing mpas-ocean non-BFB PR makes roundoff changes to answers.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants