Merge branch 'dev/gfdl' into mono_N2_depth_units

src/core/MOM.F90

@@ -1090,7 +1090,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   call cpu_clock_end(id_clock_stoch)
   call cpu_clock_begin(id_clock_varT)
   if (CS%use_stochastic_EOS) then
-    call MOM_calc_varT(G, GV, h, CS%tv, CS%stoch_eos_CS, dt)
+    call MOM_calc_varT(G, GV, US, h, CS%tv, CS%stoch_eos_CS, dt)
     if (associated(CS%tv%varT)) call pass_var(CS%tv%varT, G%Domain, clock=id_clock_pass, halo=1)
   endif
   call cpu_clock_end(id_clock_varT)
@@ -3022,7 +3022,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   new_sim = is_new_run(restart_CSp)
   if (use_temperature) then
-    CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, US, param_file, diag, CS%stoch_eos_CS, restart_CSp)
+    CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS%stoch_eos_CS, restart_CSp)
   else
     CS%use_stochastic_EOS = .false.
   endif

src/core/MOM_forcing_type.F90

@@ -29,7 +29,7 @@ module MOM_forcing_type
 
 public extractFluxes1d, extractFluxes2d, optics_type
 public MOM_forcing_chksum, MOM_mech_forcing_chksum
-public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d
+public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, find_ustar
 public forcing_accumulate, fluxes_accumulate
 public forcing_SinglePointPrint, mech_forcing_diags, forcing_diagnostics
 public register_forcing_type_diags, allocate_forcing_type, deallocate_forcing_type
@@ -53,6 +53,12 @@ module MOM_forcing_type
   module procedure allocate_mech_forcing_from_ref
 end interface allocate_mech_forcing
 
+!> Determine the friction velocity from a forcing type or a mechanical forcing type.
+interface find_ustar
+  module procedure find_ustar_fluxes
+  module procedure find_ustar_mech_forcing
+end interface find_ustar
+
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
@@ -1077,6 +1083,139 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv,
 end subroutine calculateBuoyancyFlux2d
 
 
+!> Determine the friction velocity from the contenxts of a forcing type, perhaps
+!! using the evolving surface density.
+subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units)
+  type(ocean_grid_type),   intent(in)  :: G    !< The ocean's grid structure
+  type(verticalGrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),   intent(in)  :: US   !< A dimensional unit scaling type
+  type(forcing),           intent(in)  :: fluxes !< Surface fluxes container
+  type(thermo_var_ptrs),   intent(in)  :: tv   !< Structure containing pointers to any
+                                               !! available thermodynamic fields.
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] or
+                                               !! [H T-1 ~> m s-1 or kg m-2 s-1], depending on H_T_units.
+  integer,       optional, intent(in)  :: halo !< The extra halo size to fill in, 0 by default
+  logical,       optional, intent(in)  :: H_T_units !< If present and true, return U_star in units
+                                               !! of [H T-1 ~> m s-1 or kg m-2 s-1]
+
+  ! Local variables
+  real :: I_rho        ! The inverse of the reference density times a ratio of scaling
+                       ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
+                       ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3]
+  logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is
+                       ! returned in [H T-1 ~> m s-1 or kg m-2 s-1]
+  integer :: i, j, k, is, ie, js, je, hs
+
+  hs = 0 ; if (present(halo)) hs = max(halo, 0)
+  is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs
+
+  Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units
+
+  if (.not.(associated(fluxes%ustar) .or. associated(fluxes%tau_mag))) &
+    call MOM_error(FATAL, "find_ustar_fluxes requires that either ustar or tau_mag be associated.")
+
+  if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then
+    if (Z_T_units) then
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = fluxes%ustar(i,j)
+      enddo ; enddo
+    else
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = GV%Z_to_H * fluxes%ustar(i,j)
+      enddo ; enddo
+    endif
+  elseif (allocated(tv%SpV_avg)) then
+    if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, &
+        "find_ustar_fluxes called in non-Boussinesq mode with invalid values of SpV_avg.")
+    if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, &
+        "find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.")
+    if (Z_T_units) then
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1))
+      enddo ; enddo
+    else
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
+      enddo ; enddo
+    endif
+  else
+    I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H
+    if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0
+    do j=js,je ; do i=is,ie
+      U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho)
+    enddo ; enddo
+  endif
+
+end subroutine find_ustar_fluxes
+
+
+!> Determine the friction velocity from the contenxts of a forcing type, perhaps
+!! using the evolving surface density.
+subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_units)
+  type(ocean_grid_type),   intent(in)  :: G    !< The ocean's grid structure
+  type(verticalGrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),   intent(in)  :: US   !< A dimensional unit scaling type
+  type(mech_forcing),      intent(in)  :: forces !< Surface forces container
+  type(thermo_var_ptrs),   intent(in)  :: tv   !< Structure containing pointers to any
+                                               !! available thermodynamic fields.
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1]
+  integer,       optional, intent(in)  :: halo !< The extra halo size to fill in, 0 by default
+  logical,       optional, intent(in)  :: H_T_units !< If present and true, return U_star in units
+                                               !! of [H T-1 ~> m s-1 or kg m-2 s-1]
+
+  ! Local variables
+  real :: I_rho        ! The inverse of the reference density times a ratio of scaling
+                       ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
+                       ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3]
+  logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is
+                       ! returned in [H T-1 ~> m s-1 or kg m-2 s-1]
+  integer :: i, j, k, is, ie, js, je, hs
+
+  hs = 0 ; if (present(halo)) hs = max(halo, 0)
+  is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs
+
+  Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units
+
+  if (.not.(associated(forces%ustar) .or. associated(forces%tau_mag))) &
+    call MOM_error(FATAL, "find_ustar_mech requires that either ustar or tau_mag be associated.")
+
+  if (associated(forces%ustar) .and. (GV%Boussinesq .or. .not.associated(forces%tau_mag))) then
+    if (Z_T_units) then
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = forces%ustar(i,j)
+      enddo ; enddo
+    else
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = GV%Z_to_H * forces%ustar(i,j)
+      enddo ; enddo
+    endif
+  elseif (allocated(tv%SpV_avg)) then
+    if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, &
+        "find_ustar_mech called in non-Boussinesq mode with invalid values of SpV_avg.")
+    if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, &
+        "find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.")
+    if (Z_T_units) then
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1))
+      enddo ; enddo
+    else
+      do j=js,je ; do i=is,ie
+        U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1))
+      enddo ; enddo
+    endif
+  else
+    I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H
+    if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0
+    do j=js,je ; do i=is,ie
+      U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho)
+    enddo ; enddo
+  endif
+
+end subroutine find_ustar_mech_forcing
+
+
 !> Write out chksums for thermodynamic fluxes.
 subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
   character(len=*),        intent(in) :: mesg      !< message

