Skip to content

Commit a9dd5f2

Browse files
Wendazhang33Wendazhang33
committed
Added SQG vertical structure in Varmix to provide vertical profile for 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 - to use sqg_struct for the backscatter, set BS_use_sqg=true, sqg_expo>0., and BS_EBT_power=0. - if SQG_USE_MEKE=True, use the eddy length scale from MEKE to compute sqg_struct - added eddy length scale Le in MEKE if SQG_USE_MEKE=True - added MEKE%Le into restart file if SQG_USE_MEKE=True - added MEKE in Varmix - registered N2_u and N2_v diagnostics when SQG_EXPO>0 (cherry picked from commit 6d3df0541c33d6f6d1f9fcb695f1a1eb961ec1b3)
1 parent f90b071 commit a9dd5f2

File tree

4 files changed

+153
-11
lines changed

4 files changed

+153
-11
lines changed

‎src/core/MOM.F90‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -750,7 +750,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
750750

751751
if (CS%VarMix%use_variable_mixing) then
752752
call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag)
753-
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix)
753+
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE)
754754
call calc_depth_function(G, CS%VarMix)
755755
call disable_averaging(CS%diag)
756756
endif
@@ -1898,7 +1898,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
18981898
if (.not. skip_diffusion) then
18991899
if (CS%VarMix%use_variable_mixing) then
19001900
call pass_var(CS%h, G%Domain)
1901-
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
1901+
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE)
19021902
call calc_depth_function(G, CS%VarMix)
19031903
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
19041904
endif
@@ -1925,7 +1925,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
19251925
if (.not. skip_diffusion) then
19261926
if (CS%VarMix%use_variable_mixing) then
19271927
call pass_var(CS%h, G%Domain)
1928-
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
1928+
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE)
19291929
call calc_depth_function(G, CS%VarMix)
19301930
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
19311931
endif

‎src/parameterizations/lateral/MOM_MEKE.F90‎

Lines changed: 32 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -124,6 +124,7 @@ module MOM_MEKE
124124
logical :: debug !< If true, write out checksums of data for debugging
125125
integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default),
126126
!! read in from a file, or inferred via a neural network
127+
logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure.
127128
type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output
128129
!>@{ Diagnostic handles
129130
integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
@@ -400,6 +401,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
400401
call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m)
401402
endif
402403

404+
if (allocated(MEKE%Le)) then
405+
!$OMP parallel do default(shared)
406+
do j=js,je ; do i=is,ie
407+
MEKE%Le(i,j) = LmixScale(i,j)
408+
enddo ; enddo
409+
endif
410+
403411
! Aggregate sources of MEKE (background, frictional and GM)
404412
!$OMP parallel do default(shared)
405413
do j=js,je ; do i=is,ie
@@ -757,7 +765,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
757765
enddo ; enddo
758766
endif
759767

760-
if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) then
768+
if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
769+
.or. allocated(MEKE%Le)) then
761770
call cpu_clock_begin(CS%id_clock_pass)
762771
call do_group_pass(CS%pass_Kh, G%Domain)
763772
call cpu_clock_end(CS%id_clock_pass)
@@ -1425,6 +1434,9 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
14251434
"computing beta in the expression of Rhines scale. Use 1 if full "//&
14261435
"topographic beta effect is considered; use 0 if it's completely ignored.", &
14271436
units="nondim", default=0.0)
1437+
call get_param(param_file, mdl, "SQG_USE_MEKE", CS%sqg_use_MEKE, &
1438+
"If true, the eddy scale of MEKE is used for the SQG vertical structure ",&
1439+
default=.false.)
14281440

14291441
! Nonlocal module parameters
14301442
call get_param(param_file, mdl, "CDRAG", cdrag, &
@@ -1531,13 +1543,20 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
15311543

15321544
CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE)
15331545

