Skip to content

Commit

Permalink
Solved the dimension bug for sqg_struct
Browse files Browse the repository at this point in the history
  • Loading branch information
@Wendazhang33
Wendazhang33 committed Oct 17, 2024
1 parent 8309221 commit 04e46cd
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 63 deletions.
66 changes: 49 additions & 17 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
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
94 changes: 48 additions & 46 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ module MOM_lateral_mixing_coeffs
real, allocatable :: kdgl90_struct(:,:,:) !< Vertical structure function used in GL90 diffusivity [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 :: BS_use_sqg !< If true, use sqg_stuct for backscatter vertical structure.
logical :: BS_use_sqg_struct !< If true, use sqg_stuct for backscatter vertical structure.


real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
Expand Down Expand Up @@ -282,31 +282,23 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt)
call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain)
call do_group_pass(CS%pass_cg1, G%Domain)
endif
if (CS%BS_use_sqg .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
.or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>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. .and. CS%BS_use_sqg) then
call MOM_error(FATAL, &
"calc_resoln_function: BS_EBT_POWER>0. &
and BS_USE_SQG=True cannot be set together")
elseif (CS%BS_EBT_power>0.) then
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
elseif (CS%BS_use_sqg) then
elseif (CS%BS_use_sqg_struct) 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

if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT &
and KHTH_USE_SQG_STRUCT can be true")
elseif (CS%khth_use_ebt_struct) then
if (CS%khth_use_ebt_struct) then
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
CS%khth_struct(i,j,k) = CS%ebt_struct(i,j,k)
enddo ; enddo ; enddo
Expand All @@ -316,11 +308,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt)
enddo ; enddo ; enddo
endif

if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT &
and KHTR_USE_SQG_STRUCT can be true")
elseif (CS%khtr_use_ebt_struct) then
if (CS%khtr_use_ebt_struct) then
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
CS%khtr_struct(i,j,k) = CS%ebt_struct(i,j,k)
enddo ; enddo ; enddo
Expand All @@ -330,11 +318,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt)
enddo ; enddo ; enddo
endif

if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT &
and KD_GL90_USE_SQG_STRUCT can be true")
elseif (CS%kdgl90_use_ebt_struct) then
if (CS%kdgl90_use_ebt_struct) then
do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
CS%kdgl90_struct(i,j,k) = CS%ebt_struct(i,j,k)
enddo ; enddo ; enddo
Expand Down Expand Up @@ -561,8 +545,6 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
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) :: &
Expand All @@ -576,7 +558,6 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
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

Expand All @@ -596,27 +577,21 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
if (allocated(MEKE%Le)) 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/US%T_to_s)
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 * US%s_to_T)
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*US%T_to_s)
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 * US%s_to_T)
enddo ; enddo
endif
if (CS%debug) then
call hchksum(Le, 'SQG length scale', G%HI, unscale=US%L_to_m)
call hchksum(f, 'Coriolis at h point', G%HI, unscale=US%s_to_T)
call uvchksum( 'MEKE LmixScale', dzu, dzv, G%HI, unscale=US%Z_to_m, scalar_pair=.true.)
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)*US%Z_to_L
! 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))
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)))
CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j))))
enddo ; enddo ; enddo


Expand Down Expand Up @@ -1418,7 +1393,7 @@ 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, &
call get_param(param_file, mdl, "BS_USE_SQG", CS%BS_use_sqg_struct, &
"If true, the SQG vertical structure is used for backscatter "//&
"on the condition that BS_EBT_power=0", &
default=.false.)
Expand Down Expand Up @@ -1503,8 +1478,33 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
endif


allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
CS%BS_struct(:,:,:) = 1.0
if (CS%BS_EBT_power>0. .and. CS%BS_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: BS_EBT_POWER>0. &
and BS_USE_SQG=True cannot be set together")
endif

if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT &
and KHTH_USE_SQG_STRUCT can be true")
endif

if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT &
and KHTR_USE_SQG_STRUCT can be true")
endif

if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) then
call MOM_error(FATAL, &
"calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT &
and KD_GL90_USE_SQG_STRUCT can be true")
endif

if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then
allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
endif

if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then
allocate(CS%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
Expand Down Expand Up @@ -1615,13 +1615,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)

CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, &
'Vertical structure of SQG mode', 'nondim')
if (CS%BS_use_sqg .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
.or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then
allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
endif

CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, &
'Vertical structure of backscatter', 'nondim')
if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then
CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, &
'Vertical structure of backscatter', 'nondim')
endif
if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then
CS%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, Time, &
'Vertical structure of thickness diffusivity', 'nondim')
Expand Down

0 comments on commit 04e46cd

Please sign in to comment.