+Cleanup MOM_set_diffusivity and MOM_geothermal #1228
Conversation
Merged the calculations in sfc_bkgnd_mixing into calculate_bkgnd_mixing, changed some of the arguments to calculate_bkgnd_mixing, and moved the diagnostics in MOM_bkgnd_mixing to MOM_set_diffusivity to reduce the use of elements of a control structure across modules. Also made the Kd_lay arguments to some of the internal subroutines in MOM_set_diffusivity optional. All answers are bitwise identical.
Made the Kd_lay arguments to calculate_tidal_mixing, calculate_CVMix_tidal and add_init_tide_diffusivity optional. All answers are bitwise identical.
Added () to argumentless subroutine calls and corrected some indents to the 2-space MOM6 standard. All answers are bitwise identical.
Made a series of changes to MOM_bkgnd_mixing to eliminate the unnecessary use of 3-d arrays for diagnostics. These include passing in 2-d slices of memory for use with the layer and interface diffusivities and interface viscosity to calculate_bkgnd_mixing, and the elimination of the kd_bkgnd and kv_bkgnd elements of bkgnd_mixing_cs. The diagnostics of background diffusivities and viscosities are still being made available, but they are only stored in a 3-d array if the diagnostics are in use. With these changes the contents of the control structure for MOM_bkgnd_mixing can now be made private. All answers are bitwise identical, but there are changes to public interfaces.
Moved the registration calls for the double diffusivity diagnostics from MOM_CVMix_ddiff to MOM_set_diffusivity, so there are no longer duplicated diagnostic registration calls, and there are not 3-d arrays being allocated in the CVmix double diffusion module. Also, the contents of the control structure for MOM_CVMix_ddiff are now private. There is a new optional argument to compute_ddiff_coeffs. All solutions are bitwise identical, but the order of some entries in the available_diags files may have changed in some cases.
Added extra non-logged get_parameter calls to set_diffusivity_init and the new publicly visible subroutine tidal_mixing_h_amp to the MOM_tidal_mixing module, so that the elements of tidal_mixing_cs can be declared to be private. All answers are bitwise identical.
Made the argument Kd_lay to set_diffusivity into an optional argument, and pass a 2-d (i-,k-) slice of this diffusivity to many of the subroutines called from within set_diffusivity, including calculate_tidal_mixing. All answers are bitwise identical, but there are changes to public interfaces.
Pass around and work with a 2-d (i-,k-) slice of Kd_int within set_diffusivity, including calculate_tidal_mixing. All answers are bitwise identical, but there are changes to public interfaces.
Eliminated the variable Kd_lay from diabatic_ALE and diabatic_ALE_legacy, taking advantage of the fact that Kd_lay is now an optional argument to set_diffusivity. Also eliminated unused logical branches in diabatic_ALE_legacy that had not been cleaned up when this was separated off from diabatic_ALE. All answers are bitwise identical.
Moved the init and end calls for tidal_mixing into the set_diffusivity module, simplifying connections between modules and the set_diffusivity_init interface. All answers are bitwise identical, but there are changes in the order of entries in the MOM_parameter_doc files.
Reuse the arrays with the temperature and salinity after vertical smoothing for the temperature and salinities after convective adjustment. Both are temporary arrays and do not reappear as diagnostics. This required some minor reordering of the calls setting these two adjusted temperatures. All answers are bitwise identical.
Overloaded the interface to geothermal to either call the existing version or a much simpler version when used in ALE mode with no movement of mass between layers. Also modified the diabatic_ALE routines to use this new interface. All answers are bitwise identical.
Eliminated several variables or variable copies that are no longer needed due to the simplified interface to geothermal. Also added the optional argument useALEalgorithm to geothermal_init, which is used to enable a diagnostic that only makes sense when not using the ALE mode of geothermal heating. All answers are bitwise identical.
Refactored the diagnostics in geothermal_in_place to use fewer 3-d arrays and simplified the code. All solutions are bitwise identical, but there may be changes at the level of roundoff in the internal_heat_heat_tendency and internal_heat_temp_tendency diagnostics.
Reordering of code in calculate_bkgnd_mixing for greater clarity and a less convoluted structure. Also added comments highlighting and explaining a bug that has been in the MOM6 code for about 18 months. All answers are bitwise identical.
Added the runtime parameter PRANDTL_EPBL to specify the Prandtl number used convert the diffusivities from ePBL into viscosities. With the default value (1.0), all answers are bitwise identical to the previous version of the code, but there is a new runtime parameter in some MOM_parameter_doc files.
Corrected the calculation of the diagnostics internal_heat_temp_tendency and internal_heat_h_tendency in the original version of geothermal. They had been wrong because they are based on differences between the temperature and thickness at the end of the subroutine and stored values, but the stored values were only being saved in a part of the computational domain and set to 0 elsewhere. This does not impact the model solutions, but it causes the regression testing of diagnostics to fail in several of the TCs. All solutions are bitwise identical.
|
This is a complicated PR. It will pass the pipeline regression tests apart form changes in the MOM_parameter_doc files, but as expected the Travis CI tests are failing in the regression checks of the diagnostics for three separate reasons:
If this PR is too complicated to evaluate, the various commits could be broken up into multiple stages:
Please note that several of these changes would have to be carried out in sequence if the testing and approval is to be simplified, and that some of these PRs can not even be requested until the predecessor has been approved. Note also that many of these PRs would occur in the opposite of the order in which they were developed, so replacing this one PR as a series of 6 PRs would require all of the code changes to be redone. However, if the prospective reviewers are finding this to complicated an onerous to review, and if one of the other people with write access to dev/gfdl were willing to commit to handling PRs 1-5 in this list promptly, I would be willing to to withdraw this large PR and redo it in the 6 PR sequence described above. If we were to decide to redo this one PR as a sequence of PRs, there is a potentially answer changing bug-fix in MOM_bkgnd_mixing restoring the answers from prior to March, 2018 (it does not happen to change anything in MOM6-examples) that this large PR highlights with comments but does not correct, which we might want to put into this sequence in parallel with items 1, 2, 4 or 5 in the broken-out list above. Otherwise this bug fix will follow in a small PR of its own. |
|
Concerning the |
| @@ -411,7 +338,11 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kv, j, G, GV, US, CS) | |||
|
|
|||
| ! set some parameters | |||
| deg_to_rad = atan(1.0)/45.0 ! = PI/180 | |||
| epsilon = 1.e-10 | |||
| min_sinlat = 1.e-10 | |||
There was a problem hiding this comment.
Would change answers, but should this be a limitation on the minimum latitude instead (slightly easier to interpret physically)?
There was a problem hiding this comment.
Curiously, if you are willing to work in radians instead of degrees, this value of the minimum value of sin(lat) is also the minimum value of lat to within 64-bit machine precision, since sin(x) = x * (1 - x^2/3! + x^4/5! - ...) and (1.e-10)**2/6 < 1.0e-17. So maybe we just need to revise the description in the comment (or not).
| @@ -154,14 +161,17 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo) | |||
| .or. CS%id_internal_heat_temp_tendency > 0 | |||
|
|
|||
| if (CS%id_internal_heat_heat_tendency > 0) work_3d(:,:,:) = 0.0 | |||
There was a problem hiding this comment.
- I don't think this resetting of the
work3darray is necessary - Potentially outside the scope of the PR, but refactor the variable names of the diagnostic ids to
id_internal_heatingso we don't get things likeinternal_heat_heat
There was a problem hiding this comment.
It is necessary to set work3d to 0 here if this heat_tendency diagnostic is in use, for precisely the same reason why the code storing h_old and T_old needed to be moved outside of the loops with conditional exits where the heating occurs.
I strongly agree that shorter diagnostic id names here would be less cumbersome.
There was a problem hiding this comment.
This is overall a good tidy up.
The new geothermal_in_place() is a much cleaner implementation for when we are in general coordinate mode. Stripping out the interleaving of entrainment calculations is both more efficient and makes the code easier to understand. It is a shame the reverse can't be done for isopycnal_geothermal().
The one change I request: I do not think this is useful example of overloaded interfaces so please call geothermal_in_place() and isopycnal_geothermal() directly. When following the call tree it is better to know where you are going. I've spent too much of my life reading the wrong version of a routine.
| real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. | ||
| type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. | ||
| ! Kd_T and Kd_S are intent inout because only one j-row is set here, but they are essentially outputs. | ||
| real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kd_T !< Interface double diffusion diapycnal |
There was a problem hiding this comment.
This was a good catch - lucky for use the compilers have been treating it as inout.
There was a problem hiding this comment.
I suspect that we might have detected this if the compilers had not be treating this array as intent(inout).
There was a problem hiding this comment.
The overloaded interface to geothermal() has been replaced with separate interfaces to geothermal_entraining() and geothermal_in_place().
| if (CS%ePBL_is_additive) then | ||
| Kd_add_here = Kd_ePBL(i,j,K) | ||
| visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) | ||
| visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + CS%ePBL_Prandtl*Kd_ePBL(i,j,K) |
There was a problem hiding this comment.
This is a bug fix that doesn't change answers for the M6E suite but does affect OM4.1 for which Pr=1.25 . We can forego a run-time parameter so long as no one else has a non-unity Pr.
There was a problem hiding this comment.
I introduced a new runtime parameter (PRANDTL_EPBL) that is used set ePBL_Prandtl with a default value of 1 so this can't possibly affect OM4.1. We can consider using a global Prandtl number to set the defaults throughout the code, but this PR does not do this.
There was a problem hiding this comment.
Looks good to me and definitely clears up the code flow quite a bit. I'm not sure why the diagnostic tests are failing either.
I agree with @adcroft that the geothermal routines probably don't need to be overloaded. We're already using ALE as a big-switch elsewhere and makes it abundantly clear which processes are applied differently between layer-isopycnal and ALE configurations.
| call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, G, GV, US, CS) | ||
| ! diagnose temperature, salinity, heat, and salt tendencies | ||
| ! Note: hold here refers to the thicknesses from before the dual-entraintment when using | ||
| ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will (not?) have changed |
There was a problem hiding this comment.
Since eatr(k+1)=ebtr(k) in ALE mode, this diagnostic should always be 0. Maybe don't register in ALE mode?
Rearranged the argument declarations for compute_ddiff_coeffs and calculate_bkgnd_mixing to be as close to the order of the arguments themselves as Fortran syntax permits, as suggested in a PR review. Also changed the macros for the vertical extent to use SZK_(GV) instead of SZK_(G). All answers are bitwise identical.
Changed the overloaded interface to geothermal() into geothermal_entraining() and geothermal_in_place(). All answers are bitwise identical but there are changes to a public interface.
Refactored MOM_set_diffusivity.F90, MOM_geothermal.F90, MOM_diabatic_driver.F90
and related modules to improve the code structure, and to avoid unnecessarily setting
or using 3-d arrays that are not used themselves. This includes the elimination
of at least 3 3-d arrays and a dramatic simplification of the geothermal heating
routine in ALE mode. The registration calls for some modules and diagnostics
were changed, so there are changes in the order of the entries in many
MOM_parameter_doc files. Also added the new runtime parameter PRANDTL_EPBL, and
corrected two diagnostics of the thickness and temperature changes due to
geothermal heating that were just wrong before. All tested solutions are bitwise
identical, but the corrections to diagnostics cause some automated regression
tests fail. The commits in this PR include: