Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Adds parameterizations of SGS variance to PGF and isopycnal slopes #1156

Merged
merged 37 commits into from
Jul 11, 2020
Merged

Adds parameterizations of SGS variance to PGF and isopycnal slopes #1156

merged 37 commits into from
Jul 11, 2020

Conversation

adcroft
Copy link
Collaborator

@adcroft adcroft commented Jul 8, 2020

This adds deterministic parameterizations of both i) the Brankart terms in the pressure gradient force (PGF) and ii) related effect on isopycnal slope in the GM parameterization. There is more to come but these changes add up such that maintaining compatibility with dev/gfdl is becoming a burden. The deterministic schemes in this PR are being evaluated in OM4.1. There are quite a few steps including:

  • major re-factor of the line integrals and PGF code
    • variants of the line integrals are now grouped into MOM_density_integrals where there were scattered across PGF and EOS modules before. This happens to localize the logic of where stuff is switched but the primary reason is for being able to access data type that should not be seen by the EOS.
    • the routines always and still, worked on a single layer at a time but now take 3D inputs along with "k". This allows reconstruction to be moved into this routine and avoid 3D pre-calculation (still in progress with a new monotonized PLM scheme). For now the main advantage is to avoid the addition of numerous arguments associated with each of the SGS fields.
  • adds the PPM form of Boussinesq PGF using quadrature. This was once used in 2d for published results but never extended to three dimensions. PPM has been validated by zeroing out the curvature term to bitwise reproduce the PLM results. New tests that use PPM-PGF will be coming later in MOM6-examples.
  • a less major re-factor of some of the ALE reconstruction code because they are used in PGF (includes renaming ppoly_E -> edge_values as discussed a while ago)
  • introduction of new equation of state functions that use SGS variance and covariance as inputs
  • added new fields to the tv data type. This are currently determined by parameterizations but will later be optionally prognostic fields.
  • updated TC2 to enable (test) both these new schemes

Interactions:

  • This adds new parameters to MOM6-examples will need updating once it points to this code.
  • There are no answer changes for existing configurations in MOM6-examples but tc2 and tc2.a were updated to test these schemes so a regression test will fail.

adcroft and others added 30 commits July 7, 2020 12:05
- Stanley et al., 2020, adds the Brankart modification to volume mean
  density via linear corrections involving SGS sample variances and
  covariances of T and S. This commit adds the new interfaces that
  allow a call to calculate_density to include the variances and
  covariance as arguments. This correction sits above the particular
  EOS and thus can use any EOS so long as second derivatives are
  available within the EOS module.
- The routines pressure_gradient_plm() and pressure_gradient_ppm()
  were poorly named and had comments referring to the pressure gradient
  calculation even though the only calculate edge values in the vertical
  for T/S using ALE functions. The routines are actually more general
  and can be used outside of the PGF. The comments have been shortened
  and no longer refer to the PGF.
- The integrals of density routines used (mostly) by the PGF calculation
  were part of MOM_EOS. Originally, when writing the analytic FV PGF paper,
  this was the right place to put the integrals. The additional variants
  using the ALE reconstruction functions mean that it is cleaner to have
  these routines sit in a layer above EOS and ALE.
- ppoly_E meant something to someone a while ago but we felt it
  would be better to clean up the ALE APIs. This is a pre-cursor
  to switching to a more precise description of reconstructions.
- Undid nesting of if / select statemnent
- Changed FATAL messages
- Function had been moved to density_integrals already but not
  deleted
- Also ordered public statements in MOM_EOS to help find things.
- Order was frustratingly illogical
- When converting from the two hor_index types to one (HII and HIO
  became HI) I retained the HI%isd:... declaration statements. We
  normally use `SZI_(G)` or `SZI_(HI)` so this switches to that style.
- This adds elemental functions for cellwise operations within the PLM
  construction procedure which allows the operations to be accessed
  from outside of the ALE functions on different array shapes but recover
  bitwise identical results.
- The older subroutines now use the functions and some optimizations were
  obtained in the process.
- Out of curiosity I want to be sure these weren't implicated in the long
  initialization cost.
- Soft index convention was not properly implemented in original PLM
  density integrals. Fixed to avoid confusion.
- Probably makes no difference but we re-computed weights repeatedly
  even though we'd pre-computed them (and used only once).
  - Reduces number of local variables.
