MOM_continuity_PPM.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!> Solve the layer continuity equation using the PPM method for layer fluxes.

6module mom_continuity_ppm

7

8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine

9use mom_diag_mediator, only : time_type, diag_ctrl

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

11use mom_file_parser, only : get_param, log_version, param_file_type

12use mom_grid, only : ocean_grid_type

13use mom_open_boundary, only : ocean_obc_type, obc_segment_type, obc_none

14use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s

15use mom_unit_scaling, only : unit_scale_type

16use mom_variables, only : bt_cont_type, porous_barrier_type

17use mom_verticalgrid, only : verticalgrid_type

18

19implicit none ; private

20

21#include <MOM_memory.h>

22

23public continuity_ppm, continuity_ppm_init, continuity_ppm_stencil

24public continuity_fluxes, continuity_adjust_vel

25public zonal_mass_flux, meridional_mass_flux

26public zonal_edge_thickness, meridional_edge_thickness

27public continuity_zonal_convergence, continuity_merdional_convergence

28public zonal_flux_thickness, meridional_flux_thickness

29public zonal_bt_mass_flux, meridional_bt_mass_flux

30public set_continuity_loop_bounds

31

32!>@{ CPU time clock IDs

33integer :: id_clock_reconstruct, id_clock_update, id_clock_correct

34!>@}

35

36!> Control structure for mom_continuity_ppm

37type, public :: continuity_ppm_cs ; private

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

39 type(diag_ctrl), pointer :: diag !< Diagnostics control structure.

40 logical :: upwind_1st !< If true, use a first-order upwind scheme.

41 logical :: monotonic !< If true, use the Colella & Woodward monotonic

42 !! limiter; otherwise use a simple positive

43 !! definite limiter.

44 logical :: simple_2nd !< If true, use a simple second order (arithmetic

45 !! mean) interpolation of the edge values instead

46 !! of the higher order interpolation.

47 real :: tol_eta !< The tolerance for free-surface height

48 !! discrepancies between the barotropic solution and

49 !! the sum of the layer thicknesses [H ~> m or kg m-2].

50 real :: tol_vel !< The tolerance for barotropic velocity

51 !! discrepancies between the barotropic solution and

52 !! the sum of the layer thicknesses [L T-1 ~> m s-1].

53 real :: cfl_limit_adjust !< The maximum CFL of the adjusted velocities [nondim]

54 real :: h_marg_min !< Negligible floor on h_marg, the marginal thickness

55 !! used to calculate the partial derivative of transports

56 !! with velocities [H ~> m or kg m-2]

57 logical :: aggress_adjust !< If true, allow the adjusted velocities to have a

58 !! relative CFL change up to 0.5. False by default.

59 logical :: vol_cfl !< If true, use the ratio of the open face lengths

60 !! to the tracer cell areas when estimating CFL

61 !! numbers. Without aggress_adjust, the default is

62 !! false; it is always true with.

63 logical :: better_iter !< If true, stop corrective iterations using a

64 !! velocity-based criterion and only stop if the

65 !! iteration is better than all predecessors.

66 logical :: use_visc_rem_max !< If true, use more appropriate limiting bounds

67 !! for corrections in strongly viscous columns.

68 logical :: marginal_faces !< If true, use the marginal face areas from the

69 !! continuity solver for use as the weights in the

70 !! barotropic solver. Otherwise use the transport

71 !! averaged areas.

72end type continuity_ppm_cs

73

74!> A container for loop bounds

75type, public :: cont_loop_bounds_type ; private

76 !>@{ Loop bounds

77 integer :: ish, ieh, jsh, jeh

78 !>@}

79end type cont_loop_bounds_type

80

81!> Finds the thickness fluxes from the continuity solver or their vertical sum without

82!! actually updating the layer thicknesses.

83interface continuity_fluxes

84 module procedure continuity_3d_fluxes, continuity_2d_fluxes

85end interface continuity_fluxes

86

87contains

88

89!> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme,

90!! based on Lin (1994).

91subroutine continuity_ppm(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &

92 visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont, du_cor, dv_cor)

93 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

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

95 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

96 intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

97 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

98 intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

99 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

100 intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].

101 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

102 intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2].

103 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

104 intent(out) :: uh !< Zonal volume flux, u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].

105 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

106 intent(out) :: vh !< Meridional volume flux, v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].

107 real, intent(in) :: dt !< Time increment [T ~> s].

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

109 type(continuity_ppm_cs), intent(in) :: cs !< Module's control structure.

110 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

111 type(porous_barrier_type), intent(in) :: pbv !< pointers to porous barrier fractional cell metrics

112 real, dimension(SZIB_(G),SZJ_(G)), &

113 optional, intent(in) :: uhbt !< The summed volume flux through zonal faces

114 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

115 real, dimension(SZI_(G),SZJB_(G)), &

116 optional, intent(in) :: vhbt !< The summed volume flux through meridional faces

117 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

118 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

119 optional, intent(in) :: visc_rem_u

120 !< The fraction of zonal momentum originally

121 !! in a layer that remains after a time-step of viscosity, and the

122 !! fraction of a time-step's worth of a barotropic acceleration that

123 !! a layer experiences after viscosity is applied [nondim].

124 !! Visc_rem_u is between 0 (at the bottom) and 1 (far above the bottom).

125 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

126 optional, intent(in) :: visc_rem_v

127 !< The fraction of meridional momentum originally

128 !! in a layer that remains after a time-step of viscosity, and the

129 !! fraction of a time-step's worth of a barotropic acceleration that

130 !! a layer experiences after viscosity is applied [nondim].

131 !! Visc_rem_v is between 0 (at the bottom) and 1 (far above the bottom).

132 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

133 optional, intent(out) :: u_cor

134 !< The zonal velocities that give uhbt as the depth-integrated transport [L T-1 ~> m s-1].

135 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

136 optional, intent(out) :: v_cor

137 !< The meridional velocities that give vhbt as the depth-integrated

138 !! transport [L T-1 ~> m s-1].

139 type(bt_cont_type), optional, pointer :: bt_cont !< A structure with elements that describe

140 !! the effective open face areas as a function of barotropic flow.

141 real, dimension(SZIB_(G),SZJ_(G)), &

142 optional, intent(out) :: du_cor !< The zonal velocity increments from u that give uhbt

143 !! as the depth-integrated transports [L T-1 ~> m s-1].

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

145 optional, intent(out) :: dv_cor !< The meridional velocity increments from v that give vhbt

146 !! as the depth-integrated transports [L T-1 ~> m s-1].

147

148 ! Local variables

149 real :: h_w(szi_(g),szj_(g),szk_(gv)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

150 real :: h_e(szi_(g),szj_(g),szk_(gv)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

151 real :: h_s(szi_(g),szj_(g),szk_(gv)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

152 real :: h_n(szi_(g),szj_(g),szk_(gv)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

153 real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.

154 type(cont_loop_bounds_type) :: lb ! A type indicating the loop range for a phase of the updates

155 logical :: x_first

156

157 h_min = gv%Angstrom_H

158

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

160 "MOM_continuity_PPM: Module must be initialized before it is used.")

161

162 x_first = (mod(g%first_direction,2) == 0)

163

164 if (present(visc_rem_u) .neqv. present(visc_rem_v)) call mom_error(fatal, &

165 "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither "// &

166 "one must be present in call to continuity_PPM.")

167

168 if (x_first) then

169 ! First advect zonally, with loop bounds that accomodate the subsequent meridional advection.

170 lb = set_continuity_loop_bounds(g, cs, i_stencil=.false., j_stencil=.true.)

171 call zonal_edge_thickness(hin, h_w, h_e, g, gv, us, cs, obc, lb)

172 call zonal_mass_flux(u, hin, h_w, h_e, uh, dt, g, gv, us, cs, obc, pbv%por_face_areaU, &

173 lb, uhbt, visc_rem_u, u_cor, bt_cont, du_cor)

174 call continuity_zonal_convergence(h, uh, dt, g, gv, lb, hin)

175

176 ! Now advect meridionally, using the updated thicknesses to determine the fluxes.

177 lb = set_continuity_loop_bounds(g, cs, i_stencil=.false., j_stencil=.false.)

178 call meridional_edge_thickness(h, h_s, h_n, g, gv, us, cs, obc, lb)

179 call meridional_mass_flux(v, h, h_s, h_n, vh, dt, g, gv, us, cs, obc, pbv%por_face_areaV, &

180 lb, vhbt, visc_rem_v, v_cor, bt_cont, dv_cor)

181 call continuity_merdional_convergence(h, vh, dt, g, gv, lb, hmin=h_min)

182

183 else ! .not. x_first

184 ! First advect meridionally, with loop bounds that accomodate the subsequent zonal advection.

185 lb = set_continuity_loop_bounds(g, cs, i_stencil=.true., j_stencil=.false.)

186 call meridional_edge_thickness(hin, h_s, h_n, g, gv, us, cs, obc, lb)

187 call meridional_mass_flux(v, hin, h_s, h_n, vh, dt, g, gv, us, cs, obc, pbv%por_face_areaV, &

188 lb, vhbt, visc_rem_v, v_cor, bt_cont, dv_cor)

189 call continuity_merdional_convergence(h, vh, dt, g, gv, lb, hin)

190

191 ! Now advect zonally, using the updated thicknesses to determine the fluxes.

192 lb = set_continuity_loop_bounds(g, cs, i_stencil=.false., j_stencil=.false.)

193 call zonal_edge_thickness(h, h_w, h_e, g, gv, us, cs, obc, lb)

194 call zonal_mass_flux(u, h, h_w, h_e, uh, dt, g, gv, us, cs, obc, pbv%por_face_areaU, &

195 lb, uhbt, visc_rem_u, u_cor, bt_cont, du_cor)

196 call continuity_zonal_convergence(h, uh, dt, g, gv, lb, hmin=h_min)

197 endif

198

199end subroutine continuity_ppm

200

201!> Finds the thickness fluxes from the continuity solver without actually updating the

202!! layer thicknesses. Because the fluxes in the two directions are calculated based on the

203!! input thicknesses, which are not updated between the direcitons, the fluxes returned here

204!! are not the same as those that would be returned by a call to continuity.

205subroutine continuity_3d_fluxes(u, v, h, uh, vh, dt, G, GV, US, CS, OBC, pbv)

206 type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.

207 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure.

208 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

209 intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

210 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

211 intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

212 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

213 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].

214 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

215 intent(out) :: uh !< Thickness fluxes through zonal faces,

216 !! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].

217 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

218 intent(out) :: vh !< Thickness fluxes through meridional faces,

219 !! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].

220 real, intent(in) :: dt !< Time increment [T ~> s].

221 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

222 type(continuity_ppm_cs), intent(in) :: CS !< Control structure for mom_continuity.

223 type(ocean_obc_type), pointer :: OBC !< Open boundaries control structure.

224 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics

225

226 ! Local variables

227 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

228 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

229 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

230 real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

231

232 call zonal_edge_thickness(h, h_w, h_e, g, gv, us, cs, obc)

233 call zonal_mass_flux(u, h, h_w, h_e, uh, dt, g, gv, us, cs, obc, pbv%por_face_areaU)

234

235 call meridional_edge_thickness(h, h_s, h_n, g, gv, us, cs, obc)

236 call meridional_mass_flux(v, h, h_s, h_n, vh, dt, g, gv, us, cs, obc, pbv%por_face_areaV)

237

238end subroutine continuity_3d_fluxes

239

240!> Find the vertical sum of the thickness fluxes from the continuity solver without actually

241!! updating the layer thicknesses. Because the fluxes in the two directions are calculated

242!! based on the input thicknesses, which are not updated between the directions, the fluxes

243!! returned here are not the same as those that would be returned by a call to continuity.

244subroutine continuity_2d_fluxes(u, v, h, uhbt, vhbt, dt, G, GV, US, CS, OBC, pbv)

245 type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.

246 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure.

247 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

248 intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

249 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

250 intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

251 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

252 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].

253 real, dimension(SZIB_(G),SZJ_(G)), &

254 intent(out) :: uhbt !< Vertically summed thickness flux through

255 !! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].

256 real, dimension(SZI_(G),SZJB_(G)), &

257 intent(out) :: vhbt !< Vertically summed thickness flux through

258 !! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].

259 real, intent(in) :: dt !< Time increment [T ~> s].

260 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

261 type(continuity_ppm_cs), intent(in) :: CS !< Control structure for mom_continuity.

262 type(ocean_obc_type), pointer :: OBC !< Open boundaries control structure.

263 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics

264

265 ! Local variables

266 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

267 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

268 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

269 real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

270

271 call zonal_edge_thickness(h, h_w, h_e, g, gv, us, cs, obc)

272 call zonal_bt_mass_flux(u, h, h_w, h_e, uhbt, dt, g, gv, us, cs, obc, pbv%por_face_areaU)

273

274 call meridional_edge_thickness(h, h_s, h_n, g, gv, us, cs, obc)

275 call meridional_bt_mass_flux(v, h, h_s, h_n, vhbt, dt, g, gv, us, cs, obc, pbv%por_face_areaV)

276

277end subroutine continuity_2d_fluxes

278

279!> Correct the velocities to give the specified depth-integrated transports by applying a

280!! barotropic acceleration (subject to viscous drag) to the velocities.

281subroutine continuity_adjust_vel(u, v, h, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, visc_rem_u, visc_rem_v)

282 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure.

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

284 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

285 intent(inout) :: u !< Zonal velocity, which will be adjusted to

286 !! give uhbt as the depth-integrated

287 !! transport [L T-1 ~> m s-1].

288 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

289 intent(inout) :: v !< Meridional velocity, which will be adjusted

290 !! to give vhbt as the depth-integrated

291 !! transport [L T-1 ~> m s-1].

292 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

293 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].

294 real, intent(in) :: dt !< Time increment [T ~> s].

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

296 type(continuity_ppm_cs), intent(in) :: cs !< Control structure for mom_continuity.

297 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

298 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics

299 real, dimension(SZIB_(G),SZJ_(G)), &

300 intent(in) :: uhbt !< The vertically summed thickness flux through

301 !! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].

302 real, dimension(SZI_(G),SZJB_(G)), &

303 intent(in) :: vhbt !< The vertically summed thickness flux through

304 !! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].

305 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

306 optional, intent(in) :: visc_rem_u !< Both the fraction of the zonal momentum

307 !! that remains after a time-step of viscosity, and

308 !! the fraction of a time-step's worth of a barotropic

309 !! acceleration that a layer experiences after viscosity

310 !! is applied [nondim]. This goes between 0 (at the

311 !! bottom) and 1 (far above the bottom). When this

312 !! column is under an ice shelf, this also goes to 0

313 !! at the top due to the no-slip boundary condition there.

314 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

315 optional, intent(in) :: visc_rem_v !< Both the fraction of the meridional momentum

316 !! that remains after a time-step of viscosity, and

317 !! the fraction of a time-step's worth of a barotropic

318 !! acceleration that a layer experiences after viscosity

319 !! is applied [nondim]. This goes between 0 (at the

320 !! bottom) and 1 (far above the bottom). When this

321 !! column is under an ice shelf, this also goes to 0

322 !! at the top due to the no-slip boundary condition there.

323

324 ! Local variables

325 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_in !< Input zonal velocity [L T-1 ~> m s-1]

326 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_in !< Input meridional velocity [L T-1 ~> m s-1]

327 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uh !< Volume flux through zonal faces =

328 !! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].

329 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vh !< Volume flux through meridional faces =

330 !! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].

331 real :: h_w(szi_(g),szj_(g),szk_(gv)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

332 real :: h_e(szi_(g),szj_(g),szk_(gv)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]

333 real :: h_s(szi_(g),szj_(g),szk_(gv)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

334 real :: h_n(szi_(g),szj_(g),szk_(gv)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]

335

336 ! It might not be necessary to separate the input velocity array from the adjusted velocities,

337 ! but it seems safer to do so, even if it might be less efficient.

338 u_in(:,:,:) = u(:,:,:)

339 v_in(:,:,:) = v(:,:,:)

340

341 call zonal_edge_thickness(h, h_w, h_e, g, gv, us, cs, obc)

342 call zonal_mass_flux(u_in, h, h_w, h_e, uh, dt, g, gv, us, cs, obc, pbv%por_face_areaU, &

343 uhbt=uhbt, visc_rem_u=visc_rem_u, u_cor=u)

344

345 call meridional_edge_thickness(h, h_s, h_n, g, gv, us, cs, obc)

346 call meridional_mass_flux(v_in, h, h_s, h_n, vh, dt, g, gv, us, cs, obc, pbv%por_face_areaV, &

347 vhbt=vhbt, visc_rem_v=visc_rem_v, v_cor=v)

348

349end subroutine continuity_adjust_vel

350

351

352!> Updates the thicknesses due to zonal thickness fluxes.

353subroutine continuity_zonal_convergence(h, uh, dt, G, GV, LB, hin, hmin)

354 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure

355 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure

356 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

357 intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2]

358 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

359 intent(in) :: uh !< Zonal thickness flux, u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]

360 real, intent(in) :: dt !< Time increment [T ~> s]

361 type(cont_loop_bounds_type), intent(in) :: lb !< Loop bounds structure

362 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

363 optional, intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].

364 !! If hin is absent, h is also the initial thickness.

365 real, optional, intent(in) :: hmin !< The minimum layer thickness [H ~> m or kg m-2]

366

367 real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.

368 integer :: i, j, k

369

370 call cpu_clock_begin(id_clock_update)

371

372 h_min = 0.0 ; if (present(hmin)) h_min = hmin

373

374 if (present(hin)) then

375 !$OMP parallel do default(shared)

376 do k=1,gv%ke ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh

377 h(i,j,k) = max( hin(i,j,k) - dt * g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k)), h_min )

378 enddo ; enddo ; enddo

379 else

380 !$OMP parallel do default(shared)

381 do k=1,gv%ke ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh

382 h(i,j,k) = max( h(i,j,k) - dt * g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k)), h_min )

383 enddo ; enddo ; enddo

384 endif

385

386 call cpu_clock_end(id_clock_update)

387

388end subroutine continuity_zonal_convergence

389

390!> Updates the thicknesses due to meridional thickness fluxes.

391subroutine continuity_merdional_convergence(h, vh, dt, G, GV, LB, hin, hmin)

392 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure

393 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure

394 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

395 intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2]

396 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

397 intent(in) :: vh !< Meridional thickness flux, v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]

398 real, intent(in) :: dt !< Time increment [T ~> s]

399 type(cont_loop_bounds_type), intent(in) :: lb !< Loop bounds structure

400 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

401 optional, intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].

402 !! If hin is absent, h is also the initial thickness.

403 real, optional, intent(in) :: hmin !< The minimum layer thickness [H ~> m or kg m-2]

404

405 real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.

406 integer :: i, j, k

407

408 call cpu_clock_begin(id_clock_update)

409

410 h_min = 0.0 ; if (present(hmin)) h_min = hmin

411

412 if (present(hin)) then

413 !$OMP parallel do default(shared)

414 do k=1,gv%ke ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh

415 h(i,j,k) = max( hin(i,j,k) - dt * g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k)), h_min )

416 enddo ; enddo ; enddo

417 else

418 !$OMP parallel do default(shared)

419 do k=1,gv%ke ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh

420 h(i,j,k) = max( h(i,j,k) - dt * g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k)), h_min )

421 enddo ; enddo ; enddo

422 endif

423

424 call cpu_clock_end(id_clock_update)

425

426end subroutine continuity_merdional_convergence

427

428

429!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.

430subroutine zonal_edge_thickness(h_in, h_W, h_E, G, GV, US, CS, OBC, LB_in)

431 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

432 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

433 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

434 intent(in) :: h_in !< Tracer cell layer thickness [H ~> m or kg m-2].

435 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

436 intent(out) :: h_w !< Western edge layer thickness [H ~> m or kg m-2].

437 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

438 intent(out) :: h_e !< Eastern edge layer thickness [H ~> m or kg m-2].

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

440 type(continuity_ppm_cs), intent(in) :: cs !< This module's control structure.

441 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

442 type(cont_loop_bounds_type), &

443 optional, intent(in) :: lb_in !< Loop bounds structure.

444

445 ! Local variables

446 type(cont_loop_bounds_type) :: lb

447 integer :: i, j, k, ish, ieh, jsh, jeh, nz

448

449 call cpu_clock_begin(id_clock_reconstruct)

450

451 if (present(lb_in)) then

452 lb = lb_in

453 else

454 lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec

455 endif

456 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = gv%ke

457

458 if (cs%upwind_1st) then

459 !$OMP parallel do default(shared)

460 do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh+1

461 h_w(i,j,k) = h_in(i,j,k) ; h_e(i,j,k) = h_in(i,j,k)

462 enddo ; enddo ; enddo

463 else

464 !$OMP parallel do default(shared)

465 do k=1,nz

466 call ppm_reconstruction_x(h_in(:,:,k), h_w(:,:,k), h_e(:,:,k), g, lb, &

467 2.0*gv%Angstrom_H, cs%monotonic, cs%simple_2nd, obc, k)

468 enddo

469 endif

470

471 call cpu_clock_end(id_clock_reconstruct)

472

473end subroutine zonal_edge_thickness

474

475

476!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.

477subroutine meridional_edge_thickness(h_in, h_S, h_N, G, GV, US, CS, OBC, LB_in)

478 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

479 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

480 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

481 intent(in) :: h_in !< Tracer cell layer thickness [H ~> m or kg m-2].

482 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

483 intent(out) :: h_s !< Southern edge layer thickness [H ~> m or kg m-2].

484 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

485 intent(out) :: h_n !< Northern edge layer thickness [H ~> m or kg m-2].

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

487 type(continuity_ppm_cs), intent(in) :: cs !< This module's control structure.

488 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

489 type(cont_loop_bounds_type), &

490 optional, intent(in) :: lb_in !< Loop bounds structure.

491

492 ! Local variables

493 type(cont_loop_bounds_type) :: lb

494 integer :: i, j, k, ish, ieh, jsh, jeh, nz

495

496 call cpu_clock_begin(id_clock_reconstruct)

497

498 if (present(lb_in)) then

499 lb = lb_in

500 else

501 lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec

502 endif

503 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = gv%ke

504

505 if (cs%upwind_1st) then

506 !$OMP parallel do default(shared)

507 do k=1,nz ; do j=jsh-1,jeh+1 ; do i=ish,ieh

508 h_s(i,j,k) = h_in(i,j,k) ; h_n(i,j,k) = h_in(i,j,k)

509 enddo ; enddo ; enddo

510 else

511 !$OMP parallel do default(shared)

512 do k=1,nz

513 call ppm_reconstruction_y(h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), g, lb, &

514 2.0*gv%Angstrom_H, cs%monotonic, cs%simple_2nd, obc, k)

515 enddo

516 endif

517

518 call cpu_clock_end(id_clock_reconstruct)

519

520end subroutine meridional_edge_thickness

521

522

523!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.

524subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_face_areaU, &

525 LB_in, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)

526 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

527 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

528 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

529 intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

530 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

531 intent(in) :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].

532 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

533 intent(in) :: h_w !< Western edge thicknesses [H ~> m or kg m-2].

534 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

535 intent(in) :: h_e !< Eastern edge thicknesses [H ~> m or kg m-2].

536 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

537 intent(out) :: uh !< Volume flux through zonal faces = u*h*dy

538 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

539 real, intent(in) :: dt !< Time increment [T ~> s].

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

541 type(continuity_ppm_cs), intent(in) :: cs !< This module's control structure.

542 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

543 real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &

544 intent(in) :: por_face_areau !< fractional open area of U-faces [nondim]

545 type(cont_loop_bounds_type), &

546 optional, intent(in) :: lb_in !< Loop bounds structure.

547 real, dimension(SZIB_(G),SZJ_(G)), &

548 optional, intent(in) :: uhbt !< The summed volume flux through zonal faces

549 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

550 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

551 optional, intent(in) :: visc_rem_u

552 !< The fraction of zonal momentum originally in a layer that remains after a

553 !! time-step of viscosity, and the fraction of a time-step's worth of a barotropic

554 !! acceleration that a layer experiences after viscosity is applied [nondim].

555 !! Visc_rem_u is between 0 (at the bottom) and 1 (far above the bottom).

556 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

557 optional, intent(out) :: u_cor

558 !< The zonal velocities (u with a barotropic correction)

559 !! that give uhbt as the depth-integrated transport [L T-1 ~> m s-1]

560 type(bt_cont_type), optional, pointer :: bt_cont !< A structure with elements that describe the

561 !! effective open face areas as a function of barotropic flow.

562 real, dimension(SZIB_(G),SZJ_(G)), &

563 optional, intent(out) :: du_cor !< The zonal velocity increments from u that give uhbt

564 !! as the depth-integrated transports [L T-1 ~> m s-1].

565

566 ! Local variables

567 real, dimension(SZIB_(G),SZK_(GV)) :: duhdu ! Partial derivative of uh with u [H L ~> m2 or kg m-1].

568 real, dimension(SZIB_(G)) :: &

569 du, & ! Corrective barotropic change in the velocity to give uhbt [L T-1 ~> m s-1].

570 du_min_cfl, & ! Lower limit on du correction to avoid CFL violations [L T-1 ~> m s-1]

571 du_max_cfl, & ! Upper limit on du correction to avoid CFL violations [L T-1 ~> m s-1]

572 duhdu_tot_0, & ! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].

573 uh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].

574 visc_rem_max ! The column maximum of visc_rem [nondim].

575 logical, dimension(SZIB_(G)) :: do_i

576 real, dimension(SZIB_(G),SZK_(GV)) :: &

577 visc_rem ! A 2-D copy of visc_rem_u or an array of 1's [nondim].

578 real, dimension(SZIB_(G)) :: faui ! A list of sums of zonal face areas [H L ~> m2 or kg m-1].

579 real :: fa_u ! A sum of zonal face areas [H L ~> m2 or kg m-1].

580 real :: i_vrm ! 1.0 / visc_rem_max [nondim]

581 real :: cfl_dt ! The maximum CFL ratio of the adjusted velocities divided by

582 ! the time step [T-1 ~> s-1].

583 real :: i_dt ! 1.0 / dt [T-1 ~> s-1].

584 real :: du_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].

585 real :: dx_e, dx_w ! Effective x-grid spacings to the east and west [L ~> m].

586 type(cont_loop_bounds_type) :: lb

587 integer :: i, j, k, ish, ieh, jsh, jeh, n, nz

588 integer :: l_seg ! The OBC segment number

589 logical :: use_visc_rem, set_bt_cont

590 logical :: local_specified_bc, local_flather_obc, local_open_bc, any_simple_obc ! OBC-related logicals

591 logical :: simple_obc_pt(szib_(g)) ! Indicates points in a row with specified transport OBCs

592

593 call cpu_clock_begin(id_clock_correct)

594

595 use_visc_rem = present(visc_rem_u)

596

597 set_bt_cont = .false. ; if (present(bt_cont)) set_bt_cont = (associated(bt_cont))

598

599 local_specified_bc = .false. ; local_flather_obc = .false. ; local_open_bc = .false.

600 if (associated(obc)) then ; if (obc%OBC_pe) then

601 local_specified_bc = obc%specified_u_BCs_exist_globally

602 local_flather_obc = obc%Flather_u_BCs_exist_globally

603 local_open_bc = obc%open_u_BCs_exist_globally

604 endif ; endif

605

606 if (present(du_cor)) du_cor(:,:) = 0.0

607

608 if (present(lb_in)) then

609 lb = lb_in

610 else

611 lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec

612 endif

613 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = gv%ke

614

615 cfl_dt = cs%CFL_limit_adjust / dt

616 i_dt = 1.0 / dt

617 if (cs%aggress_adjust) cfl_dt = i_dt

618

619 if (.not.use_visc_rem) visc_rem(:,:) = 1.0

620 !$OMP parallel do default(shared) private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0, &

621 !$OMP duhdu_tot_0,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W, &

622 !$OMP simple_OBC_pt,any_simple_OBC,l_seg) &

623 !$OMP firstprivate(visc_rem)

624 do j=jsh,jeh

625 do i=ish-1,ieh ; do_i(i) = .true. ; enddo

626 ! Set uh and duhdu.

627 do k=1,nz

628 if (use_visc_rem) then ; do i=ish-1,ieh

629 visc_rem(i,k) = visc_rem_u(i,j,k)

630 enddo ; endif

631 call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_w(:,j,k), h_e(:,j,k), &

632 uh(:,j,k), duhdu(:,k), visc_rem(:,k), &

633 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, por_face_areau(:,j,k), &

634 cs%h_marg_min, obc)

635 if (local_specified_bc) then

636 do i=ish-1,ieh ; if (obc%segnum_u(i,j) /= 0) then

637 l_seg = abs(obc%segnum_u(i,j))

638 if (obc%segment(l_seg)%specified) uh(i,j,k) = obc%segment(l_seg)%normal_trans(i,j,k)

639 endif ; enddo

640 endif

641 enddo

642

643 if (present(uhbt) .or. set_bt_cont) then

644 if (use_visc_rem .and. cs%use_visc_rem_max) then

645 visc_rem_max(:) = 0.0

646 do k=1,nz ; do i=ish-1,ieh

647 visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))

648 enddo ; enddo

649 else

650 visc_rem_max(:) = 1.0

651 endif

652 ! Set limits on du that will keep the CFL number between -1 and 1.

653 ! This should be adequate to keep the root bracketed in all cases.

654 do i=ish-1,ieh

655 i_vrm = 0.0

656 if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)

657 if (cs%vol_CFL) then

658 dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))

659 dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))

660 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

661 du_max_cfl(i) = 2.0* (cfl_dt * dx_w) * i_vrm

662 du_min_cfl(i) = -2.0 * (cfl_dt * dx_e) * i_vrm

663 uh_tot_0(i) = 0.0 ; duhdu_tot_0(i) = 0.0

664 enddo

665 do k=1,nz ; do i=ish-1,ieh

666 duhdu_tot_0(i) = duhdu_tot_0(i) + duhdu(i,k)

667 uh_tot_0(i) = uh_tot_0(i) + uh(i,j,k)

668 enddo ; enddo

669 if (use_visc_rem) then

670 if (cs%aggress_adjust) then

671 do k=1,nz ; do i=ish-1,ieh

672 if (cs%vol_CFL) then

673 dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))

674 dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))

675 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

676

677 du_lim = 0.499*((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k)))

678 if (du_max_cfl(i) * visc_rem(i,k) > du_lim) &

679 du_max_cfl(i) = du_lim / visc_rem(i,k)

680

681 du_lim = 0.499*((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k)))

682 if (du_min_cfl(i) * visc_rem(i,k) < du_lim) &

683 du_min_cfl(i) = du_lim / visc_rem(i,k)

684 enddo ; enddo

685 else

686 do k=1,nz ; do i=ish-1,ieh

687 if (cs%vol_CFL) then

688 dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))

689 dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))

690 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

691

692 if (du_max_cfl(i) * visc_rem(i,k) > dx_w*cfl_dt - u(i,j,k)*g%mask2dCu(i,j)) &

693 du_max_cfl(i) = (dx_w*cfl_dt - u(i,j,k)) / visc_rem(i,k)

694 if (du_min_cfl(i) * visc_rem(i,k) < -dx_e*cfl_dt - u(i,j,k)*g%mask2dCu(i,j)) &

695 du_min_cfl(i) = -(dx_e*cfl_dt + u(i,j,k)) / visc_rem(i,k)

696 enddo ; enddo

697 endif

698 else

699 if (cs%aggress_adjust) then

700 do k=1,nz ; do i=ish-1,ieh

701 if (cs%vol_CFL) then

702 dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))

703 dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))

704 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

705

706 du_max_cfl(i) = min(du_max_cfl(i), 0.499 * &

707 ((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k))) )

708 du_min_cfl(i) = max(du_min_cfl(i), 0.499 * &

709 ((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k))) )

710 enddo ; enddo

711 else

712 do k=1,nz ; do i=ish-1,ieh

713 if (cs%vol_CFL) then

714 dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))

715 dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))

716 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

717

718 du_max_cfl(i) = min(du_max_cfl(i), dx_w*cfl_dt - u(i,j,k))

719 du_min_cfl(i) = max(du_min_cfl(i), -(dx_e*cfl_dt + u(i,j,k)))

720 enddo ; enddo

721 endif

722 endif

723 do i=ish-1,ieh

724 du_max_cfl(i) = max(du_max_cfl(i),0.0)

725 du_min_cfl(i) = min(du_min_cfl(i),0.0)

726 enddo

727

728 any_simple_obc = .false.

729 if (present(uhbt) .or. set_bt_cont) then

730 if (local_specified_bc .or. local_flather_obc) then ; do i=ish-1,ieh

731 l_seg = abs(obc%segnum_u(i,j))

732

733 ! Avoid reconciling barotropic/baroclinic transports if transport is specified

734 simple_obc_pt(i) = .false.

735 if (l_seg /= obc_none) simple_obc_pt(i) = obc%segment(l_seg)%specified

736 do_i(i) = .not.simple_obc_pt(i)

737 any_simple_obc = any_simple_obc .or. simple_obc_pt(i)

738 enddo ; else ; do i=ish-1,ieh

739 do_i(i) = .true.

740 enddo ; endif

741 endif

742

743 if (present(uhbt)) then

744 ! Find du and uh.

745 call zonal_flux_adjust(u, h_in, h_w, h_e, uhbt(:,j), uh_tot_0, duhdu_tot_0, du, &

746 du_max_cfl, du_min_cfl, dt, g, gv, us, cs, visc_rem, &

747 j, ish, ieh, do_i, por_face_areau, uh, obc=obc)

748

749 if (present(u_cor)) then ; do k=1,nz

750 do i=ish-1,ieh ; u_cor(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo

751 if (any_simple_obc) then ; do i=ish-1,ieh ; if (simple_obc_pt(i)) then

752 u_cor(i,j,k) = obc%segment(abs(obc%segnum_u(i,j)))%normal_vel(i,j,k)

753 endif ; enddo ; endif

754 enddo ; endif ! u-corrected

755

756 if (present(du_cor)) then

757 do i=ish-1,ieh ; du_cor(i,j) = du(i) ; enddo

758 endif

759

760 endif

761

762 if (set_bt_cont) then

763 call set_zonal_bt_cont(u, h_in, h_w, h_e, bt_cont, uh_tot_0, duhdu_tot_0,&

764 du_max_cfl, du_min_cfl, dt, g, gv, us, cs, visc_rem, &

765 visc_rem_max, j, ish, ieh, do_i, por_face_areau)

766 if (any_simple_obc) then

767 do i=ish-1,ieh

768 if (simple_obc_pt(i)) faui(i) = gv%H_subroundoff*g%dy_Cu(i,j)

769 enddo

770 ! NOTE: simple_OBC_pt(I) should prevent access to segment OBC_NONE

771 do k=1,nz ; do i=ish-1,ieh ; if (simple_obc_pt(i)) then

772 l_seg = abs(obc%segnum_u(i,j))

773 if ((abs(obc%segment(l_seg)%normal_vel(i,j,k)) > 0.0) .and. (obc%segment(l_seg)%specified)) &

774 faui(i) = faui(i) + obc%segment(l_seg)%normal_trans(i,j,k) / obc%segment(l_seg)%normal_vel(i,j,k)

775 endif ; enddo ; enddo

776 do i=ish-1,ieh ; if (simple_obc_pt(i)) then

777 bt_cont%FA_u_W0(i,j) = faui(i) ; bt_cont%FA_u_E0(i,j) = faui(i)

778 bt_cont%FA_u_WW(i,j) = faui(i) ; bt_cont%FA_u_EE(i,j) = faui(i)

779 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0

780 endif ; enddo

781 endif

782 endif ! set_BT_cont

783

784 endif ! present(uhbt) or set_BT_cont

785

786 enddo ! j-loop

787

788 if (local_open_bc .and. set_bt_cont) then

789 do n = 1, obc%number_of_segments

790 if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then

791 i = obc%segment(n)%HI%IsdB

792 if (obc%segment(n)%direction == obc_direction_e) then

793 do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed

794 fa_u = 0.0

795 do k=1,nz ; fa_u = fa_u + h_in(i,j,k)*(g%dy_Cu(i,j)*por_face_areau(i,j,k)) ; enddo

796 bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u

797 bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u

798 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0

799 enddo

800 else

801 do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed

802 fa_u = 0.0

803 do k=1,nz ; fa_u = fa_u + h_in(i+1,j,k)*(g%dy_Cu(i,j)*por_face_areau(i,j,k)) ; enddo

804 bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u

805 bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u

806 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0

807 enddo

808 endif

809 endif

810 enddo

811 endif

812

813 if (set_bt_cont) then ; if (allocated(bt_cont%h_u)) then

814 if (present(u_cor)) then

815 call zonal_flux_thickness(u_cor, h_in, h_w, h_e, bt_cont%h_u, dt, g, gv, us, lb, &

816 cs%vol_CFL, cs%marginal_faces, obc, por_face_areau, visc_rem_u)

817 else

818 call zonal_flux_thickness(u, h_in, h_w, h_e, bt_cont%h_u, dt, g, gv, us, lb, &

819 cs%vol_CFL, cs%marginal_faces, obc, por_face_areau, visc_rem_u)

820 endif

821 endif ; endif

822

823 call cpu_clock_end(id_clock_correct)

824

825end subroutine zonal_mass_flux

826

827

828!> Calculates the vertically integrated mass or volume fluxes through the zonal faces.

829subroutine zonal_bt_mass_flux(u, h_in, h_W, h_E, uhbt, dt, G, GV, US, CS, OBC, por_face_areaU, LB_in)

830 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

831 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

832 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]

833 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to

834 !! calculate fluxes [H ~> m or kg m-2]

835 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_w !< Western edge thickness in the PPM

836 !! reconstruction [H ~> m or kg m-2].

837 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_e !< Eastern edge thickness in the PPM

838 !! reconstruction [H ~> m or kg m-2].

839 real, dimension(SZIB_(G),SZJ_(G)), intent(out) :: uhbt !< The summed volume flux through zonal

840 !! faces [H L2 T-1 ~> m3 s-1 or kg s-1].

841 real, intent(in) :: dt !< Time increment [T ~> s].

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

843 type(continuity_ppm_cs), intent(in) :: cs !< This module's control structure.G

844 type(ocean_obc_type), pointer :: obc !< Open boundary condition type

845 !! specifies whether, where, and what

846 !! open boundary conditions are used.

847 real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: por_face_areau !< fractional open area of U-faces [nondim]

848 type(cont_loop_bounds_type), optional, intent(in) :: lb_in !< Loop bounds structure.

849

850 ! Local variables

851 real :: uh(szib_(g)) ! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]

852 real :: duhdu(szib_(g)) ! Partial derivative of uh with u [H L ~> m2 or kg m-1].

853 logical, dimension(SZIB_(G)) :: do_i

854 real :: ones(szib_(g)) ! An array of 1's [nondim]

855 integer :: i, j, k, ish, ieh, jsh, jeh, nz, l_seg

856 logical :: local_specified_bc, obc_in_row

857

858 call cpu_clock_begin(id_clock_correct)

859

860 local_specified_bc = .false.

861 if (associated(obc)) then ; if (obc%OBC_pe) then

862 local_specified_bc = obc%specified_v_BCs_exist_globally

863 endif ; endif

864

865 if (present(lb_in)) then

866 ish = lb_in%ish ; ieh = lb_in%ieh ; jsh = lb_in%jsh ; jeh = lb_in%jeh ; nz = gv%ke

867 else

868 ish = g%isc ; ieh = g%iec ; jsh = g%jsc ; jeh = g%jec ; nz = gv%ke

869 endif

870

871 ones(:) = 1.0 ; do_i(:) = .true.

872

873 uhbt(:,:) = 0.0

874 !$OMP parallel do default(shared) private(uh,duhdu,OBC_in_row)

875 do j=jsh,jeh

876 ! Determining whether there are any OBC points outside of the k-loop should be more efficient.

877 obc_in_row = .false.

878 if (local_specified_bc) then ; do i=ish-1,ieh ; if (obc%segnum_u(i,j) /= 0) then

879 if (obc%segment(abs(obc%segnum_u(i,j)))%specified) obc_in_row = .true.

880 endif ; enddo ; endif

881 do k=1,nz

882 ! This sets uh and duhdu.

883 call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_w(:,j,k), h_e(:,j,k), uh, duhdu, ones, &

884 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, por_face_areau(:,j,k), &

885 cs%h_marg_min, obc)

886 if (obc_in_row) then ; do i=ish-1,ieh ; if (obc%segnum_u(i,j) /= 0) then

887 l_seg = abs(obc%segnum_u(i,j))

888 if (obc%segment(l_seg)%specified) uh(i) = obc%segment(l_seg)%normal_trans(i,j,k)

889 endif ; enddo ; endif

890

891 ! Accumulate the barotropic transport.

892 do i=ish-1,ieh

893 uhbt(i,j) = uhbt(i,j) + uh(i)

894 enddo

895 enddo ! k-loop

896 enddo ! j-loop

897 call cpu_clock_end(id_clock_correct)

898

899end subroutine zonal_bt_mass_flux

900

901

902!> Evaluates the zonal mass or volume fluxes in a layer.

903subroutine zonal_flux_layer(u, h, h_W, h_E, uh, duhdu, visc_rem, dt, G, US, j, &

904 ish, ieh, do_I, vol_CFL, por_face_areaU, h_marg_min, OBC)

905 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

906 real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

907 real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the

908 !! momentum originally in a layer that remains after a time-step

909 !! of viscosity, and the fraction of a time-step's worth of a barotropic

910 !! acceleration that a layer experiences after viscosity is applied [nondim].

911 !! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).

912 real, dimension(SZI_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].

913 real, dimension(SZI_(G)), intent(in) :: h_W !< West edge thickness [H ~> m or kg m-2].

914 real, dimension(SZI_(G)), intent(in) :: h_E !< East edge thickness [H ~> m or kg m-2].

915 real, dimension(SZIB_(G)), intent(inout) :: uh !< Zonal mass or volume

916 !! transport [H L2 T-1 ~> m3 s-1 or kg s-1].

917 real, dimension(SZIB_(G)), intent(inout) :: duhdu !< Partial derivative of uh

918 !! with u [H L ~> m2 or kg m-1].

919 real, intent(in) :: dt !< Time increment [T ~> s]

920 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

921 integer, intent(in) :: j !< Spatial index.

922 integer, intent(in) :: ish !< Start of index range.

923 integer, intent(in) :: ieh !< End of index range.

924 logical, dimension(SZIB_(G)), intent(in) :: do_I !< Which i values to work on.

925 logical, intent(in) :: vol_CFL !< If true, rescale the

926 real, dimension(SZIB_(G)), intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]

927 !! ratio of face areas to the cell areas when estimating the CFL number.

928 real, intent(in) :: h_marg_min !< Negligible floor on h_marg [H ~> m or kg m-2]

929 type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.

930 ! Local variables

931 real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]

932 real :: curv_3 ! A measure of the thickness curvature over a grid length [H ~> m or kg m-2]

933 real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].

934 integer :: i

935 logical :: local_open_BC

936

937 local_open_bc = .false.

938 if (present(obc)) then ; if (associated(obc)) then

939 local_open_bc = obc%open_u_BCs_exist_globally

940 endif ; endif

941

942 do i=ish-1,ieh ; if (do_i(i)) then

943 ! Set new values of uh and duhdu.

944 if (u(i) > 0.0) then

945 if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))

946 else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif

947 curv_3 = (h_w(i) + h_e(i)) - 2.0*h(i)

948 uh(i) = (g%dy_Cu(i,j) * por_face_areau(i)) * u(i) * &

949 (h_e(i) + cfl * (0.5*(h_w(i) - h_e(i)) + curv_3*(cfl - 1.5)))

950 h_marg = h_e(i) + cfl * ((h_w(i) - h_e(i)) + 3.0*curv_3*(cfl - 1.0))

951 elseif (u(i) < 0.0) then

952 if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))

953 else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif

954 curv_3 = (h_w(i+1) + h_e(i+1)) - 2.0*h(i+1)

955 uh(i) = (g%dy_Cu(i,j) * por_face_areau(i)) * u(i) * &

956 (h_w(i+1) + cfl * (0.5*(h_e(i+1)-h_w(i+1)) + curv_3*(cfl - 1.5)))

957 h_marg = h_w(i+1) + cfl * ((h_e(i+1)-h_w(i+1)) + 3.0*curv_3*(cfl - 1.0))

958 else

959 uh(i) = 0.0

960 h_marg = 0.5 * (h_w(i+1) + h_e(i))

961 endif

962 h_marg = max(h_marg, h_marg_min)

963 duhdu(i) = (g%dy_Cu(i,j) * por_face_areau(i)) * h_marg * visc_rem(i)

964 endif ; enddo

965

966 if (local_open_bc) then

967 do i=ish-1,ieh ; if (do_i(i)) then ; if (obc%segnum_u(i,j) /= 0) then

968 if (obc%segment(abs(obc%segnum_u(i,j)))%open) then

969 if (obc%segnum_u(i,j) > 0) then ! OBC_DIRECTION_E

970 uh(i) = (g%dy_Cu(i,j) * por_face_areau(i)) * u(i) * h(i)

971 duhdu(i) = (g%dy_Cu(i,j) * por_face_areau(i)) * max(h(i), h_marg_min) * visc_rem(i)

972 else ! OBC_DIRECTION_W

973 uh(i) = (g%dy_Cu(i,j) * por_face_areau(i)) * u(i) * h(i+1)

974 duhdu(i) = (g%dy_Cu(i,j)* por_face_areau(i)) * max(h(i+1), h_marg_min) * visc_rem(i)

975 endif

976 endif

977 endif ; endif ; enddo

978 endif

979end subroutine zonal_flux_layer

980

981!> Sets the effective interface thickness associated with the fluxes at each zonal velocity point,

982!! optionally scaling back these thicknesses to account for viscosity and fractional open areas.

983subroutine zonal_flux_thickness(u, h, h_W, h_E, h_u, dt, G, GV, US, LB, vol_CFL, &

984 marginal, OBC, por_face_areaU, visc_rem_u)

985 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

986 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

987 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

988 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to

989 !! calculate fluxes [H ~> m or kg m-2].

990 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_w !< West edge thickness in the

991 !! reconstruction [H ~> m or kg m-2].

992 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_e !< East edge thickness in the

993 !! reconstruction [H ~> m or kg m-2].

994 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_u !< Effective thickness at zonal faces,

995 !! scaled down to account for the effects of

996 !! viscosity and the fractional open area

997 !! [H ~> m or kg m-2].

998 real, intent(in) :: dt !< Time increment [T ~> s].

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

1000 type(cont_loop_bounds_type), intent(in) :: lb !< Loop bounds structure.

1001 logical, intent(in) :: vol_cfl !< If true, rescale the ratio

1002 !! of face areas to the cell areas when estimating the CFL number.

1003 logical, intent(in) :: marginal !< If true, report the

1004 !! marginal face thicknesses; otherwise report transport-averaged thicknesses.

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

1006 intent(in) :: por_face_areau !< fractional open area of U-faces [nondim]

1007 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

1008 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1009 optional, intent(in) :: visc_rem_u

1010 !< Both the fraction of the momentum originally in a layer that remains after

1011 !! a time-step of viscosity, and the fraction of a time-step's worth of a

1012 !! barotropic acceleration that a layer experiences after viscosity is applied [nondim].

1013 !! Visc_rem_u is between 0 (at the bottom) and 1 (far above the bottom).

1014

1015 ! Local variables

1016 real :: cfl ! The CFL number based on the local velocity and grid spacing [nondim]

1017 real :: curv_3 ! A measure of the thickness curvature over a grid length [H ~> m or kg m-2]

1018 real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].

1019 real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].

1020 logical :: local_open_bc

1021 integer :: i, j, k, ish, ieh, jsh, jeh, nz, n

1022 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = gv%ke

1023

1024 !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)

1025 do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh

1026 if (u(i,j,k) > 0.0) then

1027 if (vol_cfl) then ; cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))

1028 else ; cfl = u(i,j,k) * dt * g%IdxT(i,j) ; endif

1029 curv_3 = (h_w(i,j,k) + h_e(i,j,k)) - 2.0*h(i,j,k)

1030 h_avg = h_e(i,j,k) + cfl * (0.5*(h_w(i,j,k) - h_e(i,j,k)) + curv_3*(cfl - 1.5))

1031 h_marg = h_e(i,j,k) + cfl * ((h_w(i,j,k) - h_e(i,j,k)) + 3.0*curv_3*(cfl - 1.0))

1032 elseif (u(i,j,k) < 0.0) then

1033 if (vol_cfl) then ; cfl = (-u(i,j,k)*dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))

1034 else ; cfl = -u(i,j,k) * dt * g%IdxT(i+1,j) ; endif

1035 curv_3 = (h_w(i+1,j,k) + h_e(i+1,j,k)) - 2.0*h(i+1,j,k)

1036 h_avg = h_w(i+1,j,k) + cfl * (0.5*(h_e(i+1,j,k)-h_w(i+1,j,k)) + curv_3*(cfl - 1.5))

1037 h_marg = h_w(i+1,j,k) + cfl * ((h_e(i+1,j,k)-h_w(i+1,j,k)) + &

1038 3.0*curv_3*(cfl - 1.0))

1039 else

1040 h_avg = 0.5 * (h_w(i+1,j,k) + h_e(i,j,k))

1041 ! The choice to use the arithmetic mean here is somewhat arbitrarily, but

1042 ! it should be noted that h_W(i+1,j,k) and h_E(i,j,k) are usually the same.

1043 h_marg = 0.5 * (h_w(i+1,j,k) + h_e(i,j,k))

1044 ! h_marg = (2.0 * h_W(i+1,j,k) * h_E(i,j,k)) / &

1045 ! (h_W(i+1,j,k) + h_E(i,j,k) + GV%H_subroundoff)

1046 endif

1047

1048 if (marginal) then ; h_u(i,j,k) = h_marg

1049 else ; h_u(i,j,k) = h_avg ; endif

1050 enddo ; enddo ; enddo

1051 if (present(visc_rem_u)) then

1052 ! Scale back the thickness to account for the effects of viscosity and the fractional open

1053 ! thickness to give an appropriate non-normalized weight for each layer in determining the

1054 ! barotropic acceleration.

1055 !$OMP parallel do default(shared)

1056 do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh

1057 h_u(i,j,k) = h_u(i,j,k) * (visc_rem_u(i,j,k) * por_face_areau(i,j,k))

1058 enddo ; enddo ; enddo

1059 else

1060 !$OMP parallel do default(shared)

1061 do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh

1062 h_u(i,j,k) = h_u(i,j,k) * por_face_areau(i,j,k)

1063 enddo ; enddo ; enddo

1064 endif

1065

1066 local_open_bc = .false.

1067 if (associated(obc)) local_open_bc = obc%open_u_BCs_exist_globally

1068 if (local_open_bc) then

1069 do n = 1, obc%number_of_segments

1070 if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then

1071 i = obc%segment(n)%HI%IsdB

1072 if (obc%segment(n)%direction == obc_direction_e) then

1073 if (present(visc_rem_u)) then ; do k=1,nz

1074 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed

1075 h_u(i,j,k) = h(i,j,k) * (visc_rem_u(i,j,k) * por_face_areau(i,j,k))

1076 enddo

1077 enddo ; else ; do k=1,nz

1078 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed

1079 h_u(i,j,k) = h(i,j,k) * por_face_areau(i,j,k)

1080 enddo

1081 enddo ; endif

1082 else

1083 if (present(visc_rem_u)) then ; do k=1,nz

1084 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed

1085 h_u(i,j,k) = h(i+1,j,k) * (visc_rem_u(i,j,k) * por_face_areau(i,j,k))

1086 enddo

1087 enddo ; else ; do k=1,nz

1088 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed

1089 h_u(i,j,k) = h(i+1,j,k) * por_face_areau(i,j,k)

1090 enddo

1091 enddo ; endif

1092 endif

1093 endif

1094 enddo

1095 endif

1096

1097end subroutine zonal_flux_thickness

1098

1099!> Returns the barotropic velocity adjustment that gives the

1100!! desired barotropic (layer-summed) transport.

1101subroutine zonal_flux_adjust(u, h_in, h_W, h_E, uhbt, uh_tot_0, duhdu_tot_0, &

1102 du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &

1103 j, ish, ieh, do_I_in, por_face_areaU, uh_3d, OBC)

1104

1105 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

1106 type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure.

1107 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

1108 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to

1109 !! calculate fluxes [H ~> m or kg m-2].

1110 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_W !< West edge thickness in the

1111 !! reconstruction [H ~> m or kg m-2].

1112 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_E !< East edge thickness in the

1113 !! reconstruction [H ~> m or kg m-2].

1114 real, dimension(SZIB_(G),SZK_(GV)), intent(in) :: visc_rem !< Both the fraction of the

1115 !! momentum originally in a layer that remains after a time-step of viscosity, and

1116 !! the fraction of a time-step's worth of a barotropic acceleration that a layer

1117 !! experiences after viscosity is applied [nondim].

1118 !! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).

1119 real, dimension(SZIB_(G)), intent(in) :: uhbt !< The summed volume flux

1120 !! through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].

1121

1122 real, dimension(SZIB_(G)), intent(in) :: du_max_CFL !< Maximum acceptable

1123 !! value of du [L T-1 ~> m s-1].

1124 real, dimension(SZIB_(G)), intent(in) :: du_min_CFL !< Minimum acceptable

1125 !! value of du [L T-1 ~> m s-1].

1126 real, dimension(SZIB_(G)), intent(in) :: uh_tot_0 !< The summed transport

1127 !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].

1128 real, dimension(SZIB_(G)), intent(in) :: duhdu_tot_0 !< The partial derivative

1129 !! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].

1130 real, dimension(SZIB_(G)), intent(out) :: du !<

1131 !! The barotropic velocity adjustment [L T-1 ~> m s-1].

1132 real, intent(in) :: dt !< Time increment [T ~> s].

1133 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1134 type(continuity_ppm_cs), intent(in) :: CS !< This module's control structure.

1135 integer, intent(in) :: j !< Spatial index.

1136 integer, intent(in) :: ish !< Start of index range.

1137 integer, intent(in) :: ieh !< End of index range.

1138 logical, dimension(SZIB_(G)), intent(in) :: do_I_in !<

1139 !! A logical flag indicating which I values to work on.

1140 real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &

1141 intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]

1142 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), optional, intent(inout) :: uh_3d !<

1143 !! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].

1144 type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.

1145 ! Local variables

1146 real, dimension(SZIB_(G),SZK_(GV)) :: &

1147 uh_aux, & ! An auxiliary zonal volume flux [H L2 T-1 ~> m3 s-1 or kg s-1].

1148 duhdu ! Partial derivative of uh with u [H L ~> m2 or kg m-1].

1149 real, dimension(SZIB_(G)) :: &

1150 uh_err, & ! Difference between uhbt and the summed uh [H L2 T-1 ~> m3 s-1 or kg s-1].

1151 uh_err_best, & ! The smallest value of uh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].

1152 u_new, & ! The velocity with the correction added [L T-1 ~> m s-1].

1153 duhdu_tot,&! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].

1154 du_min, & ! Lower limit on du correction based on CFL limits and previous iterations [L T-1 ~> m s-1]

1155 du_max ! Upper limit on du correction based on CFL limits and previous iterations [L T-1 ~> m s-1]

1156 real :: du_prev ! The previous value of du [L T-1 ~> m s-1].

1157 real :: ddu ! The change in du from the previous iteration [L T-1 ~> m s-1].

1158 real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].

1159 real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].

1160 integer :: i, k, nz, itt, max_itts = 20

1161 logical :: domore, do_I(SZIB_(G))

1162

1163 nz = gv%ke

1164

1165 uh_aux(:,:) = 0.0 ; duhdu(:,:) = 0.0

1166

1167 if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh

1168 uh_aux(i,k) = uh_3d(i,j,k)

1169 enddo ; enddo ; endif

1170

1171 do i=ish-1,ieh

1172 du(i) = 0.0 ; do_i(i) = do_i_in(i)

1173 du_max(i) = du_max_cfl(i) ; du_min(i) = du_min_cfl(i)

1174 uh_err(i) = uh_tot_0(i) - uhbt(i) ; duhdu_tot(i) = duhdu_tot_0(i)

1175 uh_err_best(i) = abs(uh_err(i))

1176 enddo

1177

1178 do itt=1,max_itts

1179 select case (itt)

1180 case (:1) ; tol_eta = 1e-6 * cs%tol_eta

1181 case (2) ; tol_eta = 1e-4 * cs%tol_eta

1182 case (3) ; tol_eta = 1e-2 * cs%tol_eta

1183 case default ; tol_eta = cs%tol_eta

1184 end select

1185 tol_vel = cs%tol_vel

1186

1187 do i=ish-1,ieh

1188 if (uh_err(i) > 0.0) then ; du_max(i) = du(i)

1189 elseif (uh_err(i) < 0.0) then ; du_min(i) = du(i)

1190 else ; do_i(i) = .false. ; endif

1191 enddo

1192 domore = .false.

1193 do i=ish-1,ieh ; if (do_i(i)) then

1194 if ((dt * min(g%IareaT(i,j),g%IareaT(i+1,j))*abs(uh_err(i)) > tol_eta) .or. &

1195 (cs%better_iter .and. ((abs(uh_err(i)) > tol_vel * duhdu_tot(i)) .or. &

1196 (abs(uh_err(i)) > uh_err_best(i))) )) then

1197 ! Use Newton's method, provided it stays bounded. Otherwise bisect

1198 ! the value with the appropriate bound.

1199 ddu = -uh_err(i) / duhdu_tot(i)

1200 du_prev = du(i)

1201 du(i) = du(i) + ddu

1202 if (abs(ddu) < 1.0e-15*abs(du(i))) then

1203 do_i(i) = .false. ! ddu is small enough to quit.

1204 elseif (ddu > 0.0) then

1205 if (du(i) >= du_max(i)) then

1206 du(i) = 0.5*(du_prev + du_max(i))

1207 if (du_max(i) - du_prev < 1.0e-15*abs(du(i))) do_i(i) = .false.

1208 endif

1209 else ! ddu < 0.0

1210 if (du(i) <= du_min(i)) then

1211 du(i) = 0.5*(du_prev + du_min(i))

1212 if (du_prev - du_min(i) < 1.0e-15*abs(du(i))) do_i(i) = .false.

1213 endif

1214 endif

1215 if (do_i(i)) domore = .true.

1216 else

1217 do_i(i) = .false.

1218 endif

1219 endif ; enddo

1220 if (.not.domore) exit

1221

1222 if ((itt < max_itts) .or. present(uh_3d)) then ; do k=1,nz

1223 do i=ish-1,ieh ; u_new(i) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo

1224 call zonal_flux_layer(u_new, h_in(:,j,k), h_w(:,j,k), h_e(:,j,k), &

1225 uh_aux(:,k), duhdu(:,k), visc_rem(:,k), &

1226 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, por_face_areau(:,j,k), &

1227 cs%h_marg_min, obc)

1228 enddo ; endif

1229

1230 if (itt < max_itts) then

1231 do i=ish-1,ieh

1232 uh_err(i) = -uhbt(i) ; duhdu_tot(i) = 0.0

1233 enddo

1234 do k=1,nz ; do i=ish-1,ieh

1235 uh_err(i) = uh_err(i) + uh_aux(i,k)

1236 duhdu_tot(i) = duhdu_tot(i) + duhdu(i,k)

1237 enddo ; enddo

1238 do i=ish-1,ieh

1239 uh_err_best(i) = min(uh_err_best(i), abs(uh_err(i)))

1240 enddo

1241 endif

1242 enddo ! itt-loop

1243 ! If there are any faces which have not converged to within the tolerance,

1244 ! so-be-it, or else use a final upwind correction?

1245 ! This never seems to happen with 20 iterations as max_itt.

1246

1247 if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh

1248 uh_3d(i,j,k) = uh_aux(i,k)

1249 enddo ; enddo ; endif

1250

1251end subroutine zonal_flux_adjust

1252

1253!> Sets a structure that describes the zonal barotropic volume or mass fluxes as a

1254!! function of barotropic flow to agree closely with the sum of the layer's transports.

1255subroutine set_zonal_bt_cont(u, h_in, h_W, h_E, BT_cont, uh_tot_0, duhdu_tot_0, &

1256 du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &

1257 visc_rem_max, j, ish, ieh, do_I, por_face_areaU)

1258 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

1259 type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure.

1260 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].

1261 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to

1262 !! calculate fluxes [H ~> m or kg m-2].

1263 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_W !< West edge thickness in the

1264 !! reconstruction [H ~> m or kg m-2].

1265 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_E !< East edge thickness in the

1266 !! reconstruction [H ~> m or kg m-2].

1267 type(bt_cont_type), intent(inout) :: BT_cont !< A structure with elements

1268 !! that describe the effective open face areas as a function of barotropic flow.

1269 real, dimension(SZIB_(G)), intent(in) :: uh_tot_0 !< The summed transport

1270 !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].

1271 real, dimension(SZIB_(G)), intent(in) :: duhdu_tot_0 !< The partial derivative

1272 !! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].

1273 real, dimension(SZIB_(G)), intent(in) :: du_max_CFL !< Maximum acceptable

1274 !! value of du [L T-1 ~> m s-1].

1275 real, dimension(SZIB_(G)), intent(in) :: du_min_CFL !< Minimum acceptable

1276 !! value of du [L T-1 ~> m s-1].

1277 real, intent(in) :: dt !< Time increment [T ~> s].

1278 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1279 type(continuity_ppm_cs), intent(in) :: CS !< This module's control structure.

1280 real, dimension(SZIB_(G),SZK_(GV)), intent(in) :: visc_rem !< Both the fraction of the

1281 !! momentum originally in a layer that remains after a time-step of viscosity, and

1282 !! the fraction of a time-step's worth of a barotropic acceleration that a layer

1283 !! experiences after viscosity is applied [nondim].

1284 !! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).

1285 real, dimension(SZIB_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem [nondim].

1286 integer, intent(in) :: j !< Spatial index.

1287 integer, intent(in) :: ish !< Start of index range.

1288 integer, intent(in) :: ieh !< End of index range.

1289 logical, dimension(SZIB_(G)), intent(in) :: do_I !< A logical flag indicating

1290 !! which I values to work on.

1291 real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &

1292 intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]

1293 ! Local variables

1294 real, dimension(SZIB_(G)) :: &

1295 du0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].

1296 duL, duR, & ! The barotropic velocity increments that give the westerly

1297 ! (duL) and easterly (duR) test velocities [L T-1 ~> m s-1].

1298 zeros, & ! An array of full of 0 transports [H L2 T-1 ~> m3 s-1 or kg s-1]

1299 du_cfl, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].

1300 u_l, u_r, & ! The westerly (u_L), easterly (u_R), and zero-barotropic

1301 u_0, & ! transport (u_0) layer test velocities [L T-1 ~> m s-1].

1302 duhdu_l, & ! The effective layer marginal face areas with the westerly

1303 duhdu_r, & ! (_L), easterly (_R), and zero-barotropic (_0) test

1304 duhdu_0, & ! velocities [H L ~> m2 or kg m-1].

1305 uh_l, uh_r, & ! The layer transports with the westerly (_L), easterly (_R),

1306 uh_0, & ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].

1307 famt_l, famt_r, & ! The summed effective marginal face areas for the 3

1308 famt_0, & ! test velocities [H L ~> m2 or kg m-1].

1309 uhtot_l, & ! The summed transport with the westerly (uhtot_L) and

1310 uhtot_r ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].

1311 real :: FA_0 ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m-1].

1312 real :: FA_avg ! The average effective face area [L H ~> m2 or kg m-1], nominally given by

1313 ! the realized transport divided by the barotropic velocity.

1314 real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]. This

1315 ! limiting is necessary to keep the inverse of visc_rem

1316 ! from leading to large CFL numbers.

1317 real :: min_visc_rem ! The smallest permitted value for visc_rem that is used

1318 ! in finding the barotropic velocity that changes the

1319 ! flow direction [nondim]. This is necessary to keep the inverse

1320 ! of visc_rem from leading to large CFL numbers.

1321 real :: CFL_min ! A minimal increment in the CFL to try to ensure that the

1322 ! flow is truly upwind [nondim]

1323 real :: Idt ! The inverse of the time step [T-1 ~> s-1].

1324 logical :: domore

1325 integer :: i, k, nz

1326

1327 nz = gv%ke ; idt = 1.0 / dt

1328 min_visc_rem = 0.1 ; cfl_min = 1e-6

1329

1330 ! Diagnose the zero-transport correction, du0.

1331 do i=ish-1,ieh ; zeros(i) = 0.0 ; enddo

1332 call zonal_flux_adjust(u, h_in, h_w, h_e, zeros, uh_tot_0, duhdu_tot_0, du0, &

1333 du_max_cfl, du_min_cfl, dt, g, gv, us, cs, visc_rem, &

1334 j, ish, ieh, do_i, por_face_areau)

1335

1336 ! Determine the westerly- and easterly- fluxes. Choose a sufficiently

1337 ! negative velocity correction for the easterly-flux, and a sufficiently

1338 ! positive correction for the westerly-flux.

1339 domore = .false.

1340 do i=ish-1,ieh

1341 if (do_i(i)) domore = .true.

1342 du_cfl(i) = (cfl_min * idt) * g%dxCu(i,j)

1343 dur(i) = min(0.0,du0(i) - du_cfl(i))

1344 dul(i) = max(0.0,du0(i) + du_cfl(i))

1345 famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0

1346 uhtot_l(i) = 0.0 ; uhtot_r(i) = 0.0

1347 enddo

1348

1349 if (.not.domore) then

1350 do k=1,nz ; do i=ish-1,ieh

1351 bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0

1352 bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0

1353 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0

1354 enddo ; enddo

1355 return

1356 endif

1357

1358 do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then

1359 visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))

1360 if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.

1361 if (u(i,j,k) + dur(i)*visc_rem_lim > -du_cfl(i)*visc_rem(i,k)) &

1362 dur(i) = -(u(i,j,k) + du_cfl(i)*visc_rem(i,k)) / visc_rem_lim

1363 if (u(i,j,k) + dul(i)*visc_rem_lim < du_cfl(i)*visc_rem(i,k)) &

1364 dul(i) = -(u(i,j,k) - du_cfl(i)*visc_rem(i,k)) / visc_rem_lim

1365 endif

1366 endif ; enddo ; enddo

1367

1368 do k=1,nz

1369 do i=ish-1,ieh ; if (do_i(i)) then

1370 u_l(i) = u(i,j,k) + dul(i) * visc_rem(i,k)

1371 u_r(i) = u(i,j,k) + dur(i) * visc_rem(i,k)

1372 u_0(i) = u(i,j,k) + du0(i) * visc_rem(i,k)

1373 endif ; enddo

1374 call zonal_flux_layer(u_0, h_in(:,j,k), h_w(:,j,k), h_e(:,j,k), uh_0, duhdu_0, &

1375 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, &

1376 por_face_areau(:,j,k), cs%h_marg_min)

1377 call zonal_flux_layer(u_l, h_in(:,j,k), h_w(:,j,k), h_e(:,j,k), uh_l, duhdu_l, &

1378 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, &

1379 por_face_areau(:,j,k), cs%h_marg_min)

1380 call zonal_flux_layer(u_r, h_in(:,j,k), h_w(:,j,k), h_e(:,j,k), uh_r, duhdu_r, &

1381 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, &

1382 por_face_areau(:,j,k), cs%h_marg_min)

1383 do i=ish-1,ieh ; if (do_i(i)) then

1384 famt_0(i) = famt_0(i) + duhdu_0(i)

1385 famt_l(i) = famt_l(i) + duhdu_l(i)

1386 famt_r(i) = famt_r(i) + duhdu_r(i)

1387 uhtot_l(i) = uhtot_l(i) + uh_l(i)

1388 uhtot_r(i) = uhtot_r(i) + uh_r(i)

1389 endif ; enddo

1390 enddo

1391 do i=ish-1,ieh ; if (do_i(i)) then

1392 fa_0 = famt_0(i) ; fa_avg = famt_0(i)

1393 if ((dul(i) - du0(i)) /= 0.0) &

1394 fa_avg = uhtot_l(i) / (dul(i) - du0(i))

1395 if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))

1396 elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif

1397

1398 bt_cont%FA_u_W0(i,j) = fa_0 ; bt_cont%FA_u_WW(i,j) = famt_l(i)

1399 if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%uBT_WW(i,j) = 0.0 ; else

1400 bt_cont%uBT_WW(i,j) = (1.5 * (dul(i) - du0(i))) * &

1401 ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))

1402 endif

1403

1404 fa_0 = famt_0(i) ; fa_avg = famt_0(i)

1405 if ((dur(i) - du0(i)) /= 0.0) &

1406 fa_avg = uhtot_r(i) / (dur(i) - du0(i))

1407 if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))

1408 elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif

1409

1410 bt_cont%FA_u_E0(i,j) = fa_0 ; bt_cont%FA_u_EE(i,j) = famt_r(i)

1411 if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%uBT_EE(i,j) = 0.0 ; else

1412 bt_cont%uBT_EE(i,j) = (1.5 * (dur(i) - du0(i))) * &

1413 ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))

1414 endif

1415 else

1416 bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0

1417 bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0

1418 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0

1419 endif ; enddo

1420

1421end subroutine set_zonal_bt_cont

1422

1423!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.

1424subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, por_face_areaV, &

1425 LB_in, vhbt, visc_rem_v, v_cor, BT_cont, dv_cor)

1426 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

1427 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

1428 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]

1429 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to

1430 !! calculate fluxes [H ~> m or kg m-2]

1431 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_s !< South edge thickness in the

1432 !! reconstruction [H ~> m or kg m-2].

1433 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_n !< North edge thickness in the

1434 !! reconstruction [H ~> m or kg m-2].

1435 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: vh !< Volume flux through meridional

1436 !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]

1437 real, intent(in) :: dt !< Time increment [T ~> s].

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

1439 type(continuity_ppm_cs), intent(in) :: cs !< This module's control structure.G

1440 type(ocean_obc_type), pointer :: obc !< Open boundary condition type

1441 !! specifies whether, where, and what

1442 !! open boundary conditions are used.

1443 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: por_face_areav !< fractional open area of V-faces [nondim]

1444 type(cont_loop_bounds_type), optional, intent(in) :: lb_in !< Loop bounds structure.

1445 real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt !< The summed volume flux through meridional

1446 !! faces [H L2 T-1 ~> m3 s-1 or kg s-1].

1447 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1448 optional, intent(in) :: visc_rem_v !< Both the fraction of the momentum

