From df5515b292942210e34d2fb6dd3a7a82ba7237a0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 14 Aug 2023 18:38:11 -0400 Subject: [PATCH] Use GV%dZ_subroundoff Use the minimal vertical distance GV%dZ_subroundoff in 8 spots in 2 files in place of GV%H_to_Z*GV%H_subroundoff. Four internal variables were renamed for clarity of purpose in porous_widths_layer and PressureForce_Mont_Bouss. Also modified the non-Boussinesq units in a description of two integrated energy diagnostics from [H L4 T-3 ~> kg m2 s-3] to [H L4 T-3 ~> W]; these units are equivalent but the latter is more informative. Answers in all Boussinesq test cases and the existing non-Boussinesq test case in the MOM6-examples test suite are bitwise identical. --- src/core/MOM_PressureForce_Montgomery.F90 | 16 ++++++++-------- src/core/MOM_porous_barriers.F90 | 16 ++++++++-------- src/diagnostics/MOM_diagnostics.F90 | 4 ++-- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 424e9b1a32..43bbbb82c1 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -405,7 +405,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, real :: G_Rho0 ! G_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. real :: PFu_bc, PFv_bc ! The pressure gradient force due to along-layer ! compensated density gradients [L T-2 ~> m s-2] - real :: h_neglect ! A thickness that is so small it is usually lost + real :: dz_neglect ! A vertical distance that is so small it is usually lost ! in roundoff and can be neglected [Z ~> m]. logical :: use_p_atm ! If true, use the atmospheric pressure. logical :: use_EOS ! If true, density is calculated from T & S using @@ -436,7 +436,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.") endif - h_neglect = GV%H_subroundoff * GV%H_to_Z + dz_neglect = GV%dZ_subroundoff I_Rho0 = 1.0/CS%Rho0 G_Rho0 = GV%g_Earth / GV%Rho0 @@ -552,7 +552,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, !$OMP parallel do default(shared) private(h_star,PFu_bc,PFv_bc) do k=1,nz do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - h_star(i,j) = (e(i,j,K) - e(i,j,K+1)) + h_neglect + h_star(i,j) = (e(i,j,K) - e(i,j,K+1)) + dz_neglect enddo ; enddo do j=js,je ; do I=Isq,Ieq PFu_bc = -1.0*(rho_star(i+1,j,k) - rho_star(i,j,k)) * (G%IdxCu(I,j) * & @@ -636,7 +636,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) real :: Rho0xG ! g_Earth * Rho0 [R L2 Z-1 T-2 ~> kg s-2 m-2] logical :: use_EOS ! If true, density is calculated from T & S using ! an equation of state. - real :: z_neglect ! A thickness that is so small it is usually lost + real :: dz_neglect ! A vertical distance that is so small it is usually lost ! in roundoff and can be neglected [Z ~> m]. integer, dimension(2) :: EOSdom ! The computational domain for the equation of state integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k @@ -647,14 +647,14 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) Rho0xG = Rho0 * GV%g_Earth G_Rho0 = GV%g_Earth / GV%Rho0 use_EOS = associated(tv%eqn_of_state) - z_neglect = GV%H_subroundoff*GV%H_to_Z + dz_neglect = GV%dZ_subroundoff if (use_EOS) then if (present(rho_star)) then !$OMP parallel do default(shared) private(Ihtot) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect) + Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + dz_neglect) pbce(i,j,1) = GFS_scale * rho_star(i,j,1) * GV%H_to_Z enddo do k=2,nz ; do i=Isq,Ieq+1 @@ -666,7 +666,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) !$OMP parallel do default(shared) private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect) + Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + dz_neglect) press(i) = -Rho0xG*(e(i,j,1) - G%Z_ref) enddo call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, & @@ -695,7 +695,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) !$OMP parallel do default(shared) private(Ihtot) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect) + Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + dz_neglect) pbce(i,j,1) = GV%g_prime(1) * GV%H_to_Z enddo do k=2,nz ; do i=Isq,Ieq+1 diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index d73d96b242..e212581993 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -80,7 +80,7 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) logical, dimension(SZIB_(G),SZJB_(G)) :: do_I ! Booleans for calculation at u or v points ! updated while moving up layers real :: A_layer ! Integral of fractional open width from bottom to current layer [Z ~> m] - real :: h_min ! ! The minimum layer thickness [Z ~> m] + real :: dz_min ! The minimum layer thickness [Z ~> m] real :: dmask ! The depth below which porous barrier is not applied [Z ~> m] integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq @@ -100,7 +100,7 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) call calc_eta_at_uv(eta_u, eta_v, CS%eta_interp, dmask, h, tv, G, GV, US) - h_min = GV%Angstrom_H * GV%H_to_Z + dz_min = GV%Angstrom_Z ! u-points do j=js,je ; do I=Isq,Ieq ; do_I(I,j) = .False. ; enddo ; enddo @@ -125,7 +125,7 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & eta_u(I,j,K), A_layer, do_I(I,j)) - if (eta_u(I,j,K) - (eta_u(I,j,K+1)+h_min) > 0.0) then + if (eta_u(I,j,K) - (eta_u(I,j,K+1)+dz_min) > 0.0) then pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1))) else pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice @@ -157,7 +157,7 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & eta_v(i,J,K), A_layer, do_I(i,J)) - if (eta_v(i,J,K) - (eta_v(i,J,K+1)+h_min) > 0.0) then + if (eta_v(i,J,K) - (eta_v(i,J,K+1)+dz_min) > 0.0) then pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1))) else pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice @@ -286,7 +286,7 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! Layer interface heights [Z ~> m]. - real :: h_neglect ! Negligible thicknesses [Z ~> m] + real :: dz_neglect ! A negligible height difference [Z ~> m] integer :: i, j, k, nk, is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc; ie = G%iec; js = G%jsc; je = G%jec; nk = GV%ke @@ -295,7 +295,7 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) ! currently no treatment for using optional find_eta arguments if present call find_eta(h, tv, G, GV, US, eta, halo_size=1) - h_neglect = GV%H_subroundoff * GV%H_to_Z + dz_neglect = GV%dZ_subroundoff do K=1,nk+1 do j=js,je ; do I=Isq,Ieq ; eta_u(I,j,K) = dmask ; enddo ; enddo @@ -333,10 +333,10 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) case (ETA_INTERP_HARM) ! Harmonic mean do K=1,nk+1 do j=js,je ; do I=Isq,Ieq ; if (G%porous_DavgU(I,j) < dmask) then - eta_u(I,j,K) = 2.0 * (eta(i,j,K) * eta(i+1,j,K)) / (eta(i,j,K) + eta(i+1,j,K) + h_neglect) + eta_u(I,j,K) = 2.0 * (eta(i,j,K) * eta(i+1,j,K)) / (eta(i,j,K) + eta(i+1,j,K) + dz_neglect) endif ; enddo ; enddo do J=Jsq,Jeq ; do i=is,ie ; if (G%porous_DavgV(i,J) < dmask) then - eta_v(i,J,K) = 2.0 * (eta(i,j,K) * eta(i,j+1,K)) / (eta(i,j,K) + eta(i,j+1,K) + h_neglect) + eta_v(i,J,K) = 2.0 * (eta(i,j,K) * eta(i,j+1,K)) / (eta(i,j,K) + eta(i,j+1,K) + dz_neglect) endif ; enddo ; enddo enddo case default diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index c23712d8e7..253b7189e3 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -957,9 +957,9 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS real :: KE_term(SZI_(G),SZJ_(G),SZK_(GV)) ! A term in the kinetic energy budget ! [H L2 T-3 ~> m3 s-3 or W m-2] real :: KE_u(SZIB_(G),SZJ_(G)) ! The area integral of a KE term in a layer at u-points - ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] + ! [H L4 T-3 ~> m5 s-3 or W] real :: KE_v(SZI_(G),SZJB_(G)) ! The area integral of a KE term in a layer at v-points - ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] + ! [H L4 T-3 ~> m5 s-3 or W] real :: KE_h(SZI_(G),SZJ_(G)) ! A KE term contribution at tracer points ! [H L2 T-3 ~> m3 s-3 or W m-2]