Skip to content

Commit

Permalink
Added SQG vertical structure in Varmix to provide vertical profile fo…
Browse files Browse the repository at this point in the history 
…r diffusivities

- added function calc_sqg_struct in MOM_lateral_mixing_coeffs to compute sqg_struct
- added sqg_expo to set the exponent of sqg_struct
- set BS_use_sqg=true, sqg_expo>0., and BS_EBT_power=0. to use sqg_struct for the backscatter
- if sqg_use_MEKE=true, use the eddy length scale from MEKE to compute sqg_struct
- added eddy length scale Le in MEKE
- added MEKE in Varmix
- registered N2_u and N2_v diagnostics when SQG_EXPO>0

(cherry picked from commit 6d3df0541c33d6f6d1f9fcb695f1a1eb961ec1b3)
  • Loading branch information
@Wendazhang33
Wendazhang33 committed Aug 9, 2024
1 parent e27eec5 commit 1ac5a75
Show file tree
Hide file tree
Showing 5 changed files with 135 additions and 12 deletions.
6 changes: 3 additions & 3 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -750,7 +750,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS

if (CS%VarMix%use_variable_mixing) then
call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag)
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix)
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE)
call calc_depth_function(G, CS%VarMix)
call disable_averaging(CS%diag)
endif
Expand Down Expand Up @@ -1898,7 +1898,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
if (.not. skip_diffusion) then
if (CS%VarMix%use_variable_mixing) then
call pass_var(CS%h, G%Domain)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE)
call calc_depth_function(G, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
endif
Expand All @@ -1925,7 +1925,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
if (.not. skip_diffusion) then
if (CS%VarMix%use_variable_mixing) then
call pass_var(CS%h, G%Domain)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE)
call calc_depth_function(G, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
endif
Expand Down
7 changes: 7 additions & 0 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -397,6 +397,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m)
endif

!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
MEKE%Le(i,j) = LmixScale(i,j)
enddo ; enddo

! Aggregate sources of MEKE (background, frictional and GM)
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
Expand Down Expand Up @@ -1841,6 +1846,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
call MOM_mesg("MEKE_alloc_register_restart: allocating and registering", 5)
isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed
allocate(MEKE%MEKE(isd:ied,jsd:jed), source=0.0)
allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0)
call register_restart_field(MEKE%MEKE, "MEKE", .false., restart_CS, &
longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2)

Expand Down Expand Up @@ -1899,6 +1905,7 @@ subroutine MEKE_end(MEKE)
if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh)
if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src)
if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE)
if (allocated(MEKE%Le)) deallocate(MEKE%Le)
end subroutine MEKE_end

!> \namespace mom_meke
Expand Down
1 change: 1 addition & 0 deletions src/parameterizations/lateral/MOM_MEKE_types.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module MOM_MEKE_types
!! backscatter from unresolved eddies (see Jansen and Held, 2014).
real, allocatable :: Au(:,:) !< The MEKE-derived lateral biharmonic viscosity
!! coefficient [L4 T-1 ~> m4 s-1].
real, allocatable :: Le(:,:) !< Eddy length scale [L m]

! Parameters
real :: KhTh_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTh [nondim]
Expand Down
6 changes: 4 additions & 2 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -659,7 +659,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
!$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
!$OMP is_vort, ie_vort, js_vort, je_vort, &
!$OMP is_Kh, ie_Kh, js_Kh, je_Kh, &
!$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
!$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, EY24_EBT_BS, &
!$OMP BS_coeff_h, visc_limit_h, visc_limit_q, visc_limit_h_flag, visc_limit_q_flag, &
!$OMP visc_limit_h_frac, visc_limit_q_frac, &
!$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, &
!$OMP backscat_subround, GME_effic_h, GME_effic_q, &
!$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, &
Expand All @@ -678,7 +680,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
!$OMP sh_xx_sq, sh_xy_sq, meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, &
!$OMP h_min, hrat_min, visc_bound_rem, Kh_max_here, &
!$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, &
!$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, &
!$OMP Kh, Kh_BS, Ah, AhSm, AhLth, local_strain, Sh_F_pow, tmp, &
!$OMP dDel2vdx, dDel2udy, Del2vort_q, Del2vort_h, KE, &
!$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff, &
!$OMP dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, &
Expand Down
127 changes: 120 additions & 7 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@ module MOM_lateral_mixing_coeffs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init
use MOM_open_boundary, only : ocean_OBC_type
use MOM_MEKE_types, only : MEKE_type


implicit none ; private

Expand Down Expand Up @@ -112,9 +114,13 @@ module MOM_lateral_mixing_coeffs

real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [Z L-1 ~> nondim]
real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [Z L-1 ~> nondim]
real, allocatable :: ebt_struct(:,:,:) !< Vertical structure function to scale diffusivities with [nondim]
real, allocatable :: ebt_struct(:,:,:) !< EBT vertical structure to scale diffusivities with [nondim]
real, allocatable :: sqg_struct(:,:,:) !< SQG vertical structure to scale diffusivities with [nondim]
real, allocatable :: BS_struct(:,:,:) !< Vertical structure function used in backscatter [nondim]
real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0.
real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 0.0
logical :: sqg_use_MEKE !< If true, use MEKE%Le for SQG vertical structure.
logical :: BS_use_sqg !< If true, use sqg_stuct for backscatter vertical structure.


real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
Expand Down Expand Up @@ -163,6 +169,7 @@ module MOM_lateral_mixing_coeffs
integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1
integer :: id_dzu=-1, id_dzv=-1, id_dzSxN=-1, id_dzSyN=-1
integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1
integer :: id_sqg_struct=-1, id_BS_struct=-1
type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
!! timing of diagnostic output.
!>@}
Expand All @@ -173,7 +180,7 @@ module MOM_lateral_mixing_coeffs
end type VarMix_CS

public VarMix_init, VarMix_end, calc_slope_functions, calc_resoln_function
public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function
public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function, calc_sqg_struct

contains

Expand Down Expand Up @@ -214,13 +221,14 @@ subroutine calc_depth_function(G, CS)
end subroutine calc_depth_function

!> Calculates and stores the non-dimensional resolution functions
subroutine calc_resoln_function(h, tv, G, GV, US, CS)
subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
type(MEKE_type), intent(in) :: MEKE !< MEKE struct

! Local variables
! Depending on the power-function being used, dimensional rescaling may be limited, so some
Expand All @@ -230,6 +238,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
real :: cg1_v ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1].
real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2]
integer :: power_2
real :: dt !< Time increment [T ~> s]
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: i, j, k
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Expand Down Expand Up @@ -261,11 +270,20 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain)
call do_group_pass(CS%pass_cg1, G%Domain)
endif
if (CS%sqg_expo>0.0) then
call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
call pass_var(CS%sqg_struct, G%Domain)
endif

if (CS%BS_EBT_power>0.) then
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
CS%BS_struct(i,j,k) = CS%ebt_struct(i,j,k)**CS%BS_EBT_power
enddo ; enddo ; enddo
endif
elseif (CS%BS_use_sqg) then
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
CS%BS_struct(i,j,k) = CS%sqg_struct(i,j,k)
enddo ; enddo ; enddo
endif

! Calculate and store the ratio between deformation radius and grid-spacing
! at h-points [nondim].
Expand Down Expand Up @@ -460,6 +478,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)

if (query_averaging_enabled(CS%diag)) then
if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag)
if (CS%id_BS_struct > 0) call post_data(CS%id_BS_struct, CS%BS_struct, CS%diag)
endif

if (CS%debug) then
Expand All @@ -470,6 +489,75 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)

end subroutine calc_resoln_function

!> Calculates and stores functions of SQG mode
subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
type(ocean_grid_type), intent(inout) :: 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
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !<Thermodynamic variables
real, intent(in) :: dt !< Time increment [T ~> s]
type(VarMix_CS), intent(inout) :: CS !< Variable mixing control struct
type(MEKE_type), intent(in) :: MEKE !< MEKE struct
! type(ocean_OBC_type), pointer :: OBC !< Open
! boundaries control structure.

! Local variables
real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: &
e ! The interface heights relative to mean sea level [Z ~> m].
real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2]
real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2]
real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1]
real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1]
real, dimension(SZI_(G), SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1]
real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2]
real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m]
! real, dimension(SZK_(GV)) :: zs ! Stretched vertical coordinate [m]
real, dimension(SZI_(G), SZJ_(G)) :: Le ! Eddy length scale [m]
integer :: i, j, k, is, ie, js, je, nz

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
"Module must be initialized before it is used.")

