From a67ff8ea02ffb5dbb9ac339c607b774bb4a823fc Mon Sep 17 00:00:00 2001 From: He Wang Date: Wed, 24 Sep 2025 14:22:46 -0400 Subject: [PATCH] Correction on total column thickness for wetting In a number of cases, total resting column thickness is calucated as G%bathyT + G%Z_ref, which is largely correct but for wetting, i.e. G%bathyT < 0. This commit makes a correction for eight cases with this potential bug. Additionally, a correction is made in one case that G%bathyT + G%Z_ref is used as depth where negative depth does not make sense. There is no answer changes if no wetting points are used and G%Z_ref is zero. List of modules/processes affected: * MOM_barotropic * affects only surface stress when BT_NONLIN_STRESS is False. * MOM_wave_speed * h2 calculations in * subroutine internal_tides_init * subroutine int_tide_input_int * subroutine tidal_mixing_init * MOM_tracer_z_init * D_[uv] in set_viscous_BBL, which is used only when CHANNEL_DRAG = True * MOM_lateral_mixing_coeffs * MOM_MEKE --- src/ALE/MOM_hybgen_unmix.F90 | 2 +- src/core/MOM_barotropic.F90 | 11 +++++++---- src/diagnostics/MOM_wave_speed.F90 | 4 ++-- src/parameterizations/lateral/MOM_MEKE.F90 | 4 ++-- src/parameterizations/lateral/MOM_internal_tides.F90 | 2 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 6 ++++-- .../vertical/MOM_internal_tide_input.F90 | 2 +- src/parameterizations/vertical/MOM_set_viscosity.F90 | 4 ++-- src/parameterizations/vertical/MOM_tidal_mixing.F90 | 2 +- src/tracer/MOM_tracer_Z_init.F90 | 2 +- 10 files changed, 22 insertions(+), 17 deletions(-) diff --git a/src/ALE/MOM_hybgen_unmix.F90 b/src/ALE/MOM_hybgen_unmix.F90 index bb6f64c4d7..d87762b721 100644 --- a/src/ALE/MOM_hybgen_unmix.F90 +++ b/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) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index ac2e668f8e..777b45026a 100644 --- a/src/core/MOM_barotropic.F90 +++ b/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 diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 7abdab0a90..b525efd149 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/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 diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 41c98884ba..7abc24fe2d 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/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 diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 1f51cc99a9..e4a48c84c5 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/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 diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 5c9130b6f7..2a7527ce9a 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/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 diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 2384844f6e..2a7a3a2a7c 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/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 * & diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 6466b71dd5..33d38e7538 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/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 diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 9a972e6e06..faeadbb8ab 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/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) diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 8b2c12fd9b..a2f64ff5d4 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/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