diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index d1fd8619da..9d0701ec2a 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. @@ -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. @@ -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) @@ -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) @@ -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. @@ -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) @@ -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) @@ -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()") @@ -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.) 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, & 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/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 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) 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)