Skip to content

Conversation

adcroft
Copy link
Member

@adcroft adcroft commented Nov 26, 2024

The default= optional argument to get_param() only provides a uniform value to initialize an array of reals. This commit adds the optional defaults= argument that must have the same length as the values argument.

I've also added a few instances of this optional argument:

  • by adding the initialize_thickness_param() procedure, selected by THICKNESS_CONFIG = "param". The procedure was based on the "uniform" method, and uses the parameter THICKNESS_INIT_VALUES which defaults to uniform values derived from MAXIMUM_DEPTH
  • the setting of MLD_EN_VALS in MOM_diabatic_driver.F90 which was previously using a work around to set defaults to 25, 2500, 250000 J/m2.
  • two vectors of 4 values in user/user_change_diffusivity.F90

For super repos there will be some doc file changes, but no answer changes.

This PR was inspired after realizing get_param_real_array() could not be used in PR #737 where 16 run-time non-dimensional parameters are being added.

Copy link

codecov bot commented Nov 26, 2024

Codecov Report

Attention: Patch coverage is 82.05128% with 7 lines in your changes missing coverage. Please review.

Project coverage is 36.69%. Comparing base (7a9adbc) to head (9eec351).

Files with missing lines Patch % Lines
src/framework/MOM_file_parser.F90 55.55% 2 Missing and 2 partials ⚠️
src/initialization/MOM_state_initialization.F90 88.00% 1 Missing and 2 partials ⚠️
Additional details and impacted files
@@             Coverage Diff              @@
##           dev/gfdl     #760      +/-   ##
============================================
+ Coverage     36.67%   36.69%   +0.01%     
============================================
  Files           278      278              
  Lines         84269    84306      +37     
  Branches      15869    15877       +8     
============================================
+ Hits          30908    30935      +27     
- Misses        47532    47537       +5     
- Partials       5829     5834       +5     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@adcroft
Copy link
Member Author

adcroft commented Nov 26, 2024

I see that tc2.target failed, but that's because the new parameter value isn't available with the target executable. @marshallward shouldn't we be running the target executable in the target's .testing?

@marshallward
Copy link
Member

marshallward commented Nov 26, 2024

I think regression uses everthing in its own .testing for compilation (that is, .testing/build/target_codebase/.testing) but they may use the top .testing directory when running the tests.

(Edit: We build target inside of its own .testing, but it is symlinked into .testing/build/target and treated as just another executable alongside the others. It would take some additional work to completely run the tests inside of its own .testing and run the regression checks. Running and checksum diffs are not set up for arbitrary directories.)

You are probably right and it should run its own tests. I am a bit worried that we are starting to put too much distance between the two codebases, but I can't see any specific problem here.

@adcroft adcroft added the enhancement New feature or request label Nov 26, 2024
Copy link
Member

@Hallberg-NOAA Hallberg-NOAA left a comment

Choose a reason for hiding this comment

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

This change looks like it will be a valuable new addition, but before this is accepted, I am requesting one minor change to keep this new commit in line with the overall MOM6 code style guide.

@adcroft
Copy link
Member Author

adcroft commented Dec 2, 2024

FWIW, dropping the new defaults= optional argument, and making default= a real array does work and only requires a few files to be corrected. These are the files where we use default= to set a vector of constant values, and all were easily changed:

        modified:   ALE/MOM_hybgen_regrid.F90
        modified:   core/MOM_open_boundary.F90
        modified:   initialization/MOM_state_initialization.F90
        modified:   parameterizations/vertical/MOM_diabatic_driver.F90
        modified:   tracer/oil_tracer.F90
        modified:   user/MOM_wave_interface.F90
        modified:   user/user_change_diffusivity.F90

I'll ask a the dev call if we prefer this version...

@adcroft adcroft force-pushed the add-defaults-vector-to-get_param_real_array branch from 9eec351 to 8f7dd3e Compare December 2, 2024 20:50
Removed two instances of `fail_if_missing=.true., default=0.` which
are contradictory: a default value is meaningless if the parameter must
be specified.

I encountered this when adding the `defaults=` option to `get_param_real_array()`.
@adcroft adcroft force-pushed the add-defaults-vector-to-get_param_real_array branch 3 times, most recently from 9f15344 to f8dda20 Compare December 2, 2024 22:14
The `default=` optional argument to get_param() only provides a uniform
value to initialize an array of reals. This commit adds the optional
`defaults=` argument that must have the same length as the `values`
argument.

I've also added a few instances of this optional argument:
 - by adding the `initialize_thickness_param()` procedure, selected by
   `THICKNESS_CONFIG = "param"`. The procedure was based on the "uniform"
   method, and uses the parameter `THICKNESS_INIT_VALUES` which defaults
   to uniform values derived from `MAXIMUM_DEPTH`
 - the setting of MLD_EN_VALS in MOM_diabatic_driver.F90 which was
   previously using a work around to set defaults to 25, 2500, 250000 J/m2.
 - two vectors of 4 values in user/user_change_diffusivity.F90

There will be some doc file changes, but no answer changes.
@adcroft
Copy link
Member Author

adcroft commented Dec 2, 2024

I tried assumed-rank arrays for default=, following @marshallward 's suggestion, and it looked very promising, but we discovered a feature-support problem with our favourite third compiler. 🙄 So I'm settling on the original proposal, have undone the tc2 change, and addressed the other raised issues.

Copy link
Member

@Hallberg-NOAA Hallberg-NOAA left a comment

Choose a reason for hiding this comment

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

These changes now all look good to me.

We probably should add the analogous new defaults optional argument to get_param_integer_array(), but that could wait for a subsequent PR.

@adcroft
Copy link
Member Author

adcroft commented Dec 5, 2024

@adcroft adcroft merged commit 0337147 into NOAA-GFDL:dev/gfdl Dec 5, 2024
10 checks passed
Hallberg-NOAA added a commit to Hallberg-NOAA/MOM6 that referenced this pull request Jan 1, 2025
  The `default=` optional argument to get_param() only provides a uniform value
to initialize an array of integers.  This commit adds an optional `defaults=`
argument to get_param_int_array(), doc_param_int_array() and
log_param_int_array() to allow for the specification of an array of default
values.  These additions are analogous to what had previously been added for
real arrays in github.com/NOAA-GFDL/pull/760.

  This commit also adds the new internal function int_array_string(), analogous
to real_array_string(), in MOM_document.  This differs slightly from its real
array counterpart in that it only uses the syntax like `3*75` for lists of
integers that are longer than we would use to specify dates and times or pairs
of layout parameters, because "(0, 0)" seems more readily interpretable than
"(2*0)".

  The new defaults argument is now used in the get_param calls for LAYOUT and
IO_LAYOUT, and in setting the tidal reference dates.

  Several spelling errors in comments were also corrected in the files that
were being edited.

  All answers are bitwise identical, but there are minor changes in many
MOM_parameter_doc.layout files and some MOM_parameter_doc.all files.
marshallward pushed a commit that referenced this pull request Jan 2, 2025
  The `default=` optional argument to get_param() only provides a uniform value
to initialize an array of integers.  This commit adds an optional `defaults=`
argument to get_param_int_array(), doc_param_int_array() and
log_param_int_array() to allow for the specification of an array of default
values.  These additions are analogous to what had previously been added for
real arrays in github.com//pull/760.

  This commit also adds the new internal function int_array_string(), analogous
to real_array_string(), in MOM_document.  This differs slightly from its real
array counterpart in that it only uses the syntax like `3*75` for lists of
integers that are longer than we would use to specify dates and times or pairs
of layout parameters, because "(0, 0)" seems more readily interpretable than
"(2*0)".

  The new defaults argument is now used in the get_param calls for LAYOUT and
IO_LAYOUT, and in setting the tidal reference dates.

  Several spelling errors in comments were also corrected in the files that
were being edited.

  All answers are bitwise identical, but there are minor changes in many
MOM_parameter_doc.layout files and some MOM_parameter_doc.all files.
@adcroft adcroft deleted the add-defaults-vector-to-get_param_real_array branch January 21, 2025 16:15
alex-huth pushed a commit to alex-huth/MOM6 that referenced this pull request Apr 14, 2025
  The `default=` optional argument to get_param() only provides a uniform value
to initialize an array of integers.  This commit adds an optional `defaults=`
argument to get_param_int_array(), doc_param_int_array() and
log_param_int_array() to allow for the specification of an array of default
values.  These additions are analogous to what had previously been added for
real arrays in github.com/NOAA-GFDL/pull/760.

  This commit also adds the new internal function int_array_string(), analogous
to real_array_string(), in MOM_document.  This differs slightly from its real
array counterpart in that it only uses the syntax like `3*75` for lists of
integers that are longer than we would use to specify dates and times or pairs
of layout parameters, because "(0, 0)" seems more readily interpretable than
"(2*0)".

  The new defaults argument is now used in the get_param calls for LAYOUT and
IO_LAYOUT, and in setting the tidal reference dates.

  Several spelling errors in comments were also corrected in the files that
were being edited.

  All answers are bitwise identical, but there are minor changes in many
MOM_parameter_doc.layout files and some MOM_parameter_doc.all files.
iangrooms added a commit to iangrooms/MOM6 that referenced this pull request Jun 17, 2025
Fix a bug with subroutine write_energy when using a DT<2. Otherwise,
the energy outputs are written at wrong time steps.

The reason was that time type divide is essentially a floor.
So DT/2 = 0 if DT<2.
Remove extra copy of compute_global_grid_integrals

  The subroutine compute_global_grid_integrals appeared in both the
MOM_state_initialization and MOM_shared_initialization modules, but was only
being called from the latter.   This commit removes the extra copy in
MOM_state_initialization.  It also removes some unnecessary parentheses in the
copy that is being retained, in part to facilitate the review of this commit.
All answers are bitwise identical, and no publicly visible interfaces are
altered.

Make REMAPPING_USE_OM4_SUBCELLS the default

In addition to REMAPPING_USE_OM4_SUBCELLS, for ALE remapping, there are
several parameters of the form XXX_REMAPPING_USE_OM4_SUBCELLS, where
XXX identifies the target, and they all currently default to True.

To simplify setting them all to False, which is recommended, the defaults
for the XXX versions is changed to the value of REMAPPING_USE_OM4_SUBCELLS.

Answers are only changed if REMAPPING_USE_OM4_SUBCELLS is set to False
and the default (now False) is used for one or more of the other parameters.
In such cases the original behaviour can be recovered by explicitly
setting the other parameters to True.

Removed default for mandatory time scale in OBCs

Removed two instances of `fail_if_missing=.true., default=0.` which
are contradictory: a default value is meaningless if the parameter must
be specified.

I encountered this when adding the `defaults=` option to `get_param_real_array()`.

Adds a vector of default values to get_param_real_array()

The `default=` optional argument to get_param() only provides a uniform
value to initialize an array of reals. This commit adds the optional
`defaults=` argument that must have the same length as the `values`
argument.

I've also added a few instances of this optional argument:
 - by adding the `initialize_thickness_param()` procedure, selected by
   `THICKNESS_CONFIG = "param"`. The procedure was based on the "uniform"
   method, and uses the parameter `THICKNESS_INIT_VALUES` which defaults
   to uniform values derived from `MAXIMUM_DEPTH`
 - the setting of MLD_EN_VALS in MOM_diabatic_driver.F90 which was
   previously using a work around to set defaults to 25, 2500, 250000 J/m2.
 - two vectors of 4 values in user/user_change_diffusivity.F90

There will be some doc file changes, but no answer changes.

FMS API: Convert real kind of constants

Two latent heat constants are imported directly from FMS, which is built
independently of MOM6.  Previously, it was a safe assumption that both
would be built with double precision, but this is no longer the case
since FMS now supports both single and double precision.  This could
cause conflicts with mixed-precision builds.

This patch converts the values from FMS-precision to MOM-precision.

Single->double should not affect reproducibility since every
single-precision number can be exactly represented in double precision.
Double->single could affect reproducibility, but this is not an issue
since MOM6 does not run in single precision.

Inline harmonic analysis (#744)

* Inline harmonic analysis

Important bug fix:
    1) The Cholesky decomposition was operating on entries below
the main diagonal of FtF, whereas in the accumulator of FtF, only
entries along and above the main diagonal were calculated. In this
revision, I modified HA_accum_FtF so that entries below the main
diagonal are accumulated instead.
    2) In the accumulator of FtSSH, the first entry for the mean
(zero frequency) is moved out of the loop over different tidal
constituents, so that it is not accumulated multiple times within
a single time step.

* Inline harmonic analysis

Another bug fix: initial state added back to the mean state.

* Inline harmonic analysis