- Removes (now) unused functions for quadrature and polynomial evaluation
- Has been tested by setting the logical use_PPM=.false.
  - reproduces PLM mode bitwise
- Use of new PLM functions led to out of date directives
- Now passing in 3D structures and tv so that we can access other fields
  in tv%.
- Corrected calculation of curvatures s6, t6 for interpolated interior
  line integrals.
- Tested with use_PPM=.false. and PLM edge values.
  - Using PLM edge values with use_PPM=.true. does give different answers
    because t6,s6 are non-zero albeit tiny. This is due to FP truncation
    errors.
- Since the majority of code is not about the analytic integration
  I thought it time to rename this module
- Continuation of rename of module
- This routine was suffering from if's inside loops. I've removed this
  one since it was just obfuscation, pretending to care about performance
  and actually reducing it.
- Makes use of the Stanley equation of state to include effects of
  SGS temperature variance, salinity variance and T-S covariance.
- This now calculates the gradient of SGS temperature variance
  using the same discretization as used for the gradient of
  density along coordinate surfaces.
- Added run-time coefficient for Stanley parameterization
- Fixed openmp directives.
- Alters halo over which vert_fill_TS() is called.
- Add Stanley parameter to tc2 to test new code.
- Renamed STANLEY_DET_COEFF to STANLEY_PRM_DET_COEFF to indicate
  this is for the Stanley parameterization as opposed to a related
  approach by Stanley for implementing the Brankart PGF correction.
- PGF_STANLEY_T2_DET_COEFF>=0. now adds the temperature variance
  contribution to the Brankart correction using the Stanley linearization
  of the EOS and parameterization of SGS variance.
- After starting with allocatables and switching to pointers
  I'd forgotten to nullify them.
- Passing the scale= to calculate_density_second_derivs() was the
  double scaling the density contribution from SGS variance in the
  Stanley parameterization.
- select case statements were indented as if once inside a loop.
- The pressure scaling was wrong when trying to call
  calculate_density_second_derivs_array() from within
  calculate_stanley_density_array() because the latter should not do
  any scaling but the former always did. I had to call the lower level
  functions provided by WRIGHT, TEOS10, etc to avoid get the scaling
  tests to pass.
- In preparing to add Brankartc terms to PLM form of presure gradient
  I'm adapting the routine to look like the PPM routine which takes 3D
  arguments but only works on layer "k".
- The Brankart PGF terms are now implemented in the PLM recontruction
  routines, just as they were for the PPM form.
- Added documentation lines for PGF_STANLEY_T2_DET_COEFF
PLM_slope_wa = PLM_slope_wa * ( 1. - epsilon(PLM_slope_wa) )
endif

! An attempt to avoid inconsistency when the values become unrepresentable.
Copy link
Collaborator

Choose a reason for hiding this comment

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

Please modify this comment and the one at line 107 to note that these are dimensionally inconsistent expressions. It would be better to rewrite the code to handle underflow consistently, (e.g., via a du_underflow argument), but this can be added later. For now leaving a note should suffice.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Address in commit f0bab12

@@ -193,6 +202,43 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, scale)

end subroutine calculate_density_scalar

!> Calls the appropriate subroutine to calculate density of sea water for scalar inputs
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think that perhaps calculate_stanley_density_scalar should follow the same pattern with respect to the dimensions of its arguments and rescaling via the contents of EOS as does calculate_density_scalar, and not the pattern of calculate_density_array, as it currently does. This is especially true if the calculate_density_array routines are eventually going to be hidden from public exposure.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The problem actually goes deeper because second derivatives are always scaled so it is hard to build these Stanley version out of the existing functions. The existing calculate_density_1d() calls calculate_density_array() and to do so has to "unscale" rho_ref. This did not work for the second derivatives so I've had to write out longer versions of these functions. Addressed in commit e8adc48.

@Hallberg-NOAA
Copy link
Collaborator

I have visually examined all of the code change in this PR, and I am supportive of it. There are two minor points where I had previously commented on specific parts of the the commit. Specifically there are two places in the ALE/PLM_functions code where we need comments noting dimensionally inconsistent expressions to deal with underflow, and I have a query about the handling of the units in the corresponding calculate_stanley_density_scalar and calculate_density_scalar routines.

I have also manually run this through the MOM6-examples test suite, and I have verified that this PR does not change answers there. The deliberate fail on the tc2 test case has been noted.

