src/core/MOM_continuity_PPM.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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		#include "do_concurrent_compat.h"
6
7		!> Solve the layer continuity equation using the PPM method for layer fluxes.
8		module MOM_continuity_PPM
9
10		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
11		use MOM_diag_mediator, only : time_type, diag_ctrl
12		use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
13		use MOM_file_parser, only : get_param, log_version, param_file_type
14		use MOM_grid, only : ocean_grid_type
15		use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE
16		use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
17		use MOM_unit_scaling, only : unit_scale_type
18		use MOM_variables, only : BT_cont_type, porous_barrier_type
19		use MOM_verticalGrid, only : verticalGrid_type
20
21		!$ use omp_lib, only: omp_get_num_devices
22
23		implicit none ; private
24
25		#include <MOM_memory.h>
26
27		public continuity_PPM, continuity_PPM_init, continuity_PPM_stencil
28		public continuity_fluxes, continuity_adjust_vel
29		public zonal_mass_flux, meridional_mass_flux
30		public zonal_edge_thickness, meridional_edge_thickness
31		public continuity_zonal_convergence, continuity_merdional_convergence
32		public zonal_flux_thickness, meridional_flux_thickness
33		public zonal_BT_mass_flux, meridional_BT_mass_flux
34		public set_continuity_loop_bounds
35
36		!>@{ CPU time clock IDs
37		integer :: id_clock_reconstruct, id_clock_update, id_clock_correct
38		!>@}
39
40		!> Control structure for mom_continuity_ppm
41		type, public :: continuity_PPM_CS ; private
42		logical :: initialized = .false. !< True if this control structure has been initialized.
43		type(diag_ctrl), pointer :: diag !< Diagnostics control structure.
44		logical :: upwind_1st !< If true, use a first-order upwind scheme.
45		logical :: monotonic !< If true, use the Colella & Woodward monotonic
46		!! limiter; otherwise use a simple positive
47		!! definite limiter.
48		logical :: simple_2nd !< If true, use a simple second order (arithmetic
49		!! mean) interpolation of the edge values instead
50		!! of the higher order interpolation.
51		real :: tol_eta !< The tolerance for free-surface height
52		!! discrepancies between the barotropic solution and
53		!! the sum of the layer thicknesses [H ~> m or kg m-2].
54		real :: tol_vel !< The tolerance for barotropic velocity
55		!! discrepancies between the barotropic solution and
56		!! the sum of the layer thicknesses [L T-1 ~> m s-1].
57		real :: CFL_limit_adjust !< The maximum CFL of the adjusted velocities [nondim]
58		real :: h_marg_min !< Negligible floor on h_marg, the marginal thickness
59		!! used to calculate the partial derivative of transports
60		!! with velocities [H ~> m or kg m-2]
61		logical :: aggress_adjust !< If true, allow the adjusted velocities to have a
62		!! relative CFL change up to 0.5. False by default.
63		logical :: vol_CFL !< If true, use the ratio of the open face lengths
64		!! to the tracer cell areas when estimating CFL
65		!! numbers. Without aggress_adjust, the default is
66		!! false; it is always true with.
67		logical :: better_iter !< If true, stop corrective iterations using a
68		!! velocity-based criterion and only stop if the
69		!! iteration is better than all predecessors.
70		logical :: use_visc_rem_max !< If true, use more appropriate limiting bounds
71		!! for corrections in strongly viscous columns.
72		logical :: marginal_faces !< If true, use the marginal face areas from the
73		!! continuity solver for use as the weights in the
74		!! barotropic solver. Otherwise use the transport
75		!! averaged areas.
76		integer :: niblock !< The i block size used in array calculations [nondim].
77		integer :: njblock !< The j block size used in array calculations [nondim].
78		integer :: nkblock !< The k block size used in reconstruction calculations [nondim].
79
80		end type continuity_PPM_CS
81
82		!> A container for loop bounds
83		type, public :: cont_loop_bounds_type ; private
84		!>@{ Loop bounds
85		integer :: ish, ieh, jsh, jeh
86		!>@}
87		end type cont_loop_bounds_type
88
89		!> Finds the thickness fluxes from the continuity solver or their vertical sum without
90		!! actually updating the layer thicknesses.
91		interface continuity_fluxes
92		module procedure continuity_3d_fluxes, continuity_2d_fluxes
93		end interface continuity_fluxes
94
95		contains
96
97		!> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme,
98		!! based on Lin (1994).
99	96	subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &
100	240	visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont, du_cor, dv_cor)
101		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
102		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
103		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
104		intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
105		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
106		intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
107		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
108		intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].
109		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
110		intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2].
111		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
112		intent(out) :: uh !< Zonal volume flux, uhdy [H L2 T-1 ~> m3 s-1 or kg s-1].
113		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
114		intent(out) :: vh !< Meridional volume flux, vhdx [H L2 T-1 ~> m3 s-1 or kg s-1].
115		real, intent(in) :: dt !< Time increment [T ~> s].
116		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
117		type(continuity_PPM_CS), intent(in) :: CS !< Module's control structure.
118		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
119		type(porous_barrier_type), intent(in) :: pbv !< pointers to porous barrier fractional cell metrics
120		real, dimension(SZIB_(G),SZJ_(G)), &
121		optional, intent(in) :: uhbt !< The summed volume flux through zonal faces
122		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
123		real, dimension(SZI_(G),SZJB_(G)), &
124		optional, intent(in) :: vhbt !< The summed volume flux through meridional faces
125		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
126		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
127		optional, intent(in) :: visc_rem_u
128		!< The fraction of zonal momentum originally
129		!! in a layer that remains after a time-step of viscosity, and the
130		!! fraction of a time-step's worth of a barotropic acceleration that
131		!! a layer experiences after viscosity is applied [nondim].
132		!! Visc_rem_u is between 0 (at the bottom) and 1 (far above the bottom).
133		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
134		optional, intent(in) :: visc_rem_v
135		!< The fraction of meridional momentum originally
136		!! in a layer that remains after a time-step of viscosity, and the
137		!! fraction of a time-step's worth of a barotropic acceleration that
138		!! a layer experiences after viscosity is applied [nondim].
139		!! Visc_rem_v is between 0 (at the bottom) and 1 (far above the bottom).
140		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
141		optional, intent(out) :: u_cor
142		!< The zonal velocities that give uhbt as the depth-integrated transport [L T-1 ~> m s-1].
143		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
144		optional, intent(out) :: v_cor
145		!< The meridional velocities that give vhbt as the depth-integrated
146		!! transport [L T-1 ~> m s-1].
147		type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
148		!! the effective open face areas as a function of barotropic flow.
149		real, dimension(SZIB_(G),SZJ_(G)), &
150		optional, intent(out) :: du_cor !< The zonal velocity increments from u that give uhbt
151		!! as the depth-integrated transports [L T-1 ~> m s-1].
152		real, dimension(SZI_(G),SZJB_(G)), &
153		optional, intent(out) :: dv_cor !< The meridional velocity increments from v that give vhbt
154		!! as the depth-integrated transports [L T-1 ~> m s-1].
155
156		! Local variables
157	146	real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
158	146	real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
159	146	real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
160	146	real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
161		real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
162		type(cont_loop_bounds_type) :: LB ! A type indicating the loop range for a phase of the updates
163		logical :: x_first
164		integer :: niblock !< i block size for array calculations [nondim].
165		integer :: njblock !< j block size for array calculations [nondim].
166
167	73	niblock = CS%niblock
168	73	njblock = CS%njblock
169
170	73	h_min = GV%Angstrom_H
171
172	73	if (.not.CS%initialized) call MOM_error(FATAL, &
173	0	"MOM_continuity_PPM: Module must be initialized before it is used.")
174
175	73	x_first = (MOD(G%first_direction,2) == 0)
176
177	73	if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, &
178		"MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither "// &
179	0	"one must be present in call to continuity_PPM.")
180
181		!$omp target enter data map(alloc: h_W, h_E, h_S, h_N)
182
183	73	if (x_first) then
184
185		! First advect zonally, with loop bounds that accomodate the subsequent meridional advection.
186	73	LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.true.)
187		! set whole-domain block sizes when ni/jblock is 0
188	73	if (niblock == 0) niblock = LB%ieh-LB%ish+2
189	73	if (njblock == 0) njblock = LB%jeh-LB%jsh+1
190
191	73	call zonal_edge_thickness(hin, h_W, h_E, G, GV, US, CS, OBC, LB)
192		call zonal_mass_flux(u, hin, h_W, h_E, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, &
193	73	niblock, njblock, LB, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)
194	73	call continuity_zonal_convergence(h, uh, dt, G, GV, LB, hin)
195
196		! Now advect meridionally, using the updated thicknesses to determine the fluxes.
197	73	LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
198
199	73	if (niblock == 0) niblock = LB%ieh-LB%ish+1
200	73	if (njblock == 0) njblock = LB%jeh-LB%jsh+2
201
202	73	call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC, LB)
203		call meridional_mass_flux(v, h, h_S, h_N, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, &
204	73	niblock, njblock, LB, vhbt, visc_rem_v, v_cor, BT_cont, dv_cor)
205	73	call continuity_merdional_convergence(h, vh, dt, G, GV, LB, hmin=h_min)
206
207		else ! .not. x_first
208
209		! First advect meridionally, with loop bounds that accomodate the subsequent zonal advection.
210	0	LB = set_continuity_loop_bounds(G, CS, i_stencil=.true., j_stencil=.false.)
211	0	if (niblock == 0) niblock = LB%ieh-LB%ish+1
212	0	if (njblock == 0) njblock = LB%jeh-LB%jsh+2
213	0	call meridional_edge_thickness(hin, h_S, h_N, G, GV, US, CS, OBC, LB)
214		call meridional_mass_flux(v, hin, h_S, h_N, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, &
215	0	niblock, njblock, LB, vhbt, visc_rem_v, v_cor, BT_cont, dv_cor)
216	0	call continuity_merdional_convergence(h, vh, dt, G, GV, LB, hin)
217
218		! Now advect zonally, using the updated thicknesses to determine the fluxes.
219	0	LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
220	0	if (niblock == 0) niblock = LB%ieh-LB%ish+2
221	0	if (njblock == 0) njblock = LB%jeh-LB%jsh+1
222	0	call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC, LB)
223		call zonal_mass_flux(u, h, h_W, h_E, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, &
224	0	niblock, njblock, LB, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)
225	0	call continuity_zonal_convergence(h, uh, dt, G, GV, LB, hmin=h_min)
226		endif
227
228		!$omp target exit data map(delete: h_W, h_E, h_S, h_N)
229
230	409	end subroutine continuity_PPM
231
232		!> Finds the thickness fluxes from the continuity solver without actually updating the
233		!! layer thicknesses. Because the fluxes in the two directions are calculated based on the
234		!! input thicknesses, which are not updated between the direcitons, the fluxes returned here
235		!! are not the same as those that would be returned by a call to continuity.
236	0	subroutine continuity_3d_fluxes(u, v, h, uh, vh, dt, G, GV, US, CS, OBC, pbv)
237		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.
238		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
239		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
240		intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
241		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
242		intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
243		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
244		intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
245		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
246		intent(out) :: uh !< Thickness fluxes through zonal faces,
247		!! uhdy [H L2 T-1 ~> m3 s-1 or kg s-1].
248		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
249		intent(out) :: vh !< Thickness fluxes through meridional faces,
250		!! vhdx [H L2 T-1 ~> m3 s-1 or kg s-1].
251		real, intent(in) :: dt !< Time increment [T ~> s].
252		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
253		type(continuity_PPM_CS), intent(in) :: CS !< Control structure for mom_continuity.
254		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
255		type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
256
257		! Local variables
258	0	real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
259	0	real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
260	0	real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
261	0	real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
262		integer :: niblock !< i block size for array calculations [nondim].
263		integer :: njblock !< j block size for array calculations [nondim].
264
265	0	niblock = CS%niblock
266	0	njblock = CS%njblock
267
268	0	call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC)
269	0	if (niblock == 0) niblock = G%iec-G%isc+2
270	0	if (njblock == 0) njblock = G%jec-G%jsc+1
271		call zonal_mass_flux(u, h, h_W, h_E, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, &
272	0	niblock=niblock, njblock=njblock)
273
274	0	call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC)
275	0	if (niblock == 0) niblock = G%iec-G%isc+1
276	0	if (njblock == 0) njblock = G%jec-G%jsc+2
277		call meridional_mass_flux(v, h, h_S, h_N, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, &
278	0	niblock=niblock, njblock=njblock)
279
280	0	end subroutine continuity_3d_fluxes
281
282		!> Find the vertical sum of the thickness fluxes from the continuity solver without actually
283		!! updating the layer thicknesses. Because the fluxes in the two directions are calculated
284		!! based on the input thicknesses, which are not updated between the directions, the fluxes
285		!! returned here are not the same as those that would be returned by a call to continuity.
286	0	subroutine continuity_2d_fluxes(u, v, h, uhbt, vhbt, dt, G, GV, US, CS, OBC, pbv)
287		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.
288		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
289		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
290		intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
291		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
292		intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
293		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
294		intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
295		real, dimension(SZIB_(G),SZJ_(G)), &
296		intent(out) :: uhbt !< Vertically summed thickness flux through
297		!! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
298		real, dimension(SZI_(G),SZJB_(G)), &
299		intent(out) :: vhbt !< Vertically summed thickness flux through
300		!! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
301		real, intent(in) :: dt !< Time increment [T ~> s].
302		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
303		type(continuity_PPM_CS), intent(in) :: CS !< Control structure for mom_continuity.
304		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
305		type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
306
307		! Local variables
308	0	real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
309	0	real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
310	0	real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
311	0	real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
312
313	0	call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC)
314	0	call zonal_BT_mass_flux(u, h, h_W, h_E, uhbt, dt, G, GV, US, CS, OBC, pbv%por_face_areaU)
315
316	0	call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC)
317	0	call meridional_BT_mass_flux(v, h, h_S, h_N, vhbt, dt, G, GV, US, CS, OBC, pbv%por_face_areaV)
318
319	0	end subroutine continuity_2d_fluxes
320
321		!> Correct the velocities to give the specified depth-integrated transports by applying a
322		!! barotropic acceleration (subject to viscous drag) to the velocities.
323	0	subroutine continuity_adjust_vel(u, v, h, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, visc_rem_u, visc_rem_v)
324		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.
325		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
326		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
327		intent(inout) :: u !< Zonal velocity, which will be adjusted to
328		!! give uhbt as the depth-integrated
329		!! transport [L T-1 ~> m s-1].
330		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
331		intent(inout) :: v !< Meridional velocity, which will be adjusted
332		!! to give vhbt as the depth-integrated
333		!! transport [L T-1 ~> m s-1].
334		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
335		intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
336		real, intent(in) :: dt !< Time increment [T ~> s].
337		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
338		type(continuity_PPM_CS), intent(in) :: CS !< Control structure for mom_continuity.
339		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
340		type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
341		real, dimension(SZIB_(G),SZJ_(G)), &
342		intent(in) :: uhbt !< The vertically summed thickness flux through
343		!! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
344		real, dimension(SZI_(G),SZJB_(G)), &
345		intent(in) :: vhbt !< The vertically summed thickness flux through
346		!! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
347		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
348		optional, intent(in) :: visc_rem_u !< Both the fraction of the zonal momentum
349		!! that remains after a time-step of viscosity, and
350		!! the fraction of a time-step's worth of a barotropic
351		!! acceleration that a layer experiences after viscosity
352		!! is applied [nondim]. This goes between 0 (at the
353		!! bottom) and 1 (far above the bottom). When this
354		!! column is under an ice shelf, this also goes to 0
355		!! at the top due to the no-slip boundary condition there.
356		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
357		optional, intent(in) :: visc_rem_v !< Both the fraction of the meridional momentum
358		!! that remains after a time-step of viscosity, and
359		!! the fraction of a time-step's worth of a barotropic
360		!! acceleration that a layer experiences after viscosity
361		!! is applied [nondim]. This goes between 0 (at the
362		!! bottom) and 1 (far above the bottom). When this
363		!! column is under an ice shelf, this also goes to 0
364		!! at the top due to the no-slip boundary condition there.
365
366		! Local variables
367	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_in !< Input zonal velocity [L T-1 ~> m s-1]
368	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_in !< Input meridional velocity [L T-1 ~> m s-1]
369	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uh !< Volume flux through zonal faces =
370		!! uhdy [H L2 T-1 ~> m3 s-1 or kg s-1].
371	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vh !< Volume flux through meridional faces =
372		!! vhdx [H L2 T-1 ~> m3 s-1 or kg s-1].
373	0	real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
374	0	real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
375	0	real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
376	0	real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
377		integer :: niblock !< i block size for array calculations [nondim].
378		integer :: njblock !< j block size for array calculations [nondim].
379
380	0	niblock = CS%niblock
381	0	njblock = CS%njblock
382
383		! It might not be necessary to separate the input velocity array from the adjusted velocities,
384		! but it seems safer to do so, even if it might be less efficient.
385	0	u_in(:,:,:) = u(:,:,:)
386	0	v_in(:,:,:) = v(:,:,:)
387
388	0	call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC)
389	0	if (niblock == 0) niblock = G%iec-G%isc+2
390	0	if (njblock == 0) njblock = G%jec-G%jsc+1
391		call zonal_mass_flux(u_in, h, h_W, h_E, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, &
392		niblock=niblock, njblock=njblock, &
393	0	uhbt=uhbt, visc_rem_u=visc_rem_u, u_cor=u)
394
395	0	call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC)
396	0	if (niblock == 0) niblock = G%iec-G%isc+1
397	0	if (njblock == 0) njblock = G%jec-G%jsc+2
398		call meridional_mass_flux(v_in, h, h_S, h_N, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, &
399		niblock=niblock, njblock=njblock, &
400	0	vhbt=vhbt, visc_rem_v=visc_rem_v, v_cor=v)
401
402	0	end subroutine continuity_adjust_vel
403
404
405		!> Updates the thicknesses due to zonal thickness fluxes.
406	73	subroutine continuity_zonal_convergence(h, uh, dt, G, GV, LB, hin, hmin)
407		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure
408		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure
409		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
410		intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2]
411		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
412		intent(in) :: uh !< Zonal thickness flux, uhdy [H L2 T-1 ~> m3 s-1 or kg s-1]
413		real, intent(in) :: dt !< Time increment [T ~> s]
414		type(cont_loop_bounds_type), intent(in) :: LB !< Loop bounds structure
415		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
416		optional, intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].
417		!! If hin is absent, h is also the initial thickness.
418		real, optional, intent(in) :: hmin !< The minimum layer thickness [H ~> m or kg m-2]
419
420		real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
421		integer :: i, j, k, ish, ieh, jsh, jeh, nz
422
423	73	call cpu_clock_begin(id_clock_update)
424
425	73	h_min = 0.0 ; if (present(hmin)) h_min = hmin
426
427	73	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
428
429	73	if (present(hin)) then
430	73	do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
431	43948993	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 )
432		enddo
433		else
434		! untested
435	0	do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
436	0	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 )
437		enddo
438		endif
439
440	73	call cpu_clock_end(id_clock_update)
441
442	73	end subroutine continuity_zonal_convergence
443
444		!> Updates the thicknesses due to meridional thickness fluxes.
445	0	subroutine continuity_merdional_convergence(h, vh, dt, G, GV, LB, hin, hmin)
446		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure
447		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure
448		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
449		intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2]
450		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
451		intent(in) :: vh !< Meridional thickness flux, vhdx [H L2 T-1 ~> m3 s-1 or kg s-1]
452		real, intent(in) :: dt !< Time increment [T ~> s]
453		type(cont_loop_bounds_type), intent(in) :: LB !< Loop bounds structure
454		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
455		optional, intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].
456		!! If hin is absent, h is also the initial thickness.
457		real, optional, intent(in) :: hmin !< The minimum layer thickness [H ~> m or kg m-2]
458
459		real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
460		integer :: i, j, k, ish, ieh, jsh, jeh, nz
461
462	73	call cpu_clock_begin(id_clock_update)
463
464	73	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
465
466	73	h_min = 0.0 ; if (present(hmin)) h_min = hmin
467
468	73	if (present(hin)) then
469		! untested
470	0	do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
471	0	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 )
472		enddo
473		else
474	73	do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
475	39954433	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 )
476		enddo
477		endif
478
479	73	call cpu_clock_end(id_clock_update)
480
481	73	end subroutine continuity_merdional_convergence
482
483
484		!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.
485	73	subroutine zonal_edge_thickness(h_in, h_W, h_E, G, GV, US, CS, OBC, LB_in)
486		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
487		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
488		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
489		intent(in) :: h_in !< Tracer cell layer thickness [H ~> m or kg m-2].
490		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
491		intent(out) :: h_W !< Western edge layer thickness [H ~> m or kg m-2].
492		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
493		intent(out) :: h_E !< Eastern edge layer thickness [H ~> m or kg m-2].
494		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
495		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
496		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
497		type(cont_loop_bounds_type), &
498		optional, intent(in) :: LB_in !< Loop bounds structure.
499
500		! Local variables
501		type(cont_loop_bounds_type) :: LB
502		integer :: i, j, k, ish, ieh, jsh, jeh, nz, nkblock
503
504	73	call cpu_clock_begin(id_clock_reconstruct)
505
506	73	if (present(LB_in)) then
507	73	LB = LB_in
508		else
509	0	LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
510		endif
511	73	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
512	73	nkblock = CS%nkblock
513	73	if (nkblock == 0) nkblock = nz
514
515	73	if (CS%upwind_1st) then
516	0	do concurrent (k=1:nz, j=jsh:jeh, i=ish-1:ieh+1)
517	0	h_W(i,j,k) = h_in(i,j,k) ; h_E(i,j,k) = h_in(i,j,k)
518		enddo
519		else
520		call PPM_reconstruction_x(h_in, h_W, h_E, G, GV, LB, &
521	73	nkblock, 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC)
522		endif
523
524	73	call cpu_clock_end(id_clock_reconstruct)
525
526	73	end subroutine zonal_edge_thickness
527
528
529		!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.
530	73	subroutine meridional_edge_thickness(h_in, h_S, h_N, G, GV, US, CS, OBC, LB_in)
531		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
532		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
533		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
534		intent(in) :: h_in !< Tracer cell layer thickness [H ~> m or kg m-2].
535		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
536		intent(out) :: h_S !< Southern edge layer thickness [H ~> m or kg m-2].
537		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
538		intent(out) :: h_N !< Northern edge layer thickness [H ~> m or kg m-2].
539		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
540		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
541		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
542		type(cont_loop_bounds_type), &
543		optional, intent(in) :: LB_in !< Loop bounds structure.
544
545		! Local variables
546		type(cont_loop_bounds_type) :: LB
547		integer :: i, j, k, ish, ieh, jsh, jeh, nz, nkblock
548
549	73	call cpu_clock_begin(id_clock_reconstruct)
550
551	73	if (present(LB_in)) then
552	73	LB = LB_in
553		else
554	0	LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
555		endif
556	73	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
557	73	nkblock = CS%nkblock
558	73	if (nkblock == 0) nkblock = nz
559
560	73	if (CS%upwind_1st) then
561		! untested
562	0	do concurrent (k=1:nz, j=jsh-1:jeh+1, i=ish:ieh)
563	0	h_S(i,j,k) = h_in(i,j,k) ; h_N(i,j,k) = h_in(i,j,k)
564		enddo
565		else
566		call PPM_reconstruction_y(h_in, h_S, h_N, G, GV, LB, &
567	73	nkblock, 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC)
568		endif
569
570	73	call cpu_clock_end(id_clock_reconstruct)
571
572	73	end subroutine meridional_edge_thickness
573
574
575		!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.
576	219	subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_face_areaU, &
577	168	niblock, njblock, LB_in, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)
578		integer, intent(in) :: niblock !< i block size for array calculations [nondim].
579		integer, intent(in) :: njblock !< j block size for array calculations [nondim].
580		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
581		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
582		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
583		intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
584		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
585		intent(in) :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
586		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
587		intent(in) :: h_W !< Western edge thicknesses [H ~> m or kg m-2].
588		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
589		intent(in) :: h_E !< Eastern edge thicknesses [H ~> m or kg m-2].
590		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
591		intent(out) :: uh !< Volume flux through zonal faces = uhdy
592		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
593		real, intent(in) :: dt !< Time increment [T ~> s].
594		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
595		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
596		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
597		real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &
598		intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]
599		type(cont_loop_bounds_type), &
600		optional, intent(in) :: LB_in !< Loop bounds structure.
601		real, dimension(SZIB_(G),SZJ_(G)), &
602		optional, intent(in) :: uhbt !< The summed volume flux through zonal faces
603		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
604		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
605		optional, intent(in) :: visc_rem_u
606		!< The fraction of zonal momentum originally in a layer that remains after a
607		!! time-step of viscosity, and the fraction of a time-step's worth of a barotropic
608		!! acceleration that a layer experiences after viscosity is applied [nondim].
609		!! Visc_rem_u is between 0 (at the bottom) and 1 (far above the bottom).
610		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
611		optional, intent(out) :: u_cor
612		!< The zonal velocities (u with a barotropic correction)
613		!! that give uhbt as the depth-integrated transport [L T-1 ~> m s-1]
614		type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe the
615		!! effective open face areas as a function of barotropic flow.
616		real, dimension(SZIB_(G),SZJ_(G)), &
617		optional, intent(out) :: du_cor !< The zonal velocity increments from u that give uhbt
618		!! as the depth-integrated transports [L T-1 ~> m s-1].
619
620		! Local variables
621		real, dimension(niblock,njblock,SZK_(GV)) :: &
622	146	duhdu, & ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
623	146	visc_rem ! A 2-D copy of visc_rem_u or an array of 1's [nondim].
624		real, dimension(niblock,njblock) :: &
625	146	du, & ! Corrective barotropic change in the velocity to give uhbt [L T-1 ~> m s-1].
626	146	du_min_CFL, & ! Lower limit on du correction to avoid CFL violations [L T-1 ~> m s-1]
627	146	du_max_CFL, & ! Upper limit on du correction to avoid CFL violations [L T-1 ~> m s-1]
628	146	duhdu_tot_0, & ! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
629	146	uh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
630	146	visc_rem_max ! The column maximum of visc_rem [nondim].
631	146	logical, dimension(niblock,njblock) :: do_I
632		real, dimension(niblock,njblock,SZK_(GV)) :: &
633	146	uh_t ! A block-sized copy of uh [H L2 T-1 ~> m3 s-1 or kg s-1].
634		real, dimension(niblock,njblock) :: &
635	146	uhbt_t ! A block-sized copy of uhbt [H L2 T-1 ~> m3 s-1 or kg s-1].
636	146	real, dimension(niblock,njblock) :: FAuI ! A list of sums of zonal face areas [H L ~> m2 or kg m-1].
637		real :: FA_u ! A sum of zonal face areas [H L ~> m2 or kg m-1].
638		real :: I_vrm ! 1.0 / visc_rem_max [nondim]
639		real :: CFL_dt ! The maximum CFL ratio of the adjusted velocities divided by
640		! the time step [T-1 ~> s-1].
641		real :: I_dt ! 1.0 / dt [T-1 ~> s-1].
642		real :: du_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].
643		real :: dx_E, dx_W ! Effective x-grid spacings to the east and west [L ~> m].
644		type(cont_loop_bounds_type) :: LB
645		integer :: i, j, k, ish, ieh, jsh, jeh, n, nz, ii, jj
646		integer :: l_seg ! The OBC segment number
647		logical :: use_visc_rem, set_BT_cont
648		logical :: local_specified_BC, local_Flather_OBC, local_open_BC, any_simple_OBC ! OBC-related logicals
649	146	logical :: simple_OBC_pt(niblock,njblock) ! Indicates points in a row with specified transport OBCs
650		integer :: j_start !< The j-index of the start of the current tile [nondim].
651		integer :: i_start !< The i-index of the start of the current tile [nondim].
652		integer :: j_end !< The j-index of the end of the current tile [nondim].
653		integer :: i_end !< The i-index of the end of the current tile [nondim].
654		integer :: nteams
655
656	73	call cpu_clock_begin(id_clock_correct)
657
658		!$omp target enter data &
659		!$omp map(alloc:uhbt_t,uh_t,duhdu,du,du_min_CFL,du_max_CFL,duhdu_tot_0,uh_tot_0, &
660		!$omp visc_rem_max,do_I,visc_rem,simple_OBC_pt)
661
662	73	use_visc_rem = present(visc_rem_u)
663
664	73	set_BT_cont = .false. ; if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
665
666	73	local_specified_BC = .false. ; local_Flather_OBC = .false. ; local_open_BC = .false.
667	73	if (associated(OBC)) then ; if (OBC%OBC_pe) then
668	0	local_specified_BC = OBC%specified_u_BCs_exist_globally
669	0	local_Flather_OBC = OBC%Flather_u_BCs_exist_globally
670	0	local_open_BC = OBC%open_u_BCs_exist_globally
671		endif ; endif
672
673	73	if (present(du_cor)) then
674	0	do concurrent (j=G%jsd:G%jed, I=G%IsdB:G%IedB)
675	0	du_cor(I,j) = 0.0
676		enddo
677		endif
678
679	73	if (present(LB_in)) then
680	73	LB = LB_in
681		else
682	0	LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
683		endif
684	73	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
685
686	73	CFL_dt = CS%CFL_limit_adjust / dt
687	73	I_dt = 1.0 / dt
688	73	if (CS%aggress_adjust) CFL_dt = I_dt
689
690	73	if (.not.use_visc_rem) then
691	1	do concurrent (k=1:nz, jj=1:njblock, ii=1:niblock)
692	9761	visc_rem(ii,jj,k) = 1.0
693		enddo
694		endif
695
696	7519	do j_start=jsh,jeh,njblock ; do i_start=ish-1,ieh,niblock
697	4964	i_end = min(i_start+niblock-1,ieh)
698	4964	j_end = min(j_start+njblock-1,jeh)
699
700		! calculate number of teams
701		!$ nteams = ceiling(real((j_end-j_start+1)*(i_end-i_start+1))/128.)
702
703	155125	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
704	738103	do_I(ii,jj) = .true.
705		enddo
706		! Set uh and duhdu.
707		!$omp target teams num_teams(nteams)
708	377264	do k=1,nz
709	372300	if (use_visc_rem) then
710		!$omp loop collapse(2) private(ii,jj)
711	44917200	do j=j_start,j_end ; do I=i_start,i_end
712	43124400	ii=I-i_start+1 ; jj=j-j_start+1
713	44550000	visc_rem(ii,jj,k) = visc_rem_u(I,j,k)
714		enddo ; enddo
715		endif
716		!$omp loop collapse(2) private(ii,jj)
717	45541050	do j=j_start,j_end ; do i=i_start,i_end
718	43723350	ii=i-i_start+1 ; jj=j-j_start+1
719		call flux_elem(u(i,j,k),h_in(i,j,k),h_in(i+1,j,k),h_W(i,j,k),h_W(i+1,j,k),h_E(i,j,k),&
720		h_E(i+1,j,k),uh_t(ii,jj,k),duhdu(ii,jj,k),visc_rem(ii,jj,k),G%dy_Cu(i,j),&
721		G%IareaT(i,j),G%IareaT(i+1,j),G%IdxT(i,j),G%IdxT(i+1,j),dt,CS%vol_CFL,&
722	43723350	por_face_areaU(I,j,k))
723	43723350	if (local_open_BC) &
724		call flux_elem_OBC(u(i,j,k),h_in(i,j,k),h_in(i+1,j,k),uh_t(ii,jj,k),duhdu(ii,jj,k),&
725		visc_rem(ii,jj,k),por_face_areaU(i,j,k),G%dy_Cu(i,j),OBC,&
726	1445400	OBC%segnum_u(i,j))
727		enddo ; enddo
728	377264	if (local_specified_BC) then
729		!$omp loop collapse(2) private(l_seg)
730	0	do j=j_start,j_end ; do i=i_start,i_end ; if (OBC%segnum_u(I,j) /= 0) then
731	0	l_seg = abs(OBC%segnum_u(I,j))
732	0	if (OBC%segment(l_seg)%specified) uh_t(I-i_start+1,j-j_start+1,k) = OBC%segment(l_seg)%normal_trans(I,j,k)
733		endif ; enddo ; enddo
734		endif
735		enddo
736		!$omp end target teams
737
738	4964	if (present(uhbt) .or. set_BT_cont) then
739	4896	if (use_visc_rem .and. CS%use_visc_rem_max) then
740		!$omp target
741		!$omp loop collapse(2) private(ii,jj)
742	598896	do j=j_start,j_end ; do i=i_start,i_end
743	574992	ii=i-i_start+1 ; jj=j-j_start+1
744	594000	visc_rem_max(ii,jj) = 0.0
745		enddo ; enddo
746	372096	do k=1,nz
747		!$omp loop collapse(2) private(ii,jj)
748	44922096	do j=j_start,j_end ; do i=i_start,i_end
749	43124400	ii=i-i_start+1 ; jj=j-j_start+1
750	44550000	visc_rem_max(ii,jj) = max(visc_rem_max(ii,jj), visc_rem(ii,jj,k))
751		enddo ; enddo
752		enddo
753		!$omp end target
754		else
755	0	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
756	0	visc_rem_max(ii,jj) = 1.0
757		enddo
758		endif
759		! Set limits on du that will keep the CFL number between -1 and 1.
760		! This should be adequate to keep the root bracketed in all cases.
761		!$omp target
762		!$omp loop collapse(2) private(ii,jj,I_vrm,dx_W,dx_E)
763	598896	do j=j_start,j_end ; do I=i_start,i_end
764	574992	ii=I-i_start+1 ; jj=j-j_start+1
765	574992	I_vrm = 0.0
766	574992	if (visc_rem_max(ii,jj) > 0.0) I_vrm = 1.0 / visc_rem_max(ii,jj)
767	574992	if (CS%vol_CFL) then
768	0	dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
769	0	dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
770	574992	else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
771	574992	du_max_CFL(ii,jj) = 2.0* (CFL_dt * dx_W) * I_vrm
772	574992	du_min_CFL(ii,jj) = -2.0 * (CFL_dt * dx_E) * I_vrm
773	594000	uh_tot_0(ii,jj) = 0.0 ; duhdu_tot_0(ii,jj) = 0.0
774		enddo ; enddo
775	372096	do k=1,nz
776		!$omp loop collapse(2) private(ii,jj)
777	44922096	do j=j_start,j_end ; do I=i_start,i_end
778	43124400	ii=I-i_start+1 ; jj=j-j_start+1
779	43124400	duhdu_tot_0(ii,jj) = duhdu_tot_0(ii,jj) + duhdu(ii,jj,k)
780	44550000	uh_tot_0(ii,jj) = uh_tot_0(ii,jj) + uh_t(ii,jj,k)
781		enddo
782		enddo ; enddo
783		!$omp end target
784	4896	if (use_visc_rem) then
785	4896	if (CS%aggress_adjust) then
786		!$omp target
787	0	do k=1,nz
788		!$omp loop collapse(2) private(ii,jj,dx_W,dx_E,du_lim)
789	0	do j=j_start,j_end ; do I=i_start,i_end
790	0	ii=I-i_start+1 ; jj=j-j_start+1
791	0	if (CS%vol_CFL) then
792	0	dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
793	0	dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
794	0	else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
795
796	0	du_lim = 0.499((dx_WI_dt - u(I,j,k)) + MIN(0.0,u(I-1,j,k)))
797	0	if (du_max_CFL(ii,jj) * visc_rem(ii,jj,k) > du_lim) &
798	0	du_max_CFL(ii,jj) = du_lim / visc_rem(ii,jj,k)
799
800	0	du_lim = 0.499((-dx_EI_dt - u(I,j,k)) + MAX(0.0,u(I+1,j,k)))
801	0	if (du_min_CFL(ii,jj) * visc_rem(ii,jj,k) < du_lim) &
802	0	du_min_CFL(ii,jj) = du_lim / visc_rem(ii,jj,k)
803		enddo ; enddo
804		enddo
805		!$omp end target
806		else
807		!$omp target
808	372096	do k=1,nz
809		!$omp loop collapse(2) private(ii,jj,dx_E,dx_W)
810	44922096	do j=j_start,j_end ; do I=i_start,i_end
811	43124400	ii=I-i_start+1 ; jj=j-j_start+1
812	43124400	if (CS%vol_CFL) then
813	0	dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
814	0	dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
815	43124400	else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
816
817	43124400	if (du_max_CFL(ii,jj) * visc_rem(ii,jj,k) > dx_WCFL_dt - u(I,j,k)G%mask2dCu(I,j)) &
818	359309	du_max_CFL(ii,jj) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(ii,jj,k)
819	43124400	if (du_min_CFL(ii,jj) * visc_rem(ii,jj,k) < -dx_ECFL_dt - u(I,j,k)G%mask2dCu(I,j)) &
820	9646385	du_min_CFL(ii,jj) = -(dx_E*CFL_dt + u(I,j,k)) / visc_rem(ii,jj,k)
821		enddo ; enddo
822		enddo
823		!$omp end target
824		endif
825		else
826	0	if (CS%aggress_adjust) then
827		!$omp target
828	0	do k=1,nz
829		!$omp loop collapse(2) private(ii,jj,dx_W,dx_E)
830	0	do j=j_start,j_end ; do I=i_start,i_end
831	0	ii=I-i_start+1 ; jj=j-j_start+1
832	0	if (CS%vol_CFL) then
833	0	dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
834	0	dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
835	0	else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
836
837		du_max_CFL(ii,jj) = MIN(du_max_CFL(ii,jj), 0.499 * &
838	0	((dx_W*I_dt - u(I,j,k)) + MIN(0.0,u(I-1,j,k))) )
839		du_min_CFL(ii,jj) = MAX(du_min_CFL(ii,jj), 0.499 * &
840	0	((-dx_E*I_dt - u(I,j,k)) + MAX(0.0,u(I+1,j,k))) )
841		enddo ; enddo
842		enddo
843		!$omp end target
844		else
845		!$omp target
846	0	do k=1,nz
847		!$omp loop collapse(2) private(ii,jj,dx_W,dx_E)
848	0	do j=j_start,j_end ; do I=i_start,i_end
849	0	ii=I-i_start+1 ; jj=j-j_start+1
850	0	if (CS%vol_CFL) then
851	0	dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
852	0	dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
853	0	else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
854
855	0	du_max_CFL(ii,jj) = MIN(du_max_CFL(ii,jj), dx_W*CFL_dt - u(I,j,k))
856	0	du_min_CFL(ii,jj) = MAX(du_min_CFL(ii,jj), -(dx_E*CFL_dt + u(I,j,k)))
857		enddo ; enddo
858		enddo
859		!$omp end target
860		endif
861		endif
862	153000	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
863	574992	du_max_CFL(ii,jj) = max(du_max_CFL(ii,jj),0.0)
864	727992	du_min_CFL(ii,jj) = min(du_min_CFL(ii,jj),0.0)
865		enddo
866
867	4896	any_simple_OBC = .false.
868	4896	if (present(uhbt) .or. set_BT_cont) then
869	4896	if (local_specified_BC .or. local_Flather_OBC) then
870	0	do concurrent (j=j_start:j_end, I=i_start:i_end) DO_LOCALITY(reduce(.or.:any_simple_OBC))
871	0	ii=I-i_start+1 ; jj=j-j_start+1
872	0	l_seg = abs(OBC%segnum_u(I,j))
873
874		! Avoid reconciling barotropic/baroclinic transports if transport is specified
875	0	simple_OBC_pt(ii,jj) = .false.
876	0	if (l_seg /= OBC_NONE) simple_OBC_pt(ii,jj) = OBC%segment(l_seg)%specified
877	0	do_I(ii,jj) = .not.simple_OBC_pt(ii,jj)
878	0	any_simple_OBC = any_simple_OBC .or. simple_OBC_pt(ii,jj)
879		enddo ; else
880	153000	do concurrent (jj=1:j_end-j_start+1, II=1:i_end-i_start+1)
881	727992	do_I(II,jj) = .true.
882		enddo ; endif
883		endif
884
885	4896	if (present(uhbt)) then
886	3264	do concurrent (j=j_start:j_end, I=i_start:i_end)
887	485328	uhbt_t(I-i_start+1,j-j_start+1) = uhbt(I,j)
888		enddo
889		! Find du and uh.
890		call zonal_flux_adjust(u, h_in, h_W, h_E, uhbt_t, uh_tot_0, duhdu_tot_0, du, &
891		du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
892		i_start, i_end, j_start, j_end, do_I, por_face_areaU, niblock, njblock, &
893	3264	uh_t, OBC=OBC)
894
895	3264	if (present(u_cor)) then
896	3264	do concurrent (k=1:nz, j=j_start:j_end, I=i_start:i_end)
897	28749600	ii=I-i_start+1 ; jj=j-j_start+1
898	28749600	u_cor(I,j,k) = u(I,j,k) + du(ii,jj) * visc_rem(ii,jj,k)
899	29234928	if (any_simple_OBC) then ; if (simple_OBC_pt(ii,jj)) then
900	0	u_cor(I,j,k) = OBC%segment(abs(OBC%segnum_u(I,j)))%normal_vel(I,j,k)
901		endif ; endif
902		enddo
903		endif ! u-corrected
904
905	3264	if (present(du_cor)) then
906	0	do concurrent (j=j_start:j_end, I=i_start:i_end)
907	0	ii=I-i_start+1 ; jj=j-j_start+1
908	0	du_cor(I,j) = du(ii,jj)
909		enddo
910		endif
911
912		endif
913
914	4896	if (set_BT_cont) then
915		call set_zonal_BT_cont(u, h_in, h_W, h_E, BT_cont, uh_tot_0, duhdu_tot_0,&
916		du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
917		visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
918	3264	por_face_areaU, niblock, njblock)
919	3264	if (any_simple_OBC) then
920		!$omp target
921		!$omp loop collapse(2) private(ii,jj)
922	0	do j=j_start,j_end ; do I=i_start,i_end
923	0	ii=i-i_start+1 ; jj=j-j_start+1
924	0	if (simple_OBC_pt(II,jj)) FAuI(II,jj) = GV%H_subroundoff*G%dy_Cu(I,j)
925		enddo ; enddo
926		! NOTE: simple_OBC_pt should prevent access to segment OBC_NONE
927	0	do k=1,nz
928		!$omp loop collapse(2) private(ii,jj,l_seg)
929	0	do j=j_start,j_end ; do I=i_start,i_end
930	0	ii=i-i_start+1 ; jj=j-j_start+1
931	0	if (simple_OBC_pt(II,jj)) then
932	0	l_seg = abs(OBC%segnum_u(I,j))
933	0	if ((abs(OBC%segment(l_seg)%normal_vel(I,j,k)) > 0.0) .and. &
934		(OBC%segment(l_seg)%specified)) &
935		FAuI(II,jj) = FAuI(II,jj) + &
936	0	OBC%segment(l_seg)%normal_trans(I,j,k) / OBC%segment(l_seg)%normal_vel(I,j,k)
937		endif
938		enddo ; enddo
939		enddo
940		!$omp loop collapse(2) private(ii,jj)
941	0	do j=j_start,j_end ; do I=i_start,i_end
942	0	ii=i-i_start+1 ; jj=j-j_start+1
943	0	if (simple_OBC_pt(ii,jj)) then
944	0	BT_cont%FA_u_W0(I,j) = FAuI(ii,jj) ; BT_cont%FA_u_E0(I,j) = FAuI(ii,jj)
945	0	BT_cont%FA_u_WW(I,j) = FAuI(ii,jj) ; BT_cont%FA_u_EE(I,j) = FAuI(ii,jj)
946	0	BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0
947		endif
948		enddo ; enddo
949		!$omp end target
950		endif
951		endif ! set_BT_cont
952
953		endif ! present(uhbt) or set_BT_cont
954
955	1241	do concurrent (k=1:nz, j=j_start:j_end, I=i_start:i_end)
956	44461453	uh(I,j,k) = uh_t(I-i_start+1,j-j_start+1,k)
957		enddo
958
959		enddo ; enddo ! ij block loop
960
961		!$omp target exit data &
962		!$omp map(release:uhbt_t,uh_t,duhdu,du,du_min_CFL,du_max_CFL,duhdu_tot_0,uh_tot_0,&
963		!$omp visc_rem_max,do_I,visc_rem,simple_OBC_pt)
964
965	73	if (local_open_BC .and. set_BT_cont) then
966	0	do n = 1, OBC%number_of_segments
967	0	if (OBC%segment(n)%open .and. OBC%segment(n)%is_E_or_W) then
968	0	I = OBC%segment(n)%HI%IsdB
969	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
970	0	do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed
971	0	FA_u = 0.0
972	0	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
973	0	BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u
974	0	BT_cont%FA_u_WW(I,j) = FA_u ; BT_cont%FA_u_EE(I,j) = FA_u
975	0	BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0
976		enddo
977		else
978	0	do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed
979	0	FA_u = 0.0
980	0	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
981	0	BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u
982	0	BT_cont%FA_u_WW(I,j) = FA_u ; BT_cont%FA_u_EE(I,j) = FA_u
983	0	BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0
984		enddo
985		endif
986		endif
987		enddo
988		endif
989
990	73	if (set_BT_cont) then ; if (allocated(BT_cont%h_u)) then
991	48	if (present(u_cor)) then
992		call zonal_flux_thickness(u_cor, h_in, h_W, h_E, BT_cont%h_u, dt, G, GV, US, LB, &
993	24	CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU, visc_rem_u)
994		else
995		call zonal_flux_thickness(u, h_in, h_W, h_E, BT_cont%h_u, dt, G, GV, US, LB, &
996	24	CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU, visc_rem_u)
997		endif
998		endif ; endif
999
1000	73	call cpu_clock_end(id_clock_correct)
1001
1002	193	end subroutine zonal_mass_flux
1003
1004
1005		!> Calculates the vertically integrated mass or volume fluxes through the zonal faces.
1006	0	subroutine zonal_BT_mass_flux(u, h_in, h_W, h_E, uhbt, dt, G, GV, US, CS, OBC, por_face_areaU, LB_in)
1007		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1008		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1009		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
1010		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to
1011		!! calculate fluxes [H ~> m or kg m-2]
1012		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_W !< Western edge thickness in the PPM
1013		!! reconstruction [H ~> m or kg m-2].
1014		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_E !< Eastern edge thickness in the PPM
1015		!! reconstruction [H ~> m or kg m-2].
1016		real, dimension(SZIB_(G),SZJ_(G)), intent(out) :: uhbt !< The summed volume flux through zonal
1017		!! faces [H L2 T-1 ~> m3 s-1 or kg s-1].
1018		real, intent(in) :: dt !< Time increment [T ~> s].
1019		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1020		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.G
1021		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition type
1022		!! specifies whether, where, and what
1023		!! open boundary conditions are used.
1024		real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]
1025		type(cont_loop_bounds_type), optional, intent(in) :: LB_in !< Loop bounds structure.
1026
1027		! Local variables
1028	0	real :: uh(SZIB_(G),SZJ_(G),SZK_(GV)) ! Volume flux through zonal faces = uhdy [H L2 T-1 ~> m3 s-1 or kg s-1]
1029	0	real :: duhdu(SZIB_(G),SZJ_(G),SZK_(GV)) ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
1030		integer :: i, j, k, ish, ieh, jsh, jeh, nz, l_seg
1031		logical :: local_specified_BC
1032	0	logical, dimension(SZJ_(G)) :: OBC_in_row
1033
1034	0	call cpu_clock_begin(id_clock_correct)
1035
1036	0	local_specified_BC = .false.
1037	0	if (associated(OBC)) then ; if (OBC%OBC_pe) then
1038	0	local_specified_BC = OBC%specified_v_BCs_exist_globally
1039		endif ; endif
1040
1041	0	if (present(LB_in)) then
1042	0	ish = LB_in%ish ; ieh = LB_in%ieh ; jsh = LB_in%jsh ; jeh = LB_in%jeh ; nz = GV%ke
1043		else
1044	0	ish = G%isc ; ieh = G%iec ; jsh = G%jsc ; jeh = G%jec ; nz = GV%ke
1045		endif
1046
1047	0	OBC_in_row(:) = .false.
1048
1049	0	uhbt(:,:) = 0.0
1050
1051		! Determining whether there are any OBC points outside of the k-loop should be more efficient.
1052	0	if (local_specified_BC) then
1053	0	do j=jsh,jeh ; do I=ish-1,ieh ; if (OBC%segnum_u(I,j) /= 0) then
1054	0	if (OBC%segment(abs(OBC%segnum_u(I,j)))%specified) OBC_in_row(j) = .true.
1055		endif ; enddo ; enddo
1056		endif
1057
1058		! This sets uh and duhdu.
1059	0	do concurrent (k=1:nz, j=jsh:jeh, I=ish-1:ieh)
1060		call flux_elem(u(I,j,k), h_in(I,j,k), h_in(I+1,j,k), h_W(I,j,k), h_W(I+1,j,k), h_E(I,j,k), &
1061		h_E(I+1,j,k), uh(I,j,k), duhdu(I,j,k), 1.0, G%dy_Cu(I,j), G%IareaT(I,j), &
1062		G%IareaT(I+1,j), G%IdxT(I,j), G%IdxT(I+1,j), dt, CS%vol_CFL, &
1063	0	por_face_areaU(I,j,k))
1064	0	if (local_specified_BC) &
1065		call flux_elem_OBC(u(I,j,k), h_in(I,j,k), h_in(I+1,j,k), uh(I,j,k), duhdu(I,j,k), 1.0, &
1066	0	por_face_areaU(I,j,k), G%dy_Cu(I,j), OBC, OBC%segnum_u(I,j))
1067		enddo
1068
1069	0	do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
1070	0	if (OBC_in_row(j) .and. OBC%segnum_u(I,j) /= 0) then
1071	0	l_seg = abs(OBC%segnum_u(I,j))
1072	0	if (OBC%segment(l_seg)%specified) uh(I,j,k) = OBC%segment(l_seg)%normal_trans(I,j,k)
1073		endif
1074		enddo ; enddo ; enddo
1075
1076		! Accumulate the barotropic transport.
1077	0	do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh
1078	0	uhbt(I,j) = uhbt(I,j) + uh(I,j,k)
1079		enddo ; enddo ; enddo ! j-loop
1080
1081	0	call cpu_clock_end(id_clock_correct)
1082
1083	0	end subroutine zonal_BT_mass_flux
1084
1085		!> Evaluates the zonal mass or volume fluxes in an element.
1086		!DIR$ ATTRIBUTES FORCEINLINE :: flux_elem
1087	343010625	elemental subroutine flux_elem(u, h, h_p1, h_L, h_L_p1, h_R, h_R_p1, uh, duhdu, visc_rem, &
1088		G_dy_Cu, G_IareaT, G_IareaT_p1, G_IdxT, G_IdxT_p1, dt, &
1089		vol_CFL, por_face_area)
1090		real, intent(in) :: u !< Zonal or meridional velocity [L T-1 ~> m s-1].
1091		real, intent(in) :: visc_rem !< Both the fraction of the
1092		!! momentum originally in a layer that remains after a time-step
1093		!! of viscosity, and the fraction of a time-step's worth of a barotropic
1094		!! acceleration that a layer experiences after viscosity is applied [nondim].
1095		!! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).
1096		real, intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
1097		real, intent(in) :: h_p1 !< Layer thickness - offset by 1 [ H ~> m or kg m-2].
1098		real, intent(in) :: h_L !< West/South edge thickness [H ~> m or kg m-2].
1099		real, intent(in) :: h_L_p1 !< West/South edge thickness - offset by 1 [H ~> m or kg m-2].
1100		real, intent(in) :: h_R !< East/North edge thickness [H ~> m or kg m-2].
1101		real, intent(in) :: h_R_p1 !< East/North edge thickness - offset by 1 [H ~> m or kg m-2].
1102		real, intent(out) :: uh !< Zonal or meridional mass or volume transport
1103		!! [H L2 T-1 ~> m3 s-1 or kg s-1].
1104		real, intent(out) :: duhdu !< Partial derivative of uh
1105		!! with u [H L ~> m2 or kg m-1].
1106		real, intent(in) :: dt !< Time increment [T ~> s]
1107		logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face areas to the
1108		!! cell areas when estimating the CFL number.
1109		real, intent(in) :: por_face_area !< fractional open area of U/V-faces [nondim].
1110		real, intent(in) :: G_dy_Cu !< The grid cell's unblocked lengths of the u/v-faces
1111		!! of the h-cell [L ~> m].
1112		real, intent(in) :: G_IareaT !< The grid cell's 1/areaT [L-2 ~> m-2].
1113		real, intent(in) :: G_IareaT_p1 !< The grid cell's 1/areaT - offset by 1 [L-2 ~> m-2].
1114		real, intent(in) :: G_IdxT !< The grid cell's 1/dxT [L-1 ~> m-1].
1115		real, intent(in) :: G_IdxT_p1 !< The grid cell's 1/dxT - offset by 1 [L-1 ~> m-1].
1116		! Local variables
1117		real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1118		real :: curv_3 ! A measure of the thickness curvature over a grid length [H ~> m or kg m-2]
1119		real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1120		real :: tmp ! temporary variable to store precalculted values
1121		real :: dh ! h differential between E/W
1122
1123		! Set new values of uh and duhdu.
1124	343010625	tmp = G_dy_Cu * por_face_area ! precalculate things
1125	343010625	if (u > 0.0) then
1126	87385524	if (vol_CFL) then ; CFL = (u * dt) * (G_dy_Cu * G_IareaT)
1127	87385524	else ; CFL = u * dt * G_IdxT ; endif
1128	87385524	curv_3 = (h_L + h_R) - 2.0*h
1129	87385524	dh = h_L - h_R
1130		uh = tmp * u * &
1131	87385524	(h_R + CFL * (0.5dh + curv_3(CFL - 1.5)))
1132	87385524	h_marg = h_R + CFL * (dh + 3.0curv_3(CFL - 1.0))
1133	255625101	elseif (u < 0.0) then
1134	164484051	if (vol_CFL) then ; CFL = (-u * dt) * (G_dy_Cu * G_IareaT_p1)
1135	164484051	else ; CFL = -u * dt * G_IdxT_p1 ; endif
1136	164484051	curv_3 = (h_L_p1 + h_R_p1) - 2.0*h_p1
1137	164484051	dh = h_R_p1-h_L_p1
1138		uh = tmp * u * &
1139	164484051	(h_L_p1 + CFL * (0.5dh + curv_3(CFL - 1.5)))
1140	164484051	h_marg = h_L_p1 + CFL * (dh + 3.0curv_3(CFL - 1.0))
1141		else
1142	91141050	uh = 0.0
1143	91141050	h_marg = 0.5 * (h_L_p1 + h_R)
1144		endif
1145	343010625	duhdu = tmp * h_marg * visc_rem
1146
1147	343010625	end subroutine flux_elem
1148
1149		!DIR$ ATTRIBUTES FORCEINLINE :: flux_elem_OBC
1150	0	elemental subroutine flux_elem_OBC(u, h, h_p1, uh, duhdu, visc_rem, por_face_area, &
1151		G_dy_Cu, OBC, l_seg)
1152		real, intent(in) :: u !< Zonal/meridional velocity [L T-1 ~> m s-1].
1153		real, intent(in) :: visc_rem !< Both the fraction of the
1154		!! momentum originally in a layer that remains after a time-step
1155		!! of viscosity, and the fraction of a time-step's worth of a barotropic
1156		!! acceleration that a layer experiences after viscosity is applied [nondim].
1157		!! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).
1158		real, intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
1159		real, intent(in) :: h_p1 !< Layer thickness offset by 1 [H ~> m or kg m-2].
1160		real, intent(inout) :: uh !< Zonal/meridional mass or volume
1161		!! transport [H L2 T-1 ~> m3 s-1 or kg s-1].
1162		real, intent(inout) :: duhdu !< Partial derivative of uh
1163		!! with u [H L ~> m2 or kg m-1].
1164		real, intent(in) :: por_face_area !< fractional open area of U/V-faces
1165		!! [nondim].
1166		real, intent(in) :: G_dy_Cu !< The grid cell's unblocked lengths of the
1167		!! u/v-faces of the h-cell [L ~> m].
1168		!! ratio of face areas to the cell areas when estimating the CFL number.
1169		type(ocean_OBC_type), intent(in) :: OBC !< Open boundaries control structure.
1170		integer, intent(in) :: l_seg !< Segment index.
1171
1172		! untested
1173	0	if (l_seg /= 0) then
1174	0	if (OBC%segment(abs(l_seg))%open) then
1175	0	if (l_seg > 0) then ! OBC_DIRECTION_E or OBC_DIRECTION_N
1176	0	uh = (G_dy_Cu * por_face_area) * u * h
1177	0	duhdu = (G_dy_Cu * por_face_area) * h * visc_rem
1178		else ! OBC_DIRECTION_W or OBC_DIRECTION_S
1179	0	uh = (G_dy_Cu * por_face_area) * u * h_p1
1180	0	duhdu = (G_dy_Cu* por_face_area) * h_p1 * visc_rem
1181		endif
1182		endif
1183		endif
1184
1185	0	end subroutine flux_elem_OBC
1186
1187
1188		!> Sets the effective interface thickness associated with the fluxes at each zonal velocity point,
1189		!! optionally scaling back these thicknesses to account for viscosity and fractional open areas.
1190	48	subroutine zonal_flux_thickness(u, h, h_W, h_E, h_u, dt, G, GV, US, LB, vol_CFL, &
1191	48	marginal, OBC, por_face_areaU, visc_rem_u)
1192		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1193		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1194		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
1195		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to
1196		!! calculate fluxes [H ~> m or kg m-2].
1197		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_W !< West edge thickness in the
1198		!! reconstruction [H ~> m or kg m-2].
1199		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_E !< East edge thickness in the
1200		!! reconstruction [H ~> m or kg m-2].
1201		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_u !< Effective thickness at zonal faces,
1202		!! scaled down to account for the effects of
1203		!! viscosity and the fractional open area
1204		!! [H ~> m or kg m-2].
1205		real, intent(in) :: dt !< Time increment [T ~> s].
1206		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1207		type(cont_loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1208		logical, intent(in) :: vol_CFL !< If true, rescale the ratio
1209		!! of face areas to the cell areas when estimating the CFL number.
1210		logical, intent(in) :: marginal !< If true, report the
1211		!! marginal face thicknesses; otherwise report transport-averaged thicknesses.
1212		real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &
1213		intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]
1214		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
1215		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1216		optional, intent(in) :: visc_rem_u
1217		!< Both the fraction of the momentum originally in a layer that remains after
1218		!! a time-step of viscosity, and the fraction of a time-step's worth of a
1219		!! barotropic acceleration that a layer experiences after viscosity is applied [nondim].
1220		!! Visc_rem_u is between 0 (at the bottom) and 1 (far above the bottom).
1221
1222		! Local variables
1223		real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1224		real :: curv_3 ! A measure of the thickness curvature over a grid length [H ~> m or kg m-2]
1225		logical :: local_open_BC
1226		integer :: i, j, k, ish, ieh, jsh, jeh, nz, n
1227		real :: dh
1228	48	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
1229
1230	48	do concurrent (k=1:nz, j=jsh:jeh, I=ish-1:ieh)
1231	28749600	if (u(I,j,k) > 0.0) then
1232	7951478	if (vol_CFL) then ; CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j))
1233	7951478	else ; CFL = u(I,j,k) * dt * G%IdxT(i,j) ; endif
1234	7951478	curv_3 = (h_W(i,j,k) + h_E(i,j,k)) - 2.0*h(i,j,k)
1235	7951478	dh = h_W(i,j,k) - h_E(i,j,k)
1236	7951478	if (marginal) then
1237	7951478	h_u(I,j,k) = h_E(i,j,k) + CFL * (dh + 3.0curv_3(CFL - 1.0))
1238		else
1239	0	h_u(I,j,k) = h_E(i,j,k) + CFL * (0.5dh + curv_3(CFL - 1.5))
1240		endif
1241	20798122	elseif (u(I,j,k) < 0.0) then
1242	9490222	if (vol_CFL) then ; CFL = (-u(I,j,k)dt) (G%dy_Cu(I,j) * G%IareaT(i+1,j))
1243	9490222	else ; CFL = -u(I,j,k) * dt * G%IdxT(i+1,j) ; endif
1244	9490222	curv_3 = (h_W(i+1,j,k) + h_E(i+1,j,k)) - 2.0*h(i+1,j,k)
1245	9490222	dh = h_E(i+1,j,k)-h_W(i+1,j,k)
1246	9490222	if (marginal) then
1247	9490222	h_u(I,j,k) = h_W(i+1,j,k) + CFL * (dh + 3.0curv_3(CFL - 1.0))
1248		else
1249	0	h_u(I,j,k) = h_W(i+1,j,k) + CFL * (0.5dh + curv_3(CFL - 1.5))
1250		endif
1251		else
1252		! The choice to use the arithmetic mean here is somewhat arbitrarily, but
1253		! it should be noted that h_W(i+1,j,k) and h_E(i,j,k) are usually the same.
1254	11307900	h_u(I,j,k) = 0.5 * (h_W(i+1,j,k) + h_E(i,j,k))
1255		! h_marg = (2.0 * h_W(i+1,j,k) * h_E(i,j,k)) / &
1256		! (h_W(i+1,j,k) + h_E(i,j,k) + GV%H_subroundoff)
1257		endif
1258
1259	57888384	if (present(visc_rem_u)) then
1260		! Scale back the thickness to account for the effects of viscosity and the fractional open
1261		! thickness to give an appropriate non-normalized weight for each layer in determining the
1262		! barotropic acceleration.
1263	28749600	h_u(I,j,k) = h_u(I,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
1264		else
1265	0	h_u(I,j,k) = h_u(I,j,k) * por_face_areaU(I,j,k)
1266		endif
1267		enddo
1268
1269	48	local_open_BC = .false.
1270	48	if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally
1271	48	if (local_open_BC) then
1272		! untested
1273	0	do n = 1, OBC%number_of_segments
1274	0	if (OBC%segment(n)%open .and. OBC%segment(n)%is_E_or_W) then
1275	0	I = OBC%segment(n)%HI%IsdB
1276	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
1277	0	if (present(visc_rem_u)) then
1278	0	do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
1279	0	h_u(I,j,k) = h(i,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
1280		enddo
1281		else
1282	0	do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
1283	0	h_u(I,j,k) = h(i,j,k) * por_face_areaU(I,j,k)
1284		enddo
1285		endif
1286		else
1287	0	if (present(visc_rem_u)) then
1288	0	do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
1289	0	h_u(I,j,k) = h(i+1,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
1290		enddo
1291		else
1292	0	do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
1293	0	h_u(I,j,k) = h(i+1,j,k) * por_face_areaU(I,j,k)
1294		enddo
1295		endif
1296		endif
1297		endif
1298		enddo
1299		endif
1300
1301	48	end subroutine zonal_flux_thickness
1302
1303		!> Returns the barotropic velocity adjustment that gives the
1304		!! desired barotropic (layer-summed) transport.
1305	6528	subroutine zonal_flux_adjust(u, h_in, h_W, h_E, uhbt, uh_tot_0, duhdu_tot_0, &
1306	6528	du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
1307	6528	i_start, i_end, j_start, j_end, do_I_in, por_face_areaU, niblock, njblock, &
1308	6528	uh_3d, OBC)
1309		integer, intent(in) :: niblock !< i block size for array calculations [nondim].
1310		integer, intent(in) :: njblock !< j block size for array calculations [nondim].
1311		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1312		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1313		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
1314		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to
1315		!! calculate fluxes [H ~> m or kg m-2].
1316		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_W !< West edge thickness in the
1317		!! reconstruction [H ~> m or kg m-2].
1318		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_E !< East edge thickness in the
1319		!! reconstruction [H ~> m or kg m-2].
1320		real, dimension(niblock,njblock,SZK_(GV)),intent(in) :: visc_rem !< Both the fraction of the
1321		!! momentum originally in a layer that remains after a time-step of viscosity, and
1322		!! the fraction of a time-step's worth of a barotropic acceleration that a layer
1323		!! experiences after viscosity is applied [nondim].
1324		!! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).
1325		real, dimension(niblock,njblock), intent(in) :: uhbt !< The summed volume flux
1326		!! through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
1327
1328		real, dimension(niblock,njblock), intent(in) :: du_max_CFL !< Maximum acceptable
1329		!! value of du [L T-1 ~> m s-1].
1330		real, dimension(niblock,njblock), intent(in) :: du_min_CFL !< Minimum acceptable
1331		!! value of du [L T-1 ~> m s-1].
1332		real, dimension(niblock,njblock), intent(in) :: uh_tot_0 !< The summed transport
1333		!! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
1334		real, dimension(niblock,njblock), intent(in) :: duhdu_tot_0 !< The partial derivative
1335		!! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
1336		real, dimension(niblock,njblock), intent(out) :: du !<
1337		!! The barotropic velocity adjustment [L T-1 ~> m s-1].
1338		real, intent(in) :: dt !< Time increment [T ~> s].
1339		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1340		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
1341		integer, intent(in) :: i_start !< Start of i index range [nondim].
1342		integer, intent(in) :: i_end !< End of i index range [nondim].
1343		integer, intent(in) :: j_start !< Start of j index range [nondim].
1344		integer, intent(in) :: j_end !< End of j index range [nondim].
1345		logical, dimension(niblock,njblock),intent(in) :: do_I_in !<
1346		!! A logical flag indicating which I values to work on.
1347		real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &
1348		intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]
1349		real, dimension(niblock,njblock,SZK_(GV)), intent(inout) :: uh_3d !<
1350		!! Volume flux through zonal faces = uhdy [H L2 T-1 ~> m3 s-1 or kg s-1].
1351		type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1352		! Local variables
1353		real :: &
1354		u_new, & ! The velocity with the correction added [L T-1 ~> m s-1].
1355		duhdu ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
1356		real, dimension(niblock,njblock) :: &
1357	13056	uh_err, & ! Difference between uhbt and the summed uh [H L2 T-1 ~> m3 s-1 or kg s-1].
1358	13056	uh_err_best, & ! The smallest value of uh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
1359	13056	duhdu_tot,&! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
1360	13056	du_min, & ! Lower limit on du correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
1361	13056	du_max ! Upper limit on du correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
1362		real :: du_prev ! The previous value of du [L T-1 ~> m s-1].
1363		real :: ddu ! The change in du from the previous iteration [L T-1 ~> m s-1].
1364		real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
1365		real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
1366		integer :: i, j, k, nz, itt, ii, jj
1367	13056	logical :: domore, do_I(niblock,njblock)
1368		logical :: local_open_BC ! True if there are open OBC points on this PE [nondim].
1369		integer, parameter :: max_itts = 20
1370
1371		!$omp target enter data &
1372		!$omp map(alloc: uh_err,uh_err_best,duhdu_tot,du_min,du_max,do_I)
1373
1374	6528	nz = GV%ke
1375	6528	local_open_BC = .false.
1376	6528	if (present(OBC)) then ; if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally ; endif
1377
1378	6528	tol_vel = CS%tol_vel
1379
1380		!$omp target
1381
1382		!$omp loop collapse(2) private(ii,jj)
1383	798528	do j=j_start,j_end ; do I=i_start,i_end
1384	766656	ii=I-i_start+1 ; jj=j-j_start+1
1385	766656	du(ii,jj) = 0.0 ; do_I(ii,jj) = do_I_in(ii,jj)
1386	766656	du_max(ii,jj) = du_max_CFL(ii,jj) ; du_min(ii,jj) = du_min_CFL(ii,jj)
1387	766656	uh_err(ii,jj) = uh_tot_0(ii,jj) - uhbt(ii,jj) ; duhdu_tot(ii,jj) = duhdu_tot_0(ii,jj)
1388	792000	uh_err_best(ii,jj) = abs(uh_err(ii,jj))
1389		enddo ; enddo
1390
1391	18450	do itt=1,max_itts
1392	6528	select case (itt)
1393	6528	case (:1) ; tol_eta = 1e-6 * CS%tol_eta
1394	6080	case (2) ; tol_eta = 1e-4 * CS%tol_eta
1395	5413	case (3) ; tol_eta = 1e-2 * CS%tol_eta
1396	18450	case default ; tol_eta = CS%tol_eta
1397		end select
1398
1399		!$omp loop collapse(2) private(ii,jj)
1400	2300098	do j=j_start,j_end ; do I=i_start,i_end
1401	2208616	ii=I-i_start+1 ; jj=j-j_start+1
1402	2281648	if (uh_err(ii,jj) > 0.0) then ; du_max(ii,jj) = du(ii,jj)
1403	1542761	elseif (uh_err(ii,jj) < 0.0) then ; du_min(ii,jj) = du(ii,jj)
1404	894347	else ; do_I(ii,jj) = .false. ; endif
1405		enddo ; enddo
1406		#ifndef __NVCOMPILER_OPENMP_GPU
1407	18450	domore = .false.
1408		#endif
1409		!$omp loop collapse(2) private(ii,jj,ddu,du_prev)
1410	2300098	do j=j_start,j_end ; do I=i_start,i_end
1411	2208616	ii=I-i_start+1 ; jj=j-j_start+1
1412	2281648	if (do_I(ii,jj)) then
1413	1141957	if ((dt * min(G%IareaT(i,j),G%IareaT(i+1,j))*abs(uh_err(ii,jj)) > tol_eta) .or. &
1414		(CS%better_iter .and. ((abs(uh_err(ii,jj)) > tol_vel * duhdu_tot(ii,jj)) .or. &
1415		(abs(uh_err(ii,jj)) > uh_err_best(ii,jj))) )) then
1416		! Use Newton's method, provided it stays bounded. Otherwise bisect
1417		! the value with the appropriate bound.
1418	671902	ddu = -uh_err(ii,jj) / duhdu_tot(ii,jj)
1419	671902	du_prev = du(ii,jj)
1420	671902	du(ii,jj) = du(ii,jj) + ddu
1421	671902	if (abs(ddu) < 1.0e-15*abs(du(ii,jj))) then
1422	0	do_I(ii,jj) = .false. ! ddu is small enough to quit.
1423	671902	elseif (ddu > 0.0) then
1424	201417	if (du(ii,jj) >= du_max(ii,jj)) then
1425	0	du(ii,jj) = 0.5*(du_prev + du_max(ii,jj))
1426	0	if (du_max(ii,jj) - du_prev < 1.0e-15*abs(du(ii,jj))) do_I(ii,jj) = .false.
1427		endif
1428		else ! ddu < 0.0
1429	470485	if (du(ii,jj) <= du_min(ii,jj)) then
1430	0	du(ii,jj) = 0.5*(du_prev + du_min(ii,jj))
1431	0	if (du_prev - du_min(ii,jj) < 1.0e-15*abs(du(ii,jj))) do_I(ii,jj) = .false.
1432		endif
1433		endif
1434		#ifndef __NVCOMPILER_OPENMP_GPU
1435	671902	if (do_I(ii,jj)) domore = .true.
1436		#endif
1437		else
1438	470055	do_I(ii,jj) = .false.
1439		endif
1440		endif
1441		enddo ; enddo
1442		#ifndef __NVCOMPILER_OPENMP_GPU
1443	18450	if (.not.domore) exit
1444		#endif
1445
1446		!$omp loop collapse(2) private(ii,jj)
1447	1501570	do j=j_start,j_end ; do I=i_start,i_end
1448	1441960	ii=I-i_start+1 ; jj=j-j_start+1
1449	1489648	uh_err(ii,jj) = -uhbt(ii,jj) ; duhdu_tot(ii,jj) = 0.0
1450		enddo ; enddo
1451	906072	do k=1,nz
1452		!$omp loop collapse(2) private(ii,jj,duhdu,u_new)
1453	112629672	do j=j_start,j_end ; do i=i_start,i_end
1454	108147000	ii=i-i_start+1 ; jj=j-j_start+1
1455	111723600	if (do_I(ii,jj)) then
1456	50392650	u_new = u(i,j,k) + du(ii,jj) * visc_rem(ii,jj,k)
1457		call flux_elem(u_new,h_in(i,j,k),h_in(i+1,j,k),h_W(i,j,k),h_W(i+1,j,k),h_E(i,j,k),&
1458		h_E(i+1,j,k),uh_3d(ii,jj,k),duhdu,visc_rem(ii,jj,k),G%dy_Cu(i,j),&
1459		G%IareaT(i,j),G%IareaT(i+1,j),G%IdxT(i,j),G%IdxT(i+1,j),dt,CS%vol_CFL,&
1460	50392650	por_face_areaU(i,j,k))
1461	50392650	if (local_open_BC) &
1462		call flux_elem_OBC(u_new,h_in(i,j,k),h_in(i+1,j,k),uh_3d(ii,jj,k),duhdu,&
1463		visc_rem(ii,jj,k),por_face_areaU(i,j,k),G%dy_Cu(i,j),OBC,&
1464	0	OBC%segnum_u(i,j))
1465	50392650	uh_err(ii,jj) = uh_err(ii,jj) + uh_3d(ii,jj,k)
1466	50392650	duhdu_tot(ii,jj) = duhdu_tot(ii,jj) + duhdu
1467		endif
1468		enddo ; enddo
1469		enddo
1470		!$omp loop collapse(2) private(ii,jj)
1471	1501570	do j=j_start,j_end ; do I=i_start,i_end
1472	1441960	ii=I-i_start+1 ; jj=j-j_start+1
1473	1489648	uh_err_best(ii,jj) = min(uh_err_best(ii,jj), abs(uh_err(ii,jj)))
1474		enddo ; enddo
1475		enddo ! itt-loop
1476		! If there are any faces which have not converged to within the tolerance,
1477		! so-be-it, or else use a final upwind correction?
1478		! This never seems to happen with 20 iterations as max_itt.
1479		!$omp end target
1480
1481		!$omp target exit data &
1482		!$omp map(release: uh_err,uh_err_best,duhdu_tot,du_min,du_max,do_I)
1483
1484	6528	end subroutine zonal_flux_adjust
1485
1486
1487		!> Sets a structure that describes the zonal barotropic volume or mass fluxes as a
1488		!! function of barotropic flow to agree closely with the sum of the layer's transports.
1489	3264	subroutine set_zonal_BT_cont(u, h_in, h_W, h_E, BT_cont, uh_tot_0, duhdu_tot_0, &
1490	3264	du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
1491	3264	visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
1492	3264	por_face_areaU, niblock, njblock)
1493		integer, intent(in) :: niblock !< i block size for array calculations [nondim].
1494		integer, intent(in) :: njblock !< j block size for array calculations [nondim].
1495		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1496		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1497		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
1498		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to
1499		!! calculate fluxes [H ~> m or kg m-2].
1500		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_W !< West edge thickness in the
1501		!! reconstruction [H ~> m or kg m-2].
1502		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_E !< East edge thickness in the
1503		!! reconstruction [H ~> m or kg m-2].
1504		type(BT_cont_type), intent(inout) :: BT_cont !< A structure with elements
1505		!! that describe the effective open face areas as a function of barotropic flow.
1506		real, dimension(niblock,njblock), intent(in) :: uh_tot_0 !< The summed transport
1507		!! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
1508		real, dimension(niblock,njblock), intent(in) :: duhdu_tot_0 !< The partial derivative
1509		!! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
1510		real, dimension(niblock,njblock), intent(in) :: du_max_CFL !< Maximum acceptable
1511		!! value of du [L T-1 ~> m s-1].
1512		real, dimension(niblock,njblock), intent(in) :: du_min_CFL !< Minimum acceptable
1513		!! value of du [L T-1 ~> m s-1].
1514		real, intent(in) :: dt !< Time increment [T ~> s].
1515		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1516		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
1517		real, dimension(niblock,njblock,SZK_(GV)),intent(in) :: visc_rem !< Both the fraction of the
1518		!! momentum originally in a layer that remains after a time-step of viscosity, and
1519		!! the fraction of a time-step's worth of a barotropic acceleration that a layer
1520		!! experiences after viscosity is applied [nondim].
1521		!! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).
1522		real, dimension(niblock,njblock), intent(in) :: visc_rem_max !< Maximum allowable visc_rem [nondim].
1523		integer, intent(in) :: i_start !< Start of i index range [nondim].
1524		integer, intent(in) :: i_end !< End of i index range [nondim].
1525		integer, intent(in) :: j_start !< Start of j index range [nondim].
1526		integer, intent(in) :: j_end !< End of j index range [nondim].
1527		logical, dimension(niblock,njblock),intent(in) :: do_I !< A logical flag indicating
1528		!! which I values to work on.
1529		real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &
1530		intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim]
1531		! Local variables
1532		real, dimension(niblock,njblock) :: &
1533	6528	du0 ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
1534		real, dimension(niblock,njblock) :: &
1535	6528	zeros ! An array of full of 0 transports [H L2 T-1 ~> m3 s-1 or kg s-1]
1536		real, dimension(niblock,njblock) :: &
1537	6528	duL, duR, & ! The barotropic velocity increments that give the westerly
1538		! (duL) and easterly (duR) test velocities [L T-1 ~> m s-1].
1539	6528	du_CFL, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
1540	6528	FAmt_L, FAmt_R, & ! The summed effective marginal face areas for the 3
1541	6528	FAmt_0, & ! test velocities [H L ~> m2 or kg m-1].
1542	6528	uhtot_L, & ! The summed transport with the westerly (uhtot_L) and
1543	6528	uhtot_R ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
1544		real :: &
1545		u_L, u_R, & ! The westerly (u_L), easterly (u_R), and zero-barotropic
1546		u_0, & ! transport (u_0) layer test velocities [L T-1 ~> m s-1].
1547		duhdu_L, & ! The effective layer marginal face areas with the westerly
1548		duhdu_R, & ! (_L), easterly (_R), and zero-barotropic (_0) test
1549		duhdu_0, & ! velocities [H L ~> m2 or kg m-1].
1550		uh_L, uh_R, & ! The layer transports with the westerly (_L), easterly (_R),
1551		uh_0 ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
1552		real :: FA_0 ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m-1].
1553		real :: FA_avg ! The average effective face area [L H ~> m2 or kg m-1], nominally given by
1554		! the realized transport divided by the barotropic velocity.
1555		real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]. This
1556		! limiting is necessary to keep the inverse of visc_rem
1557		! from leading to large CFL numbers.
1558		real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
1559		! in finding the barotropic velocity that changes the
1560		! flow direction [nondim]. This is necessary to keep the inverse
1561		! of visc_rem from leading to large CFL numbers.
1562		real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
1563		! flow is truly upwind [nondim]
1564		real :: Idt ! The inverse of the time step [T-1 ~> s-1].
1565		integer :: i, j, k, nz, ii, jj !< Tile loop indices [nondim].
1566	6528	real, dimension(niblock,njblock,SZK_(GV)) :: uh_tmp
1567
1568		!$omp target enter data &
1569		!$omp map(alloc: du0,duL,duR,du_CFL,FAmt_L,FAmt_R,FAmt_0,uhtot_L,uhtot_R,zeros,uh_tmp)
1570
1571	3264	nz = GV%ke ; Idt = 1.0 / dt
1572	3264	min_visc_rem = 0.1 ; CFL_min = 1e-6
1573
1574		! Diagnose the zero-transport correction, du0.
1575	102000	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
1576	485328	zeros(ii,jj) = 0.0
1577		enddo
1578		call zonal_flux_adjust(u, h_in, h_W, h_E, zeros, uh_tot_0, duhdu_tot_0, du0, &
1579		du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
1580	3264	i_start, i_end, j_start, j_end, do_I, por_face_areaU, niblock, njblock, uh_tmp)
1581
1582		! Determine the westerly- and easterly- fluxes. Choose a sufficiently
1583		! negative velocity correction for the easterly-flux, and a sufficiently
1584		! positive correction for the westerly-flux.
1585		!$omp target
1586		!$omp loop collapse(2) private(ii,jj)
1587	399264	do j=j_start,j_end ; do i=i_start,i_end
1588	383328	ii=i-i_start+1 ; jj=j-j_start+1
1589	383328	du_CFL(ii,jj) = (CFL_min * Idt) * G%dxCu(i,j)
1590	383328	duR(ii,jj) = min(0.0,du0(ii,jj) - du_CFL(ii,jj))
1591	383328	duL(ii,jj) = max(0.0,du0(ii,jj) + du_CFL(ii,jj))
1592	383328	FAmt_L(ii,jj) = 0.0 ; FAmt_R(ii,jj) = 0.0 ; FAmt_0(ii,jj) = 0.0
1593	396000	uhtot_L(ii,jj) = 0.0 ; uhtot_R(ii,jj) = 0.0
1594		enddo ; enddo
1595
1596	248064	do k=1,nz
1597		!$omp loop collapse(2) private(ii,jj,visc_rem_lim)
1598	29948064	do j=j_start,j_end ; do i=i_start,i_end
1599	28749600	ii=i-i_start+1 ; jj=j-j_start+1
1600	29700000	if (do_I(ii,jj)) then
1601	28749600	visc_rem_lim = max(visc_rem(ii,jj,k), min_visc_rem*visc_rem_max(ii,jj))
1602	28749600	if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
1603	17812800	if (u(i,j,k) + duR(ii,jj)visc_rem_lim > -du_CFL(ii,jj)visc_rem(ii,jj,k)) &
1604	232568	duR(ii,jj) = -(u(i,j,k) + du_CFL(ii,jj)*visc_rem(ii,jj,k)) / visc_rem_lim
1605	17812800	if (u(i,j,k) + duL(ii,jj)visc_rem_lim < du_CFL(ii,jj)visc_rem(ii,jj,k)) &
1606	3688842	duL(ii,jj) = -(u(i,j,k) - du_CFL(ii,jj)*visc_rem(ii,jj,k)) / visc_rem_lim
1607		endif
1608		endif
1609		enddo ; enddo
1610		enddo
1611		!$omp end target
1612
1613		!$omp target
1614	248064	do k=1,nz
1615		!$omp loop collapse(2) private(ii,jj,u_L,u_R,u_0,duhdu_0,duhdu_L,duhdu_R,uh_L,uh_R)
1616	29948064	do j=j_start,j_end ; do i=i_start,i_end
1617	28749600	ii=i-i_start+1 ; jj=j-j_start+1
1618	29700000	if (do_I(ii,jj)) then
1619	28749600	u_L = u(i,j,k) + duL(ii,jj) * visc_rem(ii,jj,k)
1620	28749600	u_R = u(i,j,k) + duR(ii,jj) * visc_rem(ii,jj,k)
1621	28749600	u_0 = u(i,j,k) + du0(ii,jj) * visc_rem(ii,jj,k)
1622		call flux_elem(u_0,h_in(i,j,k),h_in(i+1,j,k),h_W(i,j,k),h_W(i+1,j,k),h_E(i,j,k),&
1623		h_E(i+1,j,k),uh_0,duhdu_0,visc_rem(ii,jj,k),G%dy_Cu(i,j),&
1624		G%IareaT(i,j),G%IareaT(i+1,j),G%IdxT(i,j),G%IdxT(i+1,j),dt,CS%vol_CFL,&
1625	28749600	por_face_areaU(i,j,k))
1626		call flux_elem(u_L,h_in(i,j,k),h_in(i+1,j,k),h_W(i,j,k),h_W(i+1,j,k),h_E(i,j,k),&
1627		h_E(i+1,j,k),uh_L,duhdu_L,visc_rem(ii,jj,k),G%dy_Cu(i,j),&
1628		G%IareaT(i,j),G%IareaT(i+1,j),G%IdxT(i,j),G%IdxT(i+1,j),dt,CS%vol_CFL,&
1629	28749600	por_face_areaU(i,j,k))
1630		call flux_elem(u_R,h_in(i,j,k),h_in(i+1,j,k),h_W(i,j,k),h_W(i+1,j,k),h_E(i,j,k),&
1631		h_E(i+1,j,k),uh_R,duhdu_R,visc_rem(ii,jj,k),G%dy_Cu(i,j),&
1632		G%IareaT(i,j),G%IareaT(i+1,j),G%IdxT(i,j),G%IdxT(i+1,j),dt,CS%vol_CFL,&
1633	28749600	por_face_areaU(i,j,k))
1634	28749600	FAmt_0(ii,jj) = FAmt_0(ii,jj) + duhdu_0
1635	28749600	FAmt_L(ii,jj) = FAmt_L(ii,jj) + duhdu_L
1636	28749600	FAmt_R(ii,jj) = FAmt_R(ii,jj) + duhdu_R
1637	28749600	uhtot_L(ii,jj) = uhtot_L(ii,jj) + uh_L
1638	28749600	uhtot_R(ii,jj) = uhtot_R(ii,jj) + uh_R
1639		endif
1640		enddo ; enddo
1641		enddo
1642		!$omp loop collapse(2) private(ii,jj,FA_0,FA_avg)
1643	399264	do j=j_start,j_end ; do i=i_start,i_end
1644	383328	ii=i-i_start+1 ; jj=j-j_start+1
1645	396000	if (do_I(ii,jj)) then
1646	383328	FA_0 = FAmt_0(ii,jj) ; FA_avg = FAmt_0(ii,jj)
1647	383328	if ((duL(ii,jj) - du0(ii,jj)) /= 0.0) &
1648	383328	FA_avg = uhtot_L(ii,jj) / (duL(ii,jj) - du0(ii,jj))
1649	383328	if (FA_avg > max(FA_0, FAmt_L(ii,jj))) then ; FA_avg = max(FA_0, FAmt_L(ii,jj))
1650	382195	elseif (FA_avg < min(FA_0, FAmt_L(ii,jj))) then ; FA_0 = FA_avg ; endif
1651
1652	383328	BT_cont%FA_u_W0(i,j) = FA_0 ; BT_cont%FA_u_WW(i,j) = FAmt_L(ii,jj)
1653	383328	if (abs(FA_0-FAmt_L(ii,jj)) <= 1e-12*FA_0) then ; BT_cont%uBT_WW(i,j) = 0.0 ; else
1654		BT_cont%uBT_WW(i,j) = (1.5 * (duL(ii,jj) - du0(ii,jj))) * &
1655	237502	((FAmt_L(ii,jj) - FA_avg) / (FAmt_L(ii,jj) - FA_0))
1656		endif
1657
1658	383328	FA_0 = FAmt_0(ii,jj) ; FA_avg = FAmt_0(ii,jj)
1659	383328	if ((duR(ii,jj) - du0(ii,jj)) /= 0.0) &
1660	383328	FA_avg = uhtot_R(ii,jj) / (duR(ii,jj) - du0(ii,jj))
1661	383328	if (FA_avg > max(FA_0, FAmt_R(ii,jj))) then ; FA_avg = max(FA_0, FAmt_R(ii,jj))
1662	383323	elseif (FA_avg < min(FA_0, FAmt_R(ii,jj))) then ; FA_0 = FA_avg ; endif
1663
1664	383328	BT_cont%FA_u_E0(i,j) = FA_0 ; BT_cont%FA_u_EE(i,j) = FAmt_R(ii,jj)
1665	383328	if (abs(FAmt_R(ii,jj) - FA_0) <= 1e-12*FA_0) then ; BT_cont%uBT_EE(i,j) = 0.0 ; else
1666		BT_cont%uBT_EE(i,j) = (1.5 * (duR(ii,jj) - du0(ii,jj))) * &
1667	237504	((FAmt_R(ii,jj) - FA_avg) / (FAmt_R(ii,jj) - FA_0))
1668		endif
1669		else
1670	0	BT_cont%FA_u_W0(i,j) = 0.0 ; BT_cont%FA_u_WW(i,j) = 0.0
1671	0	BT_cont%FA_u_E0(i,j) = 0.0 ; BT_cont%FA_u_EE(i,j) = 0.0
1672	0	BT_cont%uBT_WW(i,j) = 0.0 ; BT_cont%uBT_EE(i,j) = 0.0
1673		endif
1674		enddo ; enddo
1675		!$omp end target
1676
1677		!$omp target exit data &
1678		!$omp map(release: du0,duL,duR,du_CFL,FAmt_L,FAmt_R,FAmt_0,uhtot_L,uhtot_R,zeros,uh_tmp)
1679
1680	3264	end subroutine set_zonal_BT_cont
1681
1682		!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.
1683	219	subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, por_face_areaV, &
1684	48	niblock, njblock, LB_in, vhbt, &
1685	120	visc_rem_v, v_cor, BT_cont, dv_cor)
1686		integer, intent(in) :: niblock !< i block size for array calculations [nondim].
1687		integer, intent(in) :: njblock !< j block size for array calculations [nondim].
1688		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1689		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1690		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
1691		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to
1692		!! calculate fluxes [H ~> m or kg m-2]
1693		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_S !< South edge thickness in the
1694		!! reconstruction [H ~> m or kg m-2].
1695		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_N !< North edge thickness in the
1696		!! reconstruction [H ~> m or kg m-2].
1697		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: vh !< Volume flux through meridional
1698		!! faces = vhdx [H L2 T-1 ~> m3 s-1 or kg s-1]
1699		real, intent(in) :: dt !< Time increment [T ~> s].
1700		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1701		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.G
1702		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition type
1703		!! specifies whether, where, and what
1704		!! open boundary conditions are used.
1705		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]
1706		type(cont_loop_bounds_type), optional, intent(in) :: LB_in !< Loop bounds structure.
1707		real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt !< The summed volume flux through meridional
1708		!! faces [H L2 T-1 ~> m3 s-1 or kg s-1].
1709		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1710		optional, intent(in) :: visc_rem_v !< Both the fraction of the momentum
1711		!! originally in a layer that remains after a time-step of viscosity,
1712		!! and the fraction of a time-step's worth of a barotropic acceleration
1713		!! that a layer experiences after viscosity is applied [nondim].
1714		!! Visc_rem_v is between 0 (at the bottom) and 1 (far above the bottom).
1715		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1716		optional, intent(out) :: v_cor
1717		!< The meridional velocities (v with a barotropic correction)
1718		!! that give vhbt as the depth-integrated transport [L T-1 ~> m s-1].
1719		type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
1720		!! the effective open face areas as a function of barotropic flow.
1721		real, dimension(SZI_(G),SZJB_(G)), &
1722		optional, intent(out) :: dv_cor !< The meridional velocity increments from v
1723		!! that give vhbt as the depth-integrated
1724		!! transports [L T-1 ~> m s-1].
1725
1726		! Local variables
1727		real, dimension(niblock,njblock) :: &
1728	146	vhbt_t
1729		real, dimension(niblock,njblock,SZK_(GV)) :: &
1730	146	vh_t
1731		real, dimension(niblock,njblock,SZK_(GV)) :: &
1732	146	dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
1733		real, dimension(niblock,njblock) :: &
1734	146	dv, & ! Corrective barotropic change in the velocity to give vhbt [L T-1 ~> m s-1].
1735	146	dv_min_CFL, & ! Lower limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]
1736	146	dv_max_CFL, & ! Upper limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]
1737	146	dvhdv_tot_0, & ! Summed partial derivative of vh with v [H L ~> m2 or kg m-1].
1738	146	vh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
1739	146	visc_rem_max ! The column maximum of visc_rem [nondim]
1740	146	logical, dimension(niblock,njblock) :: do_I
1741	146	real, dimension(niblock,njblock) :: FAvi ! A list of sums of meridional face areas
1742		! [H L ~> m2 or kg m-1].
1743		real :: FA_v ! A sum of meridional face areas [H L ~> m2 or kg m-1].
1744		real, dimension(niblock,njblock,SZK_(GV)) :: &
1745	146	visc_rem ! A 2-D copy of visc_rem_v or an array of 1's [nondim]
1746		real :: I_vrm ! 1.0 / visc_rem_max [nondim]
1747		real :: CFL_dt ! The maximum CFL ratio of the adjusted velocities divided by
1748		! the time step [T-1 ~> s-1].
1749		real :: I_dt ! 1.0 / dt [T-1 ~> s-1].
1750		real :: dv_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].
1751		real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [L ~> m].
1752		type(cont_loop_bounds_type) :: LB
1753		integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1754		integer :: l_seg ! The OBC segment number
1755		logical :: use_visc_rem, set_BT_cont
1756		logical :: local_specified_BC, local_Flather_OBC ! OBC-related logicals
1757		logical :: local_open_BC, any_simple_OBC ! OBC-related logicals
1758	146	logical :: simple_OBC_pt(niblock,njblock) ! Indicates points with specified transport OBCs
1759		type(OBC_segment_type), pointer :: segment => NULL()
1760		integer :: j_start, j_end, i_start, i_end, ii, jj
1761		integer :: nteams
1762
1763	73	call cpu_clock_begin(id_clock_correct)
1764
1765		!$omp target enter data &
1766		!$omp map(alloc:vhbt_t,vh_t,dvhdv,dv,dv_min_CFL,dv_max_CFL,dvhdv_tot_0,vh_tot_0, &
1767		!$omp visc_rem_max,do_I,visc_rem,simple_OBC_pt)
1768
1769	73	use_visc_rem = present(visc_rem_v)
1770
1771	73	set_BT_cont = .false. ; if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
1772
1773	73	local_specified_BC = .false. ; local_Flather_OBC = .false. ; local_open_BC = .false.
1774	73	if (associated(OBC)) then ; if (OBC%OBC_pe) then
1775	0	local_specified_BC = OBC%specified_v_BCs_exist_globally
1776	0	local_Flather_OBC = OBC%Flather_v_BCs_exist_globally
1777	0	local_open_BC = OBC%open_v_BCs_exist_globally
1778		endif ; endif
1779
1780	73	if (present(dv_cor)) then
1781	0	do concurrent (jj=G%jsdB:G%jedB, ii=G%isd:G%ied)
1782	0	dv_cor(ii,jj) = 0.0
1783		enddo
1784		endif
1785
1786	73	if (present(LB_in)) then
1787	73	LB = LB_in
1788		else
1789	0	LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
1790		endif
1791	73	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
1792
1793	73	CFL_dt = CS%CFL_limit_adjust / dt
1794	73	I_dt = 1.0 / dt
1795	73	if (CS%aggress_adjust) CFL_dt = I_dt
1796
1797	73	if (.not.use_visc_rem) then
1798	1	do concurrent (k=1:nz, jj=1:njblock, ii=1:niblock)
1799	9761	visc_rem(ii,jj,k) = 1.0
1800		enddo
1801		endif
1802
1803	7081	do j_start = jsh-1, jeh, njblock ; do i_start = ish, ieh, niblock
1804	4672	j_end = min(j_start + njblock-1, jeh)
1805	4672	i_end = min(i_start + niblock-1, ieh)
1806		!$ nteams = ceiling(real((j_end-j_start+1)*(i_end-i_start+1))/128.)
1807	144832	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
1808	679192	do_I(ii,JJ) = .true.
1809		enddo
1810		! This sets vh and dvhdv.
1811		!$omp target teams num_teams(nteams)
1812	355072	do k=1,nz
1813	350400	if (use_visc_rem) then
1814		!$omp loop collapse(2) private(ii,jj)
1815	41191200	do J=j_start,j_end ; do i=i_start,i_end
1816	39528000	ii=i-i_start+1 ; jj=j-j_start+1
1817	40845600	visc_rem(ii,JJ,k) = visc_rem_v(i,J,k)
1818		enddo ; enddo
1819		endif
1820		!$omp loop collapse(2) private(ii,JJ)
1821	41763300	do J=j_start,j_end ; do i=i_start,i_end
1822	40077000	ii=i-i_start+1 ; jj=j-j_start+1
1823		call flux_elem(v(i,J,k),h_in(i,J,k),h_in(i,J+1,k),h_S(i,J,k),h_S(i,J+1,k),&
1824		h_N(i,J,k),h_N(i,J+1,k),vh_t(ii,JJ,k),dvhdv(ii,JJ,k),visc_rem(ii,JJ,k),&
1825		G%dx_Cv(i,J),G%IareaT(i,J),G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,&
1826	40077000	por_face_areaV(i,J,k))
1827	40077000	if (local_open_BC) &
1828		call flux_elem_OBC(v(i,J,k),h_in(i,J,k),h_in(i,J+1,k),vh_t(ii,JJ,k),dvhdv(ii,JJ,k),&
1829		visc_rem(ii,JJ,k),por_face_areaV(i,J,k),G%dx_Cv(i,J),OBC,&
1830	1335900	OBC%segnum_v(i,J))
1831		enddo ; enddo
1832	355072	if (local_specified_BC) then
1833		!$omp loop collapse(2) private(l_seg)
1834	0	do J=J_start,J_end ; do i=i_start,i_end ; if (OBC%segnum_v(i,J) /= 0) then
1835	0	l_seg = abs(OBC%segnum_v(i,J))
1836	0	if (OBC%segment(l_seg)%specified) vh_t(i-i_start+1,J-J_start+1,k) = OBC%segment(l_seg)%normal_trans(i,J,k)
1837		endif ; enddo ; enddo
1838		endif
1839		enddo ! k-loop
1840		!$omp end target teams
1841
1842	4672	if (present(vhbt) .or. set_BT_cont) then
1843	4608	if (use_visc_rem .and. CS%use_visc_rem_max) then
1844		!$omp target
1845		!$omp loop collapse(2) private(ii,jj)
1846	549216	do j=j_start,j_end ; do i=i_start,i_end
1847	527040	ii=i-i_start+1 ; jj=j-j_start+1
1848	544608	visc_rem_max(ii,JJ) = 0.0
1849		enddo ; enddo
1850	350208	do k=1,nz
1851		!$omp loop collapse(2) private(ii,jj)
1852	41195808	do j=j_start,j_end ; do i=i_start,i_end
1853	39528000	ii=i-i_start+1 ; jj=j-j_start+1
1854	40845600	visc_rem_max(ii,JJ) = max(visc_rem_max(ii,JJ), visc_rem(ii,JJ,k))
1855		enddo ; enddo
1856		enddo
1857		!$omp end target
1858		else
1859	0	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
1860	0	visc_rem_max(ii,JJ) = 1.0
1861		enddo
1862		endif
1863		! Set limits on dv that will keep the CFL number between -1 and 1.
1864		! This should be adequate to keep the root bracketed in all cases.
1865		!$omp target
1866		!$omp loop collapse(2) private(ii,jj,I_vrm,dy_S,dy_N)
1867	549216	do J=J_start,J_end ; do i=i_start,i_end
1868	527040	ii=i-i_start+1 ; jj=j-j_start+1
1869	527040	I_vrm = 0.0
1870	527040	if (visc_rem_max(ii,JJ) > 0.0) I_vrm = 1.0 / visc_rem_max(ii,JJ)
1871	527040	if (CS%vol_CFL) then
1872	0	dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(i,J), 1000.0*G%dyT(i,j))
1873	0	dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(i,J), 1000.0*G%dyT(i,j+1))
1874	527040	else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
1875	527040	dv_max_CFL(ii,JJ) = 2.0 * (CFL_dt * dy_S) * I_vrm
1876	527040	dv_min_CFL(ii,JJ) = -2.0 * (CFL_dt * dy_N) * I_vrm
1877	544608	vh_tot_0(ii,JJ) = 0.0 ; dvhdv_tot_0(ii,JJ) = 0.0
1878		enddo ; enddo
1879	350208	do k=1,nz
1880		!$omp loop collapse(2) private(ii,jj)
1881	41195808	do J=J_start,J_end ; do i=i_start,i_end
1882	39528000	ii=i-i_start+1 ; jj=j-j_start+1
1883	39528000	dvhdv_tot_0(ii,JJ) = dvhdv_tot_0(ii,JJ) + dvhdv(ii,jj,k)
1884	40845600	vh_tot_0(ii,JJ) = vh_tot_0(ii,JJ) + vh_t(ii,JJ,k)
1885		enddo
1886		enddo ; enddo
1887		!$omp end target
1888
1889	4608	if (use_visc_rem) then
1890	4608	if (CS%aggress_adjust) then
1891		!$omp target
1892	0	do k=1,nz
1893		!$omp loop collapse(2) private(ii,jj,dy_S,dy_N,dv_lim)
1894	0	do J=J_start,J_end ; do i=i_start,i_end
1895	0	ii=i-i_start+1 ; jj=j-j_start+1
1896	0	if (CS%vol_CFL) then
1897	0	dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
1898	0	dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
1899	0	else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
1900	0	dv_lim = 0.499((dy_SI_dt - v(i,J,k)) + MIN(0.0,v(i,J-1,k)))
1901	0	if (dv_max_CFL(ii,JJ) * visc_rem(ii,JJ,k) > dv_lim) &
1902	0	dv_max_CFL(ii,JJ) = dv_lim / visc_rem(ii,JJ,k)
1903
1904	0	dv_lim = 0.499((-dy_NCFL_dt - v(i,J,k)) + MAX(0.0,v(i,J+1,k)))
1905	0	if (dv_min_CFL(ii,JJ) * visc_rem(ii,JJ,k) < dv_lim) &
1906	0	dv_min_CFL(ii,JJ) = dv_lim / visc_rem(ii,JJ,k)
1907		enddo
1908		enddo ; enddo
1909		!$omp end target
1910		else
1911		!$omp target
1912	350208	do k=1,nz
1913		!$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
1914	41195808	do J=J_start,J_end ; do i=i_start,i_end
1915	39528000	ii=i-i_start+1 ; jj=j-j_start+1
1916	39528000	if (CS%vol_CFL) then
1917	0	dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
1918	0	dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
1919	39528000	else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
1920	39528000	if (dv_max_CFL(ii,JJ) * visc_rem(ii,JJ,k) > &
1921		dy_SCFL_dt - v(i,J,k)G%mask2dCv(i,J)) &
1922	569593	dv_max_CFL(ii,JJ) = (dy_S*CFL_dt - v(i,J,k)) / visc_rem(ii,JJ,k)
1923	39528000	if (dv_min_CFL(ii,JJ) * visc_rem(ii,JJ,k) < &
1924		-dy_NCFL_dt - v(i,J,k)G%mask2dCv(i,J)) &
1925	9948917	dv_min_CFL(ii,JJ) = -(dy_N*CFL_dt + v(i,J,k)) / visc_rem(ii,JJ,k)
1926		enddo ; enddo
1927		enddo
1928		!$omp end target
1929		endif
1930		else
1931	0	if (CS%aggress_adjust) then
1932		!$omp target
1933	0	do k=1,nz
1934		!$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
1935	0	do J=J_start,J_end ; do i=i_start,i_end
1936	0	ii=i-i_start+1 ; jj=j-j_start+1
1937	0	if (CS%vol_CFL) then
1938	0	dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
1939	0	dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
1940	0	else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
1941		dv_max_CFL(ii,JJ) = min(dv_max_CFL(ii,JJ), 0.499 * &
1942	0	((dy_S*I_dt - v(i,J,k)) + MIN(0.0,v(i,J-1,k))) )
1943		dv_min_CFL(ii,JJ) = max(dv_min_CFL(ii,JJ), 0.499 * &
1944	0	((-dy_N*I_dt - v(i,J,k)) + MAX(0.0,v(i,J+1,k))) )
1945		enddo ; enddo
1946		enddo
1947		!$omp end target
1948		else
1949		!$omp target
1950	0	do k=1,nz
1951		!$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
1952	0	do J=J_start,J_end ; do i=i_start,i_end
1953	0	ii=i-i_start+1 ; jj=j-j_start+1
1954	0	if (CS%vol_CFL) then
1955	0	dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
1956	0	dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
1957	0	else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
1958	0	dv_max_CFL(ii,JJ) = min(dv_max_CFL(ii,JJ), dy_S*CFL_dt - v(i,J,k))
1959	0	dv_min_CFL(ii,JJ) = max(dv_min_CFL(ii,JJ), -(dy_N*CFL_dt + v(i,J,k)))
1960		enddo ; enddo
1961		enddo
1962		!$omp end target
1963		endif
1964		endif
1965	142848	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
1966	527040	dv_max_CFL(ii,JJ) = max(dv_max_CFL(ii,JJ),0.0)
1967	669888	dv_min_CFL(ii,JJ) = min(dv_min_CFL(ii,JJ),0.0)
1968		enddo
1969
1970	4608	any_simple_OBC = .false.
1971	4608	if (present(vhbt) .or. set_BT_cont) then
1972	4608	if (local_specified_BC .or. local_Flather_OBC) then
1973	0	do concurrent (J=J_start:J_end, i=i_start:i_end) DO_LOCALITY(reduce(.or.:any_simple_OBC))
1974	0	ii=i-i_start+1 ; jj=J-j_start+1
1975	0	l_seg = abs(OBC%segnum_v(i,J))
1976
1977		! Avoid reconciling barotropic/baroclinic transports if transport is specified
1978	0	simple_OBC_pt(ii,jj) = .false.
1979	0	if (l_seg /= 0) simple_OBC_pt(ii,jj) = OBC%segment(l_seg)%specified
1980	0	do_I(ii,jj) = .not.simple_OBC_pt(ii,jj)
1981	0	any_simple_OBC = any_simple_OBC .or. simple_OBC_pt(ii,jj)
1982		enddo ; else
1983	142848	do concurrent (JJ=1:J_end-J_start+1, ii=1:i_end-i_start+1)
1984	669888	do_I(ii,JJ) = .true.
1985		enddo ; endif
1986		endif
1987
1988	4608	if (present(vhbt)) then
1989	3072	do concurrent (j=j_start:j_end, i=i_start:i_end)
1990	351360	ii=i-i_start+1 ; jj=j-j_start+1
1991	446592	vhbt_t(ii,Jj) = vhbt(i,j)
1992		enddo
1993		! Find dv and vh.
1994		call meridional_flux_adjust(v, h_in, h_S, h_N, vhbt_t, vh_tot_0, dvhdv_tot_0, dv, &
1995		dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
1996		i_start, i_end, j_start, j_end, do_I, por_face_areaV, &
1997	3072	niblock, njblock, vh_t, OBC)
1998
1999	3072	if (present(v_cor)) then
2000	3072	do concurrent (k=1:nz, J=J_start:J_end, i=i_start:i_end)
2001	26352000	ii=i-i_start+1 ; jj=j-j_start+1
2002	26352000	v_cor(i,J,k) = v(i,J,k) + dv(ii,JJ) * visc_rem(ii,JJ,k)
2003	26798592	if (any_simple_OBC) then ; if (simple_OBC_pt(ii,jj)) then
2004	0	v_cor(i,J,k) = OBC%segment(abs(OBC%segnum_v(i,J)))%normal_vel(i,J,k)
2005		endif ; endif
2006		enddo
2007		endif ! v-corrected
2008
2009	3072	if (present(dv_cor)) then
2010	0	do concurrent (J=J_start:J_end, i=i_start:i_end)
2011	0	ii=i-i_start+1 ; jj=j-j_start+1
2012	0	dv_cor(i,J) = dv(ii,JJ)
2013		enddo
2014		endif
2015		endif
2016
2017	4608	if (set_BT_cont) then
2018		call set_merid_BT_cont(v, h_in, h_S, h_N, BT_cont, vh_tot_0, dvhdv_tot_0, &
2019		dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
2020		visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
2021	3072	por_face_areaV, niblock, njblock)
2022	3072	if (any_simple_OBC) then
2023		!$omp target
2024		!$omp loop collapse(2) private(ii,jj)
2025	0	do J=J_start,J_end ; do i=i_start,i_end
2026	0	ii=i-i_start+1 ; jj=j-j_start+1
2027	0	if (simple_OBC_pt(ii,JJ)) FAvi(ii,JJ) = GV%H_subroundoff*G%dx_Cv(i,J)
2028		enddo ; enddo
2029		! NOTE: simple_OBC_pt should prevent access to segment OBC_NONE
2030	0	do k=1,nz
2031		!$omp loop collapse(2) private(ii,jj,l_seg)
2032	0	do J=J_start,J_end ; do i=i_start,i_end
2033	0	ii=i-i_start+1 ; jj=j-j_start+1
2034	0	if (simple_OBC_pt(ii,JJ)) then
2035	0	l_seg = abs(OBC%segnum_v(i,J))
2036	0	if ((abs(OBC%segment(l_seg)%normal_vel(i,J,k)) > 0.0) .and. &
2037		(OBC%segment(l_seg)%specified)) &
2038		FAvi(ii,JJ) = FAvi(ii,JJ) + &
2039	0	OBC%segment(l_seg)%normal_trans(i,J,k) / OBC%segment(l_seg)%normal_vel(i,J,k)
2040		endif
2041		enddo ; enddo
2042		enddo
2043		!$omp loop collapse(2) private(ii,jj)
2044	0	do J=J_start,J_end ; do i=i_start,i_end
2045	0	ii=i-i_start+1 ; jj=j-j_start+1
2046	0	if (simple_OBC_pt(ii,jj)) then
2047	0	BT_cont%FA_v_S0(i,J) = FAvi(ii,jj) ; BT_cont%FA_v_N0(i,J) = FAvi(ii,jj)
2048	0	BT_cont%FA_v_SS(i,J) = FAvi(ii,jj) ; BT_cont%FA_v_NN(i,J) = FAvi(ii,jj)
2049	0	BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0
2050		endif
2051		enddo ; enddo
2052		!$omp end target
2053		endif
2054		endif ! set_BT_cont
2055
2056		endif ! present(vhbt) or set_BT_cont
2057
2058	1168	do concurrent (k=1:nz, j=j_start:j_end, i=i_start:i_end)
2059	40077000	jj=j-j_start+1 ; ii=i-i_start+1
2060	40756192	vh(i,j,k) = vh_t(ii,jj,k)
2061		enddo
2062
2063		enddo ; enddo ! ij block loops
2064
2065		!$omp target exit data &
2066		!$omp map(release:vhbt_t,vh_t,dvhdv,dv,dv_min_CFL,dv_max_CFL,dvhdv_tot_0,vh_tot_0,&
2067		!$omp visc_rem_max,do_I,visc_rem,simple_OBC_pt)
2068
2069	73	if (local_open_BC .and. set_BT_cont) then
2070	0	do n = 1, OBC%number_of_segments
2071	0	if (OBC%segment(n)%open .and. OBC%segment(n)%is_N_or_S) then
2072	0	J = OBC%segment(n)%HI%JsdB
2073	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
2074	0	do i = OBC%segment(n)%HI%Isd, OBC%segment(n)%HI%Ied
2075	0	FA_v = 0.0
2076	0	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
2077	0	BT_cont%FA_v_S0(i,J) = FA_v ; BT_cont%FA_v_N0(i,J) = FA_v
2078	0	BT_cont%FA_v_SS(i,J) = FA_v ; BT_cont%FA_v_NN(i,J) = FA_v
2079	0	BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0
2080		enddo
2081		else
2082	0	do i = OBC%segment(n)%HI%Isd, OBC%segment(n)%HI%Ied
2083	0	FA_v = 0.0
2084	0	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
2085	0	BT_cont%FA_v_S0(i,J) = FA_v ; BT_cont%FA_v_N0(i,J) = FA_v
2086	0	BT_cont%FA_v_SS(i,J) = FA_v ; BT_cont%FA_v_NN(i,J) = FA_v
2087	0	BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0
2088		enddo
2089		endif
2090		endif
2091		enddo
2092		endif
2093
2094	73	if (set_BT_cont) then ; if (allocated(BT_cont%h_v)) then
2095	48	if (present(v_cor)) then
2096		call meridional_flux_thickness(v_cor, h_in, h_S, h_N, BT_cont%h_v, dt, G, GV, US, LB, &
2097	24	CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaV, visc_rem_v)
2098		else
2099		call meridional_flux_thickness(v, h_in, h_S, h_N, BT_cont%h_v, dt, G, GV, US, LB, &
2100	24	CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaV, visc_rem_v)
2101		endif
2102		endif ; endif
2103
2104	73	call cpu_clock_end(id_clock_correct)
2105
2106	193	end subroutine meridional_mass_flux
2107
2108
2109		!> Calculates the vertically integrated mass or volume fluxes through the meridional faces.
2110	0	subroutine meridional_BT_mass_flux(v, h_in, h_S, h_N, vhbt, dt, G, GV, US, CS, OBC, por_face_areaV, LB_in)
2111		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2112		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2113		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
2114		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to
2115		!! calculate fluxes [H ~> m or kg m-2]
2116		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_S !< Southern edge thickness in the PPM
2117		!! reconstruction [H ~> m or kg m-2].
2118		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_N !< Northern edge thickness in the PPM
2119		!! reconstruction [H ~> m or kg m-2].
2120		real, dimension(SZI_(G),SZJB_(G)), intent(out) :: vhbt !< The summed volume flux through meridional
2121		!! faces [H L2 T-1 ~> m3 s-1 or kg s-1].
2122		real, intent(in) :: dt !< Time increment [T ~> s].
2123		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2124		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
2125		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition type
2126		!! specifies whether, where, and what
2127		!! open boundary conditions are used.
2128		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]
2129		type(cont_loop_bounds_type), optional, intent(in) :: LB_in !< Loop bounds structure.
2130
2131		! Local variables
2132	0	real :: vh(SZI_(G),SZJB_(G),SZK_(GV)) ! Volume flux through meridional faces = vhdx [H L2 T-1 ~> m3 s-1 or kg s-1]
2133	0	real :: dvhdv(SZI_(G),SZJB_(G),SZK_(GV)) ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
2134		integer :: i, j, k, ish, ieh, jsh, jeh, nz, l_seg
2135	0	logical :: local_specified_BC, OBC_in_row(SZJB_(G))
2136
2137	0	call cpu_clock_begin(id_clock_correct)
2138
2139	0	local_specified_BC = .false.
2140	0	if (associated(OBC)) then ; if (OBC%OBC_pe) then
2141	0	local_specified_BC = OBC%specified_v_BCs_exist_globally
2142		endif ; endif
2143
2144	0	if (present(LB_in)) then
2145	0	ish = LB_in%ish ; ieh = LB_in%ieh ; jsh = LB_in%jsh ; jeh = LB_in%jeh ; nz = GV%ke
2146		else
2147	0	ish = G%isc ; ieh = G%iec ; jsh = G%jsc ; jeh = G%jec ; nz = GV%ke
2148		endif
2149
2150	0	vhbt(:,:) = 0.0
2151
2152		! Determining whether there are any OBC points outside of the k-loop should be more efficient.
2153	0	OBC_in_row(:) = .false.
2154	0	if (local_specified_BC) then
2155	0	do j=jsh-1,jeh ; do i=ish,ieh ; if (OBC%segnum_v(i,J) /= 0) then
2156	0	if (OBC%segment(abs(OBC%segnum_v(i,J)))%specified) OBC_in_row(j) = .true.
2157		endif ; enddo ; enddo
2158		endif
2159
2160		! This sets vh and dvhdv.
2161	0	do concurrent (k=1:nz, J=jsh-1:jeh, i=ish:ieh)
2162		call flux_elem(v(i,J,k), h_in(i,J,k), h_in(i,J+1,k), h_S(i,J,k), h_S(i,J+1,k), &
2163		h_N(i,J,k), h_N(i,J+1,k), vh(i,J,k), dvhdv(i,J,k), 1.0, G%dx_Cv(I,j), &
2164		G%IareaT(i,J), G%IareaT(i,J+1), G%IdyT(i,J), G%IdyT(i,J+1), dt, &
2165	0	CS%vol_CFL, por_face_areaV(i,J,k))
2166	0	if (local_specified_BC) &
2167		call flux_elem_OBC(v(i,J,k), h_in(i,J,k), h_in(i,J+1,k), vh(i,J,k), dvhdv(i,J,k), 1.0, &
2168	0	por_face_areaV(i,J,k), G%dx_Cv(i,J), OBC, OBC%segnum_v(i,J))
2169		enddo
2170
2171	0	do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
2172	0	if (OBC_in_row(j) .and. OBC%segnum_v(i,J) /= 0) then
2173	0	l_seg = abs(OBC%segnum_v(i,J))
2174	0	if (OBC%segment(l_seg)%specified) vh(i,j,k) = OBC%segment(l_seg)%normal_trans(i,J,k)
2175		endif
2176		enddo ; enddo ; enddo
2177
2178
2179		! Accumulate the barotropic transport.
2180	0	do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh
2181	0	vhbt(i,J) = vhbt(i,J) + vh(i,J,k)
2182		enddo ; enddo ; enddo
2183
2184	0	call cpu_clock_end(id_clock_correct)
2185
2186	0	end subroutine meridional_BT_mass_flux
2187
2188
2189		!> Sets the effective interface thickness associated with the fluxes at each meridional velocity point,
2190		!! optionally scaling back these thicknesses to account for viscosity and fractional open areas.
2191	48	subroutine meridional_flux_thickness(v, h, h_S, h_N, h_v, dt, G, GV, US, LB, vol_CFL, &
2192	48	marginal, OBC, por_face_areaV, visc_rem_v)
2193		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2194		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2195		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
2196		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to calculate fluxes,
2197		!! [H ~> m or kg m-2].
2198		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_S !< South edge thickness in the reconstruction,
2199		!! [H ~> m or kg m-2].
2200		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_N !< North edge thickness in the reconstruction,
2201		!! [H ~> m or kg m-2].
2202		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: h_v !< Effective thickness at meridional faces,
2203		!! scaled down to account for the effects of
2204		!! viscosity and the fractional open area
2205		!! [H ~> m or kg m-2].
2206		real, intent(in) :: dt !< Time increment [T ~> s].
2207		type(cont_loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
2208		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2209		logical, intent(in) :: vol_CFL !< If true, rescale the ratio
2210		!! of face areas to the cell areas when estimating the CFL number.
2211		logical, intent(in) :: marginal !< If true, report the marginal
2212		!! face thicknesses; otherwise report transport-averaged thicknesses.
2213		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
2214		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
2215		intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]
2216		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), optional, intent(in) :: visc_rem_v !< Both the fraction
2217		!! of the momentum originally in a layer that remains after a time-step of
2218		!! viscosity, and the fraction of a time-step's worth of a barotropic
2219		!! acceleration that a layer experiences after viscosity is applied [nondim].
2220		!! Visc_rem_v is between 0 (at the bottom) and 1 (far above the bottom).
2221
2222		! Local variables
2223		real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
2224		real :: curv_3 ! A measure of the thickness curvature over a grid length,
2225		! with the same units as h [H ~> m or kg m-2] .
2226		real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].
2227		real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
2228		logical :: local_open_BC
2229		integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
2230		real :: dh
2231	48	ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
2232
2233	48	do concurrent (k=1:nz, J=jsh-1:jeh, i=ish:ieh)
2234	26352000	if (v(i,J,k) > 0.0) then
2235	5937285	if (vol_CFL) then ; CFL = (v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j))
2236	5937285	else ; CFL = v(i,J,k) * dt * G%IdyT(i,j) ; endif
2237	5937285	curv_3 = (h_S(i,j,k) + h_N(i,j,k)) - 2.0*h(i,j,k)
2238	5937285	dh = h_S(i,J,k) - h_N(i,J,k)
2239	5937285	if (marginal) then
2240	5937285	h_v(i,J,k) = h_N(i,j,k) + CFL * (dh + 3.0curv_3(CFL - 1.0))
2241		else
2242	0	h_v(i,J,k) = h_N(i,j,k) + CFL * (0.5dh + curv_3(CFL - 1.5))
2243		endif
2244	20414715	elseif (v(i,J,k) < 0.0) then
2245	11314065	if (vol_CFL) then ; CFL = (-v(i,J,k)dt) (G%dx_Cv(i,J) * G%IareaT(i,j+1))
2246	11314065	else ; CFL = -v(i,J,k) * dt * G%IdyT(i,j+1) ; endif
2247	11314065	curv_3 = (h_S(i,j+1,k) + h_N(i,j+1,k)) - 2.0*h(i,j+1,k)
2248	11314065	dh = h_N(i,j+1,k)-h_S(i,j+1,k)
2249	11314065	if (marginal) then
2250	11314065	h_v(i,J,k) = h_S(i,j+1,k) + CFL * (dh + 3.0curv_3(CFL - 1.0))
2251		else
2252	0	h_v(i,J,k) = h_S(i,j+1,k) + CFL * (0.5dh + curv_3(CFL - 1.5))
2253		endif
2254		else
2255		! The choice to use the arithmetic mean here is somewhat arbitrarily, but
2256		! it should be noted that h_S(i+1,j,k) and h_N(i,j,k) are usually the same.
2257	9100650	h_v(i,J,k) = 0.5 * (h_S(i,j+1,k) + h_N(i,j,k))
2258		! h_marg = (2.0 * h_S(i,j+1,k) * h_N(i,j,k)) / &
2259		! (h_S(i,j+1,k) + h_N(i,j,k) + GV%H_subroundoff)
2260		endif
2261
2262	53061168	if (present(visc_rem_v)) then
2263		! Scale back the thickness to account for the effects of viscosity and the fractional open
2264		! thickness to give an appropriate non-normalized weight for each layer in determining the
2265		! barotropic acceleration.
2266	26352000	h_v(i,J,k) = h_v(i,J,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
2267		else
2268	0	h_v(i,J,k) = h_v(i,J,k) * por_face_areaV(i,J,k)
2269		endif
2270		enddo
2271
2272	48	local_open_BC = .false.
2273	48	if (associated(OBC)) local_open_BC = OBC%open_v_BCs_exist_globally
2274		! untested - will need to be refactored to be performant on GPUs
2275	48	if (local_open_BC) then
2276	0	do n = 1, OBC%number_of_segments
2277	0	if (OBC%segment(n)%open .and. OBC%segment(n)%is_N_or_S) then
2278	0	J = OBC%segment(n)%HI%JsdB
2279	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
2280	0	if (present(visc_rem_v)) then
2281	0	do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
2282	0	h_v(i,J,k) = h(i,J,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
2283		enddo
2284		else
2285	0	do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
2286	0	h_v(i,J,k) = h(i,J,k) * por_face_areaV(i,J,k)
2287		enddo
2288		endif
2289		else
2290	0	if (present(visc_rem_v)) then
2291	0	do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
2292	0	h_v(i,J,k) = h(i,J+1,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
2293		enddo
2294		else
2295	0	do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
2296	0	h_v(i,J,k) = h(i,J+1,k) * por_face_areaV(i,J,k)
2297		enddo
2298		endif
2299		endif
2300		endif
2301		enddo
2302		endif
2303
2304	48	end subroutine meridional_flux_thickness
2305
2306
2307		!> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.
2308	6144	subroutine meridional_flux_adjust(v, h_in, h_S, h_N, vhbt, vh_tot_0, dvhdv_tot_0, &
2309	6144	dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
2310	12288	i_start, i_end, j_start, j_end, do_I_in, por_face_areaV, &
2311	6144	niblock, njblock, vh_3d, OBC)
2312		integer, intent(in) :: niblock !< i block size for array calculations [nondim].
2313		integer, intent(in) :: njblock !< j block size for array calculations [nondim].
2314		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2315		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2316		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2317		intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
2318		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2319		intent(in) :: h_in !< Layer thickness used to calculate fluxes
2320		!! [H ~> m or kg m-2].
2321		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2322		intent(in) :: h_S !< South edge thickness in the reconstruction
2323		!! [H ~> m or kg m-2].
2324		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2325		intent(in) :: h_N !< North edge thickness in the reconstruction
2326		!! [H ~> m or kg m-2].
2327		real, dimension(niblock,njblock,SZK_(GV)), intent(in) :: visc_rem
2328		!< Both the fraction of the momentum originally
2329		!! in a layer that remains after a time-step of viscosity, and the
2330		!! fraction of a time-step's worth of a barotropic acceleration that
2331		!! a layer experiences after viscosity is applied [nondim].
2332		!! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).
2333		real, dimension(niblock,njblock), intent(in) :: vhbt !< The summed volume flux through
2334		!! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
2335		real, dimension(niblock,njblock), intent(in) :: dv_max_CFL !< Maximum acceptable value
2336		!! of dv [L T-1 ~> m s-1].
2337		real, dimension(niblock,njblock), intent(in) :: dv_min_CFL !< Minimum acceptable value
2338		!! of dv [L T-1 ~> m s-1].
2339		real, dimension(niblock,njblock), intent(in) :: vh_tot_0 !< The summed transport with
2340		!! 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
2341		real, dimension(niblock,njblock), intent(in) :: dvhdv_tot_0 !< The partial derivative
2342		!! of dv_err with dv at 0 adjustment [H L ~> m2 or kg m-1].
2343		real, dimension(niblock,njblock), intent(out) :: dv !< The barotropic velocity adjustment
2344		!! [L T-1 ~> m s-1].
2345		real, intent(in) :: dt !< Time increment [T ~> s].
2346		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2347		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
2348		integer, intent(in) :: j_start !< Spatial index.
2349		integer, intent(in) :: j_end !< Spatial index.
2350		integer, intent(in) :: i_start !< Start of index range.
2351		integer, intent(in) :: i_end !< End of index range.
2352		logical, dimension(niblock,njblock), &
2353		intent(in) :: do_I_in !< A flag indicating which I values to work on.
2354		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
2355		intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]
2356		real, dimension(niblock,njblock,SZK_(GV)), &
2357		optional, intent(inout) :: vh_3d !< Volume flux through meridional
2358		!! faces = vhdx [H L2 T-1 ~> m3 s-1 or kg s-1].
2359		type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
2360		! Local variables
2361		real :: &
2362		dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
2363		real, dimension(niblock,njblock) :: &
2364	12288	vh_err, & ! Difference between vhbt and the summed vh [H L2 T-1 ~> m3 s-1 or kg s-1].
2365	12288	vh_err_best, & ! The smallest value of vh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
2366	12288	dvhdv_tot,&! Summed partial derivative of vh with u [H L ~> m2 or kg m-1].
2367	12288	dv_min, & ! Lower limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
2368	12288	dv_max ! Upper limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
2369		real :: &
2370		v_new ! The velocity with the correction added [L T-1 ~> m s-1].
2371		real :: dv_prev ! The previous value of dv [L T-1 ~> m s-1].
2372		real :: ddv ! The change in dv from the previous iteration [L T-1 ~> m s-1].
2373		real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
2374		real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
2375		integer :: i, j, k, nz, itt, ii, jj
2376	12288	logical :: domore, do_I(niblock,njblock), local_open_BC
2377		integer, parameter :: max_itts = 20
2378
2379		!$omp target enter data &
2380		!$omp map(alloc:vh_err,vh_err_best,dvhdv_tot,dv_min,dv_max,do_I)
2381
2382	6144	local_open_BC = .false.
2383	6144	if (present(OBC)) then ; if (associated(OBC)) then
2384	0	local_open_BC = OBC%open_v_BCs_exist_globally
2385		endif ; endif
2386
2387	6144	nz = GV%ke
2388	6144	tol_vel = CS%tol_vel
2389
2390		!$omp target
2391
2392		!$omp loop collapse(2) private(ii,jj)
2393	732288	do J=J_start,J_end ; do i=i_start,i_end
2394	702720	ii=i-i_start+1 ; jj=j-j_start+1
2395	702720	dv(ii,JJ) = 0.0 ; do_I(ii,JJ) = do_I_in(ii,JJ)
2396	702720	dv_max(ii,JJ) = dv_max_CFL(ii,JJ) ; dv_min(ii,JJ) = dv_min_CFL(ii,JJ)
2397	702720	vh_err(ii,JJ) = vh_tot_0(ii,JJ) - vhbt(ii,JJ) ; dvhdv_tot(ii,JJ) = dvhdv_tot_0(ii,JJ)
2398	726144	vh_err_best(ii,JJ) = abs(vh_err(ii,JJ))
2399		enddo ; enddo
2400
2401	17162	do itt=1,max_itts
2402	6144	select case (itt)
2403	6144	case (:1) ; tol_eta = 1e-6 * CS%tol_eta
2404	5700	case (2) ; tol_eta = 1e-4 * CS%tol_eta
2405	4574	case (3) ; tol_eta = 1e-2 * CS%tol_eta
2406	17162	case default ; tol_eta = CS%tol_eta
2407		end select
2408
2409		!$omp loop collapse(2) private(ii,jj)
2410	2115826	do j=j_start,j_end ; do i=i_start,i_end
2411	2031168	ii=i-i_start+1 ; jj=j-j_start+1
2412	2098664	if (vh_err(ii,JJ) > 0.0) then ; dv_max(ii,JJ) = dv(ii,JJ)
2413	1427166	elseif (vh_err(ii,JJ) < 0.0) then ; dv_min(ii,JJ) = dv(ii,JJ)
2414	845812	else ; do_I(ii,JJ) = .false. ; endif
2415		enddo ; enddo
2416		#ifndef __NVCOMPILER_OPENMP_GPU
2417	17162	domore = .false.
2418		#endif
2419		!$omp loop collapse(2) private(ii,jj,ddv,dv_prev)
2420	2115826	do J=J_start,J_end ; do i=i_start,i_end
2421	2098664	ii=i-i_start+1 ; jj=j-j_start+1 ; if (do_I(ii,JJ)) then
2422	1045099	if ((dt * min(G%IareaT(i,j),G%IareaT(i,j+1))*abs(vh_err(ii,JJ)) > tol_eta) .or. &
2423		(CS%better_iter .and. ((abs(vh_err(ii,JJ)) > tol_vel * dvhdv_tot(ii,JJ)) .or. &
2424		(abs(vh_err(ii,JJ)) > vh_err_best(ii,JJ))) )) then
2425		! Use Newton's method, provided it stays bounded. Otherwise bisect
2426		! the value with the appropriate bound.
2427	580171	ddv = -vh_err(ii,JJ) / dvhdv_tot(ii,JJ)
2428	580171	dv_prev = dv(ii,JJ)
2429	580171	dv(ii,JJ) = dv(ii,JJ) + ddv
2430	580171	if (abs(ddv) < 1.0e-15*abs(dv(ii,JJ))) then
2431	0	do_I(ii,JJ) = .false. ! ddv is small enough to quit.
2432	580171	elseif (ddv > 0.0) then
2433	231787	if (dv(ii,JJ) >= dv_max(ii,JJ)) then
2434	0	dv(ii,JJ) = 0.5*(dv_prev + dv_max(ii,JJ))
2435	0	if (dv_max(ii,JJ) - dv_prev < 1.0e-15*abs(dv(ii,JJ))) do_I(ii,JJ) = .false.
2436		endif
2437		else ! dvv(i,J) < 0.0
2438	348384	if (dv(ii,JJ) <= dv_min(ii,JJ)) then
2439	0	dv(ii,JJ) = 0.5*(dv_prev + dv_min(ii,JJ))
2440	0	if (dv_prev - dv_min(ii,JJ) < 1.0e-15*abs(dv(ii,JJ))) do_I(ii,JJ) = .false.
2441		endif
2442		endif
2443		#ifndef __NVCOMPILER_OPENMP_GPU
2444	580171	if (do_I(ii,JJ)) domore = .true.
2445		#endif
2446		else
2447	464928	do_I(ii,JJ) = .false.
2448		endif
2449		endif ; enddo ; enddo
2450		#ifndef __NVCOMPILER_OPENMP_GPU
2451	17162	if (.not.domore) exit
2452		#endif
2453		!$omp loop collapse(2) private(ii,jj)
2454	1383538	do J=J_start,J_end ; do i=i_start,i_end
2455	1328448	ii=i-i_start+1 ; jj=J-j_start+1
2456	1372520	vh_err(ii,JJ) = -vhbt(ii,JJ) ; dvhdv_tot(ii,JJ) = 0.0
2457		enddo ; enddo
2458	837368	do k=1,nz
2459		!$omp loop collapse(2) private(ii,jj,dvhdv,v_new)
2460	103776368	do J=J_start,J_end ; do i=i_start,i_end
2461	99633600	ii=i-i_start+1 ; jj=j-j_start+1
2462	102939000	if (do_I(ii,JJ)) then
2463	43512825	v_new = v(i,J,k) + dv(ii,JJ) * visc_rem(ii,JJ,k)
2464		call flux_elem(v_new,h_in(i,J,k),h_in(i,J+1,k),h_S(i,J,k),h_S(i,J+1,k),h_N(i,J,k),&
2465		h_N(i,J+1,k),vh_3d(ii,JJ,k),dvhdv,visc_rem(ii,JJ,k),G%dx_Cv(i,J),&
2466		G%IareaT(i,J),G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,&
2467	43512825	por_face_areaV(i,J,k))
2468	43512825	if (local_open_BC) &
2469		call flux_elem_OBC(v_new,h_in(i,J,k),h_in(i,J+1,k),vh_3d(ii,JJ,k),dvhdv,&
2470		visc_rem(ii,JJ,k),por_face_areaV(i,J,k),G%dx_Cv(i,J),OBC,&
2471	0	OBC%segnum_v(i,J))
2472	43512825	vh_err(ii,JJ) = vh_err(ii,JJ) + vh_3d(ii,JJ,k)
2473	43512825	dvhdv_tot(ii,JJ) = dvhdv_tot(ii,JJ) + dvhdv
2474		endif
2475		enddo ; enddo
2476		enddo
2477		!$omp loop collapse(2) private(ii,jj)
2478	1383538	do J=J_start,J_end ; do i=i_start,i_end
2479	1328448	ii=i-i_start+1 ; jj=J-j_start+1
2480	1372520	vh_err_best(ii,JJ) = min(vh_err_best(ii,JJ), abs(vh_err(ii,JJ)))
2481		enddo ; enddo
2482		enddo ! itt-loop
2483		! If there are any faces which have not converged to within the tolerance,
2484		! so-be-it, or else use a final upwind correction?
2485		! This never seems to happen with 20 iterations as max_itt.
2486		!$omp end target
2487
2488		!$omp target exit data &
2489		!$omp map(release:vh_err,vh_err_best,dvhdv_tot,dv_min,dv_max,do_I)
2490
2491	6144	end subroutine meridional_flux_adjust
2492
2493
2494		!> Sets of a structure that describes the meridional barotropic volume or mass fluxes as a
2495		!! function of barotropic flow to agree closely with the sum of the layer's transports.
2496	3072	subroutine set_merid_BT_cont(v, h_in, h_S, h_N, BT_cont, vh_tot_0, dvhdv_tot_0, &
2497	3072	dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
2498	3072	visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
2499	3072	por_face_areaV, niblock, njblock)
2500		integer, intent(in) :: niblock !< i block size for array calculations [nondim].
2501		integer, intent(in) :: njblock !< j block size for array calculations [nondim].
2502		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2503		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2504		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2505		intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
2506		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to
2507		!! calculate fluxes [H ~> m or kg m-2].
2508		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_S !< South edge thickness in the
2509		!! reconstruction [H ~> m or kg m-2].
2510		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_N !< North edge thickness in the
2511		!! reconstruction [H ~> m or kg m-2].
2512		type(BT_cont_type), intent(inout) :: BT_cont !< A structure with elements
2513		!! that describe the effective open face areas as a function of barotropic flow.
2514		real, dimension(niblock,njblock), intent(in) :: vh_tot_0 !< The summed transport
2515		!! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
2516		real, dimension(niblock,njblock), intent(in) :: dvhdv_tot_0 !< The partial derivative
2517		!! of du_err with dv at 0 adjustment [H L ~> m2 or kg m-1].
2518		real, dimension(niblock,njblock), intent(in) :: dv_max_CFL !< Maximum acceptable value
2519		!! of dv [L T-1 ~> m s-1].
2520		real, dimension(niblock,njblock), intent(in) :: dv_min_CFL !< Minimum acceptable value
2521		!! of dv [L T-1 ~> m s-1].
2522		real, intent(in) :: dt !< Time increment [T ~> s].
2523		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2524		type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.
2525		real, dimension(niblock,njblock,SZK_(GV)), &
2526		intent(in) :: visc_rem !< Both the fraction of the
2527		!! momentum originally in a layer that remains after a time-step
2528		!! of viscosity, and the fraction of a time-step's worth of a barotropic
2529		!! acceleration that a layer experiences after viscosity is applied [nondim].
2530		!! Visc_rem is between 0 (at the bottom) and 1 (far above the bottom).
2531		real, dimension(niblock,njblock), intent(in) :: visc_rem_max !< Maximum allowable
2532		!! visc_rem [nondim].
2533		integer, intent(in) :: j_start !< Start of j index range.
2534		integer, intent(in) :: j_end !< End of j index range.
2535		integer, intent(in) :: i_start !< Start of i index range.
2536		integer, intent(in) :: i_end !< End of i index range.
2537		logical, dimension(niblock,njblock), &
2538		intent(in) :: do_I !< A logical flag indicating which I values to work on.
2539		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
2540		intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim]
2541		! Local variables
2542		real, dimension(niblock,njblock) :: &
2543	6144	dv0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
2544	6144	dvL, dvR, & ! The barotropic velocity increments that give the southerly
2545		! (dvL) and northerly (dvR) test velocities [L T-1 ~> m s-1].
2546	6144	dv_CFL, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
2547	6144	FAmt_L, FAmt_R, & ! The summed effective marginal face areas for the 3
2548	6144	FAmt_0, & ! test velocities [H L ~> m2 or kg m-1].
2549	6144	vhtot_L, & ! The summed transport with the southerly (vhtot_L) and
2550	6144	vhtot_R ! and northerly (vhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
2551		real, dimension(niblock,njblock) :: &
2552	6144	zeros ! An array of full of 0 transports [H L2 T-1 ~> m3 s-1 or kg s-1]
2553		real :: &
2554		vh_L, vh_R, & ! The layer transports with the southerly (_L), northerly (_R)
2555		vh_0, & ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
2556		dvhdv_L, & ! The effective layer marginal face areas with the southerly
2557		dvhdv_R, & ! (_L), northerly (_R), and zero-barotropic (_0) test
2558		dvhdv_0, & ! velocities [H L ~> m2 or kg m-1].
2559		v_L, v_R, & ! The southerly (v_L), northerly (v_R), and zero-barotropic
2560		v_0 ! transport (v_0) layer test velocities [L T-1 ~> m s-1].
2561		real :: FA_0 ! The effective face area with 0 barotropic transport [H L ~> m2 or kg m-1].
2562		real :: FA_avg ! The average effective face area [H L ~> m2 or kg m-1], nominally given by
2563		! the realized transport divided by the barotropic velocity.
2564		real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim] This
2565		! limiting is necessary to keep the inverse of visc_rem
2566		! from leading to large CFL numbers.
2567		real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
2568		! in finding the barotropic velocity that changes the
2569		! flow direction [nondim]. This is necessary to keep the inverse
2570		! of visc_rem from leading to large CFL numbers.
2571		real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
2572		! flow is truly upwind [nondim]
2573		real :: Idt ! The inverse of the time step [T-1 ~> s-1].
2574		integer :: i, j, k, nz, ii, jj
2575	6144	real :: vh_tmp(niblock,njblock,SZK_(GV))
2576
2577		!$omp target enter data &
2578		!$omp map(alloc:dv0,dvL,dvR,dv_CFL,FAmt_L,FAmt_R,FAmt_0,vhtot_L,vhtot_R,zeros,vh_tmp)
2579
2580	3072	nz = GV%ke ; Idt = 1.0 / dt
2581	3072	min_visc_rem = 0.1 ; CFL_min = 1e-6
2582
2583		! Diagnose the zero-transport correction, dv0.
2584	95232	do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
2585	446592	zeros(ii,JJ) = 0.0
2586		enddo
2587		call meridional_flux_adjust(v, h_in, h_S, h_N, zeros, vh_tot_0, dvhdv_tot_0, dv0, &
2588		dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
2589		i_start, i_end, j_start, j_end, do_I, por_face_areaV, &
2590	3072	niblock, njblock, vh_tmp)
2591
2592		! Determine the southerly- and northerly- fluxes. Choose a sufficiently
2593		! negative velocity correction for the northerly-flux, and a sufficiently
2594		! positive correction for the southerly-flux.
2595		!$omp target
2596		!$omp loop collapse(2) private(ii,jj)
2597	366144	do J=J_start,J_end ; do i=i_start,i_end
2598	351360	ii=i-i_start+1 ; jj=j-j_start+1
2599	351360	dv_CFL(ii,JJ) = (CFL_min * Idt) * G%dyCv(i,J)
2600	351360	dvR(ii,JJ) = min(0.0,dv0(ii,JJ) - dv_CFL(ii,JJ))
2601	351360	dvL(ii,JJ) = max(0.0,dv0(ii,JJ) + dv_CFL(ii,JJ))
2602	351360	FAmt_L(ii,JJ) = 0.0 ; FAmt_R(ii,JJ) = 0.0 ; FAmt_0(ii,JJ) = 0.0
2603	363072	vhtot_L(ii,JJ) = 0.0 ; vhtot_R(ii,JJ) = 0.0
2604		enddo ; enddo
2605
2606	233472	do k=1,nz
2607		!$omp loop collapse(2) private(ii,jj,visc_rem_lim)
2608	27463872	do J=J_start,J_end ; do i=i_start,i_end
2609	26352000	ii=i-i_start+1 ; jj=j-j_start+1
2610	27230400	if (do_I(ii,jj)) then
2611	26352000	visc_rem_lim = max(visc_rem(ii,JJ,k), min_visc_rem*visc_rem_max(ii,JJ))
2612	26352000	if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
2613	17618400	if (v(i,J,k) + dvR(ii,JJ)visc_rem_lim > -dv_CFL(ii,JJ)visc_rem(ii,JJ,k)) &
2614	331512	dvR(ii,JJ) = -(v(i,J,k) + dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) / visc_rem_lim
2615	17618400	if (v(i,J,k) + dvL(ii,JJ)visc_rem_lim < dv_CFL(ii,JJ)visc_rem(ii,JJ,k)) &
2616	3767578	dvL(ii,JJ) = -(v(i,J,k) - dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) / visc_rem_lim
2617		endif
2618		endif ; enddo ; enddo ; enddo
2619		!$omp end target
2620		!$omp target
2621	233472	do k=1,nz
2622		!$omp loop collapse(2) private(ii,jj,v_L,v_R,v_0,dvhdv_0,dvhdv_L,dvhdv_R,vh_0,vh_L,vh_R)
2623	27463872	do J=J_start,J_end ; do i=i_start,i_end
2624	26352000	ii=i-i_start+1 ; jj=j-j_start+1
2625	27230400	if (do_I(ii,jj)) then
2626	26352000	v_L = v(i,J,k) + dvL(ii,JJ) * visc_rem(ii,JJ,k)
2627	26352000	v_R = v(i,J,k) + dvR(ii,JJ) * visc_rem(ii,JJ,k)
2628	26352000	v_0 = v(i,J,k) + dv0(ii,JJ) * visc_rem(ii,JJ,k)
2629		call flux_elem(v_0,h_in(i,J,k),h_in(i,J+1,k),h_S(i,J,k),h_S(i,J+1,k),h_N(i,J,k),&
2630		h_N(i,J+1,k),vh_0,dvhdv_0,visc_rem(ii,JJ,k),G%dx_Cv(i,J),G%IareaT(i,J),&
2631	26352000	G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,por_face_areaV(i,J,k))
2632		call flux_elem(v_L,h_in(i,J,k),h_in(i,J+1,k),h_S(i,J,k),h_S(i,J+1,k),h_N(i,J,k),&
2633		h_N(i,J+1,k),vh_L,dvhdv_L,visc_rem(ii,JJ,k),G%dx_Cv(i,J),G%IareaT(i,J),&
2634	26352000	G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,por_face_areaV(i,J,k))
2635		call flux_elem(v_R,h_in(i,J,k),h_in(i,J+1,k),h_S(i,J,k),h_S(i,J+1,k),h_N(i,J,k),&
2636		h_N(i,J+1,k),vh_R,dvhdv_R,visc_rem(ii,JJ,k),G%dx_Cv(i,J),G%IareaT(i,J),&
2637	26352000	G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,por_face_areaV(i,J,k))
2638	26352000	FAmt_0(ii,JJ) = FAmt_0(ii,JJ) + dvhdv_0
2639	26352000	FAmt_L(ii,JJ) = FAmt_L(ii,JJ) + dvhdv_L
2640	26352000	FAmt_R(ii,JJ) = FAmt_R(ii,JJ) + dvhdv_R
2641	26352000	vhtot_L(ii,JJ) = vhtot_L(ii,JJ) + vh_L
2642	26352000	vhtot_R(ii,JJ) = vhtot_R(ii,JJ) + vh_R
2643		endif
2644		enddo ; enddo
2645		enddo
2646		!$omp loop collapse(2) private(ii,jj,FA_0,FA_avg)
2647	366144	do J=J_start,J_end ; do i=i_start,i_end
2648	351360	ii=i-i_start+1 ; jj=j-j_start+1
2649	363072	if (do_I(ii,JJ)) then
2650	351360	FA_0 = FAmt_0(ii,JJ) ; FA_avg = FAmt_0(ii,JJ)
2651	351360	if ((dvL(ii,JJ) - dv0(ii,JJ)) /= 0.0) &
2652	351360	FA_avg = vhtot_L(ii,JJ) / (dvL(ii,JJ) - dv0(ii,JJ))
2653	351360	if (FA_avg > max(FA_0, FAmt_L(ii,JJ))) then ; FA_avg = max(FA_0, FAmt_L(ii,JJ))
2654	351354	elseif (FA_avg < min(FA_0, FAmt_L(ii,JJ))) then ; FA_0 = FA_avg ; endif
2655	351360	BT_cont%FA_v_S0(i,J) = FA_0 ; BT_cont%FA_v_SS(i,J) = FAmt_L(ii,JJ)
2656	351360	if (abs(FA_0-FAmt_L(ii,JJ)) <= 1e-12*FA_0) then ; BT_cont%vBT_SS(i,J) = 0.0 ; else
2657		BT_cont%vBT_SS(i,J) = (1.5 * (dvL(ii,JJ) - dv0(ii,JJ))) * &
2658	234912	((FAmt_L(ii,JJ) - FA_avg) / (FAmt_L(ii,JJ) - FA_0))
2659		endif
2660
2661	351360	FA_0 = FAmt_0(ii,JJ) ; FA_avg = FAmt_0(ii,JJ)
2662	351360	if ((dvR(ii,JJ) - dv0(ii,JJ)) /= 0.0) &
2663	351360	FA_avg = vhtot_R(ii,JJ) / (dvR(ii,JJ) - dv0(ii,JJ))
2664	351360	if (FA_avg > max(FA_0, FAmt_R(ii,JJ))) then ; FA_avg = max(FA_0, FAmt_R(ii,JJ))
2665	351312	elseif (FA_avg < min(FA_0, FAmt_R(ii,JJ))) then ; FA_0 = FA_avg ; endif
2666	351360	BT_cont%FA_v_N0(i,J) = FA_0 ; BT_cont%FA_v_NN(i,J) = FAmt_R(ii,JJ)
2667	351360	if (abs(FAmt_R(ii,JJ) - FA_0) <= 1e-12*FA_0) then ; BT_cont%vBT_NN(i,J) = 0.0 ; else
2668		BT_cont%vBT_NN(i,J) = (1.5 * (dvR(ii,JJ) - dv0(ii,JJ))) * &
2669	234911	((FAmt_R(ii,JJ) - FA_avg) / (FAmt_R(ii,JJ) - FA_0))
2670		endif
2671		else
2672	0	BT_cont%FA_v_S0(i,J) = 0.0 ; BT_cont%FA_v_SS(i,J) = 0.0
2673	0	BT_cont%FA_v_N0(i,J) = 0.0 ; BT_cont%FA_v_NN(i,J) = 0.0
2674	0	BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0
2675		endif ; enddo ; enddo
2676		!$omp end target
2677
2678		!$omp target exit data &
2679		!$omp map(release:dv0,dvL,dvR,dv_CFL,FAmt_L,FAmt_R,FAmt_0,vhtot_L,vhtot_R,zeros,vh_tmp)
2680
2681	3072	end subroutine set_merid_BT_cont
2682
2683		!> Calculates left/right edge values for PPM reconstruction.
2684	73	subroutine PPM_reconstruction_x(h_in, h_W, h_E, G, GV, LB, nkblock, h_min, monotonic, simple_2nd, OBC)
2685		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2686		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2687		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2688		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_W !< West edge thickness in the reconstruction,
2689		!! [H ~> m or kg m-2].
2690		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_E !< East edge thickness in the reconstruction,
2691		!! [H ~> m or kg m-2].
2692		type(cont_loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
2693		integer, intent(in) :: nkblock !< k block size for reconstruction calculations [nondim].
2694		real, intent(in) :: h_min !< The minimum thickness
2695		!! that can be obtained by a concave parabolic fit [H ~> m or kg m-2]
2696		logical, intent(in) :: monotonic !< If true, use the
2697		!! Colella & Woodward monotonic limiter.
2698		!! Otherwise use a simple positive-definite limiter.
2699		logical, intent(in) :: simple_2nd !< If true, use the
2700		!! arithmetic mean thicknesses as the default edge values
2701		!! for a simple 2nd order scheme.
2702		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
2703		integer :: k, kk !< vertical grid and k-block indices
2704
2705		! Local variables with useful mnemonic names.
2706		real, dimension(SZI_(G),SZJ_(G),max(1,nkblock)) :: &
2707	146	slp ! The slopes per grid point in a k block [H ~> m or kg m-2]
2708		real, parameter :: oneSixth = 1./6. ! [nondim]
2709		real :: h_ip1, h_im1 ! Neighboring thicknesses or sensibly extrapolated values [H ~> m or kg m-2]
2710		real :: dMx, dMn ! The difference between the local thickness and the maximum (dMx) or
2711		! minimum (dMn) of the surrounding values [H ~> m or kg m-2]
2712		character(len=256) :: mesg
2713		integer :: i, j, isl, iel, jsl, jel, nz, n, stencil, ks, ke
2714		logical :: local_open_BC
2715		type(OBC_segment_type), pointer :: segment => NULL()
2716
2717	73	local_open_BC = .false.
2718	73	if (associated(OBC)) then
2719	0	local_open_BC = OBC%open_u_BCs_exist_globally
2720		endif
2721
2722	73	isl = LB%ish-1 ; iel = LB%ieh+1 ; jsl = LB%jsh ; jel = LB%jeh ; nz = GV%ke
2723
2724		! This is the stencil of the reconstruction, not the scheme overall.
2725	73	stencil = 2 ; if (simple_2nd) stencil = 1
2726
2727	73	if ((isl-stencil < G%isd) .or. (iel+stencil > G%ied)) then
2728		write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
2729		& "x-halo that needs to be increased by ",I0,".")') &
2730	0	stencil + max(G%isd-isl,iel-G%ied)
2731	0	call MOM_error(FATAL,mesg)
2732		endif
2733	73	if ((jsl < G%jsd) .or. (jel > G%jed)) then
2734		write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
2735		& "y-halo that needs to be increased by ",I0,".")') &
2736	0	max(G%jsd-jsl,jel-G%jed)
2737	0	call MOM_error(FATAL,mesg)
2738		endif
2739
2740		!$omp target enter data map(alloc: slp)
2741
2742	5548	do ks = 1, nz, nkblock
2743	5475	ke = min(ks + nkblock - 1, nz)
2744
2745	5475	if (simple_2nd) then
2746		! untested
2747	0	do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_im1,h_ip1))
2748	0	h_im1 = G%mask2dT(i-1,j) * h_in(i-1,j,k) + (1.0-G%mask2dT(i-1,j)) * h_in(i,j,k)
2749	0	h_ip1 = G%mask2dT(i+1,j) * h_in(i+1,j,k) + (1.0-G%mask2dT(i+1,j)) * h_in(i,j,k)
2750	0	h_W(i,j,k) = 0.5*( h_im1 + h_in(i,j,k) )
2751	0	h_E(i,j,k) = 0.5*( h_ip1 + h_in(i,j,k) )
2752		enddo
2753		else
2754	5475	do concurrent (k=ks:ke, j=jsl:jel, i=isl-1:iel+1) DO_LOCALITY(local(dMx,dMn,kk))
2755	44807400	kk = k - ks + 1
2756	90299175	if ((G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j)) == 0.0) then
2757	17881350	slp(i,j,kk) = 0.0
2758		else
2759		! This uses a simple 2nd order slope.
2760	26926050	slp(i,j,kk) = 0.5 * (h_in(i+1,j,k) - h_in(i-1,j,k))
2761		! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
2762	26926050	dMx = max(h_in(i+1,j,k), h_in(i-1,j,k), h_in(i,j,k)) - h_in(i,j,k)
2763	26926050	dMn = h_in(i,j,k) - min(h_in(i+1,j,k), h_in(i-1,j,k), h_in(i,j,k))
2764	26926050	slp(i,j,kk) = sign(1.,slp(i,j,kk)) * min(abs(slp(i,j,kk)), 2. * min(dMx, dMn))
2765		! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
2766		endif
2767		enddo
2768
2769	5475	if (local_open_BC) then
2770		! untested
2771	0	do n=1, OBC%number_of_segments
2772	0	segment => OBC%segment(n)
2773	0	if (.not. segment%on_pe) cycle
2774	0	if (segment%is_E_or_W) then
2775	0	I=segment%HI%IsdB
2776	0	do concurrent (k=ks:ke, j=segment%HI%jsd:segment%HI%jed) DO_LOCALITY(local(kk))
2777	0	kk = k - ks + 1
2778	0	slp(i+1,j,kk) = 0.0
2779	0	slp(i,j,kk) = 0.0
2780		enddo
2781		endif
2782		enddo
2783		endif
2784
2785	5475	do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_im1,h_ip1,kk))
2786	44084700	kk = k - ks + 1
2787		! Neighboring values should take into account any boundaries. The 3
2788		! following sets of expressions are equivalent.
2789		! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
2790		! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
2791	44084700	h_im1 = G%mask2dT(i-1,j) * h_in(i-1,j,k) + (1.0-G%mask2dT(i-1,j)) * h_in(i,j,k)
2792	44084700	h_ip1 = G%mask2dT(i+1,j) * h_in(i+1,j,k) + (1.0-G%mask2dT(i+1,j)) * h_in(i,j,k)
2793		! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2794		h_W(i,j,k) = 0.5*( h_im1 + h_in(i,j,k) ) + &
2795	44084700	oneSixth*( slp(i-1,j,kk) - slp(i,j,kk) )
2796		h_E(i,j,k) = 0.5*( h_ip1 + h_in(i,j,k) ) + &
2797	88842825	oneSixth*( slp(i,j,kk) - slp(i+1,j,kk) )
2798		enddo
2799		endif
2800
2801	5475	if (local_open_BC) then
2802		! untested
2803	0	do n=1, OBC%number_of_segments
2804	0	segment => OBC%segment(n)
2805	0	if (.not. segment%on_pe) cycle
2806	0	if (segment%direction == OBC_DIRECTION_E) then
2807	0	I=segment%HI%IsdB
2808	0	do concurrent (k=ks:ke, j=segment%HI%jsd:segment%HI%jed)
2809	0	h_W(i+1,j,k) = h_in(i,j,k)
2810	0	h_E(i+1,j,k) = h_in(i,j,k)
2811	0	h_W(i,j,k) = h_in(i,j,k)
2812	0	h_E(i,j,k) = h_in(i,j,k)
2813		enddo
2814	0	elseif (segment%direction == OBC_DIRECTION_W) then
2815	0	I=segment%HI%IsdB
2816	0	do concurrent (k=ks:ke, j=segment%HI%jsd:segment%HI%jed)
2817	0	h_W(i,j,k) = h_in(i+1,j,k)
2818	0	h_E(i,j,k) = h_in(i+1,j,k)
2819	0	h_W(i+1,j,k) = h_in(i+1,j,k)
2820	0	h_E(i+1,j,k) = h_in(i+1,j,k)
2821		enddo
2822		endif
2823		enddo
2824		endif
2825
2826	5548	if (monotonic) then
2827		! untested
2828	0	call PPM_limit_CW84(h_in, h_W, h_E, G, GV, isl, iel, jsl, jel, ks, ke)
2829		else
2830	5475	call PPM_limit_pos(h_in, h_W, h_E, h_min, G, GV, isl, iel, jsl, jel, ks, ke)
2831		endif
2832		enddo
2833
2834		!$omp target exit data map(release: slp)
2835
2836	73	end subroutine PPM_reconstruction_x
2837
2838		!> Calculates left/right edge values for PPM reconstruction.
2839	73	subroutine PPM_reconstruction_y(h_in, h_S, h_N, G, GV, LB, nkblock, h_min, monotonic, simple_2nd, OBC)
2840		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2841		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2842		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2843		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_S !< South edge thickness in the reconstruction,
2844		!! [H ~> m or kg m-2].
2845		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_N !< North edge thickness in the reconstruction,
2846		!! [H ~> m or kg m-2].
2847		type(cont_loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
2848		integer, intent(in) :: nkblock !< k block size for reconstruction calculations [nondim].
2849		real, intent(in) :: h_min !< The minimum thickness
2850		!! that can be obtained by a concave parabolic fit [H ~> m or kg m-2]
2851		logical, intent(in) :: monotonic !< If true, use the
2852		!! Colella & Woodward monotonic limiter.
2853		!! Otherwise use a simple positive-definite limiter.
2854		logical, intent(in) :: simple_2nd !< If true, use the
2855		!! arithmetic mean thicknesses as the default edge values
2856		!! for a simple 2nd order scheme.
2857		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
2858		integer :: k, kk !< vertical grid and k-block indices
2859
2860		! Local variables with useful mnemonic names.
2861		real, dimension(SZI_(G),SZJ_(G),max(1,nkblock)) :: &
2862	146	slp ! The slopes per grid point in a k block [H ~> m or kg m-2]
2863		real, parameter :: oneSixth = 1./6. ! [nondim]
2864		real :: h_jp1, h_jm1 ! Neighboring thicknesses or sensibly extrapolated values [H ~> m or kg m-2]
2865		real :: dMx, dMn ! The difference between the local thickness and the maximum (dMx) or
2866		! minimum (dMn) of the surrounding values [H ~> m or kg m-2]
2867		character(len=256) :: mesg
2868		integer :: i, j, isl, iel, jsl, jel, nz, n, stencil, ks, ke
2869		logical :: local_open_BC
2870		type(OBC_segment_type), pointer :: segment => NULL()
2871
2872	73	local_open_BC = .false.
2873	73	if (associated(OBC)) then
2874	0	local_open_BC = OBC%open_v_BCs_exist_globally
2875		endif
2876
2877	73	isl = LB%ish ; iel = LB%ieh ; jsl = LB%jsh-1 ; jel = LB%jeh+1 ; nz = G%ke
2878
2879		! This is the stencil of the reconstruction, not the scheme overall.
2880	73	stencil = 2 ; if (simple_2nd) stencil = 1
2881
2882	73	if ((isl < G%isd) .or. (iel > G%ied)) then
2883		write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
2884		& "x-halo that needs to be increased by ",I0,".")') &
2885	0	max(G%isd-isl,iel-G%ied)
2886	0	call MOM_error(FATAL,mesg)
2887		endif
2888	73	if ((jsl-stencil < G%jsd) .or. (jel+stencil > G%jed)) then
2889		write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
2890		& "y-halo that needs to be increased by ",I0,".")') &
2891	0	stencil + max(G%jsd-jsl,jel-G%jed)
2892	0	call MOM_error(FATAL,mesg)
2893		endif
2894
2895		!$omp target enter data map(alloc: slp)
2896
2897	5548	do ks = 1, nz, nkblock
2898	5475	ke = min(ks + nkblock - 1, nz)
2899
2900	5475	if (simple_2nd) then
2901		! untested
2902	0	do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_jm1,h_jp1))
2903	0	h_jm1 = G%mask2dT(i,j-1) * h_in(i,j-1,k) + (1.0-G%mask2dT(i,j-1)) * h_in(i,j,k)
2904	0	h_jp1 = G%mask2dT(i,j+1) * h_in(i,j+1,k) + (1.0-G%mask2dT(i,j+1)) * h_in(i,j,k)
2905	0	h_S(i,j,k) = 0.5*( h_jm1 + h_in(i,j,k) )
2906	0	h_N(i,j,k) = 0.5*( h_jp1 + h_in(i,j,k) )
2907		enddo
2908		else
2909	662475	do concurrent (k=ks:ke, j=jsl-1:jel+1, i=isl:iel) DO_LOCALITY(local(dMx,dMn,kk))
2910	42048000	kk = k - ks + 1
2911	84758475	if ((G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1)) == 0.0) then
2912	15910350	slp(i,j,kk) = 0.0
2913		else
2914		! This uses a simple 2nd order slope.
2915	26137650	slp(i,j,kk) = 0.5 * (h_in(i,j+1,k) - h_in(i,j-1,k))
2916		! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
2917	26137650	dMx = max(h_in(i,j+1,k), h_in(i,j-1,k), h_in(i,j,k)) - h_in(i,j,k)
2918	26137650	dMn = h_in(i,j,k) - min(h_in(i,j+1,k), h_in(i,j-1,k), h_in(i,j,k))
2919	26137650	slp(i,j,kk) = sign(1.,slp(i,j,kk)) * min(abs(slp(i,j,kk)), 2. * min(dMx, dMn))
2920		! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
2921		endif
2922		enddo
2923
2924	5475	if (local_open_BC) then
2925		! untested
2926	0	do n=1, OBC%number_of_segments
2927	0	segment => OBC%segment(n)
2928	0	if (.not. segment%on_pe) cycle
2929	0	if (segment%is_N_or_S) then
2930	0	J=segment%HI%JsdB
2931	0	do concurrent (k=ks:ke, i=segment%HI%isd:segment%HI%ied) DO_LOCALITY(local(kk))
2932	0	kk = k - ks + 1
2933	0	slp(i,j+1,kk) = 0.0
2934	0	slp(i,j,kk) = 0.0
2935		enddo
2936		endif
2937		enddo
2938		endif
2939
2940	5475	do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_jm1,h_jp1,kk))
2941	40734000	kk = k - ks + 1
2942		! Neighboring values should take into account any boundaries. The 3
2943		! following sets of expressions are equivalent.
2944	40734000	h_jm1 = G%mask2dT(i,j-1) * h_in(i,j-1,k) + (1.0-G%mask2dT(i,j-1)) * h_in(i,j,k)
2945	40734000	h_jp1 = G%mask2dT(i,j+1) * h_in(i,j+1,k) + (1.0-G%mask2dT(i,j+1)) * h_in(i,j,k)
2946		! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2947		h_S(i,j,k) = 0.5*( h_jm1 + h_in(i,j,k) ) + &
2948	40734000	oneSixth*( slp(i,j-1,kk) - slp(i,j,kk) )
2949		h_N(i,j,k) = 0.5*( h_jp1 + h_in(i,j,k) ) + &
2950	82130475	oneSixth*( slp(i,j,kk) - slp(i,j+1,kk) )
2951		enddo
2952		endif
2953
2954	5475	if (local_open_BC) then
2955		! untested
2956	0	do n=1, OBC%number_of_segments
2957	0	segment => OBC%segment(n)
2958	0	if (.not. segment%on_pe) cycle
2959	0	if (segment%direction == OBC_DIRECTION_N) then
2960	0	J=segment%HI%JsdB
2961	0	do concurrent (k=ks:ke, i=segment%HI%isd:segment%HI%ied)
2962	0	h_S(i,j+1,k) = h_in(i,j,k)
2963	0	h_N(i,j+1,k) = h_in(i,j,k)
2964	0	h_S(i,j,k) = h_in(i,j,k)
2965	0	h_N(i,j,k) = h_in(i,j,k)
2966		enddo
2967	0	elseif (segment%direction == OBC_DIRECTION_S) then
2968	0	J=segment%HI%JsdB
2969	0	do concurrent (k=ks:ke, i=segment%HI%isd:segment%HI%ied)
2970	0	h_S(i,j,k) = h_in(i,j+1,k)
2971	0	h_N(i,j,k) = h_in(i,j+1,k)
2972	0	h_S(i,j+1,k) = h_in(i,j+1,k)
2973	0	h_N(i,j+1,k) = h_in(i,j+1,k)
2974		enddo
2975		endif
2976		enddo
2977		endif
2978
2979	5548	if (monotonic) then
2980		! untested
2981	0	call PPM_limit_CW84(h_in, h_S, h_N, G, GV, isl, iel, jsl, jel, ks, ke)
2982		else
2983	5475	call PPM_limit_pos(h_in, h_S, h_N, h_min, G, GV, isl, iel, jsl, jel, ks, ke)
2984		endif
2985		enddo
2986
2987		!$omp target exit data map(release: slp)
2988
2989	73	end subroutine PPM_reconstruction_y
2990
2991		!> This subroutine limits the left/right edge values of the PPM reconstruction
2992		!! to give a reconstruction that is positive-definite. Here this is
2993		!! reinterpreted as giving a constant thickness if the mean thickness is less
2994		!! than h_min, with a minimum of h_min otherwise.
2995	10950	subroutine PPM_limit_pos(h_in, h_L, h_R, h_min, G, GV, iis, iie, jis, jie, ks, ke)
2996		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2997		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
2998		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2999		intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
3000		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3001		intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
3002		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3003		intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
3004		real, intent(in) :: h_min !< The minimum thickness
3005		!! that can be obtained by a concave parabolic fit [H ~> m or kg m-2]
3006		integer, intent(in) :: iis !< Start of i index range.
3007		integer, intent(in) :: iie !< End of i index range.
3008		integer, intent(in) :: jis !< Start of j index range.
3009		integer, intent(in) :: jie !< End of j index range.
3010		integer, intent(in) :: ks !< Start of k index range.
3011		integer, intent(in) :: ke !< End of k index range.
3012
3013		! Local variables
3014		real :: curv ! The grid-normalized curvature of the three thicknesses [H ~> m or kg m-2]
3015		real :: dh ! The difference between the edge thicknesses [H ~> m or kg m-2]
3016		real :: scale ! A scaling factor to reduce the curvature of the fit [nondim]
3017		integer :: i,j,k
3018
3019	86154600	do concurrent (k=ks:ke, j=jis:jie, i=iis:iie)
3020		! This limiter prevents undershooting minima within the domain with
3021		! values less than h_min.
3022	84818700	curv = 3.0((h_L(i,j,k) + h_R(i,j,k)) - 2.0h_in(i,j,k))
3023	170973300	if (curv > 0.0) then ! Only minima are limited.
3024	25644351	dh = h_R(i,j,k) - h_L(i,j,k)
3025	25644351	if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
3026	18775254	if (h_in(i,j,k) <= h_min) then
3027	0	h_L(i,j,k) = h_in(i,j,k) ; h_R(i,j,k) = h_in(i,j,k)
3028	18775254	elseif (12.0curv(h_in(i,j,k) - h_min) < (curv*2 + 3.0dh**2)) then
3029		! The minimum value is h_in - (curv^2 + 3dh^2)/(12curv), and must
3030		! be limited in this case. 0 < scale < 1.
3031	2160655	scale = 12.0curv(h_in(i,j,k) - h_min) / (curv*2 + 3.0dh**2)
3032	2160655	h_L(i,j,k) = h_in(i,j,k) + scale*(h_L(i,j,k) - h_in(i,j,k))
3033	2160655	h_R(i,j,k) = h_in(i,j,k) + scale*(h_R(i,j,k) - h_in(i,j,k))
3034		endif
3035		endif
3036		endif
3037		enddo
3038
3039	10950	end subroutine PPM_limit_pos
3040
3041		!> This subroutine limits the left/right edge values of the PPM reconstruction
3042		!! according to the monotonic prescription of Colella and Woodward, 1984.
3043	0	subroutine PPM_limit_CW84(h_in, h_L, h_R, G, GV, iis, iie, jis, jie, ks, ke)
3044		type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
3045		type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
3046		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
3047		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_L !< Left thickness in the reconstruction,
3048		!! [H ~> m or kg m-2].
3049		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_R !< Right thickness in the reconstruction,
3050		!! [H ~> m or kg m-2].
3051		integer, intent(in) :: iis !< Start of i index range.
3052		integer, intent(in) :: iie !< End of i index range.
3053		integer, intent(in) :: jis !< Start of j index range.
3054		integer, intent(in) :: jie !< End of j index range.
3055		integer, intent(in) :: ks !< Start of k index range.
3056		integer, intent(in) :: ke !< End of k index range.
3057
3058		! Local variables
3059		real :: h_i ! A copy of the cell-average layer thickness [H ~> m or kg m-2]
3060		real :: RLdiff ! The difference between the input edge values [H ~> m or kg m-2]
3061		real :: RLdiff2 ! The squared difference between the input edge values [H2 ~> m2 or kg2 m-4]
3062		real :: RLmean ! The average of the input edge thicknesses [H ~> m or kg m-2]
3063		real :: FunFac ! A curious product of the thickness slope and curvature [H2 ~> m2 or kg2 m-4]
3064		integer :: i, j, k
3065
3066		! untested
3067	0	do concurrent (k=ks:ke, j=jis:jie, i=iis:iie) DO_LOCALITY(local(h_i,RLdiff,RLdiff2,RLmean,FunFac))
3068		! This limiter monotonizes the parabola following
3069		! Colella and Woodward, 1984, Eq. 1.10
3070	0	h_i = h_in(i,j,k)
3071	0	if ( ( h_R(i,j,k) - h_i ) * ( h_i - h_L(i,j,k) ) <= 0. ) then
3072	0	h_L(i,j,k) = h_i ; h_R(i,j,k) = h_i
3073		else
3074	0	RLdiff = h_R(i,j,k) - h_L(i,j,k) ! Difference of edge values
3075	0	RLmean = 0.5 * ( h_R(i,j,k) + h_L(i,j,k) ) ! Mean of edge values
3076	0	FunFac = 6. * RLdiff * ( h_i - RLmean ) ! Some funny factor
3077	0	RLdiff2 = RLdiff * RLdiff ! Square of difference
3078	0	if ( FunFac > RLdiff2 ) h_L(i,j,k) = 3. * h_i - 2. * h_R(i,j,k)
3079	0	if ( FunFac < -RLdiff2 ) h_R(i,j,k) = 3. * h_i - 2. * h_L(i,j,k)
3080		endif
3081		enddo
3082
3083	0	end subroutine PPM_limit_CW84
3084
3085		!> Return the maximum ratio of a/b or maxrat.
3086	0	pure function ratio_max(a, b, maxrat) result(ratio)
3087		real, intent(in) :: a !< Numerator, in arbitrary units [A]
3088		real, intent(in) :: b !< Denominator, in arbitrary units [B]
3089		real, intent(in) :: maxrat !< Maximum value of ratio [A B-1]
3090		real :: ratio !< Return value [A B-1]
3091
3092	0	if (abs(a) > abs(maxrat*b)) then
3093	0	ratio = maxrat
3094		else
3095	0	ratio = a / b
3096		endif
3097	0	end function ratio_max
3098
3099		!> Initializes continuity_ppm_cs
3100	1	subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS, OBC)
3101		type(time_type), target, intent(in) :: Time !< The current model time.
3102		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3103		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
3104		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3105		type(param_file_type), intent(in) :: param_file !< A structure indicating
3106		!! the open file to parse for model parameter values.
3107		type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
3108		!! regulate diagnostic output.
3109		type(continuity_PPM_CS), intent(inout) :: CS !< Module's control structure.
3110		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
3111		logical :: local_open_BC, use_h_marg_min
3112		type(OBC_segment_type), pointer :: segment => NULL()
3113		integer :: n
3114
3115		!> This include declares and sets the variable "version".
3116		# include "version_variable.h"
3117		character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name.
3118		character(len=256) :: mesg
3119		character(len=10) :: niblock_dflt_str, njblock_dflt_str, nkblock_dflt_str
3120		#ifdef __NVCOMPILER_OPENMP_GPU
3121		integer, parameter :: default_niblock = 0 !< Default i block size for array calculations [nondim].
3122		integer, parameter :: default_njblock = 0 !< Default j block size for array calculations [nondim].
3123		integer, parameter :: default_nkblock = 0 !< Default k block size for reconstruction calculations [nondim].
3124		#else
3125		! These were found to give best performance in limited tests.
3126		integer, parameter :: default_niblock = 32 !< Default i block size for array calculations [nondim].
3127		integer, parameter :: default_njblock = 4 !< Default j block size for array calculations [nondim].
3128		integer, parameter :: default_nkblock = 1 !< Default k block size for reconstruction calculations [nondim].
3129		#endif
3130
3131	1	CS%initialized = .true.
3132
3133	1	local_open_BC = .false.
3134	1	if (associated(OBC)) then
3135	0	local_open_BC = OBC%open_u_BCs_exist_globally
3136		endif
3137
3138		! Read all relevant parameters and write them to the model log.
3139	1	call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.true.)
3140		call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", CS%monotonic, &
3141		"If true, CONTINUITY_PPM uses the Colella and Woodward "//&
3142		"monotonic limiter. The default (false) is to use a "//&
3143	1	"simple positive definite limiter.", default=.false.)
3144		call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", CS%simple_2nd, &
3145		"If true, CONTINUITY_PPM uses a simple 2nd order "//&
3146		"(arithmetic mean) interpolation of the edge values. "//&
3147		"This may give better PV conservation properties. While "//&
3148		"it formally reduces the accuracy of the continuity "//&
3149		"solver itself in the strongly advective limit, it does "//&
3150		"not reduce the overall order of accuracy of the dynamic "//&
3151	1	"core.", default=.false.)
3152		call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", CS%upwind_1st, &
3153		"If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
3154		"continuity solver. This scheme is highly diffusive "//&
3155		"but may be useful for debugging or in single-column "//&
3156	1	"mode where its minimal stencil is useful.", default=.false.)
3157		call get_param(param_file, mdl, "ETA_TOLERANCE", CS%tol_eta, &
3158		"The tolerance for the differences between the "//&
3159		"barotropic and baroclinic estimates of the sea surface "//&
3160		"height due to the fluxes through each face. The total "//&
3161		"tolerance for SSH is 4 times this value. The default "//&
3162		"is 0.5NKANGSTROM, and this should not be set less "//&
3163		"than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=GV%m_to_H, &
3164	1	default=0.5GV%keGV%Angstrom_m)
3165
3166		call get_param(param_file, mdl, "VELOCITY_TOLERANCE", CS%tol_vel, &
3167		"The tolerance for barotropic velocity discrepancies "//&
3168		"between the barotropic solution and the sum of the "//&
3169	1	"layer thicknesses.", units="m s-1", default=3.0e8, scale=US%m_s_to_L_T)
3170		! The speed of light is the default.
3171
3172		call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", CS%aggress_adjust,&
3173		"If true, allow the adjusted velocities to have a "//&
3174	1	"relative CFL change up to 0.5.", default=.false.)
3175	1	CS%vol_CFL = CS%aggress_adjust
3176		call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", CS%vol_CFL, &
3177		"If true, use the ratio of the open face lengths to the "//&
3178		"tracer cell areas when estimating CFL numbers. The "//&
3179		"default is set by CONT_PPM_AGGRESS_ADJUST.", &
3180	1	default=CS%aggress_adjust, do_not_read=CS%aggress_adjust)
3181		call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", CS%CFL_limit_adjust, &
3182		"The maximum CFL of the adjusted velocities.", units="nondim", &
3183	1	default=0.5)
3184		call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", CS%better_iter, &
3185		"If true, stop corrective iterations using a velocity "//&
3186		"based criterion and only stop if the iteration is "//&
3187	1	"better than all predecessors.", default=.true.)
3188		call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", CS%use_visc_rem_max, &
3189		"If true, use more appropriate limiting bounds for "//&
3190	1	"corrections in strongly viscous columns.", default=.true.)
3191		call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", CS%marginal_faces, &
3192		"If true, use the marginal face areas from the continuity "//&
3193		"solver for use as the weights in the barotropic solver. "//&
3194	1	"Otherwise use the transport averaged areas.", default=.true.)
3195		call get_param(param_file, mdl, "CONT_USE_H_MARG_MIN", use_h_marg_min, &
3196		"If true, the marginal thickness used and returned from continuity "//&
3197		"is bounded from below by a sub-roundoff value. Otherwise the "//&
3198	1	"minimum is 0.", default=.false.)
3199	1	write(niblock_dflt_str, '(I0)') default_niblock
3200	1	write(njblock_dflt_str, '(I0)') default_njblock
3201	1	write(nkblock_dflt_str, '(I0)') default_nkblock
3202		call get_param(param_file, mdl, "CONTINUITY_NIBLOCK", CS%niblock, &
3203		"The i-direction block size used in the mass and volume flux calculations. "//&
3204		"the default 0 setting is dynamic and fits the "//&
3205	1	"full computational i-domain length.", default=default_niblock, layoutParam=.true.)
3206		call get_param(param_file, mdl, "CONTINUITY_NJBLOCK", CS%njblock, &
3207		"The j-direction block size used in the mass and volume flux calculations. "//&
3208		"the default 0 setting is dynamic and fits the "//&
3209	1	"full computational j-domain length.", default=default_njblock, layoutParam=.true.)
3210		call get_param(param_file, mdl, "CONTINUITY_NKBLOCK", CS%nkblock, &
3211		"The k-direction block size used in PPM reconstruction edge value calculations. "//&
3212		"the default 0 setting is dynamic and fits the "//&
3213	1	"full vertical column.", default=default_nkblock, layoutParam=.true.)
3214	1	if (CS%niblock < 0) &
3215		call MOM_error(FATAL, "CONTINUITY_NIBLOCK must be nonnegative; "//&
3216	0	"use 0 to select the default block size.")
3217	1	if (CS%njblock < 0) &
3218		call MOM_error(FATAL, "CONTINUITY_NJBLOCK must be nonnegative; "//&
3219	0	"use 0 to select the default block size.")
3220	1	if (CS%nkblock < 0) &
3221		call MOM_error(FATAL, "CONTINUITY_NKBLOCK must be nonnegative; "//&
3222	0	"use 0 to select the default block size.")
3223	1	CS%diag => diag
3224		!$omp target update to(CS)
3225
3226	1	id_clock_reconstruct = cpu_clock_id('(Ocean continuity reconstruction)', grain=CLOCK_ROUTINE)
3227	1	id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=CLOCK_ROUTINE)
3228	1	id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=CLOCK_ROUTINE)
3229
3230	1	if (use_h_marg_min) then
3231	0	CS%h_marg_min = GV%H_subroundoff
3232		else
3233	1	CS%h_marg_min = 0.
3234		endif
3235
3236	1	if (local_open_BC) then
3237	0	do n=1, OBC%number_of_segments
3238	0	segment => OBC%segment(n)
3239	0	if (associated(segment%h_Reg)) then
3240	0	if (.not. allocated(segment%h_Reg%h_res)) then
3241		write(mesg,'("In MOM_continuity_PPM, continuity_PPM_init called with ", &
3242	0	& "badly configured h_res.")')
3243	0	call MOM_error(FATAL, mesg)
3244		endif
3245		endif
3246		enddo
3247		endif
3248
3249	1	end subroutine continuity_PPM_init
3250
3251		!> continuity_PPM_stencil returns the continuity solver stencil size
3252	171	function continuity_PPM_stencil(CS) result(stencil)
3253		type(continuity_PPM_CS), intent(in) :: CS !< Module's control structure.
3254		integer :: stencil !< The continuity solver stencil size with the current settings.
3255
3256	171	stencil = 3 ; if (CS%simple_2nd) stencil = 2 ; if (CS%upwind_1st) stencil = 1
3257
3258	171	end function continuity_PPM_stencil
3259
3260		!> Set up a structure that stores the sizes of the i- and j-loops to to work on in the continuity solver.
3261	146	function set_continuity_loop_bounds(G, CS, i_stencil, j_stencil) result(LB)
3262		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3263		type(continuity_PPM_CS), intent(in) :: CS !< Module's control structure.
3264		logical, optional, intent(in) :: i_stencil !< If present and true, extend the i-loop bounds
3265		!! by the stencil width of the continuity scheme.
3266		logical, optional, intent(in) :: j_stencil !< If present and true, extend the j-loop bounds
3267		!! by the stencil width of the continuity scheme.
3268		type(cont_loop_bounds_type) :: LB !< A type storing the array sizes to work on in the continuity routines.
3269
3270		! Local variables
3271		logical :: add_i_stencil, add_j_stencil ! Local variables set based on i_stencil and j_stensil
3272		integer :: stencil ! The continuity solver stencil size with the current continuity scheme.
3273
3274	146	add_i_stencil = .false. ; if (present(i_stencil)) add_i_stencil = i_stencil
3275	146	add_j_stencil = .false. ; if (present(j_stencil)) add_j_stencil = j_stencil
3276
3277	146	stencil = continuity_PPM_stencil(CS)
3278
3279	146	if (add_i_stencil) then
3280	0	LB%ish = G%isc-stencil ; LB%ieh = G%iec+stencil
3281		else
3282	146	LB%ish = G%isc ; LB%ieh = G%iec
3283		endif
3284
3285	146	if (add_j_stencil) then
3286	73	LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil
3287		else
3288	73	LB%jsh = G%jsc ; LB%jeh = G%jec
3289		endif
3290
3291	146	end function set_continuity_loop_bounds
3292
3293		!> \namespace mom_continuity_ppm
3294		!!
3295		!! This module contains the subroutines that advect layer
3296		!! thickness. The scheme here uses a Piecewise-Parabolic method with
3297		!! a positive definite limiter.
3298
3299	0	end module MOM_continuity_PPM