Bugfix on total column thickness for wetting

src/ALE/MOM_hybgen_unmix.F90

@@ -214,7 +214,7 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
     endif
 
     ! The following block of code is used to trigger z* stretching of the targets heights.
-    if (allocated(tv%SpV_avg)) then  ! This is the fully non-Boussiesq version
+    if (allocated(tv%SpV_avg)) then  ! This is the fully non-Boussinesq version
       dz_tot = 0.0
       do k=1,nk
         dz_tot = dz_tot + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h_col(k)

src/core/MOM_barotropic.F90

@@ -5533,6 +5533,7 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
                                        ! name in wave_drag_file.
   real :: mean_SL     ! The mean sea level that is used along with the bathymetry to estimate the
                       ! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m].
+  real :: htot        ! Total column thickness used when BT_NONLIN_STRESS is false [Z ~> m].
   real :: Z_to_H      ! A local unit conversion factor [H Z-1 ~> nondim or kg m-3]
   real :: H_to_Z      ! A local unit conversion factor [Z H-1 ~> nondim or m3 kg-1]
   real :: det_de      ! The partial derivative due to self-attraction and loading of the reference
@@ -6437,15 +6438,17 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
     Mean_SL = G%Z_ref
     Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
     do j=js,je ; do I=is-1,ie
-      if (G%OBCmaskCu(I,j) > 0.) then
-        CS%IDatu(I,j) = G%OBCmaskCu(I,j) * 2.0 / (Z_to_H * ((G%bathyT(i+1,j) + G%bathyT(i,j)) + 2.0*Mean_SL))
+      htot = max(G%bathyT(i+1,j) + G%Z_ref, 0.0) + max(G%bathyT(i,j) + G%Z_ref, 0.0)
+      if (G%OBCmaskCu(I,j) * htot > 0.) then
+        CS%IDatu(I,j) = G%OBCmaskCu(I,j) * 2.0 / (Z_to_H * htot)
       else ! Both neighboring H points are masked out or this is an OBC face so IDatu(I,j) is unused
         CS%IDatu(I,j) = 0.
       endif
     enddo ; enddo
     do J=js-1,je ; do i=is,ie
-      if (G%OBCmaskCv(i,J) > 0.) then
-        CS%IDatv(i,J) = G%OBCmaskCv(i,J) * 2.0 / (Z_to_H * ((G%bathyT(i,j+1) + G%bathyT(i,j)) + 2.0*Mean_SL))
+      htot = max(G%bathyT(i,j+1) + G%Z_ref, 0.0) + max(G%bathyT(i,j) + G%Z_ref, 0.0)
+      if (G%OBCmaskCv(i,J) * htot > 0.) then
+        CS%IDatv(i,J) = G%OBCmaskCv(i,J) * 2.0 / (Z_to_H * htot)
       else ! Both neighboring H points are masked out or this is an OBC face so IDatv(i,J) is unused
         CS%IDatv(i,J) = 0.
       endif

src/diagnostics/MOM_wave_speed.F90

@@ -550,8 +550,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N
                 ! Determine whether N2 estimates should not be allowed to increase with depth.
                 if (l_mono_N2_column_fraction>0.) then
                   if (GV%Boussinesq .or. GV%semi_Boussinesq) then
-                    below_mono_N2_frac = ((G%bathyT(i,j)+G%Z_ref) - GV%H_to_Z*sum_hc < &
-                                          l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref))
+                    below_mono_N2_frac = (max(G%bathyT(i,j)+G%Z_ref, 0.0) - GV%H_to_Z*sum_hc < &
+                                          l_mono_N2_column_fraction*max(G%bathyT(i,j)+G%Z_ref, 0.0))
                   else
                     below_mono_N2_frac = (htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i))
                   endif

src/parameterizations/lateral/MOM_MEKE.F90

@@ -372,12 +372,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       if (GV%Boussinesq) then
         !$OMP parallel do default(shared)
         do j=js-1,je+1 ; do i=is-1,ie+1
-          depth_tot(i,j) = (G%bathyT(i,j) + G%Z_ref) * GV%Z_to_H
+          depth_tot(i,j) = max(G%bathyT(i,j) + G%Z_ref, 0.0) * GV%Z_to_H
         enddo ; enddo
       else
         !$OMP parallel do default(shared)
         do j=js-1,je+1 ; do i=is-1,ie+1
-          depth_tot(i,j) = (G%bathyT(i,j) + G%Z_ref) * CS%rho_fixed_total_depth * GV%RZ_to_H
+          depth_tot(i,j) = max(G%bathyT(i,j) + G%Z_ref, 0.0) * CS%rho_fixed_total_depth * GV%RZ_to_H
         enddo ; enddo
       endif
     else

src/parameterizations/lateral/MOM_internal_tides.F90

@@ -3773,7 +3773,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
     ! Restrict RMS topographic roughness to a fraction (10 percent by default) of the column depth.
     if (RMS_roughness_frac >= 0.0) then
-      h2(i,j) = max(min((RMS_roughness_frac*(G%bathyT(i,j)+G%Z_ref))**2, h2(i,j)), 0.0)
+      h2(i,j) = max(min((RMS_roughness_frac * max(G%bathyT(i,j)+G%Z_ref, 0.0))**2, h2(i,j)), 0.0)
     else
       h2(i,j) = max(h2(i,j), 0.0)
     endif

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -224,11 +224,13 @@ subroutine calc_depth_function(G, CS)
   expo = CS%depth_scaled_khth_exp
 !$OMP do
   do j=js,je ; do I=is-1,Ieq
-    CS%Depth_fn_u(I,j) = (MIN(1.0, (0.5*(G%bathyT(i,j) + G%bathyT(i+1,j)) + G%Z_ref)/H0))**expo
+    CS%Depth_fn_u(I,j) = (MIN(1.0, &
+      (0.5 * (max(G%bathyT(i,j) + G%Z_ref, 0.0) + max(G%bathyT(i+1,j) + G%Z_ref, 0.0))) / H0))**expo
   enddo ; enddo
 !$OMP do
   do J=js-1,Jeq ; do i=is,ie
-    CS%Depth_fn_v(i,J) = (MIN(1.0, (0.5*(G%bathyT(i,j) + G%bathyT(i,j+1)) + G%Z_ref)/H0))**expo
+    CS%Depth_fn_v(i,J) = (MIN(1.0, &
+      (0.5 * (max(G%bathyT(i,j) + G%Z_ref, 0.0) + max(G%bathyT(i,j+1) + G%Z_ref, 0.0))) / H0))**expo
   enddo ; enddo
 
 end subroutine calc_depth_function

src/parameterizations/vertical/MOM_internal_tide_input.F90

@@ -557,7 +557,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
 
     ! Restrict rms topo to a fraction (often 10 percent) of the column depth.
     if (max_frac_rough >= 0.0) &
-      itide%h2(i,j) = min((max_frac_rough*(G%bathyT(i,j)+G%Z_ref))**2, itide%h2(i,j))
+      itide%h2(i,j) = min((max_frac_rough * max(G%bathyT(i,j)+G%Z_ref, 0.0))**2, itide%h2(i,j))
 
     ! Compute the fixed part of internal tidal forcing; units are [R Z4 H-1 T-2 ~> J m-2 or J m kg-1] here.
     CS%TKE_itidal_coef(i,j,fr) = 0.5*US%L_to_Z*kappa_h2_factor * GV%H_to_RZ * &

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -366,12 +366,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
 
   !$OMP parallel do default(shared)
   do J=js-1,je ; do i=is-1,ie+1
-    D_v(i,J) = 0.5*(G%bathyT(i,j) + G%bathyT(i,j+1)) + G%Z_ref
+    D_v(i,J) = 0.5 * (max(G%bathyT(i,j) + G%Z_ref, 0.0) + max(G%bathyT(i,j+1) + G%Z_ref, 0.0))
     mask_v(i,J) = G%mask2dCv(i,J)
   enddo ; enddo
   !$OMP parallel do default(shared)
   do j=js-1,je+1 ; do I=is-1,ie
-    D_u(I,j) = 0.5*(G%bathyT(i,j) + G%bathyT(i+1,j)) + G%Z_ref
+    D_u(I,j) = 0.5 * (max(G%bathyT(i,j) + G%Z_ref, 0.0) + max(G%bathyT(i+1,j) + G%Z_ref, 0.0))
     mask_u(I,j) = G%mask2dCu(I,j)
   enddo ; enddo
 

src/parameterizations/vertical/MOM_tidal_mixing.F90

@@ -516,7 +516,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di
         CS%h2(i,j) = hamp*hamp
       else
         if (max_frac_rough >= 0.0) &
-          CS%h2(i,j) = min((max_frac_rough*(G%bathyT(i,j)+G%Z_ref))**2, CS%h2(i,j))
+          CS%h2(i,j) = min((max_frac_rough * max(G%bathyT(i,j)+G%Z_ref, 0.0))**2, CS%h2(i,j))
       endif
 
       utide = CS%tideamp(i,j)

src/tracer/MOM_tracer_Z_init.F90

@@ -138,7 +138,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va
 
       do i=is,ie ; if (G%mask2dT(i,j)*htot(i) > 0.0) then
         ! Determine the z* heights of the model interfaces.
-        dilate = (G%bathyT(i,j) + G%Z_ref) / htot(i)
+        dilate = max(G%bathyT(i,j) + G%Z_ref, 0.0) / htot(i)
         e(nz+1) = -G%bathyT(i,j) - G%Z_ref
         do k=nz,1,-1 ; e(K) = e(K+1) + dilate * h(i,j,k) ; enddo