Skip to content

Commit

Permalink
Change units of slope returned from calc_isoneutral_slopes() to "Z L-1"
Browse files Browse the repository at this point in the history 
- 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.
  • Loading branch information
@adcroft
adcroft committed Mar 12, 2021
1 parent 250f007 commit 0d60fd0
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 36 deletions.
26 changes: 13 additions & 13 deletions src/core/MOM_isopycnal_slopes.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -86,15 +86,15 @@ 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].
real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
10 changes: 6 additions & 4 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
38 changes: 19 additions & 19 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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].
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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.", &
Expand Down Expand Up @@ -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', &
Expand Down

0 comments on commit 0d60fd0

Please sign in to comment.