1546+
15341547
! Detect whether this instance of MEKE_init() is at the beginning of a run
15351548
! or after a restart. If at the beginning, we will initialize MEKE to a local
15361549
! equilibrium.
15371550
CS%initialize = .not.query_initialized(MEKE%MEKE, "MEKE", restart_CS)
15381551
if (coldStart) CS%initialize = .false.
15391552
if (CS%initialize) call MOM_error(WARNING, &
15401553
"MEKE_init: Initializing MEKE with a local equilibrium balance.")
1554+
if (.not.query_initialized(MEKE%Le, "MEKE_Le", restart_CS) .and. allocated(MEKE%Le)) then
1555+
!$OMP parallel do default(shared)
1556+
do j=js,je ; do i=is,ie
1557+
MEKE%Le(i,j) = sqrt(G%areaT(i,j))
1558+
enddo ; enddo
1559+
endif
15411560

15421561
! Set up group passes. In the case of a restart, these fields need a halo update now.
15431562
if (allocated(MEKE%MEKE)) then
@@ -1548,8 +1567,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
15481567
if (allocated(MEKE%Kh)) call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain)
15491568
if (allocated(MEKE%Ku)) call create_group_pass(CS%pass_Kh, MEKE%Ku, G%Domain)
15501569
if (allocated(MEKE%Au)) call create_group_pass(CS%pass_Kh, MEKE%Au, G%Domain)
1570+
if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain)
15511571

1552-
if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) &
1572+
if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
1573+
.or. allocated(MEKE%Le)) &
15531574
call do_group_pass(CS%pass_Kh, G%Domain)
15541575

15551576
end function MEKE_init
@@ -1839,6 +1860,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
18391860
real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim]
18401861
logical :: Use_KH_in_MEKE
18411862
logical :: useMEKE
1863+
logical :: sqg_use_MEKE
18421864
integer :: isd, ied, jsd, jed
18431865

18441866
! Determine whether this module will be used
@@ -1853,6 +1875,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
18531875
MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku)
18541876
MEKE_viscCoeff_Au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au)
18551877
Use_KH_in_MEKE = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE)
1878+
sqg_use_MEKE = .false. ; call read_param(param_file,"SQG_USE_MEKE", sqg_use_MEKE)
18561879

18571880
if (.not. useMEKE) return
18581881

@@ -1884,6 +1907,12 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
18841907
longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", &
18851908
units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T)
18861909
endif
1910+
if (sqg_use_MEKE) then
1911+
allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0)
1912+
call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, &
1913+
longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", &
1914+
units="m", conversion=US%L_to_m)
1915+
endif
18871916
if (Use_Kh_in_MEKE) then
18881917
allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0)
18891918
call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, &
@@ -1918,6 +1947,7 @@ subroutine MEKE_end(MEKE)
19181947
if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh)
19191948
if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src)
19201949
if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE)
1950+
if (allocated(MEKE%Le)) deallocate(MEKE%Le)
19211951
end subroutine MEKE_end
19221952

19231953
!> \namespace mom_meke

‎src/parameterizations/lateral/MOM_MEKE_types.F90‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ module MOM_MEKE_types
2424
!! backscatter from unresolved eddies (see Jansen and Held, 2014).
2525
real, allocatable :: Au(:,:) !< The MEKE-derived lateral biharmonic viscosity
2626
!! coefficient [L4 T-1 ~> m4 s-1].
27+
real, allocatable :: Le(:,:) !< Eddy length scale [L m]
2728

2829
! Parameters
2930
real :: KhTh_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTh [nondim]

‎src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90‎

Lines changed: 117 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@ module MOM_lateral_mixing_coeffs
1818
use MOM_verticalGrid, only : verticalGrid_type
1919
use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init
2020
use MOM_open_boundary, only : ocean_OBC_type
21+
use MOM_MEKE_types, only : MEKE_type
22+
2123

2224
implicit none ; private
2325

@@ -112,9 +114,12 @@ module MOM_lateral_mixing_coeffs
112114

113115
real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [Z L-1 ~> nondim]
114116
real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [Z L-1 ~> nondim]
115-
real, allocatable :: ebt_struct(:,:,:) !< Vertical structure function to scale diffusivities with [nondim]
117+
real, allocatable :: ebt_struct(:,:,:) !< EBT vertical structure to scale diffusivities with [nondim]
118+
real, allocatable :: sqg_struct(:,:,:) !< SQG vertical structure to scale diffusivities with [nondim]
116119
real, allocatable :: BS_struct(:,:,:) !< Vertical structure function used in backscatter [nondim]
117120
real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0.
121+
real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 0.0
122+
logical :: BS_use_sqg !< If true, use sqg_stuct for backscatter vertical structure.
118123

