Skip to content

Commit

Permalink
(*)Corrected the units of VarMix_CS%slope_x
Browse files Browse the repository at this point in the history 
 VarMix_CS%slope_x was being set with units of [Z L-1 ~> nondim], but described
in comments as though it was simply [nondim], and then used in the (apparently
unused?) calculate_slopes=.false. branch in calc_slope_functions_using_just_e as
though its units actually were [nondim].  This commit corrects this
inconsistency, while also rescaling the internal slope variables in that routine
to also have the proper units of [Z L-1 ~> nondim].  In so doing, several
rescaling factors could be eliminated from the calculations.  In addition, the
slopes used in calc_QG_Leith_viscosity were also being rescaled with the wrong
factor or had dimensionally incorrect tiny values in some denominators, and this
has been corrected as well.  In testing this rescaling fix, a number of other
bugs were identified with USE_QG_LEITH_VISC=True (as described at
github.com/mom-ocean/issues/1590), so a fatal error message was added if
this option is enabled.  All answers in the MOM6-examples test suite are bitwise
identical, but the code will now give a fatal error if USE_QG_LEITH_VISC=.true.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Jan 31, 2023
1 parent b22617c commit 3348172
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 39 deletions.
8 changes: 7 additions & 1 deletion src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1936,8 +1936,14 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, &
"If true, use QG Leith nonlinear eddy viscosity.", &
default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) )
if (CS%use_QG_Leith_visc) then
call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//&
"should not be used until a number of bugs are fixed. Specifically it does not "//&
"reproduce across PE count or layout, and may use arrays that have not been properly "//&
"set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.")
endif
if (CS%use_QG_Leith_visc .and. .not. (CS%Leith_Kh .or. CS%Leith_Ah) ) then
call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//&
call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//&
"LEITH_KH or LEITH_AH must be True when USE_QG_LEITH_VISC=True.")
endif

Expand Down
90 changes: 52 additions & 38 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,8 +105,8 @@ module MOM_lateral_mixing_coeffs
!! spacing squared at v [L2 T-2 ~> m2 s-2].
real, allocatable :: Rd_dx_h(:,:) !< Deformation radius over grid spacing [nondim]

real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [nondim]
real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [nondim]
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 ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
Expand Down Expand Up @@ -142,7 +142,7 @@ module MOM_lateral_mixing_coeffs
integer :: Res_fn_power_visc !< The power of dx/Ld in the Kh resolution function. Any
!! positive integer power may be used, but even powers
!! and especially 2 are coded to be more efficient.
real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate [nondim].
real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate [Z L-1 ~> nondim].

! Leith parameters
logical :: use_QG_Leith_GM !< If true, uses the QG Leith viscosity as the GM coefficient
Expand Down Expand Up @@ -514,28 +514,33 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
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]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [nondim]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [Z L-1 ~> nondim]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency
!! at u-points [L2 Z-2 T-2 ~> s-2]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope [nondim]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope
!! [Z L-1 ~> nondim]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency
!! at v-points [L2 Z-2 T-2 ~> s-2]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure

! Local variables
real :: S2 ! Interface slope squared [nondim]
real :: N2 ! Positive buoyancy frequency or zero [T-2 ~> s-2]
real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
real :: N2 ! Positive buoyancy frequency or zero [L2 Z-2 T-2 ~> s-2]
real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
real :: H_geom ! The geometric mean of Hup and Hdn [H ~> m or kg m-2].
real :: S2max ! An upper bound on the squared slopes [nondim]
real :: S2max ! An upper bound on the squared slopes [Z2 L-2 ~> nondim]
real :: wNE, wSE, wSW, wNW ! Weights of adjacent points [nondim]
real :: H_u(SZIB_(G)), H_v(SZI_(G)) ! Layer thicknesses at u- and v-points [H ~> m or kg m-2]

! Note that at some points in the code S2_u and S2_v hold the running depth
! integrals of the squared slope [H ~> m or kg m-2] before the average is taken.
real :: S2_u(SZIB_(G),SZJ_(G)) ! The thickness-weighted depth average of the squared slope at u points [nondim].
real :: S2_v(SZI_(G),SZJB_(G)) ! The thickness-weighted depth average of the squared slope at v points [nondim].
real :: S2_u(SZIB_(G),SZJ_(G)) ! At first the thickness-weighted depth integral of the squared
! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
! squared slope [Z2 L-2 ~> nondim] at u points.
real :: S2_v(SZI_(G),SZJB_(G)) ! At first the thickness-weighted depth integral of the squared
! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
! squared slope [Z2 L-2 ~> nondim] at v points.

integer :: i, j, k, is, ie, js, je, nz, l_seg

