Skip to content

Commit

Permalink
(*)Offset G%bathyT by -G%Z_ref
Browse files Browse the repository at this point in the history 
  Added code to adjust the value of G%bathyT in the ocean_grid_type by -G%Z_ref,
with changes scattered throughout the code to compensate.  This does not impact
the value of bathyT in the dyn_hor_grid.  Using a non-zero value of
REFERENCE_HEIGHT leads to mathematically equivalent solutions that differ at
round-off but are demonstrably similar for values of REFERENCE_HEIGHT ranging
from -1e4 m to 1e4 m, and in cases that do not use the less exact
..._2018_ANSWERS algorithms the solutions are qualitatively equivalent for
values ranging from -1e8 m to 1e8 m.  (For these more extreme values, the 64-bit
roundoff in the free surface height calculation starts to become physically
significant at of order 0.1 mm.)  This new capability is useful for developing
wetting and drying and for identifying algorithmic shortcomings, but because
answers are not bitwise identical for various reference heights, it is not such
a good candidate for fully automated testing.  By default, all answers are
bitwise identical and all output files are the same.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Aug 31, 2021
1 parent df97d4f commit 060c412
Show file tree
Hide file tree
Showing 30 changed files with 219 additions and 200 deletions.
4 changes: 2 additions & 2 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
if (CS%id_T_preale > 0) call post_data(CS%id_T_preale, tv%T, CS%diag)
if (CS%id_S_preale > 0) call post_data(CS%id_S_preale, tv%S, CS%diag)
if (CS%id_e_preale > 0) then
call find_eta(h, tv, G, GV, US, eta_preale)
call find_eta(h, tv, G, GV, US, eta_preale, dZref=G%Z_ref)
call post_data(CS%id_e_preale, eta_preale, CS%diag)
endif

Expand Down Expand Up @@ -1304,7 +1304,7 @@ subroutine ALE_initThicknessToCoord( CS, G, GV, h )
integer :: i, j, k

do j = G%jsd,G%jed ; do i = G%isd,G%ied
h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j) )
h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref )
enddo ; enddo

end subroutine ALE_initThicknessToCoord
Expand Down
12 changes: 6 additions & 6 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1113,7 +1113,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)
endif

! Local depth (G%bathyT is positive downward)
nominalDepth = G%bathyT(i,j)*GV%Z_to_H
nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H

! Determine water column thickness
totalThickness = 0.0
Expand Down Expand Up @@ -1203,7 +1203,7 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface )
endif

! The rest of the model defines grids integrating up from the bottom
nominalDepth = G%bathyT(i,j)*GV%Z_to_H
nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H

! Determine water column height
totalThickness = 0.0
Expand Down Expand Up @@ -1314,7 +1314,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel


! Local depth (G%bathyT is positive downward)
nominalDepth = G%bathyT(i,j)*GV%Z_to_H
nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H

! Determine total water column thickness
totalThickness = 0.0
Expand Down Expand Up @@ -1444,7 +1444,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_she
do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
if (G%mask2dT(i,j)>0.) then

nominalDepth = G%bathyT(i,j) * GV%Z_to_H
nominalDepth = (G%bathyT(i,j)+G%Z_ref) * GV%Z_to_H

if (ice_shelf) then
totalThickness = 0.0
Expand Down Expand Up @@ -1592,7 +1592,7 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS)
do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
if (G%mask2dT(i,j)>0.) then

depth = G%bathyT(i,j) * GV%Z_to_H
depth = (G%bathyT(i,j)+G%Z_ref) * GV%Z_to_H
z_col(1) = 0. ! Work downward rather than bottom up
do K=1,nz
z_col(K+1) = z_col(K) + h(i,j,k)
Expand Down Expand Up @@ -1718,7 +1718,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS )
do i = G%isc-1,G%iec+1

! Local depth
local_depth = G%bathyT(i,j)*GV%Z_to_H
local_depth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H

! Determine water column height
total_height = 0.0
Expand Down
2 changes: 1 addition & 1 deletion src/ALE/coord_adapt.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
zNext(nz+1) = zInt(i,j,nz+1)

! local depth for scaling diffusivity
depth = G%bathyT(i,j) * GV%Z_to_H
depth = (G%bathyT(i,j) + G%Z_ref) * GV%Z_to_H

! initialize del2sigma and the thickness change response to it zero
del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0
Expand Down
18 changes: 9 additions & 9 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -850,7 +850,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
! Determining the time-average sea surface height is part of the algorithm.
! This may be eta_av if Boussinesq, or need to be diagnosed if not.
CS%time_in_cycle = CS%time_in_cycle + dt
call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0)
call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0, dZref=G%Z_ref)
do j=js,je ; do i=is,ie
CS%ssh_rint(i,j) = CS%ssh_rint(i,j) + dt*ssh(i,j)
enddo ; enddo
Expand Down Expand Up @@ -2278,7 +2278,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &

ALLOC_(CS%ssh_rint(isd:ied,jsd:jed)) ; CS%ssh_rint(:,:) = 0.0
ALLOC_(CS%ave_ssh_ibc(isd:ied,jsd:jed)) ; CS%ave_ssh_ibc(:,:) = 0.0
ALLOC_(CS%eta_av_bc(isd:ied,jsd:jed)) ; CS%eta_av_bc(:,:) = 0.0
ALLOC_(CS%eta_av_bc(isd:ied,jsd:jed)) ; CS%eta_av_bc(:,:) = 0.0 ! -G%Z_ref
CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0

! Use the Wright equation of state by default, unless otherwise specified
Expand Down Expand Up @@ -2789,9 +2789,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &

if (.not.query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then
if (CS%split) then
call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0)
call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0, dZref=G%Z_ref)
else
call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta_to_m=1.0)
call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta_to_m=1.0, dZref=G%Z_ref)
endif
endif
if (CS%split) deallocate(eta)
Expand Down Expand Up @@ -2852,7 +2852,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
restart_CSp_tmp = restart_CSp
call restart_registry_lock(restart_CSp_tmp, unlocked=.true.)
allocate(z_interface(SZI_(G),SZJ_(G),SZK_(GV)+1))
call find_eta(CS%h, CS%tv, G, GV, US, z_interface, eta_to_m=1.0)
call find_eta(CS%h, CS%tv, G, GV, US, z_interface, eta_to_m=1.0, dZref=G%Z_ref)
call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, &
"Interface heights", "meter", z_grid='i')
! NOTE: write_ic=.true. routes routine to fms2 IO write_initial_conditions interface
Expand Down Expand Up @@ -3401,10 +3401,10 @@ subroutine extract_surface_state(CS, sfc_state_in)
numberOfErrors=0 ! count number of errors
do j=js,je ; do i=is,ie
if (G%mask2dT(i,j)>0.) then
localError = sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) &
localError = sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) - G%Z_ref &
.or. sfc_state%sea_lev(i,j) >= CS%bad_val_ssh_max &
.or. sfc_state%sea_lev(i,j) <= -CS%bad_val_ssh_max &
.or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) < CS%bad_val_col_thick
.or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) + G%Z_ref < CS%bad_val_col_thick
if (use_temperature) localError = localError &
.or. sfc_state%SSS(i,j)<0. &
.or. sfc_state%SSS(i,j)>=CS%bad_val_sss_max &
Expand All @@ -3420,7 +3420,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), &
'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), &
'D=',CS%US%Z_to_m*G%bathyT(i,j), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), &
'D=',CS%US%Z_to_m*(G%bathyT(i,j)+G%Z_ref), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), &
'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), &
'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), &
'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J)
Expand All @@ -3429,7 +3429,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), &
'x=',G%gridLonT(i), 'y=',G%gridLatT(j), &
'D=',CS%US%Z_to_m*G%bathyT(i,j), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), &
'D=',CS%US%Z_to_m*(G%bathyT(i,j)+G%Z_ref), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), &
'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), &
'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J)
endif
Expand Down
42 changes: 21 additions & 21 deletions src/core/MOM_PressureForce_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,10 +90,9 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
!! or atmosphere-ocean interface [R L2 T-2 ~> Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure
!! anomaly in each layer due to eta anomalies
!! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
!! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
!! contributions or compressibility compensation.
!! [L2 T-2 H-1 ~> m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to
!! calculate PFu and PFv [H ~> kg m-2].
! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
Expand Down Expand Up @@ -301,7 +300,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
! Find and add the tidal geopotential anomaly.
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth
SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref
enddo ; enddo
call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
!$OMP parallel do default(shared)
Expand Down Expand Up @@ -430,15 +429,16 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
!! or atmosphere-ocean interface [R L2 T-2 ~> Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure
!! anomaly in each layer due to eta anomalies
!! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
!! calculate PFu and PFv [H ~> m or kg m-2], with any
!! tidal contributions or compressibility compensation.
!! [L2 T-2 H-1 ~> m s-2].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The sea-surface height used to
!! calculate PFu and PFv [H ~> m], with any
!! tidal contributions.
! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height in depth units [Z ~> m].
real, dimension(SZI_(G),SZJ_(G)) :: &
e_tidal, & ! The bottom geopotential anomaly due to tidal forces from
! astronomical sources and self-attraction and loading [Z ~> m].
SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
dM ! The barotropic adjustment to the Montgomery potential to
! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
real, dimension(SZI_(G)) :: &
Expand All @@ -451,7 +451,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
dpa, & ! The change in pressure anomaly between the top and bottom
! of a layer [R L2 T-2 ~> Pa].
intz_dpa ! The vertical integral in depth of the pressure anomaly less the
! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa or kg m-2 Pa].
! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa].
real, dimension(SZIB_(G),SZJ_(G)) :: &
intx_pa, & ! The zonal integral of the pressure anomaly along the interface
! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
Expand Down Expand Up @@ -485,7 +485,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
real :: Tl(5) ! copy and T in local stencil [degC]
real :: mn_T ! mean of T in local stencil [degC]
real :: mn_T2 ! mean of T**2 in local stencil [degC]
real :: mn_T2 ! mean of T**2 in local stencil [degC2]
real :: hl(5) ! Copy of local stencil of H [H ~> m]
real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
real, parameter :: C1_6 = 1.0/6.0
Expand Down Expand Up @@ -565,13 +565,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
e(i,j,1) = -G%bathyT(i,j)
SSH(i,j) = -G%bathyT(i,j) - G%Z_ref
enddo
do k=1,nz ; do i=Isq,Ieq+1
e(i,j,1) = e(i,j,1) + h(i,j,k)*GV%H_to_Z
SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
enddo ; enddo
enddo
call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
endif

! Here layer interface heights, e, are calculated.
Expand Down Expand Up @@ -637,13 +637,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
tv%eqn_of_state, EOSdom)
endif
do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1)
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * (e(i,j,1) - G%Z_ref)
enddo
enddo
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * GV%Rlay(1)) * e(i,j,1)
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * GV%Rlay(1)) * (e(i,j,1) - G%Z_ref)
enddo ; enddo
endif
endif
Expand All @@ -667,12 +667,12 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
if (use_p_atm) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pa(i,j) = (rho_ref*GV%g_Earth)*e(i,j,1) + p_atm(i,j)
pa(i,j) = (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref) + p_atm(i,j)
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pa(i,j) = (rho_ref*GV%g_Earth)*e(i,j,1)
pa(i,j) = (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
enddo ; enddo
endif
!$OMP parallel do default(shared)
Expand Down Expand Up @@ -700,17 +700,17 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
useMassWghtInterp=CS%useMassWghtInterp, &
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom)
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
useMassWghtInterp=CS%useMassWghtInterp)
useMassWghtInterp=CS%useMassWghtInterp, Z_0p=G%Z_ref)
endif
else
call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa, &
intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp)
intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, Z_0p=G%Z_ref)
endif
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down
19 changes: 10 additions & 9 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
optional, intent(out) :: pbce !< The baroclinic pressure anomaly in
!! each layer due to free surface height anomalies,
!! [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> kg m-1].

real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to calculate
!! PFu and PFv [H ~> kg m-2].
! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
M, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2].
Expand All @@ -104,7 +104,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
e_tidal, & ! Bottom geopotential anomaly due to tidal forces from
! astronomical sources and self-attraction and loading [Z ~> m].
geopot_bot ! Bottom geopotential relative to time-mean sea level,
geopot_bot ! Bottom geopotential relative to a temporally fixed reference value,
! including any tidal contributions [L2 T-2 ~> m2 s-2].
real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
! density [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
Expand Down Expand Up @@ -183,7 +183,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
! of self-attraction and loading.
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
SSH(i,j) = -G%bathyT(i,j)
SSH(i,j) = -G%bathyT(i,j) - G%Z_ref
enddo ; enddo
if (use_EOS) then
!$OMP parallel do default(shared)
Expand Down Expand Up @@ -393,6 +393,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
! the deepest variable density near-surface layer [R ~> kg m-3].
real :: h_star(SZI_(G),SZJ_(G)) ! Layer thickness after compensation
! for compressibility [Z ~> m].
real :: SSH(SZI_(G),SZJ_(G)) ! The sea surface height anomaly, in depth units [Z ~> m].
real :: e_tidal(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to tidal
! forces from astronomical sources and self-
! attraction and loading, in depth units [Z ~> m].
Expand Down Expand Up @@ -444,12 +445,12 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
! barotropic tides.
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1
do i=Isq,Ieq+1 ; e(i,j,1) = -G%bathyT(i,j) ; enddo
do i=Isq,Ieq+1 ; SSH(i,j) = -G%bathyT(i,j) - G%Z_ref ; enddo
do k=1,nz ; do i=Isq,Ieq+1
e(i,j,1) = e(i,j,1) + h(i,j,k)*GV%H_to_Z
SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
enddo ; enddo
enddo
call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
endif

! Here layer interface heights, e, are calculated.
Expand Down Expand Up @@ -664,7 +665,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
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)
press(i) = -Rho0xG*e(i,j,1)
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, &
tv%eqn_of_state, EOSdom)
Expand All @@ -673,7 +674,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
enddo
do k=2,nz
do i=Isq,Ieq+1
press(i) = -Rho0xG*e(i,j,K)
press(i) = -Rho0xG*(e(i,j,K) - G%Z_ref)
T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
enddo
Expand Down
Loading

0 comments on commit 060c412

Please sign in to comment.