call find_eta(h, tv, G, GV, US, e, halo_size=2)
call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v,dzu=dzu, dzv=dzv, &
dzSxN=dzSxN, dzSyN=dzSyN, halo=1)

do j=js,je ; do i=is,ie
CS%sqg_struct(i,j,1) = 1.0
enddo ; enddo
if (CS%sqg_use_MEKE) then
do j=js,je ; do i=is,ie
Le(i,j) = MEKE%Le(i,j)
f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8)
enddo ; enddo
else
do j=js,je ; do i=is,ie
Le(i,j) = sqrt(G%areaT(i,j))
f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8)
enddo ; enddo
endif
do k=2,nz ; do j=js,je ; do i=is,ie
N2 = max(0.25*(N2_u(I-1,j,k) + N2_u(I,j,k) + N2_v(i,J-1,k) + N2_v(i,J,k)),0.0)
dzc = 0.25*(dzu(I-1,j,k) + dzu(I,j,k) + dzv(i,J-1,k) + dzv(i,J,k))*N2**0.5/f(i,j)
! dzs = -N2**0.5/f(i,j)*dzc
CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1)*exp(-CS%sqg_expo*dzc/Le(i,j))
enddo ; enddo ; enddo


if (query_averaging_enabled(CS%diag)) then
if (CS%id_sqg_struct > 0) call post_data(CS%id_sqg_struct, CS%sqg_struct, CS%diag)
if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
endif

end subroutine calc_sqg_struct

!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
!! style scaling of diffusivity
subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
Expand Down Expand Up @@ -1260,6 +1348,18 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, &
"Power to raise EBT vertical structure to when backscatter "// &
"has vertical structure.", units="nondim", default=0.0)
call get_param(param_file, mdl, "BS_USE_SQG", CS%BS_use_sqg, &
"If true, the SQG vertical structure is used for backscatter "//&
"on the condition that BS_EBT_power=0", &
default=.false.)
if (CS%BS_EBT_power>0.) CS%BS_use_sqg = .false.
call get_param(param_file, mdl, "SQG_EXPO", CS%sqg_expo, &
"Nondimensional exponent coeffecient of the SQG mode "// &
"that is used for the vertical struture of diffusivities.", units="nondim", default=0.0)
if (CS%sqg_expo==0.) CS%BS_use_sqg = .false.
call get_param(param_file, mdl, "SQG_USE_MEKE", CS%sqg_use_MEKE, &
"If true, the eddy scale of MEKE is used for the SQG vertical structure ",&
default=.false.)
call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, &
"If true, uses the equivalent barotropic structure "//&
"as the vertical structure of thickness diffusivity.",&
Expand Down Expand Up @@ -1320,6 +1420,10 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
units="m", default=-1.0, scale=GV%m_to_H)
allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
endif
if (CS%sqg_expo>0.0) then
allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
endif

allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
CS%BS_struct(:,:,:) = 1.0

Expand All @@ -1335,7 +1439,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
endif
endif

if (CS%use_stored_slopes) then
if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then
! CS%calculate_Eady_growth_rate=.true.
in_use = .true.
allocate(CS%slope_x(IsdB:IedB,jsd:jed,GV%ke+1), source=0.0)
Expand Down Expand Up @@ -1418,7 +1522,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
'm2', conversion=US%L_to_m**2)
endif

if (CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) then
if (CS%sqg_expo>0.0) then
CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, &
'Vertical structure of SQG mode', 'nondim')
endif

CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, &
'Vertical structure of backscatter', 'nondim')

if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) .or. (CS%sqg_expo>0.0)) then
CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, &
'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
's-2', conversion=(US%L_to_Z*US%s_to_T)**2)
Expand Down Expand Up @@ -1672,9 +1784,10 @@ subroutine VarMix_end(CS)

if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.) &
deallocate(CS%ebt_struct)
if (CS%sqg_expo>0.0) deallocate(CS%sqg_struct)
if (allocated(CS%BS_struct)) deallocate(CS%BS_struct)

if (CS%use_stored_slopes) then
if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then
deallocate(CS%slope_x)
deallocate(CS%slope_y)
endif
Expand Down

0 comments on commit 1ac5a75

Please sign in to comment.