1449 !! originally in a layer that remains after a time-step of viscosity,

1450 !! and the fraction of a time-step's worth of a barotropic acceleration

1451 !! that a layer experiences after viscosity is applied [nondim].

1452 !! Visc_rem_v is between 0 (at the bottom) and 1 (far above the bottom).

1453 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1454 optional, intent(out) :: v_cor

1455 !< The meridional velocities (v with a barotropic correction)

1456 !! that give vhbt as the depth-integrated transport [L T-1 ~> m s-1].

1457 type(bt_cont_type), optional, pointer :: bt_cont !< A structure with elements that describe

1458 !! the effective open face areas as a function of barotropic flow.

1459 real, dimension(SZI_(G),SZJB_(G)), &

1460 optional, intent(out) :: dv_cor !< The meridional velocity increments from v

1461 !! that give vhbt as the depth-integrated

1462 !! transports [L T-1 ~> m s-1].

1463

1464 ! Local variables

1465 real, dimension(SZI_(G),SZK_(GV)) :: &

1466 dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].

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

1468 dv, & ! Corrective barotropic change in the velocity to give vhbt [L T-1 ~> m s-1].

1469 dv_min_cfl, & ! Lower limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]

1470 dv_max_cfl, & ! Upper limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]

1471 dvhdv_tot_0, & ! Summed partial derivative of vh with v [H L ~> m2 or kg m-1].

1472 vh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].

1473 visc_rem_max ! The column maximum of visc_rem [nondim]

1474 logical, dimension(SZI_(G)) :: do_i

1475 real, dimension(SZI_(G)) :: favi ! A list of sums of meridional face areas [H L ~> m2 or kg m-1].

1476 real :: fa_v ! A sum of meridional face areas [H L ~> m2 or kg m-1].

1477 real, dimension(SZI_(G),SZK_(GV)) :: &

1478 visc_rem ! A 2-D copy of visc_rem_v or an array of 1's [nondim]

1479 real :: i_vrm ! 1.0 / visc_rem_max [nondim]

1480 real :: cfl_dt ! The maximum CFL ratio of the adjusted velocities divided by

1481 ! the time step [T-1 ~> s-1].

1482 real :: i_dt ! 1.0 / dt [T-1 ~> s-1].

1483 real :: dv_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].

1484 real :: dy_n, dy_s ! Effective y-grid spacings to the north and south [L ~> m].

1485 type(cont_loop_bounds_type) :: lb

1486 integer :: i, j, k, ish, ieh, jsh, jeh, n, nz

1487 integer :: l_seg ! The OBC segment number

1488 logical :: use_visc_rem, set_bt_cont

1489 logical :: local_specified_bc, local_flather_obc, local_open_bc, any_simple_obc ! OBC-related logicals

1490 logical :: simple_obc_pt(szi_(g)) ! Indicates points in a row with specified transport OBCs

1491 type(obc_segment_type), pointer :: segment => null()

1492

1493 call cpu_clock_begin(id_clock_correct)

1494

1495 use_visc_rem = present(visc_rem_v)

1496

1497 set_bt_cont = .false. ; if (present(bt_cont)) set_bt_cont = (associated(bt_cont))

1498

1499 local_specified_bc = .false. ; local_flather_obc = .false. ; local_open_bc = .false.

1500 if (associated(obc)) then ; if (obc%OBC_pe) then

1501 local_specified_bc = obc%specified_v_BCs_exist_globally

1502 local_flather_obc = obc%Flather_v_BCs_exist_globally

1503 local_open_bc = obc%open_v_BCs_exist_globally

1504 endif ; endif

1505

1506 if (present(dv_cor)) dv_cor(:,:) = 0.0

1507

1508 if (present(lb_in)) then

1509 lb = lb_in

1510 else

1511 lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec

1512 endif

1513 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = gv%ke

1514

1515 cfl_dt = cs%CFL_limit_adjust / dt

1516 i_dt = 1.0 / dt

1517 if (cs%aggress_adjust) cfl_dt = i_dt

1518

1519 if (.not.use_visc_rem) visc_rem(:,:) = 1.0

1520 !$OMP parallel do default(shared) private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0, &

1521 !$OMP dvhdv_tot_0,FAvi,visc_rem_max,I_vrm,dv_lim,dy_N,dy_S, &

1522 !$OMP simple_OBC_pt,any_simple_OBC,l_seg) &

1523 !$OMP firstprivate(visc_rem)

1524 do j=jsh-1,jeh

1525 do i=ish,ieh ; do_i(i) = .true. ; enddo

1526 ! This sets vh and dvhdv.

1527 do k=1,nz

1528 if (use_visc_rem) then ; do i=ish,ieh

1529 visc_rem(i,k) = visc_rem_v(i,j,k)

1530 enddo ; endif

1531 call merid_flux_layer(v(:,j,k), h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), &

1532 vh(:,j,k), dvhdv(:,k), visc_rem(:,k), &

1533 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, por_face_areav(:,:,k), &

1534 cs%h_marg_min, obc)

1535 if (local_specified_bc) then

1536 do i=ish,ieh ; if (obc%segnum_v(i,j) /= 0) then

1537 l_seg = abs(obc%segnum_v(i,j))

1538 if (obc%segment(l_seg)%specified) vh(i,j,k) = obc%segment(l_seg)%normal_trans(i,j,k)

1539 endif ; enddo

1540 endif

1541 enddo ! k-loop

1542

1543 if (present(vhbt) .or. set_bt_cont) then

1544 if (use_visc_rem .and. cs%use_visc_rem_max) then

1545 visc_rem_max(:) = 0.0

1546 do k=1,nz ; do i=ish,ieh

1547 visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))

1548 enddo ; enddo

1549 else

1550 visc_rem_max(:) = 1.0

1551 endif

1552 ! Set limits on dv that will keep the CFL number between -1 and 1.

1553 ! This should be adequate to keep the root bracketed in all cases.

1554 do i=ish,ieh

1555 i_vrm = 0.0

1556 if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)

1557 if (cs%vol_CFL) then

1558 dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))

1559 dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))

1560 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif

1561 dv_max_cfl(i) = 2.0 * (cfl_dt * dy_s) * i_vrm

1562 dv_min_cfl(i) = -2.0 * (cfl_dt * dy_n) * i_vrm

1563 vh_tot_0(i) = 0.0 ; dvhdv_tot_0(i) = 0.0

1564 enddo

1565 do k=1,nz ; do i=ish,ieh

1566 dvhdv_tot_0(i) = dvhdv_tot_0(i) + dvhdv(i,k)

1567 vh_tot_0(i) = vh_tot_0(i) + vh(i,j,k)

1568 enddo ; enddo

1569

1570 if (use_visc_rem) then

1571 if (cs%aggress_adjust) then

1572 do k=1,nz ; do i=ish,ieh

1573 if (cs%vol_CFL) then

1574 dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))

1575 dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))

1576 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif

1577 dv_lim = 0.499*((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k)))

1578 if (dv_max_cfl(i) * visc_rem(i,k) > dv_lim) &

1579 dv_max_cfl(i) = dv_lim / visc_rem(i,k)

1580

1581 dv_lim = 0.499*((-dy_n*cfl_dt - v(i,j,k)) + max(0.0,v(i,j+1,k)))

1582 if (dv_min_cfl(i) * visc_rem(i,k) < dv_lim) &

1583 dv_min_cfl(i) = dv_lim / visc_rem(i,k)

1584 enddo ; enddo

1585 else

1586 do k=1,nz ; do i=ish,ieh

1587 if (cs%vol_CFL) then

1588 dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))

1589 dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))

1590 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif

1591 if (dv_max_cfl(i) * visc_rem(i,k) > dy_s*cfl_dt - v(i,j,k)*g%mask2dCv(i,j)) &

1592 dv_max_cfl(i) = (dy_s*cfl_dt - v(i,j,k)) / visc_rem(i,k)

1593 if (dv_min_cfl(i) * visc_rem(i,k) < -dy_n*cfl_dt - v(i,j,k)*g%mask2dCv(i,j)) &

1594 dv_min_cfl(i) = -(dy_n*cfl_dt + v(i,j,k)) / visc_rem(i,k)

1595 enddo ; enddo

1596 endif

1597 else

1598 if (cs%aggress_adjust) then

1599 do k=1,nz ; do i=ish,ieh

1600 if (cs%vol_CFL) then

1601 dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))

1602 dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))

1603 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif

1604 dv_max_cfl(i) = min(dv_max_cfl(i), 0.499 * &

1605 ((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k))) )

1606 dv_min_cfl(i) = max(dv_min_cfl(i), 0.499 * &

1607 ((-dy_n*i_dt - v(i,j,k)) + max(0.0,v(i,j+1,k))) )

1608 enddo ; enddo

1609 else

1610 do k=1,nz ; do i=ish,ieh

1611 if (cs%vol_CFL) then

1612 dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))

1613 dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))

1614 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif

1615 dv_max_cfl(i) = min(dv_max_cfl(i), dy_s*cfl_dt - v(i,j,k))

1616 dv_min_cfl(i) = max(dv_min_cfl(i), -(dy_n*cfl_dt + v(i,j,k)))

1617 enddo ; enddo

1618 endif

1619 endif

1620 do i=ish,ieh

1621 dv_max_cfl(i) = max(dv_max_cfl(i),0.0)

1622 dv_min_cfl(i) = min(dv_min_cfl(i),0.0)

1623 enddo

1624

1625 any_simple_obc = .false.

1626 if (present(vhbt) .or. set_bt_cont) then

1627 if (local_specified_bc .or. local_flather_obc) then ; do i=ish,ieh

1628 l_seg = abs(obc%segnum_v(i,j))

1629

1630 ! Avoid reconciling barotropic/baroclinic transports if transport is specified

1631 simple_obc_pt(i) = .false.

1632 if (l_seg /= 0) simple_obc_pt(i) = obc%segment(l_seg)%specified

1633 do_i(i) = .not.simple_obc_pt(i)

1634 any_simple_obc = any_simple_obc .or. simple_obc_pt(i)

1635 enddo ; else ; do i=ish,ieh

1636 do_i(i) = .true.

1637 enddo ; endif

1638 endif

1639

1640 if (present(vhbt)) then

1641 ! Find dv and vh.

1642 call meridional_flux_adjust(v, h_in, h_s, h_n, vhbt(:,j), vh_tot_0, dvhdv_tot_0, dv, &

1643 dv_max_cfl, dv_min_cfl, dt, g, gv, us, cs, visc_rem, &

1644 j, ish, ieh, do_i, por_face_areav, vh, obc=obc)

1645

1646 if (present(v_cor)) then ; do k=1,nz

1647 do i=ish,ieh ; v_cor(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo

1648 if (any_simple_obc) then ; do i=ish,ieh ; if (simple_obc_pt(i)) then

1649 v_cor(i,j,k) = obc%segment(abs(obc%segnum_v(i,j)))%normal_vel(i,j,k)

1650 endif ; enddo ; endif

1651 enddo ; endif ! v-corrected

1652

1653 if (present(dv_cor)) then

1654 do i=ish,ieh ; dv_cor(i,j) = dv(i) ; enddo

1655 endif

1656

1657 endif

1658

1659 if (set_bt_cont) then

1660 call set_merid_bt_cont(v, h_in, h_s, h_n, bt_cont, vh_tot_0, dvhdv_tot_0,&

1661 dv_max_cfl, dv_min_cfl, dt, g, gv, us, cs, visc_rem, &

1662 visc_rem_max, j, ish, ieh, do_i, por_face_areav)

1663 if (any_simple_obc) then

1664 do i=ish,ieh

1665 if (simple_obc_pt(i)) favi(i) = gv%H_subroundoff*g%dx_Cv(i,j)

1666 enddo

1667 ! NOTE: simple_OBC_pt(i) should prevent access to segment OBC_NONE

1668 do k=1,nz ; do i=ish,ieh ; if (simple_obc_pt(i)) then

1669 segment => obc%segment(abs(obc%segnum_v(i,j)))

1670 if ((abs(segment%normal_vel(i,j,k)) > 0.0) .and. (segment%specified)) &

1671 favi(i) = favi(i) + segment%normal_trans(i,j,k) / segment%normal_vel(i,j,k)

1672 endif ; enddo ; enddo

1673 do i=ish,ieh ; if (simple_obc_pt(i)) then

1674 bt_cont%FA_v_S0(i,j) = favi(i) ; bt_cont%FA_v_N0(i,j) = favi(i)

1675 bt_cont%FA_v_SS(i,j) = favi(i) ; bt_cont%FA_v_NN(i,j) = favi(i)

1676 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0

1677 endif ; enddo

1678 endif

1679 endif ! set_BT_cont

1680

1681 endif ! present(vhbt) or set_BT_cont

1682

1683 enddo ! j-loop

1684

1685 if (local_open_bc .and. set_bt_cont) then

1686 do n = 1, obc%number_of_segments

1687 if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then

1688 j = obc%segment(n)%HI%JsdB

1689 if (obc%segment(n)%direction == obc_direction_n) then

1690 do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied

1691 fa_v = 0.0

1692 do k=1,nz ; fa_v = fa_v + h_in(i,j,k)*(g%dx_Cv(i,j)*por_face_areav(i,j,k)) ; enddo

1693 bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v

1694 bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v

1695 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0

1696 enddo

1697 else

1698 do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied

1699 fa_v = 0.0

1700 do k=1,nz ; fa_v = fa_v + h_in(i,j+1,k)*(g%dx_Cv(i,j)*por_face_areav(i,j,k)) ; enddo

1701 bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v

1702 bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v

1703 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0

1704 enddo

1705 endif

1706 endif

1707 enddo

1708 endif

1709

1710 if (set_bt_cont) then ; if (allocated(bt_cont%h_v)) then

1711 if (present(v_cor)) then

1712 call meridional_flux_thickness(v_cor, h_in, h_s, h_n, bt_cont%h_v, dt, g, gv, us, lb, &

1713 cs%vol_CFL, cs%marginal_faces, obc, por_face_areav, visc_rem_v)

1714 else

1715 call meridional_flux_thickness(v, h_in, h_s, h_n, bt_cont%h_v, dt, g, gv, us, lb, &

1716 cs%vol_CFL, cs%marginal_faces, obc, por_face_areav, visc_rem_v)

1717 endif

1718 endif ; endif

1719

1720 call cpu_clock_end(id_clock_correct)

1721

1722end subroutine meridional_mass_flux

1723

1724

1725!> Calculates the vertically integrated mass or volume fluxes through the meridional faces.

1726subroutine meridional_bt_mass_flux(v, h_in, h_S, h_N, vhbt, dt, G, GV, US, CS, OBC, por_face_areaV, LB_in)

1727 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

1728 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

1729 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]

1730 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to

1731 !! calculate fluxes [H ~> m or kg m-2]

1732 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_s !< Southern edge thickness in the PPM

1733 !! reconstruction [H ~> m or kg m-2].

1734 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_n !< Northern edge thickness in the PPM

1735 !! reconstruction [H ~> m or kg m-2].

1736 real, dimension(SZI_(G),SZJB_(G)), intent(out) :: vhbt !< The summed volume flux through meridional

1737 !! faces [H L2 T-1 ~> m3 s-1 or kg s-1].

1738 real, intent(in) :: dt !< Time increment [T ~> s].

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

1740 type(continuity_ppm_cs), intent(in) :: cs !< This module's control structure.G

1741 type(ocean_obc_type), pointer :: obc !< Open boundary condition type

1742 !! specifies whether, where, and what

1743 !! open boundary conditions are used.

1744 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: por_face_areav !< fractional open area of V-faces [nondim]

1745 type(cont_loop_bounds_type), optional, intent(in) :: lb_in !< Loop bounds structure.

1746

1747 ! Local variables

1748 real :: vh(szi_(g)) ! Volume flux through meridional faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]

1749 real :: dvhdv(szi_(g)) ! Partial derivative of vh with v [H L ~> m2 or kg m-1].

1750 logical, dimension(SZI_(G)) :: do_i

1751 real :: ones(szi_(g)) ! An array of 1's [nondim]

1752 integer :: i, j, k, ish, ieh, jsh, jeh, nz, l_seg

1753 logical :: local_specified_bc, obc_in_row

1754

1755 call cpu_clock_begin(id_clock_correct)

1756

1757 local_specified_bc = .false.

1758 if (associated(obc)) then ; if (obc%OBC_pe) then

1759 local_specified_bc = obc%specified_v_BCs_exist_globally

1760 endif ; endif

1761

1762 if (present(lb_in)) then

1763 ish = lb_in%ish ; ieh = lb_in%ieh ; jsh = lb_in%jsh ; jeh = lb_in%jeh ; nz = gv%ke

1764 else

1765 ish = g%isc ; ieh = g%iec ; jsh = g%jsc ; jeh = g%jec ; nz = gv%ke

1766 endif

1767

1768 ones(:) = 1.0 ; do_i(:) = .true.

1769

1770 vhbt(:,:) = 0.0

1771 !$OMP parallel do default(shared) private(vh,dvhdv,OBC_in_row)

1772 do j=jsh-1,jeh

1773 ! Determining whether there are any OBC points outside of the k-loop should be more efficient.

1774 obc_in_row = .false.

1775 if (local_specified_bc) then ; do i=ish,ieh ; if (obc%segnum_v(i,j) /= 0) then

1776 if (obc%segment(abs(obc%segnum_v(i,j)))%specified) obc_in_row = .true.

1777 endif ; enddo ; endif

1778 do k=1,nz

1779 ! This sets vh and dvhdv.

1780 call merid_flux_layer(v(:,j,k), h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), vh, dvhdv, ones, &

1781 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, por_face_areav(:,:,k), &

1782 cs%h_marg_min, obc)

1783 if (obc_in_row) then ; do i=ish,ieh ; if (obc%segnum_v(i,j) /= 0) then

