From 6ad1530c3770fb400e980aff6b724e1c5639a1cf Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Tue, 2 Jul 2024 14:00:14 -0600 Subject: [PATCH] Stochastic GM+E (#280) * SKEB+GM This commit adds a Stochastic Kinetic Energy Backscatter Scheme (SKEBS) where the amplitude of the backscatter is based on either the GM work rate and/or the viscous work rate (FrictWork). Each of these can be multiplied by a coefficient so that, e.g. the backscatter rate could be 50% of the GM work rate plus 20% of the viscous work rate. The vertical structure of the backscatter rate associated with GM has vertical structure given by the EBT profile. This code was developed starting from Phil Pegion's branch, and it builds on the stochastic physics package (external) developed by Phil Pegion, Niraj Agarwal, and collaborators. This package allows the length and time scales of the backscatter to be set via namelist parameters. This commit may break the stochastic EPBL and SPPT schemes developed by P. Pegion. * Fix merge * Whitespace --- .../stochastic_physics/stochastic_physics.F90 | 11 +- src/core/MOM.F90 | 22 +- src/core/MOM_dynamics_split_RK2.F90 | 6 +- src/core/MOM_dynamics_unsplit.F90 | 8 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 6 +- .../lateral/MOM_hor_visc.F90 | 18 +- .../lateral/MOM_thickness_diffuse.F90 | 59 ++- .../stochastic/MOM_stochastics.F90 | 390 ++++++++++++++++-- 8 files changed, 453 insertions(+), 67 deletions(-) diff --git a/config_src/external/stochastic_physics/stochastic_physics.F90 b/config_src/external/stochastic_physics/stochastic_physics.F90 index 14fa1bf289..fdfd701892 100644 --- a/config_src/external/stochastic_physics/stochastic_physics.F90 +++ b/config_src/external/stochastic_physics/stochastic_physics.F90 @@ -16,7 +16,7 @@ module stochastic_physics !> Initializes the stochastic physics perturbations. subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_epbl_in, do_sppt_in, & - mpiroot, mpicomm, iret) + do_skeb_in,mpiroot, mpicomm, iret) real, intent(in) :: delt !< timestep in seconds between calls to run_stochastic_physics_ocn [s] integer, intent(in) :: nx !< number of gridpoints in the x-direction of the compute grid integer, intent(in) :: ny !< number of gridpoints in the y-direction of the compute grid @@ -25,6 +25,7 @@ subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_ real, intent(in) :: geoLatT(nx,ny) !< Latitude in degrees logical, intent(in) :: pert_epbl_in !< logical flag, if true generate random pattern for ePBL perturbations logical, intent(in) :: do_sppt_in !< logical flag, if true generate random pattern for SPPT perturbations + logical, intent(in) :: do_skeb_in !< logical flag, if true generate random pattern for SKEB perturbations integer, intent(in) :: mpiroot !< root processor integer, intent(in) :: mpicomm !< mpi communicator integer, intent(out) :: iret !< return code @@ -38,14 +39,20 @@ subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_ call MOM_error(WARNING, 'init_stochastic_physics_ocn: do_sppt needs to be false if using the stub') iret=-1 endif + if (do_skeb_in) then + call MOM_error(WARNING, 'init_stochastic_physics_ocn: do_skeb needs to be false if using the stub') + iret=-1 + endif ! This stub function does not actually do anything. return end subroutine init_stochastic_physics_ocn + !> Determines the stochastic physics perturbations. -subroutine run_stochastic_physics_ocn(sppt_wts, t_rp1, t_rp2) +subroutine run_stochastic_physics_ocn(sppt_wts, skeb_wts, t_rp1, t_rp2) real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2] + real, intent(inout) :: skeb_wts(:,:) !< array containing random weights for SKEB real, intent(inout) :: t_rp1(:,:) !< array containing random weights for ePBL !! perturbations (KE generation) range [0,2] real, intent(inout) :: t_rp2(:,:) !< array containing random weights for ePBL diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9098b245dd..7b080a5537 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -64,7 +64,7 @@ module MOM use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end use MOM_diabatic_driver, only : register_diabatic_restarts -use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS +use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS, apply_skeb use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init use MOM_diagnostics, only : register_transport_diags, post_transport_diagnostics use MOM_diagnostics, only : register_surface_diags, write_static_fields @@ -741,7 +741,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS endif endif ! advance the random pattern if stochastic physics is active - if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl) call update_stochastics(CS%stoch_CS) + if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl .OR. CS%stoch_CS%do_skeb) & + call update_stochastics(CS%stoch_CS) if (do_dyn) then if (nonblocking_p_surf_update) & @@ -1151,7 +1152,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%VarMix%use_variable_mixing) & call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, & - CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) + CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, & + CS%stoch_CS) call cpu_clock_end(id_clock_thick_diff) call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil)) if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)") @@ -1223,7 +1225,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, & - CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, waves=waves) + CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, CS%stoch_CS, waves=waves) if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)") elseif (CS%do_dynamics) then ! ------------------------------------ not SPLIT @@ -1237,11 +1239,13 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%use_RK2) then call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv, & + CS%stoch_CS) else call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, Waves=Waves) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, & + CS%stoch_CS, Waves=Waves) endif if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)") @@ -1282,7 +1286,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%VarMix%use_variable_mixing) & call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, & - CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) + CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, CS%stoch_CS) if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, scale=GV%H_to_MKS) call cpu_clock_end(id_clock_thick_diff) @@ -1584,6 +1588,10 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS, CS%OBC, Waves) fluxes%fluxes_used = .true. + if (CS%stoch_CS%do_skeb) then + call apply_skeb(CS%G,CS%GV,CS%stoch_CS,CS%u,CS%v,CS%h,CS%tv,dtdia,Time_end_thermo) + endif + if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)") ! Regridding/remapping is done here, at end of thermodynamics time step diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index debc63cb46..4bbd03a46a 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -59,6 +59,7 @@ module MOM_dynamics_split_RK2 use MOM_PressureForce, only : PressureForce, PressureForce_CS use MOM_PressureForce, only : PressureForce_init use MOM_set_visc, only : set_viscous_ML, set_visc_CS +use MOM_stochastics, only : stochastic_CS use MOM_thickness_diffuse, only : thickness_diffuse_CS use MOM_self_attr_load, only : SAL_CS use MOM_self_attr_load, only : SAL_init, SAL_end @@ -286,7 +287,7 @@ module MOM_dynamics_split_RK2 !> RK2 splitting for time stepping MOM adiabatic dynamics subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, & uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, & - MEKE, thickness_diffuse_CSp, pbv, Waves) + MEKE, thickness_diffuse_CSp, pbv, STOCH, Waves) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -326,6 +327,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s type(thickness_diffuse_CS), intent(inout) :: thickness_diffuse_CSp !< Pointer to a structure containing !! interface height diffusivities type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -851,7 +853,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & - ADp=CS%ADp) + ADp=CS%ADp, STOCH=STOCH) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index c87e6e9958..f1d3311a89 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -87,8 +87,9 @@ module MOM_dynamics_unsplit use MOM_open_boundary, only : open_boundary_zero_normal_flow use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_ML, set_visc_CS -use MOM_self_attr_load, only : SAL_init, SAL_end, SAL_CS +use MOM_stochastics, only : stochastic_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end, tidal_forcing_CS +use MOM_self_attr_load, only : SAL_init, SAL_end, SAL_CS use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_init, vertvisc_CS use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units @@ -189,7 +190,7 @@ module MOM_dynamics_unsplit !! 3rd order (for the inviscid momentum equations) order scheme subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, pbv, Waves) + VarMix, MEKE, pbv, STOCH, Waves) type(ocean_grid_type), intent(inout) :: 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 @@ -223,6 +224,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure type(MEKE_type), intent(inout) :: MEKE !< MEKE fields type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -263,7 +265,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt, Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc) + call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, STOCH=STOCH) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index b515229566..3c71fac0e4 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -87,6 +87,7 @@ module MOM_dynamics_unsplit_RK2 use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_ML, set_visc_CS use MOM_self_attr_load, only : SAL_init, SAL_end, SAL_CS +use MOM_stochastics, only : stochastic_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end, tidal_forcing_CS use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_init, vertvisc_CS @@ -192,7 +193,7 @@ module MOM_dynamics_unsplit_RK2 !> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, pbv) + VarMix, MEKE, pbv, STOCH) type(ocean_grid_type), intent(inout) :: 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 @@ -237,6 +238,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, !! fields related to the Mesoscale !! Eddy Kinetic Energy. type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av ! Averaged layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted layer thicknesses [H ~> m or kg m-2] @@ -276,7 +278,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call enable_averages(dt,Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, & - G, GV, US, CS%hor_visc) + G, GV, US, CS%hor_visc, STOCH=STOCH) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) call pass_vector(CS%diffu, CS%diffv, G%Domain, clock=id_clock_pass) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 9b1d81348e..d59e6b3871 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -20,6 +20,7 @@ module MOM_hor_visc use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE +use MOM_stochastics, only : stochastic_CS use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type use MOM_variables, only : accel_diag_ptrs @@ -241,7 +242,7 @@ module MOM_hor_visc !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT, TD, ADp) + CS, OBC, BT, TD, ADp, STOCH) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & @@ -265,6 +266,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, type(barotropic_CS), intent(in), optional :: BT !< Barotropic control structure type(thickness_diffuse_CS), intent(in), optional :: TD !< Thickness diffusion control structure type(accel_diag_ptrs), intent(in), optional :: ADp !< Acceleration diagnostics + type(stochastic_CS), intent(inout), optional :: STOCH !< Stochastic control structure ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & @@ -395,6 +397,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, logical :: apply_OBC = .false. logical :: use_MEKE_Ku logical :: use_MEKE_Au + logical :: skeb_use_frict integer :: is_vort, ie_vort, js_vort, je_vort ! Loop ranges for vorticity terms integer :: is_Kh, ie_Kh, js_Kh, je_Kh ! Loop ranges for thickness point viscosities integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -426,6 +429,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 + skeb_use_frict = .false. + if (present(STOCH)) skeb_use_frict = STOCH%skeb_use_frict + m_leithy(:,:) = 0.0 ! Initialize if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then @@ -588,12 +594,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, & !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, is_vort, ie_vort, js_vort, je_vort, & !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & - !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, & + !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, skeb_use_frict, & !$OMP backscat_subround, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & - !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt & + !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, STOCH & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, & @@ -1786,6 +1792,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) enddo ; enddo ; endif + if (skeb_use_frict) then ; do j=js,je ; do i=is,ie + ! Note that the sign convention is FrictWork < 0 means energy dissipation. + STOCH%skeb_diss(i,j,k) = STOCH%skeb_diss(i,j,k) - STOCH%skeb_frict_coef * & + FrictWork(i,j,k) / (GV%H_to_RZ * (h(i,j,k) + h_neglect)) + enddo ; enddo ; endif + ! Make a similar calculation as for FrictWork above but accumulating into ! the vertically integrated MEKE source term, and adjusting for any ! energy loss seen as a reduction in the (biharmonic) frictional source term. diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 2638ca71e1..52b55ad252 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -18,6 +18,7 @@ module MOM_thickness_diffuse use MOM_isopycnal_slopes, only : vert_fill_TS use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type +use MOM_stochastics, only : stochastic_CS use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, cont_diag_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -119,7 +120,7 @@ module MOM_thickness_diffuse !> Calculates isopycnal height diffusion coefficients and applies isopycnal height diffusion !! by modifying to the layer thicknesses, h. Diffusivities are limited to ensure stability. !! Also returns along-layer mass fluxes used in the continuity equation. -subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS) +subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS, STOCH) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -134,6 +135,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp type(VarMix_CS), target, intent(in) :: VarMix !< Variable mixing coefficients type(cont_diag_ptrs), intent(inout) :: CDp !< Diagnostics for the continuity equation type(thickness_diffuse_CS), intent(inout) :: CS !< Control structure for thickness_diffuse + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure ! Local variables real :: e(SZI_(G),SZJ_(G),SZK_(GV)+1) ! heights of interfaces, relative to mean ! sea level [Z ~> m], positive up. @@ -477,12 +479,23 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif ! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S - if (use_stored_slopes) then - call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & - int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y) + if (STOCH%skeb_use_gm) then + if (use_stored_slopes) then + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y, & + STOCH=STOCH, VarMix=VarMix) + else + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v, STOCH=STOCH, VarMix=VarMix) + endif else - call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & - int_slope_u, int_slope_v) + if (use_stored_slopes) then + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y) + else + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v) + endif endif if (VarMix%use_variable_mixing) then @@ -593,7 +606,7 @@ end subroutine thickness_diffuse !! Fluxes are limited to give positive definite thicknesses. !! Called by thickness_diffuse(). subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, & - CS, int_slope_u, int_slope_v, slope_x, slope_y) + CS, int_slope_u, int_slope_v, slope_x, slope_y, STOCH, VarMix) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -622,6 +635,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !! density gradients [nondim]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopyc. slope at u [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopyc. slope at v [Z L-1 ~> nondim] + type(stochastic_CS), optional, intent(inout) :: STOCH !< Stochastic control structure + type(VarMix_CS), target, optional, intent(in) :: VarMix !< Variable mixing coefficents ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & @@ -759,6 +774,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! [H L2 T-1 ~> m3 s-1 or kg s-1] real :: diag_sfn_unlim_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction before ! applying limiters [Z L2 T-1 ~> m3 s-1] + ! applying limiters [H L2 T-1 ~> m3 s-1 or kg s-1] + real, allocatable :: skeb_gm_work(:,:) ! Temp array to hold GM work for SKEB + real, allocatable :: skeb_ebt_norm2(:,:) ! Used to normalize EBT for SKEB + real :: h_tot ! total depth [H ~> m] + logical :: present_slope_x, present_slope_y, calc_derivatives integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of ! state calculations at u-points. @@ -766,7 +786,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! 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 + logical :: use_stanley, skeb_use_gm 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 @@ -786,6 +806,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV use_stanley = CS%use_stanley_gm + skeb_use_gm = .false. + if (present(STOCH)) skeb_use_gm = STOCH%skeb_use_gm + if (skeb_use_gm) then + allocate(skeb_gm_work(is:ie,js:je), source=0.) + allocate(skeb_ebt_norm2(is:ie,js:je), source=0.) + endif + nk_linear = max(GV%nkml, 1) Slope_x_PE(:,:,:) = 0.0 @@ -795,6 +822,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV find_work = allocated(MEKE%GM_src) find_work = (allocated(CS%GMwork) .or. find_work) + find_work = (skeb_use_gm .or. find_work) if (use_EOS) then halo = 1 ! Default halo to fill is 1 @@ -1548,8 +1576,23 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (.not. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h endif ; endif + if (skeb_use_gm) then + h_tot = sum(h(i,j,1:nz)) + skeb_gm_work(i,j) = STOCH%skeb_gm_coef * Work_h + skeb_ebt_norm2(i,j) = GV%H_to_RZ * & + (sum(h(i,j,1:nz) * VarMix%ebt_struct(i,j,1:nz)**2) + h_neglect) + endif enddo ; enddo ; endif + if (skeb_use_gm) then + ! This block spreads the GM work down through the column using the ebt vertical structure, squared. + ! Note the sign convention. + do k=1,nz ; do j=js,je ; do i=is,ie + STOCH%skeb_diss(i,j,k) = STOCH%skeb_diss(i,j,k) - skeb_gm_work(i,j) * & + VarMix%ebt_struct(i,j,k)**2 / skeb_ebt_norm2(i,j) + enddo ; enddo ; enddo + endif + if (find_work .and. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then do j=js,je ; do i=is,ie ; do k=nz,1,-1 PE_release_h = -0.25 * (GV%H_to_RZ*US%L_to_Z**2) * & diff --git a/src/parameterizations/stochastic/MOM_stochastics.F90 b/src/parameterizations/stochastic/MOM_stochastics.F90 index 04a29019fa..1bc42660d3 100644 --- a/src/parameterizations/stochastic/MOM_stochastics.F90 +++ b/src/parameterizations/stochastic/MOM_stochastics.F90 @@ -8,8 +8,12 @@ module MOM_stochastics ! particular version wraps all of the calls for MOM6 in the calls that had ! been used for MOM4. ! -use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type +use MOM_debugging, only : hchksum, uvchksum, qchksum +use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type, post_data +use MOM_diag_mediator, only : register_static_field, enable_averages, disable_averaging use MOM_grid, only : ocean_grid_type +use MOM_variables, only : thermo_var_ptrs +use MOM_domains, only : pass_var, pass_vector, CORNER, SCALAR_PAIR use MOM_verticalGrid, only : verticalGrid_type use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave @@ -18,28 +22,56 @@ module MOM_stochastics use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use MOM_domains, only : root_PE, num_PEs use MOM_coms, only : Get_PElist +use MOM_EOS, only : calculate_density, EOS_domain use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn #include implicit none ; private -public stochastics_init, update_stochastics +public stochastics_init, update_stochastics, apply_skeb !> This control structure holds parameters for the MOM_stochastics module type, public:: stochastic_CS logical :: do_sppt !< If true, stochastically perturb the diabatic + logical :: do_skeb !< If true, stochastically perturb the diabatic + logical :: skeb_use_gm !< If true, adds GM work to the amplitude of SKEBS + logical :: skeb_use_frict !< If true, adds viscous dissipation rate to the amplitude of SKEBS logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and genration terms - integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT - integer :: id_epbl1_wts = -1 !< Diagnostic id for epbl generation perturbation - integer :: id_epbl2_wts = -1 !< Diagnostic id for epbl dissipation perturbation + integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT + integer :: id_skeb_wts = -1 !< Diagnostic id for SKEB + integer :: id_skebu = -1 !< Diagnostic id for SKEB + integer :: id_skebv = -1 !< Diagnostic id for SKEB + integer :: id_diss = -1 !< Diagnostic id for SKEB + integer :: skeb_npass = -1 !< number of passes of the 9-point smoother for the dissipation estimate + integer :: id_psi = -1 !< Diagnostic id for SPPT + integer :: id_epbl1_wts = -1 !< Diagnostic id for epbl generation perturbation + integer :: id_epbl2_wts = -1 !< Diagnostic id for epbl dissipation perturbation + integer :: id_skeb_taperu = -1 !< Diagnostic id for u taper of SKEB velocity increment + integer :: id_skeb_taperv = -1 !< Diagnostic id for v taper of SKEB velocity increment + real :: skeb_gm_coef !< If skeb_use_gm is true, then skeb_gm_coef * GM_work is added to the + !! dissipation rate used to set the amplitude of SKEBS [nondim] + real :: skeb_frict_coef !< If skeb_use_frict is true, then skeb_gm_coef * GM_work is added to the + !! dissipation rate used to set the amplitude of SKEBS [nondim] + real, allocatable :: skeb_diss(:,:,:) !< Dissipation rate used to set amplitude of SKEBS [L2 T-3 ~> m2 s-2] + !! Index into this at h points. ! stochastic patterns real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT !! tendencies with a number between 0 and 2 + real, allocatable :: skeb_wts(:,:) !< Random pattern for ocean SKEB real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation - type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output type(time_type), pointer :: Time !< Pointer to model time (needed for sponges) + type(diag_ctrl), pointer :: diag=>NULL() !< A structure that is used to regulate the + + ! Taper array to smoothly zero out the SKEBS velocity increment near land + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: taperCu !< Taper applied to u component of + !! stochastic velocity increment + !! range [0,1], [nondim] + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: taperCv !< Taper applied to v component of + !! stochastic velocity increment + !! range [0,1], [nondim] + end type stochastic_CS contains @@ -62,20 +94,24 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) integer :: pe_zero ! root pe integer :: nx ! number of x-points including halo integer :: ny ! number of x-points including halo + integer :: i, j, k ! loop indices + real :: tmp(grid%isdB:grid%iedB,grid%jsdB:grid%jedB) ! Used to construct tapers + integer :: taper_width ! Width (in cells) of the taper that brings the stochastic velocity + ! increments to 0 at the boundary. ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "ocean_stochastics_init" ! This module's name. - call callTree_enter("ocean_model_stochastic_init(), MOM_stochastics.F90") + call callTree_enter("stochastic_init(), MOM_stochastics.F90") if (associated(CS)) then call MOM_error(WARNING, "MOM_stochastics_init called with an "// & "associated control structure.") return else ; allocate(CS) ; endif - CS%diag => diag CS%Time => Time + CS%diag => diag ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") @@ -83,48 +119,130 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) ! get number of processors and PE list for stochastic physics initialization call get_param(param_file, mdl, "DO_SPPT", CS%do_sppt, & "If true, then stochastically perturb the thermodynamic "//& - "tendemcies of T,S, amd h. Amplitude and correlations are "//& + "tendencies of T,S, amd h. Amplitude and correlations are "//& "controlled by the nam_stoch namelist in the UFS model only.", & default=.false.) + call get_param(param_file, mdl, "DO_SKEB", CS%do_skeb, & + "If true, then stochastically perturb the currents "//& + "using the stochastic kinetic energy backscatter scheme.",& + default=.false.) + call get_param(param_file, mdl, "SKEB_NPASS", CS%skeb_npass, & + "number of passes of a 9-point smoother of the "//& + "dissipation estimate.", default=3, do_not_log=.not.CS%do_skeb) + call get_param(param_file, mdl, "SKEB_TAPER_WIDTH", taper_width, & + "number of cells over which the stochastic velocity increment "//& + "is tapered to zero.", default=4, do_not_log=.not.CS%do_skeb) + call get_param(param_file, mdl, "SKEB_USE_GM", CS%skeb_use_gm, & + "If true, adds GM work rate to the SKEBS amplitude.", & + default=.false., do_not_log=.not.CS%do_skeb) + if ((.not. CS%do_skeb) .and. (CS%skeb_use_gm)) call MOM_error(FATAL, "If SKEB_USE_GM is True "//& + "then DO_SKEB must also be True.") + call get_param(param_file, mdl, "SKEB_GM_COEF", CS%skeb_gm_coef, & + "Fraction of GM work that is added to backscatter rate.", & + units="nondim", default=0.0, do_not_log=.not.CS%skeb_use_gm) + call get_param(param_file, mdl, "SKEB_USE_FRICT", CS%skeb_use_frict, & + "If true, adds horizontal friction dissipation rate "//& + "to the SKEBS amplitude.", default=.false., do_not_log=.not.CS%do_skeb) + if ((.not. CS%do_skeb) .and. (CS%skeb_use_frict)) call MOM_error(FATAL, "If SKEB_USE_FRICT is "//& + "True then DO_SKEB must also be True.") + call get_param(param_file, mdl, "SKEB_FRICT_COEF", CS%skeb_frict_coef, & + "Fraction of horizontal friction work that is added to backscatter rate.", & + units="nondim", default=0.0, do_not_log=.not.CS%skeb_use_frict) call get_param(param_file, mdl, "PERT_EPBL", CS%pert_epbl, & "If true, then stochastically perturb the kinetic energy "//& "production and dissipation terms. Amplitude and correlations are "//& "controlled by the nam_stoch namelist in the UFS model only.", & default=.false.) - if (CS%do_sppt .OR. CS%pert_epbl) then - num_procs = num_PEs() - allocate(pelist(num_procs)) - call Get_PElist(pelist,commID = mom_comm) - pe_zero = root_PE() - nx = grid%ied - grid%isd + 1 - ny = grid%jed - grid%jsd + 1 - call init_stochastic_physics_ocn(dt,grid%geoLonT,grid%geoLatT,nx,ny,GV%ke, & - CS%pert_epbl,CS%do_sppt,pe_zero,mom_comm,iret) - if (iret/=0) then - call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed") - endif - - if (CS%do_sppt) allocate(CS%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed), source=0.0) - if (CS%pert_epbl) then - allocate(CS%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed), source=0.0) - allocate(CS%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed), source=0.0) - endif - endif - if (CS%do_sppt) then - CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesT1, Time, & - 'random pattern for sppt', 'None') + + if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) then + num_procs = num_PEs() + allocate(pelist(num_procs)) + call Get_PElist(pelist,commID = mom_comm) + pe_zero = root_PE() + nx = grid%iedB - grid%isdB + 1 + ny = grid%jedB - grid%jsdB + 1 + call init_stochastic_physics_ocn(dt,grid%geoLonBu,grid%geoLatBu,nx,ny,GV%ke, & + CS%pert_epbl,CS%do_sppt,CS%do_skeb,pe_zero,mom_comm,iret) + if (iret/=0) then + call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed") + return + endif + + if (CS%do_sppt) allocate(CS%sppt_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + if (CS%do_skeb) allocate(CS%skeb_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + if (CS%do_skeb) allocate(CS%skeb_diss(grid%isd:grid%ied,grid%jsd:grid%jed,GV%ke), source=0.) + if (CS%pert_epbl) then + allocate(CS%epbl1_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + allocate(CS%epbl2_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + endif endif - if (CS%pert_epbl) then - CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, & - 'random pattern for KE generation', 'None') - CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, & - 'random pattern for KE dissipation', 'None') + + CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesB1, Time, & + 'random pattern for sppt', 'None') + CS%id_skeb_wts = register_diag_field('ocean_model', 'skeb_pattern', CS%diag%axesB1, Time, & + 'random pattern for skeb', 'None') + CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesB1, Time, & + 'random pattern for KE generation', 'None') + CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesB1, Time, & + 'random pattern for KE dissipation', 'None') + CS%id_skebu = register_diag_field('ocean_model', 'skebu', CS%diag%axesCuL, Time, & + 'zonal current perts', 'None') + CS%id_skebv = register_diag_field('ocean_model', 'skebv', CS%diag%axesCvL, Time, & + 'zonal current perts', 'None') + CS%id_diss = register_diag_field('ocean_model', 'skeb_amp', CS%diag%axesTL, Time, & + 'SKEB amplitude', 'm s-1') + CS%id_psi = register_diag_field('ocean_model', 'psi', CS%diag%axesBL, Time, & + 'stream function', 'None') + CS%id_skeb_taperu = register_static_field('ocean_model', 'skeb_taper_u', CS%diag%axesCu1, & + 'SKEB taper u', 'None', interp_method='none') + CS%id_skeb_taperv = register_static_field('ocean_model', 'skeb_taper_v', CS%diag%axesCv1, & + 'SKEB taper v', 'None', interp_method='none') + + ! Initialize the "taper" fields. These fields multiply the components of the stochastic + ! velocity increment in such a way as to smoothly taper them to zero at land boundaries. + if ((CS%do_skeb) .or. (CS%id_skeb_taperu > 0) .or. (CS%id_skeb_taperv > 0)) then + ALLOC_(CS%taperCu(grid%IsdB:grid%IedB,grid%jsd:grid%jed)) + ALLOC_(CS%taperCv(grid%isd:grid%ied,grid%JsdB:grid%JedB)) + ! Initialize taper from land mask + do j=grid%jsd,grid%jed ; do I=grid%isdB,grid%iedB + CS%taperCu(I,j) = grid%mask2dCu(I,j) + enddo ; enddo + do J=grid%jsdB,grid%jedB ; do i=grid%isd,grid%ied + CS%taperCv(i,J) = grid%mask2dCv(i,J) + enddo ; enddo + ! Extend taper land + do k=1,(taper_width / 2) + do j=grid%jsc-1,grid%jec+1 ; do I=grid%iscB-1,grid%iecB+1 + tmp(I,j) = minval(CS%taperCu(I-1:I+1,j-1:j+1)) + enddo ; enddo + do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB + CS%taperCu(I,j) = minval(tmp(I-1:I+1,j-1:j+1)) + enddo ; enddo + do J=grid%jscB-1,grid%jecB+1 ; do i=grid%isc-1,grid%iec+1 + tmp(i,J) = minval(CS%taperCv(i-1:i+1,J-1:J+1)) + enddo ; enddo + do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec + CS%taperCv(i,J) = minval(tmp(i-1:i+1,J-1:J+1)) + enddo ; enddo + ! Update halo + call pass_vector(CS%taperCu, CS%taperCv, grid%Domain, SCALAR_PAIR) + enddo + ! Smooth tapers. Each call smooths twice. + do k=1,(taper_width - (taper_width/2)) + call smooth_x9_uv(grid, CS%taperCu, CS%taperCv, zero_land=.true.) + call pass_vector(CS%taperCu, CS%taperCv, grid%Domain, SCALAR_PAIR) + enddo endif - if (CS%do_sppt .OR. CS%pert_epbl) & + !call uvchksum("SKEB taper [uv]", CS%taperCu, CS%taperCv, grid%HI) + + if (CS%id_skeb_taperu > 0) call post_data(CS%id_skeb_taperu, CS%taperCu, CS%diag, .true.) + if (CS%id_skeb_taperv > 0) call post_data(CS%id_skeb_taperv, CS%taperCv, CS%diag, .true.) + + if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) & call MOM_mesg(' === COMPLETED MOM STOCHASTIC INITIALIZATION =====') - call callTree_leave("ocean_model_init(") + call callTree_leave("stochastic_init(), MOM_stochastics.F90") end subroutine stochastics_init @@ -138,10 +256,202 @@ subroutine update_stochastics(CS) call callTree_enter("update_stochastics(), MOM_stochastics.F90") ! update stochastic physics patterns before running next time-step - call run_stochastic_physics_ocn(CS%sppt_wts,CS%epbl1_wts,CS%epbl2_wts) + call run_stochastic_physics_ocn(CS%sppt_wts,CS%skeb_wts,CS%epbl1_wts,CS%epbl2_wts) + + call callTree_leave("update_stochastics(), MOM_stochastics.F90") - return end subroutine update_stochastics +subroutine apply_skeb(grid,GV,CS,uc,vc,thickness,tv,dt,Time_end) + + type(ocean_grid_type), intent(in) :: grid !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid + type(stochastic_CS), intent(inout) :: CS !< stochastic control structure + + real, dimension(SZIB_(grid),SZJ_(grid),SZK_(GV)), intent(inout) :: uc !< zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(grid),SZJB_(grid),SZK_(GV)), intent(inout) :: vc !< meridional velocity [L T-1 ~> m s-1] + real, dimension(SZI_(grid),SZJ_(grid),SZK_(GV)), intent(in) :: thickness !< thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< points to thermodynamic fields + real, intent(in) :: dt !< time increment [T ~> s] + type(time_type), intent(in) :: Time_end !< Time at the end of the interval +! locals + + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_,NKMEM_) :: psi + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: ustar + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: vstar + real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: diss_tmp + + real, dimension(3,3) :: local_weights + + real :: shr,ten,tot,kh + integer :: i,j,k,iter + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + + call callTree_enter("apply_skeb(), MOM_stochastics.F90") + ALLOC_(diss_tmp(grid%isd:grid%ied,grid%jsd:grid%jed)) + ALLOC_(psi(grid%isdB:grid%iedB,grid%jsdB:grid%jedB,GV%ke)) + ALLOC_(ustar(grid%isdB:grid%iedB,grid%jsd:grid%jed,GV%ke)) + ALLOC_(vstar(grid%isd:grid%ied,grid%jsdB:grid%jedB,GV%ke)) + + if ((.not. CS%skeb_use_gm) .and. (.not. CS%skeb_use_frict)) then + ! fill in halos with zeros + do k=1,GV%ke + do j=grid%jsd,grid%jed ; do i=grid%isd,grid%ied + CS%skeb_diss(i,j,k) = 0.0 + enddo ; enddo + enddo + + !kh needs to be scaled + + kh=1!(120*111)**2 + do k=1,GV%ke + do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec + ! Shear + shr = (vc(i,J,k)-vc(i-1,J,k))*grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdxCv(i,J)+& + (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdyCu(I,j) + ! Tension + ten = (vc(i,J,k)-vc(i-1,J,k))*grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdyCv(i,J)+& + (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdxCu(I,j) + + tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j) + CS%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j)!!**2 + enddo ; enddo + enddo + endif ! Sets CS%skeb_diss without GM or FrictWork + + ! smooth dissipation skeb_npass times + do iter=1,CS%skeb_npass + if (mod(iter,2) == 1) call pass_var(CS%skeb_diss, grid%domain) + do k=1,GV%ke + do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 + ! This does not preserve rotational symmetry + local_weights = grid%mask2dT(i-1:i+1,j-1:j+1)*grid%areaT(i-1:i+1,j-1:j+1) + diss_tmp(i,j) = sum(local_weights*CS%skeb_diss(i-1:i+1,j-1:j+1,k)) / & + (sum(local_weights) + 1.E-16) + enddo ; enddo + do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 + if (grid%mask2dT(i,j)==0.) cycle + CS%skeb_diss(i,j,k) = diss_tmp(i,j) + enddo ; enddo + enddo + enddo + call pass_var(CS%skeb_diss, grid%domain) + + ! call hchksum(CS%skeb_diss, "SKEB DISS", grid%HI, haloshift=2) + ! call qchksum(CS%skeb_wts, "SKEB WTS", grid%HI, haloshift=1) + + do k=1,GV%ke + do J=grid%jscB-1,grid%jecB ; do I=grid%iscB-1,grid%iecB + psi(I,J,k) = sqrt(0.25 * dt * max((CS%skeb_diss(i ,j ,k) + CS%skeb_diss(i+1,j+1,k)) + & + (CS%skeb_diss(i ,j+1,k) + CS%skeb_diss(i+1,j ,k)), 0.) ) & + * CS%skeb_wts(I,J) + enddo ; enddo + enddo + !call qchksum(psi,"SKEB PSI", grid%HI, haloshift=1) + !call pass_var(psi, grid%domain, position=CORNER) + do k=1,GV%ke + do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB + ustar(I,j,k) = - (psi(I,J,k) - psi(I,J-1,k)) * CS%taperCu(I,j) * grid%IdyCu(I,j) + uc(I,j,k) = uc(I,j,k) + ustar(I,j,k) + enddo ; enddo + do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec + vstar(i,J,k) = (psi(I,J,k) - psi(I-1,J,k)) * CS%taperCv(i,J) * grid%IdxCv(i,J) + vc(i,J,k) = vc(i,J,k) + vstar(i,J,k) + enddo ; enddo + enddo + + !call uvchksum("SKEB increment [uv]", ustar, vstar, grid%HI) + + call enable_averages(dt, Time_end, CS%diag) + if (CS%id_diss > 0) then + call post_data(CS%id_diss, sqrt(dt * max(CS%skeb_diss(:,:,:), 0.)), CS%diag) + endif + if (CS%id_skeb_wts > 0) then + call post_data(CS%id_skeb_wts, CS%skeb_wts, CS%diag) + endif + if (CS%id_skebu > 0) then + call post_data(CS%id_skebu, ustar(:,:,:), CS%diag) + endif + if (CS%id_skebv > 0) then + call post_data(CS%id_skebv, vstar(:,:,:), CS%diag) + endif + if (CS%id_psi > 0) then + call post_data(CS%id_psi, psi(:,:,:), CS%diag) + endif + call disable_averaging(CS%diag) + DEALLOC_(diss_tmp) + DEALLOC_(ustar) + DEALLOC_(vstar) + DEALLOC_(psi) + CS%skeb_diss(:,:,:) = 0.0 ! Must zero before next time step. + + call callTree_leave("apply_skeb(), MOM_stochastics.F90") + +end subroutine apply_skeb + +!> Apply a 9-point smoothing filter twice to a pair of velocity components to reduce +!! horizontal two-grid-point noise. +!! Note that this subroutine does not conserve angular momentum, so don't use it +!! in situations where you need conservation. Also note that it assumes that the +!! input fields have valid values in the first two halo points upon entry. +subroutine smooth_x9_uv(G, field_u, field_v, zero_land) + type(ocean_grid_type), intent(in) :: G !< Ocean grid + real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed[arbitrary] + real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary] + logical, optional, intent(in) :: zero_land !< If present and false, return the average + !! of the surrounding ocean points when + !! smoothing, otherwise use a value of 0 for + !! land points and include them in the averages. + + ! Local variables. + real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary] + real :: fv_prev(SZI_(G),SZJB_(G)) ! The value of the v-point field at the previous iteration [arbitrary] + real :: Iwts ! The inverse of the sum of the weights [nondim] + logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent. + integer :: i, j, s, is, ie, js, je, Isq, Ieq, Jsq, Jeq + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + + zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land + + do s=1,0,-1 + fu_prev(:,:) = field_u(:,:) + ! apply smoothing on field_u using rotationally symmetric expressions. + do j=js-s,je+s ; do I=Isq-s,Ieq+s ; if (G%mask2dCu(I,j) > 0.0) then + Iwts = 0.0625 + if (.not. zero_land_val) & + Iwts = 1.0 / ( (4.0*G%mask2dCu(I,j) + & + ( 2.0*((G%mask2dCu(I-1,j) + G%mask2dCu(I+1,j)) + & + (G%mask2dCu(I,j-1) + G%mask2dCu(I,j+1))) + & + ((G%mask2dCu(I-1,j-1) + G%mask2dCu(I+1,j+1)) + & + (G%mask2dCu(I-1,j+1) + G%mask2dCu(I+1,j-1))) ) ) + 1.0e-16 ) + field_u(I,j) = Iwts * ( 4.0*G%mask2dCu(I,j) * fu_prev(I,j) & + + (2.0*((G%mask2dCu(I-1,j) * fu_prev(I-1,j) + G%mask2dCu(I+1,j) * fu_prev(I+1,j)) + & + (G%mask2dCu(I,j-1) * fu_prev(I,j-1) + G%mask2dCu(I,j+1) * fu_prev(I,j+1))) & + + ((G%mask2dCu(I-1,j-1) * fu_prev(I-1,j-1) + G%mask2dCu(I+1,j+1) * fu_prev(I+1,j+1)) + & + (G%mask2dCu(I-1,j+1) * fu_prev(I-1,j+1) + G%mask2dCu(I+1,j-1) * fu_prev(I-1,j-1))) )) + endif ; enddo ; enddo + + fv_prev(:,:) = field_v(:,:) + ! apply smoothing on field_v using rotationally symmetric expressions. + do J=Jsq-s,Jeq+s ; do i=is-s,ie+s ; if (G%mask2dCv(i,J) > 0.0) then + Iwts = 0.0625 + if (.not. zero_land_val) & + Iwts = 1.0 / ( (4.0*G%mask2dCv(i,J) + & + ( 2.0*((G%mask2dCv(i-1,J) + G%mask2dCv(i+1,J)) + & + (G%mask2dCv(i,J-1) + G%mask2dCv(i,J+1))) + & + ((G%mask2dCv(i-1,J-1) + G%mask2dCv(i+1,J+1)) + & + (G%mask2dCv(i-1,J+1) + G%mask2dCv(i+1,J-1))) ) ) + 1.0e-16 ) + field_v(i,J) = Iwts * ( 4.0*G%mask2dCv(i,J) * fv_prev(i,J) & + + (2.0*((G%mask2dCv(i-1,J) * fv_prev(i-1,J) + G%mask2dCv(i+1,J) * fv_prev(i+1,J)) + & + (G%mask2dCv(i,J-1) * fv_prev(i,J-1) + G%mask2dCv(i,J+1) * fv_prev(i,J+1))) & + + ((G%mask2dCv(i-1,J-1) * fv_prev(i-1,J-1) + G%mask2dCv(i+1,J+1) * fv_prev(i+1,J+1)) + & + (G%mask2dCv(i-1,J+1) * fv_prev(i-1,J+1) + G%mask2dCv(i+1,J-1) * fv_prev(i-1,J-1))) )) + endif ; enddo ; enddo + enddo + +end subroutine smooth_x9_uv + end module MOM_stochastics