From 8f97b3d2f73f80ad4dd5c0ac30cca0c1fad8e435 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Sat, 7 May 2022 12:00:36 -0600 Subject: [PATCH] Corrections to cell-averaged density computation (#213) * initial hooks for stochastic EOS modifications * remove debug statements * add documentation * Change ampltiude from 0.39 to sqrt(.39) * remove global_indexing logic from stoch_eos_init * switch to using MOM_random and add restart capability * update random sequence to update each each time-step * remove tseed0 from MOM_random (leftover from debugging) * Added necessary submodules and S^2, T^2 diagnostics to MOM_diagnostics * Added diagnostics for outputting variables related to the stochastic parameterization. * Diagnostics in MOM_PressureForce_FV updated for stochastic (rather than deterministic) Stanley SGS T variance parameterization. * Added parentheses for reproducibility. * Changed diagnostics to account for possible absence of stoch_eos_pattern in MOM_PressureForce_FV, when deterministic parameterization is on. * remove mom6_da_hooks and geoKdTree from pkg * add stochastic compoment to MOM_thickness_diffuse * fix array size declaration and post_data * Corrected indexing of loops in MOM_calc_varT * Changed how parameterization of SGS T variance (deterministic and stochastic) is switched on in PGF and thickness diffusion codes * Corrected a few typos * Cleaned up indices, redundant diagnostic, printing * Fixed diagnostic IDs * Fixed diagnostics typo * Corrected indices in calculation of tv%varT * Minor index fix * Corrected bug in pressure in Stanley diagnostics * Fixed whitespace error * Stoch eos clock (#5) *Added a clock for the Stanley parameterization Co-authored-by: jkenigson * add halo update to random pattern * Update MOM_stoch_eos.F90 Fix bug for looping over compute domain (is -> isc etc.) * Avoid unnessary computations on halo (MOM_stoch_eos) and code clean-up (MOM_thickness_diffuse) * Removed halo updates before determ param calc * Update MOM_stoch_eos.F90 Removed unnecessary code * Bug - indices are transposed * Changed Stanley stochastic coefficient from exp(X) to exp(aX) (#9) * Changed Stanley stochastic coefficient from exp(X) to exp(aX) * Extra spaces removed * Stoch eos init fix (#10) * Don't bother calculating tv%varT if stanley_coeff<0 * Missing then added * Merge Ian Grooms Tvar Discretization (#11) * Update MOM_stoch_eos.F90 In progress updating stencil for$ | dx \times \nabla T|^2$ calculation * New discretization of |dx\circ\nablaT|^2 Co-authored-by: Ian Grooms * Multiplied tvar%SGS by grid cell thickness ratio * Added limiter for tv%varT * Stoch eos ncar linear disc (#12) * Update MOM_stoch_eos.F90 In progress updating stencil for$ | dx \times \nabla T|^2$ calculation * New discretization of |dx\circ\nablaT|^2 * AR1 timescale land mask Adds land mask to the computation of the AR1 decorrelation time * Update dt in call to MOM_stoch_eos_run The call to `MOM_stoch_eos_run` (which time steps the noise) is from within `step_MOM_dynamics`. `step_MOM_dynamics` advances on time step `dt` (per line 957), but the noise is updated using `dt_thermo`. It seems more appropriate to update the noise using `dt`, since it gets called from within `step_MOM_dynamics`. * Fixed the units for r_sm_H * Remove vestigial declarations The variables `hl`, `Tl`, `mn_T`, `mn_T2`, and `r_sm_H` are no longer used, so I removed their declarations and an OMP private clause Co-authored-by: Ian Grooms * Update MOM_thickness_diffuse.F90 Changed index for soft convention * Update CVMix-src * Ensure use_varT, etc., initialized * Don't register stanley diagnostics if scheme is off * Stanley density second derivs at h pts (#15) * Change discretization of Stanley correction (drho_dT_dT at h points) * Limit Stanley noise, shrink limiting value * Revert t variance discretization * Reverted variable declarations * Stanley scheme in mixed_layer_restrat, vert_fill in stoch_eos, code cleanup (#19) * Test Stanley EOS param in mixed_layer_restrat * Fix size of TS cov, S var in Stanley calculate_density calls * Test move stanley scheme initialization * Added missing openMP directives * Revert Stanley tvar discretization (#18) * Perform vertical filling in calculation of T variance * Variable declaration syntax error, remove scaling from get_param * Fix call to vert_fill_TS * Code cleanup, whitespace cleanup Co-authored-by: Jessica Kenigson * Use Stanley (2020) variance; scheme off at coast * Comment clean-up * Remove factor of 0.5 in Tvar * Don't calculate Stanley diagnostics on halo * Change start indices in stanley_density_1d * Stanley param in MOM_isopycnal_slopes (#22) Stanley param in MOM_isopycnal_slopes and thickness diffuse index fix * Set eady flag to true if use_stored_slopes is true * Cleanup, docs, whitespace * Docs and whitespace * Docs and whitespace * Docs and whitespace * Whitespace cleanup * Whitespace cleanup * Clean up whitespace * Docs cleanup * use_stanley * Update MOM_lateral_mixing_coeffs.F90 * Adds link to another TEOS10 module * Set Stanley off for testing * Line continuation Co-authored-by: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Co-authored-by: Philip Pegion Co-authored-by: Jessica Kenigson Co-authored-by: Jessica Kenigson Co-authored-by: jkenigson Co-authored-by: jskenigson Co-authored-by: Jessica Kenigson Co-authored-by: Jessica Kenigson Co-authored-by: Philip Pegion Co-authored-by: Jessica Kenigson --- .testing/tc2/MOM_input | 4 +- src/core/MOM.F90 | 22 +- src/core/MOM_PressureForce_FV.F90 | 109 ++++----- src/core/MOM_density_integrals.F90 | 37 +-- src/core/MOM_isopycnal_slopes.F90 | 72 +++++- src/core/MOM_stoch_eos.F90 | 212 ++++++++++++++++++ src/diagnostics/MOM_diagnostics.F90 | 37 ++- src/equation_of_state/MOM_EOS.F90 | 18 +- .../TEOS10/gsw_mod_error_functions.f90 | 1 + src/framework/MOM_random.F90 | 1 + .../lateral/MOM_lateral_mixing_coeffs.F90 | 9 +- .../lateral/MOM_mixed_layer_restrat.F90 | 34 ++- .../lateral/MOM_thickness_diffuse.F90 | 141 +++++------- 13 files changed, 515 insertions(+), 182 deletions(-) create mode 100644 src/core/MOM_stoch_eos.F90 create mode 120000 src/equation_of_state/TEOS10/gsw_mod_error_functions.f90 diff --git a/.testing/tc2/MOM_input b/.testing/tc2/MOM_input index ca84d1c382..c7d2a35aa6 100644 --- a/.testing/tc2/MOM_input +++ b/.testing/tc2/MOM_input @@ -297,7 +297,7 @@ BOUND_CORIOLIS = True ! [Boolean] default = False ! v-points, and similarly at v-points. This option would ! have no effect on the SADOURNY Coriolis scheme if it ! were possible to use centered difference thickness fluxes. -PGF_STANLEY_T2_DET_COEFF = 0.5 ! [nondim] default = -1.0 +PGF_STANLEY_T2_DET_COEFF = -1.0 ! [nondim] default = -1.0 ! The coefficient correlating SGS temperature variance with the mean temperature ! gradient in the deterministic part of the Stanley form of the Brankart ! correction. Negative values disable the scheme. @@ -430,7 +430,7 @@ KHTH = 1.0 ! [m2 s-1] default = 0.0 ! The background horizontal thickness diffusivity. KHTH_MAX = 900.0 ! [m2 s-1] default = 0.0 ! The maximum horizontal thickness diffusivity. -STANLEY_PRM_DET_COEFF = 0.5 ! [nondim] default = -1.0 +STANLEY_PRM_DET_COEFF = -1.0 ! [nondim] default = -1.0 ! The coefficient correlating SGS temperature variance with the mean temperature ! gradient in the deterministic part of the Stanley parameterization. Negative ! values disable the scheme. diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e2e966e3b4..c14085c923 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -108,6 +108,7 @@ 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_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 @@ -242,6 +243,7 @@ module MOM logical :: use_ALE_algorithm !< If true, use the ALE algorithm rather than layered !! isopycnal/stacked shallow water mode. This logical is set by calling the !! function useRegridding() from the MOM_regridding module. + type(MOM_stoch_eos_CS) :: stoch_eos_CS !< structure containing random pattern for stoch EOS logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction !! updates occur first in directionally split parts of the calculation. real :: first_dir_restart = -1.0 !< A real copy of G%first_direction for use in restart files @@ -444,6 +446,8 @@ module MOM integer :: id_clock_other integer :: id_clock_offline_tracer integer :: id_clock_unit_tests +integer :: id_clock_stoch +integer :: id_clock_varT !>@} contains @@ -1063,6 +1067,15 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & showCallTree = callTree_showQuery() 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) + 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) + endif + call cpu_clock_end(id_clock_varT) if ((CS%t_dyn_rel_adv == 0.0) .and. CS%thickness_diffuse .and. CS%thickness_diffuse_first) then @@ -1229,6 +1242,9 @@ 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) call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other) @@ -2765,6 +2781,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call set_visc_init(Time, G, GV, US, param_file, diag, CS%visc, CS%set_visc_CSp, restart_CSp, CS%OBC) call thickness_diffuse_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%thickness_diffuse_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 (CS%split) then allocate(eta(SZI_(G),SZJ_(G)), source=0.0) call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, & @@ -2854,7 +2872,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif ! This subroutine initializes any tracer packages. - new_sim = is_new_run(restart_CSp) call tracer_flow_control_init(.not.new_sim, Time, G, GV, US, CS%h, param_file, & CS%diag, CS%OBC, CS%tracer_flow_CSp, CS%sponge_CSp, & CS%ALE_sponge_CSp, CS%tv) @@ -3065,6 +3082,7 @@ subroutine register_diags(Time, G, GV, US, IDs, diag) v_extensive=.true.) IDs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, & Time, 'Instantaneous Sea Surface Height', 'm', conversion=US%Z_to_m) + end subroutine register_diags !> Set up CPU clock IDs for timing various subroutines. @@ -3097,6 +3115,8 @@ subroutine MOM_timing_init(CS) if (CS%offline_tracer_mode) then id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=CLOCK_SUBCOMPONENT) endif + id_clock_stoch = cpu_clock_id('(Stochastic EOS)', grain=CLOCK_MODULE) + id_clock_varT = cpu_clock_id('(SGS Temperature Variance)', grain=CLOCK_MODULE) end subroutine MOM_timing_init diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 2a79486a5f..d50ce4b364 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -58,11 +58,12 @@ module MOM_PressureForce_FV integer :: Recon_Scheme !< Order of the polynomial of the reconstruction of T & S !! for the finite volume pressure gradient calculation. !! By the default (1) is for a piecewise linear method - real :: Stanley_T2_det_coeff !< The coefficient correlating SGS temperature variance with - !! the mean temperature gradient in the deterministic part of - !! the Stanley form of the Brankart correction. + + logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF integer :: id_e_tidal = -1 !< Diagnostic identifier - integer :: id_tvar_sgs = -1 !< Diagnostic identifier + integer :: id_rho_pgf = -1 !< Diagnostic identifier + integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier + integer :: id_p_stanley = -1 !< Diagnostic identifier type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure end type PressureForce_FV_CS @@ -167,7 +168,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.") - if (CS%Stanley_T2_det_coeff>=0.) call MOM_error(FATAL, & + if (CS%use_stanley_pgf) call MOM_error(FATAL, & "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//& "implemented in non-Boussinesq mode.") @@ -473,6 +474,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & S_t, S_b, T_t, T_b ! Top and bottom edge values for linear reconstructions ! of salinity and temperature within each layer. + 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. + 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. 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). @@ -487,11 +495,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm 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 :: Tl(5) ! copy and T in local stencil [degC] - real :: mn_T ! mean of T in local stencil [degC] - real :: mn_T2 ! mean of T**2 in local stencil [degC2] - real :: hl(5) ! Copy of local stencil of H [H ~> m] - real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1] real, parameter :: C1_6 = 1.0/6.0 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb @@ -511,49 +514,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm use_ALE = .false. if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS - if (CS%Stanley_T2_det_coeff>=0.) then - if (.not. associated(tv%varT)) call safe_alloc_ptr(tv%varT, G%isd, G%ied, G%jsd, G%jed, GV%ke) - do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 - ! Strictly speaking we should estimate the *horizontal* grid-scale variance - ! but neither of the following blocks make a rotation to the horizontal - ! and instead work along coordinate. - - ! This block calculates a simple |delta T| along coordinates and does - ! not allow vanishing layer thicknesses or layers tracking topography - !! SGS variance in i-direction [degC2] - !dTdi2 = ( ( G%mask2dCu(I ,j) * G%IdxCu(I ,j) * ( tv%T(i+1,j,k) - tv%T(i,j,k) ) & - ! + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( tv%T(i,j,k) - tv%T(i-1,j,k) ) & - ! ) * G%dxT(i,j) * 0.5 )**2 - !! SGS variance in j-direction [degC2] - !dTdj2 = ( ( G%mask2dCv(i,J ) * G%IdyCv(i,J ) * ( tv%T(i,j+1,k) - tv%T(i,j,k) ) & - ! + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( tv%T(i,j,k) - tv%T(i,j-1,k) ) & - ! ) * G%dyT(i,j) * 0.5 )**2 - !tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * 0.5 * ( dTdi2 + dTdj2 ) - - ! 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. - hl(1) = h(i,j,k) * G%mask2dT(i,j) - hl(2) = h(i-1,j,k) * G%mask2dCu(I-1,j) - hl(3) = h(i+1,j,k) * G%mask2dCu(I,j) - hl(4) = h(i,j-1,k) * G%mask2dCv(i,J-1) - hl(5) = h(i,j+1,k) * G%mask2dCv(i,J) - r_sm_H = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + GV%H_subroundoff ) - ! Mean of T - Tl(1) = tv%T(i,j,k) ; Tl(2) = tv%T(i-1,j,k) ; Tl(3) = tv%T(i+1,j,k) - Tl(4) = tv%T(i,j-1,k) ; Tl(5) = tv%T(i,j+1,k) - mn_T = ( hl(1)*Tl(1) + ( ( hl(2)*Tl(2) + hl(3)*Tl(3) ) + ( hl(4)*Tl(4) + hl(5)*Tl(5) ) ) ) * r_sm_H - ! Adjust T vectors to have zero mean - Tl(:) = Tl(:) - mn_T ; mn_T = 0. - ! Variance of T - mn_T2 = ( hl(1)*Tl(1)*Tl(1) + ( ( hl(2)*Tl(2)*Tl(2) + hl(3)*Tl(3)*Tl(3) ) & - + ( hl(4)*Tl(4)*Tl(4) + hl(5)*Tl(5)*Tl(5) ) ) ) * r_sm_H - ! Variance should be positive but round-off can violate this. Calculating - ! variance directly would fix this but requires more operations. - tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * max(0., mn_T2) - enddo ; enddo ; enddo - endif - h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff * GV%H_to_Z I_Rho0 = 1.0 / GV%Rho0 @@ -690,7 +650,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm do k=1,nz ! Calculate 4 integrals through the layer that are required in the ! subsequent calculation. - if (use_EOS) then ! The following routine computes the integrals that are needed to ! calculate the pressure gradient force. Linear profiles for T and S are @@ -701,13 +660,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if ( CS%Recon_Scheme == 1 ) then call int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & - G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, & + G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, intx_dpa, inty_dpa, & useMassWghtInterp=CS%useMassWghtInterp, & use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & - G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, & + G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, intx_dpa, inty_dpa, & useMassWghtInterp=CS%useMassWghtInterp, Z_0p=G%Z_ref) endif else @@ -797,8 +756,26 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif 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 + if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag) - if (CS%id_tvar_sgs>0) call post_data(CS%id_tvar_sgs, tv%varT, 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 @@ -859,14 +836,16 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS "boundary cells is extrapolated, rather than using PCM "//& "in these cells. If true, the same order polynomial is "//& "used as is used for the interior cells.", default=.true.) - call get_param(param_file, mdl, "PGF_STANLEY_T2_DET_COEFF", CS%Stanley_T2_det_coeff, & - "The coefficient correlating SGS temperature variance with "// & - "the mean temperature gradient in the deterministic part of "// & - "the Stanley form of the Brankart correction. "// & - "Negative values disable the scheme.", units="nondim", default=-1.0) - if (CS%Stanley_T2_det_coeff>=0.) then - CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs_pgf', diag%axesTL, & - Time, 'SGS temperature variance used in PGF', 'degC2') + call get_param(param_file, mdl, "USE_STANLEY_PGF", CS%use_stanley_pgf, & + "If true, turn on Stanley SGS T variance parameterization "// & + "in PGF code.", default=.false.) + if (CS%use_stanley_pgf) then + CS%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, & + Time, 'rho in PGF', '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') + CS%id_p_stanley = register_diag_field('ocean_model', 'p_stanley', diag%axesTL, & + Time, 'p in PGF with Stanley correction', 'Pa') endif if (CS%tides) then CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, & diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 8f26918253..c4420e0541 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -338,7 +338,7 @@ end subroutine int_density_dz_generic_pcm !> Compute pressure gradient force integrals by quadrature for the case where !! T and S are linear profiles. subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & - rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, dpa, & + rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, dpa, & intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, & use_inaccurate_form, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for @@ -365,6 +365,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] type(EOS_type), intent(in) :: EOS !< Equation of state structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + logical, intent(in) :: use_stanley_eos !< If true, turn on Stanley SGS T variance parameterization real, dimension(SZI_(HI),SZJ_(HI)), & intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa] real, dimension(SZI_(HI),SZJ_(HI)), & @@ -436,7 +437,6 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & real :: hWght ! A topographically limited thickness weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] - logical :: use_stanley_eos ! True is SGS variance fields exist in tv. logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation ! of density anomalies. logical :: use_varT, use_varS, use_covarTS ! Logicals for SGS variances fields @@ -459,10 +459,15 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form endif - use_varT = associated(tv%varT) - use_covarTS = associated(tv%covarTS) - use_varS = associated(tv%varS) - use_stanley_eos = use_varT .or. use_covarTS .or. use_varS + use_varT = .false. !ensure initialized + use_covarTS = .false. + use_varS = .false. + if (use_stanley_eos) then + use_varT = associated(tv%varT) + use_covarTS = associated(tv%covarTS) + use_varS = associated(tv%varS) + endif + T25(:) = 0. TS5(:) = 0. S25(:) = 0. @@ -489,7 +494,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & if (use_covarTS) TS5(i*5+1:i*5+5) = tv%covarTS(i,j,k) if (use_varS) S25(i*5+1:i*5+5) = tv%varS(i,j,k) enddo - if (use_Stanley_eos) then + if (use_stanley_eos) then if (rho_scale /= 1.0) then call calculate_density(T5, S5, p5, T25, TS5, S25, r5, 1, (ieq-isq+2)*5, EOS, & rho_ref=rho_ref_mks, scale=rho_scale) @@ -781,7 +786,7 @@ end subroutine int_density_dz_generic_plm !> Compute pressure gradient force integrals for layer "k" and the case where T and S !! are parabolic profiles subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & - rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, & + rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, & dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays @@ -807,6 +812,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] type(EOS_type), intent(in) :: EOS !< Equation of state structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + logical, intent(in) :: use_stanley_eos !< If true, turn on Stanley SGS T variance parameterization real, dimension(SZI_(HI),SZJ_(HI)), & intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa] real, dimension(SZI_(HI),SZJ_(HI)), & @@ -868,7 +874,6 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n logical :: use_PPM ! If false, assume zero curvature in reconstruction, i.e. PLM - logical :: use_stanley_eos ! True is SGS variance fields exist in tv. logical :: use_varT, use_varS, use_covarTS Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB @@ -888,10 +893,15 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & t6 = 0. use_PPM = .true. ! This is a place-holder to allow later re-use of this function - use_varT = associated(tv%varT) - use_covarTS = associated(tv%covarTS) - use_varS = associated(tv%varS) - use_stanley_eos = use_varT .or. use_covarTS .or. use_varS + use_varT = .false. !ensure initialized + use_covarTS = .false. + use_varS = .false. + if (use_stanley_eos) then + use_varT = associated(tv%varT) + use_covarTS = associated(tv%covarTS) + use_varS = associated(tv%varS) + endif + T25(:) = 0. TS5(:) = 0. S25(:) = 0. @@ -1003,7 +1013,6 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & S5(n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) ) T5(n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) ) enddo - if (use_stanley_eos) then if (use_varT) T25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k) if (use_covarTS) TS5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 80d94ec7fe..27a2217413 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -8,6 +8,7 @@ module MOM_isopycnal_slopes 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_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 @@ -26,7 +27,7 @@ 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, & +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) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure @@ -38,6 +39,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !! thermodynamic variables real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity !! times a smoothing timescale [Z2 ~> m2]. + logical, intent(in) :: use_stanley !< turn on stanley param in slope real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & @@ -70,12 +72,15 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & ! 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. real, dimension(SZIB_(G)) :: & drho_dT_u, & ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1]. drho_dS_u ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1]. real, dimension(SZI_(G)) :: & drho_dT_v, & ! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1]. - drho_dS_v ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1]. + drho_dS_v, & ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1]. + drho_dT_dT_h, & ! The second derivative of density with temperature at h points [R degC-2 ~> kg m-3 degC-2] + drho_dT_dT_hr ! The second derivative of density with temperature at h (+1) points [R degC-2 ~> kg m-3 degC-2] real, dimension(SZIB_(G)) :: & T_u, & ! Temperature on the interface at the u-point [degC]. S_u, & ! Salinity on the interface at the u-point [ppt]. @@ -84,6 +89,13 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & T_v, & ! Temperature on the interface at the v-point [degC]. S_v, & ! Salinity on the interface at the v-point [ppt]. pres_v ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa]. + real, dimension(SZI_(G)) :: & + T_h, & ! Temperature on the interface at the h-point [degC]. + S_h, & ! Salinity on the interface at the h-point [ppt] + pres_h, & ! Pressure on the interface at the h-point [R L2 T-2 ~> Pa]. + T_hr, & ! Temperature on the interface at the h (+1) point [degC]. + S_hr, & ! Salinity on the interface at the h (+1) point [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]. @@ -215,9 +227,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !$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) & + !$OMP 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, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & !$OMP drdx,mag_grad2,slope,slope2_Ratio,l_seg) do j=js,je ; do K=nz,2,-1 @@ -237,6 +250,19 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & tv%eqn_of_state, EOSdom_u) endif + if (use_stanley) then + do i=is-1,ie+1 + pres_h(i) = pres(i,j,K) + T_h(i) = 0.5*(T(i,j,k) + T(i,j,k-1)) + S_h(i) = 0.5*(S(i,j,k) + S(i,j,k-1)) + enddo + ! The second line below would correspond to arguments + ! 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, & + is-1, ie-is+3, tv%eqn_of_state) + endif + do I=is-1,ie if (use_EOS) then ! Estimate the horizontal density gradients along layers. @@ -251,7 +277,14 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & drdkR = (drho_dT_u(I) * (T(i+1,j,k)-T(i+1,j,k-1)) + & drho_dS_u(I) * (S(i+1,j,k)-S(i+1,j,k-1))) endif - + if (use_stanley) then + ! Correction to the horizontal density gradient due to nonlinearity in + ! the EOS rectifying SGS temperature anomalies + drdiA = drdiA + 0.5 * ((drho_dT_dT_h(i+1) * tv%varT(i+1,j,k-1)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k-1)) ) + drdiB = drdiB + 0.5 * ((drho_dT_dT_h(i+1) * tv%varT(i+1,j,k)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k)) ) + endif hg2A = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2 hg2B = h(i,j,k)*h(i+1,j,k) + h_neglect2 @@ -325,9 +358,11 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !$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_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, & - !$OMP dzv,local_open_v_BC,OBC) & + !$OMP dzv,local_open_v_BC,OBC,use_stanley) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & + !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & + !$OMP drho_dT_dT_hr,pres_hr,T_hr,S_hr, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & !$OMP drdy,mag_grad2,slope,slope2_Ratio,l_seg) do j=js-1,je ; do K=nz,2,-1 @@ -345,6 +380,25 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, tv%eqn_of_state, & EOSdom_v) endif + if (use_stanley) then + do i=is,ie + pres_h(i) = pres(i,j,K) + T_h(i) = 0.5*(T(i,j,k) + T(i,j,k-1)) + S_h(i) = 0.5*(S(i,j,k) + S(i,j,k-1)) + + pres_hr(i) = pres(i,j+1,K) + T_hr(i) = 0.5*(T(i,j+1,k) + T(i,j+1,k-1)) + S_hr(i) = 0.5*(S(i,j+1,k) + S(i,j+1,k-1)) + enddo + ! The second line below would correspond to arguments + ! 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, & + is, ie-is+1, tv%eqn_of_state) + call calculate_density_second_derivs(T_hr, S_hr, pres_hr, & + scrap, scrap, drho_dT_dT_hr, scrap, scrap, & + is, ie-is+1, tv%eqn_of_state) + endif do i=is,ie if (use_EOS) then ! Estimate the horizontal density gradients along layers. @@ -359,6 +413,14 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & drdkR = (drho_dT_v(i) * (T(i,j+1,k)-T(i,j+1,k-1)) + & drho_dS_v(i) * (S(i,j+1,k)-S(i,j+1,k-1))) endif + if (use_stanley) then + ! Correction to the horizontal density gradient due to nonlinearity in + ! the EOS rectifying SGS temperature anomalies + drdjA = drdjA + 0.5 * ((drho_dT_dT_hr(i) * tv%varT(i,j+1,k-1)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k-1)) ) + drdjB = drdjB + 0.5 * ((drho_dT_dT_hr(i) * tv%varT(i,j+1,k)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k)) ) + endif hg2A = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2 hg2B = h(i,j,k)*h(i,j+1,k) + h_neglect2 diff --git a/src/core/MOM_stoch_eos.F90 b/src/core/MOM_stoch_eos.F90 new file mode 100644 index 0000000000..0ee6d6b1be --- /dev/null +++ b/src/core/MOM_stoch_eos.F90 @@ -0,0 +1,212 @@ +!> Provides the ocean stochastic equation of state +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 + +implicit none +#include + +public MOM_stoch_eos_init +public MOM_stoch_eos_run +public MOM_calc_varT + +real,private ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: l2_inv + !< One over sum of the T cell side side lengths squared +real,private ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: rgauss !< nondimensional random Gaussian +real, parameter,private :: tfac=0.27 !< Nondimensional decorrelation time factor, ~1/3.7 +real, parameter,private :: amplitude=0.624499 !< Nondimensional std dev of Gaussian +integer ,private :: seed !< PRNG seed +type(PRNG) :: rn_CS !< PRNG control structure + +!> Describes parameters of the stochastic component of the EOS +!! correction, described in Stanley et al. JAMES 2020. +type, public :: MOM_stoch_eos_CS + real,public ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: pattern + !< Random pattern for stochastic EOS + real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: phi + !< temporal correlation stochastic EOS (deugging) + 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 ! m2 s-1] + !>@{ Diagnostic IDs + integer :: id_stoch_eos = -1, id_stoch_phi = -1, id_tvar_sgs = -1 + !>@} + +end type MOM_stoch_eos_CS + + +contains + subroutine MOM_stoch_eos_init(G,Time,param_file,stoch_eos_CS,restart_CS,diag) +! initialization subroutine called by MOM.F90, + 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) :: stoch_eos_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 + integer :: i,j + type(vardesc) :: vd + seed=0 + ! contants + !pi=2*acos(0.0) + call get_param(param_file, "MOM_stoch_eos", "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", stoch_eos_CS%stanley_coeff, & + "Coefficient correlating the temperature gradient "//& + "and SGS T variance.", default=-1.0) + call get_param(param_file, "MOM_stoch_eos", "STANLEY_A", stoch_eos_CS%stanley_a, & + "Coefficient a which scales chi in stochastic perturbation of the "//& + "SGS T variance.", default=1.0) + call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", stoch_eos_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) + + !don't run anything if STANLEY_COEFF < 0 + if (stoch_eos_CS%stanley_coeff >= 0.0) then + + ALLOC_(stoch_eos_CS%pattern(G%isd:G%ied,G%jsd:G%jed)) ; stoch_eos_CS%pattern(:,:) = 0.0 + vd = var_desc("stoch_eos_pattern","nondim","Random pattern for stoch EOS",'h','1') + call register_restart_field(stoch_eos_CS%pattern, vd, .false., restart_CS) + ALLOC_(stoch_eos_CS%phi(G%isd:G%ied,G%jsd:G%jed)) ; stoch_eos_CS%phi(:,:) = 0.0 + ALLOC_(l2_inv(G%isd:G%ied,G%jsd:G%jed)) + ALLOC_(rgauss(G%isd:G%ied,G%jsd:G%jed)) + call get_param(param_file, "MOM_stoch_eos", "SEED_STOCH_EOS", seed, & + "Specfied seed for random number sequence ", default=0) + call random_2d_constructor(rn_CS, G%HI, Time, seed) + call random_2d_norm(rn_CS, G%HI, rgauss) + ! 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 + 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 + stoch_eos_CS%pattern(i,j)=amplitude*rgauss(i,j) + enddo + enddo + endif + + !register diagnostics + stoch_eos_CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs', diag%axesTL, Time, & + 'Parameterized SGS Temperature Variance ', 'None') + if (stoch_eos_CS%use_stoch_eos) then + stoch_eos_CS%id_stoch_eos = register_diag_field('ocean_model', 'stoch_eos', diag%axesT1, Time, & + 'random pattern for EOS', 'None') + stoch_eos_CS%id_stoch_phi = register_diag_field('ocean_model', 'stoch_phi', diag%axesT1, Time, & + 'phi for EOS', 'None') + endif + endif + + end subroutine MOM_stoch_eos_init + + subroutine MOM_stoch_eos_run(G,u,v,delt,Time,stoch_eos_CS,diag) + 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]. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. + 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) :: stoch_eos_CS !< Stochastic control structure + type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics +! locals + integer :: i,j + integer :: yr,mo,dy,hr,mn,sc + real :: phi,ubar,vbar + + call random_2d_constructor(rn_CS, G%HI, Time, seed) + call random_2d_norm(rn_CS, G%HI, rgauss) + ! 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*tfac*sqrt((ubar**2+vbar**2)*l2_inv(i,j))) + stoch_eos_CS%pattern(i,j)=phi*stoch_eos_CS%pattern(i,j) + amplitude*sqrt(1-phi**2)*rgauss(i,j) + stoch_eos_CS%phi(i,j)=phi + enddo + enddo + + end subroutine MOM_stoch_eos_run + + + subroutine MOM_calc_varT(G,GV,h,tv,stoch_eos_CS,dt) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m] + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + type(MOM_stoch_eos_CS), intent(inout) :: stoch_eos_CS !< Stochastic control structure. + real, intent(in) :: dt !< Time increment [T ~> s] +! locals + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & + T, & !> The temperature (or density) [degC], with the values in + !! in massless layers filled vertically by diffusion. + S !> The filled salinity [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 [degC2] + + ! 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) + + call vert_fill_TS(h, tv%T, tv%S, stoch_eos_CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.) + + do k=1,G%ke + do j=G%jsc,G%jec + do i=G%isc,G%iec + hl(1) = h(i,j,k) * G%mask2dT(i,j) + hl(2) = h(i-1,j,k) * G%mask2dCu(I-1,j) + hl(3) = h(i+1,j,k) * G%mask2dCu(I,j) + hl(4) = h(i,j-1,k) * G%mask2dCv(i,J-1) + hl(5) = h(i,j+1,k) * G%mask2dCv(i,J) + + ! SGS variance in i-direction [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 + ! SGS variance in j-direction [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 + tv%varT(i,j,k) = stoch_eos_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)) + enddo + enddo + enddo + ! if stochastic, perturb + if (stoch_eos_CS%use_stoch_eos) then + do k=1,G%ke + do j=G%jsc,G%jec + do i=G%isc,G%iec + tv%varT(i,j,k) = exp (stoch_eos_CS%stanley_a * stoch_eos_CS%pattern(i,j)) * tv%varT(i,j,k) + enddo + enddo + enddo + endif + end subroutine MOM_calc_varT + +end module MOM_stoch_eos + diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 8d667503d7..c98fe7539a 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -106,6 +106,8 @@ module MOM_diagnostics integer :: id_rhopot0 = -1, id_rhopot2 = -1 integer :: id_drho_dT = -1, id_drho_dS = -1 integer :: id_h_pre_sync = -1 + integer :: id_tosq = -1, id_sosq = -1 + !>@} type(wave_speed_CS) :: wave_speed !< Wave speed control struct @@ -402,16 +404,28 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! Internal T&S variables are conservative temperature & absolute salinity, ! so they need to converted to potential temperature and practical salinity ! for some diagnostics using TEOS-10 function calls. - if ((CS%id_Tpot > 0) .or. (CS%id_tob > 0)) then + if ((CS%id_Tpot > 0) .or. (CS%id_tob > 0) .or. (CS%id_tosq > 0)) then do k=1,nz ; do j=js,je ; do i=is,ie work_3d(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k)) enddo ; enddo ; enddo if (CS%id_Tpot > 0) call post_data(CS%id_Tpot, work_3d, CS%diag) if (CS%id_tob > 0) call post_data(CS%id_tob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_tosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = work_3d(i,j,k)*work_3d(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_tosq, work_3d, CS%diag) + endif endif else ! Internal T&S variables are potential temperature & practical salinity if (CS%id_tob > 0) call post_data(CS%id_tob, tv%T(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_tosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = tv%T(i,j,k)*tv%T(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_tosq, work_3d, CS%diag) + endif endif ! Calculate additional, potentially derived salinity diagnostics @@ -419,16 +433,28 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! Internal T&S variables are conservative temperature & absolute salinity, ! so they need to converted to potential temperature and practical salinity ! for some diagnostics using TEOS-10 function calls. - if ((CS%id_Sprac > 0) .or. (CS%id_sob > 0)) then + if ((CS%id_Sprac > 0) .or. (CS%id_sob > 0) .or. (CS%id_sosq >0)) then do k=1,nz ; do j=js,je ; do i=is,ie work_3d(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k)) enddo ; enddo ; enddo if (CS%id_Sprac > 0) call post_data(CS%id_Sprac, work_3d, CS%diag) if (CS%id_sob > 0) call post_data(CS%id_sob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_sosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = work_3d(i,j,k)*work_3d(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_sosq, work_3d, CS%diag) + endif endif else ! Internal T&S variables are potential temperature & practical salinity if (CS%id_sob > 0) call post_data(CS%id_sob, tv%S(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_sosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = tv%S(i,j,k)*tv%S(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_sosq, work_3d, CS%diag) + endif endif ! volume mean potential temperature @@ -1619,6 +1645,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag long_name='Sea Water Salinity at Sea Floor', & standard_name='sea_water_salinity_at_sea_floor', units='psu') + CS%id_tosq = register_diag_field('ocean_model', 'tosq', diag%axesTL,& + Time, 'Square of Potential Temperature', 'degc2', & + standard_name='Potential Temperature Squared') + CS%id_sosq = register_diag_field('ocean_model', 'sosq', diag%axesTL,& + Time, 'Square of Salinity', 'psu2', & + standard_name='Salinity Squared') + CS%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', & diag%axesZL, Time, 'Layer Average Ocean Temperature', 'degC') CS%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', & diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 39b626985a..881cea329d 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -64,9 +64,9 @@ module MOM_EOS !> Calculates density of sea water from T, S and P interface calculate_density - module procedure calculate_density_scalar, calculate_density_array, calculate_density_1d - module procedure calculate_stanley_density_scalar, calculate_stanley_density_array - module procedure calculate_stanley_density_1d + module procedure calculate_density_scalar, calculate_density_array, calculate_density_1d, & + calculate_stanley_density_scalar, calculate_stanley_density_array, & + calculate_stanley_density_1d end interface calculate_density !> Calculates specific volume of sea water from T, S and P @@ -430,18 +430,18 @@ subroutine calculate_stanley_density_1d(T, S, pressure, Tvar, TScov, Svar, rho, select case (EOS%form_of_EOS) case (EOS_LINEAR) - call calculate_density_linear(T, S, pres, rho, 1, npts, & + call calculate_density_linear(T, S, pres, rho, is, npts, & EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS, rho_ref) call calculate_density_second_derivs_linear(T, S, pres, d2RdSS, d2RdST, & - d2RdTT, d2RdSp, d2RdTP, 1, npts) + d2RdTT, d2RdSp, d2RdTP, is, npts) case (EOS_WRIGHT) - call calculate_density_wright(T, S, pres, rho, 1, npts, rho_ref) + call calculate_density_wright(T, S, pres, rho, is, npts, rho_ref) call calculate_density_second_derivs_wright(T, S, pres, d2RdSS, d2RdST, & - d2RdTT, d2RdSp, d2RdTP, 1, npts) + d2RdTT, d2RdSp, d2RdTP, is, npts) case (EOS_TEOS10) - call calculate_density_teos10(T, S, pres, rho, 1, npts, rho_ref) + call calculate_density_teos10(T, S, pres, rho, is, npts, rho_ref) call calculate_density_second_derivs_teos10(T, S, pres, d2RdSS, d2RdST, & - d2RdTT, d2RdSp, d2RdTP, 1, npts) + d2RdTT, d2RdSp, d2RdTP, is, npts) case default call MOM_error(FATAL, "calculate_stanley_density_scalar: EOS is not valid.") end select diff --git a/src/equation_of_state/TEOS10/gsw_mod_error_functions.f90 b/src/equation_of_state/TEOS10/gsw_mod_error_functions.f90 new file mode 120000 index 0000000000..1c3b7bfb3c --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_error_functions.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/modules/gsw_mod_error_functions.f90 \ No newline at end of file diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index bef78a433a..38e330c4be 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -156,6 +156,7 @@ subroutine random_2d_constructor(CS, HI, Time, seed) if (.not. allocated(CS%stream2d)) allocate( CS%stream2d(HI%isd:HI%ied,HI%jsd:HI%jed) ) tseed = seed_from_time(Time) + tseed = ieor(tseed*9007, seed) do j = HI%jsd,HI%jed do i = HI%isd,HI%ied diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index da8e936642..236a2cebce 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -59,6 +59,7 @@ module MOM_lateral_mixing_coeffs !! This parameter is set depending on other parameters. logical :: calculate_Eady_growth_rate !< If true, calculate all the Eady growth rate. !! This parameter is set depending on other parameters. + logical :: use_stanley_iso !< If true, use Stanley parameterization in MOM_isopycnal_slopes logical :: use_simpler_Eady_growth_rate !< If true, use a simpler method to calculate the !! Eady growth rate that avoids division by layer thickness. !! This parameter is set depending on other parameters. @@ -471,14 +472,14 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) if (CS%calculate_Eady_growth_rate) then if (CS%use_simpler_Eady_growth_rate) then call find_eta(h, tv, G, GV, US, e, halo_size=2) - call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, & + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, dzu=dzu, dzv=dzv, & dzSxN=dzSxN, dzSyN=dzSyN, halo=1, OBC=OBC) call calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, dzSyN, CS%SN_u, CS%SN_v) else call find_eta(h, tv, G, GV, US, e, halo_size=2) if (CS%use_stored_slopes) then - call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, & + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC) call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC) else @@ -1267,6 +1268,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.) + 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%Resoln_use_ebt .or. CS%khth_use_ebt_struct) then in_use = .true. @@ -1290,6 +1294,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) endif if (CS%use_stored_slopes) then + ! CS%calculate_Eady_growth_rate=.true. in_use = .true. allocate(CS%slope_x(IsdB:IedB,jsd:jed,GV%ke+1), source=0.0) allocate(CS%slope_y(isd:ied,JsdB:JedB,GV%ke+1), source=0.0) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 04982d7171..7ebc6c0eff 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -60,6 +60,7 @@ module MOM_mixed_layer_restrat logical :: debug = .false. !< If true, calculate checksums of fields for debugging. type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. + logical :: use_stanley_ml !< If true, use the Stanley parameterization of SGS T variance real, dimension(:,:), allocatable :: & MLD_filtered, & !< Time-filtered MLD [H ~> m or kg m-2] @@ -178,6 +179,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer ! densities [R L2 T-2 ~> Pa]. real, dimension(SZI_(G)) :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2 + real, dimension(SZI_(G)) :: covTS, varS !SGS TS covariance, S variance in Stanley param; currently 0 real :: aFac, bFac ! Nondimensional ratios [nondim] real :: ddRho ! A density difference [R ~> kg m-3] real :: hAtVel, zpa, zpb, dh, res_scaling_fac @@ -198,6 +200,9 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + covTS(:)=0.0 !!Functionality not implemented yet; in future, should be passed in tv + varS(:)=0.0 + if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// & "An equation of state must be used with this module.") if (.not. allocated(VarMix%Rd_dx_h) .and. CS%front_length > 0.) & @@ -210,7 +215,12 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var EOSdom(:) = EOS_domain(G%HI, halo=1) do j = js-1, je+1 dK(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom) + if (CS%use_stanley_ml) then + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, tv%varT(:,j,1), covTS, varS, & + rhoSurf, tv%eqn_of_state, EOSdom) + else + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom) + endif deltaRhoAtK(:) = 0. MLD_fast(:,j) = 0. do k = 2, nz @@ -218,7 +228,12 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Depth of center of layer K ! Mixed-layer depth, using sigma-0 (surface reference pressure) deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, EOSdom) + if (CS%use_stanley_ml) then + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, tv%varT(:,j,k), covTS, varS, & + deltaRhoAtK, tv%eqn_of_state, EOSdom) + else + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, EOSdom) + endif do i = is-1,ie+1 deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface enddo @@ -303,6 +318,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var !$OMP h_neglect,g_Rho0,I4dt,CS,uhml,uhtr,dt,vhml,vhtr,EOSdom, & !$OMP utimescale_diag,vtimescale_diag,forces,dz_neglect, & !$OMP htot_slow,MLD_slow,Rml_av_slow,VarMix,I_LFront, & +!$OMP covTS, varS, & !$OMP res_upscale, nz,MLD_fast,uDml_diag,vDml_diag) & !$OMP private(rho_ml,h_vel,u_star,absf,mom_mixrate,timescale, & !$OMP line_is_empty, keep_going,res_scaling_fac, & @@ -320,7 +336,12 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var h_avail(i,j,k) = max(I4dt*G%areaT(i,j)*(h(i,j,k)-GV%Angstrom_H),0.0) enddo if (keep_going) then - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, EOSdom) + if (CS%use_stanley_ml) then + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, tv%varT(:,j,k), covTS, varS, & + rho_ml(:), tv%eqn_of_state, EOSdom) + else + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, EOSdom) + endif line_is_empty = .true. do i=is-1,ie+1 if (htot_fast(i,j) < MLD_fast(i,j)) then @@ -619,6 +640,10 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) if ((nkml<2) .or. (CS%ml_restrat_coef<=0.0)) return + if (CS%use_stanley_ml) call MOM_error(FATAL, & + "MOM_mixedlayer_restrat: The Stanley parameterization is not"//& + "available with the BML.") + uDml(:) = 0.0 ; vDml(:) = 0.0 I4dt = 0.25 / dt g_Rho0 = GV%g_Earth / GV%Rho0 @@ -842,6 +867,9 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, "geostrophic kinetic energy or 1 plus the square of the "//& "grid spacing over the deformation radius, as detailed "//& "by Fox-Kemper et al. (2010)", units="nondim", default=0.0) + 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.) ! We use GV%nkml to distinguish between the old and new implementation of MLE. ! The old implementation only works for the layer model with nkml>0. if (GV%nkml==0) then diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 8303d30621..14a7c30e3e 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -78,9 +78,7 @@ module MOM_thickness_diffuse !! than the streamfunction for the GM source term. logical :: use_GM_work_bug !< If true, use the incorrect sign for the !! top-level work tendency on the top layer. - real :: Stanley_det_coeff !< The coefficient correlating SGS temperature variance with the mean - !! temperature gradient in the deterministic part of the Stanley parameterization. - !! Negative values disable the scheme." [nondim] + logical :: use_stanley_gm !< If true, also use the Stanley parameterization in MOM_thickness_diffuse type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics real, allocatable :: GMwork(:,:) !< Work by thickness diffusivity [R Z L2 T-3 ~> W m-2] @@ -611,13 +609,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV h_avail_rsum ! The running sum of h_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZIB_(G)) :: & drho_dT_u, & ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1] - drho_dS_u, & ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1]. - drho_dT_dT_u ! The second derivative of density with temperature at u points [R degC-2 ~> kg m-3 degC-2] - real, dimension(SZIB_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ingored. + drho_dS_u ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1]. + real, dimension(SZI_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ingored. real, dimension(SZI_(G)) :: & drho_dT_v, & ! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1] drho_dS_v, & ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1]. - drho_dT_dT_v ! The second derivative of density with temperature at v points [R degC-2 ~> kg m-3 degC-2] + drho_dT_dT_h, & ! The second derivative of density with temperature at h points [R degC-2 ~> kg m-3 degC-2] + drho_dT_dT_hr ! The second derivative of density with temperature at h (+1) points [R degC-2 ~> kg m-3 degC-2] real :: uhtot(SZIB_(G), SZJ_(G)) ! The vertical sum of uhD [H L2 T-1 ~> m3 s-1 or kg s-1]. real :: vhtot(SZI_(G), SZJB_(G)) ! The vertical sum of vhD [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZIB_(G)) :: & @@ -628,6 +626,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV T_v, & ! Temperature on the interface at the v-point [degC]. S_v, & ! Salinity on the interface at the v-point [ppt]. pres_v ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa]. + real, dimension(SZI_(G)) :: & + T_h, & ! Temperature on the interface at the h-point [degC]. + S_h, & ! Salinity on the interface at the h-point [ppt]. + pres_h, & ! Pressure on the interface at the h-point [R L2 T-2 ~> Pa]. + T_hr, & ! Temperature on the interface at the h (+1) point [degC]. + S_hr, & ! Salinity on the interface at the h (+1) point [ppt]. + pres_hr ! Pressure on the interface at the h (+1) point [R L2 T-2 ~> Pa]. real :: Work_u(SZIB_(G), SZJ_(G)) ! The work being done by the thickness real :: Work_v(SZI_(G), SZJB_(G)) ! diffusion integrated over a cell [R Z L4 T-3 ~> W ] real :: Work_h ! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2]. @@ -683,13 +688,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real :: G_rho0 ! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2]. real :: N2_floor ! A floor for N2 to avoid degeneracy in the elliptic solver ! times unit conversion factors [T-2 L2 Z-2 ~> s-2] - real :: Tl(5) ! copy and T in local stencil [degC] - real :: mn_T ! mean of T in local stencil [degC] - real :: mn_T2 ! mean of T**2 in local stencil [degC] - real :: hl(5) ! Copy of local stencil of H [H ~> m] - real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1] - real, dimension(SZI_(G), SZJ_(G),SZK_(GV)) :: Tsgs2 ! Sub-grid temperature variance [degC2] - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: diag_sfn_x, diag_sfn_unlim_x ! Diagnostics real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics logical :: present_slope_x, present_slope_y, calc_derivatives @@ -697,7 +695,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! 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. - logical :: use_Stanley + logical :: use_stanley integer :: is, ie, js, je, nz, IsdB, halo integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke ; IsdB = G%IsdB @@ -714,7 +712,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV use_EOS = associated(tv%eqn_of_state) present_slope_x = PRESENT(slope_x) present_slope_y = PRESENT(slope_y) - use_Stanley = CS%Stanley_det_coeff >= 0. + + use_stanley = CS%use_stanley_gm nk_linear = max(GV%nkml, 1) @@ -728,17 +727,15 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (use_EOS) then halo = 1 ! Default halo to fill is 1 - if (use_Stanley) halo = 2 ! Need wider valid halo for gradients of T call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo, larger_h_denom=.true.) endif if (CS%use_FGNV_streamfn .and. .not. associated(cg1)) call MOM_error(FATAL, & "cg1 must be associated when using FGNV streamfunction.") -!$OMP parallel default(none) shared(is,ie,js,je,h_avail_rsum,pres,h_avail,I4dt, use_Stanley, & -!$OMP CS,G,GV,tv,h,h_frac,nz,uhtot,Work_u,vhtot,Work_v,Tsgs2,T, & -!$OMP diag_sfn_x, diag_sfn_y, diag_sfn_unlim_x, diag_sfn_unlim_y ) & -!$OMP private(hl,r_sm_H,Tl,mn_T,mn_T2) +!$OMP parallel default(none) shared(is,ie,js,je,h_avail_rsum,pres,h_avail,I4dt, use_stanley, & +!$OMP CS,G,GV,tv,h,h_frac,nz,uhtot,Work_u,vhtot,Work_v,T, & +!$OMP diag_sfn_x, diag_sfn_y, diag_sfn_unlim_x, diag_sfn_unlim_y ) ! Find the maximum and minimum permitted streamfunction. !$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 @@ -751,41 +748,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV h_frac(i,j,1) = 1.0 pres(i,j,2) = pres(i,j,1) + (GV%g_Earth*GV%H_to_RZ) * h(i,j,1) enddo ; enddo - if (use_Stanley) then -!$OMP do - do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 - !! SGS variance in i-direction [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 - !! SGS variance in j-direction [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 - !Tsgs2(i,j,k) = CS%Stanley_det_coeff * 0.5 * ( dTdi2 + dTdj2 ) - ! 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. - hl(1) = h(i,j,k) * G%mask2dT(i,j) - hl(2) = h(i-1,j,k) * G%mask2dCu(I-1,j) - hl(3) = h(i+1,j,k) * G%mask2dCu(I,j) - hl(4) = h(i,j-1,k) * G%mask2dCv(i,J-1) - hl(5) = h(i,j+1,k) * G%mask2dCv(i,J) - r_sm_H = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + GV%H_subroundoff ) - ! Mean of T - Tl(1) = T(i,j,k) ; Tl(2) = T(i-1,j,k) ; Tl(3) = T(i+1,j,k) - Tl(4) = T(i,j-1,k) ; Tl(5) = T(i,j+1,k) - mn_T = ( hl(1)*Tl(1) + ( ( hl(2)*Tl(2) + hl(3)*Tl(3) ) + ( hl(4)*Tl(4) + hl(5)*Tl(5) ) ) ) * r_sm_H - ! Adjust T vectors to have zero mean - Tl(:) = Tl(:) - mn_T ; mn_T = 0. - ! Variance of T - mn_T2 = ( hl(1)*Tl(1)*Tl(1) + ( ( hl(2)*Tl(2)*Tl(2) + hl(3)*Tl(3)*Tl(3) ) & - + ( hl(4)*Tl(4)*Tl(4) + hl(5)*Tl(5)*Tl(5) ) ) ) * r_sm_H - ! Variance should be positive but round-off can violate this. Calculating - ! variance directly would fix this but requires more operations. - Tsgs2(i,j,k) = CS%Stanley_det_coeff * max(0., mn_T2) - enddo ; enddo ; enddo - endif !$OMP do do j=js-1,je+1 do k=2,nz ; do i=is-1,ie+1 @@ -815,11 +777,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !$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,use_stanley, & !$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, & -!$OMP drho_dT_dT_u,scrap, & +!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & !$OMP drdx,mag_grad2,Slope,slope2_Ratio_u,hN2_u, & !$OMP Sfn_unlim_u,drdi_u,drdkDe_u,h_harm,c2_h_u, & @@ -833,7 +795,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV endif calc_derivatives = use_EOS .and. (k >= nk_linear) .and. & - (find_work .or. .not. present_slope_x .or. CS%use_FGNV_streamfn .or. use_Stanley) + (find_work .or. .not. present_slope_x .or. CS%use_FGNV_streamfn .or. use_stanley) ! Calculate the zonal fluxes and gradients. if (calc_derivatives) then @@ -845,12 +807,17 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, & tv%eqn_of_state, EOSdom_u) endif - if (use_Stanley) then + if (use_stanley) then + do i=is-1,ie+1 + pres_h(i) = pres(i,j,K) + T_h(i) = 0.5*(T(i,j,k) + T(i,j,k-1)) + S_h(i) = 0.5*(S(i,j,k) + S(i,j,k-1)) + enddo ! The second line below would correspond to arguments ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & - call calculate_density_second_derivs(T_u, S_u, pres_u, & - scrap, scrap, drho_dT_dT_u, scrap, scrap, & - (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state) + call calculate_density_second_derivs(T_h, S_h, pres_h, & + scrap, scrap, drho_dT_dT_h, scrap, scrap, & + is-1, ie-is+3, tv%eqn_of_state) endif do I=is-1,ie @@ -870,11 +837,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV elseif (find_work) then ! This is used in pure stacked SW mode drdkDe_u(I,K) = drdkR * e(i+1,j,K) - drdkL * e(i,j,K) endif - if (use_Stanley) then + if (use_stanley) then ! Correction to the horizontal density gradient due to nonlinearity in ! the EOS rectifying SGS temperature anomalies - drdiA = drdiA + drho_dT_dT_u(I) * 0.5 * ( Tsgs2(i+1,j,k-1)-Tsgs2(i,j,k-1) ) - drdiB = drdiB + drho_dT_dT_u(I) * 0.5 * ( Tsgs2(i+1,j,k)-Tsgs2(i,j,k) ) + drdiA = drdiA + 0.5 * ((drho_dT_dT_h(i+1) * tv%varT(i+1,j,k-1)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k-1)) ) + drdiB = drdiB + 0.5 * ((drho_dT_dT_h(i+1) * tv%varT(i+1,j,k)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k)) ) endif if (find_work) drdi_u(I,k) = drdiB @@ -1083,11 +1052,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & !$OMP h_neglect2,int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, & !$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1,diag_sfn_y, & -!$OMP diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley,Tsgs2, & +!$OMP diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley, & !$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & -!$OMP drho_dT_dT_v,scrap, & +!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & +!$OMP drho_dT_dT_hr,pres_hr,T_hr,S_hr, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & !$OMP drdy,mag_grad2,Slope,slope2_Ratio_v,hN2_v, & !$OMP Sfn_unlim_v,drdj_v,drdkDe_v,h_harm,c2_h_v, & @@ -1100,7 +1070,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV endif calc_derivatives = use_EOS .and. (k >= nk_linear) .and. & - (find_work .or. .not. present_slope_y .or. CS%use_FGNV_streamfn .or. use_Stanley) + (find_work .or. .not. present_slope_y .or. CS%use_FGNV_streamfn .or. use_stanley) if (calc_derivatives) then do i=is,ie @@ -1111,11 +1081,23 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, & tv%eqn_of_state, EOSdom_v) endif - if (use_Stanley) then + if (use_stanley) then + do i=is,ie + pres_h(i) = pres(i,j,K) + T_h(i) = 0.5*(T(i,j,k) + T(i,j,k-1)) + S_h(i) = 0.5*(S(i,j,k) + S(i,j,k-1)) + + pres_hr(i) = pres(i,j+1,K) + T_hr(i) = 0.5*(T(i,j+1,k) + T(i,j+1,k-1)) + S_hr(i) = 0.5*(S(i,j+1,k) + S(i,j+1,k-1)) + enddo ! The second line below would correspond to arguments ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & - call calculate_density_second_derivs(T_v, S_v, pres_v, & - scrap, scrap, drho_dT_dT_v, scrap, scrap, & + call calculate_density_second_derivs(T_h, S_h, pres_h, & + scrap, scrap, drho_dT_dT_h, scrap, scrap, & + is, ie-is+1, tv%eqn_of_state) + call calculate_density_second_derivs(T_hr, S_hr, pres_hr, & + scrap, scrap, drho_dT_dT_hr, scrap, scrap, & is, ie-is+1, tv%eqn_of_state) endif do i=is,ie @@ -1135,11 +1117,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV elseif (find_work) then ! This is used in pure stacked SW mode drdkDe_v(i,K) = drdkR * e(i,j+1,K) - drdkL * e(i,j,K) endif - if (use_Stanley) then + if (use_stanley) then ! Correction to the horizontal density gradient due to nonlinearity in ! the EOS rectifying SGS temperature anomalies - drdjA = drdjA + drho_dT_dT_v(I) * 0.5 * ( Tsgs2(i,j+1,k-1)-Tsgs2(i,j,k-1) ) - drdjB = drdjB + drho_dT_dT_v(I) * 0.5 * ( Tsgs2(i,j+1,k)-Tsgs2(i,j,k) ) + drdjA = drdjA + 0.5 * ((drho_dT_dT_hr(i) * tv%varT(i,j+1,k-1)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k-1)) ) + drdjB = drdjB + 0.5 * ((drho_dT_dT_hr(i) * tv%varT(i,j+1,k)) - & + (drho_dT_dT_h(i) * tv%varT(i,j,k)) ) endif if (find_work) drdj_v(i,k) = drdjB @@ -1969,10 +1953,9 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "streamfunction formulation, expressed as a fraction of planetary "//& "rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", & default=1.e-15, units="nondim", do_not_log=.not.CS%use_FGNV_streamfn) - call get_param(param_file, mdl, "STANLEY_PRM_DET_COEFF", CS%Stanley_det_coeff, & - "The coefficient correlating SGS temperature variance with the mean "//& - "temperature gradient in the deterministic part of the Stanley parameterization. "//& - "Negative values disable the scheme.", units="nondim", default=-1.0) + 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.) 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)