119124

120125
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
@@ -163,6 +168,7 @@ module MOM_lateral_mixing_coeffs
163168
integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1
164169
integer :: id_dzu=-1, id_dzv=-1, id_dzSxN=-1, id_dzSyN=-1
165170
integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1
171+
integer :: id_sqg_struct=-1, id_BS_struct=-1
166172
type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
167173
!! timing of diagnostic output.
168174
!>@}
@@ -173,7 +179,7 @@ module MOM_lateral_mixing_coeffs
173179
end type VarMix_CS
174180

175181
public VarMix_init, VarMix_end, calc_slope_functions, calc_resoln_function
176-
public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function
182+
public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function, calc_sqg_struct
177183

178184
contains
179185

@@ -214,13 +220,14 @@ subroutine calc_depth_function(G, CS)
214220
end subroutine calc_depth_function
215221

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

225232
! Local variables
226233
! Depending on the power-function being used, dimensional rescaling may be limited, so some
@@ -230,6 +237,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
230237
real :: cg1_v ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1].
231238
real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2]
232239
integer :: power_2
240+
real :: dt !< Time increment [T ~> s]
233241
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
234242
integer :: i, j, k
235243
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
@@ -261,10 +269,19 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
261269
call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain)
262270
call do_group_pass(CS%pass_cg1, G%Domain)
263271
endif
272+
if (CS%sqg_expo>0.0) then
273+
call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
274+
call pass_var(CS%sqg_struct, G%Domain)
275+
endif
276+
264277
if (CS%BS_EBT_power>0.) then
265278
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
266279
CS%BS_struct(i,j,k) = CS%ebt_struct(i,j,k)**CS%BS_EBT_power
267280
enddo ; enddo ; enddo
281+
elseif (CS%BS_use_sqg) then
282+
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
283+
CS%BS_struct(i,j,k) = CS%sqg_struct(i,j,k)
284+
enddo ; enddo ; enddo
268285
endif
269286

270287
! Calculate and store the ratio between deformation radius and grid-spacing
@@ -460,6 +477,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
460477

461478
if (query_averaging_enabled(CS%diag)) then
462479
if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag)
480+
if (CS%id_BS_struct > 0) call post_data(CS%id_BS_struct, CS%BS_struct, CS%diag)
463481
endif
464482

465483
if (CS%debug) then
@@ -470,6 +488,77 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
470488

471489
end subroutine calc_resoln_function
472490

