MOM_density_integrals.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!> Provides integrals of density

6module mom_density_integrals

7

8use mom_eos, only : eos_type

9use mom_eos, only : eos_quadrature, eos_domain

10use mom_eos, only : analytic_int_density_dz

11use mom_eos, only : analytic_int_specific_vol_dp

12use mom_eos, only : calculate_density

13use mom_eos, only : calculate_spec_vol

14use mom_eos, only : calculate_specific_vol_derivs

15use mom_eos, only : average_specific_vol

16use mom_error_handler, only : mom_error, fatal, warning, mom_mesg

17use mom_hor_index, only : hor_index_type

18use mom_string_functions, only : uppercase

19use mom_variables, only : thermo_var_ptrs

20use mom_unit_scaling, only : unit_scale_type

21use mom_verticalgrid, only : verticalgrid_type

22

23implicit none ; private

24

25#include <MOM_memory.h>

26

27public int_density_dz

28public int_density_dz_generic_pcm

29public int_density_dz_generic_plm

30public int_density_dz_generic_ppm

31public int_specific_vol_dp

32public int_spec_vol_dp_generic_pcm

33public int_spec_vol_dp_generic_plm

34public avg_specific_vol

35public find_depth_of_pressure_in_cell

36public diagnose_mass_weight_z, diagnose_mass_weight_p

37

38contains

39

40!> Calls the appropriate subroutine to calculate analytical and nearly-analytical

41!! integrals in z across layers of pressure anomalies, which are

42!! required for calculating the finite-volume form pressure accelerations in a

43!! Boussinesq model.

44subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &

45 intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p, &

46 MassWghtInterpVanOnly, h_nv)

47 type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structures for the arrays

48 real, dimension(SZI_(HI),SZJ_(HI)), &

49 intent(in) :: t !< Potential temperature referenced to the surface [C ~> degC]

50 real, dimension(SZI_(HI),SZJ_(HI)), &

51 intent(in) :: s !< Salinity [S ~> ppt]

52 real, dimension(SZI_(HI),SZJ_(HI)), &

53 intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]

54 real, dimension(SZI_(HI),SZJ_(HI)), &

55 intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]

56 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is

57 !! subtracted out to reduce the magnitude of each of the

58 !! integrals.

59 real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used

60 !! to calculate the pressure (as p~=-z*rho_0*G_e)

61 !! used in the equation of state.

62 real, intent(in) :: g_e !< The Earth's gravitational acceleration

63 !! [L2 Z-1 T-2 ~> m s-2]

64 type(eos_type), intent(in) :: eos !< Equation of state structure

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

66 real, dimension(SZI_(HI),SZJ_(HI)), &

67 intent(inout) :: dpa !< The change in the pressure anomaly

68 !! across the layer [R L2 T-2 ~> Pa]

69 real, dimension(SZI_(HI),SZJ_(HI)), &

70 optional, intent(inout) :: intz_dpa !< The integral through the thickness of the

71 !! layer of the pressure anomaly relative to the

72 !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]

73 real, dimension(SZIB_(HI),SZJ_(HI)), &

74 optional, intent(inout) :: intx_dpa !< The integral in x of the difference between

75 !! the pressure anomaly at the top and bottom of the

76 !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]

77 real, dimension(SZI_(HI),SZJB_(HI)), &

78 optional, intent(inout) :: inty_dpa !< The integral in y of the difference between

79 !! the pressure anomaly at the top and bottom of the

80 !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]

81 real, dimension(SZI_(HI),SZJ_(HI)), &

82 optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]

83 real, dimension(SZI_(HI),SZJ_(HI)), &

84 optional, intent(in) :: ssh !< The sea surface height [Z ~> m]

85 real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]

86 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

87 !! mass weighting to interpolate T/S in integrals

88 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

89 !! of T/S unless one side smaller than h_nv (i.e. vanished)

90 real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]

91

92 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &

93 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]

94

95 if (eos_quadrature(eos)) then

96 call int_density_dz_generic_pcm(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, eos, us, dpa, &

97 intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, dz_neglect, &

98 masswghtinterp, z_0p=z_0p, masswghtinterpvanonly=masswghtinterpvanonly, &

99 h_nv=h_nv)

100 else

101 call analytic_int_density_dz(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, eos, dpa, &

102 intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, dz_neglect, masswghtinterp, z_0p=z_0p)

103 endif

104

105end subroutine int_density_dz

106

107

108!> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which

109!! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model.

110subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &

111 EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &

112 dz_neglect, MassWghtInterp, use_inaccurate_form, Z_0p, &

113 MassWghtInterpVanOnly, h_nv)

114 type(hor_index_type), intent(in) :: hi !< Horizontal index type for input variables.

115 real, dimension(SZI_(HI),SZJ_(HI)), &

116 intent(in) :: t !< Potential temperature of the layer [C ~> degC]

117 real, dimension(SZI_(HI),SZJ_(HI)), &

118 intent(in) :: s !< Salinity of the layer [S ~> ppt]

119 real, dimension(SZI_(HI),SZJ_(HI)), &

120 intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]

121 real, dimension(SZI_(HI),SZJ_(HI)), &

122 intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]

123 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is

124 !! subtracted out to reduce the magnitude

125 !! of each of the integrals.

126 real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used

127 !! to calculate the pressure (as p~=-z*rho_0*G_e)

128 !! used in the equation of state.

129 real, intent(in) :: g_e !< The Earth's gravitational acceleration

130 !! [L2 Z-1 T-2 ~> m s-2]

131 type(eos_type), intent(in) :: eos !< Equation of state structure

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

133 real, dimension(SZI_(HI),SZJ_(HI)), &

134 intent(inout) :: dpa !< The change in the pressure anomaly

135 !! across the layer [R L2 T-2 ~> Pa]

136 real, dimension(SZI_(HI),SZJ_(HI)), &

137 optional, intent(inout) :: intz_dpa !< The integral through the thickness of the

138 !! layer of the pressure anomaly relative to the

139 !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]

140 real, dimension(SZIB_(HI),SZJ_(HI)), &

141 optional, intent(inout) :: intx_dpa !< The integral in x of the difference between

142 !! the pressure anomaly at the top and bottom of the

143 !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]

144 real, dimension(SZI_(HI),SZJB_(HI)), &

145 optional, intent(inout) :: inty_dpa !< The integral in y of the difference between

146 !! the pressure anomaly at the top and bottom of the

147 !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]

148 real, dimension(SZI_(HI),SZJ_(HI)), &

149 optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]

150 real, dimension(SZI_(HI),SZJ_(HI)), &

151 optional, intent(in) :: ssh !< The sea surface height [Z ~> m]

152 real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]

153 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

154 !! mass weighting to interpolate T/S in integrals

155 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

156 !! of T/S unless one side smaller than h_nv (i.e. vanished)

157 real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]

158 logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of

159 !! density anomalies, as was used prior to March 2018.

160 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &

161 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]

162

163 ! Local variables

164 real :: t5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Temperatures along a line of subgrid locations [C ~> degC]

165 real :: s5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Salinities along a line of subgrid locations [S ~> ppt]

166 real :: p5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]

167 real :: r5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Densities anomalies along a line of subgrid locations [R ~> kg m-3]

168 real :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Temperatures at an array of subgrid locations [C ~> degC]

169 real :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Salinities at an array of subgrid locations [S ~> ppt]

170 real :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]

171 real :: r15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Densities at an array of subgrid locations [R ~> kg m-3]

172 real :: rho_anom ! The depth averaged density anomaly [R ~> kg m-3]

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

174 real :: gxrho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]

175 real :: dz ! The layer thickness [Z ~> m]

176 real :: dz_x(5,hi%iscb:hi%iecb) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]

177 real :: dz_y(5,hi%isc:hi%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]

178 real :: z0pres(hi%isd:hi%ied,hi%jsd:hi%jed) ! The height at which the pressure is zero [Z ~> m]

179 real :: hwght ! A pressure-thickness below topography [Z ~> m]

180 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m]

181 real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2]

182 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]

183 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]

184 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]

185 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]

186 real :: intz(5) ! The gravitational acceleration times the integrals of density

187 ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]

188 logical :: do_massweight ! Indicates whether to do mass weighting near bathymetry

189 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface

190 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

191 ! if at least one side vanished (0 or 1) [nondim]

192 real :: h_nonvanished ! nonvanished height [Z ~> m]

193 logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation

194 ! of density anomalies.

195 integer, dimension(2) :: eosdom_h5 ! The 5-point h-point i-computational domain for the equation of state

196 integer, dimension(2) :: eosdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state

197 integer, dimension(2) :: eosdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state

198 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, pos

199

200 ! These array bounds work for the indexing convention of the input arrays, but

201 ! on the computational domain defined for the output arrays.

202 isq = hi%IscB ; ieq = hi%IecB

203 jsq = hi%JscB ; jeq = hi%JecB

204 is = hi%isc ; ie = hi%iec

205 js = hi%jsc ; je = hi%jec

206

207 gxrho = g_e * rho_0

208 if (present(z_0p)) then

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

210 z0pres(i,j) = z_0p(i,j)

211 enddo ; enddo

212 else

213 z0pres(:,:) = 0.0

214 endif

215 use_rho_ref = .true.

216 if (present(use_inaccurate_form)) then

217 if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form

218 endif

219

220 do_massweight = .false. ; top_massweight = .false.

221 if (present(masswghtinterp)) then

222 do_massweight = btest(masswghtinterp, 0) ! True for odd values

223 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set

224 if (do_massweight .and. .not.present(bathyt)) call mom_error(fatal, &

225 "int_density_dz_generic: bathyT must be present if near-bottom mass weighting is in use.")

226 if (top_massweight .and. .not.present(ssh)) call mom_error(fatal, &

227 "int_density_dz_generic: SSH must be present if near-surface mass weighting is in use.")

228 if ((do_massweight .or. top_massweight) .and. .not.present(dz_neglect)) call mom_error(fatal, &

229 "int_density_dz_generic: dz_neglect must be present if mass weighting is in use.")

230 endif

231 massweightnvonlytoggle = 1.

232 if (present(masswghtinterpvanonly)) then

233 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

234 endif

235 h_nonvanished = 0.

236 if (present(h_nv)) then

237 h_nonvanished = h_nv

238 endif

239

240 ! Set the loop ranges for equation of state calculations at various points.

241 eosdom_h5(1) = 1 ; eosdom_h5(2) = 5*(ieq-isq+2)

242 eosdom_q15(1) = 1 ; eosdom_q15(2) = 15*(ieq-isq+1)

243 eosdom_h15(1) = 1 ; eosdom_h15(2) = 15*(hi%iec-hi%isc+1)

244

245 do j=jsq,jeq+1

246 do i=isq,ieq+1

247 dz = z_t(i,j) - z_b(i,j)

248 do n=1,5

249 t5(i*5+n) = t(i,j) ; s5(i*5+n) = s(i,j)

250 p5(i*5+n) = -gxrho*((z_t(i,j) - z0pres(i,j)) - 0.25*real(n-1)*dz)

251 enddo

252 enddo

253

254 if (use_rho_ref) then

255 call calculate_density(t5, s5, p5, r5, eos, eosdom_h5, rho_ref=rho_ref)

256 else

257 call calculate_density(t5, s5, p5, r5, eos, eosdom_h5)

258 endif

259

260 do i=isq,ieq+1

261 ! Use Boole's rule to estimate the pressure anomaly change.

262 rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))

263 if (.not.use_rho_ref) rho_anom = rho_anom - rho_ref

264 dz = z_t(i,j) - z_b(i,j)

265 dpa(i,j) = g_e*dz*rho_anom

266 ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of

267 ! the pressure anomaly.

268 if (present(intz_dpa)) intz_dpa(i,j) = 0.5*g_e*dz**2 * &

269 (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )

270 enddo

271 enddo

272

273 if (present(intx_dpa)) then ; do j=js,je

274 do i=isq,ieq

275 ! hWght is the distance measure by which the cell is violation of

276 ! hydrostatic consistency. For large hWght we bias the interpolation of

277 ! T & S along the top and bottom integrals, akin to thickness weighting.

278 hwght = 0.0

279 if (do_massweight) &

280 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))

281 if (top_massweight) &

282 hwght = max(hwght, z_b(i+1,j)-ssh(i,j), z_b(i,j)-ssh(i+1,j))

283 ! If both sides are nonvanished, then set it back to zero.

284 if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i+1,j) - z_b(i+1,j)) > h_nonvanished)) then

285 hwght = massweightnvonlytoggle * hwght

286 endif

287 if (hwght > 0.) then

288 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect

289 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect

290 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

291 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

292 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom

293 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom

294 else

295 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0

296 endif

297

298 do m=2,4

299 ! T, S, and z are interpolated in the horizontal. The z interpolation

300 ! is linear, but for T and S it may be thickness weighted.

301 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l

302 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)

303 dz_x(m,i) = (wt_l*(z_t(i,j) - z_b(i,j))) + (wt_r*(z_t(i+1,j) - z_b(i+1,j)))

304 pos = i*15+(m-2)*5

305 t15(pos+1) = (wtt_l*t(i,j)) + (wtt_r*t(i+1,j))

306 s15(pos+1) = (wtt_l*s(i,j)) + (wtt_r*s(i+1,j))

307 p15(pos+1) = -gxrho * ((wt_l*(z_t(i,j)-z0pres(i,j))) + (wt_r*(z_t(i+1,j)-z0pres(i+1,j))))

308 do n=2,5

309 t15(pos+n) = t15(pos+1) ; s15(pos+n) = s15(pos+1)

310 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)

311 enddo

312 enddo

313 enddo

314

315 if (use_rho_ref) then

316 call calculate_density(t15, s15, p15, r15, eos, eosdom_q15, rho_ref=rho_ref)

317 else

318 call calculate_density(t15, s15, p15, r15, eos, eosdom_q15)

319 endif

320

321 do i=isq,ieq

322 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)

323 ! Use Boole's rule to estimate the pressure anomaly change.

324 if (use_rho_ref) then

325 do m=2,4

326 pos = i*15+(m-2)*5

327 intz(m) = (g_e*dz_x(m,i)*(c1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + &

328 32.0*(r15(pos+2)+r15(pos+4)) + &

329 12.0*r15(pos+3)) ))

330 enddo

331 else

332 do m=2,4

333 pos = i*15+(m-2)*5

334 intz(m) = (g_e*dz_x(m,i)*(c1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + &

335 32.0*(r15(pos+2)+r15(pos+4)) + &

336 12.0*r15(pos+3)) - rho_ref ))

337 enddo

338 endif

339 ! Use Boole's rule to integrate the bottom pressure anomaly values in x.

340 intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &

341 12.0*intz(3))

342 enddo

343 enddo ; endif

344

345 if (present(inty_dpa)) then ; do j=jsq,jeq

346 do i=is,ie

347 ! hWght is the distance measure by which the cell is violation of

348 ! hydrostatic consistency. For large hWght we bias the interpolation of

349 ! T & S along the top and bottom integrals, akin to thickness weighting.

350 hwght = 0.0

351 if (do_massweight) &

352 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))

353 if (top_massweight) &

354 hwght = max(hwght, z_b(i,j+1)-ssh(i,j), z_b(i,j)-ssh(i,j+1))

355 ! If both sides are nonvanished, then set it back to zero.

356 if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i,j+1) - z_b(i,j+1)) > h_nonvanished)) then

357 hwght = massweightnvonlytoggle * hwght

358 endif

359 if (hwght > 0.) then

360 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect

361 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect

362 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

363 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

364 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom

365 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom

366 else

367 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0

368 endif

369

370 do m=2,4

371 ! T, S, and z are interpolated in the horizontal. The z interpolation

372 ! is linear, but for T and S it may be thickness weighted.

373 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l

374 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)

375 dz_y(m,i) = (wt_l*(z_t(i,j) - z_b(i,j))) + (wt_r*(z_t(i,j+1) - z_b(i,j+1)))

376 pos = i*15+(m-2)*5

377 t15(pos+1) = (wtt_l*t(i,j)) + (wtt_r*t(i,j+1))

378 s15(pos+1) = (wtt_l*s(i,j)) + (wtt_r*s(i,j+1))

379 p15(pos+1) = -gxrho * ((wt_l*(z_t(i,j)-z0pres(i,j))) + (wt_r*(z_t(i,j+1)-z0pres(i,j+1))))

380 do n=2,5

381 t15(pos+n) = t15(pos+1) ; s15(pos+n) = s15(pos+1)

382 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i)

383 enddo

384 enddo

385 enddo

386

387 if (use_rho_ref) then

388 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

389 r15(15*hi%isc+1:), eos, eosdom_h15, rho_ref=rho_ref)

390 else

391 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

392 r15(15*hi%isc+1:), eos, eosdom_h15)

393 endif

394

395 do i=is,ie

396 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)

397 ! Use Boole's rule to estimate the pressure anomaly change.

398 do m=2,4

399 pos = i*15+(m-2)*5

400 if (use_rho_ref) then

401 intz(m) = (g_e*dz_y(m,i)*(c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &

402 32.0*(r15(pos+2)+r15(pos+4)) + &

403 12.0*r15(pos+3)) ))

404 else

405 intz(m) = (g_e*dz_y(m,i)*(c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &

406 32.0*(r15(pos+2)+r15(pos+4)) + &

407 12.0*r15(pos+3)) - rho_ref ))

408 endif

409 enddo

410 ! Use Boole's rule to integrate the values.

411 inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &

412 12.0*intz(3))

413 enddo

414 enddo ; endif

415end subroutine int_density_dz_generic_pcm

416

417

418!> Compute pressure gradient force integrals by quadrature for the case where

419!! T and S are linear profiles.

420subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &

421 rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, dpa, &

422 intz_dpa, intx_dpa, inty_dpa, MassWghtInterp, &

423 use_inaccurate_form, Z_0p, MassWghtInterpVanOnly, h_nv)

424 integer, intent(in) :: k !< Layer index to calculate integrals for

425 type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structures for the input arrays

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

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

428 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

429 intent(in) :: t_t !< Potential temperature at the cell top [C ~> degC]

430 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

431 intent(in) :: t_b !< Potential temperature at the cell bottom [C ~> degC]

432 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

433 intent(in) :: s_t !< Salinity at the cell top [S ~> ppt]

434 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

435 intent(in) :: s_b !< Salinity at the cell bottom [S ~> ppt]

436 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &

437 intent(in) :: e !< Height of interfaces [Z ~> m]

438 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted

439 !! out to reduce the magnitude of each of the integrals.

440 real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used to calculate

441 !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state.

442 real, intent(in) :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]

443 real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m]

444 real, dimension(SZI_(HI),SZJ_(HI)), &

445 intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]

446 type(eos_type), intent(in) :: eos !< Equation of state structure

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

448 logical, intent(in) :: use_stanley_eos !< If true, turn on Stanley SGS T variance parameterization

449 real, dimension(SZI_(HI),SZJ_(HI)), &

450 intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]

451 real, dimension(SZI_(HI),SZJ_(HI)), &

452 optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of

453 !! the pressure anomaly relative to the anomaly at the

454 !! top of the layer [R L2 Z T-2 ~> Pa m]

455 real, dimension(SZIB_(HI),SZJ_(HI)), &

456 optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the

457 !! pressure anomaly at the top and bottom of the layer

458 !! divided by the x grid spacing [R L2 T-2 ~> Pa]

459 real, dimension(SZI_(HI),SZJB_(HI)), &

460 optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the

461 !! pressure anomaly at the top and bottom of the layer

462 !! divided by the y grid spacing [R L2 T-2 ~> Pa]

463 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

464 !! mass weighting to interpolate T/S in integrals

465 logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of

466 !! density anomalies, as was used prior to March 2018.

467 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

468 !! of T/S unless one side smaller than h_nv (i.e. vanished)

469 real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]

470 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &

471 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]

472

473! This subroutine calculates (by numerical quadrature) integrals of

474! pressure anomalies across layers, which are required for calculating the

475! finite-volume form pressure accelerations in a Boussinesq model. The one

476! potentially dodgy assumption here is that rho_0 is used both in the denominator

477! of the accelerations, and in the pressure used to calculated density (the

478! latter being -z*rho_0*G_e). These two uses could be separated if need be.

479!

480! It is assumed that the salinity and temperature profiles are linear in the

481! vertical. The top and bottom values within each layer are provided and

482! a linear interpolation is used to compute intermediate values.

483

484 ! Local variables

485 real :: t5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Temperatures along a line of subgrid locations [C ~> degC]

486 real :: s5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Salinities along a line of subgrid locations [S ~> ppt]

487 real :: t25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temperature variance along a line of subgrid

488 ! locations [C2 ~> degC2]

489 real :: ts5((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temp-salt covariance along a line of subgrid

490 ! locations [C S ~> degC ppt]

491 real :: s25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS salinity variance along a line of subgrid locations [S2 ~> ppt2]

492 real :: p5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]

493 real :: r5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Densities anomalies along a line of subgrid

494 ! locations [R ~> kg m-3]

495 real :: u5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Densities anomalies along a line of subgrid locations

496 ! (used for inaccurate form) [R ~> kg m-3]

497 real :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Temperatures at an array of subgrid locations [C ~> degC]

498 real :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Salinities at an array of subgrid locations [S ~> ppt]

499 real :: t215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temperature variance along a line of subgrid

500 ! locations [C2 ~> degC2]

501 real :: ts15((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temp-salt covariance along a line of subgrid

502 ! locations [C S ~> degC ppt]

503 real :: s215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS salinity variance along a line of subgrid

504 ! locations [S2 ~> ppt2]

505 real :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]

506 real :: r15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Densities at an array of subgrid locations [R ~> kg m-3]

507 real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]

508 real :: rho_anom ! A density anomaly [R ~> kg m-3]

509 real :: w_left, w_right ! Left and right weights [nondim]

510 real :: intz(5) ! The gravitational acceleration times the integrals of density

511 ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]

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

513 real :: gxrho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]

514 real :: dz(hi%iscb:hi%iecb+1) ! Layer thicknesses at tracer points [Z ~> m]

515 real :: dz_x(5,hi%iscb:hi%iecb) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]

516 real :: dz_y(5,hi%isc:hi%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]

517 real :: massweighttoggle ! A non-dimensional toggle factor for near-bottom mass weighting (0 or 1) [nondim]

518 real :: topweighttoggle ! A non-dimensional toggle factor for near-surface mass weighting (0 or 1) [nondim]

519 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

520 ! if at least one side vanished (0 or 1) [nondim]

521 real :: ttl, tbl, ttr, tbr ! Temperatures at the velocity cell corners [C ~> degC]

522 real :: stl, sbl, str, sbr ! Salinities at the velocity cell corners [S ~> ppt]

523 real :: z0pres(hi%isd:hi%ied,hi%jsd:hi%jed) ! The height at which the pressure is zero [Z ~> m]

524 real :: hwght ! A topographically limited thickness weight [Z ~> m]

525 real :: hwghttop ! An ice draft limited thickness weight [Z ~> m]

526 real :: hl, hr ! Thicknesses to the left and right [Z ~> m]

527 real :: idenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]

528 real :: h_nonvanished ! nonvanished height [Z ~> m]

529 logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation

530 ! of density anomalies.

531 logical :: use_vart, use_vars, use_covarts ! Logicals for SGS variances fields

532 integer, dimension(2) :: eosdom_h5 ! The 5-point h-point i-computational domain for the equation of state

533 integer, dimension(2) :: eosdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state

534 integer, dimension(2) :: eosdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state

535 integer :: isq, ieq, jsq, jeq, i, j, m, n, pos

536

537 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

538

539 gxrho = g_e * rho_0

540 if (present(z_0p)) then

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

542 z0pres(i,j) = z_0p(i,j)

543 enddo ; enddo

544 else

545 z0pres(:,:) = 0.0

546 endif

547 massweighttoggle = 0. ; topweighttoggle = 0.

548 if (present(masswghtinterp)) then

549 if (btest(masswghtinterp, 0)) massweighttoggle = 1.

550 if (btest(masswghtinterp, 1)) topweighttoggle = 1.

551 endif

552 massweightnvonlytoggle = 1.

553 if (present(masswghtinterpvanonly)) then

554 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

555 endif

556 h_nonvanished = 0.

557 if (present(h_nv)) then

558 h_nonvanished = h_nv

559 endif

560 use_rho_ref = .true.

561 if (present(use_inaccurate_form)) use_rho_ref = .not. use_inaccurate_form

562

563 use_vart = .false. !ensure initialized

564 use_covarts = .false.

565 use_vars = .false.

566 if (use_stanley_eos) then

567 use_vart = associated(tv%varT)

568 use_covarts = associated(tv%covarTS)

569 use_vars = associated(tv%varS)

570 endif

571

572 t25(:) = 0.

573 ts5(:) = 0.

574 s25(:) = 0.

575 t215(:) = 0.

576 ts15(:) = 0.

577 s215(:) = 0.

578

579 do n = 1, 5

580 wt_t(n) = 0.25 * real(5-n)

581 wt_b(n) = 1.0 - wt_t(n)

582 enddo

583

584 ! Set the loop ranges for equation of state calculations at various points.

585 eosdom_h5(1) = 1 ; eosdom_h5(2) = 5*(ieq-isq+2)

586 eosdom_q15(1) = 1 ; eosdom_q15(2) = 15*(ieq-isq+1)

587 eosdom_h15(1) = 1 ; eosdom_h15(2) = 15*(hi%iec-hi%isc+1)

588

589 ! 1. Compute vertical integrals

590 do j=jsq,jeq+1

591 do i = isq,ieq+1

592 dz(i) = e(i,j,k) - e(i,j,k+1)

593 do n=1,5

594 p5(i*5+n) = -gxrho*((e(i,j,k) - z0pres(i,j)) - 0.25*real(n-1)*dz(i))

595 ! Salinity and temperature points are linearly interpolated

596 s5(i*5+n) = wt_t(n) * s_t(i,j,k) + wt_b(n) * s_b(i,j,k)

597 t5(i*5+n) = wt_t(n) * t_t(i,j,k) + wt_b(n) * t_b(i,j,k)

598 enddo

599 if (use_vart) t25(i*5+1:i*5+5) = tv%varT(i,j,k)

600 if (use_covarts) ts5(i*5+1:i*5+5) = tv%covarTS(i,j,k)

601 if (use_vars) s25(i*5+1:i*5+5) = tv%varS(i,j,k)

602 enddo

603 if (use_stanley_eos) then

604 call calculate_density(t5, s5, p5, t25, ts5, s25, r5, eos, eosdom_h5, rho_ref=rho_ref)

605 else

606 if (use_rho_ref) then

607 call calculate_density(t5, s5, p5, r5, eos, eosdom_h5, rho_ref=rho_ref)

608 else

609 call calculate_density(t5, s5, p5, r5, eos, eosdom_h5)

610 u5(:) = r5(:) - rho_ref

611 endif

612 endif

613

614 if (use_rho_ref) then

615 do i=isq,ieq+1

616 ! Use Boole's rule to estimate the pressure anomaly change.

617 rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))

618 dpa(i,j) = g_e*dz(i)*rho_anom

619 if (present(intz_dpa)) then

620 ! Use a Boole's-rule-like fifth-order accurate estimate of

621 ! the double integral of the pressure anomaly.

622 intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &

623 (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )

624 endif

625 enddo

626 else

627 do i=isq,ieq+1

628 ! Use Boole's rule to estimate the pressure anomaly change.

629 rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) &

630 - rho_ref

631 dpa(i,j) = g_e*dz(i)*rho_anom

632 if (present(intz_dpa)) then

633 ! Use a Boole's-rule-like fifth-order accurate estimate of

634 ! the double integral of the pressure anomaly.

635 intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &

636 (rho_anom - c1_90*(16.0*(u5(i*5+4)-u5(i*5+2)) + 7.0*(u5(i*5+5)-u5(i*5+1))) )

637 endif

638 enddo

639 endif

640 enddo ! end loops on j

641

642 ! 2. Compute horizontal integrals in the x direction

643 if (present(intx_dpa)) then ; do j=hi%jsc,hi%jec

644 do i=isq,ieq

645 ! Corner values of T and S

646 ! hWght is the distance measure by which the cell is violation of