src/core/MOM_isopycnal_slopes.F90

@@ -36,8 +36,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: e    !< Interface heights [Z ~> m]
   type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various
                                                                      !! thermodynamic variables
-  real,                                        intent(in)    :: dt_kappa_smooth !< A smoothing vertical diffusivity
-                                                                     !! times a smoothing timescale [Z2 ~> m2].
+  real,                                        intent(in)    :: dt_kappa_smooth !< A smoothing vertical
+                                                                     !! diffusivity times a smoothing
+                                                                     !! timescale [H Z ~> m2 or kg m-1]
   logical,                                     intent(in)    :: use_stanley !< turn on stanley param in slope
   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]
@@ -142,7 +143,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
 
 
   h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2
-  dz_neglect = GV%H_subroundoff * GV%H_to_Z
+  dz_neglect = GV%dZ_subroundoff
 
   local_open_u_BC = .false.
   local_open_v_BC = .false.
@@ -195,9 +196,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
 
   if (use_EOS) then
     if (present(halo)) then
-      call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, halo+1)
+      call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, halo+1)
     else
-      call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, 1)
+      call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, 1)
     endif
   endif
 
@@ -341,9 +342,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
         slope = slope * max(g%mask2dT(i,j),g%mask2dT(i+1,j))
       endif
       slope_x(I,j,K) = slope
-      if (present(dzSxN)) dzSxN(I,j,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) &
-                                                              + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N
-                                         * abs(slope) * G%mask2dCu(I,j) ! x-direction contribution to S^2
+      if (present(dzSxN)) &
+        dzSxN(I,j,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) &
+                                            + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N
+                       * abs(slope) * G%mask2dCu(I,j) ! x-direction contribution to S^2
 
     enddo ! I
   enddo ; enddo ! end of j-loop
@@ -477,9 +479,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
         slope = slope * max(g%mask2dT(i,j),g%mask2dT(i,j+1))
       endif
       slope_y(i,J,K) = slope
-      if (present(dzSyN)) dzSyN(i,J,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) &
-                                                              + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N
-                                         * abs(slope) * G%mask2dCv(i,J) ! x-direction contribution to S^2
+      if (present(dzSyN)) &
+        dzSyN(i,J,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) &
+                                             + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N
+                        * abs(slope) * G%mask2dCv(i,J) ! x-direction contribution to S^2
 
     enddo ! i
   enddo ; enddo ! end of j-loop
@@ -488,14 +491,15 @@ end subroutine calc_isoneutral_slopes
 
 !> Returns tracer arrays (nominally T and S) with massless layers filled with
 !! sensible values, by diffusing vertically with a small but constant diffusivity.
-subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, larger_h_denom)
+subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, larger_h_denom)
   type(ocean_grid_type),                     intent(in)  :: G    !< The ocean's grid structure
   type(verticalGrid_type),                   intent(in)  :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),                     intent(in)  :: US   !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: T_in !< Input temperature [C ~> degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: S_in !< Input salinity [S ~> ppt]
   real,                                      intent(in)  :: kappa_dt !< A vertical diffusivity to use for smoothing
-                                                                 !! times a smoothing timescale [Z2 ~> m2].
+                                                                 !! times a smoothing timescale [H Z ~> m2 or kg m-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T_f  !< Filled temperature [C ~> degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S_f  !< Filled salinity [S ~> ppt]
   integer,                         optional, intent(in)  :: halo_here !< Number of halo points to work on,
@@ -525,10 +529,15 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar
   is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke
 
   h_neglect = GV%H_subroundoff
-  kap_dt_x2 = (2.0*kappa_dt)*GV%Z_to_H**2
+  ! The use of the fixed rescaling factor in the next line avoids an extra call to thickness_to_dz()
+  ! and the use of an extra 3-d array of vertical distnaces across layers (dz).  This would be more
+  ! physically consistent, but it would also be more expensive, and given that this routine applies
+  ! a small (but arbitrary) amount of mixing to clean up the properties of nearly massless layers,
+  ! the added expense is hard to justify.
+  kap_dt_x2 = (2.0*kappa_dt) * (US%Z_to_m*GV%m_to_H) ! Usually the latter term is GV%Z_to_H.
   h0 = h_neglect
   if (present(larger_h_denom)) then
-    if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*GV%Z_to_H
+    if (larger_h_denom) h0 = 1.0e-16*sqrt(0.5*kap_dt_x2)
   endif
 
   if (kap_dt_x2 <= 0.0) then

src/core/MOM_stoch_eos.F90

@@ -40,7 +40,7 @@ module MOM_stoch_eos
   real :: stanley_coeff   !< Coefficient correlating the temperature gradient
                           !! and SGS T variance [nondim]; if <0, turn off scheme in all codes
   real :: stanley_a       !< a in exp(aX) in stochastic coefficient [nondim]
-  real :: kappa_smooth    !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1]
+  real :: kappa_smooth    !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
 
   !>@{ Diagnostic IDs
   integer :: id_stoch_eos  = -1, id_stoch_phi  = -1, id_tvar_sgs = -1
@@ -51,9 +51,10 @@ module MOM_stoch_eos
 contains
 
 !> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used.
-logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_CS)
+logical function MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS, restart_CS)
   type(time_type),         intent(in)    :: Time       !< Time for stochastic process
   type(ocean_grid_type),   intent(in)    :: G          !< The ocean's grid structure.
+  type(verticalGrid_type), intent(in)    :: GV         !< Vertical grid structure
   type(unit_scale_type),   intent(in)    :: US         !< A dimensional unit scaling type
   type(param_file_type),   intent(in)    :: param_file !< structure indicating parameter file to parse
   type(diag_ctrl), target, intent(inout) :: diag       !< Structure used to control diagnostics
@@ -80,7 +81,7 @@ logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_C
   call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", CS%kappa_smooth, &
                  "A diapycnal diffusivity that is used to interpolate "//&
                  "more sensible values of T & S into thin layers.", &
-                 units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s, &
+                 units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T, &
                  do_not_log=(CS%stanley_coeff<0.0))
 
   ! Don't run anything if STANLEY_COEFF < 0
@@ -193,9 +194,10 @@ subroutine post_stoch_EOS_diags(CS, tv, diag)
 end subroutine post_stoch_EOS_diags
 
 !> Computes a parameterization of the SGS temperature variance
-subroutine MOM_calc_varT(G, GV, h, tv, CS, dt)
+subroutine MOM_calc_varT(G, GV, US, h, tv, CS, dt)
   type(ocean_grid_type),   intent(in)   :: G   !< The ocean's grid structure.
   type(verticalGrid_type), intent(in)   :: GV  !< Vertical grid structure
+  type(unit_scale_type),   intent(in)   :: US  !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: h   !< Layer thickness [H ~> m]
   type(thermo_var_ptrs),  intent(inout) :: tv  !< Thermodynamics structure
@@ -219,7 +221,7 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt)
   ! extreme gradients along layers which are vanished against topography. It is
   ! still a poor approximation in the interior when coordinates are strongly tilted.
   if (.not. associated(tv%varT)) allocate(tv%varT(G%isd:G%ied, G%jsd:G%jed, GV%ke), source=0.0)
-  call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.)
+  call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo_here=1, larger_h_denom=.true.)
 
   do k=1,G%ke
     do j=G%jsc,G%jec

src/parameterizations/lateral/MOM_MEKE.F90

@@ -1578,7 +1578,8 @@ subroutine ML_MEKE_calculate_features(G, GV, US, CS, Rd_dx_h, u, v, tv, h, dt, f
     h_v(i,J,k) = 0.5*(h(i,j,k)*G%mask2dT(i,j) + h(i,j+1,k)*G%mask2dT(i,j+1)) + GV%Angstrom_H
   enddo; enddo; enddo;
   call find_eta(h, tv, G, GV, US, e, halo_size=2)
-  call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7, .false., slope_x, slope_y)
+  ! Note the hard-coded dimenisional constant in the following line.
+  call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7*GV%m2_s_to_HZ_T, .false., slope_x, slope_y)
   call pass_vector(slope_x, slope_y, G%Domain)
   do j=js-1,je+1; do i=is-1,ie+1
     slope_x_vert_avg(I,j) = vertical_average_interface(slope_x(i,j,:), h_u(i,j,:), GV%H_subroundoff)