Expand All @@ -545,7 +550,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
if (.not. CS%calculate_Eady_growth_rate) return
if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// &
"%SN_u is not associated with use_variable_mixing.")
if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:"// &
if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:R"// &
"%SN_v is not associated with use_variable_mixing.")

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Expand All @@ -563,7 +568,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
! and midlatitude deformation radii, using calc_resoln_function as a template.

!$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
do j = js,je
do j=js,je
do I=is-1,ie
CS%SN_u(I,j) = 0. ; H_u(I) = 0. ; S2_u(I,j) = 0.
enddo
Expand Down Expand Up @@ -599,7 +604,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
enddo

!$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
do J = js-1,je
do J=js-1,je
do i=is,ie
CS%SN_v(i,J) = 0. ; H_v(i) = 0. ; S2_v(i,J) = 0.
enddo
Expand Down Expand Up @@ -644,7 +649,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, &
scale=US%Z_to_L, haloshift=1)
call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, &
scale=US%L_to_Z**2 * US%s_to_T**2, scalar_pair=.true.)
scale=US%L_to_Z**2*US%s_to_T**2, scalar_pair=.true.)
call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
scale=US%s_to_T, scalar_pair=.true.)
endif
Expand Down Expand Up @@ -698,7 +703,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN,
enddo ; enddo

!$OMP parallel do default(shared) private(dnew,dz,weight,l_seg,vint_SN,sum_dz)
do j = G%jsc-1,G%jec+1
do j=G%jsc-1,G%jec+1
do I=G%isc-1,G%iec
vint_SN(I) = 0.
sum_dz(I) = dz_neglect
Expand Down Expand Up @@ -741,7 +746,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN,
enddo

!$OMP parallel do default(shared) private(dnew,dz,weight,l_seg)
do J = G%jsc-1,G%jec
do J=G%jsc-1,G%jec
do i=G%isc-1,G%iec+1
vint_SN(i) = 0.
sum_dz(i) = dz_neglect
Expand Down Expand Up @@ -782,14 +787,14 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN,
enddo
enddo

do j = G%jsc,G%jec
do j=G%jsc,G%jec
do I=G%isc-1,G%iec
CS%SN_u(I,j) = sqrt( SN_cpy(I,j)**2 &
+ 0.25*( (CS%SN_v(i,J)**2 + CS%SN_v(i+1,J-1)**2) &
+ (CS%SN_v(i+1,J)**2 + CS%SN_v(i,J-1)**2) ) )
enddo
enddo
do J = G%jsc-1,G%jec
do J=G%jsc-1,G%jec
do i=G%isc,G%iec
CS%SN_v(i,J) = sqrt( CS%SN_v(i,J)**2 &
+ 0.25*( (SN_cpy(I,j)**2 + SN_cpy(I-1,j+1)**2) &
Expand All @@ -816,13 +821,13 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
logical, intent(in) :: calculate_slopes !< If true, calculate slopes
!! internally otherwise use slopes stored in CS
! Local variables
real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics)
real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics)
real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [Z L-1 ~> nondim] (for diagnostics)
real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics)
real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: S2 ! Interface slope squared [nondim]
real :: N2 ! Brunt-Vaisala frequency squared [T-2 ~> s-2]
real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2]
real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times
Expand Down Expand Up @@ -857,16 +862,16 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
if (calculate_slopes) then
! Calculate the interface slopes E_x and E_y and u- and v- points respectively
do j=js-1,je+1 ; do I=is-1,ie
E_x(I,j) = US%Z_to_L*(e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)
E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)
! Mask slopes where interface intersects topography
if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) E_x(I,j) = 0.
enddo ; enddo
do J=js-1,je ; do i=is-1,ie+1
E_y(i,J) = US%Z_to_L*(e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)
E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)
! Mask slopes where interface intersects topography
if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0.
enddo ; enddo
else
else ! This branch is not used.
do j=js-1,je+1 ; do I=is-1,ie
E_x(I,j) = CS%slope_x(I,j,k)
if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) E_x(I,j) = 0.
Expand All @@ -880,22 +885,22 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
! Calculate N*S*h from this layer and add to the sum
do j=js,je ; do I=is-1,ie
S2 = ( E_x(I,j)**2 + 0.25*( &
(E_y(I,j)**2+E_y(I+1,j-1)**2)+(E_y(I+1,j)**2+E_y(I,j-1)**2) ) )
(E_y(I,j)**2+E_y(I+1,j-1)**2) + (E_y(I+1,j)**2+E_y(I,j-1)**2) ) )
Hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
Hup = 2.*h(i+1,j,k)*h(i+1,j,k-1) / (h(i+1,j,k) + h(i+1,j,k-1) + h_neglect)
H_geom = sqrt(Hdn*Hup)
N2 = GV%g_prime(k)*US%L_to_Z**2 / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < H_cutoff) &
S2 = 0.0
S2N2_u_local(I,j,k) = (H_geom * GV%H_to_Z) * S2 * N2
enddo ; enddo
do J=js-1,je ; do i=is,ie
S2 = ( E_y(i,J)**2 + 0.25*( &
(E_x(i,J)**2+E_x(i-1,J+1)**2)+(E_x(i,J+1)**2+E_x(i-1,J)**2) ) )
(E_x(i,J)**2+E_x(i-1,J+1)**2) + (E_x(i,J+1)**2+E_x(i-1,J)**2) ) )
Hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
Hup = 2.*h(i,j+1,k)*h(i,j+1,k-1) / (h(i,j+1,k) + h(i,j+1,k-1) + h_neglect)
H_geom = sqrt(Hdn*Hup)
N2 = GV%g_prime(k)*US%L_to_Z**2 / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < H_cutoff) &
S2 = 0.0
S2N2_v_local(i,J,k) = (H_geom * GV%H_to_Z) * S2 * N2
Expand Down Expand Up @@ -960,14 +965,14 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo
!! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
! Local variables
real, dimension(SZI_(G),SZJB_(G)) :: &
dslopey_dz, & ! z-derivative of y-slope at v-points [Z-1 ~> m-1]
dslopey_dz, & ! z-derivative of y-slope at v-points [L-1 ~> m-1]
h_at_v, & ! Thickness at v-points [H ~> m or kg m-2]
beta_v, & ! Beta at v-points [T-1 L-1 ~> s-1 m-1]
grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1]
grad_div_mag_v ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1]

