← back to index

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.
8module MOM_continuity_PPM
9
10use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
11use MOM_diag_mediator, only : time_type, diag_ctrl
12use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
13use MOM_file_parser, only : get_param, log_version, param_file_type
14use MOM_grid, only : ocean_grid_type
15use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE
16use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
17use MOM_unit_scaling, only : unit_scale_type
18use MOM_variables, only : BT_cont_type, porous_barrier_type
19use MOM_verticalGrid, only : verticalGrid_type
20
21!$ use omp_lib, only: omp_get_num_devices
22
23implicit none ; private
24
25#include <MOM_memory.h>
26
27public continuity_PPM, continuity_PPM_init, continuity_PPM_stencil
28public continuity_fluxes, continuity_adjust_vel
29public zonal_mass_flux, meridional_mass_flux
30public zonal_edge_thickness, meridional_edge_thickness
31public continuity_zonal_convergence, continuity_merdional_convergence
32public zonal_flux_thickness, meridional_flux_thickness
33public zonal_BT_mass_flux, meridional_BT_mass_flux
34public set_continuity_loop_bounds
35
36!>@{ CPU time clock IDs
37integer :: id_clock_reconstruct, id_clock_update, id_clock_correct
38!>@}
39
40!> Control structure for mom_continuity_ppm
41type, 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
80end type continuity_PPM_CS
81
82!> A container for loop bounds
83type, public :: cont_loop_bounds_type ; private
84 !>@{ Loop bounds
85 integer :: ish, ieh, jsh, jeh
86 !>@}
87end 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.
91interface continuity_fluxes
92 module procedure continuity_3d_fluxes, continuity_2d_fluxes
93end interface continuity_fluxes
94
95contains
96
97!> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme,
98!! based on Lin (1994).
9916subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &
10040 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, u*h*dy [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, v*h*dx [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
15726 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
15826 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
15926 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
16026 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
16713 niblock = CS%niblock
16813 njblock = CS%njblock
169
17013 h_min = GV%Angstrom_H
171
17213 if (.not.CS%initialized) call MOM_error(FATAL, &
1730 "MOM_continuity_PPM: Module must be initialized before it is used.")
174
17513 x_first = (MOD(G%first_direction,2) == 0)
176
17713 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 "// &
1790 "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
18313 if (x_first) then
184
185 ! First advect zonally, with loop bounds that accomodate the subsequent meridional advection.
18613 LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.true.)
187 ! set whole-domain block sizes when ni/jblock is 0
18813 if (niblock == 0) niblock = LB%ieh-LB%ish+2
18913 if (njblock == 0) njblock = LB%jeh-LB%jsh+1
190
19113 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, &
19313 niblock, njblock, LB, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)
19413 call continuity_zonal_convergence(h, uh, dt, G, GV, LB, hin)
195
196 ! Now advect meridionally, using the updated thicknesses to determine the fluxes.
19713 LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
198
19913 if (niblock == 0) niblock = LB%ieh-LB%ish+1
20013 if (njblock == 0) njblock = LB%jeh-LB%jsh+2
201
20213 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, &
20413 niblock, njblock, LB, vhbt, visc_rem_v, v_cor, BT_cont, dv_cor)
20513 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.
2100 LB = set_continuity_loop_bounds(G, CS, i_stencil=.true., j_stencil=.false.)
2110 if (niblock == 0) niblock = LB%ieh-LB%ish+1
2120 if (njblock == 0) njblock = LB%jeh-LB%jsh+2
2130 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, &
2150 niblock, njblock, LB, vhbt, visc_rem_v, v_cor, BT_cont, dv_cor)
2160 call continuity_merdional_convergence(h, vh, dt, G, GV, LB, hin)
217
218 ! Now advect zonally, using the updated thicknesses to determine the fluxes.
2190 LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
2200 if (niblock == 0) niblock = LB%ieh-LB%ish+2
2210 if (njblock == 0) njblock = LB%jeh-LB%jsh+1
2220 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, &
2240 niblock, njblock, LB, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)
2250 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
23069end 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.
2360subroutine 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 !! u*h*dy [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 !! v*h*dx [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
2580 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
2590 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
2600 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
2610 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
2650 niblock = CS%niblock
2660 njblock = CS%njblock
267
2680 call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC)
2690 if (niblock == 0) niblock = G%iec-G%isc+2
2700 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, &
2720 niblock=niblock, njblock=njblock)
273
2740 call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC)
2750 if (niblock == 0) niblock = G%iec-G%isc+1
2760 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, &
2780 niblock=niblock, njblock=njblock)
279
2800end 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.
2860subroutine 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
3080 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
3090 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
3100 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
3110 real :: h_N(SZI_(G),SZJ_(G),SZK_(GV)) ! North edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
312
3130 call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC)
3140 call zonal_BT_mass_flux(u, h, h_W, h_E, uhbt, dt, G, GV, US, CS, OBC, pbv%por_face_areaU)
315
3160 call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC)
3170 call meridional_BT_mass_flux(v, h, h_S, h_N, vhbt, dt, G, GV, US, CS, OBC, pbv%por_face_areaV)
318
3190end 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.
3230subroutine 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
3670 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_in !< Input zonal velocity [L T-1 ~> m s-1]
3680 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_in !< Input meridional velocity [L T-1 ~> m s-1]
3690 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uh !< Volume flux through zonal faces =
370 !! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
3710 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vh !< Volume flux through meridional faces =
372 !! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].
3730 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
3740 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
3750 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
3760 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
3800 niblock = CS%niblock
3810 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.
3850 u_in(:,:,:) = u(:,:,:)
3860 v_in(:,:,:) = v(:,:,:)
387
3880 call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC)
3890 if (niblock == 0) niblock = G%iec-G%isc+2
3900 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, &
3930 uhbt=uhbt, visc_rem_u=visc_rem_u, u_cor=u)
394
3950 call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC)
3960 if (niblock == 0) niblock = G%iec-G%isc+1
3970 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, &
4000 vhbt=vhbt, visc_rem_v=visc_rem_v, v_cor=v)
401
4020end subroutine continuity_adjust_vel
403
404
405!> Updates the thicknesses due to zonal thickness fluxes.
40613subroutine 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, u*h*dy [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
42313 call cpu_clock_begin(id_clock_update)
424
42513 h_min = 0.0 ; if (present(hmin)) h_min = hmin
426
42713 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
428
42913 if (present(hin)) then
43013 do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
43120025733 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
4350 do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
4360 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
44013 call cpu_clock_end(id_clock_update)
441
44213end subroutine continuity_zonal_convergence
443
444!> Updates the thicknesses due to meridional thickness fluxes.
4450subroutine 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, v*h*dx [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
46213 call cpu_clock_begin(id_clock_update)
463
46413 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
465
46613 h_min = 0.0 ; if (present(hmin)) h_min = hmin
467
46813 if (present(hin)) then
469 ! untested
4700 do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
4710 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
47413 do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
47519379893 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
47913 call cpu_clock_end(id_clock_update)
480
48113end subroutine continuity_merdional_convergence
482
483
484!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.
48513subroutine 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
50413 call cpu_clock_begin(id_clock_reconstruct)
505
50613 if (present(LB_in)) then
50713 LB = LB_in
508 else
5090 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
510 endif
51113 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
51213 nkblock = CS%nkblock
51313 if (nkblock == 0) nkblock = nz
514
51513 if (CS%upwind_1st) then
5160 do concurrent (k=1:nz, j=jsh:jeh, i=ish-1:ieh+1)
5170 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, &
52113 nkblock, 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC)
522 endif
523
52413 call cpu_clock_end(id_clock_reconstruct)
525
52613end subroutine zonal_edge_thickness
527
528
529!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.
53013subroutine 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
54913 call cpu_clock_begin(id_clock_reconstruct)
550
55113 if (present(LB_in)) then
55213 LB = LB_in
553 else
5540 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
555 endif
55613 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
55713 nkblock = CS%nkblock
55813 if (nkblock == 0) nkblock = nz
559
56013 if (CS%upwind_1st) then
561 ! untested
5620 do concurrent (k=1:nz, j=jsh-1:jeh+1, i=ish:ieh)
5630 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, &
56713 nkblock, 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC)
568 endif
569
57013 call cpu_clock_end(id_clock_reconstruct)
571
57213end subroutine meridional_edge_thickness
573
574
575!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.
57639subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_face_areaU, &
57728 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 = u*h*dy
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)) :: &
62226 duhdu, & ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
62326 visc_rem ! A 2-D copy of visc_rem_u or an array of 1's [nondim].
624 real, dimension(niblock,njblock) :: &
62526 du, & ! Corrective barotropic change in the velocity to give uhbt [L T-1 ~> m s-1].
62626 du_min_CFL, & ! Lower limit on du correction to avoid CFL violations [L T-1 ~> m s-1]
62726 du_max_CFL, & ! Upper limit on du correction to avoid CFL violations [L T-1 ~> m s-1]
62826 duhdu_tot_0, & ! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
62926 uh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
63026 visc_rem_max ! The column maximum of visc_rem [nondim].
63126 logical, dimension(niblock,njblock) :: do_I
632 real, dimension(niblock,njblock,SZK_(GV)) :: &
63326 uh_t ! A block-sized copy of uh [H L2 T-1 ~> m3 s-1 or kg s-1].
634 real, dimension(niblock,njblock) :: &
63526 uhbt_t ! A block-sized copy of uhbt [H L2 T-1 ~> m3 s-1 or kg s-1].
63626 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
64926 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
65613 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
66213 use_visc_rem = present(visc_rem_u)
663
66413 set_BT_cont = .false. ; if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
665
66613 local_specified_BC = .false. ; local_Flather_OBC = .false. ; local_open_BC = .false.
66713 if (associated(OBC)) then ; if (OBC%OBC_pe) then
6680 local_specified_BC = OBC%specified_u_BCs_exist_globally
6690 local_Flather_OBC = OBC%Flather_u_BCs_exist_globally
6700 local_open_BC = OBC%open_u_BCs_exist_globally
671 endif ; endif
672
67313 if (present(du_cor)) then
6740 do concurrent (j=G%jsd:G%jed, I=G%IsdB:G%IedB)
6750 du_cor(I,j) = 0.0
676 enddo
677 endif
678
67913 if (present(LB_in)) then
68013 LB = LB_in
681 else
6820 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
683 endif
68413 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
685
68613 CFL_dt = CS%CFL_limit_adjust / dt
68713 I_dt = 1.0 / dt
68813 if (CS%aggress_adjust) CFL_dt = I_dt
689
69013 if (.not.use_visc_rem) then
6911 do concurrent (k=1:nz, jj=1:njblock, ii=1:niblock)
6922977 visc_rem(ii,jj,k) = 1.0
693 enddo
694 endif
695
6968567 do j_start=jsh,jeh,njblock ; do i_start=ish-1,ieh,niblock
6977332 i_end = min(i_start+niblock-1,ieh)
6987332 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
703227903 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
7041100801 do_I(ii,jj) = .true.
705 enddo
706 ! Set uh and duhdu.
707 !$omp target teams num_teams(nteams)
708168636 do k=1,nz
709161304 if (use_visc_rem) then
710 !$omp loop collapse(2) private(ii,jj)
71118464688 do j=j_start,j_end ; do I=i_start,i_end
71217726544 ii=I-i_start+1 ; jj=j-j_start+1
71318315792 visc_rem(ii,jj,k) = visc_rem_u(I,j,k)
714 enddo ; enddo
715 endif
716 !$omp loop collapse(2) private(ii,jj)
71720003412 do j=j_start,j_end ; do i=i_start,i_end
71819203756 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,&
72219203756 por_face_areaU(I,j,k))
72319203756 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,&
726638352 OBC%segnum_u(i,j))
727 enddo ; enddo
728168636 if (local_specified_BC) then
729 !$omp loop collapse(2) private(l_seg)
7300 do j=j_start,j_end ; do i=i_start,i_end ; if (OBC%segnum_u(I,j) /= 0) then
7310 l_seg = abs(OBC%segnum_u(I,j))
7320 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
7387332 if (present(uhbt) .or. set_BT_cont) then
7396768 if (use_visc_rem .and. CS%use_visc_rem_max) then
740 !$omp target
741 !$omp loop collapse(2) private(ii,jj)
742839304 do j=j_start,j_end ; do i=i_start,i_end
743805752 ii=i-i_start+1 ; jj=j-j_start+1
744832536 visc_rem_max(ii,jj) = 0.0
745 enddo ; enddo
746155664 do k=1,nz
747 !$omp loop collapse(2) private(ii,jj)
74818471456 do j=j_start,j_end ; do i=i_start,i_end
74917726544 ii=i-i_start+1 ; jj=j-j_start+1
75018315792 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
7550 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
7560 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)
763839304 do j=j_start,j_end ; do I=i_start,i_end
764805752 ii=I-i_start+1 ; jj=j-j_start+1
765805752 I_vrm = 0.0
766805752 if (visc_rem_max(ii,jj) > 0.0) I_vrm = 1.0 / visc_rem_max(ii,jj)
767805752 if (CS%vol_CFL) then
7680 dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
7690 dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
770805752 else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
771805752 du_max_CFL(ii,jj) = 2.0* (CFL_dt * dx_W) * I_vrm
772805752 du_min_CFL(ii,jj) = -2.0 * (CFL_dt * dx_E) * I_vrm
773832536 uh_tot_0(ii,jj) = 0.0 ; duhdu_tot_0(ii,jj) = 0.0
774 enddo ; enddo
775155664 do k=1,nz
776 !$omp loop collapse(2) private(ii,jj)
77718471456 do j=j_start,j_end ; do I=i_start,i_end
77817726544 ii=I-i_start+1 ; jj=j-j_start+1
77917726544 duhdu_tot_0(ii,jj) = duhdu_tot_0(ii,jj) + duhdu(ii,jj,k)
78018315792 uh_tot_0(ii,jj) = uh_tot_0(ii,jj) + uh_t(ii,jj,k)
781 enddo
782 enddo ; enddo
783 !$omp end target
7846768 if (use_visc_rem) then
7856768 if (CS%aggress_adjust) then
786 !$omp target
7870 do k=1,nz
788 !$omp loop collapse(2) private(ii,jj,dx_W,dx_E,du_lim)
7890 do j=j_start,j_end ; do I=i_start,i_end
7900 ii=I-i_start+1 ; jj=j-j_start+1
7910 if (CS%vol_CFL) then
7920 dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
7930 dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
7940 else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
795
7960 du_lim = 0.499*((dx_W*I_dt - u(I,j,k)) + MIN(0.0,u(I-1,j,k)))
7970 if (du_max_CFL(ii,jj) * visc_rem(ii,jj,k) > du_lim) &
7980 du_max_CFL(ii,jj) = du_lim / visc_rem(ii,jj,k)
799
8000 du_lim = 0.499*((-dx_E*I_dt - u(I,j,k)) + MAX(0.0,u(I+1,j,k)))
8010 if (du_min_CFL(ii,jj) * visc_rem(ii,jj,k) < du_lim) &
8020 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
808155664 do k=1,nz
809 !$omp loop collapse(2) private(ii,jj,dx_E,dx_W)
81018471456 do j=j_start,j_end ; do I=i_start,i_end
81117726544 ii=I-i_start+1 ; jj=j-j_start+1
81217726544 if (CS%vol_CFL) then
8130 dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
8140 dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
81517726544 else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
816
81717726544 if (du_max_CFL(ii,jj) * visc_rem(ii,jj,k) > dx_W*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
818585021 du_max_CFL(ii,jj) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(ii,jj,k)
81917726544 if (du_min_CFL(ii,jj) * visc_rem(ii,jj,k) < -dx_E*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
8209532118 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
8260 if (CS%aggress_adjust) then
827 !$omp target
8280 do k=1,nz
829 !$omp loop collapse(2) private(ii,jj,dx_W,dx_E)
8300 do j=j_start,j_end ; do I=i_start,i_end
8310 ii=I-i_start+1 ; jj=j-j_start+1
8320 if (CS%vol_CFL) then
8330 dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
8340 dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
8350 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 * &
8380 ((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 * &
8400 ((-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
8460 do k=1,nz
847 !$omp loop collapse(2) private(ii,jj,dx_W,dx_E)
8480 do j=j_start,j_end ; do I=i_start,i_end
8490 ii=I-i_start+1 ; jj=j-j_start+1
8500 if (CS%vol_CFL) then
8510 dx_W = ratio_max(G%areaT(i,j), G%dy_Cu(I,j), 1000.0*G%dxT(i,j))
8520 dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
8530 else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
854
8550 du_max_CFL(ii,jj) = MIN(du_max_CFL(ii,jj), dx_W*CFL_dt - u(I,j,k))
8560 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
862210372 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
863805752 du_max_CFL(ii,jj) = max(du_max_CFL(ii,jj),0.0)
8641016124 du_min_CFL(ii,jj) = min(du_min_CFL(ii,jj),0.0)
865 enddo
866
8676768 any_simple_OBC = .false.
8686768 if (present(uhbt) .or. set_BT_cont) then
8696768 if (local_specified_BC .or. local_Flather_OBC) then
8700 do concurrent (j=j_start:j_end, I=i_start:i_end) DO_LOCALITY(reduce(.or.:any_simple_OBC))
8710 ii=I-i_start+1 ; jj=j-j_start+1
8720 l_seg = abs(OBC%segnum_u(I,j))
873
874 ! Avoid reconciling barotropic/baroclinic transports if transport is specified
8750 simple_OBC_pt(ii,jj) = .false.
8760 if (l_seg /= OBC_NONE) simple_OBC_pt(ii,jj) = OBC%segment(l_seg)%specified
8770 do_I(ii,jj) = .not.simple_OBC_pt(ii,jj)
8780 any_simple_OBC = any_simple_OBC .or. simple_OBC_pt(ii,jj)
879 enddo ; else
880210372 do concurrent (jj=1:j_end-j_start+1, II=1:i_end-i_start+1)
8811016124 do_I(II,jj) = .true.
882 enddo ; endif
883 endif
884
8856768 if (present(uhbt)) then
8864512 do concurrent (j=j_start:j_end, I=i_start:i_end)
887677416 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, &
8934512 uh_t, OBC=OBC)
894
8954512 if (present(u_cor)) then
8964512 do concurrent (k=1:nz, j=j_start:j_end, I=i_start:i_end)
89711817696 ii=I-i_start+1 ; jj=j-j_start+1
89811817696 u_cor(I,j,k) = u(I,j,k) + du(ii,jj) * visc_rem(ii,jj,k)
89912495112 if (any_simple_OBC) then ; if (simple_OBC_pt(ii,jj)) then
9000 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
9054512 if (present(du_cor)) then
9060 do concurrent (j=j_start:j_end, I=i_start:i_end)
9070 ii=I-i_start+1 ; jj=j-j_start+1
9080 du_cor(I,j) = du(ii,jj)
909 enddo
910 endif
911
912 endif
913
9146768 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, &
9184512 por_face_areaU, niblock, njblock)
9194512 if (any_simple_OBC) then
920 !$omp target
921 !$omp loop collapse(2) private(ii,jj)
9220 do j=j_start,j_end ; do I=i_start,i_end
9230 ii=i-i_start+1 ; jj=j-j_start+1
9240 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
9270 do k=1,nz
928 !$omp loop collapse(2) private(ii,jj,l_seg)
9290 do j=j_start,j_end ; do I=i_start,i_end
9300 ii=i-i_start+1 ; jj=j-j_start+1
9310 if (simple_OBC_pt(II,jj)) then
9320 l_seg = abs(OBC%segnum_u(I,j))
9330 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) + &
9360 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)
9410 do j=j_start,j_end ; do I=i_start,i_end
9420 ii=i-i_start+1 ; jj=j-j_start+1
9430 if (simple_OBC_pt(ii,jj)) then
9440 BT_cont%FA_u_W0(I,j) = FAuI(ii,jj) ; BT_cont%FA_u_E0(I,j) = FAuI(ii,jj)
9450 BT_cont%FA_u_WW(I,j) = FAuI(ii,jj) ; BT_cont%FA_u_EE(I,j) = FAuI(ii,jj)
9460 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
955611 do concurrent (k=1:nz, j=j_start:j_end, I=i_start:i_end)
95620304557 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
96513 if (local_open_BC .and. set_BT_cont) then
9660 do n = 1, OBC%number_of_segments
9670 if (OBC%segment(n)%open .and. OBC%segment(n)%is_E_or_W) then
9680 I = OBC%segment(n)%HI%IsdB
9690 if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
9700 do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed
9710 FA_u = 0.0
9720 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
9730 BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u
9740 BT_cont%FA_u_WW(I,j) = FA_u ; BT_cont%FA_u_EE(I,j) = FA_u
9750 BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0
976 enddo
977 else
9780 do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed
9790 FA_u = 0.0
9800 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
9810 BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u
9820 BT_cont%FA_u_WW(I,j) = FA_u ; BT_cont%FA_u_EE(I,j) = FA_u
9830 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
99013 if (set_BT_cont) then ; if (allocated(BT_cont%h_u)) then
9918 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, &
9934 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, &
9964 CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU, visc_rem_u)
997 endif
998 endif ; endif
999
100013 call cpu_clock_end(id_clock_correct)
1001
100233end subroutine zonal_mass_flux
1003
1004
1005!> Calculates the vertically integrated mass or volume fluxes through the zonal faces.
10060subroutine 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
10280 real :: uh(SZIB_(G),SZJ_(G),SZK_(GV)) ! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]
10290 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
10320 logical, dimension(SZJ_(G)) :: OBC_in_row
1033
10340 call cpu_clock_begin(id_clock_correct)
1035
10360 local_specified_BC = .false.
10370 if (associated(OBC)) then ; if (OBC%OBC_pe) then
10380 local_specified_BC = OBC%specified_v_BCs_exist_globally
1039 endif ; endif
1040
10410 if (present(LB_in)) then
10420 ish = LB_in%ish ; ieh = LB_in%ieh ; jsh = LB_in%jsh ; jeh = LB_in%jeh ; nz = GV%ke
1043 else
10440 ish = G%isc ; ieh = G%iec ; jsh = G%jsc ; jeh = G%jec ; nz = GV%ke
1045 endif
1046
10470 OBC_in_row(:) = .false.
1048
10490 uhbt(:,:) = 0.0
1050
1051 ! Determining whether there are any OBC points outside of the k-loop should be more efficient.
10520 if (local_specified_BC) then
10530 do j=jsh,jeh ; do I=ish-1,ieh ; if (OBC%segnum_u(I,j) /= 0) then
10540 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.
10590 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, &
10630 por_face_areaU(I,j,k))
10640 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, &
10660 por_face_areaU(I,j,k), G%dy_Cu(I,j), OBC, OBC%segnum_u(I,j))
1067 enddo
1068
10690 do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
10700 if (OBC_in_row(j) .and. OBC%segnum_u(I,j) /= 0) then
10710 l_seg = abs(OBC%segnum_u(I,j))
10720 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.
10770 do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh
10780 uhbt(I,j) = uhbt(I,j) + uh(I,j,k)
1079 enddo ; enddo ; enddo ! j-loop
1080
10810 call cpu_clock_end(id_clock_correct)
1082
10830end subroutine zonal_BT_mass_flux
1084
1085!> Evaluates the zonal mass or volume fluxes in an element.
1086!DIR$ ATTRIBUTES FORCEINLINE :: flux_elem
1087160752768elemental 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.
1124160752768 tmp = G_dy_Cu * por_face_area ! precalculate things
1125160752768 if (u > 0.0) then
112670037623 if (vol_CFL) then ; CFL = (u * dt) * (G_dy_Cu * G_IareaT)
112770037623 else ; CFL = u * dt * G_IdxT ; endif
112870037623 curv_3 = (h_L + h_R) - 2.0*h
112970037623 dh = h_L - h_R
1130 uh = tmp * u * &
113170037623 (h_R + CFL * (0.5*dh + curv_3*(CFL - 1.5)))
113270037623 h_marg = h_R + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
113390715145 elseif (u < 0.0) then
113451801243 if (vol_CFL) then ; CFL = (-u * dt) * (G_dy_Cu * G_IareaT_p1)
113551801243 else ; CFL = -u * dt * G_IdxT_p1 ; endif
113651801243 curv_3 = (h_L_p1 + h_R_p1) - 2.0*h_p1
113751801243 dh = h_R_p1-h_L_p1
1138 uh = tmp * u * &
113951801243 (h_L_p1 + CFL * (0.5*dh + curv_3*(CFL - 1.5)))
114051801243 h_marg = h_L_p1 + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
1141 else
114238913902 uh = 0.0
114338913902 h_marg = 0.5 * (h_L_p1 + h_R)
1144 endif
1145160752768 duhdu = tmp * h_marg * visc_rem
1146
1147160752768end subroutine flux_elem
1148
1149!DIR$ ATTRIBUTES FORCEINLINE :: flux_elem_OBC
11500elemental 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
11730 if (l_seg /= 0) then
11740 if (OBC%segment(abs(l_seg))%open) then
11750 if (l_seg > 0) then ! OBC_DIRECTION_E or OBC_DIRECTION_N
11760 uh = (G_dy_Cu * por_face_area) * u * h
11770 duhdu = (G_dy_Cu * por_face_area) * h * visc_rem
1178 else ! OBC_DIRECTION_W or OBC_DIRECTION_S
11790 uh = (G_dy_Cu * por_face_area) * u * h_p1
11800 duhdu = (G_dy_Cu* por_face_area) * h_p1 * visc_rem
1181 endif
1182 endif
1183 endif
1184
11850end 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.
11908subroutine zonal_flux_thickness(u, h, h_W, h_E, h_u, dt, G, GV, US, LB, vol_CFL, &
11918 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
12288 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
1229
12308 do concurrent (k=1:nz, j=jsh:jeh, I=ish-1:ieh)
123111817696 if (u(I,j,k) > 0.0) then
12324685842 if (vol_CFL) then ; CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j))
12334685842 else ; CFL = u(I,j,k) * dt * G%IdxT(i,j) ; endif
12344685842 curv_3 = (h_W(i,j,k) + h_E(i,j,k)) - 2.0*h(i,j,k)
12354685842 dh = h_W(i,j,k) - h_E(i,j,k)
12364685842 if (marginal) then
12374685842 h_u(I,j,k) = h_E(i,j,k) + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
1238 else
12390 h_u(I,j,k) = h_E(i,j,k) + CFL * (0.5*dh + curv_3*(CFL - 1.5))
1240 endif
12417131854 elseif (u(I,j,k) < 0.0) then
12422408476 if (vol_CFL) then ; CFL = (-u(I,j,k)*dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j))
12432408476 else ; CFL = -u(I,j,k) * dt * G%IdxT(i+1,j) ; endif
12442408476 curv_3 = (h_W(i+1,j,k) + h_E(i+1,j,k)) - 2.0*h(i+1,j,k)
12452408476 dh = h_E(i+1,j,k)-h_W(i+1,j,k)
12462408476 if (marginal) then
12472408476 h_u(I,j,k) = h_W(i+1,j,k) + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
1248 else
12490 h_u(I,j,k) = h_W(i+1,j,k) + CFL * (0.5*dh + 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.
12544723378 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
125924175456 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.
126311817696 h_u(I,j,k) = h_u(I,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k))
1264 else
12650 h_u(I,j,k) = h_u(I,j,k) * por_face_areaU(I,j,k)
1266 endif
1267 enddo
1268
12698 local_open_BC = .false.
12708 if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally
12718 if (local_open_BC) then
1272 ! untested
12730 do n = 1, OBC%number_of_segments
12740 if (OBC%segment(n)%open .and. OBC%segment(n)%is_E_or_W) then
12750 I = OBC%segment(n)%HI%IsdB
12760 if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
12770 if (present(visc_rem_u)) then
12780 do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
12790 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
12820 do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
12830 h_u(I,j,k) = h(i,j,k) * por_face_areaU(I,j,k)
1284 enddo
1285 endif
1286 else
12870 if (present(visc_rem_u)) then
12880 do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
12890 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
12920 do concurrent (k=1:nz, j = OBC%segment(n)%HI%jsd:OBC%segment(n)%HI%jed)
12930 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
13018end subroutine zonal_flux_thickness
1302
1303!> Returns the barotropic velocity adjustment that gives the
1304!! desired barotropic (layer-summed) transport.
13059024subroutine zonal_flux_adjust(u, h_in, h_W, h_E, uhbt, uh_tot_0, duhdu_tot_0, &
13069024 du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
13079024 i_start, i_end, j_start, j_end, do_I_in, por_face_areaU, niblock, njblock, &
13089024 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 = u*h*dy [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) :: &
135718048 uh_err, & ! Difference between uhbt and the summed uh [H L2 T-1 ~> m3 s-1 or kg s-1].
135818048 uh_err_best, & ! The smallest value of uh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
135918048 duhdu_tot,&! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
136018048 du_min, & ! Lower limit on du correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
136118048 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
136718048 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
13749024 nz = GV%ke
13759024 local_open_BC = .false.
13769024 if (present(OBC)) then ; if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally ; endif
1377
13789024 tol_vel = CS%tol_vel
1379
1380 !$omp target
1381
1382 !$omp loop collapse(2) private(ii,jj)
13831119072 do j=j_start,j_end ; do I=i_start,i_end
13841074336 ii=I-i_start+1 ; jj=j-j_start+1
13851074336 du(ii,jj) = 0.0 ; do_I(ii,jj) = do_I_in(ii,jj)
13861074336 du_max(ii,jj) = du_max_CFL(ii,jj) ; du_min(ii,jj) = du_min_CFL(ii,jj)
13871074336 uh_err(ii,jj) = uh_tot_0(ii,jj) - uhbt(ii,jj) ; duhdu_tot(ii,jj) = duhdu_tot_0(ii,jj)
13881110048 uh_err_best(ii,jj) = abs(uh_err(ii,jj))
1389 enddo ; enddo
1390
139121112 do itt=1,max_itts
13929024 select case (itt)
13939024 case (:1) ; tol_eta = 1e-6 * CS%tol_eta
13946585 case (2) ; tol_eta = 1e-4 * CS%tol_eta
13955352 case (3) ; tol_eta = 1e-2 * CS%tol_eta
139621112 case default ; tol_eta = CS%tol_eta
1397 end select
1398
1399 !$omp loop collapse(2) private(ii,jj)
14002687952 do j=j_start,j_end ; do I=i_start,i_end
14012582776 ii=I-i_start+1 ; jj=j-j_start+1
14022666840 if (uh_err(ii,jj) > 0.0) then ; du_max(ii,jj) = du(ii,jj)
14031878060 elseif (uh_err(ii,jj) < 0.0) then ; du_min(ii,jj) = du(ii,jj)
1404733705 else ; do_I(ii,jj) = .false. ; endif
1405 enddo ; enddo
1406#ifndef __NVCOMPILER_OPENMP_GPU
140721112 domore = .false.
1408#endif
1409 !$omp loop collapse(2) private(ii,jj,ddu,du_prev)
14102687952 do j=j_start,j_end ; do I=i_start,i_end
14112582776 ii=I-i_start+1 ; jj=j-j_start+1
14122666840 if (do_I(ii,jj)) then
14131741878 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.
14181074861 ddu = -uh_err(ii,jj) / duhdu_tot(ii,jj)
14191074861 du_prev = du(ii,jj)
14201074861 du(ii,jj) = du(ii,jj) + ddu
14211074861 if (abs(ddu) < 1.0e-15*abs(du(ii,jj))) then
14220 do_I(ii,jj) = .false. ! ddu is small enough to quit.
14231074861 elseif (ddu > 0.0) then
1424744249 if (du(ii,jj) >= du_max(ii,jj)) then
14250 du(ii,jj) = 0.5*(du_prev + du_max(ii,jj))
14260 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
1429330612 if (du(ii,jj) <= du_min(ii,jj)) then
14300 du(ii,jj) = 0.5*(du_prev + du_min(ii,jj))
14310 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
14351074861 if (do_I(ii,jj)) domore = .true.
1436#endif
1437 else
1438667017 do_I(ii,jj) = .false.
1439 endif
1440 endif
1441 enddo ; enddo
1442#ifndef __NVCOMPILER_OPENMP_GPU
144321112 if (.not.domore) exit
1444#endif
1445
1446 !$omp loop collapse(2) private(ii,jj)
14471568880 do j=j_start,j_end ; do I=i_start,i_end
14481508440 ii=I-i_start+1 ; jj=j-j_start+1
14491556792 uh_err(ii,jj) = -uhbt(ii,jj) ; duhdu_tot(ii,jj) = 0.0
1450 enddo ; enddo
1451278024 do k=1,nz
1452 !$omp loop collapse(2) private(ii,jj,duhdu,u_new)
145334527448 do j=j_start,j_end ; do i=i_start,i_end
145433185680 ii=i-i_start+1 ; jj=j-j_start+1
145534249424 if (do_I(ii,jj)) then
145623646942 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,&
146023646942 por_face_areaU(i,j,k))
146123646942 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,&
14640 OBC%segnum_u(i,j))
146523646942 uh_err(ii,jj) = uh_err(ii,jj) + uh_3d(ii,jj,k)
146623646942 duhdu_tot(ii,jj) = duhdu_tot(ii,jj) + duhdu
1467 endif
1468 enddo ; enddo
1469 enddo
1470 !$omp loop collapse(2) private(ii,jj)
14711568880 do j=j_start,j_end ; do I=i_start,i_end
14721508440 ii=I-i_start+1 ; jj=j-j_start+1
14731556792 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
14849024end 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.
14894512subroutine set_zonal_BT_cont(u, h_in, h_W, h_E, BT_cont, uh_tot_0, duhdu_tot_0, &
14904512 du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
14914512 visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
14924512 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) :: &
15339024 du0 ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
1534 real, dimension(niblock,njblock) :: &
15359024 zeros ! An array of full of 0 transports [H L2 T-1 ~> m3 s-1 or kg s-1]
1536 real, dimension(niblock,njblock) :: &
15379024 duL, duR, & ! The barotropic velocity increments that give the westerly
1538 ! (duL) and easterly (duR) test velocities [L T-1 ~> m s-1].
15399024 du_CFL, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
15409024 FAmt_L, FAmt_R, & ! The summed effective marginal face areas for the 3
15419024 FAmt_0, & ! test velocities [H L ~> m2 or kg m-1].
15429024 uhtot_L, & ! The summed transport with the westerly (uhtot_L) and
15439024 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].
15669024 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
15714512 nz = GV%ke ; Idt = 1.0 / dt
15724512 min_visc_rem = 0.1 ; CFL_min = 1e-6
1573
1574 ! Diagnose the zero-transport correction, du0.
1575140248 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
1576677416 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, &
15804512 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)
1587559536 do j=j_start,j_end ; do i=i_start,i_end
1588537168 ii=i-i_start+1 ; jj=j-j_start+1
1589537168 du_CFL(ii,jj) = (CFL_min * Idt) * G%dxCu(i,j)
1590537168 duR(ii,jj) = min(0.0,du0(ii,jj) - du_CFL(ii,jj))
1591537168 duL(ii,jj) = max(0.0,du0(ii,jj) + du_CFL(ii,jj))
1592537168 FAmt_L(ii,jj) = 0.0 ; FAmt_R(ii,jj) = 0.0 ; FAmt_0(ii,jj) = 0.0
1593555024 uhtot_L(ii,jj) = 0.0 ; uhtot_R(ii,jj) = 0.0
1594 enddo ; enddo
1595
1596103776 do k=1,nz
1597 !$omp loop collapse(2) private(ii,jj,visc_rem_lim)
159812314304 do j=j_start,j_end ; do i=i_start,i_end
159911817696 ii=i-i_start+1 ; jj=j-j_start+1
160012210528 if (do_I(ii,jj)) then
160111817696 visc_rem_lim = max(visc_rem(ii,jj,k), min_visc_rem*visc_rem_max(ii,jj))
160211817696 if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
16038107792 if (u(i,j,k) + duR(ii,jj)*visc_rem_lim > -du_CFL(ii,jj)*visc_rem(ii,jj,k)) &
1604319522 duR(ii,jj) = -(u(i,j,k) + du_CFL(ii,jj)*visc_rem(ii,jj,k)) / visc_rem_lim
16058107792 if (u(i,j,k) + duL(ii,jj)*visc_rem_lim < du_CFL(ii,jj)*visc_rem(ii,jj,k)) &
16061088028 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
1614103776 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)
161612314304 do j=j_start,j_end ; do i=i_start,i_end
161711817696 ii=i-i_start+1 ; jj=j-j_start+1
161812210528 if (do_I(ii,jj)) then
161911817696 u_L = u(i,j,k) + duL(ii,jj) * visc_rem(ii,jj,k)
162011817696 u_R = u(i,j,k) + duR(ii,jj) * visc_rem(ii,jj,k)
162111817696 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,&
162511817696 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,&
162911817696 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,&
163311817696 por_face_areaU(i,j,k))
163411817696 FAmt_0(ii,jj) = FAmt_0(ii,jj) + duhdu_0
163511817696 FAmt_L(ii,jj) = FAmt_L(ii,jj) + duhdu_L
163611817696 FAmt_R(ii,jj) = FAmt_R(ii,jj) + duhdu_R
163711817696 uhtot_L(ii,jj) = uhtot_L(ii,jj) + uh_L
163811817696 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)
1643559536 do j=j_start,j_end ; do i=i_start,i_end
1644537168 ii=i-i_start+1 ; jj=j-j_start+1
1645555024 if (do_I(ii,jj)) then
1646537168 FA_0 = FAmt_0(ii,jj) ; FA_avg = FAmt_0(ii,jj)
1647537168 if ((duL(ii,jj) - du0(ii,jj)) /= 0.0) &
1648537168 FA_avg = uhtot_L(ii,jj) / (duL(ii,jj) - du0(ii,jj))
1649537168 if (FA_avg > max(FA_0, FAmt_L(ii,jj))) then ; FA_avg = max(FA_0, FAmt_L(ii,jj))
1650531655 elseif (FA_avg < min(FA_0, FAmt_L(ii,jj))) then ; FA_0 = FA_avg ; endif
1651
1652537168 BT_cont%FA_u_W0(i,j) = FA_0 ; BT_cont%FA_u_WW(i,j) = FAmt_L(ii,jj)
1653537168 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))) * &
1655368516 ((FAmt_L(ii,jj) - FA_avg) / (FAmt_L(ii,jj) - FA_0))
1656 endif
1657
1658537168 FA_0 = FAmt_0(ii,jj) ; FA_avg = FAmt_0(ii,jj)
1659537168 if ((duR(ii,jj) - du0(ii,jj)) /= 0.0) &
1660537168 FA_avg = uhtot_R(ii,jj) / (duR(ii,jj) - du0(ii,jj))
1661537168 if (FA_avg > max(FA_0, FAmt_R(ii,jj))) then ; FA_avg = max(FA_0, FAmt_R(ii,jj))
1662537093 elseif (FA_avg < min(FA_0, FAmt_R(ii,jj))) then ; FA_0 = FA_avg ; endif
1663
1664537168 BT_cont%FA_u_E0(i,j) = FA_0 ; BT_cont%FA_u_EE(i,j) = FAmt_R(ii,jj)
1665537168 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))) * &
1667368529 ((FAmt_R(ii,jj) - FA_avg) / (FAmt_R(ii,jj) - FA_0))
1668 endif
1669 else
16700 BT_cont%FA_u_W0(i,j) = 0.0 ; BT_cont%FA_u_WW(i,j) = 0.0
16710 BT_cont%FA_u_E0(i,j) = 0.0 ; BT_cont%FA_u_EE(i,j) = 0.0
16720 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
16804512end subroutine set_zonal_BT_cont
1681
1682!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.
168339subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, por_face_areaV, &
16848 niblock, njblock, LB_in, vhbt, &
168520 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 = v*h*dx [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) :: &
172826 vhbt_t
1729 real, dimension(niblock,njblock,SZK_(GV)) :: &
173026 vh_t
1731 real, dimension(niblock,njblock,SZK_(GV)) :: &
173226 dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
1733 real, dimension(niblock,njblock) :: &
173426 dv, & ! Corrective barotropic change in the velocity to give vhbt [L T-1 ~> m s-1].
173526 dv_min_CFL, & ! Lower limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]
173626 dv_max_CFL, & ! Upper limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]
173726 dvhdv_tot_0, & ! Summed partial derivative of vh with v [H L ~> m2 or kg m-1].
173826 vh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
173926 visc_rem_max ! The column maximum of visc_rem [nondim]
174026 logical, dimension(niblock,njblock) :: do_I
174126 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)) :: &
174526 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
175826 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
176313 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
176913 use_visc_rem = present(visc_rem_v)
1770
177113 set_BT_cont = .false. ; if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
1772
177313 local_specified_BC = .false. ; local_Flather_OBC = .false. ; local_open_BC = .false.
177413 if (associated(OBC)) then ; if (OBC%OBC_pe) then
17750 local_specified_BC = OBC%specified_v_BCs_exist_globally
17760 local_Flather_OBC = OBC%Flather_v_BCs_exist_globally
17770 local_open_BC = OBC%open_v_BCs_exist_globally
1778 endif ; endif
1779
178013 if (present(dv_cor)) then
17810 do concurrent (jj=G%jsdB:G%jedB, ii=G%isd:G%ied)
17820 dv_cor(ii,jj) = 0.0
1783 enddo
1784 endif
1785
178613 if (present(LB_in)) then
178713 LB = LB_in
1788 else
17890 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
1790 endif
179113 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
1792
179313 CFL_dt = CS%CFL_limit_adjust / dt
179413 I_dt = 1.0 / dt
179513 if (CS%aggress_adjust) CFL_dt = I_dt
1796
179713 if (.not.use_visc_rem) then
17981 do concurrent (k=1:nz, jj=1:njblock, ii=1:niblock)
17992977 visc_rem(ii,jj,k) = 1.0
1800 enddo
1801 endif
1802
18038385 do j_start = jsh-1, jeh, njblock ; do i_start = ish, ieh, niblock
18047176 j_end = min(j_start + njblock-1, jeh)
18057176 i_end = min(i_start + niblock-1, ieh)
1806 !$ nteams = ceiling(real((j_end-j_start+1)*(i_end-i_start+1))/128.)
1807222456 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
18081069536 do_I(ii,JJ) = .true.
1809 enddo
1810 ! This sets vh and dvhdv.
1811 !$omp target teams num_teams(nteams)
1812165048 do k=1,nz
1813157872 if (use_visc_rem) then
1814 !$omp loop collapse(2) private(ii,jj)
181517921376 do J=j_start,j_end ; do i=i_start,i_end
181617202240 ii=i-i_start+1 ; jj=j-j_start+1
181717775648 visc_rem(ii,JJ,k) = visc_rem_v(i,J,k)
1818 enddo ; enddo
1819 endif
1820 !$omp loop collapse(2) private(ii,JJ)
182119414824 do J=j_start,j_end ; do i=i_start,i_end
182218635760 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,&
182618635760 por_face_areaV(i,J,k))
182718635760 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,&
1830621192 OBC%segnum_v(i,J))
1831 enddo ; enddo
1832165048 if (local_specified_BC) then
1833 !$omp loop collapse(2) private(l_seg)
18340 do J=J_start,J_end ; do i=i_start,i_end ; if (OBC%segnum_v(i,J) /= 0) then
18350 l_seg = abs(OBC%segnum_v(i,J))
18360 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
18427176 if (present(vhbt) .or. set_BT_cont) then
18436624 if (use_visc_rem .and. CS%use_visc_rem_max) then
1844 !$omp target
1845 !$omp loop collapse(2) private(ii,jj)
1846814608 do j=j_start,j_end ; do i=i_start,i_end
1847781920 ii=i-i_start+1 ; jj=j-j_start+1
1848807984 visc_rem_max(ii,JJ) = 0.0
1849 enddo ; enddo
1850152352 do k=1,nz
1851 !$omp loop collapse(2) private(ii,jj)
185217928000 do j=j_start,j_end ; do i=i_start,i_end
185317202240 ii=i-i_start+1 ; jj=j-j_start+1
185417775648 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
18590 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
18600 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)
1867814608 do J=J_start,J_end ; do i=i_start,i_end
1868781920 ii=i-i_start+1 ; jj=j-j_start+1
1869781920 I_vrm = 0.0
1870781920 if (visc_rem_max(ii,JJ) > 0.0) I_vrm = 1.0 / visc_rem_max(ii,JJ)
1871781920 if (CS%vol_CFL) then
18720 dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(i,J), 1000.0*G%dyT(i,j))
18730 dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(i,J), 1000.0*G%dyT(i,j+1))
1874781920 else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
1875781920 dv_max_CFL(ii,JJ) = 2.0 * (CFL_dt * dy_S) * I_vrm
1876781920 dv_min_CFL(ii,JJ) = -2.0 * (CFL_dt * dy_N) * I_vrm
1877807984 vh_tot_0(ii,JJ) = 0.0 ; dvhdv_tot_0(ii,JJ) = 0.0
1878 enddo ; enddo
1879152352 do k=1,nz
1880 !$omp loop collapse(2) private(ii,jj)
188117928000 do J=J_start,J_end ; do i=i_start,i_end
188217202240 ii=i-i_start+1 ; jj=j-j_start+1
188317202240 dvhdv_tot_0(ii,JJ) = dvhdv_tot_0(ii,JJ) + dvhdv(ii,jj,k)
188417775648 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
18896624 if (use_visc_rem) then
18906624 if (CS%aggress_adjust) then
1891 !$omp target
18920 do k=1,nz
1893 !$omp loop collapse(2) private(ii,jj,dy_S,dy_N,dv_lim)
18940 do J=J_start,J_end ; do i=i_start,i_end
18950 ii=i-i_start+1 ; jj=j-j_start+1
18960 if (CS%vol_CFL) then
18970 dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
18980 dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
18990 else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
19000 dv_lim = 0.499*((dy_S*I_dt - v(i,J,k)) + MIN(0.0,v(i,J-1,k)))
19010 if (dv_max_CFL(ii,JJ) * visc_rem(ii,JJ,k) > dv_lim) &
19020 dv_max_CFL(ii,JJ) = dv_lim / visc_rem(ii,JJ,k)
1903
19040 dv_lim = 0.499*((-dy_N*CFL_dt - v(i,J,k)) + MAX(0.0,v(i,J+1,k)))
19050 if (dv_min_CFL(ii,JJ) * visc_rem(ii,JJ,k) < dv_lim) &
19060 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
1912152352 do k=1,nz
1913 !$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
191417928000 do J=J_start,J_end ; do i=i_start,i_end
191517202240 ii=i-i_start+1 ; jj=j-j_start+1
191617202240 if (CS%vol_CFL) then
19170 dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
19180 dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
191917202240 else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
192017202240 if (dv_max_CFL(ii,JJ) * visc_rem(ii,JJ,k) > &
1921 dy_S*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
19226122188 dv_max_CFL(ii,JJ) = (dy_S*CFL_dt - v(i,J,k)) / visc_rem(ii,JJ,k)
192317202240 if (dv_min_CFL(ii,JJ) * visc_rem(ii,JJ,k) < &
1924 -dy_N*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
19251526878 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
19310 if (CS%aggress_adjust) then
1932 !$omp target
19330 do k=1,nz
1934 !$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
19350 do J=J_start,J_end ; do i=i_start,i_end
19360 ii=i-i_start+1 ; jj=j-j_start+1
19370 if (CS%vol_CFL) then
19380 dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
19390 dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
19400 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 * &
19420 ((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 * &
19440 ((-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
19500 do k=1,nz
1951 !$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
19520 do J=J_start,J_end ; do i=i_start,i_end
19530 ii=i-i_start+1 ; jj=j-j_start+1
19540 if (CS%vol_CFL) then
19550 dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
19560 dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
19570 else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
19580 dv_max_CFL(ii,JJ) = min(dv_max_CFL(ii,JJ), dy_S*CFL_dt - v(i,J,k))
19590 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
1965205344 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
1966781920 dv_max_CFL(ii,JJ) = max(dv_max_CFL(ii,JJ),0.0)
1967987264 dv_min_CFL(ii,JJ) = min(dv_min_CFL(ii,JJ),0.0)
1968 enddo
1969
19706624 any_simple_OBC = .false.
19716624 if (present(vhbt) .or. set_BT_cont) then
19726624 if (local_specified_BC .or. local_Flather_OBC) then
19730 do concurrent (J=J_start:J_end, i=i_start:i_end) DO_LOCALITY(reduce(.or.:any_simple_OBC))
19740 ii=i-i_start+1 ; jj=J-j_start+1
19750 l_seg = abs(OBC%segnum_v(i,J))
1976
1977 ! Avoid reconciling barotropic/baroclinic transports if transport is specified
19780 simple_OBC_pt(ii,jj) = .false.
19790 if (l_seg /= 0) simple_OBC_pt(ii,jj) = OBC%segment(l_seg)%specified
19800 do_I(ii,jj) = .not.simple_OBC_pt(ii,jj)
19810 any_simple_OBC = any_simple_OBC .or. simple_OBC_pt(ii,jj)
1982 enddo ; else
1983205344 do concurrent (JJ=1:J_end-J_start+1, ii=1:i_end-i_start+1)
1984987264 do_I(ii,JJ) = .true.
1985 enddo ; endif
1986 endif
1987
19886624 if (present(vhbt)) then
19894416 do concurrent (j=j_start:j_end, i=i_start:i_end)
1990521280 ii=i-i_start+1 ; jj=j-j_start+1
1991658176 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, &
19974416 niblock, njblock, vh_t, OBC)
1998
19994416 if (present(v_cor)) then
20004416 do concurrent (k=1:nz, J=J_start:J_end, i=i_start:i_end)
200111468160 ii=i-i_start+1 ; jj=j-j_start+1
200211468160 v_cor(i,J,k) = v(i,J,k) + dv(ii,JJ) * visc_rem(ii,JJ,k)
200312126336 if (any_simple_OBC) then ; if (simple_OBC_pt(ii,jj)) then
20040 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
20094416 if (present(dv_cor)) then
20100 do concurrent (J=J_start:J_end, i=i_start:i_end)
20110 ii=i-i_start+1 ; jj=j-j_start+1
20120 dv_cor(i,J) = dv(ii,JJ)
2013 enddo
2014 endif
2015 endif
2016
20176624 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, &
20214416 por_face_areaV, niblock, njblock)
20224416 if (any_simple_OBC) then
2023 !$omp target
2024 !$omp loop collapse(2) private(ii,jj)
20250 do J=J_start,J_end ; do i=i_start,i_end
20260 ii=i-i_start+1 ; jj=j-j_start+1
20270 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
20300 do k=1,nz
2031 !$omp loop collapse(2) private(ii,jj,l_seg)
20320 do J=J_start,J_end ; do i=i_start,i_end
20330 ii=i-i_start+1 ; jj=j-j_start+1
20340 if (simple_OBC_pt(ii,JJ)) then
20350 l_seg = abs(OBC%segnum_v(i,J))
20360 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) + &
20390 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)
20440 do J=J_start,J_end ; do i=i_start,i_end
20450 ii=i-i_start+1 ; jj=j-j_start+1
20460 if (simple_OBC_pt(ii,jj)) then
20470 BT_cont%FA_v_S0(i,J) = FAvi(ii,jj) ; BT_cont%FA_v_N0(i,J) = FAvi(ii,jj)
20480 BT_cont%FA_v_SS(i,J) = FAvi(ii,jj) ; BT_cont%FA_v_NN(i,J) = FAvi(ii,jj)
20490 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
2058598 do concurrent (k=1:nz, j=j_start:j_end, i=i_start:i_end)
205918635760 jj=j-j_start+1 ; ii=i-i_start+1
206019705296 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
206913 if (local_open_BC .and. set_BT_cont) then
20700 do n = 1, OBC%number_of_segments
20710 if (OBC%segment(n)%open .and. OBC%segment(n)%is_N_or_S) then
20720 J = OBC%segment(n)%HI%JsdB
20730 if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
20740 do i = OBC%segment(n)%HI%Isd, OBC%segment(n)%HI%Ied
20750 FA_v = 0.0
20760 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
20770 BT_cont%FA_v_S0(i,J) = FA_v ; BT_cont%FA_v_N0(i,J) = FA_v
20780 BT_cont%FA_v_SS(i,J) = FA_v ; BT_cont%FA_v_NN(i,J) = FA_v
20790 BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0
2080 enddo
2081 else
20820 do i = OBC%segment(n)%HI%Isd, OBC%segment(n)%HI%Ied
20830 FA_v = 0.0
20840 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
20850 BT_cont%FA_v_S0(i,J) = FA_v ; BT_cont%FA_v_N0(i,J) = FA_v
20860 BT_cont%FA_v_SS(i,J) = FA_v ; BT_cont%FA_v_NN(i,J) = FA_v
20870 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
209413 if (set_BT_cont) then ; if (allocated(BT_cont%h_v)) then
20958 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, &
20974 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, &
21004 CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaV, visc_rem_v)
2101 endif
2102 endif ; endif
2103
210413 call cpu_clock_end(id_clock_correct)
2105
210633end subroutine meridional_mass_flux
2107
2108
2109!> Calculates the vertically integrated mass or volume fluxes through the meridional faces.
21100subroutine 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
21320 real :: vh(SZI_(G),SZJB_(G),SZK_(GV)) ! Volume flux through meridional faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]
21330 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
21350 logical :: local_specified_BC, OBC_in_row(SZJB_(G))
2136
21370 call cpu_clock_begin(id_clock_correct)
2138
21390 local_specified_BC = .false.
21400 if (associated(OBC)) then ; if (OBC%OBC_pe) then
21410 local_specified_BC = OBC%specified_v_BCs_exist_globally
2142 endif ; endif
2143
21440 if (present(LB_in)) then
21450 ish = LB_in%ish ; ieh = LB_in%ieh ; jsh = LB_in%jsh ; jeh = LB_in%jeh ; nz = GV%ke
2146 else
21470 ish = G%isc ; ieh = G%iec ; jsh = G%jsc ; jeh = G%jec ; nz = GV%ke
2148 endif
2149
21500 vhbt(:,:) = 0.0
2151
2152 ! Determining whether there are any OBC points outside of the k-loop should be more efficient.
21530 OBC_in_row(:) = .false.
21540 if (local_specified_BC) then
21550 do j=jsh-1,jeh ; do i=ish,ieh ; if (OBC%segnum_v(i,J) /= 0) then
21560 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.
21610 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, &
21650 CS%vol_CFL, por_face_areaV(i,J,k))
21660 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, &
21680 por_face_areaV(i,J,k), G%dx_Cv(i,J), OBC, OBC%segnum_v(i,J))
2169 enddo
2170
21710 do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
21720 if (OBC_in_row(j) .and. OBC%segnum_v(i,J) /= 0) then
21730 l_seg = abs(OBC%segnum_v(i,J))
21740 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.
21800 do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh
21810 vhbt(i,J) = vhbt(i,J) + vh(i,J,k)
2182 enddo ; enddo ; enddo
2183
21840 call cpu_clock_end(id_clock_correct)
2185
21860end 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.
21918subroutine meridional_flux_thickness(v, h, h_S, h_N, h_v, dt, G, GV, US, LB, vol_CFL, &
21928 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
22318 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
2232
22338 do concurrent (k=1:nz, J=jsh-1:jeh, i=ish:ieh)
223411468160 if (v(i,J,k) > 0.0) then
22356350064 if (vol_CFL) then ; CFL = (v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j))
22366350064 else ; CFL = v(i,J,k) * dt * G%IdyT(i,j) ; endif
22376350064 curv_3 = (h_S(i,j,k) + h_N(i,j,k)) - 2.0*h(i,j,k)
22386350064 dh = h_S(i,J,k) - h_N(i,J,k)
22396350064 if (marginal) then
22406350064 h_v(i,J,k) = h_N(i,j,k) + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
2241 else
22420 h_v(i,J,k) = h_N(i,j,k) + CFL * (0.5*dh + curv_3*(CFL - 1.5))
2243 endif
22445118096 elseif (v(i,J,k) < 0.0) then
2245716996 if (vol_CFL) then ; CFL = (-v(i,J,k)*dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1))
2246716996 else ; CFL = -v(i,J,k) * dt * G%IdyT(i,j+1) ; endif
2247716996 curv_3 = (h_S(i,j+1,k) + h_N(i,j+1,k)) - 2.0*h(i,j+1,k)
2248716996 dh = h_N(i,j+1,k)-h_S(i,j+1,k)
2249716996 if (marginal) then
2250716996 h_v(i,J,k) = h_S(i,j+1,k) + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
2251 else
22520 h_v(i,J,k) = h_S(i,j+1,k) + CFL * (0.5*dh + 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.
22574401100 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
226223460488 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.
226611468160 h_v(i,J,k) = h_v(i,J,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k))
2267 else
22680 h_v(i,J,k) = h_v(i,J,k) * por_face_areaV(i,J,k)
2269 endif
2270 enddo
2271
22728 local_open_BC = .false.
22738 if (associated(OBC)) local_open_BC = OBC%open_v_BCs_exist_globally
2274 ! untested - will need to be refactored to be performant on GPUs
22758 if (local_open_BC) then
22760 do n = 1, OBC%number_of_segments
22770 if (OBC%segment(n)%open .and. OBC%segment(n)%is_N_or_S) then
22780 J = OBC%segment(n)%HI%JsdB
22790 if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
22800 if (present(visc_rem_v)) then
22810 do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
22820 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
22850 do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
22860 h_v(i,J,k) = h(i,J,k) * por_face_areaV(i,J,k)
2287 enddo
2288 endif
2289 else
22900 if (present(visc_rem_v)) then
22910 do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
22920 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
22950 do concurrent (k=1:nz, i = OBC%segment(n)%HI%isd:OBC%segment(n)%HI%ied)
22960 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
23048end subroutine meridional_flux_thickness
2305
2306
2307!> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.
23088832subroutine meridional_flux_adjust(v, h_in, h_S, h_N, vhbt, vh_tot_0, dvhdv_tot_0, &
23098832 dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
231017664 i_start, i_end, j_start, j_end, do_I_in, por_face_areaV, &
23118832 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 = v*h*dx [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) :: &
236417664 vh_err, & ! Difference between vhbt and the summed vh [H L2 T-1 ~> m3 s-1 or kg s-1].
236517664 vh_err_best, & ! The smallest value of vh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
236617664 dvhdv_tot,&! Summed partial derivative of vh with u [H L ~> m2 or kg m-1].
236717664 dv_min, & ! Lower limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
236817664 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
237617664 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
23828832 local_open_BC = .false.
23838832 if (present(OBC)) then ; if (associated(OBC)) then
23840 local_open_BC = OBC%open_v_BCs_exist_globally
2385 endif ; endif
2386
23878832 nz = GV%ke
23888832 tol_vel = CS%tol_vel
2389
2390 !$omp target
2391
2392 !$omp loop collapse(2) private(ii,jj)
23931086144 do J=J_start,J_end ; do i=i_start,i_end
23941042560 ii=i-i_start+1 ; jj=j-j_start+1
23951042560 dv(ii,JJ) = 0.0 ; do_I(ii,JJ) = do_I_in(ii,JJ)
23961042560 dv_max(ii,JJ) = dv_max_CFL(ii,JJ) ; dv_min(ii,JJ) = dv_min_CFL(ii,JJ)
23971042560 vh_err(ii,JJ) = vh_tot_0(ii,JJ) - vhbt(ii,JJ) ; dvhdv_tot(ii,JJ) = dvhdv_tot_0(ii,JJ)
23981077312 vh_err_best(ii,JJ) = abs(vh_err(ii,JJ))
2399 enddo ; enddo
2400
240121808 do itt=1,max_itts
24028832 select case (itt)
24038832 case (:1) ; tol_eta = 1e-6 * CS%tol_eta
24046480 case (2) ; tol_eta = 1e-4 * CS%tol_eta
24056370 case (3) ; tol_eta = 1e-2 * CS%tol_eta
240621808 case default ; tol_eta = CS%tol_eta
2407 end select
2408
2409 !$omp loop collapse(2) private(ii,jj)
24102768464 do j=j_start,j_end ; do i=i_start,i_end
24112660000 ii=i-i_start+1 ; jj=j-j_start+1
24122746656 if (vh_err(ii,JJ) > 0.0) then ; dv_max(ii,JJ) = dv(ii,JJ)
24131239078 elseif (vh_err(ii,JJ) < 0.0) then ; dv_min(ii,JJ) = dv(ii,JJ)
2414642847 else ; do_I(ii,JJ) = .false. ; endif
2415 enddo ; enddo
2416#ifndef __NVCOMPILER_OPENMP_GPU
241721808 domore = .false.
2418#endif
2419 !$omp loop collapse(2) private(ii,jj,ddv,dv_prev)
24202768464 do J=J_start,J_end ; do i=i_start,i_end
24212746656 ii=i-i_start+1 ; jj=j-j_start+1 ; if (do_I(ii,JJ)) then
24221986113 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.
24271336761 ddv = -vh_err(ii,JJ) / dvhdv_tot(ii,JJ)
24281336761 dv_prev = dv(ii,JJ)
24291336761 dv(ii,JJ) = dv(ii,JJ) + ddv
24301336761 if (abs(ddv) < 1.0e-15*abs(dv(ii,JJ))) then
24310 do_I(ii,JJ) = .false. ! ddv is small enough to quit.
24321336761 elseif (ddv > 0.0) then
2433276693 if (dv(ii,JJ) >= dv_max(ii,JJ)) then
24340 dv(ii,JJ) = 0.5*(dv_prev + dv_max(ii,JJ))
24350 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
24381060068 if (dv(ii,JJ) <= dv_min(ii,JJ)) then
24390 dv(ii,JJ) = 0.5*(dv_prev + dv_min(ii,JJ))
24400 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
24441336761 if (do_I(ii,JJ)) domore = .true.
2445#endif
2446 else
2447649352 do_I(ii,JJ) = .false.
2448 endif
2449 endif ; enddo ; enddo
2450#ifndef __NVCOMPILER_OPENMP_GPU
245121808 if (.not.domore) exit
2452#endif
2453 !$omp loop collapse(2) private(ii,jj)
24541682320 do J=J_start,J_end ; do i=i_start,i_end
24551617440 ii=i-i_start+1 ; jj=J-j_start+1
24561669344 vh_err(ii,JJ) = -vhbt(ii,JJ) ; dvhdv_tot(ii,JJ) = 0.0
2457 enddo ; enddo
2458298448 do k=1,nz
2459 !$omp loop collapse(2) private(ii,jj,dvhdv,v_new)
246037024016 do J=J_start,J_end ; do i=i_start,i_end
246135583680 ii=i-i_start+1 ; jj=j-j_start+1
246236725568 if (do_I(ii,JJ)) then
246329408742 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,&
246729408742 por_face_areaV(i,J,k))
246829408742 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,&
24710 OBC%segnum_v(i,J))
247229408742 vh_err(ii,JJ) = vh_err(ii,JJ) + vh_3d(ii,JJ,k)
247329408742 dvhdv_tot(ii,JJ) = dvhdv_tot(ii,JJ) + dvhdv
2474 endif
2475 enddo ; enddo
2476 enddo
2477 !$omp loop collapse(2) private(ii,jj)
24781682320 do J=J_start,J_end ; do i=i_start,i_end
24791617440 ii=i-i_start+1 ; jj=J-j_start+1
24801669344 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
24918832end 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.
24964416subroutine set_merid_BT_cont(v, h_in, h_S, h_N, BT_cont, vh_tot_0, dvhdv_tot_0, &
24974416 dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
24984416 visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
24994416 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) :: &
25438832 dv0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
25448832 dvL, dvR, & ! The barotropic velocity increments that give the southerly
2545 ! (dvL) and northerly (dvR) test velocities [L T-1 ~> m s-1].
25468832 dv_CFL, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
25478832 FAmt_L, FAmt_R, & ! The summed effective marginal face areas for the 3
25488832 FAmt_0, & ! test velocities [H L ~> m2 or kg m-1].
25498832 vhtot_L, & ! The summed transport with the southerly (vhtot_L) and
25508832 vhtot_R ! and northerly (vhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
2551 real, dimension(niblock,njblock) :: &
25528832 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
25758832 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
25804416 nz = GV%ke ; Idt = 1.0 / dt
25814416 min_visc_rem = 0.1 ; CFL_min = 1e-6
2582
2583 ! Diagnose the zero-transport correction, dv0.
2584136896 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
2585658176 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, &
25904416 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)
2597543072 do J=J_start,J_end ; do i=i_start,i_end
2598521280 ii=i-i_start+1 ; jj=j-j_start+1
2599521280 dv_CFL(ii,JJ) = (CFL_min * Idt) * G%dyCv(i,J)
2600521280 dvR(ii,JJ) = min(0.0,dv0(ii,JJ) - dv_CFL(ii,JJ))
2601521280 dvL(ii,JJ) = max(0.0,dv0(ii,JJ) + dv_CFL(ii,JJ))
2602521280 FAmt_L(ii,JJ) = 0.0 ; FAmt_R(ii,JJ) = 0.0 ; FAmt_0(ii,JJ) = 0.0
2603538656 vhtot_L(ii,JJ) = 0.0 ; vhtot_R(ii,JJ) = 0.0
2604 enddo ; enddo
2605
2606101568 do k=1,nz
2607 !$omp loop collapse(2) private(ii,jj,visc_rem_lim)
260811952000 do J=J_start,J_end ; do i=i_start,i_end
260911468160 ii=i-i_start+1 ; jj=j-j_start+1
261011850432 if (do_I(ii,jj)) then
261111468160 visc_rem_lim = max(visc_rem(ii,JJ,k), min_visc_rem*visc_rem_max(ii,JJ))
261211468160 if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
26138076640 if (v(i,J,k) + dvR(ii,JJ)*visc_rem_lim > -dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) &
2614727090 dvR(ii,JJ) = -(v(i,J,k) + dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) / visc_rem_lim
26158076640 if (v(i,J,k) + dvL(ii,JJ)*visc_rem_lim < dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) &
261677340 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
2621101568 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)
262311952000 do J=J_start,J_end ; do i=i_start,i_end
262411468160 ii=i-i_start+1 ; jj=j-j_start+1
262511850432 if (do_I(ii,jj)) then
262611468160 v_L = v(i,J,k) + dvL(ii,JJ) * visc_rem(ii,JJ,k)
262711468160 v_R = v(i,J,k) + dvR(ii,JJ) * visc_rem(ii,JJ,k)
262811468160 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),&
263111468160 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),&
263411468160 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),&
263711468160 G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,por_face_areaV(i,J,k))
263811468160 FAmt_0(ii,JJ) = FAmt_0(ii,JJ) + dvhdv_0
263911468160 FAmt_L(ii,JJ) = FAmt_L(ii,JJ) + dvhdv_L
264011468160 FAmt_R(ii,JJ) = FAmt_R(ii,JJ) + dvhdv_R
264111468160 vhtot_L(ii,JJ) = vhtot_L(ii,JJ) + vh_L
264211468160 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)
2647543072 do J=J_start,J_end ; do i=i_start,i_end
2648521280 ii=i-i_start+1 ; jj=j-j_start+1
2649538656 if (do_I(ii,JJ)) then
2650521280 FA_0 = FAmt_0(ii,JJ) ; FA_avg = FAmt_0(ii,JJ)
2651521280 if ((dvL(ii,JJ) - dv0(ii,JJ)) /= 0.0) &
2652521280 FA_avg = vhtot_L(ii,JJ) / (dvL(ii,JJ) - dv0(ii,JJ))
2653521280 if (FA_avg > max(FA_0, FAmt_L(ii,JJ))) then ; FA_avg = max(FA_0, FAmt_L(ii,JJ))
2654520702 elseif (FA_avg < min(FA_0, FAmt_L(ii,JJ))) then ; FA_0 = FA_avg ; endif
2655521280 BT_cont%FA_v_S0(i,J) = FA_0 ; BT_cont%FA_v_SS(i,J) = FAmt_L(ii,JJ)
2656521280 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))) * &
2658367103 ((FAmt_L(ii,JJ) - FA_avg) / (FAmt_L(ii,JJ) - FA_0))
2659 endif
2660
2661521280 FA_0 = FAmt_0(ii,JJ) ; FA_avg = FAmt_0(ii,JJ)
2662521280 if ((dvR(ii,JJ) - dv0(ii,JJ)) /= 0.0) &
2663521280 FA_avg = vhtot_R(ii,JJ) / (dvR(ii,JJ) - dv0(ii,JJ))
2664521280 if (FA_avg > max(FA_0, FAmt_R(ii,JJ))) then ; FA_avg = max(FA_0, FAmt_R(ii,JJ))
2665518430 elseif (FA_avg < min(FA_0, FAmt_R(ii,JJ))) then ; FA_0 = FA_avg ; endif
2666521280 BT_cont%FA_v_N0(i,J) = FA_0 ; BT_cont%FA_v_NN(i,J) = FAmt_R(ii,JJ)
2667521280 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))) * &
2669367120 ((FAmt_R(ii,JJ) - FA_avg) / (FAmt_R(ii,JJ) - FA_0))
2670 endif
2671 else
26720 BT_cont%FA_v_S0(i,J) = 0.0 ; BT_cont%FA_v_SS(i,J) = 0.0
26730 BT_cont%FA_v_N0(i,J) = 0.0 ; BT_cont%FA_v_NN(i,J) = 0.0
26740 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
26814416end subroutine set_merid_BT_cont
2682
2683!> Calculates left/right edge values for PPM reconstruction.
268413subroutine 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)) :: &
270726 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
271713 local_open_BC = .false.
271813 if (associated(OBC)) then
27190 local_open_BC = OBC%open_u_BCs_exist_globally
2720 endif
2721
272213 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.
272513 stencil = 2 ; if (simple_2nd) stencil = 1
2726
272713 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,".")') &
27300 stencil + max(G%isd-isl,iel-G%ied)
27310 call MOM_error(FATAL,mesg)
2732 endif
273313 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,".")') &
27360 max(G%jsd-jsl,jel-G%jed)
27370 call MOM_error(FATAL,mesg)
2738 endif
2739
2740 !$omp target enter data map(alloc: slp)
2741
2742299 do ks = 1, nz, nkblock
2743286 ke = min(ks + nkblock - 1, nz)
2744
2745286 if (simple_2nd) then
2746 ! untested
27470 do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_im1,h_ip1))
27480 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)
27490 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)
27500 h_W(i,j,k) = 0.5*( h_im1 + h_in(i,j,k) )
27510 h_E(i,j,k) = 0.5*( h_ip1 + h_in(i,j,k) )
2752 enddo
2753 else
2754286 do concurrent (k=ks:ke, j=jsl:jel, i=isl-1:iel+1) DO_LOCALITY(local(dMx,dMn,kk))
275519363344 kk = k - ks + 1
275638831078 if ((G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j)) == 0.0) then
27576206200 slp(i,j,kk) = 0.0
2758 else
2759 ! This uses a simple 2nd order slope.
276013157144 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)
276213157144 dMx = max(h_in(i+1,j,k), h_in(i-1,j,k), h_in(i,j,k)) - h_in(i,j,k)
276313157144 dMn = h_in(i,j,k) - min(h_in(i+1,j,k), h_in(i-1,j,k), h_in(i,j,k))
276413157144 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
2769286 if (local_open_BC) then
2770 ! untested
27710 do n=1, OBC%number_of_segments
27720 segment => OBC%segment(n)
27730 if (.not. segment%on_pe) cycle
27740 if (segment%is_E_or_W) then
27750 I=segment%HI%IsdB
27760 do concurrent (k=ks:ke, j=segment%HI%jsd:segment%HI%jed) DO_LOCALITY(local(kk))
27770 kk = k - ks + 1
27780 slp(i+1,j,kk) = 0.0
27790 slp(i,j,kk) = 0.0
2780 enddo
2781 endif
2782 enddo
2783 endif
2784
2785286 do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_im1,h_ip1,kk))
278619256952 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)
279119256952 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)
279219256952 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) ) + &
279519256952 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) ) + &
279738617722 oneSixth*( slp(i,j,kk) - slp(i+1,j,kk) )
2798 enddo
2799 endif
2800
2801286 if (local_open_BC) then
2802 ! untested
28030 do n=1, OBC%number_of_segments
28040 segment => OBC%segment(n)
28050 if (.not. segment%on_pe) cycle
28060 if (segment%direction == OBC_DIRECTION_E) then
28070 I=segment%HI%IsdB
28080 do concurrent (k=ks:ke, j=segment%HI%jsd:segment%HI%jed)
28090 h_W(i+1,j,k) = h_in(i,j,k)
28100 h_E(i+1,j,k) = h_in(i,j,k)
28110 h_W(i,j,k) = h_in(i,j,k)
28120 h_E(i,j,k) = h_in(i,j,k)
2813 enddo
28140 elseif (segment%direction == OBC_DIRECTION_W) then
28150 I=segment%HI%IsdB
28160 do concurrent (k=ks:ke, j=segment%HI%jsd:segment%HI%jed)
28170 h_W(i,j,k) = h_in(i+1,j,k)
28180 h_E(i,j,k) = h_in(i+1,j,k)
28190 h_W(i+1,j,k) = h_in(i+1,j,k)
28200 h_E(i+1,j,k) = h_in(i+1,j,k)
2821 enddo
2822 endif
2823 enddo
2824 endif
2825
2826299 if (monotonic) then
2827 ! untested
28280 call PPM_limit_CW84(h_in, h_W, h_E, G, GV, isl, iel, jsl, jel, ks, ke)
2829 else
2830286 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
283613end subroutine PPM_reconstruction_x
2837
2838!> Calculates left/right edge values for PPM reconstruction.
283913subroutine 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)) :: &
286226 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
287213 local_open_BC = .false.
287313 if (associated(OBC)) then
28740 local_open_BC = OBC%open_v_BCs_exist_globally
2875 endif
2876
287713 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.
288013 stencil = 2 ; if (simple_2nd) stencil = 1
2881
288213 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,".")') &
28850 max(G%isd-isl,iel-G%ied)
28860 call MOM_error(FATAL,mesg)
2887 endif
288813 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,".")') &
28910 stencil + max(G%jsd-jsl,jel-G%jed)
28920 call MOM_error(FATAL,mesg)
2893 endif
2894
2895 !$omp target enter data map(alloc: slp)
2896
2897299 do ks = 1, nz, nkblock
2898286 ke = min(ks + nkblock - 1, nz)
2899
2900286 if (simple_2nd) then
2901 ! untested
29020 do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_jm1,h_jp1))
29030 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)
29040 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)
29050 h_S(i,j,k) = 0.5*( h_jm1 + h_in(i,j,k) )
29060 h_N(i,j,k) = 0.5*( h_jp1 + h_in(i,j,k) )
2907 enddo
2908 else
2909103246 do concurrent (k=ks:ke, j=jsl-1:jel+1, i=isl:iel) DO_LOCALITY(local(dMx,dMn,kk))
291018944640 kk = k - ks + 1
291137992526 if ((G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1)) == 0.0) then
29125923060 slp(i,j,kk) = 0.0
2913 else
2914 ! This uses a simple 2nd order slope.
291513021580 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)
291713021580 dMx = max(h_in(i,j+1,k), h_in(i,j-1,k), h_in(i,j,k)) - h_in(i,j,k)
291813021580 dMn = h_in(i,j,k) - min(h_in(i,j+1,k), h_in(i,j-1,k), h_in(i,j,k))
291913021580 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
2924286 if (local_open_BC) then
2925 ! untested
29260 do n=1, OBC%number_of_segments
29270 segment => OBC%segment(n)
29280 if (.not. segment%on_pe) cycle
29290 if (segment%is_N_or_S) then
29300 J=segment%HI%JsdB
29310 do concurrent (k=ks:ke, i=segment%HI%isd:segment%HI%ied) DO_LOCALITY(local(kk))
29320 kk = k - ks + 1
29330 slp(i,j+1,kk) = 0.0
29340 slp(i,j,kk) = 0.0
2935 enddo
2936 endif
2937 enddo
2938 endif
2939
2940286 do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_jm1,h_jp1,kk))
294118738720 kk = k - ks + 1
2942 ! Neighboring values should take into account any boundaries. The 3
2943 ! following sets of expressions are equivalent.
294418738720 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)
294518738720 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) ) + &
294818738720 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) ) + &
295037580686 oneSixth*( slp(i,j,kk) - slp(i,j+1,kk) )
2951 enddo
2952 endif
2953
2954286 if (local_open_BC) then
2955 ! untested
29560 do n=1, OBC%number_of_segments
29570 segment => OBC%segment(n)
29580 if (.not. segment%on_pe) cycle
29590 if (segment%direction == OBC_DIRECTION_N) then
29600 J=segment%HI%JsdB
29610 do concurrent (k=ks:ke, i=segment%HI%isd:segment%HI%ied)
29620 h_S(i,j+1,k) = h_in(i,j,k)
29630 h_N(i,j+1,k) = h_in(i,j,k)
29640 h_S(i,j,k) = h_in(i,j,k)
29650 h_N(i,j,k) = h_in(i,j,k)
2966 enddo
29670 elseif (segment%direction == OBC_DIRECTION_S) then
29680 J=segment%HI%JsdB
29690 do concurrent (k=ks:ke, i=segment%HI%isd:segment%HI%ied)
29700 h_S(i,j,k) = h_in(i,j+1,k)
29710 h_N(i,j,k) = h_in(i,j+1,k)
29720 h_S(i,j+1,k) = h_in(i,j+1,k)
29730 h_N(i,j+1,k) = h_in(i,j+1,k)
2974 enddo
2975 endif
2976 enddo
2977 endif
2978
2979299 if (monotonic) then
2980 ! untested
29810 call PPM_limit_CW84(h_in, h_S, h_N, G, GV, isl, iel, jsl, jel, ks, ke)
2982 else
2983286 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
298913end 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.
2995572subroutine 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
301938202736 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.
302237995672 curv = 3.0*((h_L(i,j,k) + h_R(i,j,k)) - 2.0*h_in(i,j,k))
302376198408 if (curv > 0.0) then ! Only minima are limited.
302410827795 dh = h_R(i,j,k) - h_L(i,j,k)
302510827795 if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
30264494182 if (h_in(i,j,k) <= h_min) then
30273619496 h_L(i,j,k) = h_in(i,j,k) ; h_R(i,j,k) = h_in(i,j,k)
3028874686 elseif (12.0*curv*(h_in(i,j,k) - h_min) < (curv**2 + 3.0*dh**2)) then
3029 ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
3030 ! be limited in this case. 0 < scale < 1.
3031142852 scale = 12.0*curv*(h_in(i,j,k) - h_min) / (curv**2 + 3.0*dh**2)
3032142852 h_L(i,j,k) = h_in(i,j,k) + scale*(h_L(i,j,k) - h_in(i,j,k))
3033142852 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
3039572end 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.
30430subroutine 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
30670 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
30700 h_i = h_in(i,j,k)
30710 if ( ( h_R(i,j,k) - h_i ) * ( h_i - h_L(i,j,k) ) <= 0. ) then
30720 h_L(i,j,k) = h_i ; h_R(i,j,k) = h_i
3073 else
30740 RLdiff = h_R(i,j,k) - h_L(i,j,k) ! Difference of edge values
30750 RLmean = 0.5 * ( h_R(i,j,k) + h_L(i,j,k) ) ! Mean of edge values
30760 FunFac = 6. * RLdiff * ( h_i - RLmean ) ! Some funny factor
30770 RLdiff2 = RLdiff * RLdiff ! Square of difference
30780 if ( FunFac > RLdiff2 ) h_L(i,j,k) = 3. * h_i - 2. * h_R(i,j,k)
30790 if ( FunFac < -RLdiff2 ) h_R(i,j,k) = 3. * h_i - 2. * h_L(i,j,k)
3080 endif
3081 enddo
3082
30830end subroutine PPM_limit_CW84
3084
3085!> Return the maximum ratio of a/b or maxrat.
30860pure 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
30920 if (abs(a) > abs(maxrat*b)) then
30930 ratio = maxrat
3094 else
30950 ratio = a / b
3096 endif
30970end function ratio_max
3098
3099!> Initializes continuity_ppm_cs
31001subroutine 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
31311 CS%initialized = .true.
3132
31331 local_open_BC = .false.
31341 if (associated(OBC)) then
31350 local_open_BC = OBC%open_u_BCs_exist_globally
3136 endif
3137
3138! Read all relevant parameters and write them to the model log.
31391 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 "//&
31431 "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 "//&
31511 "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 "//&
31561 "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.5*NK*ANGSTROM, and this should not be set less "//&
3163 "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=GV%m_to_H, &
31641 default=0.5*GV%ke*GV%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 "//&
31691 "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 "//&
31741 "relative CFL change up to 0.5.", default=.false.)
31751 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.", &
31801 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", &
31831 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 "//&
31871 "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 "//&
31901 "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. "//&
31941 "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 "//&
31981 "minimum is 0.", default=.false.)
31991 write(niblock_dflt_str, '(I0)') default_niblock
32001 write(njblock_dflt_str, '(I0)') default_njblock
32011 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 "//&
32051 "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 "//&
32091 "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 "//&
32131 "full vertical column.", default=default_nkblock, layoutParam=.true.)
32141 if (CS%niblock < 0) &
3215 call MOM_error(FATAL, "CONTINUITY_NIBLOCK must be nonnegative; "//&
32160 "use 0 to select the default block size.")
32171 if (CS%njblock < 0) &
3218 call MOM_error(FATAL, "CONTINUITY_NJBLOCK must be nonnegative; "//&
32190 "use 0 to select the default block size.")
32201 if (CS%nkblock < 0) &
3221 call MOM_error(FATAL, "CONTINUITY_NKBLOCK must be nonnegative; "//&
32220 "use 0 to select the default block size.")
32231 CS%diag => diag
3224 !$omp target update to(CS)
3225
32261 id_clock_reconstruct = cpu_clock_id('(Ocean continuity reconstruction)', grain=CLOCK_ROUTINE)
32271 id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=CLOCK_ROUTINE)
32281 id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=CLOCK_ROUTINE)
3229
32301 if (use_h_marg_min) then
32310 CS%h_marg_min = GV%H_subroundoff
3232 else
32331 CS%h_marg_min = 0.
3234 endif
3235
32361 if (local_open_BC) then
32370 do n=1, OBC%number_of_segments
32380 segment => OBC%segment(n)
32390 if (associated(segment%h_Reg)) then
32400 if (.not. allocated(segment%h_Reg%h_res)) then
3241 write(mesg,'("In MOM_continuity_PPM, continuity_PPM_init called with ", &
32420 & "badly configured h_res.")')
32430 call MOM_error(FATAL, mesg)
3244 endif
3245 endif
3246 enddo
3247 endif
3248
32491end subroutine continuity_PPM_init
3250
3251!> continuity_PPM_stencil returns the continuity solver stencil size
325231function 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
325631 stencil = 3 ; if (CS%simple_2nd) stencil = 2 ; if (CS%upwind_1st) stencil = 1
3257
325831end 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.
326126function 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
327426 add_i_stencil = .false. ; if (present(i_stencil)) add_i_stencil = i_stencil
327526 add_j_stencil = .false. ; if (present(j_stencil)) add_j_stencil = j_stencil
3276
327726 stencil = continuity_PPM_stencil(CS)
3278
327926 if (add_i_stencil) then
32800 LB%ish = G%isc-stencil ; LB%ieh = G%iec+stencil
3281 else
328226 LB%ish = G%isc ; LB%ieh = G%iec
3283 endif
3284
328526 if (add_j_stencil) then
328613 LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil
3287 else
328813 LB%jsh = G%jsc ; LB%jeh = G%jec
3289 endif
3290
329126end 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
32990end module MOM_continuity_PPM