From 7be08832127b4669db8644cfc261ba82163d3a5d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 7 Aug 2020 08:13:03 -0400 Subject: [PATCH 1/2] (*)Set dSV_dT and dSV_dS with unassociated fluxes Set dSV_dT and dSV_dS if present in applyBoundaryFluxesInOut, even if boundary fluxes are not associated. With this change, setting BUOY_CONFIG='NONE' and BUOY_CONFIG='zero' both work and give similar (but not identical) answers in some test cases with an ePBL boundary layer parameterization, whereas before answers were tainted with uninitialized values when BUOY_CONFIG='NONE'. All answers in the existing MOM6-examples test suite are bitwise identical, but answers can change in other cases. --- .../vertical/MOM_diabatic_aux.F90 | 15 +++++++++------ .../vertical/MOM_vert_friction.F90 | 2 +- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 91085047c9..bf2e86cb80 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -822,19 +822,18 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - ! Only apply forcing if fluxes%sw is associated. - if (.not.associated(fluxes%sw)) return - -#define _OLD_ALG_ Idt = 1.0 / dt calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS)) calculate_buoyancy = present(SkinBuoyFlux) if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0 + if (present(cTKE)) cTKE(:,:,:) = 0.0 g_Hconv2 = (US%L_to_Z**2*GV%g_Earth * GV%H_to_RZ) * GV%H_to_RZ EOSdom(:) = EOS_domain(G%HI) - if (present(cTKE)) cTKE(:,:,:) = 0.0 + ! Only apply forcing if fluxes%sw is associated. + if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return + if (calculate_buoyancy) then SurfPressure(:) = 0.0 GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 @@ -874,7 +873,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t h2d(i,k) = h(i,j,k) T2d(i,k) = tv%T(i,j,k) enddo ; enddo - if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) if (calculate_energetics) then ! The partial derivatives of specific volume with temperature and @@ -898,6 +896,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t pen_TKE_2d(:,:) = 0.0 endif + ! Nothing more is done on this j-slice if there is no buoyancy forcing. + if (.not.associated(fluxes%sw)) cycle + + if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) + ! The surface forcing is contained in the fluxes type. ! We aggregate the thermodynamic forcing for a time step into the following: ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c6a6f37b39..1a4fb58e02 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1143,12 +1143,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, 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 :: z2 ! A copy of z_i [nondim] + real :: botfn ! A function that is 1 at the bottom and small far from it [nondim] real :: topfn ! A function that is 1 at the top and small far from it [nondim] real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] logical :: do_shelf, do_OBCs integer :: i, k, is, ie, max_nk integer :: nz - real :: botfn a_cpl(:,:) = 0.0 Kv_tot(:,:) = 0.0 From feed9ba82b9ebe6805ce691e77190b2f5ba4f7ee Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 7 Aug 2020 08:13:40 -0400 Subject: [PATCH 2/2] (*)Fix an indexing bug in int_density_dz_linear Corrected a horizontal indexing bug in int_density_dz_linear that caused the ISOMIP/layer test case to fail. This bug was first introduced with PR#732 on March 8, 2018. This bug fix will change answers with a linear equation of state and the finite volume pressure gradient force, however it does not change any of the verified answers in the MOM6-examples regression suite. --- src/equation_of_state/MOM_EOS_linear.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index e3a5443840..47a2bf21b0 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -473,7 +473,7 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) + intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) do m=2,4 wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR