← 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).
992880subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &
1007200 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
1574322 real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
1584322 real :: h_E(SZI_(G),SZJ_(G),SZK_(GV)) ! East edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
1594322 real :: h_S(SZI_(G),SZJ_(G),SZK_(GV)) ! South edge thicknesses in the meridional PPM reconstruction [H ~> m or kg m-2]
1604322 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
1672161 niblock = CS%niblock
1682161 njblock = CS%njblock
169
1702161 h_min = GV%Angstrom_H
171
1722161 if (.not.CS%initialized) call MOM_error(FATAL, &
1730 "MOM_continuity_PPM: Module must be initialized before it is used.")
174
1752161 x_first = (MOD(G%first_direction,2) == 0)
176
1772161 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
1832161 if (x_first) then
184
185 ! First advect zonally, with loop bounds that accomodate the subsequent meridional advection.
1862161 LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.true.)
187 ! set whole-domain block sizes when ni/jblock is 0
1882161 if (niblock == 0) niblock = LB%ieh-LB%ish+2
1892161 if (njblock == 0) njblock = LB%jeh-LB%jsh+1
190
1912161 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, &
1932161 niblock, njblock, LB, uhbt, visc_rem_u, u_cor, BT_cont, du_cor)
1942161 call continuity_zonal_convergence(h, uh, dt, G, GV, LB, hin)
195
196 ! Now advect meridionally, using the updated thicknesses to determine the fluxes.
1972161 LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
198
1992161 if (niblock == 0) niblock = LB%ieh-LB%ish+1
2002161 if (njblock == 0) njblock = LB%jeh-LB%jsh+2
201
2022161 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, &
2042161 niblock, njblock, LB, vhbt, visc_rem_v, v_cor, BT_cont, dv_cor)
2052161 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
23012241end 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.
4062161subroutine 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
4232161 call cpu_clock_begin(id_clock_update)
424
4252161 h_min = 0.0 ; if (present(hmin)) h_min = hmin
426
4272161 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
428
4292161 if (present(hin)) then
4302161 do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
43143948993 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
4402161 call cpu_clock_end(id_clock_update)
441
4422161end 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
4622161 call cpu_clock_begin(id_clock_update)
463
4642161 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
465
4662161 h_min = 0.0 ; if (present(hmin)) h_min = hmin
467
4682161 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
4742161 do concurrent (k=1:nz, j=jsh:jeh, i=ish:ieh)
47539954433 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
4792161 call cpu_clock_end(id_clock_update)
480
4812161end subroutine continuity_merdional_convergence
482
483
484!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.
4852161subroutine 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
5042161 call cpu_clock_begin(id_clock_reconstruct)
505
5062161 if (present(LB_in)) then
5072161 LB = LB_in
508 else
5090 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
510 endif
5112161 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
5122161 nkblock = CS%nkblock
5132161 if (nkblock == 0) nkblock = nz
514
5152161 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, &
5212161 nkblock, 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC)
522 endif
523
5242161 call cpu_clock_end(id_clock_reconstruct)
525
5262161end subroutine zonal_edge_thickness
527
528
529!> Set the reconstructed thicknesses at the eastern and western edges of tracer cells.
5302161subroutine 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
5492161 call cpu_clock_begin(id_clock_reconstruct)
550
5512161 if (present(LB_in)) then
5522161 LB = LB_in
553 else
5540 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
555 endif
5562161 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
5572161 nkblock = CS%nkblock
5582161 if (nkblock == 0) nkblock = nz
559
5602161 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, &
5672161 nkblock, 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC)
568 endif
569
5702161 call cpu_clock_end(id_clock_reconstruct)
571
5722161end subroutine meridional_edge_thickness
573
574
575!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.
5766483subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_face_areaU, &
5775040 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)) :: &
6224322 duhdu, & ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
6234322 visc_rem ! A 2-D copy of visc_rem_u or an array of 1's [nondim].
624 real, dimension(niblock,njblock) :: &
6254322 du, & ! Corrective barotropic change in the velocity to give uhbt [L T-1 ~> m s-1].
6264322 du_min_CFL, & ! Lower limit on du correction to avoid CFL violations [L T-1 ~> m s-1]
6274322 du_max_CFL, & ! Upper limit on du correction to avoid CFL violations [L T-1 ~> m s-1]
6284322 duhdu_tot_0, & ! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
6294322 uh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
6304322 visc_rem_max ! The column maximum of visc_rem [nondim].
6314322 logical, dimension(niblock,njblock) :: do_I
632 real, dimension(niblock,njblock,SZK_(GV)) :: &
6334322 uh_t ! A block-sized copy of uh [H L2 T-1 ~> m3 s-1 or kg s-1].
634 real, dimension(niblock,njblock) :: &
6354322 uhbt_t ! A block-sized copy of uhbt [H L2 T-1 ~> m3 s-1 or kg s-1].
6364322 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
6494322 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
6562161 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
6622161 use_visc_rem = present(visc_rem_u)
663
6642161 set_BT_cont = .false. ; if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
665
6662161 local_specified_BC = .false. ; local_Flather_OBC = .false. ; local_open_BC = .false.
6672161 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
6732161 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
6792161 if (present(LB_in)) then
6802161 LB = LB_in
681 else
6820 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
683 endif
6842161 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
685
6862161 CFL_dt = CS%CFL_limit_adjust / dt
6872161 I_dt = 1.0 / dt
6882161 if (CS%aggress_adjust) CFL_dt = I_dt
689
6902161 if (.not.use_visc_rem) then
6911 do concurrent (k=1:nz, jj=1:njblock, ii=1:niblock)
6929761 visc_rem(ii,jj,k) = 1.0
693 enddo
694 endif
695
696105889 do j_start=jsh,jeh,njblock ; do i_start=ish-1,ieh,niblock
69751864 i_end = min(i_start+niblock-1,ieh)
69851864 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
7031218804 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
7045692074 do_I(ii,jj) = .true.
705 enddo
706 ! Set uh and duhdu.
707 !$omp target teams num_teams(nteams)
708377264 do k=1,nz
709372300 if (use_visc_rem) then
710 !$omp loop collapse(2) private(ii,jj)
71144917200 do j=j_start,j_end ; do I=i_start,i_end
71243124400 ii=I-i_start+1 ; jj=j-j_start+1
71344550000 visc_rem(ii,jj,k) = visc_rem_u(I,j,k)
714 enddo ; enddo
715 endif
716 !$omp loop collapse(2) private(ii,jj)
71745541050 do j=j_start,j_end ; do i=i_start,i_end
71843723350 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,&
72243723350 por_face_areaU(I,j,k))
72343723350 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,&
7261445400 OBC%segnum_u(i,j))
727 enddo ; enddo
728377264 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
73851864 if (present(uhbt) .or. set_BT_cont) then
73951840 if (use_visc_rem .and. CS%use_visc_rem_max) then
740 !$omp target
741 !$omp loop collapse(2) private(ii,jj)
7424721760 do j=j_start,j_end ; do i=i_start,i_end
7434471200 ii=i-i_start+1 ; jj=j-j_start+1
7444669920 visc_rem_max(ii,jj) = 0.0
745 enddo ; enddo
746372096 do k=1,nz
747 !$omp loop collapse(2) private(ii,jj)
74844922096 do j=j_start,j_end ; do i=i_start,i_end
74943124400 ii=i-i_start+1 ; jj=j-j_start+1
75044550000 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)
7634721760 do j=j_start,j_end ; do I=i_start,i_end
7644471200 ii=I-i_start+1 ; jj=j-j_start+1
7654471200 I_vrm = 0.0
7664471200 if (visc_rem_max(ii,jj) > 0.0) I_vrm = 1.0 / visc_rem_max(ii,jj)
7674471200 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))
7704471200 else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
7714471200 du_max_CFL(ii,jj) = 2.0* (CFL_dt * dx_W) * I_vrm
7724471200 du_min_CFL(ii,jj) = -2.0 * (CFL_dt * dx_E) * I_vrm
7734669920 uh_tot_0(ii,jj) = 0.0 ; duhdu_tot_0(ii,jj) = 0.0
774 enddo ; enddo
775372096 do k=1,nz
776 !$omp loop collapse(2) private(ii,jj)
77744922096 do j=j_start,j_end ; do I=i_start,i_end
77843124400 ii=I-i_start+1 ; jj=j-j_start+1
77943124400 duhdu_tot_0(ii,jj) = duhdu_tot_0(ii,jj) + duhdu(ii,jj,k)
78044550000 uh_tot_0(ii,jj) = uh_tot_0(ii,jj) + uh_t(ii,jj,k)
781 enddo
782 enddo ; enddo
783 !$omp end target
78451840 if (use_visc_rem) then
78551840 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
808372096 do k=1,nz
809 !$omp loop collapse(2) private(ii,jj,dx_E,dx_W)
81044922096 do j=j_start,j_end ; do I=i_start,i_end
81143124400 ii=I-i_start+1 ; jj=j-j_start+1
81243124400 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))
81543124400 else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif
816
81743124400 if (du_max_CFL(ii,jj) * visc_rem(ii,jj,k) > dx_W*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
8183740193 du_max_CFL(ii,jj) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(ii,jj,k)
81943124400 if (du_min_CFL(ii,jj) * visc_rem(ii,jj,k) < -dx_E*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
8209646385 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
8621218240 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
8634471200 du_max_CFL(ii,jj) = max(du_max_CFL(ii,jj),0.0)
8645689440 du_min_CFL(ii,jj) = min(du_min_CFL(ii,jj),0.0)
865 enddo
866
86751840 any_simple_OBC = .false.
86851840 if (present(uhbt) .or. set_BT_cont) then
86951840 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
8801218240 do concurrent (jj=1:j_end-j_start+1, II=1:i_end-i_start+1)
8815689440 do_I(II,jj) = .true.
882 enddo ; endif
883 endif
884
88551840 if (present(uhbt)) then
88634560 do concurrent (j=j_start:j_end, I=i_start:i_end)
8873792960 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, &
89334560 uh_t, OBC=OBC)
894
89534560 if (present(u_cor)) then
89634560 do concurrent (k=1:nz, j=j_start:j_end, I=i_start:i_end)
89728749600 ii=I-i_start+1 ; jj=j-j_start+1
89828749600 u_cor(I,j,k) = u(I,j,k) + du(ii,jj) * visc_rem(ii,jj,k)
89929234928 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
90534560 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
91451840 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, &
91834560 por_face_areaU, niblock, njblock)
91934560 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
95525932 do concurrent (k=1:nz, j=j_start:j_end, I=i_start:i_end)
95644461453 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
9652161 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
9902161 if (set_BT_cont) then ; if (allocated(BT_cont%h_u)) then
9911440 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, &
993720 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, &
996720 CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU, visc_rem_u)
997 endif
998 endif ; endif
999
10002161 call cpu_clock_end(id_clock_correct)
1001
10025761end 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
1087343010625elemental 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.
1124343010625 tmp = G_dy_Cu * por_face_area ! precalculate things
1125343010625 if (u > 0.0) then
112687385524 if (vol_CFL) then ; CFL = (u * dt) * (G_dy_Cu * G_IareaT)
112787385524 else ; CFL = u * dt * G_IdxT ; endif
112887385524 curv_3 = (h_L + h_R) - 2.0*h
112987385524 dh = h_L - h_R
1130 uh = tmp * u * &
113187385524 (h_R + CFL * (0.5*dh + curv_3*(CFL - 1.5)))
113287385524 h_marg = h_R + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
1133255625101 elseif (u < 0.0) then
1134164484051 if (vol_CFL) then ; CFL = (-u * dt) * (G_dy_Cu * G_IareaT_p1)
1135164484051 else ; CFL = -u * dt * G_IdxT_p1 ; endif
1136164484051 curv_3 = (h_L_p1 + h_R_p1) - 2.0*h_p1
1137164484051 dh = h_R_p1-h_L_p1
1138 uh = tmp * u * &
1139164484051 (h_L_p1 + CFL * (0.5*dh + curv_3*(CFL - 1.5)))
1140164484051 h_marg = h_L_p1 + CFL * (dh + 3.0*curv_3*(CFL - 1.0))
1141 else
114291141050 uh = 0.0
114391141050 h_marg = 0.5 * (h_L_p1 + h_R)
1144 endif
1145343010625 duhdu = tmp * h_marg * visc_rem
1146
1147343010625end 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.
11901440subroutine zonal_flux_thickness(u, h, h_W, h_E, h_u, dt, G, GV, US, LB, vol_CFL, &
11911440 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
12281440 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
1229
12301440 do concurrent (k=1:nz, j=jsh:jeh, I=ish-1:ieh)
123128749600 if (u(I,j,k) > 0.0) then
12327951478 if (vol_CFL) then ; CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j))
12337951478 else ; CFL = u(I,j,k) * dt * G%IdxT(i,j) ; endif
12347951478 curv_3 = (h_W(i,j,k) + h_E(i,j,k)) - 2.0*h(i,j,k)
12357951478 dh = h_W(i,j,k) - h_E(i,j,k)
12367951478 if (marginal) then
12377951478 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
124120798122 elseif (u(I,j,k) < 0.0) then
12429490222 if (vol_CFL) then ; CFL = (-u(I,j,k)*dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j))
12439490222 else ; CFL = -u(I,j,k) * dt * G%IdxT(i+1,j) ; endif
12449490222 curv_3 = (h_W(i+1,j,k) + h_E(i+1,j,k)) - 2.0*h(i+1,j,k)
12459490222 dh = h_E(i+1,j,k)-h_W(i+1,j,k)
12469490222 if (marginal) then
12479490222 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.
125411307900 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
125957888384 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.
126328749600 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
12691440 local_open_BC = .false.
12701440 if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally
12711440 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
13011440end subroutine zonal_flux_thickness
1302
1303!> Returns the barotropic velocity adjustment that gives the
1304!! desired barotropic (layer-summed) transport.
130569120subroutine zonal_flux_adjust(u, h_in, h_W, h_E, uhbt, uh_tot_0, duhdu_tot_0, &
130669120 du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
130769120 i_start, i_end, j_start, j_end, do_I_in, por_face_areaU, niblock, njblock, &
130869120 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) :: &
1357138240 uh_err, & ! Difference between uhbt and the summed uh [H L2 T-1 ~> m3 s-1 or kg s-1].
1358138240 uh_err_best, & ! The smallest value of uh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
1359138240 duhdu_tot,&! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
1360138240 du_min, & ! Lower limit on du correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
1361138240 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
1367138240 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
137469120 nz = GV%ke
137569120 local_open_BC = .false.
137669120 if (present(OBC)) then ; if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally ; endif
1377
137869120 tol_vel = CS%tol_vel
1379
1380 !$omp target
1381
1382 !$omp loop collapse(2) private(ii,jj)
13836295680 do j=j_start,j_end ; do I=i_start,i_end
13845961600 ii=I-i_start+1 ; jj=j-j_start+1
13855961600 du(ii,jj) = 0.0 ; do_I(ii,jj) = do_I_in(ii,jj)
13865961600 du_max(ii,jj) = du_max_CFL(ii,jj) ; du_min(ii,jj) = du_min_CFL(ii,jj)
13875961600 uh_err(ii,jj) = uh_tot_0(ii,jj) - uhbt(ii,jj) ; duhdu_tot(ii,jj) = duhdu_tot_0(ii,jj)
13886226560 uh_err_best(ii,jj) = abs(uh_err(ii,jj))
1389 enddo ; enddo
1390
1391199097 do itt=1,max_itts
139269120 select case (itt)
139369120 case (:1) ; tol_eta = 1e-6 * CS%tol_eta
139463338 case (2) ; tol_eta = 1e-4 * CS%tol_eta
139563338 case (3) ; tol_eta = 1e-2 * CS%tol_eta
1396199097 case default ; tol_eta = CS%tol_eta
1397 end select
1398
1399 !$omp loop collapse(2) private(ii,jj)
140018705701 do j=j_start,j_end ; do I=i_start,i_end
140117721736 ii=I-i_start+1 ; jj=j-j_start+1
140218506604 if (uh_err(ii,jj) > 0.0) then ; du_max(ii,jj) = du(ii,jj)
140310516469 elseif (uh_err(ii,jj) < 0.0) then ; du_min(ii,jj) = du(ii,jj)
14043363892 else ; do_I(ii,jj) = .false. ; endif
1405 enddo ; enddo
1406#ifndef __NVCOMPILER_OPENMP_GPU
1407199097 domore = .false.
1408#endif
1409 !$omp loop collapse(2) private(ii,jj,ddu,du_prev)
141018705701 do j=j_start,j_end ; do I=i_start,i_end
141117721736 ii=I-i_start+1 ; jj=j-j_start+1
141218506604 if (do_I(ii,jj)) then
141314131011 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.
14189679941 ddu = -uh_err(ii,jj) / duhdu_tot(ii,jj)
14199679941 du_prev = du(ii,jj)
14209679941 du(ii,jj) = du(ii,jj) + ddu
14219679941 if (abs(ddu) < 1.0e-15*abs(du(ii,jj))) then
14221 do_I(ii,jj) = .false. ! ddu is small enough to quit.
14239679940 elseif (ddu > 0.0) then
14244755297 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
14294924643 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
14359679941 if (do_I(ii,jj)) domore = .true.
1436#endif
1437 else
14384451070 do_I(ii,jj) = .false.
1439 endif
1440 endif
1441 enddo ; enddo
1442#ifndef __NVCOMPILER_OPENMP_GPU
1443199097 if (.not.domore) exit
1444#endif
1445
1446 !$omp loop collapse(2) private(ii,jj)
144712410021 do j=j_start,j_end ; do I=i_start,i_end
144811760136 ii=I-i_start+1 ; jj=j-j_start+1
144912280044 uh_err(ii,jj) = -uhbt(ii,jj) ; duhdu_tot(ii,jj) = 0.0
1450 enddo ; enddo
1451906072 do k=1,nz
1452 !$omp loop collapse(2) private(ii,jj,duhdu,u_new)
1453112629672 do j=j_start,j_end ; do i=i_start,i_end
1454108147000 ii=i-i_start+1 ; jj=j-j_start+1
1455111723600 if (do_I(ii,jj)) then
145650392650 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,&
146050392650 por_face_areaU(i,j,k))
146150392650 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))
146550392650 uh_err(ii,jj) = uh_err(ii,jj) + uh_3d(ii,jj,k)
146650392650 duhdu_tot(ii,jj) = duhdu_tot(ii,jj) + duhdu
1467 endif
1468 enddo ; enddo
1469 enddo
1470 !$omp loop collapse(2) private(ii,jj)
147112410021 do j=j_start,j_end ; do I=i_start,i_end
147211760136 ii=I-i_start+1 ; jj=j-j_start+1
147312280044 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
148469120end 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.
148934560subroutine set_zonal_BT_cont(u, h_in, h_W, h_E, BT_cont, uh_tot_0, duhdu_tot_0, &
149034560 du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
149134560 visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
149234560 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) :: &
153369120 du0 ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
1534 real, dimension(niblock,njblock) :: &
153569120 zeros ! An array of full of 0 transports [H L2 T-1 ~> m3 s-1 or kg s-1]
1536 real, dimension(niblock,njblock) :: &
153769120 duL, duR, & ! The barotropic velocity increments that give the westerly
1538 ! (duL) and easterly (duR) test velocities [L T-1 ~> m s-1].
153969120 du_CFL, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
154069120 FAmt_L, FAmt_R, & ! The summed effective marginal face areas for the 3
154169120 FAmt_0, & ! test velocities [H L ~> m2 or kg m-1].
154269120 uhtot_L, & ! The summed transport with the westerly (uhtot_L) and
154369120 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].
156669120 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
157134560 nz = GV%ke ; Idt = 1.0 / dt
157234560 min_visc_rem = 0.1 ; CFL_min = 1e-6
1573
1574 ! Diagnose the zero-transport correction, du0.
1575812160 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
15763792960 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, &
158034560 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)
15873147840 do j=j_start,j_end ; do i=i_start,i_end
15882980800 ii=i-i_start+1 ; jj=j-j_start+1
15892980800 du_CFL(ii,jj) = (CFL_min * Idt) * G%dxCu(i,j)
15902980800 duR(ii,jj) = min(0.0,du0(ii,jj) - du_CFL(ii,jj))
15912980800 duL(ii,jj) = max(0.0,du0(ii,jj) + du_CFL(ii,jj))
15922980800 FAmt_L(ii,jj) = 0.0 ; FAmt_R(ii,jj) = 0.0 ; FAmt_0(ii,jj) = 0.0
15933113280 uhtot_L(ii,jj) = 0.0 ; uhtot_R(ii,jj) = 0.0
1594 enddo ; enddo
1595
1596248064 do k=1,nz
1597 !$omp loop collapse(2) private(ii,jj,visc_rem_lim)
159829948064 do j=j_start,j_end ; do i=i_start,i_end
159928749600 ii=i-i_start+1 ; jj=j-j_start+1
160029700000 if (do_I(ii,jj)) then
160128749600 visc_rem_lim = max(visc_rem(ii,jj,k), min_visc_rem*visc_rem_max(ii,jj))
160228749600 if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
160317812800 if (u(i,j,k) + duR(ii,jj)*visc_rem_lim > -du_CFL(ii,jj)*visc_rem(ii,jj,k)) &
16041173836 duR(ii,jj) = -(u(i,j,k) + du_CFL(ii,jj)*visc_rem(ii,jj,k)) / visc_rem_lim
160517812800 if (u(i,j,k) + duL(ii,jj)*visc_rem_lim < du_CFL(ii,jj)*visc_rem(ii,jj,k)) &
16063688842 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
1614248064 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)
161629948064 do j=j_start,j_end ; do i=i_start,i_end
161728749600 ii=i-i_start+1 ; jj=j-j_start+1
161829700000 if (do_I(ii,jj)) then
161928749600 u_L = u(i,j,k) + duL(ii,jj) * visc_rem(ii,jj,k)
162028749600 u_R = u(i,j,k) + duR(ii,jj) * visc_rem(ii,jj,k)
162128749600 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,&
162528749600 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,&
162928749600 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,&
163328749600 por_face_areaU(i,j,k))
163428749600 FAmt_0(ii,jj) = FAmt_0(ii,jj) + duhdu_0
163528749600 FAmt_L(ii,jj) = FAmt_L(ii,jj) + duhdu_L
163628749600 FAmt_R(ii,jj) = FAmt_R(ii,jj) + duhdu_R
163728749600 uhtot_L(ii,jj) = uhtot_L(ii,jj) + uh_L
163828749600 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)
16433147840 do j=j_start,j_end ; do i=i_start,i_end
16442980800 ii=i-i_start+1 ; jj=j-j_start+1
16453113280 if (do_I(ii,jj)) then
16462980800 FA_0 = FAmt_0(ii,jj) ; FA_avg = FAmt_0(ii,jj)
16472980800 if ((duL(ii,jj) - du0(ii,jj)) /= 0.0) &
16482980800 FA_avg = uhtot_L(ii,jj) / (duL(ii,jj) - du0(ii,jj))
16492980800 if (FA_avg > max(FA_0, FAmt_L(ii,jj))) then ; FA_avg = max(FA_0, FAmt_L(ii,jj))
16502975511 elseif (FA_avg < min(FA_0, FAmt_L(ii,jj))) then ; FA_0 = FA_avg ; endif
1651
16522980800 BT_cont%FA_u_W0(i,j) = FA_0 ; BT_cont%FA_u_WW(i,j) = FAmt_L(ii,jj)
16532980800 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))) * &
16552456859 ((FAmt_L(ii,jj) - FA_avg) / (FAmt_L(ii,jj) - FA_0))
1656 endif
1657
16582980800 FA_0 = FAmt_0(ii,jj) ; FA_avg = FAmt_0(ii,jj)
16592980800 if ((duR(ii,jj) - du0(ii,jj)) /= 0.0) &
16602980800 FA_avg = uhtot_R(ii,jj) / (duR(ii,jj) - du0(ii,jj))
16612980800 if (FA_avg > max(FA_0, FAmt_R(ii,jj))) then ; FA_avg = max(FA_0, FAmt_R(ii,jj))
16622976412 elseif (FA_avg < min(FA_0, FAmt_R(ii,jj))) then ; FA_0 = FA_avg ; endif
1663
16642980800 BT_cont%FA_u_E0(i,j) = FA_0 ; BT_cont%FA_u_EE(i,j) = FAmt_R(ii,jj)
16652980800 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))) * &
16672464069 ((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
168034560end subroutine set_zonal_BT_cont
1681
1682!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.
16836483subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, por_face_areaV, &
16841440 niblock, njblock, LB_in, vhbt, &
16853600 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) :: &
17284322 vhbt_t
1729 real, dimension(niblock,njblock,SZK_(GV)) :: &
17304322 vh_t
1731 real, dimension(niblock,njblock,SZK_(GV)) :: &
17324322 dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
1733 real, dimension(niblock,njblock) :: &
17344322 dv, & ! Corrective barotropic change in the velocity to give vhbt [L T-1 ~> m s-1].
17354322 dv_min_CFL, & ! Lower limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]
17364322 dv_max_CFL, & ! Upper limit on dv correction to avoid CFL violations [L T-1 ~> m s-1]
17374322 dvhdv_tot_0, & ! Summed partial derivative of vh with v [H L ~> m2 or kg m-1].
17384322 vh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
17394322 visc_rem_max ! The column maximum of visc_rem [nondim]
17404322 logical, dimension(niblock,njblock) :: do_I
17414322 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)) :: &
17454322 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
17584322 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
17632161 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
17692161 use_visc_rem = present(visc_rem_v)
1770
17712161 set_BT_cont = .false. ; if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
1772
17732161 local_specified_BC = .false. ; local_Flather_OBC = .false. ; local_open_BC = .false.
17742161 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
17802161 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
17862161 if (present(LB_in)) then
17872161 LB = LB_in
1788 else
17890 LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec
1790 endif
17912161 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
1792
17932161 CFL_dt = CS%CFL_limit_adjust / dt
17942161 I_dt = 1.0 / dt
17952161 if (CS%aggress_adjust) CFL_dt = I_dt
1796
17972161 if (.not.use_visc_rem) then
17981 do concurrent (k=1:nz, jj=1:njblock, ii=1:niblock)
17999761 visc_rem(ii,jj,k) = 1.0
1800 enddo
1801 endif
1802
180397245 do j_start = jsh-1, jeh, njblock ; do i_start = ish, ieh, niblock
180447542 j_end = min(j_start + njblock-1, jeh)
180547542 i_end = min(i_start + niblock-1, ieh)
1806 !$ nteams = ceiling(real((j_end-j_start+1)*(i_end-i_start+1))/128.)
18071093466 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
18084991910 do_I(ii,JJ) = .true.
1809 enddo
1810 ! This sets vh and dvhdv.
1811 !$omp target teams num_teams(nteams)
1812355072 do k=1,nz
1813350400 if (use_visc_rem) then
1814 !$omp loop collapse(2) private(ii,jj)
181541191200 do J=j_start,j_end ; do i=i_start,i_end
181639528000 ii=i-i_start+1 ; jj=j-j_start+1
181740845600 visc_rem(ii,JJ,k) = visc_rem_v(i,J,k)
1818 enddo ; enddo
1819 endif
1820 !$omp loop collapse(2) private(ii,JJ)
182141763300 do J=j_start,j_end ; do i=i_start,i_end
182240077000 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,&
182640077000 por_face_areaV(i,J,k))
182740077000 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,&
18301335900 OBC%segnum_v(i,J))
1831 enddo ; enddo
1832355072 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
184247542 if (present(vhbt) .or. set_BT_cont) then
184347520 if (use_visc_rem .and. CS%use_visc_rem_max) then
1844 !$omp target
1845 !$omp loop collapse(2) private(ii,jj)
18464121280 do j=j_start,j_end ; do i=i_start,i_end
18473896640 ii=i-i_start+1 ; jj=j-j_start+1
18484073760 visc_rem_max(ii,JJ) = 0.0
1849 enddo ; enddo
1850350208 do k=1,nz
1851 !$omp loop collapse(2) private(ii,jj)
185241195808 do j=j_start,j_end ; do i=i_start,i_end
185339528000 ii=i-i_start+1 ; jj=j-j_start+1
185440845600 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)
18674121280 do J=J_start,J_end ; do i=i_start,i_end
18683896640 ii=i-i_start+1 ; jj=j-j_start+1
18693896640 I_vrm = 0.0
18703896640 if (visc_rem_max(ii,JJ) > 0.0) I_vrm = 1.0 / visc_rem_max(ii,JJ)
18713896640 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))
18743896640 else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
18753896640 dv_max_CFL(ii,JJ) = 2.0 * (CFL_dt * dy_S) * I_vrm
18763896640 dv_min_CFL(ii,JJ) = -2.0 * (CFL_dt * dy_N) * I_vrm
18774073760 vh_tot_0(ii,JJ) = 0.0 ; dvhdv_tot_0(ii,JJ) = 0.0
1878 enddo ; enddo
1879350208 do k=1,nz
1880 !$omp loop collapse(2) private(ii,jj)
188141195808 do J=J_start,J_end ; do i=i_start,i_end
188239528000 ii=i-i_start+1 ; jj=j-j_start+1
188339528000 dvhdv_tot_0(ii,JJ) = dvhdv_tot_0(ii,JJ) + dvhdv(ii,jj,k)
188440845600 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
188947520 if (use_visc_rem) then
189047520 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
1912350208 do k=1,nz
1913 !$omp loop collapse(2) private(ii,jj,dy_S,dy_N)
191441195808 do J=J_start,J_end ; do i=i_start,i_end
191539528000 ii=i-i_start+1 ; jj=j-j_start+1
191639528000 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))
191939528000 else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
192039528000 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)
192339528000 if (dv_min_CFL(ii,JJ) * visc_rem(ii,JJ,k) < &
1924 -dy_N*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
19259948917 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
19651092960 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
19663896640 dv_max_CFL(ii,JJ) = max(dv_max_CFL(ii,JJ),0.0)
19674989600 dv_min_CFL(ii,JJ) = min(dv_min_CFL(ii,JJ),0.0)
1968 enddo
1969
197047520 any_simple_OBC = .false.
197147520 if (present(vhbt) .or. set_BT_cont) then
197247520 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
19831092960 do concurrent (JJ=1:J_end-J_start+1, ii=1:i_end-i_start+1)
19844989600 do_I(ii,JJ) = .true.
1985 enddo ; endif
1986 endif
1987
198847520 if (present(vhbt)) then
198931680 do concurrent (j=j_start:j_end, i=i_start:i_end)
19902597760 ii=i-i_start+1 ; jj=j-j_start+1
19913326400 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, &
199731680 niblock, njblock, vh_t, OBC)
1998
199931680 if (present(v_cor)) then
200031680 do concurrent (k=1:nz, J=J_start:J_end, i=i_start:i_end)
200126352000 ii=i-i_start+1 ; jj=j-j_start+1
200226352000 v_cor(i,J,k) = v(i,J,k) + dv(ii,JJ) * visc_rem(ii,JJ,k)
200326798592 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
200931680 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
201747520 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, &
202131680 por_face_areaV, niblock, njblock)
202231680 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
205823771 do concurrent (k=1:nz, j=j_start:j_end, i=i_start:i_end)
205940077000 jj=j-j_start+1 ; ii=i-i_start+1
206040756192 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
20692161 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
20942161 if (set_BT_cont) then ; if (allocated(BT_cont%h_v)) then
20951440 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, &
2097720 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, &
2100720 CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaV, visc_rem_v)
2101 endif
2102 endif ; endif
2103
21042161 call cpu_clock_end(id_clock_correct)
2105
21065761end 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.
21911440subroutine meridional_flux_thickness(v, h, h_S, h_N, h_v, dt, G, GV, US, LB, vol_CFL, &
21921440 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
22311440 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke
2232
22331440 do concurrent (k=1:nz, J=jsh-1:jeh, i=ish:ieh)
223426352000 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
224420414715 elseif (v(i,J,k) < 0.0) then
224511314065 if (vol_CFL) then ; CFL = (-v(i,J,k)*dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1))
224611314065 else ; CFL = -v(i,J,k) * dt * G%IdyT(i,j+1) ; endif
224711314065 curv_3 = (h_S(i,j+1,k) + h_N(i,j+1,k)) - 2.0*h(i,j+1,k)
224811314065 dh = h_N(i,j+1,k)-h_S(i,j+1,k)
224911314065 if (marginal) then
225011314065 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.
22579100650 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
226253061168 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.
226626352000 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
22721440 local_open_BC = .false.
22731440 if (associated(OBC)) local_open_BC = OBC%open_v_BCs_exist_globally
2274 ! untested - will need to be refactored to be performant on GPUs
22751440 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
23041440end subroutine meridional_flux_thickness
2305
2306
2307!> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.
230863360subroutine meridional_flux_adjust(v, h_in, h_S, h_N, vhbt, vh_tot_0, dvhdv_tot_0, &
230963360 dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
2310126720 i_start, i_end, j_start, j_end, do_I_in, por_face_areaV, &
231163360 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) :: &
2364126720 vh_err, & ! Difference between vhbt and the summed vh [H L2 T-1 ~> m3 s-1 or kg s-1].
2365126720 vh_err_best, & ! The smallest value of vh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
2366126720 dvhdv_tot,&! Summed partial derivative of vh with u [H L ~> m2 or kg m-1].
2367126720 dv_min, & ! Lower limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
2368126720 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
2376126720 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
238263360 local_open_BC = .false.
238363360 if (present(OBC)) then ; if (associated(OBC)) then
23840 local_open_BC = OBC%open_v_BCs_exist_globally
2385 endif ; endif
2386
238763360 nz = GV%ke
238863360 tol_vel = CS%tol_vel
2389
2390 !$omp target
2391
2392 !$omp loop collapse(2) private(ii,jj)
23935495040 do J=J_start,J_end ; do i=i_start,i_end
23945195520 ii=i-i_start+1 ; jj=j-j_start+1
23955195520 dv(ii,JJ) = 0.0 ; do_I(ii,JJ) = do_I_in(ii,JJ)
23965195520 dv_max(ii,JJ) = dv_max_CFL(ii,JJ) ; dv_min(ii,JJ) = dv_min_CFL(ii,JJ)
23975195520 vh_err(ii,JJ) = vh_tot_0(ii,JJ) - vhbt(ii,JJ) ; dvhdv_tot(ii,JJ) = dvhdv_tot_0(ii,JJ)
23985431680 vh_err_best(ii,JJ) = abs(vh_err(ii,JJ))
2399 enddo ; enddo
2400
2401179503 do itt=1,max_itts
240263360 select case (itt)
240363360 case (:1) ; tol_eta = 1e-6 * CS%tol_eta
240457580 case (2) ; tol_eta = 1e-4 * CS%tol_eta
240557556 case (3) ; tol_eta = 1e-2 * CS%tol_eta
2406179503 case default ; tol_eta = CS%tol_eta
2407 end select
2408
2409 !$omp loop collapse(2) private(ii,jj)
241016336699 do j=j_start,j_end ; do i=i_start,i_end
241115456464 ii=i-i_start+1 ; jj=j-j_start+1
241216157196 if (vh_err(ii,JJ) > 0.0) then ; dv_max(ii,JJ) = dv(ii,JJ)
24136602611 elseif (vh_err(ii,JJ) < 0.0) then ; dv_min(ii,JJ) = dv(ii,JJ)
24141224361 else ; do_I(ii,JJ) = .false. ; endif
2415 enddo ; enddo
2416#ifndef __NVCOMPILER_OPENMP_GPU
2417179503 domore = .false.
2418#endif
2419 !$omp loop collapse(2) private(ii,jj,ddv,dv_prev)
242016336699 do J=J_start,J_end ; do i=i_start,i_end
242116157196 ii=i-i_start+1 ; jj=j-j_start+1 ; if (do_I(ii,JJ)) then
242213929428 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.
24279572599 ddv = -vh_err(ii,JJ) / dvhdv_tot(ii,JJ)
24289572599 dv_prev = dv(ii,JJ)
24299572599 dv(ii,JJ) = dv(ii,JJ) + ddv
24309572599 if (abs(ddv) < 1.0e-15*abs(dv(ii,JJ))) then
24310 do_I(ii,JJ) = .false. ! ddv is small enough to quit.
24329572599 elseif (ddv > 0.0) then
24333020166 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
24386552433 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
24449572599 if (do_I(ii,JJ)) domore = .true.
2445#endif
2446 else
24474356829 do_I(ii,JJ) = .false.
2448 endif
2449 endif ; enddo ; enddo
2450#ifndef __NVCOMPILER_OPENMP_GPU
2451179503 if (.not.domore) exit
2452#endif
2453 !$omp loop collapse(2) private(ii,jj)
245410841659 do J=J_start,J_end ; do i=i_start,i_end
245510260944 ii=i-i_start+1 ; jj=J-j_start+1
245610725516 vh_err(ii,JJ) = -vhbt(ii,JJ) ; dvhdv_tot(ii,JJ) = 0.0
2457 enddo ; enddo
2458837368 do k=1,nz
2459 !$omp loop collapse(2) private(ii,jj,dvhdv,v_new)
2460103776368 do J=J_start,J_end ; do i=i_start,i_end
246199633600 ii=i-i_start+1 ; jj=j-j_start+1
2462102939000 if (do_I(ii,JJ)) then
246343512825 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,&
246743512825 por_face_areaV(i,J,k))
246843512825 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))
247243512825 vh_err(ii,JJ) = vh_err(ii,JJ) + vh_3d(ii,JJ,k)
247343512825 dvhdv_tot(ii,JJ) = dvhdv_tot(ii,JJ) + dvhdv
2474 endif
2475 enddo ; enddo
2476 enddo
2477 !$omp loop collapse(2) private(ii,jj)
247810841659 do J=J_start,J_end ; do i=i_start,i_end
247910260944 ii=i-i_start+1 ; jj=J-j_start+1
248010725516 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
249163360end 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.
249631680subroutine set_merid_BT_cont(v, h_in, h_S, h_N, BT_cont, vh_tot_0, dvhdv_tot_0, &
249731680 dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
249831680 visc_rem_max, i_start, i_end, j_start, j_end, do_I, &
249931680 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) :: &
254363360 dv0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
254463360 dvL, dvR, & ! The barotropic velocity increments that give the southerly
2545 ! (dvL) and northerly (dvR) test velocities [L T-1 ~> m s-1].
254663360 dv_CFL, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
254763360 FAmt_L, FAmt_R, & ! The summed effective marginal face areas for the 3
254863360 FAmt_0, & ! test velocities [H L ~> m2 or kg m-1].
254963360 vhtot_L, & ! The summed transport with the southerly (vhtot_L) and
255063360 vhtot_R ! and northerly (vhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
2551 real, dimension(niblock,njblock) :: &
255263360 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
257563360 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
258031680 nz = GV%ke ; Idt = 1.0 / dt
258131680 min_visc_rem = 0.1 ; CFL_min = 1e-6
2582
2583 ! Diagnose the zero-transport correction, dv0.
2584728640 do concurrent (jj=1:j_end-j_start+1, ii=1:i_end-i_start+1)
25853326400 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, &
259031680 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)
25972747520 do J=J_start,J_end ; do i=i_start,i_end
25982597760 ii=i-i_start+1 ; jj=j-j_start+1
25992597760 dv_CFL(ii,JJ) = (CFL_min * Idt) * G%dyCv(i,J)
26002597760 dvR(ii,JJ) = min(0.0,dv0(ii,JJ) - dv_CFL(ii,JJ))
26012597760 dvL(ii,JJ) = max(0.0,dv0(ii,JJ) + dv_CFL(ii,JJ))
26022597760 FAmt_L(ii,JJ) = 0.0 ; FAmt_R(ii,JJ) = 0.0 ; FAmt_0(ii,JJ) = 0.0
26032715840 vhtot_L(ii,JJ) = 0.0 ; vhtot_R(ii,JJ) = 0.0
2604 enddo ; enddo
2605
2606233472 do k=1,nz
2607 !$omp loop collapse(2) private(ii,jj,visc_rem_lim)
260827463872 do J=J_start,J_end ; do i=i_start,i_end
260926352000 ii=i-i_start+1 ; jj=j-j_start+1
261027230400 if (do_I(ii,jj)) then
261126352000 visc_rem_lim = max(visc_rem(ii,JJ,k), min_visc_rem*visc_rem_max(ii,JJ))
261226352000 if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
261317618400 if (v(i,J,k) + dvR(ii,JJ)*visc_rem_lim > -dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) &
26141097929 dvR(ii,JJ) = -(v(i,J,k) + dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) / visc_rem_lim
261517618400 if (v(i,J,k) + dvL(ii,JJ)*visc_rem_lim < dv_CFL(ii,JJ)*visc_rem(ii,JJ,k)) &
26163767578 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
2621233472 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)
262327463872 do J=J_start,J_end ; do i=i_start,i_end
262426352000 ii=i-i_start+1 ; jj=j-j_start+1
262527230400 if (do_I(ii,jj)) then
262626352000 v_L = v(i,J,k) + dvL(ii,JJ) * visc_rem(ii,JJ,k)
262726352000 v_R = v(i,J,k) + dvR(ii,JJ) * visc_rem(ii,JJ,k)
262826352000 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),&
263126352000 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),&
263426352000 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),&
263726352000 G%IareaT(i,J+1),G%IdyT(i,J),G%IdyT(i,J+1),dt,CS%vol_CFL,por_face_areaV(i,J,k))
263826352000 FAmt_0(ii,JJ) = FAmt_0(ii,JJ) + dvhdv_0
263926352000 FAmt_L(ii,JJ) = FAmt_L(ii,JJ) + dvhdv_L
264026352000 FAmt_R(ii,JJ) = FAmt_R(ii,JJ) + dvhdv_R
264126352000 vhtot_L(ii,JJ) = vhtot_L(ii,JJ) + vh_L
264226352000 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)
26472747520 do J=J_start,J_end ; do i=i_start,i_end
26482597760 ii=i-i_start+1 ; jj=j-j_start+1
26492715840 if (do_I(ii,JJ)) then
26502597760 FA_0 = FAmt_0(ii,JJ) ; FA_avg = FAmt_0(ii,JJ)
26512597760 if ((dvL(ii,JJ) - dv0(ii,JJ)) /= 0.0) &
26522597760 FA_avg = vhtot_L(ii,JJ) / (dvL(ii,JJ) - dv0(ii,JJ))
26532597760 if (FA_avg > max(FA_0, FAmt_L(ii,JJ))) then ; FA_avg = max(FA_0, FAmt_L(ii,JJ))
26542593024 elseif (FA_avg < min(FA_0, FAmt_L(ii,JJ))) then ; FA_0 = FA_avg ; endif
26552597760 BT_cont%FA_v_S0(i,J) = FA_0 ; BT_cont%FA_v_SS(i,J) = FAmt_L(ii,JJ)
26562597760 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))) * &
26582462370 ((FAmt_L(ii,JJ) - FA_avg) / (FAmt_L(ii,JJ) - FA_0))
2659 endif
2660
26612597760 FA_0 = FAmt_0(ii,JJ) ; FA_avg = FAmt_0(ii,JJ)
26622597760 if ((dvR(ii,JJ) - dv0(ii,JJ)) /= 0.0) &
26632597760 FA_avg = vhtot_R(ii,JJ) / (dvR(ii,JJ) - dv0(ii,JJ))
26642597760 if (FA_avg > max(FA_0, FAmt_R(ii,JJ))) then ; FA_avg = max(FA_0, FAmt_R(ii,JJ))
26652594208 elseif (FA_avg < min(FA_0, FAmt_R(ii,JJ))) then ; FA_0 = FA_avg ; endif
26662597760 BT_cont%FA_v_N0(i,J) = FA_0 ; BT_cont%FA_v_NN(i,J) = FAmt_R(ii,JJ)
26672597760 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))) * &
26692463605 ((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
268131680end subroutine set_merid_BT_cont
2682
2683!> Calculates left/right edge values for PPM reconstruction.
26842161subroutine 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)) :: &
27074322 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
27172161 local_open_BC = .false.
27182161 if (associated(OBC)) then
27190 local_open_BC = OBC%open_u_BCs_exist_globally
2720 endif
2721
27222161 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.
27252161 stencil = 2 ; if (simple_2nd) stencil = 1
2726
27272161 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
27332161 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
27426483 do ks = 1, nz, nkblock
27435475 ke = min(ks + nkblock - 1, nz)
2744
27455475 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
27545475 do concurrent (k=ks:ke, j=jsl:jel, i=isl-1:iel+1) DO_LOCALITY(local(dMx,dMn,kk))
275544807400 kk = k - ks + 1
275690299175 if ((G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j)) == 0.0) then
275717881350 slp(i,j,kk) = 0.0
2758 else
2759 ! This uses a simple 2nd order slope.
276026926050 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)
276226926050 dMx = max(h_in(i+1,j,k), h_in(i-1,j,k), h_in(i,j,k)) - h_in(i,j,k)
276326926050 dMn = h_in(i,j,k) - min(h_in(i+1,j,k), h_in(i-1,j,k), h_in(i,j,k))
276426926050 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
27695475 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
27855475 do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_im1,h_ip1,kk))
278644084700 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)
279144084700 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)
279244084700 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) ) + &
279544084700 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) ) + &
279788842825 oneSixth*( slp(i,j,kk) - slp(i+1,j,kk) )
2798 enddo
2799 endif
2800
28015475 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
28266483 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
28305475 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
28362161end subroutine PPM_reconstruction_x
2837
2838!> Calculates left/right edge values for PPM reconstruction.
28392161subroutine 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)) :: &
28624322 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
28722161 local_open_BC = .false.
28732161 if (associated(OBC)) then
28740 local_open_BC = OBC%open_v_BCs_exist_globally
2875 endif
2876
28772161 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.
28802161 stencil = 2 ; if (simple_2nd) stencil = 1
2881
28822161 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
28882161 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
28976483 do ks = 1, nz, nkblock
28985475 ke = min(ks + nkblock - 1, nz)
2899
29005475 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
2909662475 do concurrent (k=ks:ke, j=jsl-1:jel+1, i=isl:iel) DO_LOCALITY(local(dMx,dMn,kk))
291042048000 kk = k - ks + 1
291184758475 if ((G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1)) == 0.0) then
291215910350 slp(i,j,kk) = 0.0
2913 else
2914 ! This uses a simple 2nd order slope.
291526137650 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)
291726137650 dMx = max(h_in(i,j+1,k), h_in(i,j-1,k), h_in(i,j,k)) - h_in(i,j,k)
291826137650 dMn = h_in(i,j,k) - min(h_in(i,j+1,k), h_in(i,j-1,k), h_in(i,j,k))
291926137650 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
29245475 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
29405475 do concurrent (k=ks:ke, j=jsl:jel, i=isl:iel) DO_LOCALITY(local(h_jm1,h_jp1,kk))
294140734000 kk = k - ks + 1
2942 ! Neighboring values should take into account any boundaries. The 3
2943 ! following sets of expressions are equivalent.
294440734000 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)
294540734000 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) ) + &
294840734000 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) ) + &
295082130475 oneSixth*( slp(i,j,kk) - slp(i,j+1,kk) )
2951 enddo
2952 endif
2953
29545475 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
29796483 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
29835475 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
29892161end 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.
299510950subroutine 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
301986154600 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.
302284818700 curv = 3.0*((h_L(i,j,k) + h_R(i,j,k)) - 2.0*h_in(i,j,k))
3023170973300 if (curv > 0.0) then ! Only minima are limited.
302425644351 dh = h_R(i,j,k) - h_L(i,j,k)
302525644351 if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
302618775254 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)
302818775254 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.
30312160655 scale = 12.0*curv*(h_in(i,j,k) - h_min) / (curv**2 + 3.0*dh**2)
30322160655 h_L(i,j,k) = h_in(i,j,k) + scale*(h_L(i,j,k) - h_in(i,j,k))
30332160655 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
303910950end 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
32525043function 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
32565043 stencil = 3 ; if (CS%simple_2nd) stencil = 2 ; if (CS%upwind_1st) stencil = 1
3257
32585043end 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.
32614322function 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
32744322 add_i_stencil = .false. ; if (present(i_stencil)) add_i_stencil = i_stencil
32754322 add_j_stencil = .false. ; if (present(j_stencil)) add_j_stencil = j_stencil
3276
32774322 stencil = continuity_PPM_stencil(CS)
3278
32794322 if (add_i_stencil) then
32800 LB%ish = G%isc-stencil ; LB%ieh = G%iec+stencil
3281 else
32824322 LB%ish = G%isc ; LB%ieh = G%iec
3283 endif
3284
32854322 if (add_j_stencil) then
32862161 LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil
3287 else
32882161 LB%jsh = G%jsc ; LB%jeh = G%jec
3289 endif
3290
32914322end 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