real, dimension(SZIB_(G),SZJ_(G)) :: &
dslopex_dz, & ! z-derivative of x-slope at u-points [Z-1 ~> m-1]
dslopex_dz, & ! z-derivative of x-slope at u-points [L-1 ~> m-1]
h_at_u, & ! Thickness at u-points [H ~> m or kg m-2]
beta_u, & ! Beta at u-points [T-1 L-1 ~> s-1 m-1]
grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1]
Expand All @@ -987,41 +992,50 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo

if ((k > 1) .and. (k < nz)) then

! With USE_QG_LEITH_VISC=True, this might need to change to
! do j=js-2,je+2 ; do I=is-2,ie+1
! but other arrays used here (e.g., h and CS%slope_x) would also need to have wider valid halos.
do j=js-1,je+1 ; do I=is-2,Ieq+1
h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / &
( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
+ ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff )
+ ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff**2 )
h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / &
( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
+ ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff )
+ ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**2 )
Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_Z )
dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * Ih
h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
enddo ; enddo

! With USE_QG_LEITH_VISC=True, this might need to change to
! do J=js-2,je+1 ; do i=is-2,ie+2
do J=js-2,Jeq+1 ; do i=is-1,ie+1
h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / &
( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
+ ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff )
+ ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff**2 )
h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / &
( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
+ ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff )
+ ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**2 )
Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_Z )
dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * Ih
h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
enddo ; enddo

! With USE_QG_LEITH_VISC=True, this might need to be
! do J=js-2,je+1 ; do i=is-1,ie+1
do J=js-1,je ; do i=is-1,Ieq+1
f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) )
vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * US%L_to_Z * &
vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * &
( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) &
+ ( h_at_u(I-1,j) * dslopex_dz(I-1,j) + h_at_u(I,j+1) * dslopex_dz(I,j+1) ) ) / &
( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) + GV%H_subroundoff)
enddo ; enddo

! With USE_QG_LEITH_VISC=True, this might need to be
! do j=js-1,je+1 ; do I=is-2,ie+1
do j=js-1,Jeq+1 ; do I=is-1,ie
f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) )
vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * US%L_to_Z * &
vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * &
( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) &
+ ( h_at_v(i,J-1) * dslopey_dz(i,J-1) + h_at_v(i+1,J) * dslopey_dz(i+1,J) ) ) / &
( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) + GV%H_subroundoff)
Expand Down Expand Up @@ -1236,7 +1250,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"If non-zero, is an upper bound on slopes used in the "//&
"Visbeck formula for diffusivity. This does not affect the "//&
"isopycnal slope calculation used within thickness diffusion.", &
units="nondim", default=0.0)
units="nondim", default=0.0, scale=US%L_to_Z)
else
CS%Visbeck_S_max = 0.
endif
Expand Down

0 comments on commit 3348172

Please sign in to comment.