-
Notifications
You must be signed in to change notification settings - Fork 72
Adds a vector of default values to get_param_real_array() #760
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Adds a vector of default values to get_param_real_array() #760
Conversation
Codecov ReportAttention: Patch coverage is
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. |
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 |
I think regression uses everthing in its own (Edit: We build 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. |
There was a problem hiding this 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.
FWIW, dropping the new
I'll ask a the dev call if we prefer this version... |
9eec351
to
8f7dd3e
Compare
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()`.
9f15344
to
f8dda20
Compare
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.
I tried assumed-rank arrays for |
There was a problem hiding this 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.
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.
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.
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.
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…
* 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…
The
default=
optional argument to get_param() only provides a uniform value to initialize an array of reals. This commit adds the optionaldefaults=
argument that must have the same length as thevalues
argument.I've also added a few instances of this optional argument:
initialize_thickness_param()
procedure, selected byTHICKNESS_CONFIG = "param"
. The procedure was based on the "uniform" method, and uses the parameterTHICKNESS_INIT_VALUES
which defaults to uniform values derived fromMAXIMUM_DEPTH
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.