From 125e3a3d721b922840f4409c6d46c8a647188d7f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Dec 2022 08:54:13 -0500 Subject: [PATCH 1/6] +(*)Fix numerous issues with MOM_stoch_eos Corrected several bugs that would prevent cases with STOCH_EOS = True or that set a non-negative value of STANLEY_COEF from running at all. There was also extensive refactoring of the MOM_stoch_eos.F90 code. Specifically this commit makes the following changes: - The new routine stoch_eos_register_restarts was added to register the restart field associated with the MOM_stoch_EOS module before the restarts are read and before restart_registry_lock is called, which was causing any test case with STOCH_EOS=True or a positive value of STANLEY_COEF to have a fatal error. - Added a missing dimensional rescaling factor in the get_param call for KD_SMOOTH in MOM_stoch_eos_init, which would have caused cases exercising the MOM_stoch_eos options to fail the dimensional consistency tests. - MOM_stoch_eos_init was changed from a subroutine into a function that returns a logical value indicating whether this routine is used further. The order of the arguments was modified to match the order used in every other init call. - The new routine post_stoch_eos_diags was added to write diagnostics associated with the stoch_eos module. - The MOM_stoch_eos_CS type was made opaque. - The unused diag argument was removed from MOM_stoch_eos_run. - Unit arguments were added to the get_param calls for STANLEY_COEFF and STANLEY_A. - Four arrays in the MOM_stoch_eos_CS type that had been declared to optionally use static memory allocation were modified to be simple allocatables, as they are not used in the vast majority of MOM6 cases, and there is no reason to always assign memory to them. - The register_restart_field call for "stoch_eos_pattern" was revised to use the newer, more direct form rather than working via a vardesc type. - Return statements were added to several of the MOM_stoch_eos routines in the cases where they are not supposed to do anything. - The comments describing several real variables in the MOM_stoch_eos module were added or modified to describe their units using the standard format. - The module use statements at the start of the MOM_stoch_eos module were updated to reflect these changes. There were also parallel changes in MOM.F90: - Added use_stochastic_EOS element to the MOM_control_struct to indicate whether the stoch_eos calls are to be used. - MOM_stoch_eos_init is only called if temperature and salinities are state variables, as it makes no sense to call it otherwise. - Calls to stoch_eos routines were updates to reflect the new interfaces. - The contents of the MOM_stoch_eos_CS type are no longer used in MOM.F90. All answers are bitwise identical in cases that do not use dimensional rescaling, but answers will change (be corrected) in some cases that do use dimensional consistency tests. Several public interfaces to MOM_stoch_eos routines were altered, and there are changes to multiple MOM_parameter_doc files due to the new units, and due to the fact that stoch_EOS parameters are no longer being logged in cases where they are meaningless. --- src/core/MOM.F90 | 26 +++-- src/core/MOM_stoch_eos.F90 | 191 ++++++++++++++++++++++--------------- 2 files changed, 132 insertions(+), 85 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index d1fd8619da..8d7caf83ed 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -114,7 +114,8 @@ module MOM use MOM_shared_initialization, only : write_ocean_geometry_file use MOM_sponge, only : init_sponge_diags, sponge_CS use MOM_state_initialization, only : MOM_initialize_state -use MOM_stoch_eos, only : MOM_stoch_eos_init,MOM_stoch_eos_run,MOM_stoch_eos_CS,mom_calc_varT +use MOM_stoch_eos, only : MOM_stoch_eos_init, MOM_stoch_eos_run, MOM_stoch_eos_CS +use MOM_stoch_eos, only : stoch_EOS_register_restarts, post_stoch_EOS_diags, mom_calc_varT use MOM_sum_output, only : write_energy, accumulate_net_input use MOM_sum_output, only : MOM_sum_output_init, MOM_sum_output_end use MOM_sum_output, only : sum_output_CS @@ -288,6 +289,7 @@ module MOM logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics. logical :: mixedlayer_restrat !< If true, use submesoscale mixed layer restratifying scheme. logical :: useMEKE !< If true, call the MEKE parameterization. + logical :: use_stochastic_EOS !< If true, use the stochastic EOS parameterizations. logical :: useWaves !< If true, update Stokes drift logical :: use_p_surf_in_EOS !< If true, always include the surface pressure contributions !! in equation of state calculations. @@ -1079,12 +1081,12 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call cpu_clock_begin(id_clock_dynamics) call cpu_clock_begin(id_clock_stoch) - if (CS%stoch_eos_CS%use_stoch_eos) call MOM_stoch_eos_run(G,u,v,dt,Time_local,CS%stoch_eos_CS,CS%diag) + if (CS%use_stochastic_EOS) call MOM_stoch_eos_run(G, u, v, dt, Time_local, CS%stoch_eos_CS) call cpu_clock_end(id_clock_stoch) call cpu_clock_begin(id_clock_varT) - if (CS%stoch_eos_CS%stanley_coeff >= 0.0) then - call MOM_calc_varT(G,GV,h,CS%tv,CS%stoch_eos_CS,dt) - call pass_var(CS%tv%varT, G%Domain,clock=id_clock_pass,halo=1) + if (CS%use_stochastic_EOS) then + call MOM_calc_varT(G, GV, h, CS%tv, CS%stoch_eos_CS, dt) + if (associated(CS%tv%varT)) call pass_var(CS%tv%varT, G%Domain, clock=id_clock_pass, halo=1) endif call cpu_clock_end(id_clock_varT) @@ -1297,9 +1299,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (IDs%id_u > 0) call post_data(IDs%id_u, u, CS%diag) if (IDs%id_v > 0) call post_data(IDs%id_v, v, CS%diag) if (IDs%id_h > 0) call post_data(IDs%id_h, h, CS%diag) - if (CS%stoch_eos_CS%id_stoch_eos > 0) call post_data(CS%stoch_eos_CS%id_stoch_eos, CS%stoch_eos_CS%pattern, CS%diag) - if (CS%stoch_eos_CS%id_stoch_phi > 0) call post_data(CS%stoch_eos_CS%id_stoch_phi, CS%stoch_eos_CS%phi, CS%diag) - if (CS%stoch_eos_CS%id_tvar_sgs > 0) call post_data(CS%stoch_eos_CS%id_tvar_sgs, CS%tv%varT, CS%diag) + if (CS%use_stochastic_EOS) call post_stoch_EOS_diags(CS%stoch_eos_CS, CS%tv, CS%diag) call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other) @@ -2679,6 +2679,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call waves_register_restarts(waves_CSp, HI, GV, param_file, restart_CSp) endif + if (use_temperature) then + call stoch_EOS_register_restarts(HI, param_file, CS%stoch_eos_CS, restart_CSp) + endif + call callTree_waypoint("restart registration complete (initialize_MOM)") call restart_registry_lock(restart_CSp) @@ -2964,7 +2968,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call interface_filter_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%interface_filter_CSp) new_sim = is_new_run(restart_CSp) - call MOM_stoch_eos_init(G,Time,param_file,CS%stoch_eos_CS,restart_CSp,diag) + if (use_temperature) then + CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, US, param_file, diag, CS%stoch_eos_CS, restart_CSp) + else + CS%use_stochastic_EOS = .false. + endif if (CS%use_porbar) & call porous_barriers_init(Time, US, param_file, diag, CS%por_bar_CS) diff --git a/src/core/MOM_stoch_eos.F90 b/src/core/MOM_stoch_eos.F90 index 2f67077f1e..deb878e99c 100644 --- a/src/core/MOM_stoch_eos.F90 +++ b/src/core/MOM_stoch_eos.F90 @@ -2,46 +2,44 @@ module MOM_stoch_eos ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_grid, only : ocean_grid_type -use MOM_hor_index, only : hor_index_type -use MOM_file_parser, only : get_param, param_file_type -use MOM_random, only : PRNG,random_2d_constructor,random_2d_norm -use MOM_time_manager, only : time_type -use MOM_io, only : vardesc, var_desc -use MOM_restart, only : MOM_restart_CS,is_new_run -use MOM_diag_mediator, only : register_diag_field,post_data,diag_ctrl,safe_alloc_ptr -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type -use MOM_restart, only : register_restart_field -use MOM_isopycnal_slopes,only : vert_fill_TS -!use random_numbers_mod, only : getRandomNumbers,initializeRandomNumberStream,randomNumberStream +use MOM_diag_mediator, only : register_diag_field, post_data, diag_ctrl +use MOM_error_handler, only : MOM_error, FATAL +use MOM_file_parser, only : get_param, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_hor_index, only : hor_index_type +use MOM_isopycnal_slopes, only : vert_fill_TS +use MOM_random, only : PRNG, random_2d_constructor, random_2d_norm +use MOM_restart, only : MOM_restart_CS, register_restart_field, is_new_run, query_initialized +use MOM_time_manager, only : time_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type +!use random_numbers_mod, only : getRandomNumbers, initializeRandomNumberStream, randomNumberStream implicit none; private #include public MOM_stoch_eos_init public MOM_stoch_eos_run +public stoch_EOS_register_restarts +public post_stoch_EOS_diags public MOM_calc_varT !> Describes parameters of the stochastic component of the EOS !! correction, described in Stanley et al. JAMES 2020. -type, public :: MOM_stoch_eos_CS - real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: l2_inv - !< One over sum of the T cell side side lengths squared - real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: rgauss - !< nondimensional random Gaussian - real :: tfac=0.27 !< Nondimensional decorrelation time factor, ~1/3.7 - real :: amplitude=0.624499 !< Nondimensional std dev of Gaussian +type, public :: MOM_stoch_eos_CS ; private + real, allocatable :: l2_inv(:,:) !< One over sum of the T cell side side lengths squared [L-2 ~> m-2] + real, allocatable :: rgauss(:,:) !< nondimensional random Gaussian [nondim] + real :: tfac=0.27 !< Nondimensional decorrelation time factor, ~1/3.7 [nondim] + real :: amplitude=0.624499 !< Nondimensional standard deviation of Gaussian [nondim] integer :: seed !< PRNG seed type(PRNG) :: rn_CS !< PRNG control structure - real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: pattern - !< Random pattern for stochastic EOS [nondim] - real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: phi - !< temporal correlation stochastic EOS [nondim] + real, allocatable :: pattern(:,:) !< Random pattern for stochastic EOS [nondim] + real, allocatable :: phi(:,:) !< temporal correlation stochastic EOS [nondim] logical :: use_stoch_eos!< If true, use the stochastic equation of state (Stanley et al. 2020) real :: stanley_coeff !< Coefficient correlating the temperature gradient - !! and SGS T variance; if <0, turn off scheme in all codes - real :: stanley_a !< a in exp(aX) in stochastic coefficient + !! and SGS T variance [nondim]; if <0, turn off scheme in all codes + real :: stanley_a !< a in exp(aX) in stochastic coefficient [nondim] real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1] !>@{ Diagnostic IDs @@ -52,61 +50,64 @@ module MOM_stoch_eos contains -!> Initializes MOM_stoch_eos module. -subroutine MOM_stoch_eos_init(G, Time, param_file, CS, restart_CS, diag) - type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(time_type), intent(in) :: Time !< Time for stochastic process - type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure - type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. - type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics +!> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used. +logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_CS) + type(time_type), intent(in) :: Time !< Time for stochastic process + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse + type(diag_ctrl), target, intent(inout) :: diag !< Structure used to control diagnostics + type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. ! local variables integer :: i,j - type(vardesc) :: vd - CS%seed=0 - ! contants - !pi=2*acos(0.0) + + MOM_stoch_eos_init = .false. + + CS%seed = 0 + call get_param(param_file, "MOM_stoch_eos", "STOCH_EOS", CS%use_stoch_eos, & "If true, stochastic perturbations are applied "//& "to the EOS in the PGF.", default=.false.) call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", CS%stanley_coeff, & "Coefficient correlating the temperature gradient "//& - "and SGS T variance.", default=-1.0) + "and SGS T variance.", units="nondim", default=-1.0) call get_param(param_file, "MOM_stoch_eos", "STANLEY_A", CS%stanley_a, & "Coefficient a which scales chi in stochastic perturbation of the "//& - "SGS T variance.", default=1.0) + "SGS T variance.", units="nondim", default=1.0, & + do_not_log=((CS%stanley_coeff<0.0) .or. .not.CS%use_stoch_eos)) call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6) + units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s, & + do_not_log=(CS%stanley_coeff<0.0)) - !don't run anything if STANLEY_COEFF < 0 + ! Don't run anything if STANLEY_COEFF < 0 if (CS%stanley_coeff >= 0.0) then + if (.not.allocated(CS%pattern)) call MOM_error(FATAL, & + "MOM_stoch_eos_CS%pattern is not allocated when it should be, suggesting that "//& + "stoch_EOS_register_restarts() has not been called before MOM_stoch_eos_init().") - ALLOC_(CS%pattern(G%isd:G%ied,G%jsd:G%jed)) ; CS%pattern(:,:) = 0.0 - vd = var_desc("stoch_eos_pattern","nondim","Random pattern for stoch EOS",'h','1') - call register_restart_field(CS%pattern, vd, .false., restart_CS) - ALLOC_(CS%phi(G%isd:G%ied,G%jsd:G%jed)) ; CS%phi(:,:) = 0.0 - ALLOC_(CS%l2_inv(G%isd:G%ied,G%jsd:G%jed)) - ALLOC_(CS%rgauss(G%isd:G%ied,G%jsd:G%jed)) + allocate(CS%phi(G%isd:G%ied,G%jsd:G%jed), source=0.0) + allocate(CS%l2_inv(G%isd:G%ied,G%jsd:G%jed), source=0.0) + allocate(CS%rgauss(G%isd:G%ied,G%jsd:G%jed), source=0.0) call get_param(param_file, "MOM_stoch_eos", "SEED_STOCH_EOS", CS%seed, & "Specfied seed for random number sequence ", default=0) call random_2d_constructor(CS%rn_CS, G%HI, Time, CS%seed) call random_2d_norm(CS%rn_CS, G%HI, CS%rgauss) - ! fill array with approximation of grid area needed for decorrelation - ! time-scale calculation + ! fill array with approximation of grid area needed for decorrelation time-scale calculation do j=G%jsc,G%jec do i=G%isc,G%iec - CS%l2_inv(i,j)=1.0/(G%dxT(i,j)**2+G%dyT(i,j)**2) + CS%l2_inv(i,j) = 1.0/(G%dxT(i,j)**2+G%dyT(i,j)**2) enddo enddo - if (is_new_run(restart_CS)) then - do j=G%jsc,G%jec - do i=G%isc,G%iec - CS%pattern(i,j)=CS%amplitude*CS%rgauss(i,j) - enddo - enddo + + if (.not.query_initialized(CS%pattern, "stoch_eos_pattern", restart_CS) .or. & + is_new_run(restart_CS)) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + CS%pattern(i,j) = CS%amplitude*CS%rgauss(i,j) + enddo ; enddo endif !register diagnostics @@ -120,10 +121,32 @@ subroutine MOM_stoch_eos_init(G, Time, param_file, CS, restart_CS, diag) endif endif -end subroutine MOM_stoch_eos_init + ! This module is only used if explicitly enabled or a positive correlation coefficient is set. + MOM_stoch_eos_init = CS%use_stoch_eos .or. (CS%stanley_coeff >= 0.0) + +end function MOM_stoch_eos_init + +!> Register fields related to the stoch_EOS module for resarts +subroutine stoch_EOS_register_restarts(HI, param_file, CS, restart_CS) + type(hor_index_type), intent(in) :: HI !< Horizontal index structure + type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse + type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. + + call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", CS%stanley_coeff, & + "Coefficient correlating the temperature gradient "//& + "and SGS T variance.", units="nondim", default=-1.0, do_not_log=.true.) + + if (CS%stanley_coeff >= 0.0) then + allocate(CS%pattern(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.0) + call register_restart_field(CS%pattern, "stoch_eos_pattern", .false., restart_CS, & + "Random pattern for stoch EOS", "nondim") + endif + +end subroutine stoch_EOS_register_restarts !> Generates a pattern in space and time for the ocean stochastic equation of state -subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS, diag) +subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]. @@ -132,12 +155,14 @@ subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS, diag) real, intent(in) :: delt !< Time step size for AR1 process [T ~> s]. type(time_type), intent(in) :: Time !< Time for stochastic process type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure - type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics ! local variables - integer :: i,j - integer :: yr,mo,dy,hr,mn,sc - real :: phi,ubar,vbar + real :: ubar, vbar ! Averaged velocities [L T-1 ~> m s-1] + real :: phi ! A temporal correlation factor [nondim] + integer :: i, j + + ! Return without doing anything if this capability is not enabled. + if (.not.CS%use_stoch_eos) return call random_2d_constructor(CS%rn_CS, G%HI, Time, CS%seed) call random_2d_norm(CS%rn_CS, G%HI, CS%rgauss) @@ -145,16 +170,28 @@ subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS, diag) ! advance AR(1) do j=G%jsc,G%jec do i=G%isc,G%iec - ubar=0.5*(u(I,j,1)*G%mask2dCu(I,j)+u(I-1,j,1)*G%mask2dCu(I-1,j)) - vbar=0.5*(v(i,J,1)*G%mask2dCv(i,J)+v(i,J-1,1)*G%mask2dCv(i,J-1)) - phi=exp(-delt*CS%tfac*sqrt((ubar**2+vbar**2)*CS%l2_inv(i,j))) - CS%pattern(i,j)=phi*CS%pattern(i,j) + CS%amplitude*sqrt(1-phi**2)*CS%rgauss(i,j) - CS%phi(i,j)=phi + ubar = 0.5*(u(I,j,1)*G%mask2dCu(I,j)+u(I-1,j,1)*G%mask2dCu(I-1,j)) + vbar = 0.5*(v(i,J,1)*G%mask2dCv(i,J)+v(i,J-1,1)*G%mask2dCv(i,J-1)) + phi = exp(-delt*CS%tfac*sqrt((ubar**2+vbar**2)*CS%l2_inv(i,j))) + CS%pattern(i,j) = phi*CS%pattern(i,j) + CS%amplitude*sqrt(1-phi**2)*CS%rgauss(i,j) + CS%phi(i,j) = phi enddo enddo end subroutine MOM_stoch_eos_run +!> Write out any diagnostics related to this module. +subroutine post_stoch_EOS_diags(CS, tv, diag) + type(MOM_stoch_eos_CS), intent(in) :: CS !< Stochastic control structure + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + type(diag_ctrl), intent(inout) :: diag !< Structure to control diagnostics + + if (CS%id_stoch_eos > 0) call post_data(CS%id_stoch_eos, CS%pattern, diag) + if (CS%id_stoch_phi > 0) call post_data(CS%id_stoch_phi, CS%phi, diag) + if (CS%id_tvar_sgs > 0) call post_data(CS%id_tvar_sgs, tv%varT, diag) + +end subroutine post_stoch_EOS_diags + !> Computes a parameterization of the SGS temperature variance subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -171,15 +208,17 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) !! in massless layers filled vertically by diffusion. S !> The filled salinity [S ~> ppt], with the values in !! in massless layers filled vertically by diffusion. - integer :: i, j, k real :: hl(5) !> Copy of local stencil of H [H ~> m] real :: dTdi2, dTdj2 !> Differences in T variance [C2 ~> degC2] + integer :: i, j, k + + ! Nothing happens if a negative correlation coefficient is set. + if (CS%stanley_coeff < 0.0) return ! This block does a thickness weighted variance calculation and helps control for ! extreme gradients along layers which are vanished against topography. It is ! still a poor approximation in the interior when coordinates are strongly tilted. - if (.not. associated(tv%varT)) call safe_alloc_ptr(tv%varT, G%isd, G%ied, G%jsd, G%jed, GV%ke) - + if (.not. associated(tv%varT)) allocate(tv%varT(G%isd:G%ied, G%jsd:G%jed, GV%ke), source=0.0) call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.) do k=1,G%ke @@ -193,12 +232,12 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) ! SGS variance in i-direction [C2 ~> degC2] dTdi2 = ( ( G%mask2dCu(I ,j) * G%IdxCu(I ,j) * ( T(i+1,j,k) - T(i,j,k) ) & - + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) & - ) * G%dxT(i,j) * 0.5 )**2 + + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) & + ) * G%dxT(i,j) * 0.5 )**2 ! SGS variance in j-direction [C2 ~> degC2] dTdj2 = ( ( G%mask2dCv(i,J ) * G%IdyCv(i,J ) * ( T(i,j+1,k) - T(i,j,k) ) & - + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) & - ) * G%dyT(i,j) * 0.5 )**2 + + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) & + ) * G%dyT(i,j) * 0.5 )**2 tv%varT(i,j,k) = CS%stanley_coeff * ( dTdi2 + dTdj2 ) ! Turn off scheme near land tv%varT(i,j,k) = tv%varT(i,j,k) * (minval(hl) / (maxval(hl) + GV%H_subroundoff)) @@ -210,7 +249,7 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) do k=1,G%ke do j=G%jsc,G%jec do i=G%isc,G%iec - tv%varT(i,j,k) = exp (CS%stanley_a * CS%pattern(i,j)) * tv%varT(i,j,k) + tv%varT(i,j,k) = exp(CS%stanley_a * CS%pattern(i,j)) * tv%varT(i,j,k) enddo enddo enddo From eeb07b3c364f1ddae07b80223dfaa1138371a380 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Dec 2022 08:55:49 -0500 Subject: [PATCH 2/6] (*)Correct USE_STANLEY domain extents in EOS calls Corrected the domain extents used in EOS calls that are triggered with USE_STANLEY_GM and USE_STANLEY_ISO. These bugs were accidentally introduced when the changes adding the MOM_stoch_eos code to the main branch of MOM6 were merged with changes on the dev/gfdl branch. Also added a test for cases when USE_STANLEY_GM is set to true but STANLEY_COEF is negative to reset the internal versions of this flag to false with a sensible warning message rather than encountering segmentation faults. All solutions are bitwise identical in cases that worked before. --- src/core/MOM_isopycnal_slopes.F90 | 92 +++++++++---------- .../lateral/MOM_thickness_diffuse.F90 | 38 +++++--- 2 files changed, 71 insertions(+), 59 deletions(-) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index b5bd51d75a..07dd19b0a6 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -3,12 +3,12 @@ module MOM_isopycnal_slopes ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_grid, only : ocean_grid_type -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density_derivs -use MOM_EOS, only : calculate_density_second_derivs +use MOM_debugging, only : hchksum, uvchksum +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : calculate_density_derivs, calculate_density_second_derivs, EOS_domain use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S @@ -28,13 +28,12 @@ module MOM_isopycnal_slopes !> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2 !! and dz*S^2*g-prime used, or calculable from factors used, during the calculation. subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, & - slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC) !, eta_to_m) + slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m] or units - !! given by 1/eta_to_m) + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity @@ -61,15 +60,12 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan integer, optional, intent(in) :: halo !< Halo width over which to compute type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. - ! real, optional, intent(in) :: eta_to_m !< The conversion factor from the units - ! (This argument has been tested but for now serves no purpose.) !! of eta to m; US%Z_to_m by default. ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & T, & ! The temperature [C ~> degC], with the values in ! in massless layers filled vertically by diffusion. - S !, & ! The filled salinity [S ~> ppt], with the values in + S ! The filled salinity [S ~> ppt], with the values in ! in massless layers filled vertically by diffusion. -! Rho ! Density itself, when a nonlinear equation of state is not in use [R ~> kg m-3]. real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & pres ! The pressure at an interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ingored. @@ -96,15 +92,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan T_hr, & ! Temperature on the interface at the h (+1) point [C ~> degC]. S_hr, & ! Salinity on the interface at the h (+1) point [S ~> ppt] pres_hr ! Pressure on the interface at the h (+1) point [R L2 T-2 ~> Pa]. - real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density - real :: drdjA, drdjB ! gradients in the layers above (A) and below (B) the - ! interface times the grid spacing [R ~> kg m-3]. + real :: drdiA, drdiB ! Along layer zonal potential density gradients in the layers above (A) + ! and below (B) the interface times the grid spacing [R ~> kg m-3]. + real :: drdjA, drdjB ! Along layer meridional potential density gradients in the layers above (A) + ! and below (B) the interface times the grid spacing [R ~> kg m-3]. real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3]. real :: hg2A, hg2B ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4]. real :: hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4]. real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2]. real :: dzaL, dzaR ! Temporary thicknesses in eta units [Z ~> m]. - real :: wtA, wtB, wtL, wtR ! Unscaled weights, with various units. + real :: wtA, wtB ! Unnormalized weights of the slopes above and below [H3 ~> m3 or kg3 m-6] + real :: wtL, wtR ! Unnormalized weights of the slopes to the left and right [H3 Z ~> m4 or kg3 m-5] real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4]. real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4]. real :: slope ! The slope of density surfaces, calculated in a way @@ -117,33 +115,34 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! in roundoff and can be neglected [Z ~> m]. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. real :: G_Rho0 ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] - real :: Z_to_L ! A conversion factor between from units for e to the - ! units for lateral distances [L Z-1 ~> 1] - real :: L_to_Z ! A conversion factor between from units for lateral distances - ! to the units for e [Z L-1 ~> 1] - real :: H_to_Z ! A conversion factor from thickness units to the units of e [Z H-1 ~> 1 or m3 kg-1] logical :: present_N2_u, present_N2_v - integer, dimension(2) :: EOSdom_u, EOSdom_v ! Domains for the equation of state calculations at u and v points + logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain. + integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of + ! state calculations at u-points. + integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of + ! state calculations at v-points. + integer, dimension(2) :: EOSdom_h1 ! The shifted i-computational domain to use for equation of + ! state calculations at h points with 1 extra halo point integer :: is, ie, js, je, nz, IsdB integer :: i, j, k integer :: l_seg - logical :: local_open_u_BC, local_open_v_BC if (present(halo)) then is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo + EOSdom_h1(:) = EOS_domain(G%HI, halo=halo+1) else is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + EOSdom_h1(:) = EOS_domain(G%HI, halo=1) endif + EOSdom_u(1) = is-1 - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) + EOSdom_v(:) = EOS_domain(G%HI, halo=halo) + nz = GV%ke ; IsdB = G%IsdB + h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 - Z_to_L = US%Z_to_L ; H_to_Z = GV%H_to_Z - ! if (present(eta_to_m)) then - ! Z_to_L = eta_to_m*US%m_to_L ; H_to_Z = GV%H_to_m / eta_to_m - ! endif - L_to_Z = 1.0 / Z_to_L - dz_neglect = GV%H_subroundoff * H_to_Z + dz_neglect = GV%H_subroundoff * GV%H_to_Z local_open_u_BC = .false. local_open_v_BC = .false. @@ -221,12 +220,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan enddo ; enddo enddo - EOSdom_u(1) = is-1 - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) - !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, & - !$OMP h_neglect,dz_neglect,Z_to_L,L_to_Z,H_to_Z,h_neglect2, & - !$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,local_open_u_BC, & - !$OMP dzu,OBC,use_stanley) & + !$OMP h_neglect,dz_neglect,h_neglect2, & + !$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, & + !$OMP local_open_u_BC,dzu,OBC,use_stanley) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & @@ -259,7 +256,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & call calculate_density_second_derivs(T_h, S_h, pres_h, & scrap, scrap, drho_dT_dT_h, scrap, scrap, & - tv%eqn_of_state, dom=[is-1,ie-is+3]) + tv%eqn_of_state, dom=EOSdom_h1) endif do I=is-1,ie @@ -294,7 +291,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect haR = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect if (GV%Boussinesq) then - dzaL = haL * H_to_Z ; dzaR = haR * H_to_Z + dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z else dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect dzaR = 0.5*(e(i+1,j,K-1) - e(i+1,j,K+1)) + dz_neglect @@ -318,7 +315,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = (Z_to_L*drdx)**2 + drdz**2 + mag_grad2 = (US%Z_to_L*drdx)**2 + drdz**2 if (mag_grad2 > 0.0) then slope = drdx / sqrt(mag_grad2) else ! Just in case mag_grad2 = 0 ever. @@ -351,11 +348,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan enddo ! I enddo ; enddo ! end of j-loop - EOSdom_v(1) = is - (G%isd-1) ; EOSdom_v(2) = ie - (G%isd-1) - ! Calculate the meridional isopycnal slope. !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, & - !$OMP h,h_neglect,e,dz_neglect,Z_to_L,L_to_Z,H_to_Z, & + !$OMP h,h_neglect,e,dz_neglect, & !$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, & !$OMP dzv,local_open_v_BC,OBC,use_stanley) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & @@ -393,10 +388,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & call calculate_density_second_derivs(T_h, S_h, pres_h, & scrap, scrap, drho_dT_dT_h, scrap, scrap, & - tv%eqn_of_state, dom=[is,ie-is+1]) + tv%eqn_of_state, dom=EOSdom_v) call calculate_density_second_derivs(T_hr, S_hr, pres_hr, & scrap, scrap, drho_dT_dT_hr, scrap, scrap, & - tv%eqn_of_state, dom=[is,ie-is+1]) + tv%eqn_of_state, dom=EOSdom_v) endif do i=is,ie if (use_EOS) then @@ -430,7 +425,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect haR = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect if (GV%Boussinesq) then - dzaL = haL * H_to_Z ; dzaR = haR * H_to_Z + dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z else dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect dzaR = 0.5*(e(i,j+1,K-1) - e(i,j+1,K+1)) + dz_neglect @@ -454,7 +449,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = (Z_to_L*drdy)**2 + drdz**2 + mag_grad2 = (US%Z_to_L*drdy)**2 + drdz**2 if (mag_grad2 > 0.0) then slope = drdy / sqrt(mag_grad2) else ! Just in case mag_grad2 = 0 ever. @@ -513,8 +508,9 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar ! Local variables real :: ent(SZI_(G),SZK_(GV)+1) ! The diffusive entrainment (kappa*dt)/dz ! between layers in a timestep [H ~> m or kg m-2]. - real :: b1(SZI_(G)), d1(SZI_(G)) ! b1, c1, and d1 are variables used by the - real :: c1(SZI_(G),SZK_(GV)) ! tridiagonal solver. + real :: b1(SZI_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] + real :: d1(SZI_(G)) ! A variable used by the tridiagonal solver [nondim], d1 = 1 - c1. + real :: c1(SZI_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim]. real :: kap_dt_x2 ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4]. real :: h_neglect ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses. real :: h0 ! A negligible thickness to allow for zero thickness layers without @@ -541,7 +537,7 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar T_f(i,j,k) = T_in(i,j,k) ; S_f(i,j,k) = S_in(i,j,k) enddo ; enddo ; enddo else - !$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr) + !$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr) do j=js,je do i=is,ie ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index b30d24eeaf..b8d5e4c89c 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -637,7 +637,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real, dimension(SZIB_(G)) :: & drho_dT_u, & ! The derivative of density with temperature at u points [R C-1 ~> kg m-3 degC-1] drho_dS_u ! The derivative of density with salinity at u points [R S-1 ~> kg m-3 ppt-1]. - real, dimension(SZIB_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ignored. + real, dimension(SZIB_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() + ! with various units that will be ignored [various] real, dimension(SZI_(G)) :: & drho_dT_v, & ! The derivative of density with temperature at v points [R C-1 ~> kg m-3 degC-1] drho_dS_v, & ! The derivative of density with salinity at v points [R S-1 ~> kg m-3 ppt-1]. @@ -665,9 +666,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real :: PE_release_h ! The amount of potential energy released by GM averaged over an h-cell [L4 Z-1 T-3 ~> m3 s-3] ! The calculation is equal to h * S^2 * N^2 * kappa_GM. real :: I4dt ! 1 / 4 dt [T-1 ~> s-1]. - real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density - real :: drdjA, drdjB ! gradients in the layers above (A) and below(B) the - ! interface times the grid spacing [R ~> kg m-3]. + real :: drdiA, drdiB ! Along layer zonal potential density gradients in the layers above (A) + ! and below (B) the interface times the grid spacing [R ~> kg m-3]. + real :: drdjA, drdjB ! Along layer meridional potential density gradients in the layers above (A) + ! and below (B) the interface times the grid spacing [R ~> kg m-3]. real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3]. real :: drdi_u(SZIB_(G),SZK_(GV)) ! Copy of drdi at u-points [R ~> kg m-3]. real :: drdj_v(SZI_(G),SZK_(GV)) ! Copy of drdj at v-points [R ~> kg m-3]. @@ -729,10 +731,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real :: diag_sfn_unlim_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction before ! applying limiters [H L2 T-1 ~> m3 s-1 or kg s-1] logical :: present_slope_x, present_slope_y, calc_derivatives - integer, dimension(2) :: EOSdom_u ! The shifted i-computational domain to use for equation of + integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of ! state calculations at u-points. - integer, dimension(2) :: EOSdom_v ! The shifted I-computational domain to use for equation of + integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of ! state calculations at v-points. + integer, dimension(2) :: EOSdom_h1 ! The shifted i-computational domain to use for equation of + ! state calculations at h points with 1 extra halo point logical :: use_stanley integer :: is, ie, js, je, nz, IsdB, halo integer :: i, j, k @@ -809,12 +813,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%id_sfn_unlim_y > 0) then ; diag_sfn_unlim_y(:,:,1) = 0.0 ; diag_sfn_unlim_y(:,:,nz+1) = 0.0 ; endif EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) + EOSdom_v(:) = EOS_domain(G%HI) + EOSdom_h1(:) = EOS_domain(G%HI, halo=1) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & !$OMP h_neglect2,int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, & !$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1,diag_sfn_x, & - !$OMP diag_sfn_unlim_x,N2_floor,EOSdom_u,use_stanley, Tsgs2, & + !$OMP diag_sfn_unlim_x,N2_floor,EOSdom_u,EOSdom_h1,use_stanley,Tsgs2, & !$OMP present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & @@ -855,7 +861,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & call calculate_density_second_derivs(T_h, S_h, pres_h, & scrap, scrap, drho_dT_dT_h, scrap, scrap, & - tv%eqn_of_state, dom=[is-1,ie-is+3]) + tv%eqn_of_state, EOSdom_h1) endif do I=is-1,ie @@ -1085,7 +1091,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ! end of j-loop ! Calculate the meridional fluxes and gradients. - EOSdom_v(:) = EOS_domain(G%HI) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & @@ -1134,10 +1139,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & call calculate_density_second_derivs(T_h, S_h, pres_h, & scrap, scrap, drho_dT_dT_h, scrap, scrap, & - tv%eqn_of_state, dom=[is,ie-is+1]) + tv%eqn_of_state, EOSdom_v) call calculate_density_second_derivs(T_hr, S_hr, pres_hr, & scrap, scrap, drho_dT_dT_hr, scrap, scrap, & - tv%eqn_of_state, dom=[is,ie-is+1]) + tv%eqn_of_state, EOSdom_v) endif do i=is,ie if (calc_derivatives) then @@ -1971,6 +1976,8 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) real :: strat_floor ! A floor for buoyancy frequency in the Ferrari et al. 2010, ! streamfunction formulation, expressed as a fraction of planetary ! rotation [nondim]. + real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale + ! temperature variance [nondim] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. logical :: MEKE_GEOM_answers_2018 ! If true, use expressions in the MEKE_GEOMETRIC calculation @@ -2096,6 +2103,15 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) call get_param(param_file, mdl, "USE_STANLEY_GM", CS%use_stanley_gm, & "If true, turn on Stanley SGS T variance parameterization "// & "in GM code.", default=.false.) + if (CS%use_stanley_gm) then + call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, & + "Coefficient correlating the temperature gradient and SGS T variance.", & + units="nondim", default=-1.0, do_not_log=.true.) + if (Stanley_coeff < 0.0) then + call MOM_error(WARNING, "STANLEY_COEFF must be set >= 0 if USE_STANLEY_GM is true.") + CS%use_stanley_gm = .false. + endif + endif call get_param(param_file, mdl, "OMEGA", omega, & "The rotation rate of the earth.", & default=7.2921e-5, units="s-1", scale=US%T_to_s, do_not_log=.not.CS%use_FGNV_streamfn) From d4cb1d43e6c3ed28d2dae881d073c7dcb07bb6fb Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Dec 2022 08:56:22 -0500 Subject: [PATCH 3/6] (*)Correct scaling of USE_STANLEY_PGF diagnostics Corrected the diagnostics rho_pgf, rho_stanley_pgf and p_stanley to only be calculated if they would be written, give identical output when dimensional rescaling is applied, and be documented with the right units. Also added a test for cases when USE_STANLEY_GM is set to true but STANLEY_COEF is negative to reset the internal versions of this flag to false with a sensible warning message rather than encountering segmentation faults. All solutions are bitwise identical in cases that worked before, but there are changes in some diagnostics when they are dimensionally rescaled. --- src/core/MOM_PressureForce_FV.F90 | 86 +++++++++++++++++++++---------- 1 file changed, 59 insertions(+), 27 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 5fdc4a1182..854b6b788c 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -13,7 +13,7 @@ module MOM_PressureForce_FV use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain use MOM_density_integrals, only : int_density_dz, int_specific_vol_dp use MOM_density_integrals, only : int_density_dz_generic_plm, int_density_dz_generic_ppm use MOM_density_integrals, only : int_spec_vol_dp_generic_plm @@ -477,12 +477,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm T_t, T_b ! Top and bottom edge values for linear reconstructions ! of temperature within each layer [C ~> degC]. real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & - rho_pgf, rho_stanley_pgf ! Density [kg m-3] from EOS with and without SGS T variance - ! in Stanley parameterization. + rho_pgf, rho_stanley_pgf ! Density [R ~> kg m-3] from EOS with and without SGS T variance + ! in Stanley parameterization. real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & - p_stanley ! Pressure [Pa] estimated with Rho_0 - real :: rho_stanley_scalar ! Scalar quantity to hold density [kg m-3] in Stanley diagnostics. - real :: p_stanley_scalar ! Scalar quantity to hold pressure [Pa] in Stanley diagnostics. + p_stanley ! Pressure [R L2 T-2 ~> Pa] estimated with Rho_0 + real :: zeros(SZI_(G)) ! An array of zero values that can be used as an argument [various] real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3]. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). @@ -493,12 +492,15 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real :: G_Rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. real :: rho_ref ! The reference density [R ~> kg m-3]. real :: dz_neglect ! A minimal thickness [Z ~> m], like e. + real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure + ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]. logical :: use_p_atm ! If true, use the atmospheric pressure. logical :: use_ALE ! If true, use an ALE pressure reconstruction. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. real, parameter :: C1_6 = 1.0/6.0 ! [nondim] integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_h ! The i-computational domain for the equation of state at tracer points integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: i, j, k @@ -759,25 +761,43 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif if (CS%use_stanley_pgf) then - do j=js,je ; do i=is,ie ; - p_stanley_scalar=0.0 - do k=1, nz - p_stanley_scalar = p_stanley_scalar + 0.5 * h(i,j,k) * GV%H_to_Pa !Pressure at mid-point of layer - call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley_scalar, 0.0, 0.0, 0.0, & - rho_stanley_scalar, tv%eqn_of_state) - rho_pgf(i,j,k) = rho_stanley_scalar - call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley_scalar, tv%varT(i,j,k), 0.0, 0.0, & - rho_stanley_scalar, tv%eqn_of_state) - rho_stanley_pgf(i,j,k) = rho_stanley_scalar - p_stanley(i,j,k) = p_stanley_scalar - p_stanley_scalar = p_stanley_scalar + 0.5 * h(i,j,k) * GV%H_to_Pa !Pressure at bottom of layer - enddo; enddo; enddo - endif + ! Calculated diagnostics related to the Stanley parameterization + zeros(:) = 0.0 + EOSdom_h(:) = EOS_domain(G%HI) + if ((CS%id_p_stanley>0) .or. (CS%id_rho_pgf>0) .or. (CS%id_rho_stanley_pgf>0)) then + ! Find the pressure at the mid-point of each layer. + H_to_RL2_T2 = GV%g_Earth*GV%H_to_RZ + if (use_p_atm) then + do j=js,je ; do i=is,ie + p_stanley(i,j,1) = 0.5*h(i,j,1) * H_to_RL2_T2 + p_atm(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + p_stanley(i,j,1) = 0.5*h(i,j,1) * H_to_RL2_T2 + enddo ; enddo + endif + do k=2,nz ; do j=js,je ; do i=is,ie + p_stanley(i,j,k) = p_stanley(i,j,k-1) + 0.5*(h(i,j,k-1) + h(i,j,k)) * H_to_RL2_T2 + enddo ; enddo ; enddo + endif + if (CS%id_p_stanley>0) call post_data(CS%id_p_stanley, p_stanley, CS%diag) + if (CS%id_rho_pgf>0) then + do k=1,nz ; do j=js,je + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_stanley(:,j,k), zeros, & + zeros, zeros, rho_pgf(:,j,k), tv%eqn_of_state, EOSdom_h) + enddo ; enddo + call post_data(CS%id_rho_pgf, rho_pgf, CS%diag) + endif + if (CS%id_rho_stanley_pgf>0) then + do k=1,nz ; do j=js,je + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_stanley(:,j,k), tv%varT(:,j,k), & + zeros, zeros, rho_stanley_pgf(:,j,k), tv%eqn_of_state, EOSdom_h) + enddo ; enddo + call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag) + endif + endif if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag) - if (CS%id_rho_pgf>0) call post_data(CS%id_rho_pgf, rho_pgf, CS%diag) - if (CS%id_rho_stanley_pgf>0) call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag) - if (CS%id_p_stanley>0) call post_data(CS%id_p_stanley, p_stanley, CS%diag) end subroutine PressureForce_FV_Bouss @@ -791,10 +811,14 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(PressureForce_FV_CS), intent(inout) :: CS !< Finite volume PGF control structure type(tidal_forcing_CS), intent(in), target, optional :: tides_CSp !< Tides control structure + + ! Local variables + real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale + ! temperature variance [nondim] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl ! This module's name. - logical :: use_ALE + logical :: use_ALE ! If true, use the Vertical Lagrangian Remap algorithm CS%initialized = .true. CS%diag => diag ; CS%Time => Time @@ -842,12 +866,20 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS "If true, turn on Stanley SGS T variance parameterization "// & "in PGF code.", default=.false.) if (CS%use_stanley_pgf) then + call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, & + "Coefficient correlating the temperature gradient and SGS T variance.", & + units="nondim", default=-1.0, do_not_log=.true.) + if (Stanley_coeff < 0.0) then + call MOM_error(WARNING, "STANLEY_COEFF must be set >= 0 if USE_STANLEY_PGF is true.") + CS%use_stanley_pgf = .false. + endif + CS%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, & - Time, 'rho in PGF', 'kg m3') + Time, 'rho in PGF', 'kg m-3', conversion=US%R_to_kg_m3) CS%id_rho_stanley_pgf = register_diag_field('ocean_model', 'rho_stanley_pgf', diag%axesTL, & - Time, 'rho in PGF with Stanley correction', 'kg m3') + Time, 'rho in PGF with Stanley correction', 'kg m-3', conversion=US%R_to_kg_m3) CS%id_p_stanley = register_diag_field('ocean_model', 'p_stanley', diag%axesTL, & - Time, 'p in PGF with Stanley correction', 'Pa') + Time, 'p in PGF with Stanley correction', 'Pa', conversion=US%RL2_T2_to_Pa) endif if (CS%tides) then CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, & From e18e2dadca9c5fa5892ae91c06b325c1a561348b Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Dec 2022 08:57:43 -0500 Subject: [PATCH 4/6] (*)Test for inconsistent USE_STANLEY flags Also added tests for cases when USE_STANLEY_ISO or USE_STANLEY_ML are set to true but STANLEY_COEF is negative to reset the internal versions of these flags to false with a sensible warning message rather than encountering segmentation faults. All solutions are bitwise identical in cases that worked before. --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 11 +++++++++++ .../lateral/MOM_mixed_layer_restrat.F90 | 11 +++++++++++ 2 files changed, 22 insertions(+) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index a8928ef06c..1aede30a74 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1114,6 +1114,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) ! scaled by the resolution function. logical :: better_speed_est ! If true, use a more robust estimate of the first ! mode wave speed as the starting point for iterations. + real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale + ! temperature variance [nondim] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name. @@ -1208,6 +1210,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "USE_STANLEY_ISO", CS%use_stanley_iso, & "If true, turn on Stanley SGS T variance parameterization "// & "in isopycnal slope code.", default=.false.) + if (CS%use_stanley_iso) then + call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, & + "Coefficient correlating the temperature gradient and SGS T variance.", & + units="nondim", default=-1.0, do_not_log=.true.) + if (Stanley_coeff < 0.0) then + call MOM_error(WARNING, "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.") + CS%use_stanley_iso = .false. + endif + endif if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct) then in_use = .true. diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 08c29c6c9e..5d87761363 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -854,6 +854,8 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, real :: flux_to_kg_per_s ! A unit conversion factor for fluxes. [kg T s-1 H-1 L-2 ~> kg m-3 or 1] real :: omega ! The Earth's rotation rate [T-1 ~> s-1]. real :: ustar_min_dflt ! The default value for RESTRAT_USTAR_MIN [Z T-1 ~> m s-1] + real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale + ! temperature variance [nondim] ! This include declares and sets the variable "version". # include "version_variable.h" integer :: i, j @@ -891,6 +893,15 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_stanley_ml, & "If true, turn on Stanley SGS T variance parameterization "// & "in ML restrat code.", default=.false.) + if (CS%use_stanley_ml) then + call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, & + "Coefficient correlating the temperature gradient and SGS T variance.", & + units="nondim", default=-1.0, do_not_log=.true.) + if (Stanley_coeff < 0.0) then + call MOM_error(WARNING, "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ML is true.") + CS%use_stanley_ml = .false. + endif + endif call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, & 'The value the von Karman constant as used for mixed layer viscosity.', & units='nondim', default=0.41) From 5b2c5846abbf7dfa0ebd5faf91acf6e935f5f0d9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Dec 2022 09:13:07 -0500 Subject: [PATCH 5/6] (*)Correct scaling of DT_OBC_SEG_UPDATE_OBGC Added a missing dimensional rescaling factor in the get_param call for the recently added variable DT_OBC_SEG_UPDATE_OBGC in initialize_MOM, which would have caused certain cases to fail the dimensional consistency tests. All answers are bitwise identical in cases that do not use dimensional rescaling, but answers will change (and be corrected) in some cases that use dimensional consistency tests. --- src/core/MOM.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 8d7caf83ed..2288c007c8 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -300,7 +300,10 @@ module MOM !! calculated, and if it is 0, dtbt is calculated every step. type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period. type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated. - real :: dt_obc_seg_period !< The time interval between OBC segment updates for OBGC tracers + real :: dt_obc_seg_period !< The time interval between OBC segment updates for OBGC + !! tracers [T ~> s], or a negative value if the segment + !! data are time-invarant, or zero to update the OBGC + !! segment data with every call to update_OBC_segment_data. type(time_type) :: dt_obc_seg_interval !< A time_time representation of dt_obc_seg_period. type(time_type) :: dt_obc_seg_time !< The next time OBC segment update is applied to OBGC tracers. @@ -2186,12 +2189,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & units="s", default=default_val, do_not_read=(dtbt > 0.0)) endif - CS%dt_obc_seg_period = -1.0 call get_param(param_file, "MOM", "DT_OBC_SEG_UPDATE_OBGC", CS%dt_obc_seg_period, & "The time between OBC segment data updates for OBGC tracers. "//& "This must be an integer multiple of DT and DT_THERM. "//& "The default is set to DT.", & - units="s", default=US%T_to_s*CS%dt, do_not_log=.not.associated(CS%OBC)) + units="s", default=US%T_to_s*CS%dt, scale=US%s_to_T, do_not_log=.not.associated(CS%OBC)) ! This is here in case these values are used inappropriately. use_frazil = .false. ; bound_salinity = .false. From b8c3845d8e4068592bb9aa1fcd33ab99cadc7e4f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Dec 2022 10:08:02 -0500 Subject: [PATCH 6/6] Rescale interface heights for MOM_IC Use dimensionally rescaled units when preparing fields to write to the MOM_IC, and then use the conversion argument to the register_restart_field call to undue this scaling, following the pattern for other calls. All answers and output are bitwise identical. --- src/core/MOM.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2288c007c8..9d0701ec2a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3217,7 +3217,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing ! various unit conversion factors type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() - real, allocatable :: z_interface(:,:,:) ! Interface heights [m] + real, allocatable :: z_interface(:,:,:) ! Interface heights [Z ~> m] call cpu_clock_begin(id_clock_init) call callTree_enter("finish_MOM_initialization()") @@ -3240,9 +3240,9 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) restart_CSp_tmp = restart_CSp call restart_registry_lock(restart_CSp_tmp, unlocked=.true.) allocate(z_interface(SZI_(G),SZJ_(G),SZK_(GV)+1)) - call find_eta(CS%h, CS%tv, G, GV, US, z_interface, eta_to_m=1.0, dZref=G%Z_ref) + call find_eta(CS%h, CS%tv, G, GV, US, z_interface, dZref=G%Z_ref) call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, & - "Interface heights", "meter", z_grid='i') + "Interface heights", "meter", z_grid='i', conversion=US%Z_to_m) ! NOTE: write_ic=.true. routes routine to fms2 IO write_initial_conditions interface call save_restart(dirs%output_directory, Time, CS%G_in, & restart_CSp_tmp, filename=CS%IC_file, GV=GV, write_ic=.true.)