Skip to content

Commit

Permalink
- Added if statement for CS%sqg_struct<=0.
Browse files Browse the repository at this point in the history 
- Changed the subroundoff of Coriolis parameter to 1e-40 in calc_sqg_struct
- Added indentation spaces for a few continuation lines
- Changed parameter BS_USE_SQG to BS_USE_SQG_STRUCT
  • Loading branch information
@Wendazhang33
Wendazhang33 committed Nov 4, 2024
1 parent 04e46cd commit 4b4e10b
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 35 deletions.
4 changes: 2 additions & 2 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -766,7 +766,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
endif

if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
.or. allocated(MEKE%Le)) then
.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 @@ -1570,7 +1570,7 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
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) &
.or. allocated(MEKE%Le)) &
.or. allocated(MEKE%Le)) &
call do_group_pass(CS%pass_Kh, G%Domain)

end function MEKE_init
Expand Down
72 changes: 39 additions & 33 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ module MOM_lateral_mixing_coeffs
real, allocatable :: khtr_struct(:,:,:) !< Vertical structure function used in tracer diffusivity [nondim]
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
real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 1.0
logical :: BS_use_sqg_struct !< If true, use sqg_stuct for backscatter vertical structure.


Expand Down Expand Up @@ -556,12 +556,14 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
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 :: 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 :: f_subround ! The minimal resolved value of Coriolis parameter to prevent division by zero [T-1 ~> s-1]
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
f_subround = 1.0e-40 * US%s_to_T

if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
"Module must be initialized before it is used.")
Expand All @@ -570,29 +572,33 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
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 (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%s_to_T)
enddo ; enddo

if (CS%sqg_expo<=0.) then
CS%sqg_struct(:,:,:) = 1.
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%s_to_T)
CS%sqg_struct(i,j,1) = 1.0
enddo ; enddo
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))), f_subround)
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))), f_subround)
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)))
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
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)))
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


if (query_averaging_enabled(CS%diag)) then
Expand Down Expand Up @@ -1393,7 +1399,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_struct, &
call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", 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 @@ -1480,26 +1486,26 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)

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")
"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")
"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")
"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")
"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
Expand Down Expand Up @@ -1616,7 +1622,7 @@ 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_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
.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

Expand Down Expand Up @@ -1890,7 +1896,7 @@ subroutine VarMix_end(CS)
type(VarMix_CS), intent(inout) :: CS

if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct &
.or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct)
.or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct)
if (allocated(CS%sqg_struct)) deallocate(CS%sqg_struct)
if (allocated(CS%BS_struct)) deallocate(CS%BS_struct)
if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) deallocate(CS%khth_struct)
Expand Down

0 comments on commit 4b4e10b

Please sign in to comment.