1784 l_seg = abs(obc%segnum_v(i,j))

1785 if (obc%segment(l_seg)%specified) vh(i) = obc%segment(l_seg)%normal_trans(i,j,k)

1786 endif ; enddo ; endif

1787

1788 ! Accumulate the barotropic transport.

1789 do i=ish,ieh

1790 vhbt(i,j) = vhbt(i,j) + vh(i)

1791 enddo

1792 enddo ! k-loop

1793 enddo ! j-loop

1794

1795 call cpu_clock_end(id_clock_correct)

1796

1797end subroutine meridional_bt_mass_flux

1798

1799

1800!> Evaluates the meridional mass or volume fluxes in a layer.

1801subroutine merid_flux_layer(v, h, h_S, h_N, vh, dvhdv, visc_rem, dt, G, US, J, &

1802 ish, ieh, do_I, vol_CFL, por_face_areaV, h_marg_min, OBC)

1803 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

1804 real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

1805 real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the

1806 !! momentum originally in a layer that remains after a time-step

1807 !! of viscosity, and the fraction of a time-step's worth of a barotropic

1808 !! acceleration that a layer experiences after viscosity is applied [nondim].

1809 !! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).

1810 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,

1811 !! [H ~> m or kg m-2].

1812 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_S !< South edge thickness in the reconstruction

1813 !! [H ~> m or kg m-2].

1814 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_N !< North edge thickness in the reconstruction

1815 !! [H ~> m or kg m-2].

1816 real, dimension(SZI_(G)), intent(inout) :: vh !< Meridional mass or volume transport

1817 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

1818 real, dimension(SZI_(G)), intent(inout) :: dvhdv !< Partial derivative of vh with v

1819 !! [H L ~> m2 or kg m-1].

1820 real, intent(in) :: dt !< Time increment [T ~> s].

1821 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

1822 integer, intent(in) :: j !< Spatial index.

1823 integer, intent(in) :: ish !< Start of index range.

1824 integer, intent(in) :: ieh !< End of index range.

1825 logical, dimension(SZI_(G)), intent(in) :: do_I !< Which i values to work on.

1826 logical, intent(in) :: vol_CFL !< If true, rescale the

1827 !! ratio of face areas to the cell areas when estimating the CFL number.

1828 real, dimension(SZI_(G),SZJB_(G)), &

1829 intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]

1830 real, intent(in) :: h_marg_min !< Negligible floor on h_marg [H ~> m or kg m-2]

1831 type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.

1832 ! Local variables

1833 real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]

1834 real :: curv_3 ! A measure of the thickness curvature over a grid length,

1835 ! with the same units as h, i.e. [H ~> m or kg m-2].

1836 real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].

1837 integer :: i

1838 logical :: local_open_BC

1839

1840 local_open_bc = .false.

1841 if (present(obc)) then ; if (associated(obc)) then

1842 local_open_bc = obc%open_v_BCs_exist_globally

1843 endif ; endif

1844

1845 do i=ish,ieh ; if (do_i(i)) then

1846 if (v(i) > 0.0) then

1847 if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))

1848 else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif

1849 curv_3 = (h_s(i,j) + h_n(i,j)) - 2.0*h(i,j)

1850 vh(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * v(i) * ( h_n(i,j) + cfl * &

1851 (0.5*(h_s(i,j) - h_n(i,j)) + curv_3*(cfl - 1.5)) )

1852 h_marg = h_n(i,j) + cfl * ((h_s(i,j) - h_n(i,j)) + &

1853 3.0*curv_3*(cfl - 1.0))

1854 elseif (v(i) < 0.0) then

1855 if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))

1856 else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif

1857 curv_3 = (h_s(i,j+1) + h_n(i,j+1)) - 2.0*h(i,j+1)

1858 vh(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * v(i) * ( h_s(i,j+1) + cfl * &

1859 (0.5*(h_n(i,j+1)-h_s(i,j+1)) + curv_3*(cfl - 1.5)) )

1860 h_marg = h_s(i,j+1) + cfl * ((h_n(i,j+1)-h_s(i,j+1)) + &

1861 3.0*curv_3*(cfl - 1.0))

1862 else

1863 vh(i) = 0.0

1864 h_marg = 0.5 * (h_s(i,j+1) + h_n(i,j))

1865 endif

1866 h_marg = max(h_marg, h_marg_min)

1867 dvhdv(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * h_marg * visc_rem(i)

1868 endif ; enddo

1869

1870 if (local_open_bc) then

1871 do i=ish,ieh ; if (do_i(i)) then

1872 if (obc%segnum_v(i,j) /= 0) then

1873 if (obc%segment(abs(obc%segnum_v(i,j)))%open) then

1874 if (obc%segnum_v(i,j) > 0) then ! OBC_DIRECTION_N

1875 vh(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * v(i) * h(i,j)

1876 dvhdv(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * max(h(i,j), h_marg_min) * visc_rem(i)

1877 else

1878 vh(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * v(i) * h(i,j+1)

1879 dvhdv(i) = (g%dx_Cv(i,j)*por_face_areav(i,j)) * max(h(i,j+1), h_marg_min) * visc_rem(i)

1880 endif

1881 endif

1882 endif

1883 endif ; enddo

1884 endif

1885end subroutine merid_flux_layer

1886

1887!> Sets the effective interface thickness associated with the fluxes at each meridional velocity point,

1888!! optionally scaling back these thicknesses to account for viscosity and fractional open areas.

1889subroutine meridional_flux_thickness(v, h, h_S, h_N, h_v, dt, G, GV, US, LB, vol_CFL, &

1890 marginal, OBC, por_face_areaV, visc_rem_v)

1891 type(ocean_grid_type), intent(in) :: g !< Ocean's grid structure.

1892 type(verticalgrid_type), intent(in) :: gv !< Ocean's vertical grid structure.

1893 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

1894 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to calculate fluxes,

1895 !! [H ~> m or kg m-2].

1896 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_s !< South edge thickness in the reconstruction,

1897 !! [H ~> m or kg m-2].

1898 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_n !< North edge thickness in the reconstruction,

1899 !! [H ~> m or kg m-2].

1900 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: h_v !< Effective thickness at meridional faces,

1901 !! scaled down to account for the effects of

1902 !! viscosity and the fractional open area

1903 !! [H ~> m or kg m-2].

1904 real, intent(in) :: dt !< Time increment [T ~> s].

1905 type(cont_loop_bounds_type), intent(in) :: lb !< Loop bounds structure.

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

1907 logical, intent(in) :: vol_cfl !< If true, rescale the ratio

1908 !! of face areas to the cell areas when estimating the CFL number.

1909 logical, intent(in) :: marginal !< If true, report the marginal

1910 !! face thicknesses; otherwise report transport-averaged thicknesses.

1911 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

1912 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &

1913 intent(in) :: por_face_areav !< fractional open area of V-faces [nondim]

1914 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), optional, intent(in) :: visc_rem_v !< Both the fraction

1915 !! of the momentum originally in a layer that remains after a time-step of

1916 !! viscosity, and the fraction of a time-step's worth of a barotropic

1917 !! acceleration that a layer experiences after viscosity is applied [nondim].

1918 !! Visc_rem_v is between 0 (at the bottom) and 1 (far above the bottom).

1919

1920 ! Local variables

1921 real :: cfl ! The CFL number based on the local velocity and grid spacing [nondim]

1922 real :: curv_3 ! A measure of the thickness curvature over a grid length,

1923 ! with the same units as h [H ~> m or kg m-2] .

1924 real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].

1925 real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].

1926 logical :: local_open_bc

1927 integer :: i, j, k, ish, ieh, jsh, jeh, n, nz

1928 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = gv%ke

1929

1930 !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)

1931 do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh

1932 if (v(i,j,k) > 0.0) then

1933 if (vol_cfl) then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))

1934 else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ; endif

1935 curv_3 = (h_s(i,j,k) + h_n(i,j,k)) - 2.0*h(i,j,k)

1936 h_avg = h_n(i,j,k) + cfl * (0.5*(h_s(i,j,k) - h_n(i,j,k)) + curv_3*(cfl - 1.5))

1937 h_marg = h_n(i,j,k) + cfl * ((h_s(i,j,k) - h_n(i,j,k)) + &

1938 3.0*curv_3*(cfl - 1.0))

1939 elseif (v(i,j,k) < 0.0) then

1940 if (vol_cfl) then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))

1941 else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ; endif

1942 curv_3 = (h_s(i,j+1,k) + h_n(i,j+1,k)) - 2.0*h(i,j+1,k)

1943 h_avg = h_s(i,j+1,k) + cfl * (0.5*(h_n(i,j+1,k)-h_s(i,j+1,k)) + curv_3*(cfl - 1.5))

1944 h_marg = h_s(i,j+1,k) + cfl * ((h_n(i,j+1,k)-h_s(i,j+1,k)) + &

1945 3.0*curv_3*(cfl - 1.0))

1946 else

1947 h_avg = 0.5 * (h_s(i,j+1,k) + h_n(i,j,k))

1948 ! The choice to use the arithmetic mean here is somewhat arbitrarily, but

1949 ! it should be noted that h_S(i+1,j,k) and h_N(i,j,k) are usually the same.

1950 h_marg = 0.5 * (h_s(i,j+1,k) + h_n(i,j,k))

1951 ! h_marg = (2.0 * h_S(i,j+1,k) * h_N(i,j,k)) / &

1952 ! (h_S(i,j+1,k) + h_N(i,j,k) + GV%H_subroundoff)

1953 endif

1954

1955 if (marginal) then ; h_v(i,j,k) = h_marg

1956 else ; h_v(i,j,k) = h_avg ; endif

1957 enddo ; enddo ; enddo

1958

1959 if (present(visc_rem_v)) then

1960 ! Scale back the thickness to account for the effects of viscosity and the fractional open

1961 ! thickness to give an appropriate non-normalized weight for each layer in determining the

1962 ! barotropic acceleration.

1963 !$OMP parallel do default(shared)

1964 do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh

1965 h_v(i,j,k) = h_v(i,j,k) * (visc_rem_v(i,j,k) * por_face_areav(i,j,k))

1966 enddo ; enddo ; enddo

1967 else

1968 !$OMP parallel do default(shared)

1969 do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh

1970 h_v(i,j,k) = h_v(i,j,k) * por_face_areav(i,j,k)

1971 enddo ; enddo ; enddo

1972 endif

1973

1974 local_open_bc = .false.

1975 if (associated(obc)) local_open_bc = obc%open_v_BCs_exist_globally

1976 if (local_open_bc) then

1977 do n = 1, obc%number_of_segments

1978 if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then

1979 j = obc%segment(n)%HI%JsdB

1980 if (obc%segment(n)%direction == obc_direction_n) then

1981 if (present(visc_rem_v)) then ; do k=1,nz

1982 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied

1983 h_v(i,j,k) = h(i,j,k) * (visc_rem_v(i,j,k) * por_face_areav(i,j,k))

1984 enddo

1985 enddo ; else ; do k=1,nz

1986 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied

1987 h_v(i,j,k) = h(i,j,k) * por_face_areav(i,j,k)

1988 enddo

1989 enddo ; endif

1990 else

1991 if (present(visc_rem_v)) then ; do k=1,nz

1992 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied

1993 h_v(i,j,k) = h(i,j+1,k) * (visc_rem_v(i,j,k) * por_face_areav(i,j,k))

1994 enddo

1995 enddo ; else ; do k=1,nz

1996 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied

1997 h_v(i,j,k) = h(i,j+1,k) * por_face_areav(i,j,k)

1998 enddo

1999 enddo ; endif

2000 endif

2001 endif

2002 enddo

2003 endif

2004

2005end subroutine meridional_flux_thickness

2006

2007!> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.

2008subroutine meridional_flux_adjust(v, h_in, h_S, h_N, vhbt, vh_tot_0, dvhdv_tot_0, &

2009 dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &

2010 j, ish, ieh, do_I_in, por_face_areaV, vh_3d, OBC)

2011 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

2012 type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure.

2013 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

2014 intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

2015 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

2016 intent(in) :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].

2017 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),&

2018 intent(in) :: h_S !< South edge thickness in the reconstruction [H ~> m or kg m-2].

2019 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

2020 intent(in) :: h_N !< North edge thickness in the reconstruction [H ~> m or kg m-2].

2021 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: visc_rem

2022 !< Both the fraction of the momentum originally

2023 !! in a layer that remains after a time-step of viscosity, and the

2024 !! fraction of a time-step's worth of a barotropic acceleration that

2025 !! a layer experiences after viscosity is applied [nondim].

2026 !! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).

2027 real, dimension(SZI_(G)), intent(in) :: vhbt !< The summed volume flux through meridional faces

2028 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

2029 real, dimension(SZI_(G)), intent(in) :: dv_max_CFL !< Maximum acceptable value of dv [L T-1 ~> m s-1].

2030 real, dimension(SZI_(G)), intent(in) :: dv_min_CFL !< Minimum acceptable value of dv [L T-1 ~> m s-1].

2031 real, dimension(SZI_(G)), intent(in) :: vh_tot_0 !< The summed transport with 0 adjustment

2032 !! [H L2 T-1 ~> m3 s-1 or kg s-1].

2033 real, dimension(SZI_(G)), intent(in) :: dvhdv_tot_0 !< The partial derivative of dv_err with

2034 !! dv at 0 adjustment [H L ~> m2 or kg m-1].

2035 real, dimension(SZI_(G)), intent(out) :: dv !< The barotropic velocity adjustment [L T-1 ~> m s-1].

2036 real, intent(in) :: dt !< Time increment [T ~> s].

2037 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

2038 type(continuity_ppm_cs), intent(in) :: CS !< This module's control structure.

2039 integer, intent(in) :: j !< Spatial index.

2040 integer, intent(in) :: ish !< Start of index range.

2041 integer, intent(in) :: ieh !< End of index range.

2042 logical, dimension(SZI_(G)), &

2043 intent(in) :: do_I_in !< A flag indicating which I values to work on.

2044 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &

2045 intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]

2046 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

2047 optional, intent(inout) :: vh_3d !< Volume flux through meridional

2048 !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].

2049 type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.

2050 ! Local variables

2051 real, dimension(SZI_(G),SZK_(GV)) :: &

2052 vh_aux, & ! An auxiliary meridional volume flux [H L2 T-1 ~> m3 s-1 or kg s-1].

2053 dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].

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

2055 vh_err, & ! Difference between vhbt and the summed vh [H L2 T-1 ~> m3 s-1 or kg s-1].

2056 vh_err_best, & ! The smallest value of vh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].

2057 v_new, & ! The velocity with the correction added [L T-1 ~> m s-1].

2058 dvhdv_tot,&! Summed partial derivative of vh with u [H L ~> m2 or kg m-1].

2059 dv_min, & ! Lower limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]

2060 dv_max ! Upper limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]

2061 real :: dv_prev ! The previous value of dv [L T-1 ~> m s-1].

2062 real :: ddv ! The change in dv from the previous iteration [L T-1 ~> m s-1].

2063 real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].

2064 real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].

2065 integer :: i, k, nz, itt, max_itts = 20

2066 logical :: domore, do_I(SZI_(G))

2067

2068 nz = gv%ke

2069

2070 vh_aux(:,:) = 0.0 ; dvhdv(:,:) = 0.0

2071

2072 if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh

2073 vh_aux(i,k) = vh_3d(i,j,k)

2074 enddo ; enddo ; endif

2075

2076 do i=ish,ieh

2077 dv(i) = 0.0 ; do_i(i) = do_i_in(i)

2078 dv_max(i) = dv_max_cfl(i) ; dv_min(i) = dv_min_cfl(i)

2079 vh_err(i) = vh_tot_0(i) - vhbt(i) ; dvhdv_tot(i) = dvhdv_tot_0(i)

2080 vh_err_best(i) = abs(vh_err(i))

2081 enddo

2082

2083 do itt=1,max_itts

2084 select case (itt)

2085 case (:1) ; tol_eta = 1e-6 * cs%tol_eta

2086 case (2) ; tol_eta = 1e-4 * cs%tol_eta

2087 case (3) ; tol_eta = 1e-2 * cs%tol_eta

2088 case default ; tol_eta = cs%tol_eta

2089 end select

2090 tol_vel = cs%tol_vel

2091

2092 do i=ish,ieh

2093 if (vh_err(i) > 0.0) then ; dv_max(i) = dv(i)

2094 elseif (vh_err(i) < 0.0) then ; dv_min(i) = dv(i)

2095 else ; do_i(i) = .false. ; endif

2096 enddo

2097 domore = .false.

2098 do i=ish,ieh ; if (do_i(i)) then

2099 if ((dt * min(g%IareaT(i,j),g%IareaT(i,j+1))*abs(vh_err(i)) > tol_eta) .or. &

2100 (cs%better_iter .and. ((abs(vh_err(i)) > tol_vel * dvhdv_tot(i)) .or. &

2101 (abs(vh_err(i)) > vh_err_best(i))) )) then

2102 ! Use Newton's method, provided it stays bounded. Otherwise bisect

2103 ! the value with the appropriate bound.

2104 ddv = -vh_err(i) / dvhdv_tot(i)

2105 dv_prev = dv(i)

2106 dv(i) = dv(i) + ddv

2107 if (abs(ddv) < 1.0e-15*abs(dv(i))) then

2108 do_i(i) = .false. ! ddv is small enough to quit.

2109 elseif (ddv > 0.0) then

2110 if (dv(i) >= dv_max(i)) then

2111 dv(i) = 0.5*(dv_prev + dv_max(i))

2112 if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_i(i) = .false.

2113 endif

2114 else ! dvv(i) < 0.0

2115 if (dv(i) <= dv_min(i)) then

2116 dv(i) = 0.5*(dv_prev + dv_min(i))

2117 if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_i(i) = .false.

2118 endif

2119 endif

2120 if (do_i(i)) domore = .true.