491+
!> Calculates and stores functions of SQG mode
492+
subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
493+
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
494+
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
495+
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
496+
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
497+
type(thermo_var_ptrs), intent(in) :: tv !<Thermodynamic variables
498+
real, intent(in) :: dt !< Time increment [T ~> s]
499+
type(VarMix_CS), intent(inout) :: CS !< Variable mixing control struct
500+
type(MEKE_type), intent(in) :: MEKE !< MEKE struct
501+
! type(ocean_OBC_type), pointer :: OBC !< Open
502+
! boundaries control structure.
503+
504+
! Local variables
505+
real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: &
506+
e ! The interface heights relative to mean sea level [Z ~> m].
507+
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]
508+
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]
509+
real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
510+
real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
511+
real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1]
512+
real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1]
513+
real, dimension(SZI_(G), SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1]
514+
real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2]
515+
real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m]
516+
! real, dimension(SZK_(GV)) :: zs ! Stretched vertical coordinate [m]
517+
real, dimension(SZI_(G), SZJ_(G)) :: Le ! Eddy length scale [m]
518+
integer :: i, j, k, is, ie, js, je, nz
519+
520+
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
521+
522+
if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
523+
"Module must be initialized before it is used.")
524+
525+
call find_eta(h, tv, G, GV, US, e, halo_size=2)
526+
call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
527+
CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v,dzu=dzu, dzv=dzv, &
528+
dzSxN=dzSxN, dzSyN=dzSyN, halo=1)
529+
530+
do j=js,je ; do i=is,ie
531+
CS%sqg_struct(i,j,1) = 1.0
532+
enddo ; enddo
533+
if (allocated(MEKE%Le)) then
534+
do j=js,je ; do i=is,ie
535+
Le(i,j) = MEKE%Le(i,j)
536+
f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + &
537+
G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8)
538+
enddo ; enddo
539+
else
540+
do j=js,je ; do i=is,ie
541+
Le(i,j) = sqrt(G%areaT(i,j))
542+
f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + &
543+
G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8)
544+
enddo ; enddo
545+
endif
546+
do k=2,nz ; do j=js,je ; do i=is,ie
547+
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)
548+
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)
549+
! dzs = -N2**0.5/f(i,j)*dzc
550+
CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1)*exp(-CS%sqg_expo*dzc/Le(i,j))
551+
enddo ; enddo ; enddo
552+
553+
554+
if (query_averaging_enabled(CS%diag)) then
555+
if (CS%id_sqg_struct > 0) call post_data(CS%id_sqg_struct, CS%sqg_struct, CS%diag)
556+
if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
557+
if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
558+
endif
559+
560+
end subroutine calc_sqg_struct
561+
473562
!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
474563
!! style scaling of diffusivity
475564
subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
@@ -1260,6 +1349,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
12601349
call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, &
12611350
"Power to raise EBT vertical structure to when backscatter "// &
12621351
"has vertical structure.", units="nondim", default=0.0)
1352+
call get_param(param_file, mdl, "BS_USE_SQG", CS%BS_use_sqg, &
1353+
"If true, the SQG vertical structure is used for backscatter "//&
1354+
"on the condition that BS_EBT_power=0", &
1355+
default=.false.)
1356+
if (CS%BS_EBT_power>0.) CS%BS_use_sqg = .false.
1357+
call get_param(param_file, mdl, "SQG_EXPO", CS%sqg_expo, &
1358+
"Nondimensional exponent coeffecient of the SQG mode "// &
1359+
"that is used for the vertical struture of diffusivities.", units="nondim", default=0.0)
1360+
if (CS%sqg_expo==0.) CS%BS_use_sqg = .false.
12631361
call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, &
12641362
"If true, uses the equivalent barotropic structure "//&
12651363
"as the vertical structure of thickness diffusivity.",&
@@ -1320,6 +1418,10 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
13201418
units="m", default=-1.0, scale=GV%m_to_H)
13211419
allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
13221420
endif
1421+
if (CS%sqg_expo>0.0) then
1422+
allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
1423+
endif
1424+
13231425
allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
13241426
CS%BS_struct(:,:,:) = 1.0
13251427

@@ -1335,7 +1437,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
13351437
endif
13361438
endif
13371439

1338-
if (CS%use_stored_slopes) then
1440+
if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then
13391441
! CS%calculate_Eady_growth_rate=.true.
13401442
in_use = .true.
13411443
allocate(CS%slope_x(IsdB:IedB,jsd:jed,GV%ke+1), source=0.0)
@@ -1418,7 +1520,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
14181520
'm2', conversion=US%L_to_m**2)
14191521
endif
14201522

1421-
if (CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) then
1523+
if (CS%sqg_expo>0.0) then
1524+
CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, &
1525+
'Vertical structure of SQG mode', 'nondim')
1526+
endif
1527+
1528+
CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, &
1529+
'Vertical structure of backscatter', 'nondim')
1530+
1531+
if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) .or. (CS%sqg_expo>0.0)) then
14221532
CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, &
14231533
'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
14241534
's-2', conversion=(US%L_to_Z*US%s_to_T)**2)
@@ -1672,9 +1782,10 @@ subroutine VarMix_end(CS)
16721782

16731783
if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.) &
16741784
deallocate(CS%ebt_struct)
1785+
if (CS%sqg_expo>0.0) deallocate(CS%sqg_struct)
16751786
if (allocated(CS%BS_struct)) deallocate(CS%BS_struct)
16761787

1677-
if (CS%use_stored_slopes) then
1788+
if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then
16781789
deallocate(CS%slope_x)
16791790
deallocate(CS%slope_y)
16801791
endif

0 commit comments

Comments
 (0)