From 28b1ac9224da9fa5c81275ecb4132ed0c80121fb Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 2 Mar 2021 14:09:56 -0500 Subject: [PATCH 1/5] Update CVMix to v0.93-beta (no API change) - Tag v0.93-beta of CVMix is the last tag on their "master" branch before an API change within CVMix. - Answers reproduce with the three-year prior commit being used, as tested in MOM6-examples single_column cases. These do not cover all CVMix code. --- pkg/CVMix-src | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 534fc38e75..919a3a24ac 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 534fc38e759fcb8a2563fa0dc4c0bbf81f758db9 +Subproject commit 919a3a24acd67454df3ef8ad9854b734a472fea7 From 34dc0c848916075614eb7ec0a5b32a6f9fb6ebec Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 2 Mar 2021 14:16:39 -0500 Subject: [PATCH 2/5] Update CVMix to v0.94b-beta with API change - This tag involves an API change that requires changes to MOM_tidal_mixing.F90. - Changes are dropped arguments so presumably that data is not needed. - Answers reproduce in so far as they are covered by MOM6-examples. Not all of CVMix is covered in these tests and this code in particular is not. --- pkg/CVMix-src | 2 +- src/parameterizations/vertical/MOM_tidal_mixing.F90 | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 919a3a24ac..fee4701ac6 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 919a3a24acd67454df3ef8ad9854b734a472fea7 +Subproject commit fee4701ac61b69964850db5a4c3ebea41fb9346f diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index b870dff1af..512179445b 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -882,7 +882,6 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv ! summing q_i*TidalConstituent_i over the number of constituents. call CVMix_compute_SchmittnerCoeff( nlev = GV%ke, & energy_flux = tidal_qe_md(:), & - rho = rho_fw, & SchmittnerCoeff = Schmittner_coeff, & exp_hab_zetar = exp_hab_zetar, & CVmix_tidal_params_user = CS%CVMix_tidal_params) @@ -896,7 +895,6 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv Tdiff_out = Kd_tidal, & Nsqr = N2_int_i, & OceanDepth = -iFaceHeight(GV%ke+1), & - vert_dep = vert_dep, & nlev = GV%ke, & max_nlev = GV%ke, & SchmittnerCoeff = Schmittner_coeff, & From e56d4533587fdaac26f3f7faee54ab391cd6ad61 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 2 Mar 2021 14:22:54 -0500 Subject: [PATCH 3/5] Update CVMix to v0.98-beta - v0.98-beta is the latest tag of CVMix that is available a.t.t. - There were no API changes since the v0.94b-beta that affected MOM6. - No answer changes for MOM6-examples, but these do not exercise much of CVMix. --- pkg/CVMix-src | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index fee4701ac6..9423197f89 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit fee4701ac61b69964850db5a4c3ebea41fb9346f +Subproject commit 9423197f894112edfcb1502245f7d7b873d551f9 From 651b4673f46a83f5fc741feaafe9cc0586e4bb0b Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sat, 6 Mar 2021 01:19:08 -0500 Subject: [PATCH 4/5] Bugfix: FGNV streamfunction vertical bounds This patch fixes an issue with the vertical array bounds of the Ferrari et al. streamfunction. The array is bounded across interfaces, from 1 to nz+1, but only the interior values need to be determined due to the arbitrary boundary value (set here to zero). In the current source, the streamfunction is rescaled before calling streamfn_solver, but need not be applied to the boundary values. This is unlikely to cause errors in production, since the values are later reset to zero, but the rescaling can raise errors in more aggressive debugging builds, such as when the arrays are initialized with NaN values. --- .../lateral/MOM_thickness_diffuse.F90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 8c6a90ba9c..4d602f9a67 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -993,10 +993,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. do I=is-1,ie if (G%mask2dCu(I,j)>0.) then - Sfn_unlim_u(I,:) = ( 1. + CS%FGNV_scale ) * Sfn_unlim_u(I,:) + do K=2,nz + Sfn_unlim_u(I,K) = (1. + CS%FGNV_scale) * Sfn_unlim_u(I,K) + enddo call streamfn_solver(nz, c2_h_u(I,:), hN2_u(I,:), Sfn_unlim_u(I,:)) else - Sfn_unlim_u(I,:) = 0. + do K=2,nz + Sfn_unlim_u(I,K) = 0. + enddo endif enddo endif @@ -1259,10 +1263,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. do i=is,ie if (G%mask2dCv(i,J)>0.) then - Sfn_unlim_v(i,:) = ( 1. + CS%FGNV_scale ) * Sfn_unlim_v(i,:) + do K=2,nz + Sfn_unlim_v(i,K) = (1. + CS%FGNV_scale) * Sfn_unlim_v(i,K) + enddo call streamfn_solver(nz, c2_h_v(i,:), hN2_v(i,:), Sfn_unlim_v(i,:)) else - Sfn_unlim_v(i,:) = 0. + do K=2,nz + Sfn_unlim_v(i,K) = 0. + enddo endif enddo endif From 0d60fd0264ee271b1ce2a7b6c5e92d853c0c2769 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 12 Mar 2021 09:39:17 -0500 Subject: [PATCH 5/5] Change units of slope returned from calc_isoneutral_slopes() to "Z L-1" - Units of isoneutral or interface slope were recorded as "nondim". While true in SI units, not so for MOM6 units. MOM6 distinguishes between units of length in the vertical (Z) and horizontal (L) the slopes should have units of "Z L-1 ~> nondim". - This has consequences for other variables in calc_isoneutral_slopes(). - An internal constant, G_Rho0, was defined differently from elsewhere in the code. "g" has units of "L2 Z-1 T-2 ~ m s-2" because it is the vertical component of the gradient of geopotential in "L2 T-2 ~ m2 s-2". Everywhere else `G_Rho0 = g_Earth/Rho0` but in this routine it was different in order render N2 (the Brunt-Vaisala frequency) in units of "T-2" (s-2). - N2 is a quantity associated with dispersion relations and defined N2 = - g/Rho0 d/dz rho and either way acquires units of "L2 Z-2 T-2" and not just "T-2". In SI units L2 Z-2 = 1. So I have also changed the units of N2 in this, and connected, modules. - The changes also propagate to MOM_lateral_mixing_coeffs.F90 and MOM_thickness_diffuse.F90. - Changing the definition of G_Rho0 in calc_isoneutral_slopes(), and its units to "L2 Z-1 T-2", the slope and N2 calculations then require many less inline conversions. Many of the one-line changes in this commit remove factors like US%Z_to_L.There is one exception: - In the calculation of slope, we use in the denominator a mostly non-vanishing replacement for d/dz rho, the magnitude of grad rho from mag_grad2 = ( d/dx rho )^2 + ( d/dz rho )^2. In code this had `mag_grad2 = drdy**2 + (L_to_Z*drdz)**2` since this is mixing gradients in the horizontal and vertical. The result should be in "R2 Z-2" so now `mag_grad2 = (Z_to_L*drdy)**2 + drdz**2` - A few diagnostics needed new, or changed, conversion factors. - One run-time parameter needed a conversion parameter. - For the most part this commit moves inline conversions of units to the I/O stage, which is an indicator that it is the right thing to do. --- src/core/MOM_isopycnal_slopes.F90 | 26 ++++++------- .../lateral/MOM_lateral_mixing_coeffs.F90 | 10 +++-- .../lateral/MOM_thickness_diffuse.F90 | 38 +++++++++---------- 3 files changed, 38 insertions(+), 36 deletions(-) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index e1f573f6ea..98b5b10998 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -37,14 +37,14 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !! thermodynamic variables real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity !! times a smoothing timescale [Z2 ~> m2]. - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction [nondim] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction [nondim] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & optional, intent(inout) :: N2_u !< Brunt-Vaisala frequency squared at - !! interfaces between u-points [T-2 ~> s-2] + !! interfaces between u-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), & optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at - !! interfaces between u-points [T-2 ~> s-2] + !! interfaces between v-points [L2 Z-2 T-2 ~> s-2] integer, optional, intent(in) :: halo !< Halo width over which to compute type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -86,7 +86,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4]. real :: Slope ! The slope of density surfaces, calculated in a way ! that is always between -1 and 1. - real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8]. + real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 Z-2 ~> kg2 m-8]. real :: slope2_Ratio ! The ratio of the slope squared to slope_max squared. 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]. @@ -94,7 +94,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & real :: dz_neglect ! A change in interface heighs that is so small it is usually lost ! in roundoff and can be neglected [Z ~> m]. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. - real :: G_Rho0 ! The gravitational acceleration divided by density [Z2 T-2 R-1 ~> m5 kg-2 s-2] + real :: G_Rho0 ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] real :: Z_to_L ! A conversion factor between from units for e to the ! units for lateral distances. real :: L_to_Z ! A conversion factor between from units for lateral distances @@ -134,7 +134,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & present_N2_u = PRESENT(N2_u) present_N2_v = PRESENT(N2_v) - G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / GV%Rho0 + G_Rho0 = GV%g_Earth / GV%Rho0 if (present_N2_u) then do j=js,je ; do I=is-1,ie N2_u(I,j,1) = 0. @@ -248,17 +248,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdx**2 + (L_to_Z*drdz)**2 + mag_grad2 = (Z_to_L*drdx)**2 + drdz**2 if (mag_grad2 > 0.0) then 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 - if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy frequency [T-2 ~> s-2] + if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. - slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j) + slope_x(I,j,K) = (e(i,j,K)-e(i+1,j,K)) * G%IdxCu(I,j) endif if (local_open_u_BC) then l_seg = OBC%segnum_u(I,j) @@ -351,17 +351,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdy**2 + (L_to_Z*drdz)**2 + mag_grad2 = (Z_to_L*drdy)**2 + drdz**2 if (mag_grad2 > 0.0) then 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 - if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy frequency [T-2 ~> s-2] + if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. - slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J) + slope_y(i,J,K) = (e(i,j,K)-e(i,j+1,K)) * G%IdyCv(i,J) endif if (local_open_v_BC) then l_seg = OBC%segnum_v(i,J) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e3a6f1599e..7d95c43b98 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1125,16 +1125,18 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) then CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, & 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', & - 's-2', conversion=US%s_to_T**2) + 's-2', conversion=(US%L_to_Z*US%s_to_T)**2) CS%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, Time, & 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', & - 's-2', conversion=US%s_to_T**2) + 's-2', conversion=(US%L_to_Z*US%s_to_T)**2) endif if (CS%use_stored_slopes) then CS%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, Time, & - 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', 'nondim') + 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', & + 'nondim', conversion=US%Z_to_L**2) CS%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, Time, & - 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 'nondim') + 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', & + 'nondim', conversion=US%Z_to_L**2) endif oneOrTwo = 1.0 diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 8c6a90ba9c..99ecca9745 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -40,7 +40,7 @@ module MOM_thickness_diffuse real :: max_Khth_CFL !< Maximum value of the diffusive CFL for thickness diffusion real :: Khth_Min !< Minimum value of Khth [L2 T-1 ~> m2 s-1] real :: Khth_Max !< Maximum value of Khth [L2 T-1 ~> m2 s-1], or 0 for no max - real :: slope_max !< Slopes steeper than slope_max are limited in some way [nondim]. + real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim]. real :: kappa_smooth !< Vertical diffusivity used to interpolate more !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1]. logical :: thickness_diffuse !< If true, interfaces heights are diffused. @@ -83,8 +83,8 @@ module MOM_thickness_diffuse type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics real, pointer :: GMwork(:,:) => NULL() !< Work by thickness diffusivity [R Z L2 T-3 ~> W m-2] - real, pointer :: diagSlopeX(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [nondim] - real, pointer :: diagSlopeY(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [nondim] + real, pointer :: diagSlopeX(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [Z L-1 ~> nondim] + real, pointer :: diagSlopeY(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [Z L-1 ~> nondim] real, dimension(:,:,:), pointer :: & KH_u_GME => NULL(), & !< interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1] @@ -578,8 +578,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !! the isopycnal slopes are taken directly from !! the interface slopes without consideration of !! density gradients [nondim]. - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopyc. slope at u [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopyc. slope at v [Z L-1 ~> nondim] ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & T, & ! The temperature (or density) [degC], with the values in @@ -660,7 +660,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real :: Slope ! The slope of density surfaces, calculated in a way ! that is always between -1 and 1, nondimensional. real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8]. - real :: I_slope_max2 ! The inverse of slope_max squared, nondimensional. + real :: I_slope_max2 ! The inverse of slope_max squared [L2 Z-2 ~> nondim]. 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 :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]. @@ -919,7 +919,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdx**2 + (US%L_to_Z*drdz)**2 + mag_grad2 = (US%Z_to_L*drdx)**2 + drdz**2 if (mag_grad2 > 0.0) then Slope = drdx / sqrt(mag_grad2) slope2_Ratio_u(I,K) = Slope**2 * I_slope_max2 @@ -933,7 +933,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! that ignore density gradients along layers. if (present_int_slope_u) then Slope = (1.0 - int_slope_u(I,j,K)) * Slope + & - int_slope_u(I,j,K) * US%Z_to_L*((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) + int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) endif @@ -942,7 +942,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1]. - Sfn_unlim_u(I,K) = -((KH_u(I,j,K)*G%dy_Cu(I,j))*US%L_to_Z*Slope) + Sfn_unlim_u(I,K) = -((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) ! Avoid moving dense water upslope from below the level of ! the bottom on the receiving side. @@ -968,10 +968,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (present_slope_x) then Slope = slope_x(I,j,k) else - Slope = US%Z_to_L*((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j) + Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j) endif if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope - Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*US%L_to_Z*Slope) + Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) hN2_u(I,K) = GV%g_prime(K) endif ! if (use_EOS) else ! if (k > nk_linear) @@ -1185,7 +1185,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdy**2 + (US%L_to_Z*drdz)**2 + mag_grad2 = (US%Z_to_L*drdy)**2 + drdz**2 if (mag_grad2 > 0.0) then Slope = drdy / sqrt(mag_grad2) slope2_Ratio_v(i,K) = Slope**2 * I_slope_max2 @@ -1199,7 +1199,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! that ignore density gradients along layers. if (present_int_slope_v) then Slope = (1.0 - int_slope_v(i,J,K)) * Slope + & - int_slope_v(i,J,K) * US%Z_to_L*((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) + int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) endif @@ -1208,7 +1208,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1]. - Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*US%L_to_Z*Slope) + Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) ! Avoid moving dense water upslope from below the level of ! the bottom on the receiving side. @@ -1234,10 +1234,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (present_slope_y) then Slope = slope_y(i,J,k) else - Slope = US%Z_to_L*((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J) + Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J) endif if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope - Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*US%L_to_Z*Slope) + Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) hN2_v(i,K) = GV%g_prime(K) endif ! if (use_EOS) else ! if (k > nk_linear) @@ -1947,7 +1947,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "longer than DT, or 0 to use DT.", units="s", default=0.0, scale=US%s_to_T) call get_param(param_file, mdl, "KHTH_SLOPE_MAX", CS%slope_max, & "A slope beyond which the calculated isopycnal slope is "//& - "not reliable and is scaled away.", units="nondim", default=0.01) + "not reliable and is scaled away.", units="nondim", default=0.01, scale=US%L_to_Z) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & @@ -2065,10 +2065,10 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) CS%id_slope_x = register_diag_field('ocean_model', 'neutral_slope_x', diag%axesCui, Time, & - 'Zonal slope of neutral surface', 'nondim') + 'Zonal slope of neutral surface', 'nondim', conversion=US%Z_to_L) if (CS%id_slope_x > 0) call safe_alloc_ptr(CS%diagSlopeX,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke+1) CS%id_slope_y = register_diag_field('ocean_model', 'neutral_slope_y', diag%axesCvi, Time, & - 'Meridional slope of neutral surface', 'nondim') + 'Meridional slope of neutral surface', 'nondim', conversion=US%Z_to_L) if (CS%id_slope_y > 0) call safe_alloc_ptr(CS%diagSlopeY,G%isd,G%ied,G%JsdB,G%JedB,GV%ke+1) CS%id_sfn_x = register_diag_field('ocean_model', 'GM_sfn_x', diag%axesCui, Time, & 'Parameterized Zonal Overturning Streamfunction', &