2121 else

2122 do_i(i) = .false.

2123 endif

2124 endif ; enddo

2125 if (.not.domore) exit

2126

2127 if ((itt < max_itts) .or. present(vh_3d)) then ; do k=1,nz

2128 do i=ish,ieh ; v_new(i) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo

2129 call merid_flux_layer(v_new, h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), &

2130 vh_aux(:,k), dvhdv(:,k), visc_rem(:,k), &

2131 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, por_face_areav(:,:,k), &

2132 cs%h_marg_min, obc)

2133 enddo ; endif

2134

2135 if (itt < max_itts) then

2136 do i=ish,ieh

2137 vh_err(i) = -vhbt(i) ; dvhdv_tot(i) = 0.0

2138 enddo

2139 do k=1,nz ; do i=ish,ieh

2140 vh_err(i) = vh_err(i) + vh_aux(i,k)

2141 dvhdv_tot(i) = dvhdv_tot(i) + dvhdv(i,k)

2142 enddo ; enddo

2143 do i=ish,ieh

2144 vh_err_best(i) = min(vh_err_best(i), abs(vh_err(i)))

2145 enddo

2146 endif

2147 enddo ! itt-loop

2148 ! If there are any faces which have not converged to within the tolerance,

2149 ! so-be-it, or else use a final upwind correction?

2150 ! This never seems to happen with 20 iterations as max_itt.

2151

2152 if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh

2153 vh_3d(i,j,k) = vh_aux(i,k)

2154 enddo ; enddo ; endif

2155

2156end subroutine meridional_flux_adjust

2157

2158!> Sets of a structure that describes the meridional barotropic volume or mass fluxes as a

2159!! function of barotropic flow to agree closely with the sum of the layer's transports.

2160subroutine set_merid_bt_cont(v, h_in, h_S, h_N, BT_cont, vh_tot_0, dvhdv_tot_0, &

2161 dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &

2162 visc_rem_max, j, ish, ieh, do_I, por_face_areaV)

2163 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

2164 type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure.

2165 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].

2166 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to calculate fluxes,

2167 !! [H ~> m or kg m-2].

2168 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_S !< South edge thickness in the reconstruction,

2169 !! [H ~> m or kg m-2].

2170 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_N !< North edge thickness in the reconstruction,

2171 !! [H ~> m or kg m-2].

2172 type(bt_cont_type), intent(inout) :: BT_cont !< A structure with elements

2173 !! that describe the effective open face areas as a function of barotropic flow.

2174 real, dimension(SZI_(G)), intent(in) :: vh_tot_0 !< The summed transport

2175 !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].

2176 real, dimension(SZI_(G)), intent(in) :: dvhdv_tot_0 !< The partial derivative

2177 !! of du_err with dv at 0 adjustment [H L ~> m2 or kg m-1].

2178 real, dimension(SZI_(G)), intent(in) :: dv_max_CFL !< Maximum acceptable value

2179 !! of dv [L T-1 ~> m s-1].

2180 real, dimension(SZI_(G)), intent(in) :: dv_min_CFL !< Minimum acceptable value

2181 !! of dv [L T-1 ~> m s-1].

2182 real, intent(in) :: dt !< Time increment [T ~> s].

2183 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

2184 type(continuity_ppm_cs), intent(in) :: CS !< This module's control structure.

2185 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: visc_rem !< Both the fraction of the

2186 !! momentum originally in a layer that remains after a time-step

2187 !! of viscosity, and the fraction of a time-step's worth of a barotropic

2188 !! acceleration that a layer experiences after viscosity is applied [nondim].

2189 !! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).

2190 real, dimension(SZI_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem [nondim]

2191 integer, intent(in) :: j !< Spatial index.

2192 integer, intent(in) :: ish !< Start of index range.

2193 integer, intent(in) :: ieh !< End of index range.

2194 logical, dimension(SZI_(G)), intent(in) :: do_I !< A logical flag indicating

2195 !! which I values to work on.

2196 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &

2197 intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]

2198 ! Local variables

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

2200 dv0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].

2201 dvL, dvR, & ! The barotropic velocity increments that give the southerly

2202 ! (dvL) and northerly (dvR) test velocities [L T-1 ~> m s-1].

2203 zeros, & ! An array of full of 0 transports [H L2 T-1 ~> m3 s-1 or kg s-1]

2204 dv_cfl, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].

2205 v_l, v_r, & ! The southerly (v_L), northerly (v_R), and zero-barotropic

2206 v_0, & ! transport (v_0) layer test velocities [L T-1 ~> m s-1].

2207 dvhdv_l, & ! The effective layer marginal face areas with the southerly

2208 dvhdv_r, & ! (_L), northerly (_R), and zero-barotropic (_0) test

2209 dvhdv_0, & ! velocities [H L ~> m2 or kg m-1].

2210 vh_l, vh_r, & ! The layer transports with the southerly (_L), northerly (_R)

2211 vh_0, & ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].

2212 famt_l, famt_r, & ! The summed effective marginal face areas for the 3

2213 famt_0, & ! test velocities [H L ~> m2 or kg m-1].

2214 vhtot_l, & ! The summed transport with the southerly (vhtot_L) and

2215 vhtot_r ! and northerly (vhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].

2216 real :: FA_0 ! The effective face area with 0 barotropic transport [H L ~> m2 or kg m-1].

2217 real :: FA_avg ! The average effective face area [H L ~> m2 or kg m-1], nominally given by

2218 ! the realized transport divided by the barotropic velocity.

2219 real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim] This

2220 ! limiting is necessary to keep the inverse of visc_rem

2221 ! from leading to large CFL numbers.

2222 real :: min_visc_rem ! The smallest permitted value for visc_rem that is used

2223 ! in finding the barotropic velocity that changes the

2224 ! flow direction [nondim]. This is necessary to keep the inverse

2225 ! of visc_rem from leading to large CFL numbers.

2226 real :: CFL_min ! A minimal increment in the CFL to try to ensure that the

2227 ! flow is truly upwind [nondim]

2228 real :: Idt ! The inverse of the time step [T-1 ~> s-1].

2229 logical :: domore

2230 integer :: i, k, nz

2231

2232 nz = gv%ke ; idt = 1.0 / dt

2233 min_visc_rem = 0.1 ; cfl_min = 1e-6

2234

2235 ! Diagnose the zero-transport correction, dv0.

2236 do i=ish,ieh ; zeros(i) = 0.0 ; enddo

2237 call meridional_flux_adjust(v, h_in, h_s, h_n, zeros, vh_tot_0, dvhdv_tot_0, dv0, &

2238 dv_max_cfl, dv_min_cfl, dt, g, gv, us, cs, visc_rem, &

2239 j, ish, ieh, do_i, por_face_areav)

2240

2241 ! Determine the southerly- and northerly- fluxes. Choose a sufficiently

2242 ! negative velocity correction for the northerly-flux, and a sufficiently

2243 ! positive correction for the southerly-flux.

2244 domore = .false.

2245 do i=ish,ieh ; if (do_i(i)) then

2246 domore = .true.

2247 dv_cfl(i) = (cfl_min * idt) * g%dyCv(i,j)

2248 dvr(i) = min(0.0,dv0(i) - dv_cfl(i))

2249 dvl(i) = max(0.0,dv0(i) + dv_cfl(i))

2250 famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0

2251 vhtot_l(i) = 0.0 ; vhtot_r(i) = 0.0

2252 endif ; enddo

2253

2254 if (.not.domore) then

2255 do k=1,nz ; do i=ish,ieh

2256 bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0

2257 bt_cont%vBT_SS(i,j) = 0.0

2258 bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0

2259 bt_cont%vBT_NN(i,j) = 0.0

2260 enddo ; enddo

2261 return

2262 endif

2263

2264 do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then

2265 visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))

2266 if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.

2267 if (v(i,j,k) + dvr(i)*visc_rem_lim > -dv_cfl(i)*visc_rem(i,k)) &

2268 dvr(i) = -(v(i,j,k) + dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim

2269 if (v(i,j,k) + dvl(i)*visc_rem_lim < dv_cfl(i)*visc_rem(i,k)) &

2270 dvl(i) = -(v(i,j,k) - dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim

2271 endif

2272 endif ; enddo ; enddo

2273 do k=1,nz

2274 do i=ish,ieh ; if (do_i(i)) then

2275 v_l(i) = v(i,j,k) + dvl(i) * visc_rem(i,k)

2276 v_r(i) = v(i,j,k) + dvr(i) * visc_rem(i,k)

2277 v_0(i) = v(i,j,k) + dv0(i) * visc_rem(i,k)

2278 endif ; enddo

2279 call merid_flux_layer(v_0, h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), vh_0, dvhdv_0, &

2280 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, &

2281 por_face_areav(:,:,k), cs%h_marg_min)

2282 call merid_flux_layer(v_l, h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), vh_l, dvhdv_l, &

2283 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, &

2284 por_face_areav(:,:,k), cs%h_marg_min)

2285 call merid_flux_layer(v_r, h_in(:,:,k), h_s(:,:,k), h_n(:,:,k), vh_r, dvhdv_r, &

2286 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, &

2287 por_face_areav(:,:,k), cs%h_marg_min)

2288 do i=ish,ieh ; if (do_i(i)) then

2289 famt_0(i) = famt_0(i) + dvhdv_0(i)

2290 famt_l(i) = famt_l(i) + dvhdv_l(i)

2291 famt_r(i) = famt_r(i) + dvhdv_r(i)

2292 vhtot_l(i) = vhtot_l(i) + vh_l(i)

2293 vhtot_r(i) = vhtot_r(i) + vh_r(i)

2294 endif ; enddo

2295 enddo

2296 do i=ish,ieh ; if (do_i(i)) then

2297 fa_0 = famt_0(i) ; fa_avg = famt_0(i)

2298 if ((dvl(i) - dv0(i)) /= 0.0) &

2299 fa_avg = vhtot_l(i) / (dvl(i) - dv0(i))

2300 if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))

2301 elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif

2302 bt_cont%FA_v_S0(i,j) = fa_0 ; bt_cont%FA_v_SS(i,j) = famt_l(i)

2303 if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%vBT_SS(i,j) = 0.0 ; else

2304 bt_cont%vBT_SS(i,j) = (1.5 * (dvl(i) - dv0(i))) * &

2305 ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))

2306 endif

2307

2308 fa_0 = famt_0(i) ; fa_avg = famt_0(i)

2309 if ((dvr(i) - dv0(i)) /= 0.0) &

2310 fa_avg = vhtot_r(i) / (dvr(i) - dv0(i))

2311 if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))

2312 elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif

2313 bt_cont%FA_v_N0(i,j) = fa_0 ; bt_cont%FA_v_NN(i,j) = famt_r(i)

2314 if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%vBT_NN(i,j) = 0.0 ; else

2315 bt_cont%vBT_NN(i,j) = (1.5 * (dvr(i) - dv0(i))) * &

2316 ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))

2317 endif

2318 else

2319 bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0

2320 bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0

2321 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0

2322 endif ; enddo

2323

2324end subroutine set_merid_bt_cont

2325

2326!> Calculates left/right edge values for PPM reconstruction.

2327subroutine ppm_reconstruction_x(h_in, h_W, h_E, G, LB, h_min, monotonic, simple_2nd, OBC, k)

2328 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

2329 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].

2330 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_W !< West edge thickness in the reconstruction,

2331 !! [H ~> m or kg m-2].

2332 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_E !< East edge thickness in the reconstruction,

2333 !! [H ~> m or kg m-2].

2334 type(cont_loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.

2335 real, intent(in) :: h_min !< The minimum thickness

2336 !! that can be obtained by a concave parabolic fit [H ~> m or kg m-2]

2337 logical, intent(in) :: monotonic !< If true, use the

2338 !! Colella & Woodward monotonic limiter.

2339 !! Otherwise use a simple positive-definite limiter.

2340 logical, intent(in) :: simple_2nd !< If true, use the

2341 !! arithmetic mean thicknesses as the default edge values

2342 !! for a simple 2nd order scheme.

2343 type(ocean_obc_type), pointer :: OBC !< Open boundaries control structure.

2344 integer :: k !< vertical grid index

2345

2346 ! Local variables with useful mnemonic names.

2347 real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes per grid point [H ~> m or kg m-2]

2348 real, parameter :: oneSixth = 1./6. ! [nondim]

2349 real :: h_ip1, h_im1 ! Neighboring thicknesses or sensibly extrapolated values [H ~> m or kg m-2]

2350 real :: dMx, dMn ! The difference between the local thickness and the maximum (dMx) or

2351 ! minimum (dMn) of the surrounding values [H ~> m or kg m-2]

2352 character(len=256) :: mesg

2353 integer :: i, j, isl, iel, jsl, jel, n, stencil

2354 logical :: local_open_BC

2355 type(obc_segment_type), pointer :: segment => null()

2356

2357 local_open_bc = .false.

2358 if (associated(obc)) then

2359 local_open_bc = obc%open_u_BCs_exist_globally

2360 endif

2361

2362 isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh

2363

2364 ! This is the stencil of the reconstruction, not the scheme overall.

2365 stencil = 2 ; if (simple_2nd) stencil = 1

2366

2367 if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then

2368 write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &

2369 & "x-halo that needs to be increased by ",I0,".")') &

2370 stencil + max(g%isd-isl,iel-g%ied)

2371 call mom_error(fatal,mesg)

2372 endif

2373 if ((jsl < g%jsd) .or. (jel > g%jed)) then

2374 write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &

2375 & "y-halo that needs to be increased by ",I0,".")') &

2376 max(g%jsd-jsl,jel-g%jed)

2377 call mom_error(fatal,mesg)

2378 endif

2379

2380 if (simple_2nd) then

2381 do j=jsl,jel ; do i=isl,iel

2382 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)

2383 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)

2384 h_w(i,j) = 0.5*( h_im1 + h_in(i,j) )

2385 h_e(i,j) = 0.5*( h_ip1 + h_in(i,j) )

2386 enddo ; enddo

2387 else

2388 do j=jsl,jel ; do i=isl-1,iel+1

2389 if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then

2390 slp(i,j) = 0.0

2391 else

2392 ! This uses a simple 2nd order slope.

2393 slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))

2394 ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)

2395 dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)

2396 dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))

2397 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))

2398 ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))

2399 endif

2400 enddo ; enddo

2401

2402 if (local_open_bc) then

2403 do n=1, obc%number_of_segments

2404 segment => obc%segment(n)

2405 if (.not. segment%on_pe) cycle

2406 if (segment%is_E_or_W) then

2407 i=segment%HI%IsdB

2408 do j=segment%HI%jsd,segment%HI%jed

2409 slp(i+1,j) = 0.0

2410 slp(i,j) = 0.0

2411 enddo

2412 endif

2413 enddo

2414 endif

2415

2416 do j=jsl,jel ; do i=isl,iel

2417 ! Neighboring values should take into account any boundaries. The 3

2418 ! following sets of expressions are equivalent.

2419 ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)

2420 ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)

2421 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)

2422 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)

2423 ! Left/right values following Eq. B2 in Lin 1994, MWR (132)

2424 h_w(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )

2425 h_e(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )

2426 enddo ; enddo

2427 endif

2428

2429 if (local_open_bc) then

2430 do n=1, obc%number_of_segments

2431 segment => obc%segment(n)

2432 if (.not. segment%on_pe) cycle

2433 if (segment%direction == obc_direction_e) then

2434 i=segment%HI%IsdB

2435 if (associated(segment%h_Reg)) then

2436 do j=segment%HI%jsd,segment%HI%jed

2437 h_w(i+1,j) = segment%h_Reg%h_res(i,j,k)

2438 h_e(i+1,j) = segment%h_Reg%h_res(i,j,k)

2439 h_w(i,j) = segment%h_Reg%h_res(i,j,k)

2440 h_e(i,j) = segment%h_Reg%h_res(i,j,k)

2441 enddo

2442 else

2443 do j=segment%HI%jsd,segment%HI%jed

2444 h_w(i+1,j) = h_in(i,j)

2445 h_e(i+1,j) = h_in(i,j)

2446 h_w(i,j) = h_in(i,j)

2447 h_e(i,j) = h_in(i,j)

2448 enddo

2449 endif

2450 elseif (segment%direction == obc_direction_w) then

2451 i=segment%HI%IsdB

2452 if (associated(segment%h_Reg)) then

2453 do j=segment%HI%jsd,segment%HI%jed

2454 h_w(i,j) = segment%h_Reg%h_res(i,j,k)

2455 h_e(i,j) = segment%h_Reg%h_res(i,j,k)

2456 h_w(i+1,j) = segment%h_Reg%h_res(i,j,k)

2457 h_e(i+1,j) = segment%h_Reg%h_res(i,j,k)

2458 enddo

2459 else

2460 do j=segment%HI%jsd,segment%HI%jed

2461 h_w(i,j) = h_in(i+1,j)

2462 h_e(i,j) = h_in(i+1,j)

2463 h_w(i+1,j) = h_in(i+1,j)

2464 h_e(i+1,j) = h_in(i+1,j)

2465 enddo

2466 endif

2467 endif

2468 enddo

2469 endif

2470

2471 if (monotonic) then

2472 call ppm_limit_cw84(h_in, h_w, h_e, g, isl, iel, jsl, jel)

2473 else

2474 call ppm_limit_pos(h_in, h_w, h_e, h_min, g, isl, iel, jsl, jel)

2475 endif

2476

2477 return

2478end subroutine ppm_reconstruction_x

2479

2480!> Calculates left/right edge values for PPM reconstruction.

2481subroutine ppm_reconstruction_y(h_in, h_S, h_N, G, LB, h_min, monotonic, simple_2nd, OBC, k)

2482 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

2483 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].

2484 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_S !< South edge thickness in the reconstruction,

2485 !! [H ~> m or kg m-2].

2486 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_N !< North edge thickness in the reconstruction,

2487 !! [H ~> m or kg m-2].

2488 type(cont_loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.

2489 real, intent(in) :: h_min !< The minimum thickness

2490 !! that can be obtained by a concave parabolic fit [H ~> m or kg m-2]

2491 logical, intent(in) :: monotonic !< If true, use the

2492 !! Colella & Woodward monotonic limiter.

2493 !! Otherwise use a simple positive-definite limiter.

2494 logical, intent(in) :: simple_2nd !< If true, use the

2495 !! arithmetic mean thicknesses as the default edge values

2496 !! for a simple 2nd order scheme.

2497 type(ocean_obc_type), pointer :: OBC !< Open boundaries control structure.

2498 integer :: k !< vertical grid index

2499

2500 ! Local variables with useful mnemonic names.

2501 real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes per grid point [H ~> m or kg m-2]

2502 real, parameter :: oneSixth = 1./6. ! [nondim]

2503 real :: h_jp1, h_jm1 ! Neighboring thicknesses or sensibly extrapolated values [H ~> m or kg m-2]

2504 real :: dMx, dMn ! The difference between the local thickness and the maximum (dMx) or

2505 ! minimum (dMn) of the surrounding values [H ~> m or kg m-2]

2506 character(len=256) :: mesg

2507 integer :: i, j, isl, iel, jsl, jel, n, stencil

2508 logical :: local_open_BC

2509 type(obc_segment_type), pointer :: segment => null()

2510

2511 local_open_bc = .false.

2512 if (associated(obc)) then

2513 local_open_bc = obc%open_v_BCs_exist_globally

2514 endif

2515

2516 isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1

2517

2518 ! This is the stencil of the reconstruction, not the scheme overall.

2519 stencil = 2 ; if (simple_2nd) stencil = 1

2520

2521 if ((isl < g%isd) .or. (iel > g%ied)) then

2522 write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &

2523 & "x-halo that needs to be increased by ",I0,".")') &

2524 max(g%isd-isl,iel-g%ied)

2525 call mom_error(fatal,mesg)

2526 endif

2527 if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then

2528 write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &

2529 & "y-halo that needs to be increased by ",I0,".")') &

2530 stencil + max(g%jsd-jsl,jel-g%jed)

2531 call mom_error(fatal,mesg)

2532 endif

2533

2534 if (simple_2nd) then

2535 do j=jsl,jel ; do i=isl,iel

2536 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)

2537 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)

2538 h_s(i,j) = 0.5*( h_jm1 + h_in(i,j) )

2539 h_n(i,j) = 0.5*( h_jp1 + h_in(i,j) )

2540 enddo ; enddo

2541 else

2542 do j=jsl-1,jel+1 ; do i=isl,iel

2543 if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then

2544 slp(i,j) = 0.0

2545 else

2546 ! This uses a simple 2nd order slope.

2547 slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))

2548 ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)

2549 dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)

2550 dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))

2551 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))

2552 ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))

2553 endif

2554 enddo ; enddo

2555

2556 if (local_open_bc) then

2557 do n=1, obc%number_of_segments

2558 segment => obc%segment(n)

2559 if (.not. segment%on_pe) cycle

2560 if (segment%is_N_or_S) then

2561 j=segment%HI%JsdB

2562 do i=segment%HI%isd,segment%HI%ied

2563 slp(i,j+1) = 0.0

2564 slp(i,j) = 0.0

2565 enddo

2566 endif

2567 enddo

2568 endif

2569

2570 do j=jsl,jel ; do i=isl,iel

2571 ! Neighboring values should take into account any boundaries. The 3

2572 ! following sets of expressions are equivalent.

2573 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)

2574 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)

2575 ! Left/right values following Eq. B2 in Lin 1994, MWR (132)

2576 h_s(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )

2577 h_n(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )

2578 enddo ; enddo

2579 endif

2580

2581 if (local_open_bc) then

2582 do n=1, obc%number_of_segments

2583 segment => obc%segment(n)

2584 if (.not. segment%on_pe) cycle

2585 if (segment%direction == obc_direction_n) then

2586 j=segment%HI%JsdB

2587 if (associated(segment%h_Reg)) then

2588 do i=segment%HI%isd,segment%HI%ied

2589 h_s(i,j+1) = segment%h_Reg%h_res(i,j,k)

2590 h_n(i,j+1) = segment%h_Reg%h_res(i,j,k)

2591 h_s(i,j) = segment%h_Reg%h_res(i,j,k)

2592 h_n(i,j) = segment%h_Reg%h_res(i,j,k)

2593 enddo

2594 else

2595 do i=segment%HI%isd,segment%HI%ied

2596 h_s(i,j+1) = h_in(i,j)

2597 h_n(i,j+1) = h_in(i,j)

2598 h_s(i,j) = h_in(i,j)

2599 h_n(i,j) = h_in(i,j)

2600 enddo

2601 endif

2602 elseif (segment%direction == obc_direction_s) then

2603 j=segment%HI%JsdB

2604 if (associated(segment%h_Reg)) then

2605 do i=segment%HI%isd,segment%HI%ied

2606 h_s(i,j) = segment%h_Reg%h_res(i,j,k)

2607 h_n(i,j) = segment%h_Reg%h_res(i,j,k)

2608 h_s(i,j+1) = segment%h_Reg%h_res(i,j,k)

2609 h_n(i,j+1) = segment%h_Reg%h_res(i,j,k)

2610 enddo

2611 else

2612 do i=segment%HI%isd,segment%HI%ied

2613 h_s(i,j) = h_in(i,j+1)

2614 h_n(i,j) = h_in(i,j+1)

2615 h_s(i,j+1) = h_in(i,j+1)

2616 h_n(i,j+1) = h_in(i,j+1)

2617 enddo

2618 endif

2619 endif

2620 enddo

2621 endif

2622

2623 if (monotonic) then

2624 call ppm_limit_cw84(h_in, h_s, h_n, g, isl, iel, jsl, jel)

2625 else

2626 call ppm_limit_pos(h_in, h_s, h_n, h_min, g, isl, iel, jsl, jel)

2627 endif

2628

2629 return

2630end subroutine ppm_reconstruction_y

2631

2632!> This subroutine limits the left/right edge values of the PPM reconstruction

2633!! to give a reconstruction that is positive-definite. Here this is

2634!! reinterpreted as giving a constant thickness if the mean thickness is less

2635!! than h_min, with a minimum of h_min otherwise.

2636subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)

2637 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

2638 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].

2639 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].

2640 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].

2641 real, intent(in) :: h_min !< The minimum thickness

2642 !! that can be obtained by a concave parabolic fit [H ~> m or kg m-2]

2643 integer, intent(in) :: iis !< Start of i index range.

2644 integer, intent(in) :: iie !< End of i index range.

2645 integer, intent(in) :: jis !< Start of j index range.

2646 integer, intent(in) :: jie !< End of j index range.

2647

2648! Local variables

2649 real :: curv ! The grid-normalized curvature of the three thicknesses [H ~> m or kg m-2]

2650 real :: dh ! The difference between the edge thicknesses [H ~> m or kg m-2]

2651 real :: scale ! A scaling factor to reduce the curvature of the fit [nondim]

2652 integer :: i,j

2653

2654 do j=jis,jie ; do i=iis,iie

2655 ! This limiter prevents undershooting minima within the domain with

2656 ! values less than h_min.

2657 curv = 3.0*((h_l(i,j) + h_r(i,j)) - 2.0*h_in(i,j))

2658 if (curv > 0.0) then ! Only minima are limited.

2659 dh = h_r(i,j) - h_l(i,j)

2660 if (abs(dh) < curv) then ! The parabola's minimum is within the cell.

2661 if (h_in(i,j) <= h_min) then

2662 h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)

2663 elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then

2664 ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must

2665 ! be limited in this case. 0 < scale < 1.

2666 scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)

2667 h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))

2668 h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))

2669 endif

2670 endif

2671 endif

2672 enddo ; enddo

2673

2674end subroutine ppm_limit_pos

2675

2676!> This subroutine limits the left/right edge values of the PPM reconstruction

2677!! according to the monotonic prescription of Colella and Woodward, 1984.

2678subroutine ppm_limit_cw84(h_in, h_L, h_R, G, iis, iie, jis, jie)

2679 type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.

2680 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].

2681 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction,

2682 !! [H ~> m or kg m-2].

2683 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction,

2684 !! [H ~> m or kg m-2].

2685 integer, intent(in) :: iis !< Start of i index range.

2686 integer, intent(in) :: iie !< End of i index range.

2687 integer, intent(in) :: jis !< Start of j index range.

2688 integer, intent(in) :: jie !< End of j index range.

2689

2690 ! Local variables

2691 real :: h_i ! A copy of the cell-average layer thickness [H ~> m or kg m-2]

2692 real :: RLdiff ! The difference between the input edge values [H ~> m or kg m-2]

2693 real :: RLdiff2 ! The squared difference between the input edge values [H2 ~> m2 or kg2 m-4]

2694 real :: RLmean ! The average of the input edge thicknesses [H ~> m or kg m-2]

2695 real :: FunFac ! A curious product of the thickness slope and curvature [H2 ~> m2 or kg2 m-4]

2696 integer :: i, j

2697

2698 do j=jis,jie ; do i=iis,iie

2699 ! This limiter monotonizes the parabola following

2700 ! Colella and Woodward, 1984, Eq. 1.10

2701 h_i = h_in(i,j)

2702 if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. ) then

2703 h_l(i,j) = h_i ; h_r(i,j) = h_i

2704 else

2705 rldiff = h_r(i,j) - h_l(i,j) ! Difference of edge values

2706 rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) ) ! Mean of edge values

2707 funfac = 6. * rldiff * ( h_i - rlmean ) ! Some funny factor

2708 rldiff2 = rldiff * rldiff ! Square of difference

2709 if ( funfac > rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)

2710 if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)

2711 endif

2712 enddo ; enddo

2713

2714 return

2715end subroutine ppm_limit_cw84

2716

2717!> Return the maximum ratio of a/b or maxrat.

2718function ratio_max(a, b, maxrat) result(ratio)

2719 real, intent(in) :: a !< Numerator, in arbitrary units [A]

2720 real, intent(in) :: b !< Denominator, in arbitrary units [B]

2721 real, intent(in) :: maxrat !< Maximum value of ratio [A B-1]

2722 real :: ratio !< Return value [A B-1]

2723

2724 if (abs(a) > abs(maxrat*b)) then

2725 ratio = maxrat

2726 else

2727 ratio = a / b

2728 endif

2729end function ratio_max

2730

2731!> Initializes continuity_ppm_cs

2732subroutine continuity_ppm_init(Time, G, GV, US, param_file, diag, CS, OBC)

2733 type(time_type), target, intent(in) :: time !< The current model time.

2734 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

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

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

2737 type(param_file_type), intent(in) :: param_file !< A structure indicating

2738 !! the open file to parse for model parameter values.

2739 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to

2740 !! regulate diagnostic output.

2741 type(continuity_ppm_cs), intent(inout) :: cs !< Module's control structure.

2742 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure.

2743 logical :: local_open_bc, use_h_marg_min

2744 type(obc_segment_type), pointer :: segment => null()

2745 integer :: n

2746

2747 !> This include declares and sets the variable "version".

2748# include "version_variable.h"

2749 character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name.

2750 character(len=256) :: mesg

2751

2752 cs%initialized = .true.

2753

2754 local_open_bc = .false.

2755 if (associated(obc)) then

2756 local_open_bc = obc%open_u_BCs_exist_globally

2757 endif

2758

2759! Read all relevant parameters and write them to the model log.

2760 call log_version(param_file, mdl, version, "")

2761 call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", cs%monotonic, &

2762 "If true, CONTINUITY_PPM uses the Colella and Woodward "//&

2763 "monotonic limiter. The default (false) is to use a "//&

2764 "simple positive definite limiter.", default=.false.)

2765 call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &

2766 "If true, CONTINUITY_PPM uses a simple 2nd order "//&

2767 "(arithmetic mean) interpolation of the edge values. "//&

2768 "This may give better PV conservation properties. While "//&

2769 "it formally reduces the accuracy of the continuity "//&

2770 "solver itself in the strongly advective limit, it does "//&

2771 "not reduce the overall order of accuracy of the dynamic "//&

2772 "core.", default=.false.)

2773 call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &

2774 "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&

2775 "continuity solver. This scheme is highly diffusive "//&

2776 "but may be useful for debugging or in single-column "//&

2777 "mode where its minimal stencil is useful.", default=.false.)

2778 call get_param(param_file, mdl, "ETA_TOLERANCE", cs%tol_eta, &

2779 "The tolerance for the differences between the "//&

2780 "barotropic and baroclinic estimates of the sea surface "//&

2781 "height due to the fluxes through each face. The total "//&

2782 "tolerance for SSH is 4 times this value. The default "//&

2783 "is 0.5*NK*ANGSTROM, and this should not be set less "//&

2784 "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=gv%m_to_H, &

2785 default=0.5*gv%ke*gv%Angstrom_m)

2786

2787 call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &

2788 "The tolerance for barotropic velocity discrepancies "//&

2789 "between the barotropic solution and the sum of the "//&

2790 "layer thicknesses.", units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)

2791 ! The speed of light is the default.

2792

2793 call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&

2794 "If true, allow the adjusted velocities to have a "//&

2795 "relative CFL change up to 0.5.", default=.false.)

2796 cs%vol_CFL = cs%aggress_adjust

2797 call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &

2798 "If true, use the ratio of the open face lengths to the "//&

2799 "tracer cell areas when estimating CFL numbers. The "//&

2800 "default is set by CONT_PPM_AGGRESS_ADJUST.", &

2801 default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)

2802 call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &

2803 "The maximum CFL of the adjusted velocities.", units="nondim", &

2804 default=0.5)

2805 call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", cs%better_iter, &

2806 "If true, stop corrective iterations using a velocity "//&

2807 "based criterion and only stop if the iteration is "//&

2808 "better than all predecessors.", default=.true.)

2809 call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", cs%use_visc_rem_max, &

2810 "If true, use more appropriate limiting bounds for "//&

2811 "corrections in strongly viscous columns.", default=.true.)

2812 call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &

2813 "If true, use the marginal face areas from the continuity "//&

2814 "solver for use as the weights in the barotropic solver. "//&

2815 "Otherwise use the transport averaged areas.", default=.true.)

2816 call get_param(param_file, mdl, "CONT_USE_H_MARG_MIN", use_h_marg_min, &

2817 "If true, the marginal thickness used and returned from continuity "//&

2818 "is bounded from below by a sub-roundoff value. Otherwise the "//&

2819 "minimum is 0.", default=.false.)

2820 cs%diag => diag

2821

2822 id_clock_reconstruct = cpu_clock_id('(Ocean continuity reconstruction)', grain=clock_routine)

2823 id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=clock_routine)

2824 id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=clock_routine)

2825

2826 if (use_h_marg_min) then

2827 cs%h_marg_min = gv%H_subroundoff

2828 else

2829 cs%h_marg_min = 0.

2830 endif

2831

2832 if (local_open_bc) then

2833 do n=1, obc%number_of_segments

2834 segment => obc%segment(n)

2835 if (associated(segment%h_Reg)) then

2836 if (.not. allocated(segment%h_Reg%h_res)) then

2837 write(mesg,'("In MOM_continuity_PPM, continuity_PPM_init called with ", &

2838 & "badly configured h_res.")')

2839 call mom_error(fatal, mesg)

2840 endif

2841 endif

2842 enddo

2843 endif

2844

2845end subroutine continuity_ppm_init

2846

2847!> continuity_PPM_stencil returns the continuity solver stencil size

2848function continuity_ppm_stencil(CS) result(stencil)

2849 type(continuity_ppm_cs), intent(in) :: cs !< Module's control structure.

2850 integer :: stencil !< The continuity solver stencil size with the current settings.

2851

2852 stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1

2853

2854end function continuity_ppm_stencil

2855

2856!> Set up a structure that stores the sizes of the i- and j-loops to to work on in the continuity solver.

2857function set_continuity_loop_bounds(G, CS, i_stencil, j_stencil) result(LB)

2858 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.

2859 type(continuity_ppm_cs), intent(in) :: cs !< Module's control structure.

2860 logical, optional, intent(in) :: i_stencil !< If present and true, extend the i-loop bounds

2861 !! by the stencil width of the continuity scheme.

2862 logical, optional, intent(in) :: j_stencil !< If present and true, extend the j-loop bounds

2863 !! by the stencil width of the continuity scheme.

2864 type(cont_loop_bounds_type) :: lb !< A type storing the array sizes to work on in the continuity routines.

2865

2866 ! Local variables

2867 logical :: add_i_stencil, add_j_stencil ! Local variables set based on i_stencil and j_stensil

2868 integer :: stencil ! The continuity solver stencil size with the current continuity scheme.

2869

2870 add_i_stencil = .false. ; if (present(i_stencil)) add_i_stencil = i_stencil

2871 add_j_stencil = .false. ; if (present(j_stencil)) add_j_stencil = j_stencil

2872

2873 stencil = continuity_ppm_stencil(cs)

2874

2875 if (add_i_stencil) then

2876 lb%ish = g%isc-stencil ; lb%ieh = g%iec+stencil

2877 else

2878 lb%ish = g%isc ; lb%ieh = g%iec

2879 endif

2880

2881 if (add_j_stencil) then

2882 lb%jsh = g%jsc-stencil ; lb%jeh = g%jec+stencil

2883 else

2884 lb%jsh = g%jsc ; lb%jeh = g%jec

2885 endif

2886

2887end function set_continuity_loop_bounds

2888

2889!> \namespace mom_continuity_ppm

2890!!

2891!! This module contains the subroutines that advect layer

2892!! thickness. The scheme here uses a Piecewise-Parabolic method with

2893!! a positive definite limiter.

2894

2895end module mom_continuity_ppm