Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

Add SQG vertical structure to eddy diffusivities #738

Merged
merged 4 commits into from
Nov 13, 2024
Merged

Add SQG vertical structure to eddy diffusivities #738

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a9dd5f2
Added SQG vertical structure in Varmix to provide vertical profile fo…
Wendazhang33 Jul 27, 2024
8309221
Compute vertical structures for khth, khtr, backscatter, and kdgl90 a…
Wendazhang33 Oct 3, 2024
4ed3dca
Bug fixes
Wendazhang33 Oct 11, 2024
65328e6
Merge branch 'dev/gfdl' into add_sqg_struct
adcroft Nov 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -751,7 +751,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, dt)
call calc_depth_function(G, CS%VarMix)
call disable_averaging(CS%diag)
endif
Expand Down Expand Up @@ -1899,7 +1899,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, dt_offline)
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 @@ -1926,7 +1926,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, dt_offline)
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
34 changes: 32 additions & 2 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ module MOM_MEKE
logical :: debug !< If true, write out checksums of data for debugging
integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default),
!! read in from a file, or inferred via a neural network
logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure.
type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output
!>@{ Diagnostic handles
integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
Expand Down Expand Up @@ -400,6 +401,13 @@ 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

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

! 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 @@ -757,7 +765,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
enddo ; enddo
endif

if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) then
if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
.or. allocated(MEKE%Le)) then
call cpu_clock_begin(CS%id_clock_pass)
call do_group_pass(CS%pass_Kh, G%Domain)
call cpu_clock_end(CS%id_clock_pass)
Expand Down Expand Up @@ -1425,6 +1434,9 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
"computing beta in the expression of Rhines scale. Use 1 if full "//&
"topographic beta effect is considered; use 0 if it's completely ignored.", &
units="nondim", default=0.0)
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.)

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

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


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

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

if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) &
if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
.or. allocated(MEKE%Le)) &
call do_group_pass(CS%pass_Kh, G%Domain)

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

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

if (.not. useMEKE) return

Expand Down Expand Up @@ -1884,6 +1907,12 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", &
units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T)
endif
if (sqg_use_MEKE) then
allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0)
call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, &
longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", &
units="m", conversion=US%L_to_m)
endif
if (Use_Kh_in_MEKE) then
allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0)
call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, &
Expand Down Expand Up @@ -1918,6 +1947,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
Open in desktop
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
66 changes: 49 additions & 17 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -446,6 +446,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
logical :: use_MEKE_Ku
logical :: use_MEKE_Au
logical :: use_cont_huv
logical :: use_kh_struct
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
Expand Down Expand Up @@ -502,6 +503,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
if (CS%id_FrictWorkIntz > 0) find_FrictWork = .true.

if (allocated(MEKE%mom_src)) find_FrictWork = .true.
use_kh_struct = allocated(VarMix%BS_struct)
backscat_subround = 0.0
if (find_FrictWork .and. allocated(MEKE%mom_src) .and. (MEKE%backscatter_Ro_c > 0.0) .and. &
(MEKE%backscatter_Ro_Pow /= 0.0)) &
Expand Down Expand Up @@ -663,7 +665,7 @@ 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, use_kh_struct, &
!$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 Down Expand Up @@ -1181,14 +1183,26 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,

if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
! *Add* the MEKE contribution (which might be negative)
if (CS%res_scale_MEKE) then
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k)
enddo ; enddo
if (use_kh_struct) then
if (CS%res_scale_MEKE) then
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k)
enddo ; enddo
else
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
enddo ; enddo
endif
else
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
enddo ; enddo
if (CS%res_scale_MEKE) then
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j)
enddo ; enddo
else
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j)
enddo ; enddo
endif
endif
endif

Expand Down Expand Up @@ -1443,7 +1457,11 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
if (visc_limit_h_flag(i,j,k) > 0) then
Kh_BS(i,j) = 0.
else
Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
if (use_kh_struct) then
Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
else
Kh_BS(i,j) = MEKE%Ku(i,j)
endif
endif
enddo ; enddo

Expand Down Expand Up @@ -1618,10 +1636,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,

if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
! *Add* the MEKE contribution (might be negative)
Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
(MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
(MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn
if (use_kh_struct) then
Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
(MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
(MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn
else
Kh(I,J) = Kh(I,J) + 0.25*( (MEKE%Ku(i,j) + &
MEKE%Ku(i+1,j+1)) + &
(MEKE%Ku(i+1,j) + &
MEKE%Ku(i,j+1)) ) * meke_res_fn
endif
endif

if (CS%anisotropic) &
Expand Down Expand Up @@ -1789,10 +1814,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
if (visc_limit_q_flag(I,J,k) > 0) then
Kh_BS(I,J) = 0.
else
Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
(MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
(MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) )
if (use_kh_struct) then
Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
(MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
(MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) )
else
Kh_BS(I,J) = 0.25*( (MEKE%Ku(i,j) + &
MEKE%Ku(i+1,j+1)) + &
(MEKE%Ku(i+1,j) + &
MEKE%Ku(i,j+1)) )
endif
endif
enddo ; enddo

Expand Down
Loading