From efaf90336b852d11b04a7a59861e0e1de237b617 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 10:47:25 -0400 Subject: [PATCH 1/8] *Land masking isopycnal slopes - These updates apply masking of C-grid slopes and Brunt-Vaisala frequency in order to avoid dependency on land values --- src/core/MOM_isopycnal_slopes.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 9f716c0297..bd6848a342 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -203,15 +203,15 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & ! to be between -1 and 1. mag_grad2 = drdx**2 + drdz**2 if (mag_grad2 > 0.0) then - slope_x(I,j,K) = drdx / sqrt(mag_grad2) + slope_x(I,j,K) = drdx / sqrt(mag_grad2) * G%mask2dCu(I,j) else ! Just in case mag_grad2 = 0 ever. slope_x(I,j,K) = 0.0 endif - - if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz ! Square of Brunt-Vaisala frequency (s-2) + + if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency (s-2) else ! With .not.use_EOS, the layers are constant density. - slope_x(I,j,K) = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * GV%m_to_H + slope_x(I,j,K) = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * GV%m_to_H * G%mask2dCu(I,j) endif enddo ! I @@ -287,15 +287,15 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & ! to be between -1 and 1. mag_grad2 = drdy**2 + drdz**2 if (mag_grad2 > 0.0) then - slope_y(i,J,K) = drdy / sqrt(mag_grad2) + slope_y(i,J,K) = drdy / sqrt(mag_grad2) * G%mask2dCv(i,J) else ! Just in case mag_grad2 = 0 ever. slope_y(i,J,K) = 0.0 endif - if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz ! Square of Brunt-Vaisala frequency (s-2) + if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of Brunt-Vaisala frequency (s-2) else ! With .not.use_EOS, the layers are constant density. - slope_y(i,J,K) = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * GV%m_to_H + slope_y(i,J,K) = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * GV%m_to_H * G%mask2dCv(i,J) endif enddo ! i From 8bdf928073665212178b2246ace12f5665722acd Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 11:09:32 -0400 Subject: [PATCH 2/8] *Apply land mask to MEKE mass --- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index bcc56ab381..b52c0030b6 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -238,7 +238,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv) do j=js-1,je+1 do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo do k=1,nz ; do i=is-1,ie+1 - mass(i,j) = mass(i,j) + GV%H_to_kg_m2 * h(i,j,k) + mass(i,j) = mass(i,j) + GV%H_to_kg_m2 * h(i,j,k) * G%mask2dT(i,j) enddo ; enddo do i=is-1,ie+1 I_mass(i,j) = 0.0 From 6078cca773e15edd1ea2aa6ede76d2f31614bb22 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 11:24:18 -0400 Subject: [PATCH 3/8] *Additional masking for SN_uv - avoids dependency on land values in the calculation of lateral mixing coefficients. --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 42 ++++++++++++------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 94872fc4ef..ca27bbe6b4 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -3,6 +3,7 @@ module MOM_lateral_mixing_coeffs ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_debugging, only : hchksum, uvchksum use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, post_data use MOM_diag_mediator, only : diag_ctrl, time_type, query_averaging_enabled @@ -106,7 +107,7 @@ module MOM_lateral_mixing_coeffs type(wave_speed_CS), pointer :: wave_speed_CSp => NULL() !< Wave speed control structure type(group_pass_type) :: pass_cg1 !< For group halo pass - + logical :: debug !< If true, write out checksums of data for debugging end type VarMix_CS public VarMix_init, calc_slope_functions, calc_resoln_function @@ -437,12 +438,12 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) !$OMP private(E_x,E_y,S2,H_u,H_v,Hdn,Hup,H_geom,N2, & !$OMP wNE, wSE, wSW, wNW) !$OMP do + do j=js-1,je+1 ; do i=is-1,ie+1 CS%SN_u(i,j) = 0.0 CS%SN_v(i,j) = 0.0 enddo ; enddo - ! To set the length scale based on the deformation radius, use wave_speed to ! calculate the first-mode gravity wave speed and then blend the equatorial ! and midlatitude deformation radii, using calc_resoln_function as a template. @@ -469,10 +470,10 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - wSE = h(i+1,j,k)*h(i+1,j-1,k) * h(i+1,j,k)*h(i+1,j-1,k-1) - wNW = h(i ,j,k)*h(i ,j+1,k) * h(i ,j,k)*h(i ,j+1,k-1) - wNE = h(i+1,j,k)*h(i+1,j+1,k) * h(i+1,j,k)*h(i+1,j+1,k-1) - wSW = h(i ,j,k)*h(i ,j-1,k) * h(i ,j,k)*h(i ,j-1,k-1) + wSE = h(i+1,j,k)*h(i+1,j-1,k) * h(i+1,j,k)*h(i+1,j-1,k-1) * G%mask2dCv(i+1,J-1) + wNW = h(i ,j,k)*h(i ,j+1,k) * h(i ,j,k)*h(i ,j+1,k-1) * G%mask2dCv(i,J) + wNE = h(i+1,j,k)*h(i+1,j+1,k) * h(i+1,j,k)*h(i+1,j+1,k-1) * G%mask2dCv(i+1,J) + wSW = h(i ,j,k)*h(i ,j-1,k) * h(i ,j,k)*h(i ,j-1,k-1) * G%mask2dCv(i,J-1) S2 = slope_x(I,j,K)**2 + ( & (wNW*slope_y(i,J,K)**2+wSE*slope_y(i+1,J-1,K)**2) & +(wNE*slope_y(i+1,J,K)**2+wSW*slope_y(i,J-1,K)**2) ) / & @@ -485,8 +486,8 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) enddo ; enddo do I=is-1,ie if (H_u(I)>0.) then - CS%SN_u(I,j) = CS%SN_u(I,j) / H_u(I) - S2_u(I,j) = S2_u(I,j) / H_u(I) + CS%SN_u(I,j) = CS%SN_u(I,j) / H_u(I) * G%mask2dCu(I,j) + S2_u(I,j) = S2_u(I,j) / H_u(I) * G%mask2dCu(I,j) else CS%SN_u(I,j) = 0. endif @@ -504,10 +505,10 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - wSE = h(i,j ,k)*h(i+1,j ,k) * h(i,j ,k)*h(i+1,j ,k-1) - wNW = h(i,j+1,k)*h(i-1,j+1,k) * h(i,j+1,k)*h(i-1,j+1,k-1) - wNE = h(i,j+1,k)*h(i+1,j+1,k) * h(i,j+1,k)*h(i+1,j+1,k-1) - wSW = h(i,j ,k)*h(i-1,j ,k) * h(i,j ,k)*h(i-1,j ,k-1) + wSE = h(i,j ,k)*h(i+1,j ,k) * h(i,j ,k)*h(i+1,j ,k-1) * G%mask2dCu(I,j) + wNW = h(i,j+1,k)*h(i-1,j+1,k) * h(i,j+1,k)*h(i-1,j+1,k-1) * G%mask2dCu(I-1,j+1) + wNE = h(i,j+1,k)*h(i+1,j+1,k) * h(i,j+1,k)*h(i+1,j+1,k-1) * G%mask2dCu(I,j+1) + wSW = h(i,j ,k)*h(i-1,j ,k) * h(i,j ,k)*h(i-1,j ,k-1) * G%mask2dCu(I-1,j) S2 = slope_y(i,J,K)**2 + ( & (wSE*slope_x(I,j,K)**2+wNW*slope_x(I-1,j+1,K)**2) & +(wNE*slope_x(I,j+1,K)**2+wSW*slope_x(I-1,j,K)**2) ) / & @@ -520,8 +521,8 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) enddo ; enddo do i=is,ie if (H_v(i)>0.) then - CS%SN_v(i,J) = CS%SN_v(i,J) / H_v(i) - S2_v(i,J) = S2_v(i,J) / H_v(i) + CS%SN_v(i,J) = CS%SN_v(i,J) / H_v(i) * G%mask2dCv(i,J) + S2_v(i,J) = S2_v(i,J) / H_v(i) * G%mask2dCv(i,J) else CS%SN_v(i,J) = 0. endif @@ -536,6 +537,12 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) if (CS%id_S2_v > 0) call post_data(CS%id_S2_v, S2_v, CS%diag) endif + if (CS%debug) then + call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) + call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", N2_u, N2_v, G%HI) + call uvchksum("calc_Visbeck_coeffs SN_[uv]", CS%SN_u, CS%SN_v, G%HI) + endif + end subroutine calc_Visbeck_coeffs !> The original calc_slope_function() that calculated slopes using @@ -663,7 +670,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + GV%Angstrom ) ) !The code below behaves better than the line above. Not sure why? AJA if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then - CS%SN_u(I,j) = sqrt( CS%SN_u(I,j) / max(G%bathyT(I,j), G%bathyT(I+1,j)) ) + CS%SN_u(I,j) = sqrt( CS%SN_u(I,j) / max(G%bathyT(I,j), G%bathyT(I+1,j)) ) * G%mask2dCu(I,j) else CS%SN_u(I,j) = 0.0 endif @@ -678,7 +685,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + GV%Angstrom ) ) !The code below behaves better than the line above. Not sure why? AJA if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then - CS%SN_v(i,J) = sqrt( CS%SN_v(i,J) / max(G%bathyT(i,J), G%bathyT(i,J+1)) ) + CS%SN_v(i,J) = sqrt( CS%SN_v(i,J) / max(G%bathyT(i,J), G%bathyT(i,J+1)) ) * G%mask2dCv(i,J) else CS%SN_v(I,j) = 0.0 endif @@ -762,6 +769,8 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "stored for re-use. This uses more memory but avoids calling\n"//& "the equation of state more times than should be necessary.", & default=.false.) + + if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) use_variable_mixing = .true. if (use_variable_mixing .or. Resoln_scaled_Kh .or. Resoln_scaled_KhTh .or. & @@ -775,6 +784,7 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) CS%khth_use_ebt_struct = khth_use_ebt_struct CS%use_variable_mixing = use_variable_mixing CS%use_stored_slopes = use_stored_slopes + call get_param(param_file, mod, "DEBUG", CS%debug, default=.false., do_not_log=.true.) else return endif From 9d444cecb263813aeab65122cd741070605df00b Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 14:58:45 -0400 Subject: [PATCH 4/8] remove trailing spaces --- src/core/MOM_isopycnal_slopes.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index bd6848a342..68ad0973f6 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -207,11 +207,11 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & else ! Just in case mag_grad2 = 0 ever. slope_x(I,j,K) = 0.0 endif - + if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency (s-2) else ! With .not.use_EOS, the layers are constant density. - slope_x(I,j,K) = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * GV%m_to_H * G%mask2dCu(I,j) + slope_x(I,j,K) = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * GV%m_to_H * G%mask2dCu(I,j) endif enddo ! I From d2831d55f2363edeffc2c4e25177e55c610487b4 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 16:38:18 -0400 Subject: [PATCH 5/8] Revert slope masking at coastlines. - This does not change answers for MEKE relative to the previous commit. --- src/core/MOM_isopycnal_slopes.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 68ad0973f6..24d8fe85a9 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -203,7 +203,7 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & ! to be between -1 and 1. mag_grad2 = drdx**2 + drdz**2 if (mag_grad2 > 0.0) then - slope_x(I,j,K) = drdx / sqrt(mag_grad2) * G%mask2dCu(I,j) + slope_x(I,j,K) = drdx / sqrt(mag_grad2) else ! Just in case mag_grad2 = 0 ever. slope_x(I,j,K) = 0.0 endif @@ -211,7 +211,7 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency (s-2) else ! With .not.use_EOS, the layers are constant density. - slope_x(I,j,K) = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * GV%m_to_H * G%mask2dCu(I,j) + slope_x(I,j,K) = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * GV%m_to_H endif enddo ! I @@ -287,7 +287,7 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & ! to be between -1 and 1. mag_grad2 = drdy**2 + drdz**2 if (mag_grad2 > 0.0) then - slope_y(i,J,K) = drdy / sqrt(mag_grad2) * G%mask2dCv(i,J) + slope_y(i,J,K) = drdy / sqrt(mag_grad2) else ! Just in case mag_grad2 = 0 ever. slope_y(i,J,K) = 0.0 endif @@ -295,7 +295,7 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of Brunt-Vaisala frequency (s-2) else ! With .not.use_EOS, the layers are constant density. - slope_y(i,J,K) = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * GV%m_to_H * G%mask2dCv(i,J) + slope_y(i,J,K) = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * GV%m_to_H endif enddo ! i From f07a2ea9734c95acedb26648590a687eabd6bcc9 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 17:11:15 -0400 Subject: [PATCH 6/8] Add parentheses and re-order multiplication --- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index b52c0030b6..354f25f5dd 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -238,7 +238,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv) do j=js-1,je+1 do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo do k=1,nz ; do i=is-1,ie+1 - mass(i,j) = mass(i,j) + GV%H_to_kg_m2 * h(i,j,k) * G%mask2dT(i,j) + mass(i,j) = mass(i,j) + G%mask2dT(i,j) * (GV%H_to_kg_m2 * h(i,j,k)) enddo ; enddo do i=is-1,ie+1 I_mass(i,j) = 0.0 From 943ab1fccf62749bf8c65db773321b6b94c9d23d Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 17:16:39 -0400 Subject: [PATCH 7/8] *New parameter use_Visbeck_slope_bug This is False be default and is retained for legacy purposes. Setting this to True fixes issues with erroneous k-indexing in the calulation of harmonic mean thicknesses which are applied to neutral slopes in order to compute SN_uv. This will impact answers and needs to be evaluated in existing runs with MEKE turned on. --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 80 +++++++++++++------ 1 file changed, 55 insertions(+), 25 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index ca27bbe6b4..80fbbc8504 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -40,6 +40,11 @@ module MOM_lateral_mixing_coeffs !! of first baroclinic wave for calculating the resolution fn. logical :: khth_use_ebt_struct !< If true, uses the equivalent barotropic structure !! as the vertical structure of thickness diffusivity. + logical :: use_Visbeck_slope_bug !< If true, then retain a legacy bug in the calculation of weights + !! applied to isoneutral slopes. There was an erroneous k-indexing + !! for layer thicknesses. In addition, masking at coastlines was not + !! used which introduced potential restart issues. This flag will be + !! deprecated in a future release. real, dimension(:,:), pointer :: & SN_u => NULL(), & !< S*N at u-points (s^-1) SN_v => NULL(), & !< S*N at v-points (s^-1) @@ -438,7 +443,6 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) !$OMP private(E_x,E_y,S2,H_u,H_v,Hdn,Hup,H_geom,N2, & !$OMP wNE, wSE, wSW, wNW) !$OMP do - do j=js-1,je+1 ; do i=is-1,ie+1 CS%SN_u(i,j) = 0.0 CS%SN_v(i,j) = 0.0 @@ -470,14 +474,25 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - wSE = h(i+1,j,k)*h(i+1,j-1,k) * h(i+1,j,k)*h(i+1,j-1,k-1) * G%mask2dCv(i+1,J-1) - wNW = h(i ,j,k)*h(i ,j+1,k) * h(i ,j,k)*h(i ,j+1,k-1) * G%mask2dCv(i,J) - wNE = h(i+1,j,k)*h(i+1,j+1,k) * h(i+1,j,k)*h(i+1,j+1,k-1) * G%mask2dCv(i+1,J) - wSW = h(i ,j,k)*h(i ,j-1,k) * h(i ,j,k)*h(i ,j-1,k-1) * G%mask2dCv(i,J-1) - S2 = slope_x(I,j,K)**2 + ( & - (wNW*slope_y(i,J,K)**2+wSE*slope_y(i+1,J-1,K)**2) & - +(wNE*slope_y(i+1,J,K)**2+wSW*slope_y(i,J-1,K)**2) ) / & - ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**2 ) !### This should be **4 for consistent units. + if (CS%use_Visbeck_slope_bug) then + wSE = h(i+1,j,k)*h(i+1,j-1,k) * h(i+1,j,k)*h(i+1,j-1,k-1) + wNW = h(i ,j,k)*h(i ,j+1,k) * h(i ,j,k)*h(i ,j+1,k-1) + wNE = h(i+1,j,k)*h(i+1,j+1,k) * h(i+1,j,k)*h(i+1,j+1,k-1) + wSW = h(i ,j,k)*h(i ,j-1,k) * h(i ,j,k)*h(i ,j-1,k-1) + S2 = slope_x(I,j,K)**2 + ( & + (wNW*slope_y(i,J,K)**2+wSE*slope_y(i+1,J-1,K)**2) & + +(wNE*slope_y(i+1,J,K)**2+wSW*slope_y(i,J-1,K)**2) ) / & + ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**2 ) !### This should be **4 for consistent units. + else + wSE = G%mask2dCv(i+1,J-1) * ( (h(i+1,j,k)*h(i+1,j-1,k)) * (h(i+1,j,k-1)*h(i+1,j-1,k-1)) ) + wNW = G%mask2dCv(i ,J ) * ( (h(i ,j,k)*h(i ,j+1,k)) * (h(i ,j,k-1)*h(i ,j+1,k-1)) ) + wNE = G%mask2dCv(i+1,J ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) ) + wSW = G%mask2dCv(i ,J-1) * ( (h(i ,j,k)*h(i ,j-1,k)) * (h(i ,j,k-1)*h(i ,j-1,k-1)) ) + S2 = slope_x(I,j,K)**2 + ( & + (wNW*slope_y(i,J,K)**2+wSE*slope_y(i+1,J-1,K)**2) & + +(wNE*slope_y(i+1,J,K)**2+wSW*slope_y(i,J-1,K)**2) ) / & + ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**4 ) + endif if (S2max>0.) S2 = S2 * S2max / (S2 + S2max) ! Limit S2 N2 = max(0., N2_u(I,j,k)) CS%SN_u(I,j) = CS%SN_u(I,j) + sqrt( S2*N2 )*H_geom @@ -486,8 +501,8 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) enddo ; enddo do I=is-1,ie if (H_u(I)>0.) then - CS%SN_u(I,j) = CS%SN_u(I,j) / H_u(I) * G%mask2dCu(I,j) - S2_u(I,j) = S2_u(I,j) / H_u(I) * G%mask2dCu(I,j) + CS%SN_u(I,j) = G%mask2dCu(I,j) * CS%SN_u(I,j) / H_u(I) + S2_u(I,j) = G%mask2dCu(I,j) * S2_u(I,j) / H_u(I) else CS%SN_u(I,j) = 0. endif @@ -505,14 +520,25 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - wSE = h(i,j ,k)*h(i+1,j ,k) * h(i,j ,k)*h(i+1,j ,k-1) * G%mask2dCu(I,j) - wNW = h(i,j+1,k)*h(i-1,j+1,k) * h(i,j+1,k)*h(i-1,j+1,k-1) * G%mask2dCu(I-1,j+1) - wNE = h(i,j+1,k)*h(i+1,j+1,k) * h(i,j+1,k)*h(i+1,j+1,k-1) * G%mask2dCu(I,j+1) - wSW = h(i,j ,k)*h(i-1,j ,k) * h(i,j ,k)*h(i-1,j ,k-1) * G%mask2dCu(I-1,j) - S2 = slope_y(i,J,K)**2 + ( & - (wSE*slope_x(I,j,K)**2+wNW*slope_x(I-1,j+1,K)**2) & - +(wNE*slope_x(I,j+1,K)**2+wSW*slope_x(I-1,j,K)**2) ) / & - ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**2 ) !### This should be **4 for consistent units. + if (CS%use_Visbeck_slope_bug) then + wSE = h(i,j ,k)*h(i+1,j ,k) * h(i,j ,k)*h(i+1,j ,k-1) + wNW = h(i,j+1,k)*h(i-1,j+1,k) * h(i,j+1,k)*h(i-1,j+1,k-1) + wNE = h(i,j+1,k)*h(i+1,j+1,k) * h(i,j+1,k)*h(i+1,j+1,k-1) + wSW = h(i,j ,k)*h(i-1,j ,k) * h(i,j ,k)*h(i-1,j ,k-1) + S2 = slope_y(i,J,K)**2 + ( & + (wSE*slope_x(I,j,K)**2+wNW*slope_x(I-1,j+1,K)**2) & + +(wNE*slope_x(I,j+1,K)**2+wSW*slope_x(I-1,j,K)**2) ) / & + ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**2 ) !### This should be **4 for consistent units. + else + wSE = G%mask2dCu(I,j) * ( (h(i,j ,k)*h(i+1,j ,k)) * (h(i,j ,k-1)*h(i+1,j ,k-1)) ) + wNW = G%mask2dCu(I-1,j+1) * ( (h(i,j+1,k)*h(i-1,j+1,k)) * (h(i,j+1,k-1)*h(i-1,j+1,k-1)) ) + wNE = G%mask2dCu(I,j+1) * ( (h(i,j+1,k)*h(i+1,j+1,k)) * (h(i,j+1,k-1)*h(i+1,j+1,k-1)) ) + wSW = G%mask2dCu(I-1,j) * ( (h(i,j ,k)*h(i-1,j ,k)) * (h(i,j ,k-1)*h(i-1,j ,k-1)) ) + S2 = slope_y(i,J,K)**2 + ( & + (wSE*slope_x(I,j,K)**2+wNW*slope_x(I-1,j+1,K)**2) & + +(wNE*slope_x(I,j+1,K)**2+wSW*slope_x(I-1,j,K)**2) ) / & + ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**4 ) !### This should be **4 for consistent units. + endif if (S2max>0.) S2 = S2 * S2max / (S2 + S2max) ! Limit S2 N2 = max(0., N2_v(i,J,K)) CS%SN_v(i,J) = CS%SN_v(i,J) + sqrt( S2*N2 )*H_geom @@ -521,8 +547,8 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) enddo ; enddo do i=is,ie if (H_v(i)>0.) then - CS%SN_v(i,J) = CS%SN_v(i,J) / H_v(i) * G%mask2dCv(i,J) - S2_v(i,J) = S2_v(i,J) / H_v(i) * G%mask2dCv(i,J) + CS%SN_v(i,J) = G%mask2dCv(i,J) * CS%SN_v(i,J) / H_v(i) + S2_v(i,J) = G%mask2dCv(i,J) * S2_v(i,J) / H_v(i) else CS%SN_v(i,J) = 0. endif @@ -670,7 +696,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + GV%Angstrom ) ) !The code below behaves better than the line above. Not sure why? AJA if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then - CS%SN_u(I,j) = sqrt( CS%SN_u(I,j) / max(G%bathyT(I,j), G%bathyT(I+1,j)) ) * G%mask2dCu(I,j) + CS%SN_u(I,j) = G%mask2dCu(I,j) * sqrt( CS%SN_u(I,j) / max(G%bathyT(I,j), G%bathyT(I+1,j)) ) else CS%SN_u(I,j) = 0.0 endif @@ -685,7 +711,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + GV%Angstrom ) ) !The code below behaves better than the line above. Not sure why? AJA if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff ) then - CS%SN_v(i,J) = sqrt( CS%SN_v(i,J) / max(G%bathyT(i,J), G%bathyT(i,J+1)) ) * G%mask2dCv(i,J) + CS%SN_v(i,J) = G%mask2dCv(i,J) * sqrt( CS%SN_v(i,J) / max(G%bathyT(i,J), G%bathyT(i,J+1)) ) else CS%SN_v(I,j) = 0.0 endif @@ -769,8 +795,6 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "stored for re-use. This uses more memory but avoids calling\n"//& "the equation of state more times than should be necessary.", & default=.false.) - - if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) use_variable_mixing = .true. if (use_variable_mixing .or. Resoln_scaled_Kh .or. Resoln_scaled_KhTh .or. & @@ -906,6 +930,12 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "velocity points from the thickness points; otherwise \n"//& "interpolate the wave speed and calculate the resolution \n"//& "function independently at each point.", default=.true.) + call get_param(param_file, mod, "USE_VISBECK_SLOPE_BUG", CS%use_Visbeck_slope_bug, & + "If true, then retain a legacy bug in the calculation of weights \n"//& + "applied to isoneutral slopes. There was an erroneous k-indexing \n"//& + "for layer thicknesses. In addition, masking at coastlines was not \n"//& + "used which introduced potential restart issues. This flag will be \m"//& + "deprecated in a future release.", default=.false.) if (CS%interpolate_Res_fn) then if (CS%Res_coef_visc .ne. CS%Res_coef_khth) call MOM_error(FATAL, & "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//& From 16c68f159b25e5472ba146eabbf090d6138873f6 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Fri, 2 Jun 2017 17:24:01 -0400 Subject: [PATCH 8/8] remove trailing whitespace --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 80fbbc8504..6635d5f550 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -474,7 +474,7 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - if (CS%use_Visbeck_slope_bug) then + if (CS%use_Visbeck_slope_bug) then wSE = h(i+1,j,k)*h(i+1,j-1,k) * h(i+1,j,k)*h(i+1,j-1,k-1) wNW = h(i ,j,k)*h(i ,j+1,k) * h(i ,j,k)*h(i ,j+1,k-1) wNE = h(i+1,j,k)*h(i+1,j+1,k) * h(i+1,j,k)*h(i+1,j+1,k-1) @@ -486,8 +486,8 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) else wSE = G%mask2dCv(i+1,J-1) * ( (h(i+1,j,k)*h(i+1,j-1,k)) * (h(i+1,j,k-1)*h(i+1,j-1,k-1)) ) wNW = G%mask2dCv(i ,J ) * ( (h(i ,j,k)*h(i ,j+1,k)) * (h(i ,j,k-1)*h(i ,j+1,k-1)) ) - wNE = G%mask2dCv(i+1,J ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) ) - wSW = G%mask2dCv(i ,J-1) * ( (h(i ,j,k)*h(i ,j-1,k)) * (h(i ,j,k-1)*h(i ,j-1,k-1)) ) + wNE = G%mask2dCv(i+1,J ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) ) + wSW = G%mask2dCv(i ,J-1) * ( (h(i ,j,k)*h(i ,j-1,k)) * (h(i ,j,k-1)*h(i ,j-1,k-1)) ) S2 = slope_x(I,j,K)**2 + ( & (wNW*slope_y(i,J,K)**2+wSE*slope_y(i+1,J-1,K)**2) & +(wNE*slope_y(i+1,J,K)**2+wSW*slope_y(i,J-1,K)**2) ) / & @@ -520,7 +520,7 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - if (CS%use_Visbeck_slope_bug) then + if (CS%use_Visbeck_slope_bug) then wSE = h(i,j ,k)*h(i+1,j ,k) * h(i,j ,k)*h(i+1,j ,k-1) wNW = h(i,j+1,k)*h(i-1,j+1,k) * h(i,j+1,k)*h(i-1,j+1,k-1) wNE = h(i,j+1,k)*h(i+1,j+1,k) * h(i,j+1,k)*h(i+1,j+1,k-1) @@ -530,7 +530,7 @@ subroutine calc_Visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS) +(wNE*slope_x(I,j+1,K)**2+wSW*slope_x(I-1,j,K)**2) ) / & ( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**2 ) !### This should be **4 for consistent units. else - wSE = G%mask2dCu(I,j) * ( (h(i,j ,k)*h(i+1,j ,k)) * (h(i,j ,k-1)*h(i+1,j ,k-1)) ) + wSE = G%mask2dCu(I,j) * ( (h(i,j ,k)*h(i+1,j ,k)) * (h(i,j ,k-1)*h(i+1,j ,k-1)) ) wNW = G%mask2dCu(I-1,j+1) * ( (h(i,j+1,k)*h(i-1,j+1,k)) * (h(i,j+1,k-1)*h(i-1,j+1,k-1)) ) wNE = G%mask2dCu(I,j+1) * ( (h(i,j+1,k)*h(i+1,j+1,k)) * (h(i,j+1,k-1)*h(i+1,j+1,k-1)) ) wSW = G%mask2dCu(I-1,j) * ( (h(i,j ,k)*h(i-1,j ,k)) * (h(i,j ,k-1)*h(i-1,j ,k-1)) )