Once these minor points have been resolved, I will move ahead with manually merging this pull request, along with appropriate changes to the MOM_parameter_doc files.

@@ -28,10 +28,10 @@ module MOM_PressureForce
type, public :: PressureForce_CS ; private
logical :: Analytic_FV_PGF !< If true, use the analytic finite volume form
!! (Adcroft et al., Ocean Mod. 2008) of the PGF.
logical :: blocked_AFV !< If true, used the blocked version of the ANALYTIC_FV_PGF
logical :: blocked_FV !< If true, used the blocked version of the ANALYTIC_FV_PGF
Copy link
Collaborator

Choose a reason for hiding this comment

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

blocked_FV is left over from an obsolete option and it not used anywhere in this module, so rather than being renamed, I think that it should simply be deleted.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Addressed in ec0946c

@@ -186,7 +159,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
"MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")

Copy link
Collaborator

Choose a reason for hiding this comment

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

Because the Stanley option is not yet implemented for a non-Boussinesq code, it might be advisable to add an error message here if CS%Stanley_T2_det_coeff>=0. If there are plans to add this in the very near-future this might not be necessary, but it would serve both as a reminder of a missing option and to avoid users thinking that they are using something that they are not.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Addressed in 4307fa5

real, intent(in) :: u_c !< Value of center cell [units of u]
real, intent(in) :: u_r !< Value of right cell [units of u]
! Local variables
real :: sigma_l, sigma_c, sigma_r, u_min, u_max
Copy link
Collaborator

Choose a reason for hiding this comment

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

We should get into the habit of documenting the units of every real variable, and eventually enforce this for all code commits if possible. In this case it would be easy as all have [units of u].

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Addressed in 553166a

real, intent(in) :: u_c !< Value of center cell [units of u]
real, intent(in) :: u_r !< Value of right cell [units of u]
! Local variables
real :: sigma_l, sigma_c, sigma_r, u_min, u_max, h_cn
Copy link
Collaborator

Choose a reason for hiding this comment

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

These local variables should have comments indicating their units, and variables with different units (i.e. h_cn) should be on a separate line if the variable units are being documented collectively.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Addressed in 553166a

real, intent(in) :: s_c !< PLM slope of center cell [units of u]
real, intent(in) :: s_r !< PLM slope of right cell [units of u]
! Local variables
real :: e_r, e_l, edge, almost_two, slp
Copy link
Collaborator

Choose a reason for hiding this comment

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

Please document the units of these variables.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Addressed in 553166a

@codecov-commenter
Copy link

Codecov Report

Merging #1156 into dev/gfdl will decrease coverage by 0.23%.
The diff coverage is 32.27%.

Impacted file tree graph

@@             Coverage Diff              @@
##           dev/gfdl    #1156      +/-   ##
============================================
- Coverage     46.08%   45.85%   -0.24%     
============================================
  Files           214      224      +10     
  Lines         69399    70082     +683     
============================================
+ Hits          31984    32137     +153     
- Misses        37415    37945     +530     
Impacted Files Coverage Δ
...g_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 0.00% <0.00%> (ø)
...fig_src/external/GFDL_ocean_BGC/generic_tracer.F90 0.00% <0.00%> (ø)
...c/external/GFDL_ocean_BGC/generic_tracer_utils.F90 0.00% <0.00%> (ø)
config_src/external/ODA_hooks/kdtree.f90 0.00% <0.00%> (ø)
config_src/external/ODA_hooks/ocean_da_core.F90 0.00% <0.00%> (ø)
config_src/external/ODA_hooks/ocean_da_types.F90 0.00% <0.00%> (ø)
config_src/external/ODA_hooks/write_ocean_obs.F90 0.00% <0.00%> (ø)
config_src/solo_driver/MESO_surface_forcing.F90 0.00% <0.00%> (ø)
config_src/solo_driver/MOM_driver.F90 68.72% <ø> (ø)
config_src/solo_driver/user_surface_forcing.F90 0.00% <0.00%> (ø)
... and 139 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 5aae165...553166a. Read the comment docs.

Copy link
Collaborator

@Hallberg-NOAA Hallberg-NOAA left a comment

Choose a reason for hiding this comment

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

This is a valuable new contribution to MOM6, of which I enthusiastically approve.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants