Skip to content

Use GV%dZ_subroundoff #457

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marshallward merged 2 commits into from
Aug 30, 2023
Merged

Use GV%dZ_subroundoff #457

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
df5515b
Use GV%dZ_subroundoff
Hallberg-NOAA Aug 14, 2023
e85d8fc
Merge branch 'dev/gfdl' into use_dz_subroundoff
marshallward Aug 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -428,7 +428,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
Expand Down Expand Up @@ -459,7 +459,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

Expand Down Expand Up @@ -582,7 +582,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) * &
Expand Down Expand Up @@ -676,7 +676,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
Expand All @@ -687,14 +687,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
Expand All @@ -706,7 +706,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, &
Expand Down Expand Up @@ -735,7 +735,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
Expand Down
16 changes: 8 additions & 8 deletions src/core/MOM_porous_barriers.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/diagnostics/MOM_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand Down