From bd5fe0c2ca70b2226dd3bbae7f611ec364ed2063 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 14 May 2023 12:39:38 -0400 Subject: [PATCH] +(*)Use tv%SpV in MOM_sponge code Use tv%SpV convert thicknesses to vertical distances in apply_sponge when it is allocated to replace multiplication by GV%H_to_Z, thereby eliminating the dependence on GV%RHo_0 in this modue in non-Boussinesq mode. The new internal array dz_to_h is used to reduce the code duplication as a result of these changes. All answers in Boussinesq test cases are bitwise identical, but answers change in fully non-Boussinesq mode. In semi-Boussinesq mode answers are mathematically equivalent but change at roundoff unless RHO_0 is an integer power of 2. --- src/parameterizations/vertical/MOM_sponge.F90 | 76 ++++++++++++++----- 1 file changed, 58 insertions(+), 18 deletions(-) diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90 index fce1eb493d..4bdf610a24 100644 --- a/src/parameterizations/vertical/MOM_sponge.F90 +++ b/src/parameterizations/vertical/MOM_sponge.F90 @@ -347,10 +347,14 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml) ! give 0 at the surface [nondim]. real :: e(SZK_(GV)+1) ! The interface heights [Z ~> m], usually negative. + real :: dz_to_h(SZK_(GV)+1) ! Factors used to convert interface height movement + ! to thickness fluxes [H Z-1 ~> nondim or kg m-3] real :: e0 ! The height of the free surface [Z ~> m]. real :: e_str ! A nondimensional amount by which the reference ! profile must be stretched for the free surfaces ! heights in the two profiles to agree [nondim]. + real :: w_mean ! The vertical displacement of water moving upward through an + ! interface within 1 timestep [Z ~> m]. real :: w ! The thickness of water moving upward through an ! interface within 1 timestep [H ~> m or kg m-2]. real :: wm ! wm is w if w is negative and 0 otherwise [H ~> m or kg m-2]. @@ -384,9 +388,15 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml) "work properly with i-mean sponges and a bulk mixed layer.") do j=js,je ; do i=is,ie ; e_D(i,j,nz+1) = -G%bathyT(i,j) ; enddo ; enddo - do k=nz,1,-1 ; do j=js,je ; do i=is,ie - e_D(i,j,K) = e_D(i,j,K+1) + h(i,j,k)*GV%H_to_Z - enddo ; enddo ; enddo + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do k=nz,1,-1 ; do j=js,je ; do i=is,ie + e_D(i,j,K) = e_D(i,j,K+1) + GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) + enddo ; enddo ; enddo + else + do k=nz,1,-1 ; do j=js,je ; do i=is,ie + e_D(i,j,K) = e_D(i,j,K+1) + h(i,j,k)*GV%H_to_Z + enddo ; enddo ; enddo + endif do j=js,je do i=is,ie dilate(i) = (G%bathyT(i,j) + G%Z_ref) / (e_D(i,j,1) + G%bathyT(i,j)) @@ -424,20 +434,39 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml) do K=2,nz+1 ; do i=is,ie h_above(i,K) = h_above(i,K-1) + max(h(i,j,k-1)-GV%Angstrom_H, 0.0) enddo ; enddo - do K=2,nz - ! w is positive for an upward (lightward) flux of mass, resulting - ! in the downward movement of an interface. - w = damp_1pdamp * eta_mean_anom(j,K) * GV%Z_to_H - do i=is,ie + + ! In both blocks below, w is positive for an upward (lightward) flux of mass, + ! resulting in the downward movement of an interface. + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do K=2,nz + w_mean = damp_1pdamp * eta_mean_anom(j,K) + do i=is,ie + w = w_mean * 2.0*GV%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) + if (w > 0.0) then + w_int(i,j,K) = min(w, h_below(i,K)) + eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K) + else + w_int(i,j,K) = max(w, -h_above(i,K)) + ea(i,j,k) = ea(i,j,k) - w_int(i,j,K) + endif + enddo + enddo + else + do K=2,nz + w = damp_1pdamp * eta_mean_anom(j,K) * GV%Z_to_H if (w > 0.0) then - w_int(i,j,K) = min(w, h_below(i,K)) - eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K) + do i=is,ie + w_int(i,j,K) = min(w, h_below(i,K)) + eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K) + enddo else - w_int(i,j,K) = max(w, -h_above(i,K)) - ea(i,j,k) = ea(i,j,k) - w_int(i,j,K) + do i=is,ie + w_int(i,j,K) = max(w, -h_above(i,K)) + ea(i,j,k) = ea(i,j,k) - w_int(i,j,K) + enddo endif enddo - enddo + endif do k=1,nz ; do i=is,ie ea_k = max(0.0, -w_int(i,j,K)) eb_k = max(0.0, w_int(i,j,K+1)) @@ -462,9 +491,20 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml) damp = dt * CS%Iresttime_col(c) e(1) = 0.0 ; e0 = 0.0 - do K=1,nz - e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z - enddo + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do K=1,nz + e(K+1) = e(K) - GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) + enddo + dz_to_h(1) = GV%RZ_to_H / tv%SpV_avg(i,j,1) + do K=2,nz + dz_to_h(K) = 2.0*GV%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) + enddo + else + do K=1,nz + e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z + dz_to_h(K) = GV%Z_to_H + enddo + endif e_str = e(nz+1) / CS%Ref_eta(nz+1,c) if ( CS%bulkmixedlayer ) then @@ -481,7 +521,7 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml) wpb = 0.0; wb = 0.0 do k=nz,nkmb+1,-1 if (GV%Rlay(k) > Rcv_ml(i,j)) then - w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*GV%Z_to_H, & + w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*dz_to_h(K), & ((wb + h(i,j,k)) - GV%Angstrom_H)) wm = 0.5*(w-ABS(w)) do m=1,CS%fldno @@ -537,7 +577,7 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml) wpb = 0.0 wb = 0.0 do k=nz,1,-1 - w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*GV%Z_to_H, & + w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*dz_to_h(K), & ((wb + h(i,j,k)) - GV%Angstrom_H)) wm = 0.5*(w - ABS(w)) do m=1,CS%fldno