(+)Better units documentation in tracer_hordiff

  This commit includes minor cleanup of the units documentation in several of
the process parameterization routines.  These changes include:

- Change the units of the (apparently unused) optional read_khdt_[xy] arguments to
  tracer_hordiff from [m2] to [L2 ~> m2]

- Correct a tiny number used for safety in the denominator of a normalization
  factor based on the (nondimensional) masks from H_subroundoff to a tiny
  nondimensional number, giving an expression that is dimensionally consistent.
  Without this change, the diagnostic "KHTR_h" will vary with large enough
  negative values of H_RESCALE_POWER.

- Replace a local variable in  bulkmixedlayer_init with units of time with
  distinct get_param calls setting the default for BUFFER_LAY_DETRAIN_TIME
  depending on which buffer layer detrainment scheme is used.

- Correct the dimensions of the "land_mask" diagnostic internal_tides_init
  from "logical" to "nondim", because it is a real multiplicable mask, not a
  boolean logical mask.

- Eliminate the internal variable L2_to_Z2 in the calculate_projected_state
  routine in MOM_kappa_shear which had been described incorrectly and ultimately
  was not helpful.

- Corrected the allocated size declaration for three 3-d v-face arrays

- Add comments highlighting a non-stride-1 loop structure for openMP threading
  in tracer_epipycnal_ML_diff that required the addition of 3 extra 3-d arrays
  and might be a performance liability in non-openMP configurations.

- Extensive modification of the comments describing the variables in
  MOM_tracer_hordiff to clearly indicate the units of each real variable and
  similarly of one variable in MOM_geothermal.

All answers are bitwise identical unless very large rescaling is done for
thickness units, and the only change in model output is to the documented units of
one diagnostic that is not used often, and to the diagnostic "KHTR_h" when large
enough negative values of H_RESCALE_POWER are used.

src/parameterizations/lateral/MOM_internal_tides.F90

@@ -2346,8 +2346,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
                  default=.false.)
   call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
                  "CDRAG is the drag coefficient relating the magnitude of "//&
-                 "the velocity field to the bottom stress.", units="nondim", &
-                 default=0.003)
+                 "the velocity field to the bottom stress.", &
+                 units="nondim", default=0.003)
   call get_param(param_file, mdl, "INTERNAL_TIDE_ENERGIZED_ANGLE", CS%energized_angle, &
                  "If positive, only one angular band of the internal tides "//&
                  "gets all of the energy.  (This is for debugging.)", default=-1)
@@ -2518,7 +2518,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
   CS%id_dy_Cu = register_diag_field('ocean_model', 'dy_Cu', diag%axesT1, &
                  Time, 'East face unblocked width', 'm', conversion=US%L_to_m)
   CS%id_land_mask = register_diag_field('ocean_model', 'land_mask', diag%axesT1, &
-                 Time, 'Land mask', 'logical')            ! used if overriding (BDM)
+                 Time, 'Land mask', 'nondim')
   ! Output reflection parameters as diags here (not needed every timestep)
   if (CS%id_refl_ang > 0)   call post_data(CS%id_refl_ang, CS%refl_angle, CS%diag)
   if (CS%id_refl_pref > 0)  call post_data(CS%id_refl_pref, CS%refl_pref, CS%diag)

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -3360,7 +3360,6 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_mixed_layer"  ! This module's name.
-  real :: BL_detrain_time_dflt ! The default value for BUFFER_LAY_DETRAIN_TIME [s]
   real :: omega_frac_dflt  ! The default value for ML_OMEGA_FRAC [nondim]
   real :: ustar_min_dflt   ! The default value for BML_USTAR_MIN [Z T-1 ~> m s-1]
   real :: Hmix_min_m       ! The unscaled value of HMIX_MIN [m]
@@ -3454,10 +3453,15 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
                  "The minimum buffer layer thickness relative to the combined mixed "//&
                  "land buffer ayer thicknesses when they are thin.", &
                  units="nondim", default=0.1/CS%nkbl)
-  BL_detrain_time_dflt = 4.0*3600.0 ; if (CS%nkbl==1) BL_detrain_time_dflt = 86400.0*30.0
-  call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", CS%BL_detrain_time, &
+  if (CS%nkbl==1) then
+    call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", CS%BL_detrain_time, &
                  "A timescale that characterizes buffer layer detrainment events.", &
-                 units="s", default=BL_detrain_time_dflt, scale=US%s_to_T)
+                 units="s", default=86400.0*30.0, scale=US%s_to_T)
+  else
+    call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", CS%BL_detrain_time, &
+                 "A timescale that characterizes buffer layer detrainment events.", &
+                 units="s", default=4.0*3600.0, scale=US%s_to_T)
+  endif
   call get_param(param_file, mdl, "BUFFER_SPLIT_RHO_TOL", CS%BL_split_rho_tol, &
                  "The fractional tolerance for matching layer target densities when splitting "//&
                  "layers to deal with massive interior layers that are lighter than one of the "//&

src/parameterizations/vertical/MOM_geothermal.F90

@@ -107,7 +107,7 @@ subroutine geothermal_entraining(h, tv, dt, ea, eb, G, GV, US, CS, halo)
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
     T_old, & ! Temperature of each layer before any heat is added, for diagnostics [C ~> degC]
     h_old, & ! Thickness of each layer before any heat is added, for diagnostics [H ~> m or kg m-2]
-    work_3d ! Scratch variable used to calculate changes due to geothermal
+    work_3d ! Scratch variable used to calculate changes due to geothermal [various]
   real :: Idt           ! inverse of the timestep [T-1 ~> s-1]
 
   logical :: do_i(SZI_(G))
@@ -407,7 +407,7 @@ subroutine geothermal_in_place(h, tv, dt, G, GV, US, CS, halo)
   if (.not.associated(tv%T)) call MOM_error(FATAL, "MOM geothermal_in_place: "//&
       "Geothermal heating can only be applied if T & S are state variables.")
 
-!  do i=is,ie ; do j=js,je
+!  do j=js,je ; do i=is,ie
 !    resid(i,j) = tv%internal_heat(i,j)
 !  enddo ; enddo
 
@@ -573,17 +573,17 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith
   if (id > 0) call post_data(id, CS%geo_heat, diag, .true.)
 
   ! Diagnostic for tendencies due to internal heat (in 3d)
-  CS%id_internal_heat_heat_tendency=register_diag_field('ocean_model', &
+  CS%id_internal_heat_heat_tendency = register_diag_field('ocean_model', &
         'internal_heat_heat_tendency', diag%axesTL, Time,              &
         'Heat tendency (in 3D) due to internal (geothermal) sources',  &
         'W m-2', conversion=US%QRZ_T_to_W_m2, v_extensive=.true.)
-  CS%id_internal_heat_temp_tendency=register_diag_field('ocean_model', &
+  CS%id_internal_heat_temp_tendency = register_diag_field('ocean_model', &
         'internal_heat_temp_tendency', diag%axesTL, Time,              &
         'Temperature tendency (in 3D) due to internal (geothermal) sources', &
         'degC s-1', conversion=US%C_to_degC*US%s_to_T, v_extensive=.true.)
   if (.not.useALEalgorithm) then
     ! Do not offer this diagnostic if heating will be in place.
-    CS%id_internal_heat_h_tendency=register_diag_field('ocean_model',    &
+    CS%id_internal_heat_h_tendency = register_diag_field('ocean_model',    &
         'internal_heat_h_tendency', diag%axesTL, Time,                &
         'Thickness tendency (in 3D) due to internal (geothermal) sources', &
         trim(thickness_units)//' s-1', conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)

src/parameterizations/vertical/MOM_kappa_shear.F90

@@ -1061,8 +1061,6 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int
 
   ! Local variables
   real, dimension(nz+1) :: c1 ! A tridiagonal variable [nondim]
-  real :: L2_to_Z2   ! A conversion factor from horizontal length units to vertical depth
-                     ! units squared [Z2 s2 T-2 m-2 ~> 1].
   real :: a_a, a_b   ! Tridiagonal coupling coefficients [Z ~> m]
   real :: b1, b1nz_0 ! Tridiagonal variables [Z-1 ~> m-1]
   real :: bd1        ! A term in the denominator of b1 [Z ~> m]
@@ -1134,16 +1132,14 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int
   endif
 
   ! Store the squared shear at interfaces
-  ! L2_to_Z2 = US%m_to_Z**2 * US%T_to_s**2
-  L2_to_Z2 = US%L_to_Z**2
   S2(1) = 0.0 ; S2(nz+1) = 0.0
   if (ks > 1) &
-    S2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (L2_to_Z2*I_dz_int(ks)**2)
+    S2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (US%L_to_Z*I_dz_int(ks))**2
   do K=ks+1,ke
-    S2(K) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (L2_to_Z2*I_dz_int(K)**2)
+    S2(K) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (US%L_to_Z*I_dz_int(K))**2
   enddo
   if (ke<nz) &
-    S2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * (L2_to_Z2*I_dz_int(ke+1)**2)
+    S2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * (US%L_to_Z*I_dz_int(ke+1))**2
 
   ! Store the buoyancy frequency at interfaces
   N2(1) = 0.0 ; N2(nz+1) = 0.0