647 ! hydrostatic consistency. For large hWght we bias the interpolation

648 ! of T,S along the top and bottom integrals, almost like thickness

649 ! weighting.

650 ! Note: To work in terrain following coordinates we could offset

651 ! this distance by the layer thickness to replicate other models.

652 hwght = massweighttoggle * &

653 max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k))

654 ! CY: The below code just uses top interface, which may be bad in high res open ocean

655 ! We want something like if (pa(i+1,k+1)<pa(i,1)) or (pa(i+1,1) <pa(i,k+1)) then...

656 ! but pressures are not passed through to this submodule, and tv just has surface press.

657 !if ((p(i+1,j,k+1)<p(i,j,1)).or.(tv%p(i+1,j,k+1)<tv%p(i,j,1))) then

658 hwghttop = topweighttoggle * &

659 max(0., e(i+1,j,k+1)-e(i,j,1), e(i,j,k+1)-e(i+1,j,1))

660 !else ! pressure criteria not activated

661 ! hWghtTop = 0.

662 !endif

663 ! Set it to be max of the bottom and top hWghts:

664 hwght = max(hwght, hwghttop)

665 ! If both sides are nonvanished, then set it back to zero.

666 if (((e(i,j,k) - e(i,j,k+1)) > h_nonvanished) .and. ((e(i+1,j,k) - e(i+1,j,k+1)) > h_nonvanished)) then

667 hwght = massweightnvonlytoggle * hwght

668 endif

669 if (hwght > 0.) then

670 hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff

671 hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz_subroundoff

672 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

673 idenom = 1./( hwght*(hr + hl) + hl*hr )

674 ttl = ( (hwght*hr)*t_t(i+1,j,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom

675 ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i+1,j,k) ) * idenom

676 tbl = ( (hwght*hr)*t_b(i+1,j,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom

677 tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i+1,j,k) ) * idenom

678 stl = ( (hwght*hr)*s_t(i+1,j,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom

679 str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i+1,j,k) ) * idenom

680 sbl = ( (hwght*hr)*s_b(i+1,j,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom

681 sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i+1,j,k) ) * idenom

682 else

683 ttl = t_t(i,j,k) ; tbl = t_b(i,j,k) ; ttr = t_t(i+1,j,k) ; tbr = t_b(i+1,j,k)

684 stl = s_t(i,j,k) ; sbl = s_b(i,j,k) ; str = s_t(i+1,j,k) ; sbr = s_b(i+1,j,k)

685 endif

686

687 do m=2,4

688 w_left = wt_t(m) ; w_right = wt_b(m)

689 dz_x(m,i) = (w_left*(e(i,j,k) - e(i,j,k+1))) + (w_right*(e(i+1,j,k) - e(i+1,j,k+1)))

690

691 ! Salinity and temperature points are linearly interpolated in

692 ! the horizontal. The subscript (1) refers to the top value in

693 ! the vertical profile while subscript (5) refers to the bottom

694 ! value in the vertical profile.

695 pos = i*15+(m-2)*5

696 t15(pos+1) = (w_left*ttl) + (w_right*ttr)

697 t15(pos+5) = (w_left*tbl) + (w_right*tbr)

698

699 s15(pos+1) = (w_left*stl) + (w_right*str)

700 s15(pos+5) = (w_left*sbl) + (w_right*sbr)

701

702 p15(pos+1) = -gxrho * ((w_left*(e(i,j,k)-z0pres(i,j))) + (w_right*(e(i+1,j,k)-z0pres(i+1,j))))

703

704 ! Pressure

705 do n=2,5

706 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)

707 enddo

708

709 ! Salinity and temperature (linear interpolation in the vertical)

710 do n=2,4

711 s15(pos+n) = wt_t(n) * s15(pos+1) + wt_b(n) * s15(pos+5)

712 t15(pos+n) = wt_t(n) * t15(pos+1) + wt_b(n) * t15(pos+5)

713 enddo

714 if (use_vart) t215(pos+1:pos+5) = (w_left*tv%varT(i,j,k)) + (w_right*tv%varT(i+1,j,k))

715 if (use_covarts) ts15(pos+1:pos+5) = (w_left*tv%covarTS(i,j,k)) + (w_right*tv%covarTS(i+1,j,k))

716 if (use_vars) s215(pos+1:pos+5) = (w_left*tv%varS(i,j,k)) + (w_right*tv%varS(i+1,j,k))

717 enddo

718 enddo

719

720 if (use_stanley_eos) then

721 call calculate_density(t15, s15, p15, t215, ts15, s215, r15, eos, eosdom_q15, rho_ref=rho_ref)

722 else

723 if (use_rho_ref) then

724 call calculate_density(t15, s15, p15, r15, eos, eosdom_q15, rho_ref=rho_ref)

725 else

726 call calculate_density(t15, s15, p15, r15, eos, eosdom_q15)

727 endif

728 endif

729

730 do i=isq,ieq

731 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)

732

733 ! Use Boole's rule to estimate the pressure anomaly change.

734 if (use_rho_ref) then

735 do m = 2,4

736 pos = i*15+(m-2)*5

737 intz(m) = (g_e*dz_x(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &

738 12.0*r15(pos+3)) ))

739 enddo

740 else

741 do m = 2,4

742 pos = i*15+(m-2)*5

743 intz(m) = (g_e*dz_x(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &

744 12.0*r15(pos+3)) - rho_ref ))

745 enddo

746 endif

747 ! Use Boole's rule to integrate the bottom pressure anomaly values in x.

748 intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &

749 12.0*intz(3))

750 enddo

751 enddo ; endif

752

753 ! 3. Compute horizontal integrals in the y direction

754 if (present(inty_dpa)) then ; do j=jsq,jeq

755 do i=hi%isc,hi%iec

756 ! Corner values of T and S

757 ! hWght is the distance measure by which the cell is violation of

758 ! hydrostatic consistency. For large hWght we bias the interpolation

759 ! of T,S along the top and bottom integrals, almost like thickness

760 ! weighting.

761 ! Note: To work in terrain following coordinates we could offset

762 ! this distance by the layer thickness to replicate other models.

763 hwght = massweighttoggle * &

764 max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k))

765 ! CY: The below code just uses top interface, which may be bad in high res open ocean

766 ! We want something like if (pa(j+1,k+1)<pa(j,1)) or (pa(j+1,1) <pa(i,j,k+1)) then...

767 ! but pressures are not passed through to this submodule, and tv just has surface press.

768 !if ((p(i,j+1,k+1)<p(i,j,1)).or.(tv%p(i,j+1,k+1)<tv%p(i,j,1))) then

769 hwghttop = topweighttoggle * &

770 max(0., e(i,j+1,k+1)-e(i,j,1), e(i,j,k+1)-e(i,j+1,1))

771 !else ! pressure criteria not activated

772 ! hWghtTop = 0.

773 !endif

774 ! Set it to be max of the bottom and top hWghts:

775 hwght = max(hwght, hwghttop)

776 ! If both sides are nonvanished, then set it back to zero.

777 if (((e(i,j,k) - e(i,j,k+1)) > h_nonvanished) .and. ((e(i,j+1,k) - e(i,j+1,k+1)) > h_nonvanished)) then

778 hwght = massweightnvonlytoggle * hwght

779 endif

780

781 if (hwght > 0.) then

782 hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff

783 hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz_subroundoff

784 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

785 idenom = 1./( hwght*(hr + hl) + hl*hr )

786 ttl = ( (hwght*hr)*t_t(i,j+1,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom

787 ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i,j+1,k) ) * idenom

788 tbl = ( (hwght*hr)*t_b(i,j+1,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom

789 tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i,j+1,k) ) * idenom

790 stl = ( (hwght*hr)*s_t(i,j+1,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom

791 str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i,j+1,k) ) * idenom

792 sbl = ( (hwght*hr)*s_b(i,j+1,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom

793 sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i,j+1,k) ) * idenom

794 else

795 ttl = t_t(i,j,k) ; tbl = t_b(i,j,k) ; ttr = t_t(i,j+1,k) ; tbr = t_b(i,j+1,k)

796 stl = s_t(i,j,k) ; sbl = s_b(i,j,k) ; str = s_t(i,j+1,k) ; sbr = s_b(i,j+1,k)

797 endif

798

799 do m=2,4

800 w_left = wt_t(m) ; w_right = wt_b(m)

801 dz_y(m,i) = (w_left*(e(i,j,k) - e(i,j,k+1))) + (w_right*(e(i,j+1,k) - e(i,j+1,k+1)))

802

803 ! Salinity and temperature points are linearly interpolated in

804 ! the horizontal. The subscript (1) refers to the top value in

805 ! the vertical profile while subscript (5) refers to the bottom

806 ! value in the vertical profile.

807 pos = i*15+(m-2)*5

808 t15(pos+1) = (w_left*ttl) + (w_right*ttr)

809 t15(pos+5) = (w_left*tbl) + (w_right*tbr)

810

811 s15(pos+1) = (w_left*stl) + (w_right*str)

812 s15(pos+5) = (w_left*sbl) + (w_right*sbr)

813

814 p15(pos+1) = -gxrho * ((w_left*(e(i,j,k)-z0pres(i,j))) + (w_right*(e(i,j+1,k)-z0pres(i,j+1))))

815

816 ! Pressure

817 do n=2,5

818 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i)

819 enddo

820

821 ! Salinity and temperature (linear interpolation in the vertical)

822 do n=2,4

823 s15(pos+n) = wt_t(n) * s15(pos+1) + wt_b(n) * s15(pos+5)

824 t15(pos+n) = wt_t(n) * t15(pos+1) + wt_b(n) * t15(pos+5)

825 enddo

826 if (use_vart) t215(pos+1:pos+5) = (w_left*tv%varT(i,j,k)) + (w_right*tv%varT(i,j+1,k))

827 if (use_covarts) ts15(pos+1:pos+5) = (w_left*tv%covarTS(i,j,k)) + (w_right*tv%covarTS(i,j+1,k))

828 if (use_vars) s215(pos+1:pos+5) = (w_left*tv%varS(i,j,k)) + (w_right*tv%varS(i,j+1,k))

829 enddo

830 enddo

831

832 if (use_stanley_eos) then

833 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

834 t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), &

835 r15(15*hi%isc+1:), eos, eosdom_h15, rho_ref=rho_ref)

836 else

837 if (use_rho_ref) then

838 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

839 r15(15*hi%isc+1:), eos, eosdom_h15, rho_ref=rho_ref)

840 else

841 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

842 r15(15*hi%isc+1:), eos, eosdom_h15)

843 endif

844 endif

845

846 do i=hi%isc,hi%iec

847 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)

848

849 ! Use Boole's rule to estimate the pressure anomaly change.

850 if (use_rho_ref) then

851 do m = 2,4

852 pos = i*15+(m-2)*5

853 intz(m) = (g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &

854 32.0*(r15(pos+2)+r15(pos+4)) + &

855 12.0*r15(pos+3)) ))

856 enddo

857 else

858 do m = 2,4

859 pos = i*15+(m-2)*5

860 intz(m) = (g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &

861 32.0*(r15(pos+2)+r15(pos+4)) + &

862 12.0*r15(pos+3)) - rho_ref ))

863 enddo

864 endif

865 ! Use Boole's rule to integrate the values.

866 inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &

867 12.0*intz(3))

868 enddo

869 enddo ; endif

870

871end subroutine int_density_dz_generic_plm

872

873

874!> Compute pressure gradient force integrals for layer "k" and the case where T and S

875!! are parabolic profiles

876subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &

877 rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, &

878 dpa, intz_dpa, intx_dpa, inty_dpa, MassWghtInterp, Z_0p, &

879 MassWghtInterpVanOnly, h_nv)

880 integer, intent(in) :: k !< Layer index to calculate integrals for

881 type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structures for the input arrays

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

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

884 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

885 intent(in) :: t_t !< Potential temperature at the cell top [C ~> degC]

886 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

887 intent(in) :: t_b !< Potential temperature at the cell bottom [C ~> degC]

888 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

889 intent(in) :: s_t !< Salinity at the cell top [S ~> ppt]

890 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &

891 intent(in) :: s_b !< Salinity at the cell bottom [S ~> ppt]

892 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &

893 intent(in) :: e !< Height of interfaces [Z ~> m]

894 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is

895 !! subtracted out to reduce the magnitude of each of the integrals.

896 real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used to calculate

897 !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state.

898 real, intent(in) :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]

899 real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m]

900 real, dimension(SZI_(HI),SZJ_(HI)), &

901 intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]

902 type(eos_type), intent(in) :: eos !< Equation of state structure

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

904 logical, intent(in) :: use_stanley_eos !< If true, turn on Stanley SGS T variance parameterization

905 real, dimension(SZI_(HI),SZJ_(HI)), &

906 intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]

907 real, dimension(SZI_(HI),SZJ_(HI)), &

908 optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of

909 !! the pressure anomaly relative to the anomaly at the

910 !! top of the layer [R L2 Z T-2 ~> Pa m]

911 real, dimension(SZIB_(HI),SZJ_(HI)), &

912 optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the

913 !! pressure anomaly at the top and bottom of the layer

914 !! divided by the x grid spacing [R L2 T-2 ~> Pa]

915 real, dimension(SZI_(HI),SZJB_(HI)), &

916 optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the

917 !! pressure anomaly at the top and bottom of the layer

918 !! divided by the y grid spacing [R L2 T-2 ~> Pa]

919 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

920 !! mass weighting to interpolate T/S in integrals

921 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

922 !! of T/S unless one side smaller than h_nv (i.e. vanished)

923 real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]

924

925 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &

926 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]

927

928! This subroutine calculates (by numerical quadrature) integrals of

929! pressure anomalies across layers, which are required for calculating the

930! finite-volume form pressure accelerations in a Boussinesq model. The one

931! potentially dodgy assumption here is that rho_0 is used both in the denominator

932! of the accelerations, and in the pressure used to calculated density (the

933! latter being -z*rho_0*G_e). These two uses could be separated if need be.

934!

935! It is assumed that the salinity and temperature profiles are parabolic in the

936! vertical. The top and bottom values within each layer are provided and

937! a parabolic interpolation is used to compute intermediate values.

938

939 ! Local variables

940 real :: t5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Temperatures along a line of subgrid locations [C ~> degC]

941 real :: s5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Salinities along a line of subgrid locations [S ~> ppt]

942 real :: t25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temperature variance along a line of subgrid

943 ! locations [C2 ~> degC2]

944 real :: ts5((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temp-salt covariance along a line of subgrid

945 ! locations [C S ~> degC ppt]

946 real :: s25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS salinity variance along a line of subgrid locations [S2 ~> ppt2]

947 real :: p5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]

948 real :: r5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Densities anomalies along a line of subgrid

949 ! locations [R ~> kg m-3]

950 real :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Temperatures at an array of subgrid locations [C ~> degC]

951 real :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Salinities at an array of subgrid locations [S ~> ppt]

952 real :: t215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temperature variance along a line of subgrid

953 ! locations [C2 ~> degC2]

954 real :: ts15((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temp-salt covariance along a line of subgrid

955 ! locations [C S ~> degC ppt]

956 real :: s215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS salinity variance along a line of subgrid

957 ! locations [S2 ~> ppt2]

958 real :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]

959 real :: r15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Densities at an array of subgrid locations [R ~> kg m-3]

960 real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]

961 real :: rho_anom ! The integrated density anomaly [R ~> kg m-3]

962 real :: w_left, w_right ! Left and right weights [nondim]

963 real :: intz(5) ! The gravitational acceleration times the integrals of density

964 ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]

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

966 real :: gxrho ! The gravitational acceleration times density [R L2 Z-1 T-2 ~> kg m-2 s-2]

967 real :: dz ! Layer thicknesses at tracer points [Z ~> m]

968 real :: dz_x(5,hi%iscb:hi%iecb) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]

969 real :: dz_y(5,hi%isc:hi%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]

970 real :: massweighttoggle ! A non-dimensional toggle factor for near-bottom mass weighting (0 or 1) [nondim]

971 real :: topweighttoggle ! A non-dimensional toggle factor for near-surface mass weighting (0 or 1) [nondim]

972 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

973 ! if at least one side vanished (0 or 1) [nondim]

974 real :: ttl, tbl, tml, ttr, tbr, tmr ! Temperatures at the velocity cell corners [C ~> degC]

975 real :: stl, sbl, sml, str, sbr, smr ! Salinities at the velocity cell corners [S ~> ppt]

976 real :: s6 ! PPM curvature coefficient for S [S ~> ppt]

977 real :: t6 ! PPM curvature coefficient for T [C ~> degC]

978 real :: t_top, t_mn, t_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T [C ~> degC]

979 real :: s_top, s_mn, s_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S [S ~> ppt]

980 real :: z0pres(hi%isd:hi%ied,hi%jsd:hi%jed) ! The height at which the pressure is zero [Z ~> m]

981 real :: hwght ! A topographically limited thickness weight [Z ~> m]

982 real :: hwghttop ! A surface displacement limited thickness weight [Z ~> m]

983 real :: hl, hr ! Thicknesses to the left and right [Z ~> m]

984 real :: idenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]

985 real :: h_nonvanished ! nonvanished height [Z ~> m]

986 integer, dimension(2) :: eosdom_h5 ! The 5-point h-point i-computational domain for the equation of state

987 integer, dimension(2) :: eosdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state

988 integer, dimension(2) :: eosdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state

989 integer :: isq, ieq, jsq, jeq, i, j, m, n, pos

990 logical :: use_ppm ! If false, assume zero curvature in reconstruction, i.e. PLM

991 logical :: use_vart, use_vars, use_covarts ! Logicals for SGS variances fields

992

993 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

994

995 gxrho = g_e * rho_0

996 if (present(z_0p)) then

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

998 z0pres(i,j) = z_0p(i,j)

999 enddo ; enddo

1000 else

1001 z0pres(:,:) = 0.0

1002 endif

1003 massweighttoggle = 0. ; topweighttoggle = 0.

1004 if (present(masswghtinterp)) then

1005 if (btest(masswghtinterp, 0)) massweighttoggle = 1.

1006 if (btest(masswghtinterp, 1)) topweighttoggle = 1.

1007 endif

1008 massweightnvonlytoggle = 1.

1009 if (present(masswghtinterpvanonly)) then

1010 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

1011 endif

1012 h_nonvanished = 0.

1013 if (present(h_nv)) then

1014 h_nonvanished = h_nv

1015 endif

1016

1017 ! In event PPM calculation is bypassed with use_PPM=False

1018 s6 = 0.

1019 t6 = 0.

1020 use_ppm = .true. ! This is a place-holder to allow later re-use of this function

1021

1022 use_vart = .false. !ensure initialized

1023 use_covarts = .false.

1024 use_vars = .false.

1025 if (use_stanley_eos) then

1026 use_vart = associated(tv%varT)

1027 use_covarts = associated(tv%covarTS)

1028 use_vars = associated(tv%varS)

1029 endif

1030

1031 t25(:) = 0.

1032 ts5(:) = 0.

1033 s25(:) = 0.

1034 t215(:) = 0.

1035 ts15(:) = 0.

1036 s215(:) = 0.

1037

1038 do n = 1, 5

1039 wt_t(n) = 0.25 * real(5-n)

1040 wt_b(n) = 1.0 - wt_t(n)

1041 enddo

1042

1043 ! Set the loop ranges for equation of state calculations at various points.

1044 eosdom_h5(1) = 1 ; eosdom_h5(2) = 5*(ieq-isq+2)

1045 eosdom_q15(1) = 1 ; eosdom_q15(2) = 15*(ieq-isq+1)

1046 eosdom_h15(1) = 1 ; eosdom_h15(2) = 15*(hi%iec-hi%isc+1)

1047

1048 ! 1. Compute vertical integrals

1049 do j=jsq,jeq+1

1050 do i=isq,ieq+1

1051 if (use_ppm) then

1052 ! Curvature coefficient of the parabolas

1053 s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( s_t(i,j,k) + s_b(i,j,k) ) )

1054 t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( t_t(i,j,k) + t_b(i,j,k) ) )

1055 endif

1056 dz = e(i,j,k) - e(i,j,k+1)

1057 do n=1,5

1058 p5(i*5+n) = -gxrho*((e(i,j,k) - z0pres(i,j)) - 0.25*real(n-1)*dz)

1059 ! Salinity and temperature points are reconstructed with PPM

1060 s5(i*5+n) = wt_t(n) * s_t(i,j,k) + wt_b(n) * ( s_b(i,j,k) + s6 * wt_t(n) )

1061 t5(i*5+n) = wt_t(n) * t_t(i,j,k) + wt_b(n) * ( t_b(i,j,k) + t6 * wt_t(n) )

1062 enddo

1063 if (use_stanley_eos) then

1064 if (use_vart) t25(i*5+1:i*5+5) = tv%varT(i,j,k)

1065 if (use_covarts) ts5(i*5+1:i*5+5) = tv%covarTS(i,j,k)

1066 if (use_vars) s25(i*5+1:i*5+5) = tv%varS(i,j,k)

1067 endif

1068 enddo

1069

1070 if (use_stanley_eos) then

1071 call calculate_density(t5, s5, p5, t25, ts5, s25, r5, eos, eosdom_h5, rho_ref=rho_ref)

1072 else

1073 call calculate_density(t5, s5, p5, r5, eos, eosdom_h5, rho_ref=rho_ref)

1074 endif

1075

1076 do i=isq,ieq+1

1077 dz = e(i,j,k) - e(i,j,k+1)

1078 ! Use Boole's rule to estimate the pressure anomaly change.

1079 rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))

1080 dpa(i,j) = g_e*dz*rho_anom

1081 if (present(intz_dpa)) then

1082 ! Use a Boole's-rule-like fifth-order accurate estimate of

1083 ! the double integral of the pressure anomaly.

1084 intz_dpa(i,j) = 0.5*g_e*dz**2 * &

1085 (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )

1086 endif

1087 enddo ! end loop on i

1088 enddo ! end loop on j

1089

1090 ! 2. Compute horizontal integrals in the x direction

1091 if (present(intx_dpa)) then ; do j=hi%jsc,hi%jec

1092 do i=isq,ieq

1093 ! Corner values of T and S

1094 ! hWght is the distance measure by which the cell is violation of

1095 ! hydrostatic consistency. For large hWght we bias the interpolation

1096 ! of T,S along the top and bottom integrals, almost like thickness

1097 ! weighting.

1098 ! Note: To work in terrain following coordinates we could offset

1099 ! this distance by the layer thickness to replicate other models.

1100 hwght = massweighttoggle * &

1101 max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k))

1102 hwghttop = topweighttoggle * &

1103 max(0., e(i+1,j,k+1)-e(i,j,1), e(i,j,k+1)-e(i+1,j,1))

1104 hwght = max(hwght, hwghttop)

1105 ! If both sides are nonvanished, then set it back to zero.

1106 if (((e(i,j,k) - e(i,j,k+1)) > h_nonvanished) .and. ((e(i+1,j,k) - e(i+1,j,k+1)) > h_nonvanished)) then

1107 hwght = massweightnvonlytoggle * hwght

1108 endif

1109 if (hwght > 0.) then

1110 hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff

1111 hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz_subroundoff

1112 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

1113 idenom = 1./( hwght*(hr + hl) + hl*hr )

1114 ttl = ( (hwght*hr)*t_t(i+1,j,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom

1115 tbl = ( (hwght*hr)*t_b(i+1,j,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom

1116 tml = ( (hwght*hr)*tv%T(i+1,j,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom

1117 ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i+1,j,k) ) * idenom

1118 tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i+1,j,k) ) * idenom

1119 tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i+1,j,k) ) * idenom

1120 stl = ( (hwght*hr)*s_t(i+1,j,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom

1121 sbl = ( (hwght*hr)*s_b(i+1,j,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom

1122 sml = ( (hwght*hr)*tv%S(i+1,j,k) + (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom

1123 str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i+1,j,k) ) * idenom

1124 sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i+1,j,k) ) * idenom

1125 smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i+1,j,k) ) * idenom

1126 else

1127 ttl = t_t(i,j,k) ; tbl = t_b(i,j,k) ; ttr = t_t(i+1,j,k) ; tbr = t_b(i+1,j,k)

1128 tml = tv%T(i,j,k) ; tmr = tv%T(i+1,j,k)

1129 stl = s_t(i,j,k) ; sbl = s_b(i,j,k) ; str = s_t(i+1,j,k) ; sbr = s_b(i+1,j,k)

1130 sml = tv%S(i,j,k) ; smr = tv%S(i+1,j,k)

1131 endif

1132

1133 do m=2,4

1134 w_left = wt_t(m) ; w_right = wt_b(m)

1135

1136 ! Salinity and temperature points are linearly interpolated in

1137 ! the horizontal. The subscript (1) refers to the top value in

1138 ! the vertical profile while subscript (5) refers to the bottom

1139 ! value in the vertical profile.

1140 t_top = (w_left*ttl) + (w_right*ttr)

1141 t_mn = (w_left*tml) + (w_right*tmr)

1142 t_bot = (w_left*tbl) + (w_right*tbr)

1143

1144 s_top = (w_left*stl) + (w_right*str)

1145 s_mn = (w_left*sml) + (w_right*smr)

1146 s_bot = (w_left*sbl) + (w_right*sbr)

1147

1148 ! Pressure

1149 dz_x(m,i) = (w_left*(e(i,j,k) - e(i,j,k+1))) + (w_right*(e(i+1,j,k) - e(i+1,j,k+1)))

1150

1151 pos = i*15+(m-2)*5

1152 p15(pos+1) = -gxrho * ((w_left*(e(i,j,k)-z0pres(i,j))) + (w_right*(e(i+1,j,k)-z0pres(i+1,j))))

1153 do n=2,5

1154 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)

1155 enddo

1156

1157 ! Parabolic reconstructions in the vertical for T and S

1158 if (use_ppm) then

1159 ! Coefficients of the parabolas

1160 s6 = 3.0 * ( 2.0*s_mn - ( s_top + s_bot ) )

1161 t6 = 3.0 * ( 2.0*t_mn - ( t_top + t_bot ) )

1162 endif

1163 do n=1,5

1164 s15(pos+n) = wt_t(n) * s_top + wt_b(n) * ( s_bot + s6 * wt_t(n) )

1165 t15(pos+n) = wt_t(n) * t_top + wt_b(n) * ( t_bot + t6 * wt_t(n) )

1166 enddo

1167 if (use_stanley_eos) then

1168 if (use_vart) t215(pos+1:pos+5) = (w_left*tv%varT(i,j,k)) + (w_right*tv%varT(i+1,j,k))

1169 if (use_covarts) ts15(pos+1:pos+5) = (w_left*tv%covarTS(i,j,k)) + (w_right*tv%covarTS(i+1,j,k))

1170 if (use_vars) s215(pos+1:pos+5) = (w_left*tv%varS(i,j,k)) + (w_right*tv%varS(i+1,j,k))

1171 endif

1172 if (use_stanley_eos) then

1173 call calculate_density(t5, s5, p5, t25, ts5, s25, r5, eos, rho_ref=rho_ref)

1174 else

1175 call calculate_density(t5, s5, p5, r5, eos, rho_ref=rho_ref)

1176 endif

1177 enddo

1178 enddo

1179

1180 if (use_stanley_eos) then

1181 call calculate_density(t15, s15, p15, t215, ts15, s215, r15, eos, eosdom_q15, rho_ref=rho_ref)

1182 else

1183 call calculate_density(t15, s15, p15, r15, eos, eosdom_q15, rho_ref=rho_ref)

1184 endif

1185

1186 do i=isq,ieq

1187 do m=2,4

1188 pos = i*15+(m-2)*5

1189 ! Use Boole's rule to estimate the pressure anomaly change.

1190 intz(m) = (g_e*dz_x(m,i)*(c1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + &

1191 32.0*(r15(pos+2)+r15(pos+4)) + &

1192 12.0*r15(pos+3)) ))

1193 enddo ! m

1194 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)

1195

1196 ! Use Boole's rule to integrate the bottom pressure anomaly values in x.

1197 intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))

1198

1199 enddo

1200 enddo ; endif

1201

1202 ! 3. Compute horizontal integrals in the y direction

1203 if (present(inty_dpa)) then ; do j=jsq,jeq

1204 do i=hi%isc,hi%iec

1205 ! Corner values of T and S

1206 ! hWght is the distance measure by which the cell is violation of

1207 ! hydrostatic consistency. For large hWght we bias the interpolation

1208 ! of T,S along the top and bottom integrals, almost like thickness

1209 ! weighting.

1210 ! Note: To work in terrain following coordinates we could offset

1211 ! this distance by the layer thickness to replicate other models.

1212 hwght = massweighttoggle * &

1213 max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k))

1214 hwghttop = topweighttoggle * &

1215 max(0., e(i,j+1,k+1)-e(i,j,1), e(i,j,k+1)-e(i,j+1,1))

1216 hwght = max(hwght, hwghttop)

1217 ! If both sides are nonvanished, then set it back to zero.

1218 if (((e(i,j,k) - e(i,j,k+1)) > h_nonvanished) .and. ((e(i,j+1,k) - e(i,j+1,k+1)) > h_nonvanished)) then

1219 hwght = massweightnvonlytoggle * hwght

1220 endif

1221 if (hwght > 0.) then

1222 hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff

1223 hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz_subroundoff

1224 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

1225 idenom = 1./( hwght*(hr + hl) + hl*hr )

1226 ttl = ( (hwght*hr)*t_t(i,j+1,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom

1227 tbl = ( (hwght*hr)*t_b(i,j+1,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom

1228 tml = ( (hwght*hr)*tv%T(i,j+1,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom

1229 ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i,j+1,k) ) * idenom

1230 tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i,j+1,k) ) * idenom

1231 tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i,j+1,k) ) * idenom

1232 stl = ( (hwght*hr)*s_t(i,j+1,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom

1233 sbl = ( (hwght*hr)*s_b(i,j+1,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom

1234 sml = ( (hwght*hr)*tv%S(i,j+1,k)+ (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom

1235 str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i,j+1,k) ) * idenom

1236 sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i,j+1,k) ) * idenom

1237 smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i,j+1,k) ) * idenom

1238 else

1239 ttl = t_t(i,j,k) ; tbl = t_b(i,j,k) ; ttr = t_t(i,j+1,k) ; tbr = t_b(i,j+1,k)

1240 tml = tv%T(i,j,k) ; tmr = tv%T(i,j+1,k)

1241 stl = s_t(i,j,k) ; sbl = s_b(i,j,k) ; str = s_t(i,j+1,k) ; sbr = s_b(i,j+1,k)

1242 sml = tv%S(i,j,k) ; smr = tv%S(i,j+1,k)

1243 endif

1244

1245 do m=2,4

1246 w_left = wt_t(m) ; w_right = wt_b(m)

1247

1248 ! Salinity and temperature points are linearly interpolated in

1249 ! the horizontal. The subscript (1) refers to the top value in

1250 ! the vertical profile while subscript (5) refers to the bottom

1251 ! value in the vertical profile.

1252 t_top = (w_left*ttl) + (w_right*ttr)

1253 t_mn = (w_left*tml) + (w_right*tmr)

1254 t_bot = (w_left*tbl) + (w_right*tbr)

1255

1256 s_top = (w_left*stl) + (w_right*str)

1257 s_mn = (w_left*sml) + (w_right*smr)

1258 s_bot = (w_left*sbl) + (w_right*sbr)

1259

1260 ! Pressure

1261 dz_y(m,i) = (w_left*(e(i,j,k) - e(i,j,k+1))) + (w_right*(e(i,j+1,k) - e(i,j+1,k+1)))

1262

1263 pos = i*15+(m-2)*5

1264 p15(pos+1) = -gxrho * ((w_left*(e(i,j,k)-z0pres(i,j))) + (w_right*(e(i,j+1,k)-z0pres(i,j+1))))

1265 do n=2,5

1266 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i)

1267 enddo

1268

1269 ! Parabolic reconstructions in the vertical for T and S

1270 if (use_ppm) then

1271 ! Coefficients of the parabolas

1272 s6 = 3.0 * ( 2.0*s_mn - ( s_top + s_bot ) )

1273 t6 = 3.0 * ( 2.0*t_mn - ( t_top + t_bot ) )

1274 endif

1275 do n=1,5

1276 s15(pos+n) = wt_t(n) * s_top + wt_b(n) * ( s_bot + s6 * wt_t(n) )

1277 t15(pos+n) = wt_t(n) * t_top + wt_b(n) * ( t_bot + t6 * wt_t(n) )

1278 enddo

1279

1280 if (use_stanley_eos) then

1281 if (use_vart) t215(pos+1:pos+5) = (w_left*tv%varT(i,j,k)) + (w_right*tv%varT(i,j+1,k))

1282 if (use_covarts) ts15(pos+1:pos+5) = (w_left*tv%covarTS(i,j,k)) + (w_right*tv%covarTS(i,j+1,k))

1283 if (use_vars) s215(pos+1:pos+5) = (w_left*tv%varS(i,j,k)) + (w_right*tv%varS(i,j+1,k))

1284 endif

1285 enddo

1286 enddo

1287

1288 if (use_stanley_eos) then

1289 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

1290 t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), &

1291 r15(15*hi%isc+1:), eos, eosdom_h15, rho_ref=rho_ref)

1292 else

1293 call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

1294 r15(15*hi%isc+1:), eos, eosdom_h15, rho_ref=rho_ref)

1295 endif

1296

1297 do i=hi%isc,hi%iec

1298 do m=2,4

1299 ! Use Boole's rule to estimate the pressure anomaly change.

1300 pos = i*15+(m-2)*5

1301 intz(m) = (g_e*dz_y(m,i)*(c1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + &

1302 32.0*(r15(pos+2)+r15(pos+4)) + &

1303 12.0*r15(pos+3)) ))

1304 enddo ! m

1305 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)

1306

1307 ! Use Boole's rule to integrate the bottom pressure anomaly values in y.

1308 inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))

1309 enddo

1310 enddo ; endif

1311

1312end subroutine int_density_dz_generic_ppm

1313

1314!> Calls the appropriate subroutine to calculate analytical and nearly-analytical

1315!! integrals in pressure across layers of geopotential anomalies, which are

1316!! required for calculating the finite-volume form pressure accelerations in a

1317!! non-Boussinesq model. There are essentially no free assumptions, apart from the

1318!! use of Boole's rule to do the horizontal integrals, and from a truncation in the

1319!! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.

1320subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, &

1321 dza, intp_dza, intx_dza, inty_dza, halo_size, &

1322 bathyP, P_surf, dP_tiny, MassWghtInterp, &

1323 MassWghtInterpVanOnly, p_nv)

1324 type(hor_index_type), intent(in) :: hi !< The horizontal index structure

1325 real, dimension(SZI_(HI),SZJ_(HI)), &

1326 intent(in) :: t !< Potential temperature referenced to the surface [C ~> degC]

1327 real, dimension(SZI_(HI),SZJ_(HI)), &

1328 intent(in) :: s !< Salinity [S ~> ppt]

1329 real, dimension(SZI_(HI),SZJ_(HI)), &

1330 intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]

1331 real, dimension(SZI_(HI),SZJ_(HI)), &

1332 intent(in) :: p_b !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa]

1333 real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out

1334 !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]

1335 !! The calculation is mathematically identical with different values of

1336 !! alpha_ref, but this reduces the effects of roundoff.

1337 type(eos_type), intent(in) :: eos !< Equation of state structure

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

1339 real, dimension(SZI_(HI),SZJ_(HI)), &

1340 intent(inout) :: dza !< The change in the geopotential anomaly across

1341 !! the layer [L2 T-2 ~> m2 s-2]

1342 real, dimension(SZI_(HI),SZJ_(HI)), &

1343 optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the

1344 !! geopotential anomaly relative to the anomaly at the bottom of the

1345 !! layer [R L4 T-4 ~> Pa m2 s-2]

1346 real, dimension(SZIB_(HI),SZJ_(HI)), &

1347 optional, intent(inout) :: intx_dza !< The integral in x of the difference between the

1348 !! geopotential anomaly at the top and bottom of the layer divided by

1349 !! the x grid spacing [L2 T-2 ~> m2 s-2]

1350 real, dimension(SZI_(HI),SZJB_(HI)), &

1351 optional, intent(inout) :: inty_dza !< The integral in y of the difference between the

1352 !! geopotential anomaly at the top and bottom of the layer divided by

1353 !! the y grid spacing [L2 T-2 ~> m2 s-2]

1354 integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.

1355 real, dimension(SZI_(HI),SZJ_(HI)), &

1356 optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]

1357 real, dimension(SZI_(HI),SZJ_(HI)), &

1358 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]

1359 real, optional, intent(in) :: dp_tiny !< A minuscule pressure change with

1360 !! the same units as p_t [R L2 T-2 ~> Pa]

1361 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

1362 !! mass weighting to interpolate T/S in integrals

1363 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

1364 !! of T/S unless one side smaller than h_nv (i.e. vanished)

1365 real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]

1366

1367

1368 if (eos_quadrature(eos)) then

1369 call int_spec_vol_dp_generic_pcm(t, s, p_t, p_b, alpha_ref, hi, eos, us, &

1370 dza, intp_dza, intx_dza, inty_dza, halo_size, &

1371 bathyp, p_surf, dp_tiny, masswghtinterp, &

1372 masswghtinterpvanonly, p_nv)

1373 else

1374 call analytic_int_specific_vol_dp(t, s, p_t, p_b, alpha_ref, hi, eos, &

1375 dza, intp_dza, intx_dza, inty_dza, halo_size, &

1376 bathyp, p_surf, dp_tiny, masswghtinterp)

1377 endif

1378

1379end subroutine int_specific_vol_dp

1380

1381

1382!> This subroutine calculates integrals of specific volume anomalies in

1383!! pressure across layers, which are required for calculating the finite-volume

1384!! form pressure accelerations in a non-Boussinesq model. There are essentially

1385!! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.

1386subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, dza, &

1387 intp_dza, intx_dza, inty_dza, halo_size, &

1388 bathyP, P_surf, dP_neglect, MassWghtInterp, &

1389 MassWghtInterpVanOnly, p_nv)

1390 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.

1391 real, dimension(SZI_(HI),SZJ_(HI)), &

1392 intent(in) :: t !< Potential temperature of the layer [C ~> degC]

1393 real, dimension(SZI_(HI),SZJ_(HI)), &

1394 intent(in) :: s !< Salinity of the layer [S ~> ppt]

1395 real, dimension(SZI_(HI),SZJ_(HI)), &

1396 intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa]

1397 real, dimension(SZI_(HI),SZJ_(HI)), &

1398 intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa]

1399 real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out

1400 !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]

1401 !! The calculation is mathematically identical with different values of

1402 !! alpha_ref, but alpha_ref alters the effects of roundoff, and

1403 !! answers do change.

1404 type(eos_type), intent(in) :: eos !< Equation of state structure

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

1406 real, dimension(SZI_(HI),SZJ_(HI)), &

1407 intent(inout) :: dza !< The change in the geopotential anomaly

1408 !! across the layer [L2 T-2 ~> m2 s-2]

1409 real, dimension(SZI_(HI),SZJ_(HI)), &

1410 optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of

1411 !! the geopotential anomaly relative to the anomaly at the bottom of the

1412 !! layer [R L4 T-4 ~> Pa m2 s-2]

1413 real, dimension(SZIB_(HI),SZJ_(HI)), &

1414 optional, intent(inout) :: intx_dza !< The integral in x of the difference between

1415 !! the geopotential anomaly at the top and bottom of the layer divided

1416 !! by the x grid spacing [L2 T-2 ~> m2 s-2]

1417 real, dimension(SZI_(HI),SZJB_(HI)), &

1418 optional, intent(inout) :: inty_dza !< The integral in y of the difference between

1419 !! the geopotential anomaly at the top and bottom of the layer divided

1420 !! by the y grid spacing [L2 T-2 ~> m2 s-2]

1421 integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.

1422 real, dimension(SZI_(HI),SZJ_(HI)), &

1423 optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]

1424 real, dimension(SZI_(HI),SZJ_(HI)), &

1425 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]

1426 real, optional, intent(in) :: dp_neglect !< A minuscule pressure change with

1427 !! the same units as p_t [R L2 T-2 ~> Pa]

1428 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

1429 !! mass weighting to interpolate T/S in integrals

1430 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

1431 !! of T/S unless one side smaller than h_nv (i.e. vanished)

1432 real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]

1433

1434! This subroutine calculates analytical and nearly-analytical integrals in

1435! pressure across layers of geopotential anomalies, which are required for

1436! calculating the finite-volume form pressure accelerations in a non-Boussinesq

1437! model. There are essentially no free assumptions, apart from the use of

1438! Boole's rule to do the horizontal integrals, and from a truncation in the

1439! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.

1440

1441 ! Local variables

1442 real :: t5((5*hi%isd+1):(5*(hi%ied+2))) ! Temperatures along a line of subgrid locations [C ~> degC]

1443 real :: s5((5*hi%isd+1):(5*(hi%ied+2))) ! Salinities along a line of subgrid locations [S ~> ppt]

1444 real :: p5((5*hi%isd+1):(5*(hi%ied+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]

1445 real :: a5((5*hi%isd+1):(5*(hi%ied+2))) ! Specific volumes anomalies along a line of subgrid

1446 ! locations [R-1 ~> m3 kg-1]

1447 real :: t15((15*hi%isd+1):(15*(hi%ied+1))) ! Temperatures at an array of subgrid locations [C ~> degC]

1448 real :: s15((15*hi%isd+1):(15*(hi%ied+1))) ! Salinities at an array of subgrid locations [S ~> ppt]

1449 real :: p15((15*hi%isd+1):(15*(hi%ied+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]

1450 real :: a15((15*hi%isd+1):(15*(hi%ied+1))) ! Specific volumes at an array of subgrid locations [R ~> kg m-3]

1451 real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]

1452 real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]

1453 real :: dp_x(5,szib_(hi)) ! The pressure change through a layer along an x-line of subgrid locations [Z ~> m]

1454 real :: dp_y(5,szi_(hi)) ! The pressure change through a layer along a y-line of subgrid locations [Z ~> m]

1455 real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]

1456 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]

1457 real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]

1458 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]

1459 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]

1460 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]

1461 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]

1462 real :: intp(5) ! The integrals of specific volume with pressure at the

1463 ! 5 sub-column locations [L2 T-2 ~> m2 s-2]

1464 logical :: do_massweight ! Indicates whether to do mass weighting.

1465 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface

1466 logical :: massweight_bug ! If true, use an incorrect expression to determine where to apply mass weighting

1467 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

1468 ! if at least one side vanished (0 or 1) [nondim]

1469 real :: p_nonvanished ! nonvanished pressure [R L2 T-2 ~> Pa]

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

1471 integer, dimension(2) :: eosdom_h5 ! The 5-point h-point i-computational domain for the equation of state

1472 integer, dimension(2) :: eosdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state

1473 integer, dimension(2) :: eosdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state

1474 integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, n, pos, halo

1475

1476 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

1477 halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)

1478 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo

1479 if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh) ; endif

1480 if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh) ; endif

1481

1482 do_massweight = .false. ; massweight_bug = .false. ; top_massweight = .false.

1483 if (present(masswghtinterp)) then

1484 do_massweight = btest(masswghtinterp, 0) ! True for odd values

1485 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set

1486 massweight_bug = btest(masswghtinterp, 3) ! True if the 8 bit is set

1487 if (do_massweight .and. .not.present(bathyp)) call mom_error(fatal, &

1488 "int_spec_vol_dp_generic_pcm: bathyP must be present if near-bottom mass weighting is in use.")

1489 if (top_massweight .and. .not.present(p_surf)) call mom_error(fatal, &

1490 "int_spec_vol_dp_generic_pcm: P_surf must be present if near-surface mass weighting is in use.")

1491 if ((do_massweight .or. top_massweight) .and. .not.present(dp_neglect)) call mom_error(fatal, &

1492 "int_spec_vol_dp_generic_pcm: dP_neglect must be present if mass weighting is in use.")

1493 endif

1494 massweightnvonlytoggle = 1.

1495 if (present(masswghtinterpvanonly)) then

1496 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

1497 endif

1498 p_nonvanished = 0.

1499 if (present(p_nv)) then

1500 p_nonvanished = p_nv

1501 endif

1502

1503

1504 ! Set the loop ranges for equation of state calculations at various points.

1505 eosdom_h5(1) = 1 ; eosdom_h5(2) = 5*(ieh-ish+1)

1506 eosdom_q15(1) = 1 ; eosdom_q15(2) = 15*(ieq-isq+1)

1507 eosdom_h15(1) = 1 ; eosdom_h15(2) = 15*(hi%iec-hi%isc+1)

1508

1509 do j=jsh,jeh

1510 do i=ish,ieh

1511 dp = p_b(i,j) - p_t(i,j)

1512 pos = 5*i

1513 do n=1,5

1514 t5(pos+n) = t(i,j) ; s5(pos+n) = s(i,j)

1515 p5(pos+n) = p_b(i,j) - 0.25*real(n-1)*dp

1516 enddo

1517 enddo

1518

1519 call calculate_spec_vol(t5(5*ish+1:), s5(5*ish+1:), p5(5*ish+1:), a5(5*ish+1:), eos, &

1520 eosdom_h5, spv_ref=alpha_ref)

1521

1522 do i=ish,ieh

1523 dp = p_b(i,j) - p_t(i,j)

1524 ! Use Boole's rule to estimate the interface height anomaly change.

1525 pos = 5*i

1526 alpha_anom = c1_90*(7.0*(a5(pos+1)+a5(pos+5)) + 32.0*(a5(pos+2)+a5(pos+4)) + 12.0*a5(pos+3))

1527 dza(i,j) = dp*alpha_anom

1528 ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of

1529 ! the interface height anomaly.

1530 if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &

1531 (alpha_anom - c1_90*(16.0*(a5(pos+4)-a5(pos+2)) + 7.0*(a5(pos+5)-a5(pos+1))) )

1532 enddo

1533 enddo

1534

1535 if (present(intx_dza)) then ; do j=hi%jsc,hi%jec

1536 do i=isq,ieq

1537 ! hWght is the distance measure by which the cell is violation of

1538 ! hydrostatic consistency. For large hWght we bias the interpolation of

1539 ! T & S along the top and bottom integrals, akin to thickness weighting.

1540 hwght = 0.0

1541 if (do_massweight .and. massweight_bug) then

1542 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))

1543 elseif (do_massweight) then

1544 hwght = max(0., p_t(i+1,j)-bathyp(i,j), p_t(i,j)-bathyp(i+1,j))

1545 endif

1546 if (top_massweight) &

1547 hwght = max(hwght, p_surf(i,j)-p_b(i+1,j), p_surf(i+1,j)-p_b(i,j))

1548 ! If both sides are nonvanished, then set it back to zero.

1549 if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i+1,j) - p_t(i+1,j)) > p_nonvanished)) then

1550 hwght = massweightnvonlytoggle * hwght

1551 endif

1552

1553 if (hwght > 0.) then

1554 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect

1555 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect

1556 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

1557 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

1558 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom

1559 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom

1560 else

1561 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0

1562 endif

1563

1564 do m=2,4

1565 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l

1566 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)

1567 pos = i*15+(m-2)*5

1568

1569 ! T, S, and p are interpolated in the horizontal. The p interpolation

1570 ! is linear, but for T and S it may be thickness weighted.

1571 p15(pos+1) = (wt_l*p_b(i,j)) + (wt_r*p_b(i+1,j))

1572 dp_x(m,i) = (wt_l*(p_b(i,j) - p_t(i,j))) + (wt_r*(p_b(i+1,j) - p_t(i+1,j)))

1573 t15(pos+1) = (wtt_l*t(i,j)) + (wtt_r*t(i+1,j))

1574 s15(pos+1) = (wtt_l*s(i,j)) + (wtt_r*s(i+1,j))

1575

1576 do n=2,5

1577 t15(pos+n) = t15(pos+1) ; s15(pos+n) = s15(pos+1)

1578 p15(pos+n) = p15(pos+n-1) - 0.25*dp_x(m,i)

1579 enddo

1580 enddo

1581 enddo

1582

1583 call calculate_spec_vol(t15(15*isq+1:), s15(15*isq+1:), p15(15*isq+1:), &

1584 a15(15*isq+1:), eos, eosdom_q15, spv_ref=alpha_ref)

1585

1586 do i=isq,ieq

1587 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)

1588 ! Use Boole's rule to estimate the interface height anomaly change.

1589 do m=2,4

1590 pos = i*15+(m-2)*5

1591 intp(m) = (dp_x(m,i)*( c1_90*(7.0*(a15(pos+1)+a15(pos+5)) + 32.0*(a15(pos+2)+a15(pos+4)) + &

1592 12.0*a15(pos+3)) ))

1593 enddo

1594 ! Use Boole's rule to integrate the interface height anomaly values in x.

1595 intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &

1596 12.0*intp(3))

1597 enddo

1598 enddo ; endif

1599

1600 if (present(inty_dza)) then ; do j=jsq,jeq

1601 do i=hi%isc,hi%iec

1602 ! hWght is the distance measure by which the cell is violation of

1603 ! hydrostatic consistency. For large hWght we bias the interpolation of

1604 ! T & S along the top and bottom integrals, akin to thickness weighting.

1605 hwght = 0.0

1606 if (do_massweight .and. massweight_bug) then

1607 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))

1608 elseif (do_massweight) then

1609 hwght = max(0., p_t(i,j+1)-bathyp(i,j), p_t(i,j)-bathyp(i,j+1))

1610 endif

1611 if (top_massweight) &

1612 hwght = max(hwght, p_surf(i,j)-p_b(i,j+1), p_surf(i,j+1)-p_b(i,j))

1613 ! If both sides are nonvanished, then set it back to zero.

1614 if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i,j+1) - p_t(i,j+1)) > p_nonvanished)) then

1615 hwght = massweightnvonlytoggle * hwght

1616 endif

1617 if (hwght > 0.) then

1618 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect

1619 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect

1620 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

1621 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

1622 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom

1623 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom

1624 else

1625 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0

1626 endif

1627

1628 do m=2,4

1629 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l

1630 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)

1631 pos = i*15+(m-2)*5

1632

1633 ! T, S, and p are interpolated in the horizontal. The p interpolation

1634 ! is linear, but for T and S it may be thickness weighted.

1635 p15(pos+1) = (wt_l*p_b(i,j)) + (wt_r*p_b(i,j+1))

1636 dp_y(m,i) = (wt_l*(p_b(i,j) - p_t(i,j))) + (wt_r*(p_b(i,j+1) - p_t(i,j+1)))

1637 t15(pos+1) = (wtt_l*t(i,j)) + (wtt_r*t(i,j+1))

1638 s15(pos+1) = (wtt_l*s(i,j)) + (wtt_r*s(i,j+1))

1639 do n=2,5

1640 t15(pos+n) = t15(pos+1) ; s15(pos+n) = s15(pos+1)

1641 p15(pos+n) = p15(pos+n-1) - 0.25*dp_y(m,i)

1642 enddo

1643 enddo

1644 enddo

1645

1646 call calculate_spec_vol(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

1647 a15(15*hi%isc+1:), eos, eosdom_h15, spv_ref=alpha_ref)

1648

1649 do i=hi%isc,hi%iec

1650

1651 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)

1652 ! Use Boole's rule to estimate the interface height anomaly change.

1653 do m=2,4

1654 pos = i*15+(m-2)*5

1655 intp(m) = (dp_y(m,i)*( c1_90*(7.0*(a15(pos+1)+a15(pos+5)) + 32.0*(a15(pos+2)+a15(pos+4)) + &

1656 12.0*a15(pos+3)) ))

1657 enddo

1658 ! Use Boole's rule to integrate the interface height anomaly values in y.

1659 inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &

1660 12.0*intp(3))

1661 enddo

1662 enddo ; endif

1663

1664end subroutine int_spec_vol_dp_generic_pcm

1665

1666!> This subroutine calculates integrals of specific volume anomalies in

1667!! pressure across layers, which are required for calculating the finite-volume

1668!! form pressure accelerations in a non-Boussinesq model. There are essentially

1669!! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.

1670subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, &

1671 dP_neglect, bathyP, HI, EOS, US, dza, &

1672 intp_dza, intx_dza, inty_dza, P_surf, MassWghtInterp, &

1673 MassWghtInterpVanOnly, p_nv)

1674 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.

1675 real, dimension(SZI_(HI),SZJ_(HI)), &

1676 intent(in) :: t_t !< Potential temperature at the top of the layer [C ~> degC]

1677 real, dimension(SZI_(HI),SZJ_(HI)), &

1678 intent(in) :: t_b !< Potential temperature at the bottom of the layer [C ~> degC]

1679 real, dimension(SZI_(HI),SZJ_(HI)), &

1680 intent(in) :: s_t !< Salinity at the top the layer [S ~> ppt]

1681 real, dimension(SZI_(HI),SZJ_(HI)), &

1682 intent(in) :: s_b !< Salinity at the bottom the layer [S ~> ppt]

1683 real, dimension(SZI_(HI),SZJ_(HI)), &

1684 intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa]

1685 real, dimension(SZI_(HI),SZJ_(HI)), &

1686 intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa]

1687 real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out

1688 !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]

1689 !! The calculation is mathematically identical with different values of

1690 !! alpha_ref, but alpha_ref alters the effects of roundoff, and

1691 !! answers do change.

1692 real, intent(in) :: dp_neglect !<!< A miniscule pressure change with

1693 !! the same units as p_t [R L2 T-2 ~> Pa]

1694 real, dimension(SZI_(HI),SZJ_(HI)), &

1695 intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]

1696 type(eos_type), intent(in) :: eos !< Equation of state structure

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

1698 real, dimension(SZI_(HI),SZJ_(HI)), &

1699 intent(inout) :: dza !< The change in the geopotential anomaly

1700 !! across the layer [L2 T-2 ~> m2 s-2]

1701 real, dimension(SZI_(HI),SZJ_(HI)), &

1702 optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of

1703 !! the geopotential anomaly relative to the anomaly at the bottom of the

1704 !! layer [R L4 T-4 ~> Pa m2 s-2]

1705 real, dimension(SZIB_(HI),SZJ_(HI)), &

1706 optional, intent(inout) :: intx_dza !< The integral in x of the difference between

1707 !! the geopotential anomaly at the top and bottom of the layer divided

1708 !! by the x grid spacing [L2 T-2 ~> m2 s-2]

1709 real, dimension(SZI_(HI),SZJB_(HI)), &

1710 optional, intent(inout) :: inty_dza !< The integral in y of the difference between

1711 !! the geopotential anomaly at the top and bottom of the layer divided

1712 !! by the y grid spacing [L2 T-2 ~> m2 s-2]

1713 real, dimension(SZI_(HI),SZJ_(HI)), &

1714 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]

1715 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

1716 !! mass weighting to interpolate T/S in integrals

1717 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

1718 !! of T/S unless one side smaller than h_nv (i.e. vanished)

1719 real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]

1720

1721! This subroutine calculates analytical and nearly-analytical integrals in

1722! pressure across layers of geopotential anomalies, which are required for

1723! calculating the finite-volume form pressure accelerations in a non-Boussinesq

1724! model. There are essentially no free assumptions, apart from the use of

1725! Boole's rule to do the horizontal integrals, and from a truncation in the

1726! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.

1727

1728 real :: t5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Temperatures along a line of subgrid locations [C ~> degC]

1729 real :: s5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Salinities along a line of subgrid locations [S ~> ppt]

1730 real :: p5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]

1731 real :: a5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Specific volumes anomalies along a line of subgrid

1732 ! locations [R-1 ~> m3 kg-1]

1733 real :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Temperatures at an array of subgrid locations [C ~> degC]

1734 real :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Salinities at an array of subgrid locations [S ~> ppt]

1735 real :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]

1736 real :: a15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Specific volumes at an array of subgrid locations [R ~> kg m-3]

1737 real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim]

1738 real :: t_top, t_bot ! Horizontally interpolated temperature at the cell top and bottom [C ~> degC]

1739 real :: s_top, s_bot ! Horizontally interpolated salinity at the cell top and bottom [S ~> ppt]

1740 real :: p_top, p_bot ! Horizontally interpolated pressure at the cell top and bottom [R L2 T-2 ~> Pa]

1741

1742 real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]

1743 real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]

1744 real :: dp_90(2:4,szib_(hi)) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa]

1745 real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]

1746 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]

1747 real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]

1748 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]

1749 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]

1750 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]

1751 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]

1752 real :: intp(5) ! The integrals of specific volume with pressure at the

1753 ! 5 sub-column locations [L2 T-2 ~> m2 s-2]

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

1755 logical :: do_massweight ! Indicates whether to do mass weighting.

1756 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface

1757 logical :: massweight_bug ! If true, use an incorrect expression to determine where to apply mass weighting

1758 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

1759 ! if at least one side vanished (0 or 1) [nondim]

1760 real :: p_nonvanished ! nonvanished pressure [R L2 T-2 ~> Pa]

1761 integer, dimension(2) :: eosdom_h5 ! The 5-point h-point i-computational domain for the equation of state

1762 integer, dimension(2) :: eosdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state

1763 integer, dimension(2) :: eosdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state

1764 integer :: isq, ieq, jsq, jeq, i, j, m, n, pos

1765

1766 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

1767

1768 do_massweight = .false. ; massweight_bug = .false. ; top_massweight = .false.

1769 if (present(masswghtinterp)) then

1770 do_massweight = btest(masswghtinterp, 0) ! True for odd values

1771 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set

1772 massweight_bug = btest(masswghtinterp, 3) ! True if the 8 bit is set

1773 if (top_massweight .and. .not.present(p_surf)) call mom_error(fatal, &

1774 "int_spec_vol_dp_generic_plm: P_surf must be present if near-surface mass weighting is in use.")

1775 endif

1776 massweightnvonlytoggle = 1.

1777 if (present(masswghtinterpvanonly)) then

1778 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

1779 endif

1780 p_nonvanished = 0.

1781 if (present(p_nv)) then

1782 p_nonvanished = p_nv

1783 endif

1784

1785 do n = 1, 5 ! Note that these are reversed from int_density_dz.

1786 wt_t(n) = 0.25 * real(n-1)

1787 wt_b(n) = 1.0 - wt_t(n)

1788 enddo

1789

1790 ! Set the loop ranges for equation of state calculations at various points.

1791 eosdom_h5(1) = 1 ; eosdom_h5(2) = 5*(ieq-isq+2)

1792 eosdom_q15(1) = 1 ; eosdom_q15(2) = 15*(ieq-isq+1)

1793 eosdom_h15(1) = 1 ; eosdom_h15(2) = 15*(hi%iec-hi%isc+1)

1794

1795 ! 1. Compute vertical integrals

1796 do j=jsq,jeq+1

1797 do i=isq,ieq+1

1798 do n=1,5 ! T, S and p are linearly interpolated in the vertical.

1799 p5(i*5+n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j)

1800 s5(i*5+n) = wt_t(n) * s_t(i,j) + wt_b(n) * s_b(i,j)

1801 t5(i*5+n) = wt_t(n) * t_t(i,j) + wt_b(n) * t_b(i,j)

1802 enddo

1803 enddo

1804 call calculate_spec_vol(t5, s5, p5, a5, eos, eosdom_h5, spv_ref=alpha_ref)

1805 do i=isq,ieq+1

1806 ! Use Boole's rule to estimate the interface height anomaly change.

1807 dp = p_b(i,j) - p_t(i,j)

1808 alpha_anom = c1_90*((7.0*(a5(i*5+1)+a5(i*5+5)) + 32.0*(a5(i*5+2)+a5(i*5+4))) + 12.0*a5(i*5+3))

1809 dza(i,j) = dp*alpha_anom

1810 ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of

1811 ! the interface height anomaly.

1812 if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &

1813 (alpha_anom - c1_90*(16.0*(a5(i*5+4)-a5(i*5+2)) + 7.0*(a5(i*5+5)-a5(i*5+1))) )

1814 enddo

1815 enddo

1816

1817 ! 2. Compute horizontal integrals in the x direction

1818 if (present(intx_dza)) then ; do j=hi%jsc,hi%jec

1819 do i=isq,ieq

1820 ! hWght is the distance measure by which the cell is violation of

1821 ! hydrostatic consistency. For large hWght we bias the interpolation

1822 ! of T,S along the top and bottom integrals, almost like thickness

1823 ! weighting. Note: To work in terrain following coordinates we could

1824 ! offset this distance by the layer thickness to replicate other models.

1825 hwght = 0.0

1826 if (do_massweight .and. massweight_bug) then

1827 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))

1828 elseif (do_massweight) then

1829 hwght = max(0., p_t(i+1,j)-bathyp(i,j), p_t(i,j)-bathyp(i+1,j))

1830 endif

1831 if (top_massweight) &

1832 hwght = max(hwght, p_surf(i,j)-p_b(i+1,j), p_surf(i+1,j)-p_b(i,j))

1833 ! If both sides are nonvanished, then set it back to zero.

1834 if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i+1,j) - p_t(i+1,j)) > p_nonvanished)) then

1835 hwght = massweightnvonlytoggle * hwght

1836 endif

1837 if (hwght > 0.) then

1838 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect

1839 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect

1840 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

1841 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

1842 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom

1843 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom

1844 else

1845 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0

1846 endif

1847

1848 do m=2,4

1849 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l

1850 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)

1851

1852 ! T, S, and p are interpolated in the horizontal. The p interpolation

1853 ! is linear, but for T and S it may be thickness weighted.

1854 p_top = (wt_l*p_t(i,j)) + (wt_r*p_t(i+1,j))

1855 p_bot = (wt_l*p_b(i,j)) + (wt_r*p_b(i+1,j))

1856 t_top = (wtt_l*t_t(i,j)) + (wtt_r*t_t(i+1,j))

1857 t_bot = (wtt_l*t_b(i,j)) + (wtt_r*t_b(i+1,j))

1858 s_top = (wtt_l*s_t(i,j)) + (wtt_r*s_t(i+1,j))

1859 s_bot = (wtt_l*s_b(i,j)) + (wtt_r*s_b(i+1,j))

1860 dp_90(m,i) = c1_90*(p_bot - p_top)

1861

1862 ! Salinity, temperature and pressure with linear interpolation in the vertical.

1863 pos = i*15+(m-2)*5

1864 do n=1,5

1865 p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot

1866 s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot

1867 t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot

1868 enddo

1869 enddo

1870 enddo

1871

1872 call calculate_spec_vol(t15, s15, p15, a15, eos, eosdom_q15, spv_ref=alpha_ref)

1873

1874 do i=isq,ieq

1875 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)

1876 do m=2,4

1877 ! Use Boole's rule to estimate the interface height anomaly change.

1878 ! The integrals at the ends of the segment are already known.

1879 pos = i*15+(m-2)*5

1880 intp(m) = (dp_90(m,i)*((7.0*(a15(pos+1)+a15(pos+5)) + &

1881 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3) ))

1882 enddo

1883 ! Use Boole's rule to integrate the interface height anomaly values in x.

1884 intx_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &

1885 12.0*intp(3))

1886 enddo

1887 enddo ; endif

1888

1889 ! 3. Compute horizontal integrals in the y direction

1890 if (present(inty_dza)) then ; do j=jsq,jeq

1891 do i=hi%isc,hi%iec

1892 ! hWght is the distance measure by which the cell is violation of

1893 ! hydrostatic consistency. For large hWght we bias the interpolation

1894 ! of T,S along the top and bottom integrals, like thickness weighting.

1895 hwght = 0.0

1896 if (do_massweight .and. massweight_bug) then

1897 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))

1898 elseif (do_massweight) then

1899 hwght = max(0., p_t(i,j+1)-bathyp(i,j), p_t(i,j)-bathyp(i,j+1))

1900 endif

1901 if (top_massweight) &

1902 hwght = max(hwght, p_surf(i,j)-p_b(i,j+1), p_surf(i,j+1)-p_b(i,j))

1903 ! If both sides are nonvanished, then set it back to zero.

1904 if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i,j+1) - p_t(i,j+1)) > p_nonvanished)) then

1905 hwght = massweightnvonlytoggle * hwght

1906 endif

1907 if (hwght > 0.) then

1908 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect

1909 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect

1910 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

1911 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

1912 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom

1913 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom

1914 else

1915 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0

1916 endif

1917

1918 do m=2,4

1919 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l

1920 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)

1921

1922 ! T, S, and p are interpolated in the horizontal. The p interpolation

1923 ! is linear, but for T and S it may be thickness weighted.

1924 p_top = (wt_l*p_t(i,j)) + (wt_r*p_t(i,j+1))

1925 p_bot = (wt_l*p_b(i,j)) + (wt_r*p_b(i,j+1))

1926 t_top = (wtt_l*t_t(i,j)) + (wtt_r*t_t(i,j+1))

1927 t_bot = (wtt_l*t_b(i,j)) + (wtt_r*t_b(i,j+1))

1928 s_top = (wtt_l*s_t(i,j)) + (wtt_r*s_t(i,j+1))

1929 s_bot = (wtt_l*s_b(i,j)) + (wtt_r*s_b(i,j+1))

1930 dp_90(m,i) = c1_90*(p_bot - p_top)

1931

1932 ! Salinity, temperature and pressure with linear interpolation in the vertical.

1933 pos = i*15+(m-2)*5

1934 do n=1,5

1935 p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot

1936 s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot

1937 t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot

1938 enddo

1939 enddo

1940 enddo

1941

1942 call calculate_spec_vol(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &

1943 a15(15*hi%isc+1:), eos, eosdom_h15, spv_ref=alpha_ref)

1944

1945 do i=hi%isc,hi%iec

1946 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)

1947 do m=2,4

1948 ! Use Boole's rule to estimate the interface height anomaly change.

1949 ! The integrals at the ends of the segment are already known.

1950 pos = i*15+(m-2)*5

1951 intp(m) = (dp_90(m,i) * ((7.0*(a15(pos+1)+a15(pos+5)) + &

1952 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)))

1953 enddo

1954 ! Use Boole's rule to integrate the interface height anomaly values in x.

1955 inty_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &

1956 12.0*intp(3))

1957 enddo

1958 enddo ; endif

1959

1960end subroutine int_spec_vol_dp_generic_plm

1961

1962

1963!> Diagnose the fractional mass weighting in a layer that might be used with a Boussinesq calculation.

1964subroutine diagnose_mass_weight_z(z_t, z_b, bathyT, SSH, dz_neglect, MassWghtInterp, HI, &

1965 MassWt_u, MassWt_v, MassWghtInterpVanOnly, h_nv)

1966 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.

1967 real, dimension(SZI_(HI),SZJ_(HI)), &

1968 intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]

1969 real, dimension(SZI_(HI),SZJ_(HI)), &

1970 intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]

1971 real, dimension(SZI_(HI),SZJ_(HI)), &

1972 intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]

1973 real, dimension(SZI_(HI),SZJ_(HI)), &

1974 intent(in) :: ssh !< The sea surface height [Z ~> m]

1975 real, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]

1976 integer, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

1977 !! mass weighting to interpolate T/S in integrals

1978 real, dimension(SZIB_(HI),SZJ_(HI)), &

1979 intent(inout) :: masswt_u !< The fractional mass weighting at u-points [nondim]

1980 real, dimension(SZI_(HI),SZJB_(HI)), &

1981 intent(inout) :: masswt_v !< The fractional mass weighting at v-points [nondim]

1982 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

1983 !! of T/S unless one side smaller than h_nv (i.e. vanished)

1984 real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]

1985

1986 ! Local variables

1987 real :: hwght ! A pressure-thickness below topography [Z ~> m]

1988 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m]

1989 real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2]

1990 logical :: do_massweight ! Indicates whether to do mass weighting near bathymetry

1991 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface

1992 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

1993 real :: h_nonvanished ! nonvanished height [Z ~> m]

1994 integer :: isq, ieq, jsq, jeq, i, j

1995

1996 isq = hi%IscB ; ieq = hi%IecB

1997 jsq = hi%JscB ; jeq = hi%JecB

1998

1999 do_massweight = btest(masswghtinterp, 0) ! True for odd values

2000 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set

2001 massweightnvonlytoggle = 1.

2002 if (present(masswghtinterpvanonly)) then

2003 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

2004 endif

2005 h_nonvanished = 0.

2006 if (present(h_nv)) then

2007 h_nonvanished = h_nv

2008 endif

2009

2010 ! Calculate MassWt_u

2011 do j=hi%jsc,hi%jec ; do i=isq,ieq

2012 ! hWght is the distance measure by which the cell is violation of

2013 ! hydrostatic consistency. For large hWght we bias the interpolation

2014 ! of T,S along the top and bottom integrals, like thickness weighting.

2015 hwght = 0.0

2016 if (do_massweight) &

2017 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))

2018 if (top_massweight) &

2019 hwght = max(hwght, z_b(i+1,j)-ssh(i,j), z_b(i,j)-ssh(i+1,j))

2020 ! If both sides are nonvanished, then set it back to zero.

2021 if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i+1,j) - z_b(i+1,j)) > h_nonvanished)) then

2022 hwght = massweightnvonlytoggle * hwght

2023 endif

2024 if (hwght > 0.) then

2025 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect

2026 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect

2027 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

2028 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

2029 masswt_u(i,j) = (hwght*hr + hwght*hl) * idenom

2030 else

2031 masswt_u(i,j) = 0.0

2032 endif

2033 enddo ; enddo

2034

2035 ! Calculate MassWt_v

2036 do j=jsq,jeq ; do i=hi%isc,hi%iec

2037 ! hWght is the distance measure by which the cell is violation of

2038 ! hydrostatic consistency. For large hWght we bias the interpolation

2039 ! of T,S along the top and bottom integrals, like thickness weighting.

2040 hwght = 0.0

2041 if (do_massweight) &

2042 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))

2043 if (top_massweight) &

2044 hwght = max(hwght, z_b(i,j+1)-ssh(i,j), z_b(i,j)-ssh(i,j+1))

2045 ! If both sides are nonvanished, then set it back to zero.

2046 if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i,j+1) - z_b(i,j+1)) > h_nonvanished)) then

2047 hwght = massweightnvonlytoggle * hwght

2048 endif

2049 if (hwght > 0.) then

2050 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect

2051 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect

2052 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

2053 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

2054 masswt_v(i,j) = (hwght*hr + hwght*hl) * idenom

2055 else

2056 masswt_v(i,j) = 0.0

2057 endif

2058 enddo ; enddo

2059

2060end subroutine diagnose_mass_weight_z

2061

2062

2063!> Diagnose the fractional mass weighting in a layer that might be used with a non-Boussinesq calculation.

2064subroutine diagnose_mass_weight_p(p_t, p_b, bathyP, P_surf, dP_neglect, MassWghtInterp, HI, &

2065 MassWt_u, MassWt_v, MassWghtInterpVanOnly, p_nv)

2066 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.

2067 real, dimension(SZI_(HI),SZJ_(HI)), &

2068 intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa]

2069 real, dimension(SZI_(HI),SZJ_(HI)), &

2070 intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa]

2071 real, intent(in) :: dp_neglect !<!< A miniscule pressure change with

2072 !! the same units as p_t [R L2 T-2 ~> Pa]

2073 real, dimension(SZI_(HI),SZJ_(HI)), &

2074 intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]

2075 real, dimension(SZI_(HI),SZJ_(HI)), &

2076 intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]

2077 integer, intent(in) :: masswghtinterp !< A flag indicating whether and how to use

2078 !! mass weighting to interpolate T/S in integrals

2079 real, dimension(SZIB_(HI),SZJ_(HI)), &

2080 intent(inout) :: masswt_u !< The fractional mass weighting at u-points [nondim]

2081 real, dimension(SZI_(HI),SZJB_(HI)), &

2082 intent(inout) :: masswt_v !< The fractional mass weighting at v-points [nondim]

2083 logical, optional, intent(in) :: masswghtinterpvanonly !< If true, does not do mass weighting

2084 !! of T/S unless one side smaller than h_nv (i.e. vanished)

2085 real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]

2086

2087 ! Local variables

2088 real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]

2089 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]

2090 real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]

2091 logical :: do_massweight ! Indicates whether to do mass weighting.

2092 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface

2093 logical :: massweight_bug ! If true, use an incorrect expression to determine where to apply mass weighting

2094 real :: massweightnvonlytoggle ! A non-dimensional toggle factor for only using mass weighting

2095 ! if at least one side vanished (0 or 1) [nondim]

2096 real :: p_nonvanished ! nonvanished pressure [R L2 T-2 ~> Pa]

2097

2098 integer :: isq, ieq, jsq, jeq, i, j

2099

2100 isq = hi%IscB ; ieq = hi%IecB

2101 jsq = hi%JscB ; jeq = hi%JecB

2102

2103 do_massweight = btest(masswghtinterp, 0) ! True for odd values

2104 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set

2105 massweight_bug = btest(masswghtinterp, 3) ! True if the 8 bit is set

2106 massweightnvonlytoggle = 1.

2107 if (present(masswghtinterpvanonly)) then

2108 if (masswghtinterpvanonly) massweightnvonlytoggle = 0.

2109 endif

2110 p_nonvanished = 0.

2111 if (present(p_nv)) then

2112 p_nonvanished = p_nv

2113 endif

2114

2115 ! Calculate MassWt_u

2116 do j=hi%jsc,hi%jec ; do i=isq,ieq

2117 ! hWght is the distance measure by which the cell is violation of

2118 ! hydrostatic consistency. For large hWght we bias the interpolation

2119 ! of T,S along the top and bottom integrals, like thickness weighting.

2120 hwght = 0.0

2121 if (do_massweight .and. massweight_bug) then

2122 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))

2123 elseif (do_massweight) then

2124 hwght = max(0., p_t(i+1,j)-bathyp(i,j), p_t(i,j)-bathyp(i+1,j))

2125 endif

2126 if (top_massweight) &

2127 hwght = max(hwght, p_surf(i,j)-p_b(i+1,j), p_surf(i+1,j)-p_b(i,j))

2128 ! If both sides are nonvanished, then set it back to zero.

2129 if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i+1,j) - p_t(i+1,j)) > p_nonvanished)) then

2130 hwght = massweightnvonlytoggle * hwght

2131 endif

2132 if (hwght > 0.) then

2133 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect

2134 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect

2135 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

2136 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

2137 masswt_u(i,j) = (hwght*hr + hwght*hl) * idenom

2138 else

2139 masswt_u(i,j) = 0.0

2140 endif

2141 enddo ; enddo

2142

2143 ! Calculate MassWt_v

2144 do j=jsq,jeq ; do i=hi%isc,hi%iec

2145 ! hWght is the distance measure by which the cell is violation of

2146 ! hydrostatic consistency. For large hWght we bias the interpolation

2147 ! of T,S along the top and bottom integrals, like thickness weighting.

2148 hwght = 0.0

2149 if (do_massweight .and. massweight_bug) then

2150 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))

2151 elseif (do_massweight) then

2152 hwght = max(0., p_t(i,j+1)-bathyp(i,j), p_t(i,j)-bathyp(i,j+1))

2153 endif

2154 if (top_massweight) &

2155 hwght = max(hwght, p_surf(i,j)-p_b(i,j+1), p_surf(i,j+1)-p_b(i,j))

2156 ! If both sides are nonvanished, then set it back to zero.

2157 if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i,j+1) - p_t(i,j+1)) > p_nonvanished)) then

2158 hwght = massweightnvonlytoggle * hwght

2159 endif

2160 if (hwght > 0.) then

2161 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect

2162 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect

2163 hwght = hwght * ( (hl-hr)/(hl+hr) )**2

2164 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )

2165 masswt_v(i,j) = (hwght*hr + hwght*hl) * idenom

2166 else

2167 masswt_v(i,j) = 0.0

2168 endif

2169 enddo ; enddo

2170

2171end subroutine diagnose_mass_weight_p

2172

2173!> Find the depth at which the reconstructed pressure matches P_tgt

2174subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &

2175 rho_ref, G_e, EOS, US, P_b, z_out, z_tol, frac_dp_bugfix)

2176 real, intent(in) :: t_t !< Potential temperature at the cell top [C ~> degC]

2177 real, intent(in) :: t_b !< Potential temperature at the cell bottom [C ~> degC]

2178 real, intent(in) :: s_t !< Salinity at the cell top [S ~> ppt]

2179 real, intent(in) :: s_b !< Salinity at the cell bottom [S ~> ppt]

2180 real, intent(in) :: z_t !< Absolute height of top of cell [Z ~> m] (Boussinesq ????)

2181 real, intent(in) :: z_b !< Absolute height of bottom of cell [Z ~> m]

2182 real, intent(in) :: p_t !< Anomalous pressure of top of cell, relative

2183 !! to g*rho_ref*z_t [R L2 T-2 ~> Pa]

2184 real, intent(in) :: p_tgt !< Target pressure at height z_out, relative

2185 !! to g*rho_ref*z_out [R L2 T-2 ~> Pa]

2186 real, intent(in) :: rho_ref !< Reference density with which calculation

2187 !! are anomalous to [R ~> kg m-3]

2188 real, intent(in) :: g_e !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]

2189 type(eos_type), intent(in) :: eos !< Equation of state structure

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

2191 real, intent(out) :: p_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]

2192 real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m]

2193 real, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m]

2194 logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos

2195

2196 ! Local variables

2197 real :: dp ! Pressure thickness of the layer [R L2 T-2 ~> Pa]

2198 real :: f_guess, f_l, f_r ! Fractional positions [nondim]

2199 real :: gxrho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]

2200 real :: pa, pa_left, pa_right, pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz [R L2 T-2 ~> Pa]

2201 integer :: m ! A counter for how many iterations have been done in the while loop

2202 character(len=240) :: msg

2203

2204 gxrho = g_e * rho_ref

2205

2206 ! Anomalous pressure difference across whole cell

2207 dp = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, 1.0, eos, frac_dp_bugfix)

2208

2209 p_b = p_t + dp ! Anomalous pressure at bottom of cell

2210

2211 if (p_tgt <= p_t ) then

2212 z_out = z_t

2213 return

2214 endif

2215

2216 if (p_tgt >= p_b) then

2217 z_out = z_b

2218 return

2219 endif

2220

2221 f_l = 0.

2222 pa_left = p_t - p_tgt ! Pa_left < 0

2223 f_r = 1.

2224 pa_right = p_b - p_tgt ! Pa_right > 0

2225 pa_tol = gxrho * z_tol

2226

2227 f_guess = f_l - pa_left / (pa_right - pa_left) * (f_r - f_l)

2228 pa = pa_right - pa_left ! To get into iterative loop

2229 m = 0 ! Reset the counter for the loop to be zero

2230 do while ( abs(pa) > pa_tol )

2231

2232 m = m + 1

2233 if (m > 30) then ! Call an error, because convergence to the tolerance has not been achieved

2234 write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt

2235 call mom_error(fatal, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg)

2236 endif

2237 z_out = z_t + ( z_b - z_t ) * f_guess

2238 pa = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, f_guess, eos, frac_dp_bugfix) - ( p_tgt - p_t )

2239

2240 if (pa<pa_left) then

2241 write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt

2242 call mom_error(fatal, 'find_depth_of_pressure_in_cell out of bounds negative: /n'//msg)

2243 elseif (pa<0.) then

2244 pa_left = pa

2245 f_l = f_guess

2246 elseif (pa>pa_right) then

2247 write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt

2248 call mom_error(fatal, 'find_depth_of_pressure_in_cell out of bounds positive: /n'//msg)

2249 elseif (pa>0.) then

2250 pa_right = pa

2251 f_r = f_guess

2252 else ! Pa == 0

2253 return

2254 endif

2255 f_guess = f_l - pa_left / (pa_right - pa_left) * (f_r - f_l)

2256

2257 enddo

2258

2259end subroutine find_depth_of_pressure_in_cell

2260

2261!> Calculate the average in situ specific volume across layers

2262subroutine avg_specific_vol(T, S, p_t, dp, HI, EOS, SpV_avg, halo_size)

2263 type(hor_index_type), intent(in) :: hi !< The horizontal index structure

2264 real, dimension(SZI_(HI),SZJ_(HI)), &

2265 intent(in) :: t !< Potential temperature of the layer [C ~> degC]

2266 real, dimension(SZI_(HI),SZJ_(HI)), &

2267 intent(in) :: s !< Salinity of the layer [S ~> ppt]

2268 real, dimension(SZI_(HI),SZJ_(HI)), &

2269 intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]

2270 real, dimension(SZI_(HI),SZJ_(HI)), &

2271 intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa]

2272 type(eos_type), intent(in) :: eos !< Equation of state structure

2273 real, dimension(SZI_(HI),SZJ_(HI)), &

2274 intent(inout) :: spv_avg !< The vertical average specific volume

2275 !! in the layer [R-1 ~> m3 kg-1]

2276 integer, optional, intent(in) :: halo_size !< The number of halo points in which to work.

2277

2278 ! Local variables

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

2280 integer :: jsh, jeh, j, halo

2281

2282 halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)

2283 jsh = hi%jsc-halo ; jeh = hi%jec+halo

2284

2285 eosdom(:) = eos_domain(hi, halo_size)

2286 do j=jsh,jeh

2287 call average_specific_vol(t(:,j), s(:,j), p_t(:,j), dp(:,j), spv_avg(:,j), eos, eosdom)

2288 enddo

2289

2290end subroutine avg_specific_vol

2291

2292!> Returns change in anomalous pressure change from top to non-dimensional

2293!! position pos between z_t and z_b [R L2 T-2 ~> Pa]

2294real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS, frac_dp_bugfix)

2295 real, intent(in) :: t_t !< Potential temperature at the cell top [C ~> degC]

2296 real, intent(in) :: t_b !< Potential temperature at the cell bottom [C ~> degC]

2297 real, intent(in) :: s_t !< Salinity at the cell top [S ~> ppt]

2298 real, intent(in) :: s_b !< Salinity at the cell bottom [S ~> ppt]

2299 real, intent(in) :: z_t !< The geometric height at the top of the layer [Z ~> m]

2300 real, intent(in) :: z_b !< The geometric height at the bottom of the layer [Z ~> m]

2301 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted out to

2302 !! reduce the magnitude of each of the integrals.

2303 real, intent(in) :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]

2304 real, intent(in) :: pos !< The fractional vertical position, 0 to 1 [nondim]

2305 type(eos_type), intent(in) :: eos !< Equation of state structure

2306 logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos

2307

2308 ! Local variables

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

2310 real :: dz ! Distance from the layer top [Z ~> m]

2311 real :: top_weight, bottom_weight ! Fractional weights at quadrature points [nondim]

2312 real :: rho_ave ! Average density [R ~> kg m-3]

2313 real, dimension(5) :: t5 ! Temperatures at quadrature points [C ~> degC]

2314 real, dimension(5) :: s5 ! Salinities at quadrature points [S ~> ppt]

2315 real, dimension(5) :: p5 ! Pressures at quadrature points [R L2 T-2 ~> Pa]

2316 real, dimension(5) :: rho5 ! Densities at quadrature points [R ~> kg m-3]

2317 integer :: n

2318

2319 do n=1,5

2320 ! Evaluate density at five quadrature points

2321 bottom_weight = 0.25*real(n-1) * pos

2322 top_weight = 1.0 - bottom_weight

2323 ! Salinity and temperature points are linearly interpolated

2324 s5(n) = top_weight * s_t + bottom_weight * s_b

2325 t5(n) = top_weight * t_t + bottom_weight * t_b

2326 if (frac_dp_bugfix) then

2327 p5(n) = (-1) * ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )

2328 else

2329 p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )

2330 endif !bugfix

2331 enddo

2332 call calculate_density(t5, s5, p5, rho5, eos)

2333 rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref

2334

2335 ! Use Boole's rule to estimate the average density

2336 rho_ave = c1_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3))

2337

2338 dz = ( z_t - z_b ) * pos

2339 frac_dp_at_pos = g_e * dz * rho_ave

2340end function frac_dp_at_pos

2341

2342end module mom_density_integrals

2343

2344!> \namespace mom_density_integrals

2345!!