MOM_PressureForce_FV.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.

2! See the LICENSE file for licensing information.

3! SPDX-License-Identifier: Apache-2.0

4

5!> Finite volume pressure gradient (integrated by quadrature or analytically)

6module mom_pressureforce_fv

7

8use mom_debugging, only : hchksum, uvchksum

9use mom_diag_mediator, only : post_data, register_diag_field

10use mom_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type

11use mom_error_handler, only : mom_error, fatal, warning, is_root_pe

12use mom_file_parser, only : get_param, log_param, log_version, param_file_type

13use mom_grid, only : ocean_grid_type

14use mom_pressureforce_mont, only : set_pbce_bouss, set_pbce_nonbouss

15use mom_self_attr_load, only : calc_sal, sal_cs

16use mom_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_cs

17use mom_tidal_forcing, only : calc_tidal_forcing_legacy

18use mom_unit_scaling, only : unit_scale_type

19use mom_variables, only : thermo_var_ptrs, accel_diag_ptrs

20use mom_verticalgrid, only : verticalgrid_type

21use mom_eos, only : calculate_density, calculate_spec_vol, eos_domain

22use mom_density_integrals, only : int_density_dz, int_specific_vol_dp

23use mom_density_integrals, only : int_density_dz_generic_plm, int_density_dz_generic_ppm

24use mom_density_integrals, only : int_spec_vol_dp_generic_plm

25use mom_density_integrals, only : int_density_dz_generic_pcm, int_spec_vol_dp_generic_pcm

26use mom_density_integrals, only : diagnose_mass_weight_z, diagnose_mass_weight_p

27use mom_ale, only : ts_plm_edge_values, ts_ppm_edge_values, ts_plm_wls_edge_values, ale_cs

28

29implicit none ; private

30

31#include <MOM_memory.h>

32

33public pressureforce_fv_init

34public pressureforce_fv_bouss, pressureforce_fv_nonbouss

35

36! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional

37! consistency testing. These are noted in comments with units like Z, H, L, and T, along with

38! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units

39! vary with the Boussinesq approximation, the Boussinesq variant is given first.

40

41!> Finite volume pressure gradient control structure

42type, public :: pressureforce_fv_cs ; private

43 logical :: initialized = .false. !< True if this control structure has been initialized.

44 logical :: calculate_sal = .false. !< If true, calculate self-attraction and loading.

45 logical :: sal_use_bpa = .false. !< If true, use bottom pressure anomaly instead of SSH

46 !! to calculate SAL.

47 logical :: tides = .false. !< If true, apply tidal momentum forcing.

48 real :: rho_ref !< The reference density that is subtracted off when calculating pressure

49 !! gradient forces [R ~> kg m-3].

50 logical :: rho_ref_bug !< If true, recover a bug that mixes GV%Rho0 and CS%rho_ref in Boussinesq mode.

51 real :: gfs_scale !< A scaling of the surface pressure gradients to

52 !! allow the use of a reduced gravity model [nondim].

53 type(time_type), pointer :: time !< A pointer to the ocean model's clock.

54 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the

55 !! timing of diagnostic output.

56 integer :: masswghtinterp !< A flag indicating whether and how to use mass weighting in T/S interpolation

57 logical :: correction_intxpa !< If true, apply a correction to the value of intxpa at a selected

58 !! interface under ice, using matching at the end values along with a

59 !! 5-point quadrature integral of the hydrostatic pressure or height

60 !! changes along that interface. The selected interface is either at the

61 !! ocean's surface or in the interior, depending on reset_intxpa_integral.

62 logical :: reset_intxpa_integral !< If true and the surface displacement between adjacent cells

63 !! exceeds the vertical grid spacing, reset intxpa at the interface below

64 !! a trusted interior cell. (This often applies in ice shelf cavities.)

65 logical :: masswghtinterpvanonly !< If true, don't do mass weighting of T/S interpolation unless vanished

66 logical :: reset_intxpa_flattest !< If true, use flattest interface rather than top for reset integral

67 !! in cases where no best nonvanished interface

68 real :: h_nonvanished !< A minimal layer thickness that indicates that a layer is thick enough

69 !! to usefully reestimate the pressure integral across the interface

70 !! below it [H ~> m or kg m-2]

71 logical :: use_inaccurate_pgf_rho_anom !< If true, uses the older and less accurate

72 !! method to calculate density anomalies, as used prior to

73 !! March 2018.

74 logical :: boundary_extrap !< Indicate whether high-order boundary

75 !! extrapolation should be used within boundary cells

76

77 logical :: reconstruct !< If true, polynomial profiles of T & S will be

78 !! reconstructed and used in the integrals for the

79 !! finite volume pressure gradient calculation.

80 !! The default depends on whether regridding is being used.

81

82 integer :: recon_scheme !< Order of the polynomial of the reconstruction of T & S

83 !! for the finite volume pressure gradient calculation.

84 !! By the default (1) is for a piecewise linear method

85

86 logical :: debug !< If true, write verbose checksums for debugging purposes.

87 logical :: use_ssh_in_z0p !< If true, adjust the height at which the pressure used in the

88 !! equation of state is 0 to account for the displacement of the sea

89 !! surface including adjustments for atmospheric or sea-ice pressure.

90 logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF

91 logical :: bq_sal_tides = .false. !< If true, use an alternative method for SAL and tides

92 !! in Boussinesq mode

93 integer :: tides_answer_date = 99991231 !< Recover old answers with tides

94 integer :: id_e_tide = -1 !< Diagnostic identifier

95 integer :: id_e_tidal_eq = -1 !< Diagnostic identifier

96 integer :: id_e_tidal_sal = -1 !< Diagnostic identifier

97 integer :: id_e_sal = -1 !< Diagnostic identifier

98 integer :: id_rho_pgf = -1 !< Diagnostic identifier

99 integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier

100 integer :: id_p_stanley = -1 !< Diagnostic identifier

101 integer :: id_masswt_u = -1 !< Diagnostic identifier

102 integer :: id_masswt_v = -1 !< Diagnostic identifier

103 integer :: id_sal_u = -1 !< Diagnostic identifier

104 integer :: id_sal_v = -1 !< Diagnostic identifier

105 integer :: id_tides_u = -1 !< Diagnostic identifier

106 integer :: id_tides_v = -1 !< Diagnostic identifier

107 type(sal_cs), pointer :: sal_csp => null() !< SAL control structure

108 type(tidal_forcing_cs), pointer :: tides_csp => null() !< Tides control structure

109end type pressureforce_fv_cs

110

111contains

112

113!> \brief Non-Boussinesq analytically-integrated finite volume form of pressure gradient

114!!

115!! Determines the acceleration due to hydrostatic pressure forces, using

116!! the analytic finite volume form of the Pressure gradient, and does not

117!! make the Boussinesq approximation.

118!!

119!! To work, the following fields must be set outside of the usual (is:ie,js:je)

120!! range before this subroutine is called:

121!! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).

122subroutine pressureforce_fv_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp, p_atm, pbce, eta)

123 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

124 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure

125 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

126 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> kg m-2]

127 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables

128 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: pfu !< Zonal acceleration [L T-2 ~> m s-2]

129 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: pfv !< Meridional acceleration [L T-2 ~> m s-2]

130 type(pressureforce_fv_cs), intent(in) :: cs !< Finite volume PGF control structure

131 type(ale_cs), pointer :: ale_csp !< ALE control structure

132 type(accel_diag_ptrs), pointer :: adp !< Acceleration diagnostic pointers

133 real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean

134 !! or atmosphere-ocean interface [R L2 T-2 ~> Pa].

135 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure

136 !! anomaly in each layer due to eta anomalies

137 !! [L2 T-2 H-1 ~> m4 s-2 kg-1].

138 real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to

139 !! calculate PFu and PFv [H ~> kg m-2].

140 ! Local variables

141 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa].

142 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &

143 t_tmp, & ! Temporary array of temperatures where layers that are lighter

144 ! than the mixed layer have the mixed layer's properties [C ~> degC].

145 s_tmp ! Temporary array of salinities where layers that are lighter

146 ! than the mixed layer have the mixed layer's properties [S ~> ppt].

147 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

148 s_t, s_b, & ! Top and bottom edge values for linear reconstructions

149 ! of salinity within each layer [S ~> ppt].

150 t_t, t_b ! Top and bottom edge values for linear reconstructions

151 ! of temperature within each layer [C ~> degC].

152 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

153 dza, & ! The change in geopotential anomaly between the top and bottom

154 ! of a layer [L2 T-2 ~> m2 s-2].

155 intp_dza ! The vertical integral in depth of the pressure anomaly less

156 ! the pressure anomaly at the top of the layer [R L4 T-4 ~> Pa m2 s-2].

157 real, dimension(SZI_(G),SZJ_(G)) :: &

158 dp, & ! The (positive) change in pressure across a layer [R L2 T-2 ~> Pa].

159 ssh, & ! Sea surfae height anomaly for self-attraction and loading. Used if

160 ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m].

161 pbot, & ! Total bottom pressure for self-attraction and loading. Used if

162 ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa].

163 e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].

164 e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m].

165 e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading

166 ! specific to tides [Z ~> m].

167 e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing, used for recovering

168 ! old answers only [Z ~> m].

169 dm ! The barotropic adjustment to the Montgomery potential to

170 ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].

171 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

172 za ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the

173 ! interfaces [L2 T-2 ~> m2 s-2].

174

175 real, dimension(SZI_(G)) :: rho_cv_bl ! The coordinate potential density in the deepest variable

176 ! density near-surface layer [R ~> kg m-3].

177 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &

178 intx_za ! The zonal integral of the geopotential anomaly along the

179 ! interfaces, divided by the grid spacing [L2 T-2 ~> m2 s-2].

180 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &

181 intx_dza ! The change in intx_za through a layer [L2 T-2 ~> m2 s-2].

182 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &

183 inty_za ! The meridional integral of the geopotential anomaly along the

184 ! interfaces, divided by the grid spacing [L2 T-2 ~> m2 s-2].

185 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &

186 inty_dza ! The change in inty_za through a layer [L2 T-2 ~> m2 s-2].

187 real, dimension(SZI_(G),SZJ_(G)) :: &

188 t_top, & ! Temperature of top layer used with correction_intxpa [C ~> degC]

189 s_top ! Salinity of top layer used with correction_intxpa [S ~> ppt]

190 real, dimension(SZIB_(G),SZJ_(G)) :: &

191 intx_za_cor ! Correction for curvature in intx_za [L2 T-2 ~> m2 s-2]

192 real, dimension(SZI_(G),SZJB_(G)) :: &

193 inty_za_cor ! Correction for curvature in inty_za [L2 T-2 ~> m2 s-2]

194

195 ! These variables are used with reset_intxpa_integral. The values are taken from different

196 ! interfaces as a function of position.

197 real, dimension(SZIB_(G),SZJ_(G)) :: &

198 t_int_w, t_int_e, & ! Temperatures on the reference interface to the east and west of a u-point [C ~> degC]

199 s_int_w, s_int_e, & ! Salinities on the reference interface to the east and west of a u-point [S ~> ppt]

200 p_int_w, p_int_e, & ! Pressures on the reference interface to the east and west of a u-point [R L2 T-2 ~> Pa]

201 intx_za_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface

202 ! from the value that would be obtained from assuming that pressure varies

203 ! linearly with depth along that interface [R L2 T-2 ~> Pa].

204 dp_int_x, & ! The change in x in pressure along the reference interface [R L2 T-2 ~> Pa]

205 intx_za_cor_ri ! The correction to intx_za based on the reference interface calculations [L2 T-2 ~> m2 s-2]

206 real, dimension(SZI_(G),SZJB_(G)) :: &

207 t_int_s, t_int_n, & ! Temperatures on the reference interface to the north and south of a v-point [C ~> degC]

208 s_int_s, s_int_n, & ! Salinities on the reference interface to the north and south of a v-point [S ~> ppt]

209 p_int_s, p_int_n, & ! Pressures on the reference interface to the north and south of a v-point [R L2 T-2 ~> Pa]

210 inty_za_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface

211 ! from the value that would be obtained from assuming that pressure varies

212 ! linearly with depth along that interface [L2 T-2 ~> m2 s-2].

213 dp_int_y, & ! The change in y in geopotenial height along the reference interface [R L2 T-2 ~> Pa]

214 inty_za_cor_ri ! The correction to inty_za based on the reference interface calculations [L2 T-2 ~> m2 s-2]

215 logical, dimension(SZIB_(G),SZJ_(G)) :: &

216 seek_x_cor ! If true, try to find a u-point interface that would provide a better estimate

217 ! of the curvature terms in the intx_pa.

218 logical, dimension(SZI_(G),SZJB_(G)) :: &

219 seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate

220 ! of the curvature terms in the inty_pa.

221 real, dimension(SZIB_(G),SZJ_(G)) :: &

222 delta_p_x ! If using flattest interface for reset integral, store x interface

223 ! differences [R L2 T-2 ~> Pa]

224 real, dimension(SZI_(G),SZJB_(G)) :: &

225 delta_p_y ! If using flattest interface for reset integral, store y interface

226 ! differences [R L2 T-2 ~> Pa]

227 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &

228 masswt_u ! The fractional mass weighting at a u-point [nondim].

229 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &

230 masswt_v ! The fractional mass weighting at a v-point [nondim].

231 real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate

232 ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).

233 real :: dp_sfc ! The change in surface pressure between adjacent cells [R L2 T-2 ~> Pa]

234

235 real :: dp_neglect ! A thickness that is so small it is usually lost

236 ! in roundoff and can be neglected [R L2 T-2 ~> Pa].

237 real :: p_nonvanished ! nonvanshed pressure [R L2 T-2 ~> Pa]

238 real :: i_gearth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1]

239 real :: alpha_anom ! The in-situ specific volume, averaged over a

240 ! layer, less alpha_ref [R-1 ~> m3 kg-1].

241 logical :: use_p_atm ! If true, use the atmospheric pressure.

242 logical :: use_ale ! If true, use an ALE pressure reconstruction.

243 logical :: use_eos ! If true, density is calculated from T & S using an equation of state.

244 logical :: do_more_k ! If true, there are still points where a flatter interface remains to be found.

245 type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.

246

247 real :: alpha_ref ! A reference specific volume [R-1 ~> m3 kg-1] that is used

248 ! to reduce the impact of truncation errors.

249 real :: rho_in_situ(szi_(g)) ! The in situ density [R ~> kg m-3].

250 real :: pa_to_h ! A factor to convert from Pa to the thickness units (H)

251 ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1].

252 real :: h_to_rl2_t2 ! A factor to convert from thickness units (H) to pressure

253 ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].

254 real :: t5(5) ! Temperatures and salinities at five quadrature points [C ~> degC]

255 real :: s5(5) ! Salinities at five quadrature points [S ~> ppt]

256 real :: p5(5) ! Pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa]

257 real :: spv5(5) ! Specific volume anomalies at five quadrature points [R-1 ~> m3 kg-1]

258 real :: wt_r ! A weighting factor [nondim]

259

260 ! real :: oneatm ! 1 standard atmosphere of pressure in [R L2 T-2 ~> Pa]

261 real, parameter :: c1_6 = 1.0/6.0 ! [nondim]

262 real, parameter :: c1_90 = 1.0/90.0 ! A rational constant [nondim]

263 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb

264 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state

265 integer, dimension(2) :: eosdom_u ! The i-computational domain for the equation of state at u-velocity points

266 integer, dimension(2) :: eosdom_v ! The i-computational domain for the equation of state at v-velocity points

267 integer :: i, j, k, m

268

269 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

270 nkmb=gv%nk_rho_varies

271 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

272 eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)

273 eosdom_u(1) = isq - (g%IsdB-1) ; eosdom_u(2) = ieq - (g%IsdB-1)

274 eosdom_v(1) = is - (g%isd-1) ; eosdom_v(2) = ie - (g%isd-1)

275

276 if (.not.cs%initialized) call mom_error(fatal, &

277 "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")

278

279 if (cs%use_stanley_pgf) call mom_error(fatal, &

280 "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet "//&

281 "implemented in non-Boussinesq mode.")

282

283 use_p_atm = associated(p_atm)

284 use_eos = associated(tv%eqn_of_state)

285 use_ale = .false.

286 if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos

287

288 h_to_rl2_t2 = gv%g_Earth*gv%H_to_RZ

289 dp_neglect = gv%g_Earth*gv%H_to_RZ * gv%H_subroundoff

290 alpha_ref = 1.0 / cs%rho_ref

291 i_gearth = 1.0 / gv%g_Earth

292 p_nonvanished = gv%g_Earth*gv%H_to_RZ*cs%h_nonvanished

293

294 if ((cs%id_MassWt_u > 0) .or. (cs%id_MassWt_v > 0)) then

295 masswt_u(:,:,:) = 0.0 ; masswt_v(:,:,:) = 0.0

296 endif

297

298 if (use_p_atm) then

299 !$OMP parallel do default(shared)

300 do j=jsq,jeq+1 ; do i=isq,ieq+1

301 p(i,j,1) = p_atm(i,j)

302 enddo ; enddo

303 else

304 ! oneatm = 101325.0 * US%Pa_to_RL2_T2 ! 1 atm scaled to [R L2 T-2 ~> Pa]

305 !$OMP parallel do default(shared)

306 do j=jsq,jeq+1 ; do i=isq,ieq+1

307 p(i,j,1) = 0.0 ! or oneatm

308 enddo ; enddo

309 endif

310 !$OMP parallel do default(shared)

311 do j=jsq,jeq+1 ; do k=2,nz+1 ; do i=isq,ieq+1

312 p(i,j,k) = p(i,j,k-1) + h_to_rl2_t2 * h(i,j,k-1)

313 enddo ; enddo ; enddo

314

315 if (use_eos) then

316 ! With a bulk mixed layer, replace the T & S of any layers that are

317 ! lighter than the buffer layer with the properties of the buffer

318 ! layer. These layers will be massless anyway, and it avoids any

319 ! formal calculations with hydrostatically unstable profiles.

320 if (nkmb>0) then

321 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp

322 tv_tmp%eqn_of_state => tv%eqn_of_state

323 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo

324 !$OMP parallel do default(shared) private(Rho_cv_BL)

325 do j=jsq,jeq+1

326 do k=1,nkmb ; do i=isq,ieq+1

327 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)

328 enddo ; enddo

329 call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, rho_cv_bl(:), &

330 tv%eqn_of_state, eosdom)

331 do k=nkmb+1,nz ; do i=isq,ieq+1

332 if (gv%Rlay(k) < rho_cv_bl(i)) then

333 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)

334 else

335 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)

336 endif

337 enddo ; enddo

338 enddo

339 else

340 tv_tmp%T => tv%T ; tv_tmp%S => tv%S

341 tv_tmp%eqn_of_state => tv%eqn_of_state

342 endif

343 endif

344

345 ! If regridding is activated, do a linear reconstruction of salinity

346 ! and temperature across each layer. The subscripts 't' and 'b' refer

347 ! to top and bottom values within each layer (these are the only degrees

348 ! of freedom needed to know the linear profile).

349 if ( use_ale .and. (cs%Recon_Scheme == 1) ) then

350 call ts_plm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)

351 elseif ( use_ale .and. (cs%Recon_Scheme == 2) ) then

352 call ts_ppm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)

353 elseif ( use_ale .and. (cs%Recon_Scheme == 3) ) then

354 call ts_plm_wls_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h)

355 elseif (cs%reset_intxpa_integral) then

356 do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1

357 t_b(i,j,k) = tv%T(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)

358 enddo ; enddo ; enddo

359 endif

360

361 !$OMP parallel do default(shared) private(alpha_anom,dp)

362 do k=1,nz

363 ! Calculate 4 integrals through the layer that are required in the

364 ! subsequent calculation.

365 if (use_eos) then

366 if ( use_ale .and. cs%Recon_Scheme > 0 ) then

367 if ( cs%Recon_Scheme == 1 .or. cs%Recon_Scheme == 3 ) then

368 call int_spec_vol_dp_generic_plm( t_t(:,:,k), t_b(:,:,k), s_t(:,:,k), s_b(:,:,k), &

369 p(:,:,k), p(:,:,k+1), alpha_ref, dp_neglect, p(:,:,nz+1), g%HI, &

370 tv%eqn_of_state, us, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), &

371 p_surf=p(:,:,1), masswghtinterp=cs%MassWghtInterp, &

372 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, p_nv=p_nonvanished)

373 elseif ( cs%Recon_Scheme == 2 ) then

374 call mom_error(fatal, "PressureForce_FV_nonBouss: "//&

375 "int_spec_vol_dp_generic_ppm does not exist yet.")

376 ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &

377 ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), &

378 ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &

379 ! intx_dza(:,:,k), inty_dza(:,:,k), P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp)

380 endif

381 else

382 call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,k), &

383 p(:,:,k+1), alpha_ref, g%HI, tv%eqn_of_state, &

384 us, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &

385 inty_dza(:,:,k), bathyp=p(:,:,nz+1), p_surf=p(:,:,1), dp_tiny=dp_neglect, &

386 masswghtinterp=cs%MassWghtInterp, &

387 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, p_nv=p_nonvanished)

388 endif

389 if ((cs%id_MassWt_u > 0) .or. (cs%id_MassWt_v > 0)) &

390 call diagnose_mass_weight_p(p(:,:,k), p(:,:,k+1), p(:,:,nz+1), p(:,:,1), dp_neglect, cs%MassWghtInterp, &

391 g%HI, masswt_u(:,:,k), masswt_v(:,:,k), &

392 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, p_nv=p_nonvanished)

393 else

394 alpha_anom = 1.0 / gv%Rlay(k) - alpha_ref

395 do j=jsq,jeq+1 ; do i=isq,ieq+1

396 dp(i,j) = h_to_rl2_t2 * h(i,j,k)

397 dza(i,j,k) = alpha_anom * dp(i,j)

398 intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2

399 enddo ; enddo

400 do j=js,je ; do i=isq,ieq

401 intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))

402 enddo ; enddo

403 do j=jsq,jeq ; do i=is,ie

404 inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))

405 enddo ; enddo

406 endif

407 enddo

408

409 ! The bottom geopotential anomaly is calculated first so that the increments

410 ! to the geopotential anomalies can be reused. Alternately, the surface

411 ! geopotential could be calculated directly with separate calls to

412 ! int_specific_vol_dp with alpha_ref=0, and the anomalies used going

413 ! downward, which would relieve the need for dza, intp_dza, intx_dza, and

414 ! inty_dza to be 3-D arrays.

415

416 ! Sum vertically to determine the surface geopotential anomaly.

417 !$OMP parallel do default(shared)

418 do j=jsq,jeq+1

419 do i=isq,ieq+1

420 za(i,j,nz+1) = alpha_ref*p(i,j,nz+1) - gv%g_Earth*g%bathyT(i,j)

421 enddo

422 do k=nz,1,-1 ; do i=isq,ieq+1

423 za(i,j,k) = za(i,j,k+1) + dza(i,j,k)

424 enddo ; enddo

425 enddo

426

427 ! Calculate and add self-attraction and loading (SAL) geopotential height anomaly to interface height.

428 if (cs%calculate_SAL) then

429 if (cs%sal_use_bpa) then

430 !$OMP parallel do default(shared)

431 do j=jsq,jeq+1 ; do i=isq,ieq+1

432 pbot(i,j) = p(i,j,nz+1)

433 enddo ; enddo

434 call calc_sal(pbot, e_sal, g, cs%SAL_CSp, tmp_scale=us%Z_to_m)

435 else

436 !$OMP parallel do default(shared)

437 do j=jsq,jeq+1 ; do i=isq,ieq+1

438 ssh(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * i_gearth - g%Z_ref

439 ! Remove above sea level topography at floodable cells

440 ssh(i,j) = ssh(i,j) - max(-g%bathyT(i,j)-g%meanSL(i,j), 0.0)

441 enddo ; enddo

442 call calc_sal(ssh, e_sal, g, cs%SAL_CSp, tmp_scale=us%Z_to_m)

443 endif

444

445 ! This gives new answers after the change of separating SAL from tidal forcing module.

446 if ((cs%tides_answer_date>20230630) .or. (.not.gv%semi_Boussinesq) .or. (.not.cs%tides)) then

447 !$OMP parallel do default(shared)

448 do j=jsq,jeq+1 ; do i=isq,ieq+1

449 za(i,j,1) = za(i,j,1) - gv%g_Earth * e_sal(i,j)

450 enddo ; enddo

451 endif

452 endif

453

454 ! Calculate and add tidal geopotential height anomaly to interface height.

455 if (cs%tides) then

456 if ((cs%tides_answer_date>20230630) .or. (.not.gv%semi_Boussinesq)) then

457 call calc_tidal_forcing(cs%Time, e_tidal_eq, e_tidal_sal, g, us, cs%tides_CSp)

458 !$OMP parallel do default(shared)

459 do j=jsq,jeq+1 ; do i=isq,ieq+1

460 za(i,j,1) = za(i,j,1) - gv%g_Earth * (e_tidal_eq(i,j) + e_tidal_sal(i,j))

461 enddo ; enddo

462 else ! This block recreates older answers with tides.

463 if (.not.cs%calculate_SAL) e_sal(:,:) = 0.0

464 call calc_tidal_forcing_legacy(cs%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, &

465 g, us, cs%tides_CSp)

466 !$OMP parallel do default(shared)

467 do j=jsq,jeq+1 ; do i=isq,ieq+1

468 za(i,j,1) = za(i,j,1) - gv%g_Earth * e_sal_and_tide(i,j)

469 enddo ; enddo

470 endif

471 endif

472

473 ! Find the height anomalies at the interfaces. If there are no tides and no SAL,

474 ! there is no need to correct za, but omitting this changes answers at roundoff.

475 do k=1,nz

476 !$OMP parallel do default(shared)

477 do j=jsq,jeq+1 ; do i=isq,ieq+1

478 za(i,j,k+1) = za(i,j,k) - dza(i,j,k)

479 enddo ; enddo

480 enddo

481

482 if (cs%debug) then

483 call hchksum(za, "Pre-correction za", g%HI, haloshift=1, unscale=us%L_T_to_m_s**2)

484 call hchksum(p, "Pre-correction p", g%HI, haloshift=1, unscale=us%RL2_T2_to_Pa)

485 endif

486

487 ! With an ice-shelf or icebergs, this linearity condition might need to be applied

488 ! to a sub-surface interface.

489 if (cs%correction_intxpa .or. cs%reset_intxpa_integral) then

490 ! Determine surface temperature and salinity for use in the pressure gradient corrections

491 if (use_ale .and. (cs%Recon_Scheme > 0)) then

492 do j=jsq,jeq+1 ; do i=isq,ieq+1

493 t_top(i,j) = t_t(i,j,1) ; s_top(i,j) = s_t(i,j,1)

494 enddo ; enddo

495 else

496 do j=jsq,jeq+1 ; do i=isq,ieq+1

497 t_top(i,j) = tv%T(i,j,1) ; s_top(i,j) = tv%S(i,j,1)

498 enddo ; enddo

499 endif

500 endif

501

502 if (cs%correction_intxpa) then

503 ! This version makes a 5 point quadrature correction for hydrostatic variations in surface

504 ! pressure under ice.

505 !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5)

506 do j=js,je ; do i=isq,ieq

507 intx_za_cor(i,j) = 0.0

508 dp_sfc = (p(i+1,j,1) - p(i,j,1))

509 ! If the changes in pressure and height anomaly were explicable by just a hydrostatic balance,

510 ! the implied specific volume would be SpV_implied = alpha_ref - (dza_x / dp_x)

511 if (dp_sfc * (alpha_ref*dp_sfc - (za(i+1,j,1)-za(i,j,1))) > 0.0) then

512 t5(1) = t_top(i,j) ; t5(5) = t_top(i+1,j)

513 s5(1) = s_top(i,j) ; s5(5) = s_top(i+1,j)

514 p5(1) = p(i,j,1) ; p5(5) = p(i+1,j,1)

515 do m=2,4

516 wt_r = 0.25*real(m-1)

517 t5(m) = t5(1) + (t5(5)-t5(1))*wt_r

518 s5(m) = s5(1) + (s5(5)-s5(1))*wt_r

519 p5(m) = p5(1) + (p5(5)-p5(1))*wt_r

520 enddo !m

521 call calculate_spec_vol(t5, s5, p5, spv5, tv%eqn_of_state, spv_ref=alpha_ref)

522 ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression.

523 intx_za_cor(i,j) = c1_90 * (4.75*(spv5(5)-spv5(1)) + 5.5*(spv5(4)-spv5(2))) * dp_sfc

524 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12

525 endif

526 intx_za(i,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) + intx_za_cor(i,j)

527 enddo ; enddo

528 !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5)

529 do j=jsq,jeq ; do i=is,ie

530 inty_za_cor(i,j) = 0.0

531 dp_sfc = (p(i,j+1,1) - p(i,j,1))

532 if (dp_sfc * (alpha_ref*dp_sfc - (za(i,j+1,1)-za(i,j,1))) > 0.0) then

533 ! The pressure/depth relationship has a positive implied specific volume.

534 t5(1) = t_top(i,j) ; t5(5) = t_top(i,j+1)

535 s5(1) = s_top(i,j) ; s5(5) = s_top(i,j+1)

536 p5(1) = p(i,j,1) ; p5(5) = p(i,j+1,1)

537 do m=2,4

538 wt_r = 0.25*real(m-1)

539 t5(m) = t5(1) + (t5(5)-t5(1))*wt_r

540 s5(m) = s5(1) + (s5(5)-s5(1))*wt_r

541 p5(m) = p5(1) + (p5(5)-p5(1))*wt_r

542 enddo !m

543 call calculate_spec_vol(t5, s5, p5, spv5, tv%eqn_of_state, spv_ref=alpha_ref)

544 ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression.

545 inty_za_cor(i,j) = c1_90 * (4.75*(spv5(5)-spv5(1)) + 5.5*(spv5(4)-spv5(2))) * dp_sfc

546 endif

547 inty_za(i,j,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) + inty_za_cor(i,j)

548 enddo ; enddo

549 else

550 ! This order of integrating upward and then downward again is necessary with

551 ! a nonlinear equation of state, so that the surface geopotentials will go

552 ! linearly between the values at thickness points, but the bottom geopotentials

553 ! will not now be linear at the sub-grid-scale. Doing this ensures no motion

554 ! with flat isopycnals, even with a nonlinear equation of state.

555 !$OMP parallel do default(shared)

556 do j=js,je ; do i=isq,ieq

557 intx_za(i,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1))

558 enddo ; enddo

559 !$OMP parallel do default(shared)

560 do j=jsq,jeq ; do i=is,ie

561 inty_za(i,j,1) = 0.5*(za(i,j,1) + za(i,j+1,1))

562 enddo ; enddo

563 endif

564

565 do k=1,nz

566 !$OMP parallel do default(shared)

567 do j=js,je ; do i=isq,ieq

568 intx_za(i,j,k+1) = intx_za(i,j,k) - intx_dza(i,j,k)

569 enddo ; enddo

570 enddo

571 do k=1,nz

572 !$OMP parallel do default(shared)

573 do j=jsq,jeq ; do i=is,ie

574 inty_za(i,j,k+1) = inty_za(i,j,k) - inty_dza(i,j,k)

575 enddo ; enddo

576 enddo

577

578 if (cs%debug) then

579 call uvchksum("Prelim int[xy]_za", intx_za, inty_za, g%HI, haloshift=0, &

580 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%L_T_to_m_s**2)

581 call uvchksum("Prelim int[xy]_dza", intx_dza, inty_dza, g%HI, haloshift=0, &

582 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%L_T_to_m_s**2)

583 endif

584

585 if (cs%reset_intxpa_integral) then

586 ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is

587 ! reset intx_za there, then adjust intx_za throughout the water column.

588

589 ! Zero out the 2-d arrays that will be set from various reference interfaces.

590 t_int_w(:,:) = 0.0 ; s_int_w(:,:) = 0.0 ; p_int_w(:,:) = 0.0

591 t_int_e(:,:) = 0.0 ; s_int_e(:,:) = 0.0 ; p_int_e(:,:) = 0.0

592 intx_za_nonlin(:,:) = 0.0 ; intx_za_cor_ri(:,:) = 0.0 ; dp_int_x(:,:) = 0.0

593 do j=js,je ; do i=isq,ieq

594 seek_x_cor(i,j) = (g%mask2dCu(i,j) > 0.)

595 delta_p_x(i,j) = 0.0

596 enddo ; enddo

597

598 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

599 if ((p(i+1,j,2) >= p(i,j,1)) .and. (p(i,j,2) >= p(i+1,j,1))) then

600 ! This is the typical case in the open ocean, so use the topmost interface.

601 t_int_w(i,j) = t_top(i,j) ; t_int_e(i,j) = t_top(i+1,j)

602 s_int_w(i,j) = s_top(i,j) ; s_int_e(i,j) = s_top(i+1,j)

603 p_int_w(i,j) = p(i,j,1) ; p_int_e(i,j) = p(i+1,j,1)

604 intx_za_nonlin(i,j) = intx_za(i,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))

605 dp_int_x(i,j) = p(i+1,j,1)-p(i,j,1)

606 seek_x_cor(i,j) = .false.

607 endif

608 endif ; enddo ; enddo

609

610 do k=1,nz

611 do_more_k = .false.

612 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

613 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not

614 ! activated in the subgrid interpolation.

615 if (((h(i,j,k) > cs%h_nonvanished) .and. (h(i+1,j,k) > cs%h_nonvanished)) .and. &

616 (max(0., p(i,j,1)-p(i+1,j,k+1), p(i+1,j,1)-p(i,j,k+1)) <= 0.0)) then

617 ! Store properties at the bottom of this cell to get a "good estimate" for intxpa at

618 ! the interface below this cell (it might have quadratic pressure dependence if sloped)

619 t_int_w(i,j) = t_b(i,j,k) ; t_int_e(i,j) = t_b(i+1,j,k)

620 s_int_w(i,j) = s_b(i,j,k) ; s_int_e(i,j) = s_b(i+1,j,k)

621 p_int_w(i,j) = p(i,j,k+1) ; p_int_e(i,j) = p(i+1,j,k+1)

622

623 intx_za_nonlin(i,j) = intx_za(i,j,k+1) - 0.5*(za(i,j,k+1) + za(i+1,j,k+1))

624 dp_int_x(i,j) = p(i+1,j,k+1)-p(i,j,k+1)

625 seek_x_cor(i,j) = .false.

626 else

627 do_more_k = .true.

628 endif

629 endif ; enddo ; enddo

630 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.

631 enddo

632

633 if (do_more_k) then

634 if (cs%reset_intxpa_flattest) then

635 ! There are still points where a correction is needed, so use flattest interface

636 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

637 ! choose top layer first

638 t_int_w(i,j) = t_top(i,j) ; t_int_e(i,j) = t_top(i+1,j)

639 s_int_w(i,j) = s_top(i,j) ; s_int_e(i,j) = s_top(i+1,j)

640 p_int_w(i,j) = p(i,j,1) ; p_int_e(i,j) = p(i+1,j,1)

641 intx_za_nonlin(i,j) = intx_za(i,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))

642 dp_int_x(i,j) = p(i+1,j,1)-p(i,j,1)

643 delta_p_x(i,j) = abs(p(i+1,j,1)-p(i,j,1))

644 do k=1,nz

645 if (abs(p(i+1,j,k+1)-p(i,j,k+1)) < delta_p_x(i,j)) then

646 ! bottom of layer is less sloped than top. Use this layer

647 delta_p_x(i,j) = abs(p(i+1,j,k+1)-p(i,j,k+1))

648 t_int_w(i,j) = t_b(i,j,k) ; t_int_e(i,j) = t_b(i+1,j,k)

649 s_int_w(i,j) = s_b(i,j,k) ; s_int_e(i,j) = s_b(i+1,j,k)

650 p_int_w(i,j) = p(i,j,k+1) ; p_int_e(i,j) = p(i+1,j,k+1)

651 intx_za_nonlin(i,j) = intx_za(i,j,k+1) - 0.5*(za(i,j,k+1) + za(i+1,j,k+1))

652 dp_int_x(i,j) = p(i+1,j,k+1)-p(i,j,k+1)

653 endif

654 enddo

655 seek_x_cor(i,j) = .false.

656 endif ; enddo ; enddo

657 else

658 ! There are still points where a correction is needed, so use the top interface.

659 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

660 t_int_w(i,j) = t_top(i,j) ; t_int_e(i,j) = t_top(i+1,j)

661 s_int_w(i,j) = s_top(i,j) ; s_int_e(i,j) = s_top(i+1,j)

662 p_int_w(i,j) = p(i,j,1) ; p_int_e(i,j) = p(i+1,j,1)

663 intx_za_nonlin(i,j) = intx_za(i,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))

664 dp_int_x(i,j) = p(i+1,j,1)-p(i,j,1)

665 seek_x_cor(i,j) = .false.

666 endif ; enddo ; enddo

667 endif

668 endif

669

670 do j=js,je

671 do i=isq,ieq

672 ! This expression assumes that temperature and salinity vary linearly with pressure

673 ! between the corners of the reference interfaces found above to get a correction to

674 ! intx_pa that takes nonlinearities in the equation of state into account.

675 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale

676 ! linear correction so that the pressures and heights match at the end-points. It turns

677 ! out that this linear correction cancels out the mid-point specific volume.

678 ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land.

679 t5(1) = t_int_w(i,j) ; s5(1) = s_int_w(i,j) ; p5(1) = p_int_w(i,j)

680 t5(5) = t_int_e(i,j) ; s5(5) = s_int_e(i,j) ; p5(5) = p_int_e(i,j)

681 t5(2) = 0.25*(3.0*t5(1) + t5(5)) ; t5(4) = 0.25*(3.0*t5(5) + t5(1)) ; t5(3) = 0.5*(t5(5) + t5(1))

682 s5(2) = 0.25*(3.0*s5(1) + s5(5)) ; s5(4) = 0.25*(3.0*s5(5) + s5(1)) ; s5(3) = 0.5*(s5(5) + s5(1))

683 p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1))

684 call calculate_spec_vol(t5, s5, p5, spv5, tv%eqn_of_state, spv_ref=alpha_ref)

685

686 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12

687 intx_za_cor_ri(i,j) = c1_90 * (4.75*(spv5(5)-spv5(1)) + 5.5*(spv5(4)-spv5(2))) * &

688 dp_int_x(i,j) - intx_za_nonlin(i,j)

689 enddo

690 enddo

691

692 ! Repeat the calculations above for v-velocity points.

693 t_int_s(:,:) = 0.0 ; s_int_s(:,:) = 0.0 ; p_int_s(:,:) = 0.0

694 t_int_n(:,:) = 0.0 ; s_int_n(:,:) = 0.0 ; p_int_n(:,:) = 0.0

695 inty_za_nonlin(:,:) = 0.0 ; inty_za_cor_ri(:,:) = 0.0 ; dp_int_y(:,:) = 0.0

696 do j=jsq,jeq ; do i=is,ie

697 seek_y_cor(i,j) = (g%mask2dCv(i,j) > 0.)

698 delta_p_y(i,j) = 0.0

699 enddo ; enddo

700

701 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

702 if ((p(i,j+1,2) >= p(i,j,1)) .and. (p(i,j,2) >= p(i,j+1,1))) then

703 ! This is the typical case in the open ocean, so use the topmost interface.

704 t_int_s(i,j) = t_top(i,j) ; t_int_n(i,j) = t_top(i,j+1)

705 s_int_s(i,j) = s_top(i,j) ; s_int_n(i,j) = s_top(i,j+1)

706 p_int_s(i,j) = p(i,j,1) ; p_int_n(i,j) = p(i,j+1,1)

707 inty_za_nonlin(i,j) = inty_za(i,j,1) - 0.5*(za(i,j,1) + za(i,j+1,1))

708 dp_int_y(i,j) = p(i,j+1,1) - p(i,j,1)

709 seek_y_cor(i,j) = .false.

710 endif

711 endif ; enddo ; enddo

712

713 do k=1,nz

714 do_more_k = .false.

715 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

716 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not

717 ! activated in the subgrid interpolation.

718 if (((h(i,j,k) > cs%h_nonvanished) .and. (h(i,j+1,k) > cs%h_nonvanished)) .and. &

719 (max(0., p(i,j,1)-p(i,j+1,k+1), p(i,j+1,1)-p(i,j,k+1)) <= 0.0)) then

720 ! Store properties at the bottom of this cell to get a "good estimate" for intypa at

721 ! the interface below this cell (it might have quadratic pressure dependence if sloped)

722 t_int_s(i,j) = t_b(i,j,k) ; t_int_n(i,j) = t_b(i,j+1,k)

723 s_int_s(i,j) = s_b(i,j,k) ; s_int_n(i,j) = s_b(i,j+1,k)

724 p_int_s(i,j) = p(i,j,k+1) ; p_int_n(i,j) = p(i,j+1,k+1)

725 inty_za_nonlin(i,j) = inty_za(i,j,k+1) - 0.5*(za(i,j,k+1) + za(i,j+1,k+1))

726 dp_int_y(i,j) = p(i,j+1,k+1) - p(i,j,k+1)

727 seek_y_cor(i,j) = .false.

728 else

729 do_more_k = .true.

730 endif

731 endif ; enddo ; enddo

732 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.

733 enddo

734

735 if (do_more_k) then

736 if (cs%reset_intxpa_flattest) then

737 ! There are still points where a correction is needed, so use flattest interface.

738 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

739 ! choose top interface first

740 t_int_s(i,j) = t_top(i,j) ; t_int_n(i,j) = t_top(i,j+1)

741 s_int_s(i,j) = s_top(i,j) ; s_int_n(i,j) = s_top(i,j+1)

742 p_int_s(i,j) = p(i,j,1) ; p_int_n(i,j) = p(i,j+1,1)

743 inty_za_nonlin(i,j) = inty_za(i,j,1) - 0.5*(za(i,j,1) + za(i,j+1,1))

744 dp_int_y(i,j) = p(i,j+1,1) - p(i,j,1)

745 delta_p_y(i,j) = abs(p(i,j+1,1)-p(i,j,1))

746 do k=1,nz

747 if (abs(p(i,j+1,k+1)-p(i,j,k+1)) < delta_p_y(i,j)) then

748 ! bottom of layer is less sloped than top. Use this layer

749 delta_p_y(i,j) = abs(p(i,j+1,k+1)-p(i,j,k+1))

750 t_int_s(i,j) = t_b(i,j,k) ; t_int_n(i,j) = t_b(i,j+1,k)

751 s_int_s(i,j) = s_b(i,j,k) ; s_int_n(i,j) = s_b(i,j+1,k)

752 p_int_s(i,j) = p(i,j,k+1) ; p_int_n(i,j) = p(i,j+1,k+1)

753 inty_za_nonlin(i,j) = inty_za(i,j,k+1) - 0.5*(za(i,j,k+1) + za(i,j+1,k+1))

754 dp_int_y(i,j) = p(i,j+1,k+1) - p(i,j,k+1)

755 endif

756 enddo

757 seek_y_cor(i,j) = .false.

758 endif ; enddo ; enddo

759 else

760 ! There are still points where a correction is needed, so use the top interface.

761 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

762 t_int_s(i,j) = t_top(i,j) ; t_int_n(i,j) = t_top(i,j+1)

763 s_int_s(i,j) = s_top(i,j) ; s_int_n(i,j) = s_top(i,j+1)

764 p_int_s(i,j) = p(i,j,1) ; p_int_n(i,j) = p(i,j+1,1)

765 inty_za_nonlin(i,j) = inty_za(i,j,1) - 0.5*(za(i,j,1) + za(i,j+1,1))

766 dp_int_y(i,j) = p(i,j+1,1) - p(i,j,1)

767 seek_y_cor(i,j) = .false.

768 endif ; enddo ; enddo

769 endif

770 endif

771

772 do j=jsq,jeq

773 do i=is,ie

774 ! This expression assumes that temperature and salinity vary linearly with pressure

775 ! between the corners of the reference interfaces found above to get a correction to

776 ! intx_pa that takes nonlinearities in the equation of state into account.

777 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale

778 ! linear correction so that the pressures and heights match at the end-points. It turns

779 ! out that this linear correction cancels out the mid-point specific volume.

780 ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land.

781 t5(1) = t_int_s(i,j) ; s5(1) = s_int_s(i,j) ; p5(1) = p_int_s(i,j)

782 t5(5) = t_int_n(i,j) ; s5(5) = s_int_n(i,j) ; p5(5) = p_int_n(i,j)

783 t5(2) = 0.25*(3.0*t5(1) + t5(5)) ; t5(4) = 0.25*(3.0*t5(5) + t5(1)) ; t5(3) = 0.5*(t5(5) + t5(1))

784 s5(2) = 0.25*(3.0*s5(1) + s5(5)) ; s5(4) = 0.25*(3.0*s5(5) + s5(1)) ; s5(3) = 0.5*(s5(5) + s5(1))

785 p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1))

786 call calculate_spec_vol(t5, s5, p5, spv5, tv%eqn_of_state, spv_ref=alpha_ref)

787

788 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12

789 inty_za_cor_ri(i,j) = c1_90 * (4.75*(spv5(5)-spv5(1)) + 5.5*(spv5(4)-spv5(2))) * &

790 dp_int_y(i,j) - inty_za_nonlin(i,j)

791 enddo

792 enddo

793

794 if (cs%debug) then

795 call uvchksum("Pre-reset int[xy]_za", intx_za, inty_za, g%HI, haloshift=0, &

796 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%L_T_to_m_s**2)

797 call uvchksum("int[xy]_za_cor", intx_za_cor_ri, inty_za_cor_ri, g%HI, haloshift=0, &

798 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%L_T_to_m_s**2)

799 call uvchksum("int[xy]_za_nonlin", intx_za_nonlin, inty_za_nonlin, g%HI, haloshift=0, &

800 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%L_T_to_m_s**2)

801 call uvchksum("dp_int_[xy]", dp_int_x, dp_int_y, g%HI, haloshift=0, &

802 symmetric=g%Domain%symmetric, unscale=us%RL2_T2_to_Pa)

803 endif

804

805 ! Correct intx_pa and inty_pa at each interface using vertically constant corrections.

806 do k=1,nz+1 ; do j=js,je ; do i=isq,ieq

807 intx_za(i,j,k) = intx_za(i,j,k) + intx_za_cor_ri(i,j)

808 enddo ; enddo ; enddo

809

810 do k=1,nz+1 ; do j=jsq,jeq ; do i=is,ie

811 inty_za(i,j,k) = inty_za(i,j,k) + inty_za_cor_ri(i,j)

812 enddo ; enddo ; enddo

813

814 if (cs%debug) then

815 call uvchksum("Post-reset int[xy]_za", intx_za, inty_za, g%HI, haloshift=0, &

816 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%L_T_to_m_s**2)

817 endif

818

819 endif ! intx_za and inty_za have now been reset to reflect the properties of an unimpeded interface.

820

821 !$OMP parallel do default(shared) private(dp)

822 do k=1,nz

823 do j=jsq,jeq+1 ; do i=isq,ieq+1

824 dp(i,j) = h_to_rl2_t2 * h(i,j,k)

825 enddo ; enddo

826

827 ! Find the horizontal pressure gradient accelerations.

828 ! These expressions for the accelerations have been carefully checked in

829 ! a set of idealized cases, and should be bug-free.

830 do j=js,je ; do i=isq,ieq

831 pfu(i,j,k) = ( ((za(i,j,k+1)*dp(i,j) + intp_dza(i,j,k)) - &

832 (za(i+1,j,k+1)*dp(i+1,j) + intp_dza(i+1,j,k))) + &

833 ((dp(i+1,j) - dp(i,j)) * intx_za(i,j,k+1) - &

834 (p(i+1,j,k) - p(i,j,k)) * intx_dza(i,j,k)) ) * &

835 (2.0*g%IdxCu(i,j) / ((dp(i,j) + dp(i+1,j)) + dp_neglect))

836 enddo ; enddo

837

838 do j=jsq,jeq ; do i=is,ie

839 pfv(i,j,k) = (((za(i,j,k+1)*dp(i,j) + intp_dza(i,j,k)) - &

840 (za(i,j+1,k+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + &

841 ((dp(i,j+1) - dp(i,j)) * inty_za(i,j,k+1) - &

842 (p(i,j+1,k) - p(i,j,k)) * inty_dza(i,j,k))) * &

843 (2.0*g%IdyCv(i,j) / ((dp(i,j) + dp(i,j+1)) + dp_neglect))

844 enddo ; enddo

845 enddo

846

847 if (cs%GFS_scale < 1.0) then

848 ! Adjust the Montgomery potential to make this a reduced gravity model.

849 if (use_eos) then

850 !$OMP parallel do default(shared) private(rho_in_situ)

851 do j=jsq,jeq+1

852 call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, &

853 tv%eqn_of_state, eosdom)

854

855 do i=isq,ieq+1

856 dm(i,j) = (cs%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j,1))

857 enddo

858 enddo

859 else

860 !$OMP parallel do default(shared)

861 do j=jsq,jeq+1 ; do i=isq,ieq+1

862 dm(i,j) = (cs%GFS_scale - 1.0) * (p(i,j,1)*(1.0/gv%Rlay(1) - alpha_ref) + za(i,j,1))

863 enddo ; enddo

864 endif

865

866 !$OMP parallel do default(shared)

867 do k=1,nz

868 do j=js,je ; do i=isq,ieq

869 pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)

870 enddo ; enddo

871 do j=jsq,jeq ; do i=is,ie

872 pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)

873 enddo ; enddo

874 enddo

875 endif

876

877 if (present(pbce)) then

878 call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce)

879 endif

880

881 if (present(eta)) then

882 pa_to_h = 1.0 / (gv%g_Earth * gv%H_to_RZ)

883 if (use_p_atm) then

884 !$OMP parallel do default(shared)

885 do j=jsq,jeq+1 ; do i=isq,ieq+1

886 eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h ! eta has the same units as h.

887 enddo ; enddo

888 else

889 !$OMP parallel do default(shared)

890 do j=jsq,jeq+1 ; do i=isq,ieq+1

891 eta(i,j) = p(i,j,nz+1)*pa_to_h ! eta has the same units as h.

892 enddo ; enddo

893 endif

894 endif

895

896 if (cs%id_MassWt_u>0) call post_data(cs%id_MassWt_u, masswt_u, cs%diag)

897 if (cs%id_MassWt_v>0) call post_data(cs%id_MassWt_v, masswt_v, cs%diag)

898

899 ! Diagnostics for tidal forcing and SAL height anomaly

900 if (cs%id_e_tide>0) then

901 ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL.

902 ! New diagnostics are given for each individual field.

903 if (cs%tides_answer_date>20230630) then ; do j=jsq,jeq+1 ; do i=isq,ieq+1

904 e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j)

905 enddo ; enddo ; endif

906 call post_data(cs%id_e_tide, e_sal_and_tide, cs%diag)

907 endif

908 if (cs%id_e_sal>0) call post_data(cs%id_e_sal, e_sal, cs%diag)

909 if (cs%id_e_tidal_eq>0) call post_data(cs%id_e_tidal_eq, e_tidal_eq, cs%diag)

910 if (cs%id_e_tidal_sal>0) call post_data(cs%id_e_tidal_sal, e_tidal_sal, cs%diag)

911

912 ! Diagnostics for tidal forcing and SAL horizontal gradients

913 if (cs%calculate_SAL .and. (associated(adp%sal_u) .or. associated(adp%sal_v))) then

914 if (cs%tides) then ; do j=jsq,jeq+1 ; do i=isq,ieq+1

915 e_sal(i,j) = e_sal(i,j) + e_tidal_sal(i,j)

916 enddo ; enddo ; endif

917 if (associated(adp%sal_u)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq

918 adp%sal_u(i,j,k) = (e_sal(i+1,j) - e_sal(i,j)) * gv%g_Earth * g%IdxCu(i,j)

919 enddo ; enddo ; enddo ; endif

920 if (associated(adp%sal_v)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie

921 adp%sal_v(i,j,k) = (e_sal(i,j+1) - e_sal(i,j)) * gv%g_Earth * g%IdyCv(i,j)

922 enddo ; enddo ; enddo ; endif

923 if (cs%id_sal_u>0) call post_data(cs%id_sal_u, adp%sal_u, cs%diag)

924 if (cs%id_sal_v>0) call post_data(cs%id_sal_v, adp%sal_v, cs%diag)

925 endif

926

927 if (cs%tides .and. (associated(adp%tides_u) .or. associated(adp%tides_v))) then

928 if (associated(adp%tides_u)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq

929 adp%tides_u(i,j,k) = (e_tidal_eq(i+1,j) - e_tidal_eq(i,j)) * gv%g_Earth * g%IdxCu(i,j)

930 enddo ; enddo ; enddo ; endif

931 if (associated(adp%tides_v)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie

932 adp%tides_v(i,j,k) = (e_tidal_eq(i,j+1) - e_tidal_eq(i,j)) * gv%g_Earth * g%IdyCv(i,j)

933 enddo ; enddo ; enddo ; endif

934 if (cs%id_tides_u>0) call post_data(cs%id_tides_u, adp%tides_u, cs%diag)

935 if (cs%id_tides_v>0) call post_data(cs%id_tides_v, adp%tides_v, cs%diag)

936 endif

937end subroutine pressureforce_fv_nonbouss

938

939!> \brief Boussinesq analytically-integrated finite volume form of pressure gradient

940!!

941!! Determines the acceleration due to hydrostatic pressure forces, using

942!! the finite volume form of the terms and analytic integrals in depth.

943!!

944!! To work, the following fields must be set outside of the usual (is:ie,js:je)

945!! range before this subroutine is called:

946!! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).

947subroutine pressureforce_fv_bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp, p_atm, pbce, eta)

948 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

949 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure

950 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

951 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m]

952 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables

953 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: pfu !< Zonal acceleration [L T-2 ~> m s-2]

954 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: pfv !< Meridional acceleration [L T-2 ~> m s-2]

955 type(pressureforce_fv_cs), intent(in) :: cs !< Finite volume PGF control structure

956 type(ale_cs), pointer :: ale_csp !< ALE control structure

957 type(accel_diag_ptrs), pointer :: adp !< Acceleration diagnostic pointers

958 real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean

959 !! or atmosphere-ocean interface [R L2 T-2 ~> Pa].

960 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure

961 !! anomaly in each layer due to eta anomalies

962 !! [L2 T-2 H-1 ~> m s-2].

963 real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The sea-surface height used to

964 !! calculate PFu and PFv [H ~> m], with any

965 !! tidal contributions.

966 ! Local variables

967 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height in depth units [Z ~> m].

968 real, dimension(SZI_(G),SZJ_(G)) :: &

969 e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m].

970 e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].

971 e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources

972 ! [Z ~> m].

973 e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading

974 ! specific to tides [Z ~> m].

975 z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m]

976 ssh, & ! Sea surfae height anomaly for self-attraction and loading. Used if

977 ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m].

978 pbot, & ! Total bottom pressure for self-attraction and loading. Used if

979 ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa].

980 dm ! The barotropic adjustment to the Montgomery potential to

981 ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].

982 real, dimension(SZI_(G)) :: &

983 rho_cv_bl ! The coordinate potential density in the deepest variable

984 ! density near-surface layer [R ~> kg m-3].

985 real, dimension(SZI_(G),SZJ_(G)) :: &

986 dz_geo ! The change in geopotential thickness through a layer [L2 T-2 ~> m2 s-2].

987 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &

988 pa ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the

989 ! the interface atop a layer [R L2 T-2 ~> Pa].

990 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

991 dpa, & ! The change in pressure anomaly between the top and bottom

992 ! of a layer [R L2 T-2 ~> Pa].

993 intz_dpa ! The vertical integral in depth of the pressure anomaly less the

994 ! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa].

995 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &

996 intx_pa ! The zonal integral of the pressure anomaly along the interface

997 ! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].

998 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &

999 intx_dpa ! The change in intx_pa through a layer [R L2 T-2 ~> Pa].

1000 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &

1001 inty_pa ! The meridional integral of the pressure anomaly along the

1002 ! interface atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].

1003 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &

1004 inty_dpa ! The change in inty_pa through a layer [R L2 T-2 ~> Pa].

1005 real, dimension(SZIB_(G),SZJ_(G)) :: &

1006 intx_pa_cor ! Correction for curvature in intx_pa [R L2 T-2 ~> Pa]

1007 real, dimension(SZI_(G),SZJB_(G)) :: &

1008 inty_pa_cor ! Correction for curvature in inty_pa [R L2 T-2 ~> Pa]

1009

1010 ! These variables are used with reset_intxpa_integral. The values are taken from different

1011 ! interfaces as a function of position.

1012 real, dimension(SZIB_(G),SZJ_(G)) :: &

1013 t_int_w, t_int_e, & ! Temperatures on the reference interface to the east and west of a u-point [C ~> degC]

1014 s_int_w, s_int_e, & ! Salinities on the reference interface to the east and west of a u-point [S ~> ppt]

1015 p_int_w, p_int_e, & ! Pressures on the reference interface to the east and west of a u-point [R L2 T-2 ~> Pa]

1016 intx_pa_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface

1017 ! from the value that would be obtained from assuming that pressure varies

1018 ! linearly with depth along that interface [R L2 T-2 ~> Pa].

1019 dgeo_x, & ! The change in x in geopotenial height along the reference interface [L2 T-2 ~> m2 s-2]

1020 intx_pa_cor_ri ! The correction to intx_pa based on the reference interface calculations [R L2 T-2 ~> Pa]

1021 real, dimension(SZI_(G),SZJB_(G)) :: &

1022 t_int_s, t_int_n, & ! Temperatures on the reference interface to the north and south of a v-point [C ~> degC]

1023 s_int_s, s_int_n, & ! Salinities on the reference interface to the north and south of a v-point [S ~> ppt]

1024 p_int_s, p_int_n, & ! Pressures on the reference interface to the north and south of a v-point [R L2 T-2 ~> Pa]

1025 inty_pa_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface

1026 ! from the value that would be obtained from assuming that pressure varies

1027 ! linearly with depth along that interface [R L2 T-2 ~> Pa].

1028 dgeo_y, & ! The change in y in geopotenial height along the reference interface [L2 T-2 ~> m2 s-2]

1029 inty_pa_cor_ri ! The correction to inty_pa based on the reference interface calculations [R L2 T-2 ~> Pa]

1030 logical, dimension(SZIB_(G),SZJ_(G)) :: &

1031 seek_x_cor ! If true, try to find a u-point interface that would provide a better estimate

1032 ! of the curvature terms in the intx_pa.

1033 logical, dimension(SZI_(G),SZJB_(G)) :: &

1034 seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate

1035 ! of the curvature terms in the inty_pa.

1036 real, dimension(SZIB_(G),SZJ_(G)) :: &

1037 delta_z_x ! If using flattest interface for reset integral, store x interface differences [Z ~> m]

1038 real, dimension(SZI_(G),SZJB_(G)) :: &

1039 delta_z_y ! If using flattest interface for reset integral, store y interface differences [Z ~> m]

1040

1041 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &

1042 t_tmp, & ! Temporary array of temperatures where layers that are lighter

1043 ! than the mixed layer have the mixed layer's properties [C ~> degC].

1044 s_tmp ! Temporary array of salinities where layers that are lighter

1045 ! than the mixed layer have the mixed layer's properties [S ~> ppt].

1046 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

1047 s_t, s_b, & ! Top and bottom edge values for linear reconstructions

1048 ! of salinity within each layer [S ~> ppt].

1049 t_t, t_b ! Top and bottom edge values for linear reconstructions

1050 ! of temperature within each layer [C ~> degC].

1051 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &

1052 masswt_u ! The fractional mass weighting at a u-point [nondim].

1053 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &

1054 masswt_v ! The fractional mass weighting at a v-point [nondim].

1055 real, dimension(SZI_(G),SZJ_(G)) :: &

1056 t_top, & ! Temperature of top layer used with correction_intxpa [C ~> degC]

1057 s_top, & ! Salinity of top layer used with correction_intxpa [S ~> ppt]

1058 rho_top ! Density anomaly of top layer used in calculating intx_pa_cor and inty_pa_cor

1059 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &

1060 rho_pgf, rho_stanley_pgf ! Density [R ~> kg m-3] from EOS with and without SGS T variance

1061 ! in Stanley parameterization.

1062 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &

1063 p_stanley ! Pressure [R L2 T-2 ~> Pa] estimated with Rho_0

1064 real :: zeros(szi_(g)) ! An array of zero values that can be used as an argument [various]

1065 real :: rho_in_situ(szi_(g)) ! The in situ density [R ~> kg m-3].

1066 real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate

1067 ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).

1068 real :: p_surf_eos(szi_(g)) ! The pressure at the ocean surface determined from the surface height,

1069 ! consistent with what is used in the density integral routines [R L2 T-2 ~> Pa]

1070 real :: p0(szi_(g)) ! An array of zeros to use for pressure [R L2 T-2 ~> Pa].

1071 real :: dz_geo_sfc ! The change in surface geopotential height between adjacent cells [L2 T-2 ~> m2 s-2]

1072 real :: gxrho0 ! The gravitational acceleration times mean ocean density [R L2 Z-1 T-2 ~> Pa m-1]

1073 real :: gxrho_ref ! The gravitational acceleration times reference density [R L2 Z-1 T-2 ~> Pa m-1]

1074 real :: rho0_int_density ! Rho0 used in int_density_dz_* subroutines [R ~> kg m-3]

1075 real :: rho0_set_pbce ! Rho0 used in set_pbce_Bouss subroutine [R ~> kg m-3]

1076 real :: h_neglect ! A thickness that is so small it is usually lost

1077 ! in roundoff and can be neglected [H ~> m].

1078 real :: i_rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].

1079 real :: g_rho0 ! G_Earth / Rho_0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].

1080 real :: i_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1]

1081 real :: rho_ref ! The reference density [R ~> kg m-3].

1082 real :: dz_neglect ! A minimal thickness [Z ~> m], like e.

1083 real :: dz_nonvanished ! A small thickness considered to be vanished for mass weighting [Z ~> m]

1084 real :: h_to_rl2_t2 ! A factor to convert from thickness units (H) to pressure

1085 ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].

1086 real :: t5(5) ! Temperatures and salinities at five quadrature points [C ~> degC]

1087 real :: s5(5) ! Salinities at five quadrature points [S ~> ppt]

1088 real :: p5(5) ! Full pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa]

1089 real :: pa5(5) ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at five quadrature points [R L2 T-2 ~> Pa].

1090 real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3]

1091 real :: wt_r ! A weighting factor [nondim]

1092 real, parameter :: c1_6 = 1.0/6.0 ! A rational constant [nondim]

1093 real, parameter :: c1_90 = 1.0/90.0 ! A rational constant [nondim]

1094 logical :: use_p_atm ! If true, use the atmospheric pressure.

1095 logical :: use_ale ! If true, use an ALE pressure reconstruction.

1096 logical :: use_eos ! If true, density is calculated from T & S using an equation of state.

1097 logical :: do_more_k ! If true, there are still points where a flatter interface remains to be found.

1098 type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.

1099 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state

1100 integer, dimension(2) :: eosdom_h ! The i-computational domain for the equation of state at tracer points

1101 integer, dimension(2) :: eosdom_u ! The i-computational domain for the equation of state at u-velocity points

1102 integer, dimension(2) :: eosdom_v ! The i-computational domain for the equation of state at v-velocity points

1103 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb

1104 integer :: i, j, k, m

1105

1106 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

1107 nkmb=gv%nk_rho_varies

1108 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

1109 eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)

1110 eosdom_u(1) = isq - (g%IsdB-1) ; eosdom_u(2) = ieq - (g%IsdB-1)

1111 eosdom_v(1) = is - (g%isd-1) ; eosdom_v(2) = ie - (g%isd-1)

1112

1113 if (.not.cs%initialized) call mom_error(fatal, &

1114 "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.")

1115

1116 use_p_atm = associated(p_atm)

1117 use_eos = associated(tv%eqn_of_state)

1118 do i=isq,ieq+1 ; p0(i) = 0.0 ; enddo

1119 use_ale = .false.

1120 if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos

1121

1122 h_neglect = gv%H_subroundoff

1123 dz_neglect = gv%dZ_subroundoff

1124 dz_nonvanished = gv%H_to_Z*cs%h_nonvanished

1125 i_rho0 = 1.0 / gv%Rho0

1126 g_rho0 = gv%g_Earth / gv%Rho0

1127 gxrho0 = gv%g_Earth * gv%Rho0

1128 rho_ref = cs%rho_ref

1129

1130 if (cs%rho_ref_bug) then

1131 rho0_int_density = rho_ref

1132 rho0_set_pbce = rho_ref

1133 gxrho_ref = gxrho0

1134 i_g_rho = 1.0 / (rho_ref * gv%g_Earth)

1135 else

1136 rho0_int_density = gv%Rho0

1137 rho0_set_pbce = gv%Rho0

1138 gxrho_ref = gv%g_Earth * rho_ref

1139 i_g_rho = 1.0 / (gv%rho0 * gv%g_Earth)

1140 endif

1141

1142 if ((cs%id_MassWt_u > 0) .or. (cs%id_MassWt_v > 0)) then

1143 masswt_u(:,:,:) = 0.0 ; masswt_v(:,:,:) = 0.0

1144 endif

1145

1146 do j=jsq,jeq+1 ; do i=isq,ieq+1

1147 e(i,j,nz+1) = -g%bathyT(i,j)

1148 enddo ; enddo

1149

1150 ! The following two if-blocks are used to recover old answers for self-attraction and loading

1151 ! (SAL) and tides only. The old algorithm moves interface heights before density calculations,

1152 ! and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024).

1153 ! See the code right after Pa calculation loop for the new algorithm.

1154

1155 ! Calculate and add SAL geopotential anomaly to interface height (old answers)

1156 if (cs%calculate_SAL .and. cs%tides_answer_date<=20250131) then

1157 !$OMP parallel do default(shared)

1158 do j=jsq,jeq+1

1159 do i=isq,ieq+1

1160 ssh(i,j) = min(-g%bathyT(i,j) - g%meanSL(i,j), 0.0)

1161 enddo

1162 do k=1,nz ; do i=isq,ieq+1

1163 ssh(i,j) = ssh(i,j) + h(i,j,k)*gv%H_to_Z

1164 enddo ; enddo

1165 enddo

1166 call calc_sal(ssh, e_sal, g, cs%SAL_CSp, tmp_scale=us%Z_to_m)

1167

1168 if (cs%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131]

1169 !$OMP parallel do default(shared)

1170 do j=jsq,jeq+1 ; do i=isq,ieq+1

1171 e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)

1172 enddo ; enddo

1173 endif

1174 endif

1175

1176 ! Calculate and add tidal geopotential anomaly to interface height (old answers)

1177 if (cs%tides .and. cs%tides_answer_date<=20250131) then

1178 if (cs%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131]

1179 call calc_tidal_forcing(cs%Time, e_tidal_eq, e_tidal_sal, g, us, cs%tides_CSp)

1180 !$OMP parallel do default(shared)

1181 do j=jsq,jeq+1 ; do i=isq,ieq+1

1182 e(i,j,nz+1) = e(i,j,nz+1) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))

1183 enddo ; enddo

1184 else ! answers_date before 20230701

1185 if (.not.cs%calculate_SAL) e_sal(:,:) = 0.0

1186 call calc_tidal_forcing_legacy(cs%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, &

1187 g, us, cs%tides_CSp)

1188 !$OMP parallel do default(shared)

1189 do j=jsq,jeq+1 ; do i=isq,ieq+1

1190 e(i,j,nz+1) = e(i,j,nz+1) - e_sal_and_tide(i,j)

1191 enddo ; enddo

1192 endif

1193 endif

1194

1195 !$OMP parallel do default(shared)

1196 do j=jsq,jeq+1 ; do k=nz,1,-1 ; do i=isq,ieq+1

1197 e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z

1198 enddo ; enddo ; enddo

1199

1200 if (use_eos) then

1201 if (nkmb>0) then

1202 ! With a bulk mixed layer, replace the T & S of any layers that are lighter than the buffer

1203 ! layer with the properties of the buffer layer. These layers will be massless anyway, and

1204 ! it avoids any formal calculations with hydrostatically unstable profiles.

1205 tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp

1206 tv_tmp%eqn_of_state => tv%eqn_of_state

1207

1208 do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo

1209 !$OMP parallel do default(shared) private(Rho_cv_BL)

1210 do j=jsq,jeq+1

1211 do k=1,nkmb ; do i=isq,ieq+1

1212 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)

1213 enddo ; enddo

1214 call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, rho_cv_bl(:), &

1215 tv%eqn_of_state, eosdom)

1216

1217 do k=nkmb+1,nz ; do i=isq,ieq+1

1218 if (gv%Rlay(k) < rho_cv_bl(i)) then

1219 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)

1220 else

1221 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)

1222 endif

1223 enddo ; enddo

1224 enddo

1225 else

1226 tv_tmp%T => tv%T ; tv_tmp%S => tv%S

1227 tv_tmp%eqn_of_state => tv%eqn_of_state

1228 endif

1229 endif

1230

1231 ! If regridding is activated, do a linear reconstruction of salinity

1232 ! and temperature across each layer. The subscripts 't' and 'b' refer

1233 ! to top and bottom values within each layer (these are the only degrees

1234 ! of freedom needed to know the linear profile).

1235 if ( use_ale .and. (cs%Recon_Scheme == 1) ) then

1236 call ts_plm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)

1237 elseif ( use_ale .and. (cs%Recon_Scheme == 2) ) then

1238 call ts_ppm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)

1239 elseif ( use_ale .and. (cs%Recon_Scheme == 3) ) then

1240 call ts_plm_wls_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h)

1241 elseif (cs%reset_intxpa_integral) then

1242 do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1

1243 t_b(i,j,k) = tv%T(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)

1244 enddo ; enddo ; enddo

1245 endif

1246

1247 ! Set the surface boundary conditions on pressure anomaly and its horizontal

1248 ! integrals, assuming that the surface pressure anomaly varies linearly

1249 ! in x and y.

1250 if (use_p_atm) then

1251 !$OMP parallel do default(shared)

1252 do j=jsq,jeq+1 ; do i=isq,ieq+1

1253 pa(i,j,1) = gxrho_ref * (e(i,j,1) - g%Z_ref) + p_atm(i,j)

1254 enddo ; enddo

1255 else

1256 !$OMP parallel do default(shared)

1257 do j=jsq,jeq+1 ; do i=isq,ieq+1

1258 pa(i,j,1) = gxrho_ref * (e(i,j,1) - g%Z_ref)

1259 enddo ; enddo

1260 endif

1261

1262 if (cs%use_SSH_in_Z0p .and. use_p_atm) then

1263 do j=jsq,jeq+1 ; do i=isq,ieq+1

1264 z_0p(i,j) = e(i,j,1) + p_atm(i,j) * i_g_rho

1265 enddo ; enddo

1266 elseif (cs%use_SSH_in_Z0p) then

1267 do j=jsq,jeq+1 ; do i=isq,ieq+1

1268 z_0p(i,j) = e(i,j,1)

1269 enddo ; enddo

1270 else

1271 do j=jsq,jeq+1 ; do i=isq,ieq+1

1272 z_0p(i,j) = g%meanSL(i,j)

1273 enddo ; enddo

1274 endif

1275

1276 do k=1,nz

1277 ! Calculate 4 integrals through the layer that are required in the

1278 ! subsequent calculation.

1279 if (use_eos) then

1280 ! The following routine computes the integrals that are needed to

1281 ! calculate the pressure gradient force. Linear profiles for T and S are

1282 ! assumed when regridding is activated. Otherwise, the previous version

1283 ! is used, whereby densities within each layer are constant no matter

1284 ! where the layers are located.

1285 if ( use_ale .and. cs%Recon_Scheme > 0 ) then

1286 if ( cs%Recon_Scheme == 1 .or. cs%Recon_Scheme == 3 ) then

1287 call int_density_dz_generic_plm(k, tv, t_t, t_b, s_t, s_b, e, &

1288 rho_ref, rho0_int_density, gv%g_Earth, dz_neglect, g%bathyT, &

1289 g%HI, gv, tv%eqn_of_state, us, cs%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &

1290 intx_dpa(:,:,k), inty_dpa(:,:,k), &

1291 masswghtinterp=cs%MassWghtInterp, &

1292 use_inaccurate_form=cs%use_inaccurate_pgf_rho_anom, z_0p=z_0p, &

1293 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, h_nv=dz_nonvanished)

1294 elseif ( cs%Recon_Scheme == 2 ) then

1295 call int_density_dz_generic_ppm(k, tv, t_t, t_b, s_t, s_b, e, &

1296 rho_ref, rho0_int_density, gv%g_Earth, dz_neglect, g%bathyT, &

1297 g%HI, gv, tv%eqn_of_state, us, cs%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &

1298 intx_dpa(:,:,k), inty_dpa(:,:,k), &

1299 masswghtinterp=cs%MassWghtInterp, z_0p=z_0p, &

1300 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, h_nv=dz_nonvanished)

1301 endif

1302 else

1303 call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,k), e(:,:,k+1), &

1304 rho_ref, rho0_int_density, gv%g_Earth, g%HI, tv%eqn_of_state, us, dpa(:,:,k), &

1305 intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), g%bathyT, e(:,:,1), dz_neglect, &

1306 cs%MassWghtInterp, z_0p=z_0p, &

1307 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, h_nv=dz_nonvanished)

1308 endif

1309 if (gv%Z_to_H /= 1.0) then

1310 !$OMP parallel do default(shared)

1311 do j=jsq,jeq+1 ; do i=isq,ieq+1

1312 intz_dpa(i,j,k) = intz_dpa(i,j,k)*gv%Z_to_H

1313 enddo ; enddo

1314 endif

1315 if ((cs%id_MassWt_u > 0) .or. (cs%id_MassWt_v > 0)) &

1316 call diagnose_mass_weight_z(e(:,:,k), e(:,:,k+1), g%bathyT, e(:,:,1), dz_neglect, cs%MassWghtInterp, &

1317 g%HI, masswt_u(:,:,k), masswt_v(:,:,k), &

1318 masswghtinterpvanonly=cs%MassWghtInterpVanOnly, h_nv=cs%h_nonvanished)

1319 else

1320 !$OMP parallel do default(shared)

1321 do j=jsq,jeq+1 ; do i=isq,ieq+1

1322 dz_geo(i,j) = gv%g_Earth * gv%H_to_Z*h(i,j,k)

1323 dpa(i,j,k) = (gv%Rlay(k) - rho_ref) * dz_geo(i,j)

1324 intz_dpa(i,j,k) = 0.5*(gv%Rlay(k) - rho_ref) * dz_geo(i,j)*h(i,j,k)

1325 enddo ; enddo

1326 !$OMP parallel do default(shared)

1327 do j=js,je ; do i=isq,ieq

1328 intx_dpa(i,j,k) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i+1,j))

1329 enddo ; enddo

1330 !$OMP parallel do default(shared)

1331 do j=jsq,jeq ; do i=is,ie

1332 inty_dpa(i,j,k) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i,j+1))

1333 enddo ; enddo

1334 endif

1335 enddo

1336

1337 ! Set the pressure anomalies at the interfaces.

1338 do k=1,nz

1339 !$OMP parallel do default(shared)

1340 do j=jsq,jeq+1 ; do i=isq,ieq+1

1341 pa(i,j,k+1) = pa(i,j,k) + dpa(i,j,k)

1342 enddo ; enddo

1343 enddo

1344

1345 ! Calculate and add SAL geopotential anomaly to interface height (new answers)

1346 if (cs%calculate_SAL .and. cs%tides_answer_date>20250131) then

1347 if (cs%sal_use_bpa) then

1348 !$OMP parallel do default(shared)

1349 do j=jsq,jeq+1 ; do i=isq,ieq+1

1350 pbot(i,j) = pa(i,j,nz+1) - gxrho_ref * (e(i,j,nz+1) - g%Z_ref)

1351 enddo ; enddo

1352 call calc_sal(pbot, e_sal, g, cs%SAL_CSp, tmp_scale=us%Z_to_m)

1353 else

1354 !$OMP parallel do default(shared)

1355 do j=jsq,jeq+1 ; do i=isq,ieq+1

1356 ssh(i,j) = e(i,j,1) - g%Z_ref

1357 ! Remove above sea level topography at floodable cells

1358 ssh(i,j) = ssh(i,j) - max(-g%bathyT(i,j)-g%meanSL(i,j), 0.0)

1359 enddo ; enddo

1360 call calc_sal(ssh, e_sal, g, cs%SAL_CSp, tmp_scale=us%Z_to_m)

1361 endif

1362 if (.not.cs%bq_sal_tides) then ; do k=1,nz+1

1363 !$OMP parallel do default(shared)

1364 do j=jsq,jeq+1 ; do i=isq,ieq+1

1365 e(i,j,k) = e(i,j,k) - e_sal(i,j)

1366 pa(i,j,k) = pa(i,j,k) - gxrho_ref * e_sal(i,j)

1367 enddo ; enddo

1368 enddo ; endif

1369 endif

1370

1371 ! Calculate and add tidal geopotential anomaly to interface height (new answers)

1372 if (cs%tides .and. cs%tides_answer_date>20250131) then

1373 call calc_tidal_forcing(cs%Time, e_tidal_eq, e_tidal_sal, g, us, cs%tides_CSp)

1374 if (.not.cs%bq_sal_tides) then ; do k=1,nz+1

1375 !$OMP parallel do default(shared)

1376 do j=jsq,jeq+1 ; do i=isq,ieq+1

1377 e(i,j,k) = e(i,j,k) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))

1378 pa(i,j,k) = pa(i,j,k) - gxrho_ref * (e_tidal_eq(i,j) + e_tidal_sal(i,j))

1379 enddo ; enddo

1380 enddo ; endif

1381 endif

1382

1383 if (cs%correction_intxpa .or. cs%reset_intxpa_integral) then

1384 ! Determine surface temperature and salinity for use in the pressure gradient corrections

1385 if (use_ale .and. (cs%Recon_Scheme > 0)) then

1386 do j=jsq,jeq+1 ; do i=isq,ieq+1

1387 t_top(i,j) = t_t(i,j,1) ; s_top(i,j) = s_t(i,j,1)

1388 enddo ; enddo

1389 else

1390 do j=jsq,jeq+1 ; do i=isq,ieq+1

1391 t_top(i,j) = tv%T(i,j,1) ; s_top(i,j) = tv%S(i,j,1)

1392 enddo ; enddo

1393 endif

1394 endif

1395

1396 if (cs%correction_intxpa) then

1397 ! Determine surface density for use in the pressure gradient corrections

1398 !$OMP parallel do default(shared) private(p_surf_EOS)

1399 do j=jsq,jeq+1

1400 ! P_surf_EOS here is consistent with the pressure that is used in the int_density_dz routines.

1401 do i=isq,ieq+1 ; p_surf_eos(i) = -gxrho0*(e(i,j,1) - z_0p(i,j)) ; enddo

1402 call calculate_density(t_top(:,j), s_top(:,j), p_surf_eos, rho_top(:,j), &

1403 tv%eqn_of_state, eosdom, rho_ref=rho_ref)

1404 enddo

1405

1406 if (cs%debug) then

1407 call hchksum(rho_top, "intx_pa rho_top", g%HI, haloshift=1, unscale=us%R_to_kg_m3)

1408 call hchksum(e(:,:,1), "intx_pa e(1)", g%HI, haloshift=1, unscale=us%Z_to_m)

1409 call hchksum(pa(:,:,1), "intx_pa pa(1)", g%HI, haloshift=1, unscale=us%RL2_T2_to_Pa)

1410 endif

1411

1412 ! This version attempts to correct for hydrostatic variations in surface pressure under ice.

1413 !$OMP parallel do default(shared) private(dz_geo_sfc)

1414 do j=js,je ; do i=isq,ieq

1415 intx_pa_cor(i,j) = 0.0

1416 dz_geo_sfc = gv%g_Earth * (e(i+1,j,1)-e(i,j,1))

1417 if ((dz_geo_sfc * rho_ref - (pa(i+1,j,1)-pa(i,j,1)))*dz_geo_sfc > 0.0) then

1418 ! The pressure/depth relationship has a positive implied density given by

1419 ! rho_implied = rho_ref - (pa(i+1,j,1)-pa(i,j,1)) / dz_geo_sfc

1420 if (-dz_geo_sfc * (pa(i+1,j,1)-pa(i,j,1)) > &

1421 0.25*((rho_top(i+1,j)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then

1422 ! The pressure difference is at least half the size of the difference expected by hydrostatic

1423 ! balance. This test gets rid of pressure differences that are small, e.g. open ocean.

1424 ! Use 5 point quadrature to calculate intxpa

1425 t5(1) = t_top(i,j) ; t5(5) = t_top(i+1,j)

1426 s5(1) = s_top(i,j) ; s5(5) = s_top(i+1,j)

1427 pa5(1) = pa(i,j,1) ; pa5(5) = pa(i+1,j,1)

1428 ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines.

1429 p5(1) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1430 p5(5) = -gxrho0*(e(i+1,j,1) - z_0p(i,j))

1431 do m=2,4

1432 wt_r = 0.25*real(m-1)

1433 t5(m) = t5(1) + (t5(5)-t5(1))*wt_r

1434 s5(m) = s5(1) + (s5(5)-s5(1))*wt_r

1435 p5(m) = p5(1) + (p5(5)-p5(1))*wt_r

1436 enddo !m

1437 call calculate_density(t5, s5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)

1438

1439 ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure

1440 ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule

1441 ! quadrature to find the integrated correction to the integral of pressure along the interface.

1442 ! The derivation for this expression is shown below in the y-direction version.

1443 intx_pa_cor(i,j) = c1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc

1444 ! Note that (4.75 + 5.5/2) / 90 = 1/12, so this is consistent with the linear result below.

1445 endif

1446 endif

1447 intx_pa(i,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1)) + intx_pa_cor(i,j)

1448 enddo ; enddo

1449 !$OMP parallel do default(shared) private(dz_geo_sfc)

1450 do j=jsq,jeq ; do i=is,ie

1451 inty_pa_cor(i,j) = 0.0

1452 dz_geo_sfc = gv%g_Earth * (e(i,j+1,1)-e(i,j,1))

1453 if ((dz_geo_sfc * rho_ref - (pa(i,j+1,1)-pa(i,j,1)))*dz_geo_sfc > 0.0) then

1454 ! The pressure/depth relationship has a positive implied density

1455 if (-dz_geo_sfc * (pa(i,j+1,1)-pa(i,j,1)) > &

1456 0.25*((rho_top(i,j+1)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then

1457 ! The pressure difference is at least half the size of the difference expected by hydrostatic

1458 ! balance. This test gets rid of pressure differences that are small, e.g. open ocean.

1459 ! Use 5 point quadrature to calculate intypa

1460 t5(1) = t_top(i,j) ; t5(5) = t_top(i,j+1)

1461 s5(1) = s_top(i,j) ; s5(5) = s_top(i,j+1)

1462 pa5(1) = pa(i,j,1) ; pa5(5) = pa(i,j+1,1)

1463 ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines.

1464 p5(1) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1465 p5(5) = -gxrho0*(e(i,j+1,1) - z_0p(i,j))

1466

1467 do m=2,4

1468 wt_r = 0.25*real(m-1)

1469 t5(m) = t5(1) + (t5(5)-t5(1))*wt_r

1470 s5(m) = s5(1) + (s5(5)-s5(1))*wt_r

1471 p5(m) = p5(1) + (p5(5)-p5(1))*wt_r

1472 enddo !m

1473 call calculate_density(t5, s5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)

1474

1475 ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure

1476 ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule

1477 ! quadrature to find the integrated correction to the integral of pressure along the interface.

1478 inty_pa_cor(i,j) = c1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc

1479

1480 ! The derivation of this correction follows:

1481

1482 ! Make pressure curvature a difference from the linear fit of pressure between the two points

1483 ! (which is equivalent to taking 4 trapezoidal rule integrals of the hydrostatic equation on

1484 ! sub-segments), with a constant slope that is chosen so that the pressure anomalies at the

1485 ! two ends of the segment agree with their known values.

1486 ! d_geo_8 = 0.125*dz_geo_sfc

1487 ! dpa_subseg = 0.25*(pa5(5)-pa5(1)) + &

1488 ! 0.25*d_geo_8 * ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3)))

1489 ! do m=2,4

1490 ! pa5(m) = pa5(m-1) + dpa_subseg - d_geo_8*(r5(m)+r5(m-1)))

1491 ! enddo

1492

1493 ! Explicitly finding expressions for the incremental pressures from the recursion relation above:

1494 ! pa5(2) = 0.25*(3.*pa5(1) + pa5(5)) + 0.25*d_geo_8 * ( (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) )

1495 ! ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + 0.25*d_geo_8 * &

1496 ! ! ( (r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3)) + &

1497 ! ! (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) - 4.*(r5(3)+r5(2)) )

1498 ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + d_geo_8 * (0.5*(r5(5)-r5(1)) + (r5(4)-r5(2)) )

1499 ! ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * &

1500 ! ! (2.0*(r5(5)-r5(1)) + 4.0*(r5(4)-r5(2)) + (r5(5)+r5(1)) + &

1501 ! ! 2.0*(r5(4)+r5(2)) + 2.0*r5(3) - 4.*(r5(4)+r5(3)))

1502 ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) )

1503 ! ! pa5(5) = pa5(5) + 0.25*d_geo_8 * &

1504 ! ! ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) + &

1505 ! ! ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3))) - 4.*(r5(5)+r5(4)) )

1506 ! pa5(5) = pa5(5) ! As it should.

1507

1508 ! From these:

1509 ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + 0.25*d_geo_8 * &

1510 ! ( (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) + (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3))

1511 ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + d_geo_8 * ( (r5(5)-r5(1)) + (r5(4)-r5(2)) )

1512

1513 ! Get the correction from the difference between the 5-point quadrature integral of pa5 and

1514 ! its trapezoidal rule integral as:

1515 ! inty_pa_cor(i,J) = C1_90*(7.0*(pa5(1)+pa5(5)) + 32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 0.5*(pa5(1)+pa5(5)))

1516 ! inty_pa_cor(i,J) = C1_90*((32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 38.0*(pa5(1)+pa5(5)))

1517 ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ((32.0*( (r5(5)-r5(1)) + (r5(4)-r5(2)) ) + &

1518 ! (6.*(r5(5)-r5(1)) + 12.0*(r5(4)-r5(2)) ))

1519 ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ( 38.0*(r5(5)-r5(1)) + 44.0*(r5(4)-r5(2)) )

1520 endif

1521 endif

1522 inty_pa(i,j,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) + inty_pa_cor(i,j)

1523 enddo ; enddo

1524

1525 if (cs%debug) then

1526 call uvchksum("int[xy]_pa_cor", intx_pa_cor, inty_pa_cor, g%HI, haloshift=0, &

1527 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%RL2_T2_to_Pa)

1528 call uvchksum("int[xy]_pa(1)", intx_pa(:,:,1), inty_pa(:,:,1), g%HI, haloshift=0, &

1529 symmetric=g%Domain%symmetric, scalar_pair=.true., unscale=us%RL2_T2_to_Pa)

1530 endif

1531

1532 else

1533 ! Set the surface boundary conditions on the horizontally integrated pressure anomaly,

1534 ! assuming that the surface pressure anomaly varies linearly in x and y.

1535 ! If there is an ice-shelf or icebergs, this linear variation would need to be applied

1536 ! to an interior interface.

1537 !$OMP parallel do default(shared)

1538 do j=js,je ; do i=isq,ieq

1539 intx_pa(i,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1))

1540 enddo ; enddo

1541 !$OMP parallel do default(shared)

1542 do j=jsq,jeq ; do i=is,ie

1543 inty_pa(i,j,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1))

1544 enddo ; enddo

1545 endif

1546

1547 do k=1,nz

1548 !$OMP parallel do default(shared)

1549 do j=js,je ; do i=isq,ieq

1550 intx_pa(i,j,k+1) = intx_pa(i,j,k) + intx_dpa(i,j,k)

1551 enddo ; enddo

1552 enddo

1553 do k=1,nz

1554 !$OMP parallel do default(shared)

1555 do j=jsq,jeq ; do i=is,ie

1556 inty_pa(i,j,k+1) = inty_pa(i,j,k) + inty_dpa(i,j,k)

1557 enddo ; enddo

1558 enddo

1559

1560 if (cs%reset_intxpa_integral) then

1561 ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is

1562 ! reset intxpa there, then adjust intxpa throughout the water column.

1563

1564 ! Zero out the 2-d arrays that will be set from various reference interfaces.

1565 t_int_w(:,:) = 0.0 ; s_int_w(:,:) = 0.0 ; p_int_w(:,:) = 0.0

1566 t_int_e(:,:) = 0.0 ; s_int_e(:,:) = 0.0 ; p_int_e(:,:) = 0.0

1567 intx_pa_nonlin(:,:) = 0.0 ; dgeo_x(:,:) = 0.0 ; intx_pa_cor_ri(:,:) = 0.0

1568 do j=js,je ; do i=isq,ieq

1569 seek_x_cor(i,j) = (g%mask2dCu(i,j) > 0.)

1570 delta_z_x(i,j) = 0.0

1571 enddo ; enddo

1572

1573 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

1574 if ((e(i+1,j,2) <= e(i,j,1)) .and. (e(i,j,2) <= e(i+1,j,1))) then

1575 ! This is a typical case in the open ocean, so use the topmost interface.

1576 t_int_w(i,j) = t_top(i,j) ; t_int_e(i,j) = t_top(i+1,j)

1577 s_int_w(i,j) = s_top(i,j) ; s_int_e(i,j) = s_top(i+1,j)

1578 p_int_w(i,j) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1579 p_int_e(i,j) = -gxrho0*(e(i+1,j,1) - z_0p(i,j))

1580 intx_pa_nonlin(i,j) = intx_pa(i,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))

1581 dgeo_x(i,j) = gv%g_Earth * (e(i+1,j,1)-e(i,j,1))

1582 seek_x_cor(i,j) = .false.

1583 endif

1584 endif ; enddo ; enddo

1585

1586 do k=1,nz

1587 do_more_k = .false.

1588 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

1589 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not

1590 ! activated in the subgrid interpolation.

1591 if (((h(i,j,k) > cs%h_nonvanished) .and. (h(i+1,j,k) > cs%h_nonvanished)) .and. &

1592 (max(0., e(i+1,j,k+1)-e(i,j,1), e(i,j,k+1)-e(i+1,j,1)) <= 0.0)) then

1593 ! Store properties at the bottom of this cell to get a "good estimate" for intxpa at

1594 ! the interface below this cell (it might have quadratic pressure dependence if sloped)

1595 t_int_w(i,j) = t_b(i,j,k) ; t_int_e(i,j) = t_b(i+1,j,k)

1596 s_int_w(i,j) = s_b(i,j,k) ; s_int_e(i,j) = s_b(i+1,j,k)

1597 ! These pressures are only used for the equation of state, and are only a function of

1598 ! height, consistent with the expressions in the int_density_dz routines.

1599 p_int_w(i,j) = -gxrho0*(e(i,j,k+1) - z_0p(i,j))

1600 p_int_e(i,j) = -gxrho0*(e(i+1,j,k+1) - z_0p(i,j))

1601

1602 intx_pa_nonlin(i,j) = intx_pa(i,j,k+1) - 0.5*(pa(i,j,k+1) + pa(i+1,j,k+1))

1603 dgeo_x(i,j) = gv%g_Earth * (e(i+1,j,k+1)-e(i,j,k+1))

1604 seek_x_cor(i,j) = .false.

1605 else

1606 do_more_k = .true.

1607 endif

1608 endif ; enddo ; enddo

1609 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.

1610 enddo

1611

1612 if (do_more_k) then

1613 if (cs%reset_intxpa_flattest) then

1614 ! There are still points where a correction is needed, so use flattest interface

1615 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

1616 ! choose top layer first

1617 t_int_w(i,j) = t_top(i,j) ; t_int_e(i,j) = t_top(i+1,j)

1618 s_int_w(i,j) = s_top(i,j) ; s_int_e(i,j) = s_top(i+1,j)

1619 p_int_w(i,j) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1620 p_int_e(i,j) = -gxrho0*(e(i+1,j,1) - z_0p(i,j))

1621 intx_pa_nonlin(i,j) = intx_pa(i,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))

1622 dgeo_x(i,j) = gv%g_Earth * (e(i+1,j,1)-e(i,j,1))

1623 delta_z_x(i,j) = abs(e(i+1,j,1)-e(i,j,1))

1624 do k=1,nz

1625 if (abs(e(i+1,j,k+1)-e(i,j,k+1)) < delta_z_x(i,j)) then

1626 ! bottom of layer is less sloped than top. Use this layer

1627 delta_z_x(i,j) = abs(e(i+1,j,k+1)-e(i,j,k+1))

1628 t_int_w(i,j) = t_b(i,j,k) ; t_int_e(i,j) = t_b(i+1,j,k)

1629 s_int_w(i,j) = s_b(i,j,k) ; s_int_e(i,j) = s_b(i+1,j,k)

1630 p_int_w(i,j) = -gxrho0*(e(i,j,k+1) - z_0p(i,j))

1631 p_int_e(i,j) = -gxrho0*(e(i+1,j,k+1) - z_0p(i,j))

1632 intx_pa_nonlin(i,j) = intx_pa(i,j,k+1) - 0.5*(pa(i,j,k+1) + pa(i+1,j,k+1))

1633 dgeo_x(i,j) = gv%g_Earth * (e(i+1,j,k+1)-e(i,j,k+1))

1634 endif

1635 enddo

1636 seek_x_cor(i,j) = .false.

1637 endif ; enddo ; enddo

1638 else

1639 ! There are still points where a correction is needed, so use the top interface for lack of a better idea?

1640 do j=js,je ; do i=isq,ieq ; if (seek_x_cor(i,j)) then

1641 t_int_w(i,j) = t_top(i,j) ; t_int_e(i,j) = t_top(i+1,j)

1642 s_int_w(i,j) = s_top(i,j) ; s_int_e(i,j) = s_top(i+1,j)

1643 p_int_w(i,j) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1644 p_int_e(i,j) = -gxrho0*(e(i+1,j,1) - z_0p(i,j))

1645 intx_pa_nonlin(i,j) = intx_pa(i,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))

1646 dgeo_x(i,j) = gv%g_Earth * (e(i+1,j,1)-e(i,j,1))

1647 seek_x_cor(i,j) = .false.

1648 endif ; enddo ; enddo

1649 endif

1650 endif

1651

1652 do j=js,je

1653 do i=isq,ieq

1654 ! This expression assumes that temperature and salinity vary linearly with hieght

1655 ! between the corners of the reference interfaces found above to get a correction to

1656 ! intx_pa that takes nonlinearities in the equation of state into account.

1657 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale

1658 ! linear correction so that the pressures and heights match at the end-points. It turns

1659 ! out that this linear correction cancels out the mid-point density anomaly.

1660 ! This can be used without masking because dgeo_x and intx_pa_nonlin are 0 over land.

1661 t5(1) = t_int_w(i,j) ; s5(1) = s_int_w(i,j) ; p5(1) = p_int_w(i,j)

1662 t5(5) = t_int_e(i,j) ; s5(5) = s_int_e(i,j) ; p5(5) = p_int_e(i,j)

1663 t5(2) = 0.25*(3.0*t5(1) + t5(5)) ; t5(4) = 0.25*(3.0*t5(5) + t5(1)) ; t5(3) = 0.5*(t5(5) + t5(1))

1664 s5(2) = 0.25*(3.0*s5(1) + s5(5)) ; s5(4) = 0.25*(3.0*s5(5) + s5(1)) ; s5(3) = 0.5*(s5(5) + s5(1))

1665 p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1))

1666 call calculate_density(t5, s5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)

1667

1668 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12

1669 intx_pa_cor_ri(i,j) = c1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_x(i,j) - &

1670 intx_pa_nonlin(i,j)

1671 enddo

1672 enddo

1673

1674 ! Repeat the calculations above for v-velocity points.

1675 t_int_s(:,:) = 0.0 ; s_int_s(:,:) = 0.0 ; p_int_s(:,:) = 0.0

1676 t_int_n(:,:) = 0.0 ; s_int_n(:,:) = 0.0 ; p_int_n(:,:) = 0.0

1677 inty_pa_nonlin(:,:) = 0.0 ; dgeo_y(:,:) = 0.0 ; inty_pa_cor_ri(:,:) = 0.0

1678 do j=jsq,jeq ; do i=is,ie

1679 seek_y_cor(i,j) = (g%mask2dCv(i,j) > 0.)

1680 delta_z_y(i,j) = 0.0

1681 enddo ; enddo

1682

1683 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

1684 if ((e(i,j+1,2) <= e(i,j,1)) .and. (e(i,j,2) <= e(i,j+1,1))) then

1685 ! This is a typical case in the open ocean, so use the topmost interface.

1686 t_int_s(i,j) = t_top(i,j) ; t_int_n(i,j) = t_top(i,j+1)

1687 s_int_s(i,j) = s_top(i,j) ; s_int_n(i,j) = s_top(i,j+1)

1688 p_int_s(i,j) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1689 p_int_n(i,j) = -gxrho0*(e(i,j+1,1) - z_0p(i,j))

1690 inty_pa_nonlin(i,j) = inty_pa(i,j,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))

1691 dgeo_y(i,j) = gv%g_Earth * (e(i,j+1,1)-e(i,j,1))

1692 seek_y_cor(i,j) = .false.

1693 endif

1694 endif ; enddo ; enddo

1695

1696 do k=1,nz

1697 do_more_k = .false.

1698 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

1699 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not

1700 ! activated in the subgrid interpolation.

1701 if (((h(i,j,k) > cs%h_nonvanished) .and. (h(i,j+1,k) > cs%h_nonvanished)) .and. &

1702 (max(0., e(i,j+1,k+1)-e(i,j,1), e(i,j,k+1)-e(i,j+1,1)) <= 0.0)) then

1703 ! Store properties at the bottom of this cell to get a "good estimate" for intypa at

1704 ! the interface below this cell (it might have quadratic pressure dependence if sloped)

1705 t_int_s(i,j) = t_b(i,j,k) ; t_int_n(i,j) = t_b(i,j+1,k)

1706 s_int_s(i,j) = s_b(i,j,k) ; s_int_n(i,j) = s_b(i,j+1,k)

1707 ! These pressures are only used for the equation of state, and are only a function of

1708 ! height, consistent with the expressions in the int_density_dz routines.

1709 p_int_s(i,j) = -gxrho0*(e(i,j,k+1) - z_0p(i,j))

1710 p_int_n(i,j) = -gxrho0*(e(i,j+1,k+1) - z_0p(i,j))

1711 inty_pa_nonlin(i,j) = inty_pa(i,j,k+1) - 0.5*(pa(i,j,k+1) + pa(i,j+1,k+1))

1712 dgeo_y(i,j) = gv%g_Earth * (e(i,j+1,k+1)-e(i,j,k+1))

1713 seek_y_cor(i,j) = .false.

1714 else

1715 do_more_k = .true.

1716 endif

1717 endif ; enddo ; enddo

1718 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.

1719 enddo

1720

1721 if (do_more_k) then

1722 if (cs%reset_intxpa_flattest) then

1723 ! There are still points where a correction is needed, so use flattest interface.

1724 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

1725 ! choose top interface first

1726 t_int_s(i,j) = t_top(i,j) ; t_int_n(i,j) = t_top(i,j+1)

1727 s_int_s(i,j) = s_top(i,j) ; s_int_n(i,j) = s_top(i,j+1)

1728 p_int_s(i,j) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1729 p_int_n(i,j) = -gxrho0*(e(i,j+1,1) - z_0p(i,j))

1730 inty_pa_nonlin(i,j) = inty_pa(i,j,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))

1731 dgeo_y(i,j) = gv%g_Earth * (e(i,j+1,1)-e(i,j,1))

1732 delta_z_y(i,j) = abs(e(i,j+1,1)-e(i,j,1))

1733 do k=1,nz

1734 if (abs(e(i,j+1,k+1)-e(i,j,k+1)) < delta_z_y(i,j)) then

1735 ! bottom of layer is less sloped than top. Use this layer

1736 delta_z_y(i,j) = abs(e(i,j+1,k+1)-e(i,j,k+1))

1737 t_int_s(i,j) = t_b(i,j,k) ; t_int_n(i,j) = t_b(i,j+1,k)

1738 s_int_s(i,j) = s_b(i,j,k) ; s_int_n(i,j) = s_b(i,j+1,k)

1739 p_int_s(i,j) = -gxrho0*(e(i,j,k+1) - z_0p(i,j))

1740 p_int_n(i,j) = -gxrho0*(e(i,j+1,k+1) - z_0p(i,j))

1741 inty_pa_nonlin(i,j) = inty_pa(i,j,k+1) - 0.5*(pa(i,j,k+1) + pa(i,j+1,k+1))

1742 dgeo_y(i,j) = gv%g_Earth * (e(i,j+1,k+1)-e(i,j,k+1))

1743 endif

1744 enddo

1745 seek_y_cor(i,j) = .false.

1746 endif ; enddo ; enddo

1747 else

1748 ! There are still points where a correction is needed, so use the top interface for lack of a better idea?

1749 do j=jsq,jeq ; do i=is,ie ; if (seek_y_cor(i,j)) then

1750 t_int_s(i,j) = t_top(i,j) ; t_int_n(i,j) = t_top(i,j+1)

1751 s_int_s(i,j) = s_top(i,j) ; s_int_n(i,j) = s_top(i,j+1)

1752 p_int_s(i,j) = -gxrho0*(e(i,j,1) - z_0p(i,j))

1753 p_int_n(i,j) = -gxrho0*(e(i,j+1,1) - z_0p(i,j))

1754 inty_pa_nonlin(i,j) = inty_pa(i,j,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))

1755 dgeo_y(i,j) = gv%g_Earth * (e(i,j+1,1)-e(i,j,1))

1756 seek_y_cor(i,j) = .false.

1757 endif ; enddo ; enddo

1758 endif

1759 endif

1760

1761 do j=jsq,jeq

1762 do i=is,ie

1763 ! This expression assumes that temperature and salinity vary linearly with hieght

1764 ! between the corners of the reference interfaces found above to get a correction to

1765 ! intx_pa that takes nonlinearities in the equation of state into account.

1766 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale

1767 ! linear correction so that the pressures and heights match at the end-points. It turns

1768 ! out that this linear correction cancels out the mid-point density anomaly.

1769 ! This can be used without masking because dgeo_y and inty_pa_nonlin are 0 over land.

1770 t5(1) = t_int_s(i,j) ; s5(1) = s_int_s(i,j) ; p5(1) = p_int_s(i,j)

1771 t5(5) = t_int_n(i,j) ; s5(5) = s_int_n(i,j) ; p5(5) = p_int_n(i,j)

1772 t5(2) = 0.25*(3.0*t5(1) + t5(5)) ; t5(4) = 0.25*(3.0*t5(5) + t5(1)) ; t5(3) = 0.5*(t5(5) + t5(1))

1773 s5(2) = 0.25*(3.0*s5(1) + s5(5)) ; s5(4) = 0.25*(3.0*s5(5) + s5(1)) ; s5(3) = 0.5*(s5(5) + s5(1))

1774 p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1))

1775 call calculate_density(t5, s5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)

1776

1777 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12

1778 inty_pa_cor_ri(i,j) = c1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_y(i,j) - &

1779 inty_pa_nonlin(i,j)

1780 enddo

1781 enddo

1782

1783 ! Correct intx_pa and inty_pa at each interface using vertically constant corrections.

1784 do k=1,nz+1 ; do j=js,je ; do i=isq,ieq

1785 intx_pa(i,j,k) = intx_pa(i,j,k) + intx_pa_cor_ri(i,j)

1786 enddo ; enddo ; enddo

1787

1788 do k=1,nz+1 ; do j=jsq,jeq ; do i=is,ie

1789 inty_pa(i,j,k) = inty_pa(i,j,k) + inty_pa_cor_ri(i,j)

1790 enddo ; enddo ; enddo

1791 endif ! intx_pa and inty_pa have now been reset to reflect the properties of an unimpeded interface.

1792

1793 ! Compute pressure gradient in x direction

1794 !$OMP parallel do default(shared)

1795 do k=1,nz ; do j=js,je ; do i=isq,ieq

1796 pfu(i,j,k) = (((pa(i,j,k)*h(i,j,k) + intz_dpa(i,j,k)) - &

1797 (pa(i+1,j,k)*h(i+1,j,k) + intz_dpa(i+1,j,k))) + &

1798 ((h(i+1,j,k) - h(i,j,k)) * intx_pa(i,j,k) - &

1799 (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa(i,j,k) * gv%Z_to_H)) * &

1800 ((2.0*i_rho0*g%IdxCu(i,j)) / &

1801 ((h(i,j,k) + h(i+1,j,k)) + h_neglect))

1802 enddo ; enddo ; enddo

1803

1804 ! Compute pressure gradient in y direction

1805 !$OMP parallel do default(shared)

1806 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

1807 pfv(i,j,k) = (((pa(i,j,k)*h(i,j,k) + intz_dpa(i,j,k)) - &

1808 (pa(i,j+1,k)*h(i,j+1,k) + intz_dpa(i,j+1,k))) + &

1809 ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,j,k) - &

1810 (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa(i,j,k) * gv%Z_to_H)) * &

1811 ((2.0*i_rho0*g%IdyCv(i,j)) / &

1812 ((h(i,j,k) + h(i,j+1,k)) + h_neglect))

1813 enddo ; enddo ; enddo

1814

1815 ! Calculate SAL geopotential anomaly and add its gradient to pressure gradient force

1816 if (cs%calculate_SAL .and. cs%tides_answer_date>20230630 .and. cs%bq_sal_tides) then

1817 !$OMP parallel do default(shared)

1818 do k=1,nz

1819 do j=js,je ; do i=isq,ieq

1820 pfu(i,j,k) = pfu(i,j,k) + (e_sal(i+1,j) - e_sal(i,j)) * gv%g_Earth * g%IdxCu(i,j)

1821 enddo ; enddo

1822 do j=jsq,jeq ; do i=is,ie

1823 pfv(i,j,k) = pfv(i,j,k) + (e_sal(i,j+1) - e_sal(i,j)) * gv%g_Earth * g%IdyCv(i,j)

1824 enddo ; enddo

1825 enddo

1826 endif

1827

1828 ! Calculate tidal geopotential anomaly and add its gradient to pressure gradient force

1829 if (cs%tides .and. cs%tides_answer_date>20230630 .and. cs%bq_sal_tides) then

1830 !$OMP parallel do default(shared)

1831 do k=1,nz

1832 do j=js,je ; do i=isq,ieq

1833 pfu(i,j,k) = pfu(i,j,k) + ((e_tidal_eq(i+1,j) + e_tidal_sal(i+1,j)) &

1834 - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * gv%g_Earth * g%IdxCu(i,j)

1835 enddo ; enddo

1836 do j=jsq,jeq ; do i=is,ie

1837 pfv(i,j,k) = pfv(i,j,k) + ((e_tidal_eq(i,j+1) + e_tidal_sal(i,j+1)) &

1838 - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * gv%g_Earth * g%IdyCv(i,j)

1839 enddo ; enddo

1840 enddo

1841 endif

1842

1843 if (cs%GFS_scale < 1.0) then

1844 ! Adjust the Montgomery potential to make this a reduced gravity model.

1845 if (use_eos) then

1846 !$OMP parallel do default(shared)

1847 do j=jsq,jeq+1

1848 if (use_p_atm) then

1849 call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, &

1850 tv%eqn_of_state, eosdom)

1851 else

1852 call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, &

1853 tv%eqn_of_state, eosdom)

1854 endif

1855 do i=isq,ieq+1

1856 dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * rho_in_situ(i)) * (e(i,j,1) - g%Z_ref)

1857 enddo

1858 enddo

1859 else

1860 !$OMP parallel do default(shared)

1861 do j=jsq,jeq+1 ; do i=isq,ieq+1

1862 dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * gv%Rlay(1)) * (e(i,j,1) - g%Z_ref)

1863 enddo ; enddo

1864 endif

1865

1866 !$OMP parallel do default(shared)

1867 do k=1,nz

1868 do j=js,je ; do i=isq,ieq

1869 pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)

1870 enddo ; enddo

1871 do j=jsq,jeq ; do i=is,ie

1872 pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)

1873 enddo ; enddo

1874 enddo

1875 endif

1876

1877 if (present(pbce)) then

1878 call set_pbce_bouss(e, tv_tmp, g, gv, us, rho0_set_pbce, cs%GFS_scale, pbce)

1879 endif

1880

1881 if (present(eta)) then

1882 ! eta is the sea surface height relative to a time-invariant geoid, for comparison with

1883 ! what is used for eta in btstep. See how e was calculated about 200 lines above.

1884 !$OMP parallel do default(shared)

1885 do j=jsq,jeq+1 ; do i=isq,ieq+1

1886 eta(i,j) = e(i,j,1)*gv%Z_to_H

1887 enddo ; enddo

1888 if (cs%tides .and. (.not.cs%bq_sal_tides)) then

1889 if (cs%tides_answer_date>20230630) then

1890 !$OMP parallel do default(shared)

1891 do j=jsq,jeq+1 ; do i=isq,ieq+1

1892 eta(i,j) = eta(i,j) + (e_tidal_eq(i,j)+e_tidal_sal(i,j))*gv%Z_to_H

1893 enddo ; enddo

1894 else

1895 !$OMP parallel do default(shared)

1896 do j=jsq,jeq+1 ; do i=isq,ieq+1

1897 eta(i,j) = eta(i,j) + e_sal_and_tide(i,j)*gv%Z_to_H

1898 enddo ; enddo

1899 endif

1900 endif

1901 if (cs%calculate_SAL .and. (cs%tides_answer_date>20230630) .and. (.not.cs%bq_sal_tides)) then

1902 !$OMP parallel do default(shared)

1903 do j=jsq,jeq+1 ; do i=isq,ieq+1

1904 eta(i,j) = eta(i,j) + e_sal(i,j)*gv%Z_to_H

1905 enddo ; enddo

1906 endif

1907 endif

1908

1909 if (cs%use_stanley_pgf) then

1910 ! Calculated diagnostics related to the Stanley parameterization

1911 zeros(:) = 0.0

1912 eosdom_h(:) = eos_domain(g%HI)

1913 if ((cs%id_p_stanley>0) .or. (cs%id_rho_pgf>0) .or. (cs%id_rho_stanley_pgf>0)) then

1914 ! Find the pressure at the mid-point of each layer.

1915 h_to_rl2_t2 = gv%g_Earth*gv%H_to_RZ

1916 if (use_p_atm) then

1917 do j=js,je ; do i=is,ie

1918 p_stanley(i,j,1) = 0.5*h(i,j,1) * h_to_rl2_t2 + p_atm(i,j)

1919 enddo ; enddo

1920 else

1921 do j=js,je ; do i=is,ie

1922 p_stanley(i,j,1) = 0.5*h(i,j,1) * h_to_rl2_t2

1923 enddo ; enddo

1924 endif

1925 do k=2,nz ; do j=js,je ; do i=is,ie

1926 p_stanley(i,j,k) = p_stanley(i,j,k-1) + 0.5*(h(i,j,k-1) + h(i,j,k)) * h_to_rl2_t2

1927 enddo ; enddo ; enddo

1928 endif

1929 if (cs%id_p_stanley>0) call post_data(cs%id_p_stanley, p_stanley, cs%diag)

1930 if (cs%id_rho_pgf>0) then

1931 do k=1,nz ; do j=js,je

1932 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_stanley(:,j,k), zeros, &

1933 zeros, zeros, rho_pgf(:,j,k), tv%eqn_of_state, eosdom_h)

1934 enddo ; enddo

1935 call post_data(cs%id_rho_pgf, rho_pgf, cs%diag)

1936 endif

1937 if (cs%id_rho_stanley_pgf>0) then

1938 do k=1,nz ; do j=js,je

1939 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_stanley(:,j,k), tv%varT(:,j,k), &

1940 zeros, zeros, rho_stanley_pgf(:,j,k), tv%eqn_of_state, eosdom_h)

1941 enddo ; enddo

1942 call post_data(cs%id_rho_stanley_pgf, rho_stanley_pgf, cs%diag)

1943 endif

1944 endif

1945

1946 if (cs%id_MassWt_u>0) call post_data(cs%id_MassWt_u, masswt_u, cs%diag)

1947 if (cs%id_MassWt_v>0) call post_data(cs%id_MassWt_v, masswt_v, cs%diag)

1948

1949 if (cs%id_rho_pgf>0) call post_data(cs%id_rho_pgf, rho_pgf, cs%diag)

1950 if (cs%id_rho_stanley_pgf>0) call post_data(cs%id_rho_stanley_pgf, rho_stanley_pgf, cs%diag)

1951 if (cs%id_p_stanley>0) call post_data(cs%id_p_stanley, p_stanley, cs%diag)

1952

1953 ! Diagnostics for tidal forcing and SAL height anomaly

1954 if (cs%id_e_tide>0) then

1955 ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL.

1956 ! New diagnostics are given for each individual field.

1957 if (cs%tides_answer_date>20230630) then ; do j=jsq,jeq+1 ; do i=isq,ieq+1

1958 e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j)

1959 enddo ; enddo ; endif

1960 call post_data(cs%id_e_tide, e_sal_and_tide, cs%diag)

1961 endif

1962 if (cs%id_e_sal>0) call post_data(cs%id_e_sal, e_sal, cs%diag)

1963 if (cs%id_e_tidal_eq>0) call post_data(cs%id_e_tidal_eq, e_tidal_eq, cs%diag)

1964 if (cs%id_e_tidal_sal>0) call post_data(cs%id_e_tidal_sal, e_tidal_sal, cs%diag)

1965

1966 ! Diagnostics for tidal forcing and SAL horizontal gradients

1967 if (cs%calculate_SAL .and. ((associated(adp%sal_u) .or. associated(adp%sal_v)))) then

1968 if (cs%tides) then ; do j=jsq,jeq+1 ; do i=isq,ieq+1

1969 e_sal(i,j) = e_sal(i,j) + e_tidal_sal(i,j)

1970 enddo ; enddo ; endif

1971 if (cs%bq_sal_tides) then

1972 ! sal_u = ( e(i+1) - e(i) ) * g / dx

1973 if (associated(adp%sal_u)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq

1974 adp%sal_u(i,j,k) = (e_sal(i+1,j) - e_sal(i,j)) * gv%g_Earth * g%IdxCu(i,j)

1975 enddo ; enddo ; enddo ; endif

1976 if (associated(adp%sal_v)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie

1977 adp%sal_v(i,j,k) = (e_sal(i,j+1) - e_sal(i,j)) * gv%g_Earth * g%IdyCv(i,j)

1978 enddo ; enddo ; enddo ; endif

1979 else

1980 ! sal_u = ( e(i+1) - e(i) ) * g / dx * (rho(k) / rho0)

1981 if (associated(adp%sal_u)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq

1982 adp%sal_u(i,j,k) = (e_sal(i+1,j) - e_sal(i,j)) * g%IdxCu(i,j) * i_rho0 * &

1983 (2.0 * intx_dpa(i,j,k) * gv%Z_to_H / ((h(i,j,k) + h(i+1,j,k)) + h_neglect) + gxrho_ref)

1984 enddo ; enddo ; enddo ; endif

1985 if (associated(adp%sal_v)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie

1986 adp%sal_v(i,j,k) = (e_sal(i,j+1) - e_sal(i,j)) * g%IdyCv(i,j) * i_rho0 * &

1987 (2.0 * inty_dpa(i,j,k) * gv%Z_to_H / ((h(i,j,k) + h(i,j+1,k)) + h_neglect) + gxrho_ref)

1988 enddo ; enddo ; enddo ; endif

1989 endif

1990 if (cs%id_sal_u>0) call post_data(cs%id_sal_u, adp%sal_u, cs%diag)

1991 if (cs%id_sal_v>0) call post_data(cs%id_sal_v, adp%sal_v, cs%diag)

1992 endif

1993

1994 if (cs%tides .and. ((associated(adp%tides_u) .or. associated(adp%tides_v)))) then

1995 if (cs%bq_sal_tides) then

1996 ! tides_u = ( e(i+1) - e(i) ) * g / dx

1997 if (associated(adp%tides_u)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq

1998 adp%tides_u(i,j,k) = (e_tidal_eq(i+1,j) - e_tidal_eq(i,j)) * gv%g_Earth * g%IdxCu(i,j)

1999 enddo ; enddo ; enddo ; endif

2000 if (associated(adp%tides_v)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie

2001 adp%tides_v(i,j,k) = (e_tidal_eq(i,j+1) - e_tidal_eq(i,j)) * gv%g_Earth * g%IdyCv(i,j)

2002 enddo ; enddo ; enddo ; endif

2003 else

2004 ! tides_u = ( e(i+1) - e(i) ) * g / dx * (rho(k) / rho0)

2005 if (associated(adp%tides_u)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq

2006 adp%tides_u(i,j,k) = (e_tidal_eq(i+1,j) - e_tidal_eq(i,j)) * g%IdxCu(i,j) * i_rho0 * &

2007 (2.0 * intx_dpa(i,j,k) * gv%Z_to_H / ((h(i,j,k) + h(i+1,j,k)) + h_neglect) + gxrho_ref)

2008 enddo ; enddo ; enddo ; endif

2009 if (associated(adp%tides_v)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie

2010 adp%tides_v(i,j,k) = (e_tidal_eq(i,j+1) - e_tidal_eq(i,j)) * g%IdyCv(i,j) * i_rho0 * &

2011 (2.0 * inty_dpa(i,j,k) * gv%Z_to_H / ((h(i,j,k) + h(i,j+1,k)) + h_neglect) + gxrho_ref)

2012 enddo ; enddo ; enddo ; endif

2013 endif

2014 if (cs%id_tides_u>0) call post_data(cs%id_tides_u, adp%tides_u, cs%diag)

2015 if (cs%id_tides_v>0) call post_data(cs%id_tides_v, adp%tides_v, cs%diag)

2016 endif

2017end subroutine pressureforce_fv_bouss

2018

2019!> Initializes the finite volume pressure gradient control structure

2020subroutine pressureforce_fv_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL_CSp, tides_CSp)

2021 type(time_type), target, intent(in) :: time !< Current model time

2022 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

2023 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure

2024 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2025 type(param_file_type), intent(in) :: param_file !< Parameter file handles

2026 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure

2027 type(pressureforce_fv_cs), intent(inout) :: cs !< Finite volume PGF control structure

2028 type(accel_diag_ptrs), pointer :: adp !< Acceleration diagnostic pointers

2029 type(sal_cs), intent(in), target, optional :: sal_csp !< SAL control structure

2030 type(tidal_forcing_cs), intent(in), target, optional :: tides_csp !< Tides control structure

2031

2032 ! Local variables

2033 real :: stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale

2034 ! temperature variance [nondim]

2035 integer :: default_answer_date ! Global answer date

2036 logical :: use_temperature ! If true, temperature and salinity are used as state variables.

2037 logical :: use_eos ! If true, density calculated from T & S using an equation of state.

2038 logical :: usemasswghtinterp ! If true, use near-bottom mass weighting for T and S

2039 logical :: masswghtinterptop ! If true, use near-surface mass weighting for T and S under ice shelves

2040 logical :: masswghtinterp_nonbous_bug ! If true, use a buggy mass weighting when non-Boussinesq

2041 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to

2042 ! recreate the bugs, or if false bugs are only used if actively selected.

2043 ! This include declares and sets the variable "version".

2044# include "version_variable.h"

2045 character(len=40) :: mdl ! This module's name.

2046 logical :: use_ale ! If true, use the Vertical Lagrangian Remap algorithm

2047 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz

2048

2049 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke

2050 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

2051

2052 cs%initialized = .true.

2053 cs%diag => diag ; cs%Time => time

2054 if (present(tides_csp)) &

2055 cs%tides_CSp => tides_csp

2056 if (present(sal_csp)) &

2057 cs%SAL_CSp => sal_csp

2058

2059 mdl = "MOM_PressureForce_FV"

2060 call log_version(param_file, mdl, version, "")

2061 call get_param(param_file, mdl, "DEBUG", cs%debug, &

2062 "If true, write out verbose debugging data.", &

2063 default=.false., debuggingparam=.true., do_not_log=.true.)

2064 call get_param(param_file, mdl, "RHO_PGF_REF", cs%rho_ref, &

2065 "The reference density that is subtracted off when calculating pressure "//&

2066 "gradient forces. Its inverse is subtracted off of specific volumes when "//&

2067 "in non-Boussinesq mode. The default is RHO_0.", &

2068 units="kg m-3", default=gv%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R)

2069 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &

2070 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.

2071 call get_param(param_file, mdl, "RHO_PGF_REF_BUG", cs%rho_ref_bug, &

2072 "If true, recover a bug that RHO_0 (the mean seawater density in Boussinesq mode) "//&

2073 "and RHO_PGF_REF (the subtracted reference density in finite volume pressure "//&

2074 "gradient forces) are incorrectly interchanged in several instances in Boussinesq mode.", &

2075 default=enable_bugs)

2076 call get_param(param_file, mdl, "TIDES", cs%tides, &

2077 "If true, apply tidal momentum forcing.", default=.false.)

2078 if (cs%tides) then

2079 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &

2080 "This sets the default value for the various _ANSWER_DATE parameters.", &

2081 default=99991231)

2082 call get_param(param_file, mdl, "TIDES_ANSWER_DATE", cs%tides_answer_date, "The vintage of "//&

2083 "self-attraction and loading (SAL) and tidal forcing calculations. Setting "//&

2084 "dates before 20230701 recovers old answers (Boussinesq and non-Boussinesq "//&

2085 "modes) when SAL is part of the tidal forcing calculation. The answer "//&

2086 "difference is only at bit level and due to a reordered summation. Setting "//&

2087 "dates before 20250201 recovers answers (Boussinesq mode) that interface "//&

2088 "heights are modified before pressure force integrals are calculated.", &

2089 default=default_answer_date, do_not_log=(.not.cs%tides))

2090 endif

2091 call get_param(param_file, mdl, "CALCULATE_SAL", cs%calculate_SAL, &

2092 "If true, calculate self-attraction and loading.", default=cs%tides)

2093 if (cs%calculate_SAL) &

2094 call get_param(param_file, '', "SAL_USE_BPA", cs%sal_use_bpa, default=.false., &

2095 do_not_log=.true.)

2096 if ((cs%tides .or. cs%calculate_SAL) .and. gv%Boussinesq) &

2097 call get_param(param_file, mdl, "BOUSSINESQ_SAL_TIDES", cs%bq_sal_tides, "If true, "//&

2098 "in Boussinesq mode, use an alternative method to include self-attraction "//&

2099 "and loading (SAL) and tidal forcings in pressure gradient, in which their "//&

2100 "gradients are calculated separately, instead of adding geopotential "//&

2101 "anomalies as corrections to the interface height. This alternative method "//&

2102 "elimates a baroclinic component of the SAL and tidal forcings.", &

2103 default=.false.)

2104 call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &

2105 "If true, Temperature and salinity are used as state variables.", &

2106 default=.true., do_not_log=.true.)

2107 call get_param(param_file, "MOM", "USE_EOS", use_eos, &

2108 "If true, density is calculated from temperature and "//&

2109 "salinity with an equation of state. If USE_EOS is "//&

2110 "true, ENABLE_THERMODYNAMICS must be true as well.", &

2111 default=use_temperature, do_not_log=.true.)

2112

2113 call get_param(param_file, mdl, "SSH_IN_EOS_PRESSURE_FOR_PGF", cs%use_SSH_in_Z0p, &

2114 "If true, include contributions from the sea surface height in the height-based "//&

2115 "pressure used in the equation of state calculations for the Boussinesq pressure "//&

2116 "gradient forces, including adjustments for atmospheric or sea-ice pressure.", &

2117 default=.false., do_not_log=.not.gv%Boussinesq)

2118 if (cs%tides .and. cs%tides_answer_date<=20250131 .and. cs%use_SSH_in_Z0p) &

2119 call mom_error(fatal, trim(mdl) // ", PressureForce_FV_init: SSH_IN_EOS_PRESSURE_FOR_PGF "//&

2120 "needs to be FALSE to recover tide answers before 20250131.")

2121

2122 call get_param(param_file, "MOM", "USE_REGRIDDING", use_ale, &

2123 "If True, use the ALE algorithm (regridding/remapping). "//&

2124 "If False, use the layered isopycnal algorithm.", default=.false. )

2125 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", usemasswghtinterp, &

2126 "If true, use mass weighting when interpolating T/S for integrals "//&

2127 "near the bathymetry in FV pressure gradient calculations.", &

2128 default=.false.)

2129 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP", masswghtinterptop, &

2130 "If true and MASS_WEIGHT_IN_PRESSURE_GRADIENT is true, use mass weighting when "//&

2131 "interpolating T/S for integrals near the top of the water column in FV "//&

2132 "pressure gradient calculations. ", &

2133 default=usemasswghtinterp)

2134 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_NONBOUS_BUG", masswghtinterp_nonbous_bug, &

2135 "If true, use a masking bug in non-Boussinesq calculations with mass weighting "//&

2136 "when interpolating T/S for integrals near the bathymetry in FV pressure "//&

2137 "gradient calculations.", &

2138 default=.false., do_not_log=(gv%Boussinesq .or. (.not.usemasswghtinterp)))

2139 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_VANISHED_ONLY", cs%MassWghtInterpVanOnly, &

2140 "If true, use mass weighting when interpolating T/S for integrals "//&

2141 "only if one side is vanished according to RESET_INTXPA_H_NONVANISHED. ", &

2142 default=.false.)

2143

2144 cs%MassWghtInterp = 0

2145 if (usemasswghtinterp) &

2146 cs%MassWghtInterp = ibset(cs%MassWghtInterp, 0) ! Same as CS%MassWghtInterp + 1

2147 if (masswghtinterptop) &

2148 cs%MassWghtInterp = ibset(cs%MassWghtInterp, 1) ! Same as CS%MassWghtInterp + 2

2149 if ((.not.gv%Boussinesq) .and. masswghtinterp_nonbous_bug) &

2150 cs%MassWghtInterp = ibset(cs%MassWghtInterp, 3) ! Same as CS%MassWghtInterp + 8

2151

2152 call get_param(param_file, mdl, "CORRECTION_INTXPA", cs%correction_intxpa, &

2153 "If true, use a correction for surface pressure curvature in intx_pa.", &

2154 default=.false., do_not_log=.not.use_eos)

2155 call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL", cs%reset_intxpa_integral, &

2156 "If true, reset INTXPA to match pressures at first nonvanished cell. "//&

2157 "Includes pressure correction.", default=.false., do_not_log=.not.use_eos)

2158 call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL_FLATTEST", cs%reset_intxpa_flattest, &

2159 "If true, use flattest interface as reference interface where there is no "//&

2160 "better choice for RESET_INTXPA_INTEGRAL. Otherwise, use surface interface.", &

2161 default=.false., do_not_log=.not.use_eos)

2162 if (.not.use_eos) then ! These options do nothing without an equation of state.

2163 cs%correction_intxpa = .false.

2164 cs%reset_intxpa_integral = .false.

2165 cs%reset_intxpa_flattest = .false.

2166 endif

2167 call get_param(param_file, mdl, "RESET_INTXPA_H_NONVANISHED", cs%h_nonvanished, &

2168 "A minimal layer thickness that indicates that a layer is thick enough to usefully "//&

2169 "reestimate the pressure integral across the interface below.", &

2170 default=1.0e-6, units="m", scale=gv%m_to_H, do_not_log=.not.cs%reset_intxpa_integral)

2171 call get_param(param_file, mdl, "USE_INACCURATE_PGF_RHO_ANOM", cs%use_inaccurate_pgf_rho_anom, &

2172 "If true, use a form of the PGF that uses the reference density "//&

2173 "in an inaccurate way. This is not recommended.", default=.false.)

2174 call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &

2175 "If True, use vertical reconstruction of T & S within "//&

2176 "the integrals of the FV pressure gradient calculation. "//&

2177 "If False, use the constant-by-layer algorithm. "//&

2178 "The default is set by USE_REGRIDDING.", &

2179 default=use_ale )

2180 call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &

2181 "Order of vertical reconstruction of T/S to use in the "//&

2182 "integrals within the FV pressure gradient calculation.\n"//&

2183 " 0: PCM or no reconstruction.\n"//&

2184 " 1: PLM reconstruction.\n"//&

2185 " 2: PPM reconstruction.\n"//&

2186 " 3: PLM with least squares slope.", default=1)

2187 call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &

2188 "If true, the reconstruction of T & S for pressure in "//&

2189 "boundary cells is extrapolated, rather than using PCM "//&

2190 "in these cells. If true, the same order polynomial is "//&

2191 "used as is used for the interior cells.", default=.true.)

2192 call get_param(param_file, mdl, "USE_STANLEY_PGF", cs%use_stanley_pgf, &

2193 "If true, turn on Stanley SGS T variance parameterization "// &

2194 "in PGF code.", default=.false.)

2195 if (cs%use_stanley_pgf) then

2196 call get_param(param_file, mdl, "STANLEY_COEFF", stanley_coeff, &

2197 "Coefficient correlating the temperature gradient and SGS T variance.", &

2198 units="nondim", default=-1.0, do_not_log=.true.)

2199 if (stanley_coeff < 0.0) call mom_error(fatal, &

2200 "STANLEY_COEFF must be set >= 0 if USE_STANLEY_PGF is true.")

2201

2202 cs%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, &

2203 time, 'rho in PGF', 'kg m-3', conversion=us%R_to_kg_m3)

2204 cs%id_rho_stanley_pgf = register_diag_field('ocean_model', 'rho_stanley_pgf', diag%axesTL, &

2205 time, 'rho in PGF with Stanley correction', 'kg m-3', conversion=us%R_to_kg_m3)

2206 cs%id_p_stanley = register_diag_field('ocean_model', 'p_stanley', diag%axesTL, &

2207 time, 'p in PGF with Stanley correction', 'Pa', conversion=us%RL2_T2_to_Pa)

2208 endif

2209 if (cs%calculate_SAL) then

2210 cs%id_e_sal = register_diag_field('ocean_model', 'e_sal', diag%axesT1, time, &

2211 'Self-attraction and loading height anomaly', 'meter', conversion=us%Z_to_m)

2212 cs%id_sal_u = register_diag_field('ocean_model', 'SAL_u', diag%axesCuL, time, &

2213 'Zonal Acceleration due to self-attraction and loading', 'm s-2', conversion=us%L_T2_to_m_s2)

2214 cs%id_sal_v = register_diag_field('ocean_model', 'SAL_v', diag%axesCvL, time, &

2215 'Meridional Acceleration due to self-attraction and loading', 'm s-2', conversion=us%L_T2_to_m_s2)

2216 if (cs%id_sal_u > 0) &

2217 call safe_alloc_ptr(adp%sal_u, isdb, iedb, jsd, jed, nz)

2218 if (cs%id_sal_v > 0) &

2219 call safe_alloc_ptr(adp%sal_v, isd, ied, jsdb, jedb, nz)

2220 endif

2221 if (cs%tides) then

2222 cs%id_e_tide = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, time, &

2223 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)

2224 cs%id_e_tidal_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, time, &

2225 'Equilibrium tides height anomaly', 'meter', conversion=us%Z_to_m)

2226 cs%id_e_tidal_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, time, &

2227 'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=us%Z_to_m)

2228 cs%id_tides_u = register_diag_field('ocean_model', 'tides_u', diag%axesCuL, time, &

2229 'Zonal Acceleration due to tidal forcing', 'm s-2', conversion=us%L_T2_to_m_s2)

2230 cs%id_tides_v = register_diag_field('ocean_model', 'tides_v', diag%axesCvL, time, &

2231 'Meridional Acceleration due to tidal forcing', 'm s-2', conversion=us%L_T2_to_m_s2)

2232 if (cs%id_tides_u > 0) &

2233 call safe_alloc_ptr(adp%tides_u, isdb, iedb, jsd, jed, nz)

2234 if (cs%id_tides_v > 0) &

2235 call safe_alloc_ptr(adp%tides_v, isd, ied, jsdb, jedb, nz)

2236 endif

2237

2238 cs%id_MassWt_u = register_diag_field('ocean_model', 'MassWt_u', diag%axesCuL, time, &

2239 'The fractional mass weighting at u-point PGF calculations', 'nondim')

2240 cs%id_MassWt_v = register_diag_field('ocean_model', 'MassWt_v', diag%axesCvL, time, &

2241 'The fractional mass weighting at v-point PGF calculations', 'nondim')

2242

2243 cs%GFS_scale = 1.0

2244 if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth

2245

2246 call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale, units="nondim")

2247

2248end subroutine pressureforce_fv_init

2249

2250!> \namespace mom_pressureforce_fv

2251!!

2252!! Provides the Boussinesq and non-Boussinesq forms of horizontal accelerations

2253!! due to pressure gradients using a vertically integrated finite volume form,

2254!! as described by Adcroft et al., 2008. Integration in the vertical is made

2255!! either by quadrature or analytically.

2256!!

2257!! This form eliminates the thermobaric instabilities that had been a problem with

2258!! previous forms of the pressure gradient force calculation, as described by

2259!! Hallberg, 2005.

2260!!

2261!! Adcroft, A., R. Hallberg, and M. Harrison, 2008: A finite volume discretization

2262!! of the pressure gradient force using analytic integration. Ocean Modelling, 22,

2263!! 106-113. http://doi.org/10.1016/j.ocemod.2008.02.001

2264!!

2265!! Hallberg, 2005: A thermobaric instability of Lagrangian vertical coordinate

2266!! ocean models. Ocean Modelling, 8, 279-300.

2267!! http://dx.doi.org/10.1016/j.ocemod.2004.01.001

2268

2269end module mom_pressureforce_fv