Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

+(*)Fix numerous issues with MOM_stoch_eos and Stanley parameterizations #271

Merged
merged 6 commits into from
Dec 19, 2022
Merged

+(*)Fix numerous issues with MOM_stoch_eos and Stanley parameterizations #271

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
125e3a3
+(*)Fix numerous issues with MOM_stoch_eos
Hallberg-NOAA Dec 8, 2022
eeb07b3
(*)Correct USE_STANLEY domain extents in EOS calls
Hallberg-NOAA Dec 8, 2022
d4cb1d4
(*)Correct scaling of USE_STANLEY_PGF diagnostics
Hallberg-NOAA Dec 8, 2022
e18e2da
(*)Test for inconsistent USE_STANLEY flags
Hallberg-NOAA Dec 8, 2022
5b2c584
(*)Correct scaling of DT_OBC_SEG_UPDATE_OBGC
Hallberg-NOAA Dec 8, 2022
b8c3845
Rescale interface heights for MOM_IC
Hallberg-NOAA Dec 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 25 additions & 15 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -298,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.

Expand Down Expand Up @@ -1079,12 +1084,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)

Expand Down Expand Up @@ -1297,9 +1302,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)

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -2679,6 +2681,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)

Expand Down Expand Up @@ -2964,7 +2970,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)
Expand Down Expand Up @@ -3207,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()")
Expand All @@ -3230,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.)
Expand Down
86 changes: 59 additions & 27 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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).
Expand All @@ -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

Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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, &
Expand Down
Loading