Minor update to HA_solver
Tidal angular frequency has units [rad s-1] (#764)

* Tidal angular frequency has units [rad s-1]

Tidal frequencies are always angular frequencies to simplify applying
sine and cosine.  These have MKS units [rad s-1] but they are all
currently listed as [s-1].

Updated dOxygen comments for variables, e.g. [T-1 ~> s-1] becomes
[rad T-1 ~> rad s-1].  Updated get_param units. e.g. units="s-1"
becomes units="rad s-1".

No answers are changed, but the logged parameter units are different.

There are frequencies in MOM_internal_tides.F90 but these have not
been updated because they may be specified incorrectly.  They are
used as if they are [T-1] but they are calculated as 2PI/period [rad T-1].

  real, allocatable, dimension(:) :: frequency  !< The frequency of each band [T-1 ~> s-1].

  real    :: period             ! A tidal period read from namelist [T ~> s]

  ! The periods of the tidal constituents for internal tides raytracing
  call read_param(param_file, "TIDAL_PERIODS", periods)

  do fr=1,num_freq
    period = US%s_to_T*extract_real(periods, " ,", fr, 0.)
    CS%frequency(fr) = 8.0*atan(1.0)/period
  enddo

All MOM6-examples cases have INTERNAL_TIDES=False and so can't
resolve this issue.

* fixed too-long line
Refactor of vertical reconstruction adding six new schemes (#741)

add diagnostic for Kd_Work related to the Kd_add part (#765)

* add diagnostic for Kd_Work related to the Kd_add part

* also corrects some capitalization typos

* fix id init

---------

Co-authored-by: Raphael Dussin <Raphael.Dussin@noaa.gov>
Co-authored-by: Alistair Adcroft <adcroft@users.noreply.github.com>
+*Obsolete WIND_CONFIG = "SCM_ideal_hurr" (#770)

Changed the code to issue a fatal error message when WIND_CONFIG =
"SCM_ideal_hurr", with the message including instructions on how to recover
mathematically equivalent solutions, and eliminated the subroutine
SCM_idealized_hurricane_wind_forcing() from the Idealized_hurricane module.  The
ocean_only MOM_parameter_doc files have also been modified to reflect that
"SCM_ideal_hurr" is no longer a valid setting for WIND_CONFIG.   All answers are
bitwise identical in any cases that run, but some cases may fail during
initialization with instructions on how to fix them.
*Fix a bug when EPBL_ORIGINAL_PE_CALC is false

  Corrected a bug that causes ePBL column to set the wrong variable and then use
an uninitialized variable when EPBL_ORIGINAL_PE_CALC is set to false.  This bug
was present when the EPBL_ORIGINAL_PE_CALC was first added on Sept. 30, 2016,
but it was not detected because only the default case with EPBL_ORIGINAL_PE_CALC
= True appears to being used or tested.  Any runs that used this code with
debugging compile options would have trapped it immediately.  This will change
answers when EPBL_ORIGINAL_PE_CALC is false.

+Add EPBL_MLD_ITER_BUG runtime parameter

  Added the new runtime parameter EPBL_MLD_ITER_BUG that can be set to false to
correct buggy logic that gives the wrong bounds for the next iteration  when
USE_MLD_ITERATION is true and successive guesses increase by exactly
EPBL_MLD_TOLERANCE.  By default all answers are bitwise identical, but there is
a new runtime parameter in some MOM_parameter_doc files.

+Add run-time ability to debug ePBL sensitivities

  Added the ability to passively run ePBL_column twice in a diagnostic mode and
then provide diagnostics of the differences in the diffusivities and boundary
layer depths that are generated with the two options.  This is controlled by the
new runtime parameter EPBL_OPTIONS_DIFF, which is an integer that specifies
which options to change or 0 (the default) to disable this capability.
Associated with this are the new diagnostics ePBL_opt_diff_Kd_ePBL,
ePBL_opt_maxdiff_Kd_ePBL and ePBL_opt_diff_h_ML, which only are registered when
the differencing is enabled.  For now, only changes associated with the settings
of EPBL_ORIGINAL_PE_CALC and EPBL_ANSWER_DATE can be evaluated, but this list
will grow as new options are added.

  As a part of these changes, there were some other reforms to the way that
MOM_energetic_PBL handles diagnostics, with the 2-d and 3-d arrays changed from
allocatable arrays with enduring memory commitments in the energetic_PBL_CS type
into simple arrays in energetic_PBL, relying on the fact that compilers will be
smart enough to avoid actually allocating this memory when it is unused to avoid
expanding the overall memory requirements of MOM6.  A number of allocate and
deallocate calls were eliminated by these changes.

  In addition, explicit 'units=' descriptors were added to numerous
register_diag_field calls to help identify inconsistent conversion factors.

  Also, the diagnostic of the number of boundary layer depth iterations was
added to the ePBL_column_diags, which enabled the intent of the CS and G
arguments to ePBL_column to be changed to intent(in).  The unnecessary extra
logic associated with the use of the OBL_converged variable in ePBL_column was
eliminated, but the results are uneffected.

  Because the MOM_parameter_doc files were changing anyway, and because this
commit is only a diagnostic change, about 6 spelling errors were corrected in
ePBL parameter descriptions as a part of this commit.

  All answers are bitwise identical, but there is a new runtime parameter and
there may be new diagnostics depending on the choice of a non-default value for
that new parameter.

+Add and test find_Kd_from_PE_chg

  Added the new internal subroutine find_Kd_from_PE_chg inside of the
MOM_energetic_PBL module to directly calculate an increment in the diapycnal
diffusivity from an energy input.  This can be used when ePBL does not convert
released mean kinetic energy into turbulent kinetic energy (i.e., if
MKE_TO_TKE_EFFIC = 0.) and is more efficient than the more general iterative
approach.  To preserve old answers, this new option is only enabled for the
surface boundary layer when the new runtime parameter DIRECT_EPBL_MIXING_CALC is
set to true.  This new option can be tested passively by setting
EPBL_OPTIONS_DIFF to 3 in a run that uses ePBL.  By default all answers are
bitwise identical, but there is a new runtime parameter in some of the
MOM_parameter_doc files.

+(*)Add ePBL bottom boundary mixing option

  Add the option to do energetically consistent bottom boundary layer mixing
with the new routine ePBL_BBL_column.  ePBL_BBL_column is closely based on the
surface-focused ePBL mixing in ePBL_column, but without adding convective
instability driven mixing or mean-TKE driven mixing to avoid possible
double-counting.  This new option is enabled by setting the new runtime
parameter EPBL_BBL_EFFIC to be positive.

  If both EPBL_BBL_EFFIC and BBL_EFFIC are set to positive values, there is a
risk of double-counting, but this case is not being trapped for now.

  The changes include the addition of a new mandatory vertvisc_type argument to
the publicly visible routine energetic_PBL.

  When this new ePBL bottom boundary layer mixing option is enabled, there are
several new diagnostics available that are related to bottom boundary layer
mixing.  Several new checksum calls were also added with this new option when
DEBUG = True. The MOM_parameter_doc files are altered by the addition of two new
runtime parameters, and by the correction of several spelling errors in the
descriptions of other ePBL parameters.  By default, all answers are bitwise
identical.

+Add DECAY_ADJUSTED_BBL_TKE option for ePBL BBL

  Added the ability to modify the bottom boundary layer TKE budget to account
for an exponential decay of TKE away from the boundary and the fact that when
the diffusivity is increased at an interface, it causes an increased buoyancy
flux that varies linearly throughout a well mixed bottom boundary layer and
through the layer above.  This new capability is enabled by the new runtime
parameter DECAY_ADJUSTED_BBL_TKE and is implemented via the new internal
function exp_decay_TKE_adjust.

  In addition, this commit adds 9 bottom-boundary layer specific run-time
parameters that are analogous to parameters that set the properties of the
surface-driven mixing, and take their defaults from them, but can now be set
independently for the bottom boundary layer.  The inefficient option to use
bisection to estimate the bottom boundary layer depth was eliminated, as it
certain that it is not being used yet in any cases, and it should be deprecated
for the estimation of the surface boundary layer.

  By default, all answers are bitwise identical, but there are up to 10 new
runtime parameters that will appear in some MOM_parameter_doc files.

+Add optional unscale argument to reproducing_sum

  Added the new optional unscale argument to reproducing_sum() and
reproducing_sum_EFP().  When this is used, the resulting real sum is restored to
the same scaling as the input arrays, but when there is an EFP_type returned it
is the same as it would have been had the input array been rescaled before it
was passed in.  All answers are bitwise identical, but thre is a new optional
argument to a two publicly visible interfaces.

+Add RZL2_to_kg element to unit_scale_type

  Added the new element RZL2_to_kg to the unit_scale_type for convenience when
rescaling masses and other globally integrated variables.  All answers are
bitwise identical, but there is a new element in a transparent type.

Work in rescaled units in write_energy

  Revised write_energy() and accumulate_net_input()  to work more extensively in
dimensionally rescaled variables by using the new unscale arguments to the
reproducing_sum functions. As a result of these changes, 15 rescaling factors
were eliminated or moved toward the end of write_energy().  All answers are
bitwise identical.

Use reproducing_sum with unscale in 5 places

  Use reproducing_sum with unscale to calculate the global mass diagnostic,
globally integrated ocean surface area, offline tracer transport residuals and
in calculating the OBC inflow area in tidal_bay_set_OBC_data().  This change allows
for the elimination or replacement of 6 rescaling factors and one added instance
with a consistent conversion factor and diagnostic units occur on the same line.
All answers are bitwise identical.

Bug fix related to directional internal wave drag (#774)

Co-authored-by: Robert Hallberg <Robert.Hallberg@noaa.gov>
Fix indexing error in extract_surface_state

Correct a bug that G%gridLonT/G%gridLatT were not using global indexing
in outputing extreme surface information.

Minor change of SSH check in extract_surface_state

Previously extreme surface message is triggered when `sea_lev` is
smaller than or **equal** to ocean topography:

sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) - G%Z_ref

The equality is innocuous for a non-zero minimum thickness
(Angstrom/=0), but it can be problematic for true zero thickness.

Besides, there is the fourth criterion in this check:

sfc_state%sea_lev(i,j) + G%bathyT(i,j) + G%Z_ref < CS%bad_val_col_thick

This is supposed to allow a user-specified tolerance (default=0.0),
which should be a more stringent check when the tolerance is larger than
zero. The equality from the first check contradicts this logic.

Recover diagnostic "SSH_inst"

Fix a bug that sea surface height averaged over every dynamic time step
(SSH_inst) is not outputted correctly.

Update MOM_wave_interface.F90 (#784)

* Update MOM_wave_interface.F90

The index at the interface has been changed from I-1 ->i+1 and J-1 -> j+1, which helps fixed model error in 3D simulations while keeping halo_size to 1.

* Update MOM_wave_interface.F90

Corrected the index for thickness which fixed the issue of 3D wave coupled simulations while keeping the halo_size in thickness as 1
Updates to use EPBL_BBL_EFFIC
- The present code only tests BBL_EFFIC when deciding whether to set the bottom TKE and ustar, but this means they are all zero when EPBL_BBL_EFFIC is non-zero and BBL_EFFIC is zero.
- Adds logic to also check EPBL_BBL_EFFIC, thereby allowing non-zero ustar and TKE for EPBL_BBL_EFFIC>0.0
- Fixed BBL_TKE diagnostic in EPBL that was not populated.
- Will change answers when EPBL_BBL_EFFIC>0.0, but won't change answers in any of our existing configurations.

Move MOM_generic_tracer with no changes (#22)

* Move MOM_generic_tracer with no changes

Renames src/tracer/MOM_generic_tracer.F90
to config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90

* Remove MOM_generic_tracer

Git rm MOM_generic_tracer.F90

---------

Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea68.ncrc.gov>
Corrected the rescaling of 5 KPP diagnostics

  Added missing factors to the conversion arguments in the register_diag_field
calls for the diagnostics of the KPP non-local transport tendency and the net
surface tracer fluxes, and corrected the dimensional scaling that is being
applied to the KPP_Vt2 diagnostic as calculated in KPP_compute_BLD.  All
solutions are bitwise identical, but now output files with these 3 sets of KPP
diagnostics are invariant to dimensional rescaling.  This can be verified with
the visc.nc file generated by the single_column/KPP test case in MOM6-examples.
Whereas previously the diagnostics KPP_Vt2, KPP_QminusSW, KPP_netSalt,
KPP_NLT_dTdt and KPP_NLT_dSdt in that file would change when dimensional
rescaling was applied, now they do not.  No output is changed unless dimensional
rescaling is used.

+Add a vector of defaults to get_param_array_int

  The `default=` optional argument to get_param() only provides a uniform value
to initialize an array of integers.  This commit adds an optional `defaults=`
argument to get_param_int_array(), doc_param_int_array() and
log_param_int_array() to allow for the specification of an array of default
values.  These additions are analogous to what had previously been added for
real arrays in github.com/NOAA-GFDL/MOM6/pull/760.

  This commit also adds the new internal function int_array_string(), analogous
to real_array_string(), in MOM_document.  This differs slightly from its real
array counterpart in that it only uses the syntax like `3*75` for lists of
integers that are longer than we would use to specify dates and times or pairs
of layout parameters, because "(0, 0)" seems more readily interpretable than
"(2*0)".

  The new defaults argument is now used in the get_param calls for LAYOUT and
IO_LAYOUT, and in setting the tidal reference dates.

  Several spelling errors in comments were also corrected in the files that
were being edited.

  All answers are bitwise identical, but there are minor changes in many
MOM_parameter_doc.layout files and some MOM_parameter_doc.all files.

set descale in reproducing_sum_3d
Fix ALE_sponge tendency diagnostic units

  Corrected the units and conversion factor in init_ALE_sponge_diags for the
various sp_tendency_... diagnostics.  Previously they had only been correct for
the tendencies of nondimensional quantities.  The code also now stores the
scaling factor that is set in set_up_ALE_sponge_field_fixed for later use in
registering the sponge tendency diagnostics.  Several instances of unusual
spacing around semicolons in MOM_ALE_sponge were also standardized.  The
documented units and conversion factors for some diagnostics were corrected,
but all solutions are bitwise identical.

Refactor horizontally_average_field

  Refactored the horizontally_average_field() routine in MOM_diag_remap to work
in rescaled units by making use of the unscale arguments to the
reproducing_sum() routines.  A total of 9 rescaling variables were moved into
unscale arguments.  All answers and diagnostics are bitwise identical, and no
interfaces are changed.

+Refactor homogenize_field and revise its interface

  Refactored the homogenize_field routine in MOM_horizontal_regridding to make
use of the unscale argument to reproducing_sum(), and revised its interface to
make it more nearly consistent with the interface to homogenize_field_t() in
MOM_forcing_type.

  The interface changes include revising the order of the arguments, making the weight
argument options, replacing the scale argument with an optional tmp_scale
argument that is the inverse of the previous scale, and making the default for
the use of reproducing sums to be true when the answer_date argument is absent.
The two homogenize_field routines now give equivalent behavior when none of the
optional arguments to homogenize_field() are absent.  The homogenize_field calls
in MOM_temp_salt_initialize_from_Z() and the horiz_interp_and_extrap_tracer()
routines have been modified in accordance with the interface changes.

  All answers are bitwise identical, but the interface to a publicly visible
routine has been substantially changed to the point where any calls using the
previous interface will not compile.

Revise the ice_shelf dimensional rescaling

  Refactored the volume_above_floatation, write_ice_shelf_energy and
ice_shelf_solve_outer routines to work in rescaled units by making use of the
unscale arguments to reproducing_sum().  Also added or corrected comments
documenting the units of 11 real variables in these routines.

  The routine integrate_over_Ice_sheet_area was converted into a function and
var_scale was renamed to unscale for more consistency with the rest of the MOM6
code.  Additionally, add_shelf_flux and update_shelf_mass were modified to use
the scale arguments to time_interp_external.

  A total of 12 rescaling variables were eliminated or moved into unscale
arguments, and 2 blocks of code that scale input variables were eliminated.  All
answers and diagnostics are bitwise identical, and no interfaces are changed.

(*)Conversion arguments for ice shelf scalar diags

  Revised volume_above_floatation and integrate_over_ice_sheet_area to return
values in scaled units, and added conversion factors to the
register_scalar_field calls in MOM_ice_shelf so that all unit unscaling of
diagnostics occurs via  the diag mediator and not in the code itself.  After
this commit, all of the dimensional scaling factors for diagnostics in the
ice_shelf code occur via conversion factors that are immediately adjacent to the
declaration of the units of those diagnostics, facilitating comparison for
consistency.

  The declared units of one diagnostic, "taub_beta", were revised from "MPa s
m-1" to "MPa yr m-1" to be consistent with its conversion factor, which was also
corrected (essentially by a factor of [Z L-1 ~> 1].  Several missing unit
conversion factors for rates of mass change were also added, and about 15
missing or incorrect units in lines documenting real variables were added or
fixed.

  The input scale factor for the (perhaps unused) variable INPUT_VEL_ICE_SHELF
was corrected from `US%m_s_to_L_T*US%m_to_Z` to just `US%m_s_to_L_T`.

  With these changes, all solutions are bitwise identical, apart from regional
cases with a specified inflow when run with dimensional rescaling.  The ice
shelf diagnostics should now be invariant when dimensional rescaling is applied,
and the units of the ice shelf diagnostics are now all consistent with the
required rescaling factors.

Users set solo driver forcing time records per day

  Refactored the solo_driver version of MOM_surface_forcing to modify how the
time levels to read from forcing files, with a series of 11 new runtime
parameters (WIND_FILE_DAYS_PER_RECORD, etc.) giving the user to specify how many
days each time record spans (for positive values) or the number of time records
in the file per day (for negative values).  If these parameters are set to 31,
the month returned by the FMS `get_date()` function sets the time level.  If
they are set to 0 (the default) the previous behavior is obtained, in which the
time level is determined from the number of records in the file.  To always take
the first record in the file, set these parameters to large values (greater than
1000000 will always work, but in practice smaller values do the same thing).  By
default all answers are bitwise identical, but there are up to 10 new runtime
parameters in some MOM_parameter_doc files for cases that have `WIND_CONFIG =
"file"` or `BUOY_CONFIG = "file"`.

  The list of new runtime parameters is `WIND_FILE_DAYS_PER_RECORD`,
`LONGWAVE_FILE_DAYS_PER_RECORD`, `SHORTWAVE_FILE_DAYS_PER_RECORD`,
`EVAPORATION_FILE_DAYS_PER_RECORD`, `LATENTHEAT_FILE_DAYS_PER_RECORD`,
`SENSIBLEHEAT_FILE_DAYS_PER_RECORD`, `RAIN_FILE_DAYS_PER_RECORD`,
`SHORTWAVE_FILE_DAYS_PER_RECORD`, `RUNOFF_FILE_DAYS_PER_RECORD`,
`SSTRESTORE_FILE_DAYS_PER_RECORD` and  `SALINITYRESTORE_FILE_DAYS_PER_RECORD`.

  The problem identified in github.com/NOAA-GFDL/MOM6/issues/631 is addressed by
this commit, and that issue can be closed once this commit is merged onto the
main branch of MOM6.

Merge branch 'main' into dev/gfdl

Fix dimensional rescaling in predict_MEKE

  Refactored predict_MEKE() to include two missing dimensional rescaling factors
and correct a third and to clearly indicate the units and the nature of the
variables where they are being used, in coordination with Andrew Shao.  Some
white spacing around semicolons in the same MOM_MEKE file was also modified to
follow the patterns used elsewhere in the MOM6 code.  All answers are bitwise
identical in cases that do not use dimensional rescaling.

Refactor SAL in MOM_PressureForce_FV

Self-attraction and loading calculation in Boussinesq pressure gradient
force is refactored to be consistent with the algorithm in
non-Boussinesq version.

Add option for SAL to use bottom pressure anomaly

Using SSH to calculate self-attraction and loading (SAL) is only
accurate for barotropic flows. Bottom pressure anomaly should really be
used for general purposes.

* New runtime parameter
A runtime parameter SAL_USE_BPA is added to use bottom pressure anomaly.
The option requires an input file for long term mean reference bottom
pressure. The reference bottom pressure field is stored with SAL_CS.

* Refactor SAL and tides in Boussinesq mode
As the total bottom pressure is needed for bottom pressure anomaly, SAL
calculation in Boussinesq mode needs to be refactored. In addition,
there is a longstanding bug in Boussinesq mode that interface height is
modified by SAL and tides, and the modified interface height is
erroneously used for density and pressure calculation later on.
Therefore the SAL and tides are refactored by moving their calculations
after pressure is calculated. Tide answers before 20230630 is retained.
Answers after 20230630 are changed for Boussinesq mode.

Rename variables with tides in Pressure Force

* Local variable SSH is renamed to pbot_anom.
* Tides related local variable names are changed to be less confusing.
Notably, `e_sal_tide` is renamed to `e_sal_and_tide` (the summation of
sal and tide), not to be confused with `e_tide_sal`, which is renamed to
`e_tidal_sal` (sal from tides).
* Fix e_tidal diagnostics in Boussinesq mode. The term was not assigned
if tide answer date is >20230630.

Alternative method for SAL and tides in Boussinesq

Add an alternative method for SAL and tides in Boussinesq mode. The
current method adjusts interface heights with geopotential height
anomaly for SAL and tides. For non-Boussinesq mode, the current method
is algebraically the same as taking the gradient of SAL and tide
geopotential (body forcing). For Boussinesq mode, there is a baroclinic
component of tidal forcing and SAL. The alternative method is added to
calculate the gradient of tidal forcing and SAL directly at the cost of
additional multiplications. The new method is controlled by runtime
parameter BOUSSINESQ_SAL_TIDES.

Remove rescaling bottom pressure in SAL

* Instead of rescaling bottom pressure to height unit,
calc_Loving_scaling is modified to be conditionally dimensional. When
calculating self-attraction and loading, Love numbers are now
dimensional when bottom pressure anomaly is used as an input. This
change eliminate Love numbers' dependence on mean seawater density.

* A new coefficient called linear_scaling is added to SAL CS for the
same purpose, although to use bottom pressure anomaly for scalar
approximation is not quite justifiable. A WARNING is given when users
try to do that.

* Two separate field, SSH (sea surface height anomaly) and pbot (total
bottom pressure), are used for calculating self attraction and loading
when SAL_USE_BPA = False and SAL_USE_BPA = True, respectively. The
change is to enhance code readability.

Add answer date flag for tide/SAL

A new date for TIDES_ANSWER_DATE is added to recover answers for tides
in Boussinesq mode before 20250201.

Refactor Love_scaling calculation in SAL module

* Precalcualte a local field `coef_rhoE` to avoid in-loop division and
if-blocks. The unit of coef_rhoE depends on use_bpa.
* Fix a few unit description typos in SAL module and two other files.

+Refactor MOM_opacity to replace hard-coded params

  Refactored MOM_opacity to replace hard-coded dimensional parameters for the
Manizza and Morel opacity fits with run-time parameters, and also added the
runtime parameter OPACITY_BAND_WAVELENGTHS to provide the ability to set the
wavelengths of the bands, even though these are not actually used in MOM6.  By
default, these parameters are all set to the previous hard-coded values, using
the recently added defaults argument to get_param_real_array().  The bounds of
the frequency band label arrays with the MANIZZA_05 opacity scheme were also
corrected when PEN_SW_NBANDS > 3, but it would not be typical to use so many
bands for no purpose and these labeling arrays (optics%min_wavelength_band and
optics%max_wavelength_band) do not appear to be used anywhere.  In addition, the
unused publicly visible routines opacity_manizza and opacity_morel were
eliminated or made private.  All answers are bitwise identical, but there are
new entries in some MOM_parameter_doc files.

implement spatially varying decay rates for internal tides leakage (#794)

* add option to specify horizontally varying decay

* extend to vertical modes

* fix comments

* corrected description of decay_rate array

---------

Co-authored-by: Raphael Dussin <Raphael.Dussin@noaa.gov>
Add ice shelf pressure initialisation bug fix (#800)

* Add ice shelf pressure initialisation bug fix

This commit adds a new parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX
and associated bug fix. This bug occurs when TRIM_IC_FOR_P_SURF
pressure initialisation is used for ice shelf, whilst in
Boussinesq mode and a nonlinear equation of state. The subroutine
trim_for_ice calls cut_off_column_top. This routine in Boussinesq
mode calls find_depth_of_pressure_in_cell in
MOM_density_integrals.F90. This subsequently calls the function
frac_dp_at_pos which calls the density equation of state based on
given salinity, temperature and depth, which is incorrectly converted
into pressure by z position (which is negative) multiplied by g and
rho0. The bug results in incorrect densities from the equation of state
and therefore an imperfect initialisation of layer thicknesses and sea
surface height due to the displacement of water by the ice.

The bug is fixed by multiplying the pressure by -1.

This commit introduces parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX with
default value False to preserve answers. If true, the ice shelf
initialisation is accurate to a high degree with a nonlinear
equation of state.

Answers should not change, except for the added parameter in
MOM_parameter_doc.

The same functions are called by a unit test in MOM_state_initialization;
in this commit I set the MOM_state_initialization unit test to
use the existing bug (with a false flag).

* Resolve indenting and white space inconsitencies with MOM6 style
for ice shelf pressure initialisation bug fix FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX
+Add the optional argument old_name to get_param

  Added the new optional argument old_name to the 8 get_param() routines. This
new capability allows for an archaic parameter name to be specified and for
appropriate warnings encouraging the user to migrate to using the new name
while still setting the parameter as intended, or error messages in the case of
inconsistent setting via the archaic name and the correct name.  The logging
inside of the MOM_parameter_doc files only uses the correct parameter name.

  Also added the new optional argument set to the 8 read_param routines, to
indicate whether a parameter has been found and successfully set.  The new set
argument is now being used in read_param() calls in obsolete_int(),
obsolete_real(), obsolete_char() and obsolete_logical().  Obsolete_logical() in
particular was substantially simplified by the use of this new argument, and is
now only about half as long as it was.  The read_param() set argument is also
used in all of the get_param() routines when they are given an old_name
argument.

  The new old_name argument to get_param() is not yet being used in the version
of the MOM6 code that is being checked in, but it has been tested extensively
by adding or modifying get_param calls in a variant of the initialization code,
and it will be used in an updated version of github.com/NOAA-GFDL/MOM6/pull/725
to gracefully handle the deprecation of 4 parameter names.

  All answers are bitwise identical, but there are new optional arguments to
two widely used interfaces.

+Add grid_unit_to_L to the ocean_grid_type

  Add the new element grid_unit_to_L to the ocean_grid_type and the
dyn_horgrid_type, which can be used to convert the units of the geoLat and
geoLon variables to rescaled horizontal distance units ([L ~> m]) when they are
Cartesian coordinates.  When Cartesian coordinates are not in use,
G%grid_unit_to_L is set to 0.

  This new element of the grid type is used to test for inconsistent grids or to
eliminate rescaling variables in set_rotation_beta_plane(),
initialize_velocity_circular(), DOME_initialize_topography(),
DOME_initialize_sponges(), DOME_set_OBC_data(), ISOMIP_initialize_topography(),
idealized_hurricane_wind_forcing(), Kelvin_set_OBC_data(),
Rossby_front_initialize_velocity(), soliton_initialize_thickness(), and
soliton_initialize_velocity().   These are the instances where this new
variable could be used and bitwise identical answers are recovered.  There are
a few other places where they should be used, but where answers would change,
and these will be deferred to a subsequent commit.

  All answers are bitwise identical, but there are new elements in two
transparent and widely used types.

Update particles_run call to allow accurate particle advection using u, v or uh, vh

The drifters package is able to run in 2 different modes: one where the particles are advected using u and v, and one where they are advected using uh and vh. When u and v are used (i.e. the drifters are advected using the resolved velocities only, with no sub-grd component), the particles package needs the layer thickness that is consistent with these velocities, so particles_run needs to be called before any subgrid scale parameterizations are applied. This was not implemented correctly in my last PR, and is corrected here.

The particles package also needs the correct timestep for particle advection. This is added here.

+Remove dyn_horgrid_type%Rad_Earth

  Remove the unused element Rad_Earth from ocean_grid_type and dyn_horgrid_type.
The dimensionally rescaled equivalent element G%Rad_Earth_L is extensively used,
and it will continue to exist.  G%Rad_Earth_L was introduced in November 27,
2021 as a dimensionally rescaled version of G%Rad_Earth, while the unscaled
version was retained because at the time, the Rad_Earth element of the
dyn_horgrid_type is also used in SIS2.  However, SIS2 has not used G%Rad_Earth
since December 23, 2021, so after 3 years we can now safely remove this unused
element.

  Any cases on other branches that might be impacted by this change will not
compile.  All answers are bitwise identical, but there is one less element in
two transparent types.

Nodal modulation

This commit fixes a few (potential) inconsistencies between open
boundary tidal forcing and astronomical tidal forcing.

1. There was an inconsistency in the code that the nodal modulation
can be calculated in OBC tidal forcing but not in the astronomical
tidal forcing. This commit fixes this bug so that nodal modulation
is applied to both forcings.

2. In the previous version of MOM_open_boundary.F90, a different
set of tidal parameters can be set for open boundary tidal forcing,
leading to potential inconsistency with astronomical tidal forcing.
This commit obsoletes those for open boundary tidal forcing.

3. Another important bug fix is that the equilibrium phase is added
to the SAL term, which was missing in the original code.

Correct unit conversion for BS_coeff_h diagnostic

  Added missing conversion arguments for the register_diag_field calls for the
recently added diagnostics BS_coeff_h and BS_coeff_q.  All answers are bitwise
identical, but two diagnostics will have corrected dimensional rescaling when
EY24_EBT_BS is true.

Streaming Filter

The filters and their target frequencies are no longer hard-coded.
Instead, multiple filters with tidal or arbitrary frequencies as
their target frequencies can be turned on. The filter names are
specified in MOM_input and must consist of two letters/numbers. If
the name of a filter is the same as the name of a tidal constituent,
then the corresponding tidal frequency will be used as its target
frequency. Otherwise, the user must specify the target frequency by
"FILTER_${FILTER_NAME}_OMEGA" in MOM_input.

The restarting capability has also been implemented. Because the
filtering is a point-wise operation, all variables are considered
as fields, even if they are velocity components.

Merge branch 'dev/gfdl' into rescale_predict_MEKE
+Refactor the spatial mean calculations

  Refactored the spatial mean calculations in the functions global_area_mean(),
global_area_mean_u(), global_area_mean_v(), global_layer_mean(),
global_volume_mean() and adjust_area_mean_to_zero() to work in rescaled units by
making use of the unscale arguments to the reproducing_sum routines.  Comments
were also added to explain what the code does when both tmp_scale and unscale
arguments are provided, which effectively acts as to allow for changes in the
rescaled units.

  Global_area_integral() and global_mass_integral() were also similarly
refactored, but with the added difference that when a tmp_scale argument is
provided, the areas or masses in the integrals have units of [L2 A ~> m2 a] or
[R Z L2 A ~> kg a] instead of the mixed units of [m2 A ~> m2 a] or [kg A ~> kg
a].  As a result the code surrounding the 4 instances where global_area_integral
or global_mass_integral were being called with tmp_scale arguments (in MOM.F90
and MOM_ice_shelf.F90) also had to be modified in this same commit.  When the
optional var argument is absent from a call to global_mass_integral() so that it
is the mass itself that is returned, the values (if any) of scale, unscale
and tmp_scale are ignored, but the presence of tmp_scale determines whether the
mass is returned in [R Z L2 ~> kg] or [kg], consistent with the behavior that
would have been obtained if var had been an array of nondimensional 1s.

  This commit also includes a rescaling in the units of the areaT_global and
IareaT_global elements of the ocean_grid_type and dyn_horgrid_type to [L2 ~> m2]
and [L-2 ~> m-2], respectively.  Although the dyn_horgrid_type is shared between
MOM6 and SIS2, these elements are not used in SIS2.

  A total of 12 rescaling factors were eliminated or moved into unscale arguments
as a result of these changes.

  All answers are bitwise identical, but there are changes in the rescaled units
of two elements each in two transparent types, and changes to the rescaling
behavior of two publicly visible routines when they are called with tmp_scale
arguments.

Diagnose area integrated fluxes in rescaled units

  Keep 36 area integrated surface mass, heat or salt flux diagnostics in
forcing_diagnostics in rescaled units by replacing an unscale argument with a
tmp_scale argument to global_area_integral.  Also added the corresponding
conversion arguments to the register_scalar_field calls for each of these
diagnostics, so that the documented units and conversion factors can be used to
confirm the correctness of the documented units for these variables.  The
internal variables used for the integrated or averaged heat, mass or salt fluxes
were also separated for clarity about the units of these variables.  All answers
and diagnostics are bitwise identical.

Fix rescaling of internal tide of debugging code

  Corrected the internal tide rescaling arguments so that some of the debugging
variables have units that are consistent with their documented units.  This
involved changing the scale arguments to global_area_integral to tmp_scale
arguments in 8 places so that the units of the output retain the scaling of the
input variable.  The multiplication by the reverse of the scaling factor was
also added to 4 debugging output messages.

  The documented units of internal wave energies in 5 register_restart_field
calls were previously given as "m3 s-2" in non-Boussinesq mode and "J m-2" in
Boussinesq mode when the reverse was actually true.  This has now been
corrected.

  Also replaced the scale argument to 35 chksum calls with the equivalent but
preferred unscale argument, following the pattern elsewhere in the MOM6 code.

  All answers are bitwise identical, and only debugging code was modified.

Corrected many unit descriptions in comments

  Corrected the descriptions of the units of about 113 variables spread across
30 files.  These were either inconsistent unit descriptions or descriptions
using a non-standard syntax.  Only comments are changed, and all answers are
bitwise identical.

Fix rotation in set_coupler_type_data

`rotate_array` in `set_coupler_type_data` was trying to rotate an array
to another of different size when `idim` and `jdim` are present.  Some
compilers seemed to let this through, but it raised a
double-deallocation error in GCC.

I'm guessing it's because the rotation was implicitly allocating to a
new array which was automatically deallocated.  But I did not confirm
this.

This was modified to rotate onto a new array of the same size.  The
`idim` and `jdim` are passed to `CT_set_data`, which (hopefully) sorts
out the indexing.

Remove implicit copies in CT_extract_data rotation

Following on the previous commit, the implicit copies in
`extract_coupler_type_data`'s `allocate_rotated_array` and
`rotate_array` are replaced with the full-sized arrays, with halos
managed by `idim` and `jdim` arguments to `CT_extract_data`.

I tested this in a rotated `Baltic` test and saw no answer changes.  The
`ocean.stats` and CFC restart files agree, but there are still known
rotation reordering and negative-zero errors in `MOM.res.nc` output.

Refactor spherical_harmonics_forward

  Refactored the spherical_harmonics_forward routine in MOM_spherical_harmonics
to work in rescaled units by making use of the unscale arguments to
reproducing_sum().  A total of 4 rescaling variables were moved into unscale
arguments, and a block of code that restores the scaling of the output variables
was eliminated.  The comments describing the units of several variables in this
module were made more explicit.  The two temporary arrays that store un-summed
contributions to the transforms were also moved out of the control structure and
made into local variables in the spherical_harmonics_forward routine to allow
for the reuse of that memory.  All answers and diagnostics are bitwise
identical, and no interfaces are changed.

Rescale 7 ice shelf variables

  Changed the rescaling of 7 ice shelf variables to cancel out common
conversion factors that appear in several expressions.

  The C_basal_friction argument to initialize_ice_C_basal_friction is now in
partially rescaled units, reflecting the portion that does not cancel out
fractional-power units from a power law fit with an arbitrary power.  The
C_basal_friction array in ice_shelf_dyn_CS and the C_friction variable in
initialize_ice_C_basal_friction were similarly rescaled.  There are new
scale factors in a get_param_call and a MOM_read_data call and a conversion
factor in register_restart_field call that reflect these changes.

  The KE_tot and mass_tot variables in write_ice_shelf_energy,
are kept in scaled units until they are written.

  The internal variable fN in calc_shelf_taub is kept in scaled units, but there
is now a scaling factor of US%L_to_N in the expression for fB.  This latter
could be folded into the CF_Max element of ice_shelf_dyn_CS, but I am unsure
whether this would be physically sensible.

  All answers should be bitwise identical and no output should change, but this
has not been extensively tested yet.

+Rename visc%TKE_BBL to visc%BBL_meanKE_loss

  Revised the name of visc%TKE_BBL to visc%BBL_meanKE_loss to better reflect
what these variables contain, and the fact that the efficiency of the conversion
from mean kinetic energy loss to turbulent kinetic energy has not been applied
yet.  The units of visc%BBL_meanKE_loss are [H L2 T-3 ~> m3 s-3 or W m-2],
whereas those of visc%TKE_BBL were [H Z2 T-3 ~> m3 s-3 or W m-2].  The factor
rescaling between the units of mean kinetic energy and those of turbulent
kinetic energy have been incorporated into set_diffusivity_CS%BBL_effic and
energetic_PBL_CS%ePBL_BBL_effic.  All answers are bitwise identical.

+Rename fluxes%TKE_tidal to fluxes%BBL_tidal_dis

  Revised the name of fluxes%TKE_tidal to fluxes%BBL_tidal_dis to better reflect
what this field holds,and the fact that the efficiency of the conversion from
mean kinetic energy loss to turbulent kinetic energy has not been applied yet.
The units of fluxes%BBL_tidal_dis are [R Z L2 T-3 ~> W m-2], whereas those of
fluxes%TKE_tidal were [R Z3 T-3 ~> W m-2].  The factor rescaling between the
units of mean kinetic energy and those of turbulent kinetic energy were already
in set_diffusivity_CS%BBL_effic, and these have been cancelled out by this
change, but this is offset by the addition of rescaling factors in the term
setting this array in the NUOPC and mct convert_IOB_to_fluxes routines.  In the
FMS_cap version of convert_IOB_to_fluxes, the extra rescaling factor is rolled
into the scaling factor used in the get_param call for rho_TKE_tidal.  All
answers are bitwise identical, but there is a change in the name and rescaled
units of an element of a transparent type.

+Add EPBL_BBL_TIDAL_EFFIC

  Added the option to use the energy source in fluxes%BBL_tidal_dis to drive
mixing in ePBL_BBL_column(), similarly to what is done in add_drag_diffusivity()
and add_LOTW_BBL_diffusivity().  Because this was omitted when ePBL_BBL_column()
was first added, a separate runtime parameter, EPBL_BBL_TIDAL_EFFIC, is used to
determine the extent to which this tidal energy source is used.  The default for
EPBL_BBL_TIDAL_EFFIC is currently set to 0 to avoid changing answers, but
perhaps it should be changed follow the value of EPBL_BBL_EFFIC.  By default all
answers are bitwise identical, but there is a new runtime parameter in the
MOM_parameter_doc files for cases that use ePBL.

+(*)EPBL_BBL_EFFIC_BUG & DRAG_DIFFUSIVITY_ANSWER_DATE

  Previously, visc%BBL_meanKE_loss (which was called visc%TKE_BBL before it was
renamed) was missing a factor of the square root of the drag coefficient
compared with the net loss of kinetic energy.  This was compensated for by
multiplication by factors of the square root of the drag coefficient in
add_drag_diffusivity and add_LOTW_BBL_diffusivity, but when an equivalent
expression was added in the ePBL BBL code this was erroneously omitted.
Moreover, Alistair has had a comment questioning this added factor in
add_LOTW_BBL_diffusivity for a decade without adequate resolution.

  To correct this confusing situation, visc%BBL_meanKE_loss has been changed to
include the missing factor of the square root of the drag coefficient, while the
new variable visc%BBL_meanKE_loss_sqrtCd was added to allow for the previous
answers to be recovered when the new runtime parameter EPBL_BBL_EFFIC_BUG is set
to true and DRAG_DIFFUSIVITY_ANSWER_DATE is set below 20250302.  Because the
ePBL bottom boundary layer was only added to dev/gfdl a month ago and has not
yet been merged into main, we can be confident that it has only received very
limited use as yet, so the default for EPBL_BBL_EFFIC_BUG is false but this
default will change answers when EPBL_BBL_EFFIC > 0.  The default for
DRAG_DIFFUSIVITY_ANSWER_DATE is 20250101, which will preserve the previous
answers, but the default should later be taken from DEFAULT_ANSWER_DATE.  By
default, answers are unchanged in any cases that are more than a month old, but
answers can change by default in a few very recent experiments.  There are two
new runtime parameters in some MOM_parameter_doc files.

+Add alternate gravity variable GV%g_Earth_Z_T2

  Added the new gravitational acceleration element GV%g_Earth_Z_T2 to the
verticalGrid_type as a copy of GV%g_Earth that uses dimensional rescaling of
[Z T-2 ~> m s-2] instead of [L2 Z-1 T-2 ~> m s-2].   GV%g_Earth_Z_T2 is more
convenient for single-column energy calculations, but the two will only differ
by an integer power of two when dimensional rescaling is applied.  In addition,
the apparently unused element GV%mks_g_Earth was commented out, and it will be
eliminated in several months if there are no concerns raised by MOM6 users whose
code is not on any of the primary development branches of MOM6.

  Explicitly rescaled versions of GV%g_Earth were replaced by GV%g_Earth_Z_T2 in
38 places in 14 files, leading to the elimination of 38 US%L_to_Z**2
dimensional rescaling factors.

  Another US%L_to_Z**2 rescaling factor was eliminated from the code in
ePBL_column by folding it into MKE_to_TKE_effic.

  All answers are bitwise identical, but there are changes to the contents of a
transparent type.

Revise rescaling in bulkmixedlayer

  Revised the internal rescaling of turbulent kinetic energies in
MOM_bulk_mixed_layer to work in units of [H Z2 T-2] instead of [H L2 T-2], for
greater consistency with the rescaling of TKE elsewhere.  These changes included
using GV%g_Earth_Z_T2 instead of GV%g_Earth in 16 places that contribute to
potential energy calculations.  In 5 lines, the unit scaling factors were
eliminated, while in 2 others the rescaling factors were revised.  In addition,
the rescaling factors to go from the units of mean kinetic energy to those of
turbulent kinetic energy were folded into the internal representation of the
bulk Richardson numbers.  The units in comments describing 58 variables
and the conversion arguments for 10 diagnostics were updated accordingly.  All
answers and output are bitwise identical.

Rescale strat_floor

  Rescaled the internal representation of the FGNV_STRAT_FLOOR variable and
thickness_diffuse%N2_floor to include appropriate factors of US%Z_to_L to
reflect the scaling of aspect ratios.  This leads to the elimination of one
rescaling factor in thickness_diffuse_full() and its replacement by another
in a scale factor for a get_param call.  All answers are bitwise identical.

Add Option to Specify Tracer Advection Time Step (#757)

* Add option to set tracer advection timestep

Add an option to seperate the tracer advection from the diabatic
processes in MOM_step_thermo. Previously the advection and
diabatic codes were seperated but called with the same timestep.
If DT_TADVECT is not specified, then DT_THERM is used.
The comments describing DT_THERM and DT_TADVECT have been modified
to refelect this change.

* Add error warning if diabatic_first is true

Currently, this option is not intended to be used if diabatic_first is true,
so a fatal error flag has been added.
This will be addressed in a future pull request.

* Change variable, add do advection before thermo step

Change DT_TADVECT to DT_TRACER_ADVECT to be clearer. dt_tadvect
has been changed to dt_tr_adv and tadvect_spans_coupling to
tradv_spans_coupling.

An additional check has been added before the advection step.
If thermo_spans_coupling is true, do_advection is true if
t_dyn_rel_thermo is ~DT_THERM so that the thermodynamics step would
happen at the next possible time. This ensures that the
thermodynamics step is not prevented because advection has not been
resolved.

* Add logic to do_diabatic check

The do_diabatic check should only be done if do_thermo is true. This
is relevant when do_dyn or do_thermo are being used from outside
step_MOM to order the updates of the thermodynamics and dynamics,
such as when the interspesed coupling is used.

* Initialize do_diabatic

* Change TRADV_SPANS_COUPLING default

Change the default for TRADV_SPANS_COUPLING to be the same as
THERMO_SPANS_COUPLING, which is read first, instead of false.
Initialize do_advection.

* Doxygen fix

---------

Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea58.ncrc.gov>
Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea68.ncrc.gov>
Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea57.ncrc.gov>
+Add optional conversion argument to register_field

  Added the option to rescale variables as they are written out via
MOM_io_file.  These involved adding optional conversion arguments to
register_field_infra and register_field_nc, which are then stored in a new
element in the MOM_field type, and use the conversion factors to unscale
variables before they are written in the ten write_field routines in
MOM_io_file.

  The new optional arguments to register_field are used in MOM_create_file,
taking their values from the vardesc types sent to this routine.  This commit
also alters modify_vardesc to store the value of the conversion optional
argument in the conversion element of the vardesc type.  Also modified
query_vardesc so that the conversion factor is returned via the conversion
optional argument.  These steps had been intended when these optional arguments
were first added, but for some reason they had not actually been used.

  The conversion values stored in a vardesc type are also now used in the
register_diag_field call in ocean_register_diag.  However, it does not appear
that ocean_register_diag is actually used anymore, so it might be a candidate
for deletion.

  All answers are bitwise identical, but there are new optional arguments to
publicly visible routines.

Specify conversion factors in write_energy

  Revised write_energy to use conversion arguments to var_desc to unscale
variables.  Their units are also documented in the same calls, so this is now
analogous to what is done the register_diag_field calls for diagnostics that
are handled by the MOM_diag_mediator.  All calculations in write_energy and
almost all internal variables are now in rescaled units.  All answers and output
are bitwise identical.

Fix 4 conversion arguments to var_desc calls

  Corrected 4 conversion arguments in calls to var_desc for temperatures and
salinities, so that they are consistent with the units of these variables and
the described purpose of the conversion element of the var_desc type.  Until the
conversion arguments to modify_vardesc and query_vardesc, these incorrect values
were inconsequential, but now they need to be fixed before they are
inadvertently used.  All answers are bitwise identical.

Deprecate TIDE_SAL_SCALAR_VALUE (#819)

Using the recently added `old_name` option in `get_params`.
Calculate volo in scaled units

  Calculate the volo diagnostic in scaled units using the tmp_scale argument to
global_area_integral and undo this scaling with a conversion argument on its
register_scalar_field call.  Also corrected the units in the comments describing
the mass_celland masso variables in calculate_diagnostic_fields.  All answers
are bitwise identical.

Rename masked_area in compute_global_grid_integral

  Renamed the internal variable tmpForSumming to masked_area in
compute_global_grid_integrals and corrected the description of its units from
[m2] to [L2 ~> m2].  All answers are bitwise identical.

Fix the syntax or substance in 10 comment units

  Corrected incorrect syntax in the descriptions of 4 positions variables in
write energy, 1 reused position variable in set_grid_metrics_from_mosaic and 1
position variable in bkgnd_mixing_CS.  Also added the unscaled units to the
comments documenting two lines of calculations in wave_speeds, and added the
scaled units to two variables related to the ice-shelf surface mass balance.
Only comments are changed and all answers are bitwise identical.

Switched the placeholder element of file_OBC_CS

  Replaced the unused real tide_flow element of the file_OBC_CS type with the
logical OBC_file_used element to more clearly reflect that this placeholder
element is only here to avoid having a completely empty type.  The actual value
of this element is irrelevant, but some compilers require that all Fortran types
have at least one element.  This change eliminates a meaningless hard-coded real
variable with undefined units and a misleading name and replaces it with a
logical variable named to reflect its actual purpose.  All answers are bitwise
identical.

Merge branch 'main' into dev/gfdl

Correct bug in kappa shear viscosity with vertex shear option. (#824)

* Correct bug in kappa shear viscosity with vertex shear option.

- Viscosities at vertices along the coast were incorrectly zero'd out.  This commit removes that mask so the non-zero shear driven viscosities can be interpolated from in the model.  This bug caused diffusivities to be very large in channels and potentially near coastlines.
- A bugfix flag is added with an option to use the current behavior for legacy purposes.

* Fix missing paranthesis in previous commit (VS_viscosity_bug)

* Update logging of vertex shear viscosity bug parameter
+Add PHILLIPS_ANSWER_DATE runtime parameter

  Add the new runtime parameter PHILLIPS_ANSWER_DATE to enable the option to use
mathematically equivalent expressions in Phillips_initialize_velocity() that
exactly specify the arithmetic to be used, avoid excess divisions and permit
full rescaling of the internal variables and the elimination of rescaling
variables.  This new slightly answer-changing option is enabled by setting
PHILLIPS_ANSWER_DATE >= 20250101.  For now, the default for PHILLIPS_ANSWER_DATE
is set to 20241231 to avoid changing answers without explicitly setting it.

  This commit also introduces code to use G%grid_unit_to_L to detect and handle
various choices for the units of the G%geolat and G%geolon variables.

  By default, all answers are bitwise identical, but there is a new runtime
parameter in some MOM_parameter_doc files.  This commit changes answers at
roundoff when there is an explicit setting of PHILLIPS_ANSWER_DATE >= 20250101,

+Add zero_zeros optional arg to MOM_write_field

  Added the new zero_zeros optional argument to the 10 MOM_write_field routines
in MOM_io and the 6 rescale_comp_data routines in MOM_domain_infra to cause
negative zeros to replaced with ordinary signless zeros before they are written
out to files.  This has no impact at all on answers, but it does help with
comparisons between rotated restart files, in which meaningless differences
between positive and negative zeros were leading to false differences between
files.  All answers are bitwise identical, and all output is equivalent, but
there are new optional arguments to 16 routines covered by 2 publicly visible
interfaces.

+Add turns argument to MOM_read_data

  Added a new turns optional argument to 6 versions of the MOM_read_data
routines to allow for the reading to override the number of turns in the
MOM_domain that is passed into these routines.  Several internal turns variables
in the same routines were renamed qturns to allow for the new optional
arguments. Also check for whether the MOM_domain%domain_in pointer is associated
before it is used, avoiding a segmentation fault that was occurring when a
restart file is read and ROTATE_INDEX is true.  Also added rotate_array calls to
ensure that the halo values are retained while reading data into a rotated
array.  These changes are necessary to allow for the model to be initialized
from a restart files with rotated grids.  Several instances of continuation line
indentation that do not follow the typical 4-space pattern used elsewhere in the
MOM6 code and documented in the MOM6 style guide were also altered to follow the
standard.  All answers that previously worked are bitwise identical, but there
are new optional arguments to publicly visible interfaces.

(*)+Modified MOM_restart to fix rotated restarts

  Modified MOM_restart so that restart files generated by rotated runs match
unrotated restart files, and the model can be properly initialized from a
restart file when the grid is rotated.  Also added runtime options to convert
negative zeros into ordinary zeros before they are written to restart files
(selected with RESTART_UNSIGNED_ZEROS) and to properly do the checksums on the
velocity points on all of the faces (selected with RESTART_SYMMETRIC_CHECKSUMS).

  Also added the new interfaces copy_restart_var and copy_restart vector to use
the names of restart variables and the pointers stored in the restart control
structure to obtain a copy of the variables as the restart variables with the
option to undo the rotation.  These routines are necessary because the reading
of restart files occurs during a phase of the model initialization that works on
an unrotated grid, and they are called from inside of MOM_initialize_state.

  The ranges for the checksums are now set correctly for each variable,
depending on where it is discretized, but when RESTART_SYMMETRIC_CHECKSUMS is
false, the previous ranges are still used so answers do not change in unrotated
test case.

  The conversion factors used for the pair of register_restart_field calls in
register_restart_pair now include the necessary sign changes for the rotation,
as set in the new internal routine set_conversion_pair.  There is also now a
scalar_pair optional argument to the register_restart_pair routines to
accommodate the rotation of pairs of scalars that do not change sign when
rotated (e.g., grid-lengths).

  Instead of working with the hor_grid character string, the restart code has
been modified to instead use the encoded integer position argument returned
from query_vardesc.  This avoids several redundant blocks of code that translate
the hor_grid strings into positions.

  All answers are bitwise identical when there is no grid rotation, but with
grid rotation the restart files that are created are modified to have the
correct signs and replicate the restart fields with no rotation.  Also, cases
with grid rotation can now be reinitialized from restart files, while previously this
simply did not work, either giving an incorrect reinitialized state or a
segmentation fault.  There are also two new runtime parameters in some
MOM_parameter_doc files.

Describe the units of 33 real function results

  Added a description of the units of the return values to the comments
describing 33 real functions in 8 modules.  Only comments are changed and all
answers are bitwise identical.

Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting in PGF (#810)

* Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting

This commit introduces the runtime variable `MASS_WEIGHT_IN_PGF_VANISHED_ONLY`
which has default False. If true, then the `MASS_WEIGHT_IN_PRESSURE_GRADIENT`
and `MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP` effect of weighting T/S integrals
in slanted grid cell FV PGF calculation is turned off if both sides of the
grid cell are nonvanished, where nonvanished means thickness greater than
`RESET_INTXPA_H_NONVANISHED` which defaults to 1e-6 m. Since the benefit of
`MASS_WEIGHT_IN_PRESSURE_GRADIENT` happens in vanished layers (creating a
fake PGF away from vanished layer, which is arrested by upwinded viscosity)
the benefit is still there, but now we can use `MASS_WEIGHT_IN_PRESSURE_GRADIENT`
for coordinates that also have slanted layers in the open ocean that are
not vanished, e.g. sigma coordinates or SIGMA_SHELF_ZSTAR coordinates in the
ice shelf where we DO trust T and S values. Additionally, this is required
near a grounding line in a 3D z-coord ice shelf as some strange looking
slanted non-vanished cells can emerge, and MWIPG being on would create fake PGFs
in non-vanished cells (and therefore generating spurious currents).

Reccommend `MASS_WEIGHT_IN_PGF_VANISHED_ONLY` to be set to True, as well as
`MASS_WEIGHT_IN_PRESSURE_GRADIENT` and `MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP`
if you have vanishing layers with min thickness < 0.1m.

Also modifies MassWt_u and MassWt_v diagnostics to reflect usage
of MASS_WEIGHT_IN_PRESSURE_GRADIENT.

This commit should not change answers since it defaults to False. However,
my implementation is not very efficient and should probably be optimised.

* Minor style updates to previous commit of Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting

* Address comments from Bob by adding scaling to dimensional numbers, replacing Boussinesq rho in nonBoussinesq code, and removing white space to follow 2-space indenting.
+Obsoleted INTERNAL_TIDE_CORNER_ADVECT

  This commit cleans up several aspects of the MOM_internal_tides code, including
the elimination of one runtime option, the correction of some units in comments
and a start toward enabling the propagation directions to alt…
alperaltuntas pushed a commit to NCAR/MOM6 that referenced this pull request Jul 1, 2025
* Leith+E Smoothing Efficiency

This commit replaces the smoothing in Leith+E (which applies to
the velocity components, the coefficient `m_leithy`, and the
biharmonic coefficient `Ah`). The old way applied a 3x3 weighted
moving average twice. The new version applies a 5x5 weighted
moving average once. In exact arithmetic the results are
equivalent, but the new way reduces number of floating point
operations by more than a factor of 3.

* Bug fix for write_energy with short dt (#749)

Fix a bug with subroutine write_energy when using a DT<2. Otherwise,
the energy outputs are written at wrong time steps.

The reason was that time type divide is essentially a floor.
So DT/2 = 0 if DT<2.
Remove extra copy of compute_global_grid_integrals

  The subroutine compute_global_grid_integrals appeared in both the
MOM_state_initialization and MOM_shared_initialization modules, but was only
being called from the latter.   This commit removes the extra copy in
MOM_state_initialization.  It also removes some unnecessary parentheses in the
copy that is being retained, in part to facilitate the review of this commit.
All answers are bitwise identical, and no publicly visible interfaces are
altered.

Make REMAPPING_USE_OM4_SUBCELLS the default

In addition to REMAPPING_USE_OM4_SUBCELLS, for ALE remapping, there are
several parameters of the form XXX_REMAPPING_USE_OM4_SUBCELLS, where
XXX identifies the target, and they all currently default to True.

To simplify setting them all to False, which is recommended, the defaults
for the XXX versions is changed to the value of REMAPPING_USE_OM4_SUBCELLS.

Answers are only changed if REMAPPING_USE_OM4_SUBCELLS is set to False
and the default (now False) is used for one or more of the other parameters.
In such cases the original behaviour can be recovered by explicitly
setting the other parameters to True.

Removed default for mandatory time scale in OBCs

Removed two instances of `fail_if_missing=.true., default=0.` which
are contradictory: a default value is meaningless if the parameter must
be specified.

I encountered this when adding the `defaults=` option to `get_param_real_array()`.

Adds a vector of default values to get_param_real_array()

The `default=` optional argument to get_param() only provides a uniform
value to initialize an array of reals. This commit adds the optional
`defaults=` argument that must have the same length as the `values`
argument.

I've also added a few instances of this optional argument:
 - by adding the `initialize_thickness_param()` procedure, selected by
   `THICKNESS_CONFIG = "param"`. The procedure was based on the "uniform"
   method, and uses the parameter `THICKNESS_INIT_VALUES` which defaults
   to uniform values derived from `MAXIMUM_DEPTH`
 - the setting of MLD_EN_VALS in MOM_diabatic_driver.F90 which was
   previously using a work around to set defaults to 25, 2500, 250000 J/m2.
 - two vectors of 4 values in user/user_change_diffusivity.F90

There will be some doc file changes, but no answer changes.

FMS API: Convert real kind of constants

Two latent heat constants are imported directly from FMS, which is built
independently of MOM6.  Previously, it was a safe assumption that both
would be built with double precision, but this is no longer the case
since FMS now supports both single and double precision.  This could
cause conflicts with mixed-precision builds.

This patch converts the values from FMS-precision to MOM-precision.

Single->double should not affect reproducibility since every
single-precision number can be exactly represented in double precision.
Double->single could affect reproducibility, but this is not an issue
since MOM6 does not run in single precision.

Inline harmonic analysis (#744)

* Inline harmonic analysis

Important bug fix:
    1) The Cholesky decomposition was operating on entries below
the main diagonal of FtF, whereas in the accumulator of FtF, only
entries along and above the main diagonal were calculated. In this
revision, I modified HA_accum_FtF so that entries below the main
diagonal are accumulated instead.
    2) In the accumulator of FtSSH, the first entry for the mean
(zero frequency) is moved out of the loop over different tidal
constituents, so that it is not accumulated multiple times within
a single time step.

* Inline harmonic analysis

Another bug fix: initial state added back to the mean state.

* Inline harmonic analysis

Minor update to HA_solver
Tidal angular frequency has units [rad s-1] (#764)

* Tidal angular frequency has units [rad s-1]

Tidal frequencies are always angular frequencies to simplify applying
sine and cosine.  These have MKS units [rad s-1] but they are all
currently listed as [s-1].

Updated dOxygen comments for variables, e.g. [T-1 ~> s-1] becomes
[rad T-1 ~> rad s-1].  Updated get_param units. e.g. units="s-1"
becomes units="rad s-1".

No answers are changed, but the logged parameter units are different.

There are frequencies in MOM_internal_tides.F90 but these have not
been updated because they may be specified incorrectly.  They are
used as if they are [T-1] but they are calculated as 2PI/period [rad T-1].

  real, allocatable, dimension(:) :: frequency  !< The frequency of each band [T-1 ~> s-1].

  real    :: period             ! A tidal period read from namelist [T ~> s]

  ! The periods of the tidal constituents for internal tides raytracing
  call read_param(param_file, "TIDAL_PERIODS", periods)

  do fr=1,num_freq
    period = US%s_to_T*extract_real(periods, " ,", fr, 0.)
    CS%frequency(fr) = 8.0*atan(1.0)/period
  enddo

All MOM6-examples cases have INTERNAL_TIDES=False and so can't
resolve this issue.

* fixed too-long line
Refactor of vertical reconstruction adding six new schemes (#741)

add diagnostic for Kd_Work related to the Kd_add part (#765)

* add diagnostic for Kd_Work related to the Kd_add part

* also corrects some capitalization typos

* fix id init

---------

Co-authored-by: Raphael Dussin <Raphael.Dussin@noaa.gov>
Co-authored-by: Alistair Adcroft <adcroft@users.noreply.github.com>
+*Obsolete WIND_CONFIG = "SCM_ideal_hurr" (#770)

Changed the code to issue a fatal error message when WIND_CONFIG =
"SCM_ideal_hurr", with the message including instructions on how to recover
mathematically equivalent solutions, and eliminated the subroutine
SCM_idealized_hurricane_wind_forcing() from the Idealized_hurricane module.  The
ocean_only MOM_parameter_doc files have also been modified to reflect that
"SCM_ideal_hurr" is no longer a valid setting for WIND_CONFIG.   All answers are
bitwise identical in any cases that run, but some cases may fail during
initialization with instructions on how to fix them.
*Fix a bug when EPBL_ORIGINAL_PE_CALC is false

  Corrected a bug that causes ePBL column to set the wrong variable and then use
an uninitialized variable when EPBL_ORIGINAL_PE_CALC is set to false.  This bug
was present when the EPBL_ORIGINAL_PE_CALC was first added on Sept. 30, 2016,
but it was not detected because only the default case with EPBL_ORIGINAL_PE_CALC
= True appears to being used or tested.  Any runs that used this code with
debugging compile options would have trapped it immediately.  This will change
answers when EPBL_ORIGINAL_PE_CALC is false.

+Add EPBL_MLD_ITER_BUG runtime parameter

  Added the new runtime parameter EPBL_MLD_ITER_BUG that can be set to false to
correct buggy logic that gives the wrong bounds for the next iteration  when
USE_MLD_ITERATION is true and successive guesses increase by exactly
EPBL_MLD_TOLERANCE.  By default all answers are bitwise identical, but there is
a new runtime parameter in some MOM_parameter_doc files.

+Add run-time ability to debug ePBL sensitivities

  Added the ability to passively run ePBL_column twice in a diagnostic mode and
then provide diagnostics of the differences in the diffusivities and boundary
layer depths that are generated with the two options.  This is controlled by the
new runtime parameter EPBL_OPTIONS_DIFF, which is an integer that specifies
which options to change or 0 (the default) to disable this capability.
Associated with this are the new diagnostics ePBL_opt_diff_Kd_ePBL,
ePBL_opt_maxdiff_Kd_ePBL and ePBL_opt_diff_h_ML, which only are registered when
the differencing is enabled.  For now, only changes associated with the settings
of EPBL_ORIGINAL_PE_CALC and EPBL_ANSWER_DATE can be evaluated, but this list
will grow as new options are added.

  As a part of these changes, there were some other reforms to the way that
MOM_energetic_PBL handles diagnostics, with the 2-d and 3-d arrays changed from
allocatable arrays with enduring memory commitments in the energetic_PBL_CS type
into simple arrays in energetic_PBL, relying on the fact that compilers will be
smart enough to avoid actually allocating this memory when it is unused to avoid
expanding the overall memory requirements of MOM6.  A number of allocate and
deallocate calls were eliminated by these changes.

  In addition, explicit 'units=' descriptors were added to numerous
register_diag_field calls to help identify inconsistent conversion factors.

  Also, the diagnostic of the number of boundary layer depth iterations was
added to the ePBL_column_diags, which enabled the intent of the CS and G
arguments to ePBL_column to be changed to intent(in).  The unnecessary extra
logic associated with the use of the OBL_converged variable in ePBL_column was
eliminated, but the results are uneffected.

  Because the MOM_parameter_doc files were changing anyway, and because this
commit is only a diagnostic change, about 6 spelling errors were corrected in
ePBL parameter descriptions as a part of this commit.

  All answers are bitwise identical, but there is a new runtime parameter and
there may be new diagnostics depending on the choice of a non-default value for
that new parameter.

+Add and test find_Kd_from_PE_chg

  Added the new internal subroutine find_Kd_from_PE_chg inside of the
MOM_energetic_PBL module to directly calculate an increment in the diapycnal
diffusivity from an energy input.  This can be used when ePBL does not convert
released mean kinetic energy into turbulent kinetic energy (i.e., if
MKE_TO_TKE_EFFIC = 0.) and is more efficient than the more general iterative
approach.  To preserve old answers, this new option is only enabled for the
surface boundary layer when the new runtime parameter DIRECT_EPBL_MIXING_CALC is
set to true.  This new option can be tested passively by setting
EPBL_OPTIONS_DIFF to 3 in a run that uses ePBL.  By default all answers are
bitwise identical, but there is a new runtime parameter in some of the
MOM_parameter_doc files.

+(*)Add ePBL bottom boundary mixing option

  Add the option to do energetically consistent bottom boundary layer mixing
with the new routine ePBL_BBL_column.  ePBL_BBL_column is closely based on the
surface-focused ePBL mixing in ePBL_column, but without adding convective
instability driven mixing or mean-TKE driven mixing to avoid possible
double-counting.  This new option is enabled by setting the new runtime
parameter EPBL_BBL_EFFIC to be positive.

  If both EPBL_BBL_EFFIC and BBL_EFFIC are set to positive values, there is a
risk of double-counting, but this case is not being trapped for now.

  The changes include the addition of a new mandatory vertvisc_type argument to
the publicly visible routine energetic_PBL.

  When this new ePBL bottom boundary layer mixing option is enabled, there are
several new diagnostics available that are related to bottom boundary layer
mixing.  Several new checksum calls were also added with this new option when
DEBUG = True. The MOM_parameter_doc files are altered by the addition of two new
runtime parameters, and by the correction of several spelling errors in the
descriptions of other ePBL parameters.  By default, all answers are bitwise
identical.

+Add DECAY_ADJUSTED_BBL_TKE option for ePBL BBL

  Added the ability to modify the bottom boundary layer TKE budget to account
for an exponential decay of TKE away from the boundary and the fact that when
the diffusivity is increased at an interface, it causes an increased buoyancy
flux that varies linearly throughout a well mixed bottom boundary layer and
through the layer above.  This new capability is enabled by the new runtime
parameter DECAY_ADJUSTED_BBL_TKE and is implemented via the new internal
function exp_decay_TKE_adjust.

  In addition, this commit adds 9 bottom-boundary layer specific run-time
parameters that are analogous to parameters that set the properties of the
surface-driven mixing, and take their defaults from them, but can now be set
independently for the bottom boundary layer.  The inefficient option to use
bisection to estimate the bottom boundary layer depth was eliminated, as it
certain that it is not being used yet in any cases, and it should be deprecated
for the estimation of the surface boundary layer.

  By default, all answers are bitwise identical, but there are up to 10 new
runtime parameters that will appear in some MOM_parameter_doc files.

+Add optional unscale argument to reproducing_sum

  Added the new optional unscale argument to reproducing_sum() and
reproducing_sum_EFP().  When this is used, the resulting real sum is restored to
the same scaling as the input arrays, but when there is an EFP_type returned it
is the same as it would have been had the input array been rescaled before it
was passed in.  All answers are bitwise identical, but thre is a new optional
argument to a two publicly visible interfaces.

+Add RZL2_to_kg element to unit_scale_type

  Added the new element RZL2_to_kg to the unit_scale_type for convenience when
rescaling masses and other globally integrated variables.  All answers are
bitwise identical, but there is a new element in a transparent type.

Work in rescaled units in write_energy

  Revised write_energy() and accumulate_net_input()  to work more extensively in
dimensionally rescaled variables by using the new unscale arguments to the
reproducing_sum functions. As a result of these changes, 15 rescaling factors
were eliminated or moved toward the end of write_energy().  All answers are
bitwise identical.

Use reproducing_sum with unscale in 5 places

  Use reproducing_sum with unscale to calculate the global mass diagnostic,
globally integrated ocean surface area, offline tracer transport residuals and
in calculating the OBC inflow area in tidal_bay_set_OBC_data().  This change allows
for the elimination or replacement of 6 rescaling factors and one added instance
with a consistent conversion factor and diagnostic units occur on the same line.
All answers are bitwise identical.

Bug fix related to directional internal wave drag (#774)

Co-authored-by: Robert Hallberg <Robert.Hallberg@noaa.gov>
Fix indexing error in extract_surface_state

Correct a bug that G%gridLonT/G%gridLatT were not using global indexing
in outputing extreme surface information.

Minor change of SSH check in extract_surface_state

Previously extreme surface message is triggered when `sea_lev` is
smaller than or **equal** to ocean topography:

sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) - G%Z_ref

The equality is innocuous for a non-zero minimum thickness
(Angstrom/=0), but it can be problematic for true zero thickness.

Besides, there is the fourth criterion in this check:

sfc_state%sea_lev(i,j) + G%bathyT(i,j) + G%Z_ref < CS%bad_val_col_thick

This is supposed to allow a user-specified tolerance (default=0.0),
which should be a more stringent check when the tolerance is larger than
zero. The equality from the first check contradicts this logic.

Recover diagnostic "SSH_inst"

Fix a bug that sea surface height averaged over every dynamic time step
(SSH_inst) is not outputted correctly.

Update MOM_wave_interface.F90 (#784)

* Update MOM_wave_interface.F90

The index at the interface has been changed from I-1 ->i+1 and J-1 -> j+1, which helps fixed model error in 3D simulations while keeping halo_size to 1.

* Update MOM_wave_interface.F90

Corrected the index for thickness which fixed the issue of 3D wave coupled simulations while keeping the halo_size in thickness as 1
Updates to use EPBL_BBL_EFFIC
- The present code only tests BBL_EFFIC when deciding whether to set the bottom TKE and ustar, but this means they are all zero when EPBL_BBL_EFFIC is non-zero and BBL_EFFIC is zero.
- Adds logic to also check EPBL_BBL_EFFIC, thereby allowing non-zero ustar and TKE for EPBL_BBL_EFFIC>0.0
- Fixed BBL_TKE diagnostic in EPBL that was not populated.
- Will change answers when EPBL_BBL_EFFIC>0.0, but won't change answers in any of our existing configurations.

Move MOM_generic_tracer with no changes (#22)

* Move MOM_generic_tracer with no changes

Renames src/tracer/MOM_generic_tracer.F90
to config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90

* Remove MOM_generic_tracer

Git rm MOM_generic_tracer.F90

---------

Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea68.ncrc.gov>
Corrected the rescaling of 5 KPP diagnostics

  Added missing factors to the conversion arguments in the register_diag_field
calls for the diagnostics of the KPP non-local transport tendency and the net
surface tracer fluxes, and corrected the dimensional scaling that is being
applied to the KPP_Vt2 diagnostic as calculated in KPP_compute_BLD.  All
solutions are bitwise identical, but now output files with these 3 sets of KPP
diagnostics are invariant to dimensional rescaling.  This can be verified with
the visc.nc file generated by the single_column/KPP test case in MOM6-examples.
Whereas previously the diagnostics KPP_Vt2, KPP_QminusSW, KPP_netSalt,
KPP_NLT_dTdt and KPP_NLT_dSdt in that file would change when dimensional
rescaling was applied, now they do not.  No output is changed unless dimensional
rescaling is used.

+Add a vector of defaults to get_param_array_int

  The `default=` optional argument to get_param() only provides a uniform value
to initialize an array of integers.  This commit adds an optional `defaults=`
argument to get_param_int_array(), doc_param_int_array() and
log_param_int_array() to allow for the specification of an array of default
values.  These additions are analogous to what had previously been added for
real arrays in github.com/NOAA-GFDL/MOM6/pull/760.

  This commit also adds the new internal function int_array_string(), analogous
to real_array_string(), in MOM_document.  This differs slightly from its real
array counterpart in that it only uses the syntax like `3*75` for lists of
integers that are longer than we would use to specify dates and times or pairs
of layout parameters, because "(0, 0)" seems more readily interpretable than
"(2*0)".

  The new defaults argument is now used in the get_param calls for LAYOUT and
IO_LAYOUT, and in setting the tidal reference dates.

  Several spelling errors in comments were also corrected in the files that
were being edited.

  All answers are bitwise identical, but there are minor changes in many
MOM_parameter_doc.layout files and some MOM_parameter_doc.all files.

set descale in reproducing_sum_3d
Fix ALE_sponge tendency diagnostic units

  Corrected the units and conversion factor in init_ALE_sponge_diags for the
various sp_tendency_... diagnostics.  Previously they had only been correct for
the tendencies of nondimensional quantities.  The code also now stores the
scaling factor that is set in set_up_ALE_sponge_field_fixed for later use in
registering the sponge tendency diagnostics.  Several instances of unusual
spacing around semicolons in MOM_ALE_sponge were also standardized.  The
documented units and conversion factors for some diagnostics were corrected,
but all solutions are bitwise identical.

Refactor horizontally_average_field

  Refactored the horizontally_average_field() routine in MOM_diag_remap to work
in rescaled units by making use of the unscale arguments to the
reproducing_sum() routines.  A total of 9 rescaling variables were moved into
unscale arguments.  All answers and diagnostics are bitwise identical, and no
interfaces are changed.

+Refactor homogenize_field and revise its interface

  Refactored the homogenize_field routine in MOM_horizontal_regridding to make
use of the unscale argument to reproducing_sum(), and revised its interface to
make it more nearly consistent with the interface to homogenize_field_t() in
MOM_forcing_type.

  The interface changes include revising the order of the arguments, making the weight
argument options, replacing the scale argument with an optional tmp_scale
argument that is the inverse of the previous scale, and making the default for
the use of reproducing sums to be true when the answer_date argument is absent.
The two homogenize_field routines now give equivalent behavior when none of the
optional arguments to homogenize_field() are absent.  The homogenize_field calls
in MOM_temp_salt_initialize_from_Z() and the horiz_interp_and_extrap_tracer()
routines have been modified in accordance with the interface changes.

  All answers are bitwise identical, but the interface to a publicly visible
routine has been substantially changed to the point where any calls using the
previous interface will not compile.

Revise the ice_shelf dimensional rescaling

  Refactored the volume_above_floatation, write_ice_shelf_energy and
ice_shelf_solve_outer routines to work in rescaled units by making use of the
unscale arguments to reproducing_sum().  Also added or corrected comments
documenting the units of 11 real variables in these routines.

  The routine integrate_over_Ice_sheet_area was converted into a function and
var_scale was renamed to unscale for more consistency with the rest of the MOM6
code.  Additionally, add_shelf_flux and update_shelf_mass were modified to use
the scale arguments to time_interp_external.

  A total of 12 rescaling variables were eliminated or moved into unscale
arguments, and 2 blocks of code that scale input variables were eliminated.  All
answers and diagnostics are bitwise identical, and no interfaces are changed.

(*)Conversion arguments for ice shelf scalar diags

  Revised volume_above_floatation and integrate_over_ice_sheet_area to return
values in scaled units, and added conversion factors to the
register_scalar_field calls in MOM_ice_shelf so that all unit unscaling of
diagnostics occurs via  the diag mediator and not in the code itself.  After
this commit, all of the dimensional scaling factors for diagnostics in the
ice_shelf code occur via conversion factors that are immediately adjacent to the
declaration of the units of those diagnostics, facilitating comparison for
consistency.

  The declared units of one diagnostic, "taub_beta", were revised from "MPa s
m-1" to "MPa yr m-1" to be consistent with its conversion factor, which was also
corrected (essentially by a factor of [Z L-1 ~> 1].  Several missing unit
conversion factors for rates of mass change were also added, and about 15
missing or incorrect units in lines documenting real variables were added or
fixed.

  The input scale factor for the (perhaps unused) variable INPUT_VEL_ICE_SHELF
was corrected from `US%m_s_to_L_T*US%m_to_Z` to just `US%m_s_to_L_T`.

  With these changes, all solutions are bitwise identical, apart from regional
cases with a specified inflow when run with dimensional rescaling.  The ice
shelf diagnostics should now be invariant when dimensional rescaling is applied,
and the units of the ice shelf diagnostics are now all consistent with the
required rescaling factors.

Users set solo driver forcing time records per day

  Refactored the solo_driver version of MOM_surface_forcing to modify how the
time levels to read from forcing files, with a series of 11 new runtime
parameters (WIND_FILE_DAYS_PER_RECORD, etc.) giving the user to specify how many
days each time record spans (for positive values) or the number of time records
in the file per day (for negative values).  If these parameters are set to 31,
the month returned by the FMS `get_date()` function sets the time level.  If
they are set to 0 (the default) the previous behavior is obtained, in which the
time level is determined from the number of records in the file.  To always take
the first record in the file, set these parameters to large values (greater than
1000000 will always work, but in practice smaller values do the same thing).  By
default all answers are bitwise identical, but there are up to 10 new runtime
parameters in some MOM_parameter_doc files for cases that have `WIND_CONFIG =
"file"` or `BUOY_CONFIG = "file"`.

  The list of new runtime parameters is `WIND_FILE_DAYS_PER_RECORD`,
`LONGWAVE_FILE_DAYS_PER_RECORD`, `SHORTWAVE_FILE_DAYS_PER_RECORD`,
`EVAPORATION_FILE_DAYS_PER_RECORD`, `LATENTHEAT_FILE_DAYS_PER_RECORD`,
`SENSIBLEHEAT_FILE_DAYS_PER_RECORD`, `RAIN_FILE_DAYS_PER_RECORD`,
`SHORTWAVE_FILE_DAYS_PER_RECORD`, `RUNOFF_FILE_DAYS_PER_RECORD`,
`SSTRESTORE_FILE_DAYS_PER_RECORD` and  `SALINITYRESTORE_FILE_DAYS_PER_RECORD`.

  The problem identified in github.com/NOAA-GFDL/MOM6/issues/631 is addressed by
this commit, and that issue can be closed once this commit is merged onto the
main branch of MOM6.

Merge branch 'main' into dev/gfdl

Fix dimensional rescaling in predict_MEKE

  Refactored predict_MEKE() to include two missing dimensional rescaling factors
and correct a third and to clearly indicate the units and the nature of the
variables where they are being used, in coordination with Andrew Shao.  Some
white spacing around semicolons in the same MOM_MEKE file was also modified to
follow the patterns used elsewhere in the MOM6 code.  All answers are bitwise
identical in cases that do not use dimensional rescaling.

Refactor SAL in MOM_PressureForce_FV

Self-attraction and loading calculation in Boussinesq pressure gradient
force is refactored to be consistent with the algorithm in
non-Boussinesq version.

Add option for SAL to use bottom pressure anomaly

Using SSH to calculate self-attraction and loading (SAL) is only
accurate for barotropic flows. Bottom pressure anomaly should really be
used for general purposes.

* New runtime parameter
A runtime parameter SAL_USE_BPA is added to use bottom pressure anomaly.
The option requires an input file for long term mean reference bottom
pressure. The reference bottom pressure field is stored with SAL_CS.

* Refactor SAL and tides in Boussinesq mode
As the total bottom pressure is needed for bottom pressure anomaly, SAL
calculation in Boussinesq mode needs to be refactored. In addition,
there is a longstanding bug in Boussinesq mode that interface height is
modified by SAL and tides, and the modified interface height is
erroneously used for density and pressure calculation later on.
Therefore the SAL and tides are refactored by moving their calculations
after pressure is calculated. Tide answers before 20230630 is retained.
Answers after 20230630 are changed for Boussinesq mode.

Rename variables with tides in Pressure Force

* Local variable SSH is renamed to pbot_anom.
* Tides related local variable names are changed to be less confusing.
Notably, `e_sal_tide` is renamed to `e_sal_and_tide` (the summation of
sal and tide), not to be confused with `e_tide_sal`, which is renamed to
`e_tidal_sal` (sal from tides).
* Fix e_tidal diagnostics in Boussinesq mode. The term was not assigned
if tide answer date is >20230630.

Alternative method for SAL and tides in Boussinesq

Add an alternative method for SAL and tides in Boussinesq mode. The
current method adjusts interface heights with geopotential height
anomaly for SAL and tides. For non-Boussinesq mode, the current method
is algebraically the same as taking the gradient of SAL and tide
geopotential (body forcing). For Boussinesq mode, there is a baroclinic
component of tidal forcing and SAL. The alternative method is added to
calculate the gradient of tidal forcing and SAL directly at the cost of
additional multiplications. The new method is controlled by runtime
parameter BOUSSINESQ_SAL_TIDES.

Remove rescaling bottom pressure in SAL

* Instead of rescaling bottom pressure to height unit,
calc_Loving_scaling is modified to be conditionally dimensional. When
calculating self-attraction and loading, Love numbers are now
dimensional when bottom pressure anomaly is used as an input. This
change eliminate Love numbers' dependence on mean seawater density.

* A new coefficient called linear_scaling is added to SAL CS for the
same purpose, although to use bottom pressure anomaly for scalar
approximation is not quite justifiable. A WARNING is given when users
try to do that.

* Two separate field, SSH (sea surface height anomaly) and pbot (total
bottom pressure), are used for calculating self attraction and loading
when SAL_USE_BPA = False and SAL_USE_BPA = True, respectively. The
change is to enhance code readability.

Add answer date flag for tide/SAL

A new date for TIDES_ANSWER_DATE is added to recover answers for tides
in Boussinesq mode before 20250201.

Refactor Love_scaling calculation in SAL module

* Precalcualte a local field `coef_rhoE` to avoid in-loop division and
if-blocks. The unit of coef_rhoE depends on use_bpa.
* Fix a few unit description typos in SAL module and two other files.

+Refactor MOM_opacity to replace hard-coded params

  Refactored MOM_opacity to replace hard-coded dimensional parameters for the
Manizza and Morel opacity fits with run-time parameters, and also added the
runtime parameter OPACITY_BAND_WAVELENGTHS to provide the ability to set the
wavelengths of the bands, even though these are not actually used in MOM6.  By
default, these parameters are all set to the previous hard-coded values, using
the recently added defaults argument to get_param_real_array().  The bounds of
the frequency band label arrays with the MANIZZA_05 opacity scheme were also
corrected when PEN_SW_NBANDS > 3, but it would not be typical to use so many
bands for no purpose and these labeling arrays (optics%min_wavelength_band and
optics%max_wavelength_band) do not appear to be used anywhere.  In addition, the
unused publicly visible routines opacity_manizza and opacity_morel were
eliminated or made private.  All answers are bitwise identical, but there are
new entries in some MOM_parameter_doc files.

implement spatially varying decay rates for internal tides leakage (#794)

* add option to specify horizontally varying decay

* extend to vertical modes

* fix comments

* corrected description of decay_rate array

---------

Co-authored-by: Raphael Dussin <Raphael.Dussin@noaa.gov>
Add ice shelf pressure initialisation bug fix (#800)

* Add ice shelf pressure initialisation bug fix

This commit adds a new parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX
and associated bug fix. This bug occurs when TRIM_IC_FOR_P_SURF
pressure initialisation is used for ice shelf, whilst in
Boussinesq mode and a nonlinear equation of state. The subroutine
trim_for_ice calls cut_off_column_top. This routine in Boussinesq
mode calls find_depth_of_pressure_in_cell in
MOM_density_integrals.F90. This subsequently calls the function
frac_dp_at_pos which calls the density equation of state based on
given salinity, temperature and depth, which is incorrectly converted
into pressure by z position (which is negative) multiplied by g and
rho0. The bug results in incorrect densities from the equation of state
and therefore an imperfect initialisation of layer thicknesses and sea
surface height due to the displacement of water by the ice.

The bug is fixed by multiplying the pressure by -1.

This commit introduces parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX with
default value False to preserve answers. If true, the ice shelf
initialisation is accurate to a high degree with a nonlinear
equation of state.

Answers should not change, except for the added parameter in
MOM_parameter_doc.

The same functions are called by a unit test in MOM_state_initialization;
in this commit I set the MOM_state_initialization unit test to
use the existing bug (with a false flag).

* Resolve indenting and white space inconsitencies with MOM6 style
for ice shelf pressure initialisation bug fix FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX
+Add the optional argument old_name to get_param

  Added the new optional argument old_name to the 8 get_param() routines. This
new capability allows for an archaic parameter name to be specified and for
appropriate warnings encouraging the user to migrate to using the new name
while still setting the parameter as intended, or error messages in the case of
inconsistent setting via the archaic name and the correct name.  The logging
inside of the MOM_parameter_doc files only uses the correct parameter name.

  Also added the new optional argument set to the 8 read_param routines, to
indicate whether a parameter has been found and successfully set.  The new set
argument is now being used in read_param() calls in obsolete_int(),
obsolete_real(), obsolete_char() and obsolete_logical().  Obsolete_logical() in
particular was substantially simplified by the use of this new argument, and is
now only about half as long as it was.  The read_param() set argument is also
used in all of the get_param() routines when they are given an old_name
argument.

  The new old_name argument to get_param() is not yet being used in the version
of the MOM6 code that is being checked in, but it has been tested extensively
by adding or modifying get_param calls in a variant of the initialization code,
and it will be used in an updated version of github.com/NOAA-GFDL/MOM6/pull/725
to gracefully handle the deprecation of 4 parameter names.

  All answers are bitwise identical, but there are new optional arguments to
two widely used interfaces.

+Add grid_unit_to_L to the ocean_grid_type

  Add the new element grid_unit_to_L to the ocean_grid_type and the
dyn_horgrid_type, which can be used to convert the units of the geoLat and
geoLon variables to rescaled horizontal distance units ([L ~> m]) when they are
Cartesian coordinates.  When Cartesian coordinates are not in use,
G%grid_unit_to_L is set to 0.

  This new element of the grid type is used to test for inconsistent grids or to
eliminate rescaling variables in set_rotation_beta_plane(),
initialize_velocity_circular(), DOME_initialize_topography(),
DOME_initialize_sponges(), DOME_set_OBC_data(), ISOMIP_initialize_topography(),
idealized_hurricane_wind_forcing(), Kelvin_set_OBC_data(),
Rossby_front_initialize_velocity(), soliton_initialize_thickness(), and
soliton_initialize_velocity().   These are the instances where this new
variable could be used and bitwise identical answers are recovered.  There are
a few other places where they should be used, but where answers would change,
and these will be deferred to a subsequent commit.

  All answers are bitwise identical, but there are new elements in two
transparent and widely used types.

Update particles_run call to allow accurate particle advection using u, v or uh, vh

The drifters package is able to run in 2 different modes: one where the particles are advected using u and v, and one where they are advected using uh and vh. When u and v are used (i.e. the drifters are advected using the resolved velocities only, with no sub-grd component), the particles package needs the layer thickness that is consistent with these velocities, so particles_run needs to be called before any subgrid scale parameterizations are applied. This was not implemented correctly in my last PR, and is corrected here.

The particles package also needs the correct timestep for particle advection. This is added here.

+Remove dyn_horgrid_type%Rad_Earth

  Remove the unused element Rad_Earth from ocean_grid_type and dyn_horgrid_type.
The dimensionally rescaled equivalent element G%Rad_Earth_L is extensively used,
and it will continue to exist.  G%Rad_Earth_L was introduced in November 27,
2021 as a dimensionally rescaled version of G%Rad_Earth, while the unscaled
version was retained because at the time, the Rad_Earth element of the
dyn_horgrid_type is also used in SIS2.  However, SIS2 has not used G%Rad_Earth
since December 23, 2021, so after 3 years we can now safely remove this unused
element.

  Any cases on other branches that might be impacted by this change will not
compile.  All answers are bitwise identical, but there is one less element in
two transparent types.

Nodal modulation

This commit fixes a few (potential) inconsistencies between open
boundary tidal forcing and astronomical tidal forcing.

1. There was an inconsistency in the code that the nodal modulation
can be calculated in OBC tidal forcing but not in the astronomical
tidal forcing. This commit fixes this bug so that nodal modulation
is applied to both forcings.

2. In the previous version of MOM_open_boundary.F90, a different
set of tidal parameters can be set for open boundary tidal forcing,
leading to potential inconsistency with astronomical tidal forcing.
This commit obsoletes those for open boundary tidal forcing.

3. Another important bug fix is that the equilibrium phase is added
to the SAL term, which was missing in the original code.

Correct unit conversion for BS_coeff_h diagnostic

  Added missing conversion arguments for the register_diag_field calls for the
recently added diagnostics BS_coeff_h and BS_coeff_q.  All answers are bitwise
identical, but two diagnostics will have corrected dimensional rescaling when
EY24_EBT_BS is true.

Streaming Filter

The filters and their target frequencies are no longer hard-coded.
Instead, multiple filters with tidal or arbitrary frequencies as
their target frequencies can be turned on. The filter names are
specified in MOM_input and must consist of two letters/numbers. If
the name of a filter is the same as the name of a tidal constituent,
then the corresponding tidal frequency will be used as its target
frequency. Otherwise, the user must specify the target frequency by
"FILTER_${FILTER_NAME}_OMEGA" in MOM_input.

The restarting capability has also been implemented. Because the
filtering is a point-wise operation, all variables are considered
as fields, even if they are velocity components.

Merge branch 'dev/gfdl' into rescale_predict_MEKE
+Refactor the spatial mean calculations

  Refactored the spatial mean calculations in the functions global_area_mean(),
global_area_mean_u(), global_area_mean_v(), global_layer_mean(),
global_volume_mean() and adjust_area_mean_to_zero() to work in rescaled units by
making use of the unscale arguments to the reproducing_sum routines.  Comments
were also added to explain what the code does when both tmp_scale and unscale
arguments are provided, which effectively acts as to allow for changes in the
rescaled units.

  Global_area_integral() and global_mass_integral() were also similarly
refactored, but with the added difference that when a tmp_scale argument is
provided, the areas or masses in the integrals have units of [L2 A ~> m2 a] or
[R Z L2 A ~> kg a] instead of the mixed units of [m2 A ~> m2 a] or [kg A ~> kg
a].  As a result the code surrounding the 4 instances where global_area_integral
or global_mass_integral were being called with tmp_scale arguments (in MOM.F90
and MOM_ice_shelf.F90) also had to be modified in this same commit.  When the
optional var argument is absent from a call to global_mass_integral() so that it
is the mass itself that is returned, the values (if any) of scale, unscale
and tmp_scale are ignored, but the presence of tmp_scale determines whether the
mass is returned in [R Z L2 ~> kg] or [kg], consistent with the behavior that
would have been obtained if var had been an array of nondimensional 1s.

  This commit also includes a rescaling in the units of the areaT_global and
IareaT_global elements of the ocean_grid_type and dyn_horgrid_type to [L2 ~> m2]
and [L-2 ~> m-2], respectively.  Although the dyn_horgrid_type is shared between
MOM6 and SIS2, these elements are not used in SIS2.

  A total of 12 rescaling factors were eliminated or moved into unscale arguments
as a result of these changes.

  All answers are bitwise identical, but there are changes in the rescaled units
of two elements each in two transparent types, and changes to the rescaling
behavior of two publicly visible routines when they are called with tmp_scale
arguments.

Diagnose area integrated fluxes in rescaled units

  Keep 36 area integrated surface mass, heat or salt flux diagnostics in
forcing_diagnostics in rescaled units by replacing an unscale argument with a
tmp_scale argument to global_area_integral.  Also added the corresponding
conversion arguments to the register_scalar_field calls for each of these
diagnostics, so that the documented units and conversion factors can be used to
confirm the correctness of the documented units for these variables.  The
internal variables used for the integrated or averaged heat, mass or salt fluxes
were also separated for clarity about the units of these variables.  All answers
and diagnostics are bitwise identical.

Fix rescaling of internal tide of debugging code

  Corrected the internal tide rescaling arguments so that some of the debugging
variables have units that are consistent with their documented units.  This
involved changing the scale arguments to global_area_integral to tmp_scale
arguments in 8 places so that the units of the output retain the scaling of the
input variable.  The multiplication by the reverse of the scaling factor was
also added to 4 debugging output messages.

  The documented units of internal wave energies in 5 register_restart_field
calls were previously given as "m3 s-2" in non-Boussinesq mode and "J m-2" in
Boussinesq mode when the reverse was actually true.  This has now been
corrected.

  Also replaced the scale argument to 35 chksum calls with the equivalent but
preferred unscale argument, following the pattern elsewhere in the MOM6 code.

  All answers are bitwise identical, and only debugging code was modified.

Corrected many unit descriptions in comments

  Corrected the descriptions of the units of about 113 variables spread across
30 files.  These were either inconsistent unit descriptions or descriptions
using a non-standard syntax.  Only comments are changed, and all answers are
bitwise identical.

Fix rotation in set_coupler_type_data

`rotate_array` in `set_coupler_type_data` was trying to rotate an array
to another of different size when `idim` and `jdim` are present.  Some
compilers seemed to let this through, but it raised a
double-deallocation error in GCC.

I'm guessing it's because the rotation was implicitly allocating to a
new array which was automatically deallocated.  But I did not confirm
this.

This was modified to rotate onto a new array of the same size.  The
`idim` and `jdim` are passed to `CT_set_data`, which (hopefully) sorts
out the indexing.

Remove implicit copies in CT_extract_data rotation

Following on the previous commit, the implicit copies in
`extract_coupler_type_data`'s `allocate_rotated_array` and
`rotate_array` are replaced with the full-sized arrays, with halos
managed by `idim` and `jdim` arguments to `CT_extract_data`.

I tested this in a rotated `Baltic` test and saw no answer changes.  The
`ocean.stats` and CFC restart files agree, but there are still known
rotation reordering and negative-zero errors in `MOM.res.nc` output.

Refactor spherical_harmonics_forward

  Refactored the spherical_harmonics_forward routine in MOM_spherical_harmonics
to work in rescaled units by making use of the unscale arguments to
reproducing_sum().  A total of 4 rescaling variables were moved into unscale
arguments, and a block of code that restores the scaling of the output variables
was eliminated.  The comments describing the units of several variables in this
module were made more explicit.  The two temporary arrays that store un-summed
contributions to the transforms were also moved out of the control structure and
made into local variables in the spherical_harmonics_forward routine to allow
for the reuse of that memory.  All answers and diagnostics are bitwise
identical, and no interfaces are changed.

Rescale 7 ice shelf variables

  Changed the rescaling of 7 ice shelf variables to cancel out common
conversion factors that appear in several expressions.

  The C_basal_friction argument to initialize_ice_C_basal_friction is now in
partially rescaled units, reflecting the portion that does not cancel out
fractional-power units from a power law fit with an arbitrary power.  The
C_basal_friction array in ice_shelf_dyn_CS and the C_friction variable in
initialize_ice_C_basal_friction were similarly rescaled.  There are new
scale factors in a get_param_call and a MOM_read_data call and a conversion
factor in register_restart_field call that reflect these changes.

  The KE_tot and mass_tot variables in write_ice_shelf_energy,
are kept in scaled units until they are written.

  The internal variable fN in calc_shelf_taub is kept in scaled units, but there
is now a scaling factor of US%L_to_N in the expression for fB.  This latter
could be folded into the CF_Max element of ice_shelf_dyn_CS, but I am unsure
whether this would be physically sensible.

  All answers should be bitwise identical and no output should change, but this
has not been extensively tested yet.

+Rename visc%TKE_BBL to visc%BBL_meanKE_loss

  Revised the name of visc%TKE_BBL to visc%BBL_meanKE_loss to better reflect
what these variables contain, and the fact that the efficiency of the conversion
from mean kinetic energy loss to turbulent kinetic energy has not been applied
yet.  The units of visc%BBL_meanKE_loss are [H L2 T-3 ~> m3 s-3 or W m-2],
whereas those of visc%TKE_BBL were [H Z2 T-3 ~> m3 s-3 or W m-2].  The factor
rescaling between the units of mean kinetic energy and those of turbulent
kinetic energy have been incorporated into set_diffusivity_CS%BBL_effic and
energetic_PBL_CS%ePBL_BBL_effic.  All answers are bitwise identical.

+Rename fluxes%TKE_tidal to fluxes%BBL_tidal_dis

  Revised the name of fluxes%TKE_tidal to fluxes%BBL_tidal_dis to better reflect
what this field holds,and the fact that the efficiency of the conversion from
mean kinetic energy loss to turbulent kinetic energy has not been applied yet.
The units of fluxes%BBL_tidal_dis are [R Z L2 T-3 ~> W m-2], whereas those of
fluxes%TKE_tidal were [R Z3 T-3 ~> W m-2].  The factor rescaling between the
units of mean kinetic energy and those of turbulent kinetic energy were already
in set_diffusivity_CS%BBL_effic, and these have been cancelled out by this
change, but this is offset by the addition of rescaling factors in the term
setting this array in the NUOPC and mct convert_IOB_to_fluxes routines.  In the
FMS_cap version of convert_IOB_to_fluxes, the extra rescaling factor is rolled
into the scaling factor used in the get_param call for rho_TKE_tidal.  All
answers are bitwise identical, but there is a change in the name and rescaled
units of an element of a transparent type.

+Add EPBL_BBL_TIDAL_EFFIC

  Added the option to use the energy source in fluxes%BBL_tidal_dis to drive
mixing in ePBL_BBL_column(), similarly to what is done in add_drag_diffusivity()
and add_LOTW_BBL_diffusivity().  Because this was omitted when ePBL_BBL_column()
was first added, a separate runtime parameter, EPBL_BBL_TIDAL_EFFIC, is used to
determine the extent to which this tidal energy source is used.  The default for
EPBL_BBL_TIDAL_EFFIC is currently set to 0 to avoid changing answers, but
perhaps it should be changed follow the value of EPBL_BBL_EFFIC.  By default all
answers are bitwise identical, but there is a new runtime parameter in the
MOM_parameter_doc files for cases that use ePBL.

+(*)EPBL_BBL_EFFIC_BUG & DRAG_DIFFUSIVITY_ANSWER_DATE

  Previously, visc%BBL_meanKE_loss (which was called visc%TKE_BBL before it was
renamed) was missing a factor of the square root of the drag coefficient
compared with the net loss of kinetic energy.  This was compensated for by
multiplication by factors of the square root of the drag coefficient in
add_drag_diffusivity and add_LOTW_BBL_diffusivity, but when an equivalent
expression was added in the ePBL BBL code this was erroneously omitted.
Moreover, Alistair has had a comment questioning this added factor in
add_LOTW_BBL_diffusivity for a decade without adequate resolution.

  To correct this confusing situation, visc%BBL_meanKE_loss has been changed to
include the missing factor of the square root of the drag coefficient, while the
new variable visc%BBL_meanKE_loss_sqrtCd was added to allow for the previous
answers to be recovered when the new runtime parameter EPBL_BBL_EFFIC_BUG is set
to true and DRAG_DIFFUSIVITY_ANSWER_DATE is set below 20250302.  Because the
ePBL bottom boundary layer was only added to dev/gfdl a month ago and has not
yet been merged into main, we can be confident that it has only received very
limited use as yet, so the default for EPBL_BBL_EFFIC_BUG is false but this
default will change answers when EPBL_BBL_EFFIC > 0.  The default for
DRAG_DIFFUSIVITY_ANSWER_DATE is 20250101, which will preserve the previous
answers, but the default should later be taken from DEFAULT_ANSWER_DATE.  By
default, answers are unchanged in any cases that are more than a month old, but
answers can change by default in a few very recent experiments.  There are two
new runtime parameters in some MOM_parameter_doc files.

+Add alternate gravity variable GV%g_Earth_Z_T2

  Added the new gravitational acceleration element GV%g_Earth_Z_T2 to the
verticalGrid_type as a copy of GV%g_Earth that uses dimensional rescaling of
[Z T-2 ~> m s-2] instead of [L2 Z-1 T-2 ~> m s-2].   GV%g_Earth_Z_T2 is more
convenient for single-column energy calculations, but the two will only differ
by an integer power of two when dimensional rescaling is applied.  In addition,
the apparently unused element GV%mks_g_Earth was commented out, and it will be
eliminated in several months if there are no concerns raised by MOM6 users whose
code is not on any of the primary development branches of MOM6.

  Explicitly rescaled versions of GV%g_Earth were replaced by GV%g_Earth_Z_T2 in
38 places in 14 files, leading to the elimination of 38 US%L_to_Z**2
dimensional rescaling factors.

  Another US%L_to_Z**2 rescaling factor was eliminated from the code in
ePBL_column by folding it into MKE_to_TKE_effic.

  All answers are bitwise identical, but there are changes to the contents of a
transparent type.

Revise rescaling in bulkmixedlayer

  Revised the internal rescaling of turbulent kinetic energies in
MOM_bulk_mixed_layer to work in units of [H Z2 T-2] instead of [H L2 T-2], for
greater consistency with the rescaling of TKE elsewhere.  These changes included
using GV%g_Earth_Z_T2 instead of GV%g_Earth in 16 places that contribute to
potential energy calculations.  In 5 lines, the unit scaling factors were
eliminated, while in 2 others the rescaling factors were revised.  In addition,
the rescaling factors to go from the units of mean kinetic energy to those of
turbulent kinetic energy were folded into the internal representation of the
bulk Richardson numbers.  The units in comments describing 58 variables
and the conversion arguments for 10 diagnostics were updated accordingly.  All
answers and output are bitwise identical.

Rescale strat_floor

  Rescaled the internal representation of the FGNV_STRAT_FLOOR variable and
thickness_diffuse%N2_floor to include appropriate factors of US%Z_to_L to
reflect the scaling of aspect ratios.  This leads to the elimination of one
rescaling factor in thickness_diffuse_full() and its replacement by another
in a scale factor for a get_param call.  All answers are bitwise identical.

Add Option to Specify Tracer Advection Time Step (#757)

* Add option to set tracer advection timestep

Add an option to seperate the tracer advection from the diabatic
processes in MOM_step_thermo. Previously the advection and
diabatic codes were seperated but called with the same timestep.
If DT_TADVECT is not specified, then DT_THERM is used.
The comments describing DT_THERM and DT_TADVECT have been modified
to refelect this change.

* Add error warning if diabatic_first is true

Currently, this option is not intended to be used if diabatic_first is true,
so a fatal error flag has been added.
This will be addressed in a future pull request.

* Change variable, add do advection before thermo step

Change DT_TADVECT to DT_TRACER_ADVECT to be clearer. dt_tadvect
has been changed to dt_tr_adv and tadvect_spans_coupling to
tradv_spans_coupling.

An additional check has been added before the advection step.
If thermo_spans_coupling is true, do_advection is true if
t_dyn_rel_thermo is ~DT_THERM so that the thermodynamics step would
happen at the next possible time. This ensures that the
thermodynamics step is not prevented because advection has not been
resolved.

* Add logic to do_diabatic check

The do_diabatic check should only be done if do_thermo is true. This
is relevant when do_dyn or do_thermo are being used from outside
step_MOM to order the updates of the thermodynamics and dynamics,
such as when the interspesed coupling is used.

* Initialize do_diabatic

* Change TRADV_SPANS_COUPLING default

Change the default for TRADV_SPANS_COUPLING to be the same as
THERMO_SPANS_COUPLING, which is read first, instead of false.
Initialize do_advection.

* Doxygen fix

---------

Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea58.ncrc.gov>
Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea68.ncrc.gov>
Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea57.ncrc.gov>
+Add optional conversion argument to register_field

  Added the option to rescale variables as they are written out via
MOM_io_file.  These involved adding optional conversion arguments to
register_field_infra and register_field_nc, which are then stored in a new
element in the MOM_field type, and use the conversion factors to unscale
variables before they are written in the ten write_field routines in
MOM_io_file.

  The new optional arguments to register_field are used in MOM_create_file,
taking their values from the vardesc types sent to this routine.  This commit
also alters modify_vardesc to store the value of the conversion optional
argument in the conversion element of the vardesc type.  Also modified
query_vardesc so that the conversion factor is returned via the conversion
optional argument.  These steps had been intended when these optional arguments
were first added, but for some reason they had not actually been used.

  The conversion values stored in a vardesc type are also now used in the
register_diag_field call in ocean_register_diag.  However, it does not appear
that ocean_register_diag is actually used anymore, so it might be a candidate
for deletion.

  All answers are bitwise identical, but there are new optional arguments to
publicly visible routines.

Specify conversion factors in write_energy

  Revised write_energy to use conversion arguments to var_desc to unscale
variables.  Their units are also documented in the same calls, so this is now
analogous to what is done the register_diag_field calls for diagnostics that
are handled by the MOM_diag_mediator.  All calculations in write_energy and
almost all internal variables are now in rescaled units.  All answers and output
are bitwise identical.

Fix 4 conversion arguments to var_desc calls

  Corrected 4 conversion arguments in calls to var_desc for temperatures and
salinities, so that they are consistent with the units of these variables and
the described purpose of the conversion element of the var_desc type.  Until the
conversion arguments to modify_vardesc and query_vardesc, these incorrect values
were inconsequential, but now they need to be fixed before they are
inadvertently used.  All answers are bitwise identical.

Deprecate TIDE_SAL_SCALAR_VALUE (#819)

Using the recently added `old_name` option in `get_params`.
Calculate volo in scaled units

  Calculate the volo diagnostic in scaled units using the tmp_scale argument to
global_area_integral and undo this scaling with a conversion argument on its
register_scalar_field call.  Also corrected the units in the comments describing
the mass_celland masso variables in calculate_diagnostic_fields.  All answers
are bitwise identical.

Rename masked_area in compute_global_grid_integral

  Renamed the internal variable tmpForSumming to masked_area in
compute_global_grid_integrals and corrected the description of its units from
[m2] to [L2 ~> m2].  All answers are bitwise identical.

Fix the syntax or substance in 10 comment units

  Corrected incorrect syntax in the descriptions of 4 positions variables in
write energy, 1 reused position variable in set_grid_metrics_from_mosaic and 1
position variable in bkgnd_mixing_CS.  Also added the unscaled units to the
comments documenting two lines of calculations in wave_speeds, and added the
scaled units to two variables related to the ice-shelf surface mass balance.
Only comments are changed and all answers are bitwise identical.

Switched the placeholder element of file_OBC_CS

  Replaced the unused real tide_flow element of the file_OBC_CS type with the
logical OBC_file_used element to more clearly reflect that this placeholder
element is only here to avoid having a completely empty type.  The actual value
of this element is irrelevant, but some compilers require that all Fortran types
have at least one element.  This change eliminates a meaningless hard-coded real
variable with undefined units and a misleading name and replaces it with a
logical variable named to reflect its actual purpose.  All answers are bitwise
identical.

Merge branch 'main' into dev/gfdl

Correct bug in kappa shear viscosity with vertex shear option. (#824)

* Correct bug in kappa shear viscosity with vertex shear option.

- Viscosities at vertices along the coast were incorrectly zero'd out.  This commit removes that mask so the non-zero shear driven viscosities can be interpolated from in the model.  This bug caused diffusivities to be very large in channels and potentially near coastlines.
- A bugfix flag is added with an option to use the current behavior for legacy purposes.

* Fix missing paranthesis in previous commit (VS_viscosity_bug)

* Update logging of vertex shear viscosity bug parameter
+Add PHILLIPS_ANSWER_DATE runtime parameter

  Add the new runtime parameter PHILLIPS_ANSWER_DATE to enable the option to use
mathematically equivalent expressions in Phillips_initialize_velocity() that
exactly specify the arithmetic to be used, avoid excess divisions and permit
full rescaling of the internal variables and the elimination of rescaling
variables.  This new slightly answer-changing option is enabled by setting
PHILLIPS_ANSWER_DATE >= 20250101.  For now, the default for PHILLIPS_ANSWER_DATE
is set to 20241231 to avoid changing answers without explicitly setting it.

  This commit also introduces code to use G%grid_unit_to_L to detect and handle
various choices for the units of the G%geolat and G%geolon variables.

  By default, all answers are bitwise identical, but there is a new runtime
parameter in some MOM_parameter_doc files.  This commit changes answers at
roundoff when there is an explicit setting of PHILLIPS_ANSWER_DATE >= 20250101,

+Add zero_zeros optional arg to MOM_write_field

  Added the new zero_zeros optional argument to the 10 MOM_write_field routines
in MOM_io and the 6 rescale_comp_data routines in MOM_domain_infra to cause
negative zeros to replaced with ordinary signless zeros before they are written
out to files.  This has no impact at all on answers, but it does help with
comparisons between rotated restart files, in which meaningless differences
between positive and negative zeros were leading to false differences between
files.  All answers are bitwise identical, and all output is equivalent, but
there are new optional arguments to 16 routines covered by 2 publicly visible
interfaces.

+Add turns argument to MOM_read_data

  Added a new turns optional argument to 6 versions of the MOM_read_data
routines to allow for the reading to override the number of turns in the
MOM_domain that is passed into these routines.  Several internal turns variables
in the same routines were renamed qturns to allow for the new optional
arguments. Also check for whether the MOM_domain%domain_in pointer is associated
before it is used, avoiding a segmentation fault that was occurring when a
restart file is read and ROTATE_INDEX is true.  Also added rotate_array calls to
ensure that the halo values are retained while reading data into a rotated
array.  These changes are necessary to allow for the model to be initialized
from a restart files with rotated grids.  Several instances of continuation line
indentation that do not follow the typical 4-space pattern used elsewhere in the
MOM6 code and documented in the MOM6 style guide were also altered to follow the
standard.  All answers that previously worked are bitwise identical, but there
are new optional arguments to publicly visible interfaces.

(*)+Modified MOM_restart to fix rotated restarts

  Modified MOM_restart so that restart files generated by rotated runs match
unrotated restart files, and the model can be properly initialized from a
restart file when the grid is rotated.  Also added runtime options to convert
negative zeros into ordinary zeros before they are written to restart files
(selected with RESTART_UNSIGNED_ZEROS) and to properly do the checksums on the
velocity points on all of the faces (selected with RESTART_SYMMETRIC_CHECKSUMS).

  Also added the new interfaces copy_restart_var and copy_restart vector to use
the names of restart variables and the pointers stored in the restart control
structure to obtain a copy of the variables as the restart variables with the
option to undo the rotation.  These routines are necessary because the reading
of restart files occurs during a phase of the model initialization that works on
an unrotated grid, and they are called from inside of MOM_initialize_state.

  The ranges for the checksums are now set correctly for each variable,
depending on where it is discretized, but when RESTART_SYMMETRIC_CHECKSUMS is
false, the previous ranges are still used so answers do not change in unrotated
test case.

  The conversion factors used for the pair of register_restart_field calls in
register_restart_pair now include the necessary sign changes for the rotation,
as set in the new internal routine set_conversion_pair.  There is also now a
scalar_pair optional argument to the register_restart_pair routines to
accommodate the rotation of pairs of scalars that do not change sign when
rotated (e.g., grid-lengths).

  Instead of working with the hor_grid character string, the restart code has
been modified to instead use the encoded integer position argument returned
from query_vardesc.  This avoids several redundant blocks of code that translate
the hor_grid strings into positions.

  All answers are bitwise identical when there is no grid rotation, but with
grid rotation the restart files that are created are modified to have the
correct signs and replicate the restart fields with no rotation.  Also, cases
with grid rotation can now be reinitialized from restart files, while previously this
simply did not work, either giving an incorrect reinitialized state or a
segmentation fault.  There are also two new runtime parameters in some
MOM_parameter_doc files.

Describe the units of 33 real function results

  Added a description of the units of the return values to the comments
describing 33 real functions in 8 modules.  Only comments are changed and all
answers are bitwise identical.

Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting in PGF (#810)

* Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting

This commit introduces the runtime variable `MASS_WEIGHT_IN_PGF_VANISHED_ONLY`
which has default False. If true, then the `MASS_WEIGHT_IN_PRESSURE_GRADIENT`
and `MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP` effect of weighting T/S integrals
in slanted grid cell FV PGF calculation is turned off if both sides of the
grid cell are nonvanished, where nonvanished means thickness greater than
`RESET_INTXPA_H_NONVANISHED` which defaults to 1e-6 m. Since the benefit of
`MASS_WEIGHT_IN_PRESSURE_GRADIENT` happens in vanished layers (creating a
fake PGF away from vanished layer, which is arrested by upwinded viscosity)
the benefit is still there, but now we can use `MASS_WEIGHT_IN_PRESSURE_GRADIENT`
for coordinates that also have slanted layers in the open ocean that are
not vanished, e.g. sigma coordinates or SIGMA_SHELF_ZSTAR coordinates in the
ice shelf where we DO trust T and S values. Additionally, this is required
near a grounding line in a 3D z-coord ice shelf as some strange looking
slanted non-vanished cells can emerge, and MWIPG being on would create fake PGFs
in non-vanished cells (and therefore generating spurious currents).

Reccommend `MASS_WEIGHT_IN_PGF_VANISHED_ONLY` to be set to True, as well as
`MASS_WEIGHT_IN_PRESSURE_GRADIENT` and `MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP`
if you have vanishing layers with min thickness < 0.1m.

Also modifies MassWt_u and MassWt_v diagnostics to reflect usage
of MASS_WEIGHT_IN_PRESSURE_GRADIENT.

This commit should not change answers since it defaults to False. However,
my implementation is not very efficient and should probably be optimised.

* Minor style updates to previo…
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants