MOM6
MOM_continuity_PPM.F90
1 !> Solve the layer continuity equation using the PPM method for layer fluxes.
2 module mom_continuity_ppm
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
7 use mom_diag_mediator, only : time_type, diag_ctrl
8 use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
9 use mom_file_parser, only : get_param, log_version, param_file_type
10 use mom_grid, only : ocean_grid_type
11 use mom_open_boundary, only : ocean_obc_type, obc_segment_type, obc_none
12 use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s
13 use mom_unit_scaling, only : unit_scale_type
14 use mom_variables, only : bt_cont_type
15 use mom_verticalgrid, only : verticalgrid_type
16 
17 implicit none ; private
18 
19 #include <MOM_memory.h>
20 
21 public continuity_ppm, continuity_ppm_init, continuity_ppm_end, continuity_ppm_stencil
22 
23 !>@{ CPU time clock IDs
24 integer :: id_clock_update, id_clock_correct
25 !!@}
26 
27 !> Control structure for mom_continuity_ppm
28 type, public :: continuity_ppm_cs ; private
29  type(diag_ctrl), pointer :: diag !< Diagnostics control structure.
30  logical :: upwind_1st !< If true, use a first-order upwind scheme.
31  logical :: monotonic !< If true, use the Colella & Woodward monotonic
32  !! limiter; otherwise use a simple positive
33  !! definite limiter.
34  logical :: simple_2nd !< If true, use a simple second order (arithmetic
35  !! mean) interpolation of the edge values instead
36  !! of the higher order interpolation.
37  real :: tol_eta !< The tolerance for free-surface height
38  !! discrepancies between the barotropic solution and
39  !! the sum of the layer thicknesses [H ~> m or kg m-2].
40  real :: tol_vel !< The tolerance for barotropic velocity
41  !! discrepancies between the barotropic solution and
42  !! the sum of the layer thicknesses [m s-1].
43  real :: tol_eta_aux !< The tolerance for free-surface height
44  !! discrepancies between the barotropic solution and
45  !! the sum of the layer thicknesses when calculating
46  !! the auxiliary corrected velocities [H ~> m or kg m-2].
47  real :: cfl_limit_adjust !< The maximum CFL of the adjusted velocities [nondim]
48  logical :: aggress_adjust !< If true, allow the adjusted velocities to have a
49  !! relative CFL change up to 0.5. False by default.
50  logical :: vol_cfl !< If true, use the ratio of the open face lengths
51  !! to the tracer cell areas when estimating CFL
52  !! numbers. Without aggress_adjust, the default is
53  !! false; it is always true with.
54  logical :: better_iter !< If true, stop corrective iterations using a
55  !! velocity-based criterion and only stop if the
56  !! iteration is better than all predecessors.
57  logical :: use_visc_rem_max !< If true, use more appropriate limiting bounds
58  !! for corrections in strongly viscous columns.
59  logical :: marginal_faces !< If true, use the marginal face areas from the
60  !! continuity solver for use as the weights in the
61  !! barotropic solver. Otherwise use the transport
62  !! averaged areas.
63 end type continuity_ppm_cs
64 
65 !> A container for loop bounds
66 type :: loop_bounds_type ; private
67  !>@{ Loop bounds
68  integer :: ish, ieh, jsh, jeh
69  !!@}
70 end type loop_bounds_type
71 
72 contains
73 
74 !> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme,
75 !! based on Lin (1994).
76 subroutine continuity_ppm(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, &
77  visc_rem_u, visc_rem_v, u_cor, v_cor, &
78  uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont)
79  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
80  type(continuity_ppm_cs), pointer :: cs !< Module's control structure.
81  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
82  intent(in) :: u !< Zonal velocity [m s-1].
83  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
84  intent(in) :: v !< Meridional velocity [m s-1].
85  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
86  intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].
87  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
88  intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2].
89  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
90  intent(out) :: uh !< Zonal volume flux, u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
91  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
92  intent(out) :: vh !< Meridional volume flux, v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
93  real, intent(in) :: dt !< Time increment [s].
94  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
95  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
96  real, dimension(SZIB_(G),SZJ_(G)), &
97  optional, intent(in) :: uhbt !< The summed volume flux through zonal faces
98  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
99  real, dimension(SZI_(G),SZJB_(G)), &
100  optional, intent(in) :: vhbt !< The summed volume flux through meridional faces
101  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
102  type(ocean_obc_type), &
103  optional, pointer :: obc !< Open boundaries control structure.
104  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
105  optional, intent(in) :: visc_rem_u
106  !< The fraction of zonal momentum originally
107  !! in a layer that remains after a time-step of viscosity, and the
108  !! fraction of a time-step's worth of a barotropic acceleration that
109  !! a layer experiences after viscosity is applied.
110  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
111  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
112  optional, intent(in) :: visc_rem_v
113  !< The fraction of meridional momentum originally
114  !! in a layer that remains after a time-step of viscosity, and the
115  !! fraction of a time-step's worth of a barotropic acceleration that
116  !! a layer experiences after viscosity is applied.
117  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
118  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
119  optional, intent(out) :: u_cor
120  !< The zonal velocities that give uhbt as the depth-integrated transport [m s-1].
121  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
122  optional, intent(out) :: v_cor
123  !< The meridional velocities that give vhbt as the depth-integrated transport [m s-1].
124  real, dimension(SZIB_(G),SZJ_(G)), &
125  optional, intent(in) :: uhbt_aux
126  !< A second set of summed volume fluxes through zonal faces
127  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
128  real, dimension(SZI_(G),SZJB_(G)), &
129  optional, intent(in) :: vhbt_aux
130  !< A second set of summed volume fluxes through meridional faces
131  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
132  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
133  optional, intent(out) :: u_cor_aux
134  !< The zonal velocities that give uhbt_aux as the depth-integrated
135  !! transports [m s-1].
136  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
137  optional, intent(out) :: v_cor_aux
138  !< The meridional velocities that give vhbt_aux as the depth-integrated
139  !! transports [m s-1].
140  type(bt_cont_type), optional, pointer :: bt_cont !< A structure with elements that describe
141  !! the effective open face areas as a function of barotropic flow.
142 
143  ! Local variables
144  real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
145  type(loop_bounds_type) :: lb
146  integer :: is, ie, js, je, nz, stencil
147  integer :: i, j, k
148 
149  logical :: x_first
150  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
151 
152  h_min = gv%Angstrom_H
153 
154  if (.not.associated(cs)) call mom_error(fatal, &
155  "MOM_continuity_PPM: Module must be initialized before it is used.")
156  x_first = (mod(g%first_direction,2) == 0)
157 
158  if (present(visc_rem_u) .neqv. present(visc_rem_v)) call mom_error(fatal, &
159  "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// &
160  " one must be present in call to continuity_PPM.")
161 
162  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1
163 
164  if (x_first) then
165  ! First, advect zonally.
166  lb%ish = g%isc ; lb%ieh = g%iec
167  lb%jsh = g%jsc-stencil ; lb%jeh = g%jec+stencil
168  call zonal_mass_flux(u, hin, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, &
169  u_cor, uhbt_aux, u_cor_aux, bt_cont)
170 
171  call cpu_clock_begin(id_clock_update)
172  !$OMP parallel do default(shared)
173  do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
174  h(i,j,k) = hin(i,j,k) - dt* g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
175  ! Uncomment this line to prevent underflow.
176  ! if (h(i,j,k) < h_min) h(i,j,k) = h_min
177  enddo ; enddo ; enddo
178  call cpu_clock_end(id_clock_update)
179 
180  lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
181 
182  ! Now advect meridionally, using the updated thicknesses to determine
183  ! the fluxes.
184  call meridional_mass_flux(v, h, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, &
185  v_cor, vhbt_aux, v_cor_aux, bt_cont)
186 
187  call cpu_clock_begin(id_clock_update)
188  !$OMP parallel do default(shared)
189  do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
190  h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
191  ! This line prevents underflow.
192  if (h(i,j,k) < h_min) h(i,j,k) = h_min
193  enddo ; enddo ; enddo
194  call cpu_clock_end(id_clock_update)
195 
196  else ! .not. x_first
197  ! First, advect meridionally, so set the loop bounds accordingly.
198  lb%ish = g%isc-stencil ; lb%ieh = g%iec+stencil
199  lb%jsh = g%jsc ; lb%jeh = g%jec
200 
201  call meridional_mass_flux(v, hin, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, &
202  v_cor, vhbt_aux, v_cor_aux, bt_cont)
203 
204  call cpu_clock_begin(id_clock_update)
205  !$OMP parallel do default(shared)
206  do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
207  h(i,j,k) = hin(i,j,k) - dt*g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
208  enddo ; enddo ; enddo
209  call cpu_clock_end(id_clock_update)
210 
211  ! Now advect zonally, using the updated thicknesses to determine
212  ! the fluxes.
213  lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
214  call zonal_mass_flux(u, h, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, &
215  u_cor, uhbt_aux, u_cor_aux, bt_cont)
216 
217  call cpu_clock_begin(id_clock_update)
218  !$OMP parallel do default(shared)
219  do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
220  h(i,j,k) = h(i,j,k) - dt* g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
221  ! This line prevents underflow.
222  if (h(i,j,k) < h_min) h(i,j,k) = h_min
223  enddo ; enddo ; enddo
224  call cpu_clock_end(id_clock_update)
225 
226  endif
227 
228 end subroutine continuity_ppm
229 
230 !> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.
231 subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, &
232  visc_rem_u, u_cor, uhbt_aux, u_cor_aux, BT_cont)
233  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
234  type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
235  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
236  intent(in) :: u !< Zonal velocity [m s-1].
237  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
238  intent(in) :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
239  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
240  intent(out) :: uh !< Volume flux through zonal faces = u*h*dy
241  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
242  real, intent(in) :: dt !< Time increment [s].
243  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
244  type(continuity_ppm_cs), pointer :: CS !< This module's control structure.
245  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
246  type(ocean_obc_type), &
247  optional, pointer :: OBC !< Open boundaries control structure.
248  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
249  optional, intent(in) :: visc_rem_u
250  !< The fraction of zonal momentum originally in a layer that remains after a
251  !! time-step of viscosity, and the fraction of a time-step's worth of a barotropic
252  !! acceleration that a layer experiences after viscosity is applied.
253  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
254  real, dimension(SZIB_(G),SZJ_(G)), &
255  optional, intent(in) :: uhbt !< The summed volume flux through zonal faces
256  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
257  real, dimension(SZIB_(G),SZJ_(G)), &
258  optional, intent(in) :: uhbt_aux
259  !< A second set of summed volume fluxes through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
260  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
261  optional, intent(out) :: u_cor
262  !< The zonal velocitiess (u with a barotropic correction)
263  !! that give uhbt as the depth-integrated transport, m s-1.
264  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
265  optional, intent(out) :: u_cor_aux
266  !< The zonal velocities (u with a barotropic correction)
267  !! that give uhbt_aux as the depth-integrated transports [m s-1].
268  type(bt_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe the
269  !! effective open face areas as a function of barotropic flow.
270 
271  ! Local variables
272  real, dimension(SZIB_(G),SZK_(G)) :: duhdu ! Partial derivative of uh with u [H m ~> m2 or kg m-1].
273  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_L, h_R ! Left and right face thicknesses [H ~> m or kg m-2].
274  real, dimension(SZIB_(G)) :: &
275  du, & ! Corrective barotropic change in the velocity [m s-1].
276  du_min_CFL, & ! Min/max limits on du correction
277  du_max_CFL, & ! to avoid CFL violations
278  duhdu_tot_0, & ! Summed partial derivative of uh with u [H m ~> m2 or kg m-1].
279  uh_tot_0, & ! Summed transport with no barotropic correction [H m2 s-1 ~> m3 s-1 or kg s-1].
280  visc_rem_max ! The column maximum of visc_rem.
281  logical, dimension(SZIB_(G)) :: do_I
282  real, dimension(SZIB_(G),SZK_(G)) :: &
283  visc_rem ! A 2-D copy of visc_rem_u or an array of 1's.
284  real, dimension(SZIB_(G)) :: FAuI ! A list of sums of zonal face areas [H m ~> m2 or kg m-1].
285  real :: FA_u ! A sum of zonal face areas [H m ~> m2 or kg m-1].
286  real :: I_vrm ! 1.0 / visc_rem_max, nondim.
287  real :: CFL_dt ! The maximum CFL ratio of the adjusted velocities divided by
288  ! the time step [s-1].
289  real :: I_dt ! 1.0 / dt [s-1].
290  real :: du_lim ! The velocity change that give a relative CFL of 1 [m s-1].
291  real :: dx_E, dx_W ! Effective x-grid spacings to the east and west [m].
292  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
293  logical :: do_aux, local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
294  logical :: local_Flather_OBC, local_open_BC, is_simple
295  type(obc_segment_type), pointer :: segment => null()
296 
297  do_aux = (present(uhbt_aux) .and. present(u_cor_aux))
298  use_visc_rem = present(visc_rem_u)
299  local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
300  local_open_bc = .false.
301  if (present(bt_cont)) set_bt_cont = (associated(bt_cont))
302  if (present(obc)) then ; if (associated(obc)) then
303  local_specified_bc = obc%specified_u_BCs_exist_globally
304  local_flather_obc = obc%Flather_u_BCs_exist_globally
305  local_open_bc = obc%open_u_BCs_exist_globally
306  endif ; endif
307  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
308 
309  cfl_dt = cs%CFL_limit_adjust / dt
310  i_dt = 1.0 / dt
311  if (cs%aggress_adjust) cfl_dt = i_dt
312 
313  call cpu_clock_begin(id_clock_update)
314 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,CS,h_L,h_in,h_R,G,GV,LB,visc_rem,OBC)
315  do k=1,nz
316  ! This sets h_L and h_R.
317  if (cs%upwind_1st) then
318  do j=jsh,jeh ; do i=ish-1,ieh+1
319  h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
320  enddo ; enddo
321  else
322  call ppm_reconstruction_x(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
323  2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
324  endif
325  do i=ish-1,ieh ; visc_rem(i,k) = 1.0 ; enddo
326  enddo
327  call cpu_clock_end(id_clock_update)
328 
329  call cpu_clock_begin(id_clock_correct)
330 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,u,h_in,h_L,h_R,use_visc_rem,visc_rem_u, &
331 !$OMP uh,dt,G,GV,CS,local_specified_BC,OBC,uhbt,do_aux,set_BT_cont, &
332 !$OMP CFL_dt,I_dt,u_cor,uhbt_aux,u_cor_aux,BT_cont, local_Flather_OBC) &
333 !$OMP private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0,duhdu_tot_0, &
334 !$OMP is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W,any_simple_OBC ) &
335 !$OMP firstprivate(visc_rem)
336  do j=jsh,jeh
337  do i=ish-1,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ; enddo
338  ! Set uh and duhdu.
339  do k=1,nz
340  if (use_visc_rem) then ; do i=ish-1,ieh
341  visc_rem(i,k) = visc_rem_u(i,j,k)
342  visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
343  enddo ; endif
344  call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
345  uh(:,j,k), duhdu(:,k), visc_rem(:,k), &
346  dt, g, j, ish, ieh, do_i, cs%vol_CFL, obc)
347  if (local_specified_bc) then
348  do i=ish-1,ieh
349  if (obc%segment(obc%segnum_u(i,j))%specified) &
350  uh(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k)
351  enddo
352  endif
353  enddo
354 
355  if ((.not.use_visc_rem).or.(.not.cs%use_visc_rem_max)) then ; do i=ish-1,ieh
356  visc_rem_max(i) = 1.0
357  enddo ; endif
358 
359  if (present(uhbt) .or. do_aux .or. set_bt_cont) then
360  ! Set limits on du that will keep the CFL number between -1 and 1.
361  ! This should be adequate to keep the root bracketed in all cases.
362  do i=ish-1,ieh
363  i_vrm = 0.0
364  if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
365  if (cs%vol_CFL) then
366  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
367  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
368  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
369  du_max_cfl(i) = 2.0* (cfl_dt * dx_w) * i_vrm
370  du_min_cfl(i) = -2.0 * (cfl_dt * dx_e) * i_vrm
371  uh_tot_0(i) = 0.0 ; duhdu_tot_0(i) = 0.0
372  enddo
373  do k=1,nz ; do i=ish-1,ieh
374  duhdu_tot_0(i) = duhdu_tot_0(i) + duhdu(i,k)
375  uh_tot_0(i) = uh_tot_0(i) + uh(i,j,k)
376  enddo ; enddo
377  if (use_visc_rem) then
378  if (cs%aggress_adjust) then
379  do k=1,nz ; do i=ish-1,ieh
380  if (cs%vol_CFL) then
381  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
382  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
383  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
384 
385  du_lim = 0.499*((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k)))
386  if (du_max_cfl(i) * visc_rem(i,k) > du_lim) &
387  du_max_cfl(i) = du_lim / visc_rem(i,k)
388 
389  du_lim = 0.499*((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k)))
390  if (du_min_cfl(i) * visc_rem(i,k) < du_lim) &
391  du_min_cfl(i) = du_lim / visc_rem(i,k)
392  enddo ; enddo
393  else
394  do k=1,nz ; do i=ish-1,ieh
395  if (cs%vol_CFL) then
396  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
397  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
398  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
399 
400  if (du_max_cfl(i) * visc_rem(i,k) > dx_w*cfl_dt - u(i,j,k)) &
401  du_max_cfl(i) = (dx_w*cfl_dt - u(i,j,k)) / visc_rem(i,k)
402  if (du_min_cfl(i) * visc_rem(i,k) < -dx_e*cfl_dt - u(i,j,k)) &
403  du_min_cfl(i) = -(dx_e*cfl_dt + u(i,j,k)) / visc_rem(i,k)
404  enddo ; enddo
405  endif
406  else
407  if (cs%aggress_adjust) then
408  do k=1,nz ; do i=ish-1,ieh
409  if (cs%vol_CFL) then
410  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
411  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
412  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
413 
414  du_max_cfl(i) = min(du_max_cfl(i), 0.499 * &
415  ((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k))) )
416  du_min_cfl(i) = max(du_min_cfl(i), 0.499 * &
417  ((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k))) )
418  enddo ; enddo
419  else
420  do k=1,nz ; do i=ish-1,ieh
421  if (cs%vol_CFL) then
422  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
423  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
424  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
425 
426  du_max_cfl(i) = min(du_max_cfl(i), dx_w*cfl_dt - u(i,j,k))
427  du_min_cfl(i) = max(du_min_cfl(i), -(dx_e*cfl_dt + u(i,j,k)))
428  enddo ; enddo
429  endif
430  endif
431  do i=ish-1,ieh
432  du_max_cfl(i) = max(du_max_cfl(i),0.0)
433  du_min_cfl(i) = min(du_min_cfl(i),0.0)
434  enddo
435 
436  ! Up to this point, everything is shared between uhbt and uhbt_aux.
437 
438  any_simple_obc = .false.
439  if (present(uhbt) .or. do_aux .or. set_bt_cont) then
440  if (local_specified_bc .or. local_flather_obc) then ; do i=ish-1,ieh
441  ! Avoid reconciling barotropic/baroclinic transports if transport is specified
442  is_simple = obc%segment(obc%segnum_u(i,j))%specified
443  do_i(i) = .not.(obc%segnum_u(i,j) /= obc_none .and. is_simple)
444  any_simple_obc = any_simple_obc .or. is_simple
445  enddo ; else ; do i=ish-1,ieh
446  do_i(i) = .true.
447  enddo ; endif
448  endif
449 
450  if (present(uhbt)) then
451  call zonal_flux_adjust(u, h_in, h_l, h_r, uhbt(:,j), uh_tot_0, &
452  duhdu_tot_0, du, du_max_cfl, du_min_cfl, dt, g, &
453  cs, visc_rem, j, ish, ieh, do_i, .true., uh, obc=obc)
454 
455  if (present(u_cor)) then ; do k=1,nz
456  do i=ish-1,ieh ; u_cor(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
457  if (local_specified_bc) then ; do i=ish-1,ieh
458  if (obc%segment(obc%segnum_u(i,j))%specified) &
459  u_cor(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
460  enddo ; endif
461  enddo ; endif ! u-corrected
462 
463  endif
464 
465  if (do_aux) then
466  call zonal_flux_adjust(u, h_in, h_l, h_r, uhbt_aux(:,j), uh_tot_0, &
467  duhdu_tot_0, du, du_max_cfl, du_min_cfl, dt, g, &
468  cs, visc_rem, j, ish, ieh, do_i, .false., obc=obc)
469 
470  do k=1,nz
471  do i=ish-1,ieh ; u_cor_aux(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
472  if (local_specified_bc) then ; do i=ish-1,ieh
473  if (obc%segment(obc%segnum_u(i,j))%specified) &
474  u_cor_aux(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
475  enddo ; endif
476  enddo
477  endif ! do_aux
478 
479  if (set_bt_cont) then
480  call set_zonal_bt_cont(u, h_in, h_l, h_r, bt_cont, uh_tot_0, duhdu_tot_0,&
481  du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
482  visc_rem_max, j, ish, ieh, do_i)
483  if (any_simple_obc) then
484  do i=ish-1,ieh
485  do_i(i) = obc%segment(obc%segnum_u(i,j))%specified
486  if (do_i(i)) faui(i) = gv%H_subroundoff*g%dy_Cu(i,j)
487  enddo
488  do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
489  if ((abs(obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
490  (obc%segment(obc%segnum_u(i,j))%specified)) &
491  faui(i) = faui(i) + obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k) / &
492  obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
493  endif ; enddo ; enddo
494  do i=ish-1,ieh ; if (do_i(i)) then
495  bt_cont%FA_u_W0(i,j) = us%m_to_L*faui(i) ; bt_cont%FA_u_E0(i,j) = us%m_to_L*faui(i)
496  bt_cont%FA_u_WW(i,j) = us%m_to_L*faui(i) ; bt_cont%FA_u_EE(i,j) = us%m_to_L*faui(i)
497  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
498  endif ; enddo
499  endif
500  endif ! set_BT_cont
501 
502  endif ! present(uhbt) or do_aux or set_BT_cont
503  enddo ! j-loop
504  if (local_open_bc .and. set_bt_cont) then
505  do n = 1, obc%number_of_segments
506  if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
507  i = obc%segment(n)%HI%IsdB
508  if (obc%segment(n)%direction == obc_direction_e) then
509  do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
510  fa_u = 0.0
511  do k=1,nz ; fa_u = fa_u + h_in(i,j,k)*us%m_to_L*g%dy_Cu(i,j) ; enddo
512  bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
513  bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
514  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
515  enddo
516  else
517  do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
518  fa_u = 0.0
519  do k=1,nz ; fa_u = fa_u + h_in(i+1,j,k)*us%m_to_L*g%dy_Cu(i,j) ; enddo
520  bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
521  bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
522  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
523  enddo
524  endif
525  endif
526  enddo
527  endif
528  call cpu_clock_end(id_clock_correct)
529 
530  if (set_bt_cont) then ; if (allocated(bt_cont%h_u)) then
531  if (present(u_cor)) then
532  call zonal_face_thickness(u_cor, h_in, h_l, h_r, bt_cont%h_u, dt, g, lb, &
533  cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
534  else
535  call zonal_face_thickness(u, h_in, h_l, h_r, bt_cont%h_u, dt, g, lb, &
536  cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
537  endif
538  endif ; endif
539 
540 end subroutine zonal_mass_flux
541 
542 !> Evaluates the zonal mass or volume fluxes in a layer.
543 subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, j, &
544  ish, ieh, do_I, vol_CFL, OBC)
545  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
546  real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity [m s-1].
547  real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the
548  !! momentum originally in a layer that remains after a time-step
549  !! of viscosity, and the fraction of a time-step's worth of a barotropic
550  !! acceleration that a layer experiences after viscosity is applied.
551  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
552  real, dimension(SZI_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
553  real, dimension(SZI_(G)), intent(in) :: h_L !< Left thickness [H ~> m or kg m-2].
554  real, dimension(SZI_(G)), intent(in) :: h_R !< Right thickness [H ~> m or kg m-2].
555  real, dimension(SZIB_(G)), intent(inout) :: uh !< Zonal mass or volume
556  !! transport [H m2 s-1 ~> m3 s-1 or kg s-1].
557  real, dimension(SZIB_(G)), intent(inout) :: duhdu !< Partial derivative of uh
558  !! with u [H m ~> m2 or kg m-1].
559  real, intent(in) :: dt !< Time increment [s].
560  integer, intent(in) :: j !< Spatial index.
561  integer, intent(in) :: ish !< Start of index range.
562  integer, intent(in) :: ieh !< End of index range.
563  logical, dimension(SZIB_(G)), intent(in) :: do_I !< Which i values to work on.
564  logical, intent(in) :: vol_CFL !< If true, rescale the
565  !! ratio of face areas to the cell areas when estimating the CFL number.
566  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
567  ! Local variables
568  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
569  real :: curv_3 ! A measure of the thickness curvature over a grid length,
570  ! with the same units as h_in.
571  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
572  integer :: i
573  logical :: local_open_BC
574 
575  local_open_bc = .false.
576  if (present(obc)) then ; if (associated(obc)) then
577  local_open_bc = obc%open_u_BCs_exist_globally
578  endif ; endif
579 
580  do i=ish-1,ieh ; if (do_i(i)) then
581  ! Set new values of uh and duhdu.
582  if (u(i) > 0.0) then
583  if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
584  else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
585  curv_3 = h_l(i) + h_r(i) - 2.0*h(i)
586  uh(i) = g%dy_Cu(i,j) * u(i) * &
587  (h_r(i) + cfl * (0.5*(h_l(i) - h_r(i)) + curv_3*(cfl - 1.5)))
588  h_marg = h_r(i) + cfl * ((h_l(i) - h_r(i)) + 3.0*curv_3*(cfl - 1.0))
589  elseif (u(i) < 0.0) then
590  if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
591  else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
592  curv_3 = h_l(i+1) + h_r(i+1) - 2.0*h(i+1)
593  uh(i) = g%dy_Cu(i,j) * u(i) * &
594  (h_l(i+1) + cfl * (0.5*(h_r(i+1)-h_l(i+1)) + curv_3*(cfl - 1.5)))
595  h_marg = h_l(i+1) + cfl * ((h_r(i+1)-h_l(i+1)) + 3.0*curv_3*(cfl - 1.0))
596  else
597  uh(i) = 0.0
598  h_marg = 0.5 * (h_l(i+1) + h_r(i))
599  endif
600  duhdu(i) = g%dy_Cu(i,j) * h_marg * visc_rem(i)
601  endif ; enddo
602 
603  if (local_open_bc) then
604  do i=ish-1,ieh ; if (do_i(i)) then
605  if (obc%segment(obc%segnum_u(i,j))%open) then
606  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
607  uh(i) = g%dy_Cu(i,j) * u(i) * h(i)
608  duhdu(i) = g%dy_Cu(i,j) * h(i) * visc_rem(i)
609  else
610  uh(i) = g%dy_Cu(i,j) * u(i) * h(i+1)
611  duhdu(i) = g%dy_Cu(i,j) * h(i+1) * visc_rem(i)
612  endif
613  endif
614  endif ; enddo
615  endif
616 end subroutine zonal_flux_layer
617 
618 !> Sets the effective interface thickness at each zonal velocity point.
619 subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, LB, vol_CFL, &
620  marginal, visc_rem_u, OBC)
621  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
622  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [m s-1].
623  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness used to
624  !! calculate fluxes [H ~> m or kg m-2].
625  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the
626  !! reconstruction [H ~> m or kg m-2].
627  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the
628  !! reconstruction [H ~> m or kg m-2].
629  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_u !< Thickness at zonal faces [H ~> m or kg m-2].
630  real, intent(in) :: dt !< Time increment [s].
631  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
632  logical, intent(in) :: vol_CFL !< If true, rescale the ratio
633  !! of face areas to the cell areas when estimating the CFL number.
634  logical, intent(in) :: marginal !< If true, report the
635  !! marginal face thicknesses; otherwise report transport-averaged thicknesses.
636  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
637  optional, intent(in) :: visc_rem_u
638  !< Both the fraction of the momentum originally in a layer that remains after
639  !! a time-step of viscosity, and the fraction of a time-step's worth of a
640  !! barotropic acceleration that a layer experiences after viscosity is applied.
641  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
642  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
643 
644  ! Local variables
645  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
646  real :: curv_3 ! A measure of the thickness curvature over a grid length,
647  ! with the same units as h_in.
648  real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].
649  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
650  logical :: local_open_BC
651  integer :: i, j, k, ish, ieh, jsh, jeh, nz, n
652  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
653 
654  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
655  do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
656  if (u(i,j,k) > 0.0) then
657  if (vol_cfl) then ; cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
658  else ; cfl = u(i,j,k) * dt * g%IdxT(i,j) ; endif
659  curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
660  h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
661  h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + 3.0*curv_3*(cfl - 1.0))
662  elseif (u(i,j,k) < 0.0) then
663  if (vol_cfl) then ; cfl = (-u(i,j,k)*dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
664  else ; cfl = -u(i,j,k) * dt * g%IdxT(i+1,j) ; endif
665  curv_3 = h_l(i+1,j,k) + h_r(i+1,j,k) - 2.0*h(i+1,j,k)
666  h_avg = h_l(i+1,j,k) + cfl * (0.5*(h_r(i+1,j,k)-h_l(i+1,j,k)) + curv_3*(cfl - 1.5))
667  h_marg = h_l(i+1,j,k) + cfl * ((h_r(i+1,j,k)-h_l(i+1,j,k)) + &
668  3.0*curv_3*(cfl - 1.0))
669  else
670  h_avg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
671  ! The choice to use the arithmetic mean here is somewhat arbitrariy, but
672  ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
673  h_marg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
674  ! h_marg = (2.0 * h_L(i+1,j,k) * h_R(i,j,k)) / &
675  ! (h_L(i+1,j,k) + h_R(i,j,k) + GV%H_subroundoff)
676  endif
677 
678  if (marginal) then ; h_u(i,j,k) = h_marg
679  else ; h_u(i,j,k) = h_avg ; endif
680  enddo ; enddo ; enddo
681  if (present(visc_rem_u)) then
682  !$OMP parallel do default(shared)
683  do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
684  h_u(i,j,k) = h_u(i,j,k) * visc_rem_u(i,j,k)
685  enddo ; enddo ; enddo
686  endif
687 
688  local_open_bc = .false.
689  if (present(obc)) then ; if (associated(obc)) then
690  local_open_bc = obc%open_u_BCs_exist_globally
691  endif ; endif
692  if (local_open_bc) then
693  do n = 1, obc%number_of_segments
694  if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
695  i = obc%segment(n)%HI%IsdB
696  if (obc%segment(n)%direction == obc_direction_e) then
697  if (present(visc_rem_u)) then ; do k=1,nz
698  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
699  h_u(i,j,k) = h(i,j,k) * visc_rem_u(i,j,k)
700  enddo
701  enddo ; else ; do k=1,nz
702  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
703  h_u(i,j,k) = h(i,j,k)
704  enddo
705  enddo ; endif
706  else
707  if (present(visc_rem_u)) then ; do k=1,nz
708  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
709  h_u(i,j,k) = h(i+1,j,k) * visc_rem_u(i,j,k)
710  enddo
711  enddo ; else ; do k=1,nz
712  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
713  h_u(i,j,k) = h(i+1,j,k)
714  enddo
715  enddo ; endif
716  endif
717  endif
718  enddo
719  endif
720 
721 end subroutine zonal_face_thickness
722 
723 !> Returns the barotropic velocity adjustment that gives the
724 !! desired barotropic (layer-summed) transport.
725 subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, &
726  du, du_max_CFL, du_min_CFL, dt, G, CS, visc_rem, &
727  j, ish, ieh, do_I_in, full_precision, uh_3d, OBC)
728  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
729  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [m s-1].
730  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to
731  !! calculate fluxes [H ~> m or kg m-2].
732  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the
733  !! reconstruction [H ~> m or kg m-2].
734  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the
735  !! reconstruction [H ~> m or kg m-2].
736  real, dimension(SZIB_(G),SZK_(G)), intent(in) :: visc_rem !< Both the fraction of the
737  !! momentum originally in a layer that remains after a time-step of viscosity, and
738  !! the fraction of a time-step's worth of a barotropic acceleration that a layer
739  !! experiences after viscosity is applied.
740  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
741  real, dimension(SZIB_(G)), optional, intent(in) :: uhbt !< The summed volume flux
742  !! through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
743 
744  real, dimension(SZIB_(G)), intent(in) :: du_max_CFL !< Maximum acceptable
745  !! value of du [m s-1].
746  real, dimension(SZIB_(G)), intent(in) :: du_min_CFL !< Minimum acceptable
747  !! value of du [m s-1].
748  real, dimension(SZIB_(G)), intent(in) :: uh_tot_0 !< The summed transport
749  !! with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
750  real, dimension(SZIB_(G)), intent(in) :: duhdu_tot_0 !< The partial derivative
751  !! of du_err with du at 0 adjustment [H m ~> m2 or kg m-1].
752  real, dimension(SZIB_(G)), intent(out) :: du !<
753  !! The barotropic velocity adjustment [m s-1].
754  real, intent(in) :: dt !< Time increment [s].
755  type(continuity_ppm_cs), pointer :: CS !< This module's control structure.
756  integer, intent(in) :: j !< Spatial index.
757  integer, intent(in) :: ish !< Start of index range.
758  integer, intent(in) :: ieh !< End of index range.
759  logical, dimension(SZIB_(G)), intent(in) :: do_I_in !<
760  !! A logical flag indicating which I values to work on.
761  logical, optional, intent(in) :: full_precision !<
762  !! A flag indicating how carefully to iterate. The
763  !! default is .true. (more accurate).
764  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(inout) :: uh_3d !<
765  !! Volume flux through zonal faces = u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
766  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
767  ! Local variables
768  real, dimension(SZIB_(G),SZK_(G)) :: &
769  uh_aux, & ! An auxiliary zonal volume flux [H m s-1 ~> m2 s-1 or kg m-1 s-1].
770  duhdu ! Partial derivative of uh with u [H m ~> m2 or kg m-1].
771  real, dimension(SZIB_(G)) :: &
772  uh_err, & ! Difference between uhbt and the summed uh [H m2 s-1 ~> m3 s-1 or kg s-1].
773  uh_err_best, & ! The smallest value of uh_err found so far [H m2 s-1 ~> m3 s-1 or kg s-1].
774  u_new, & ! The velocity with the correction added [m s-1].
775  duhdu_tot,&! Summed partial derivative of uh with u [H m ~> m2 or kg m-1].
776  du_min, & ! Min/max limits on du correction based on CFL limits
777  du_max ! and previous iterations [m s-1].
778  real :: du_prev ! The previous value of du [m s-1].
779  real :: ddu ! The change in du from the previous iteration [m s-1].
780  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
781  real :: tol_vel ! The tolerance for velocity in the current iteration [m s-1].
782  integer :: i, k, nz, itt, max_itts = 20
783  logical :: full_prec, domore, do_I(SZIB_(G))
784 
785  nz = g%ke
786  full_prec = .true. ; if (present(full_precision)) full_prec = full_precision
787 
788  uh_aux(:,:) = 0.0 ; duhdu(:,:) = 0.0
789 
790  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
791  uh_aux(i,k) = uh_3d(i,j,k)
792  enddo ; enddo ; endif
793 
794  do i=ish-1,ieh
795  du(i) = 0.0 ; do_i(i) = do_i_in(i)
796  du_max(i) = du_max_cfl(i) ; du_min(i) = du_min_cfl(i)
797  uh_err(i) = uh_tot_0(i) - uhbt(i) ; duhdu_tot(i) = duhdu_tot_0(i)
798  uh_err_best(i) = abs(uh_err(i))
799  enddo
800 
801  do itt=1,max_itts
802  if (full_prec) then
803  select case (itt)
804  case (:1) ; tol_eta = 1e-6 * cs%tol_eta
805  case (2) ; tol_eta = 1e-4 * cs%tol_eta
806  case (3) ; tol_eta = 1e-2 * cs%tol_eta
807  case default ; tol_eta = cs%tol_eta
808  end select
809  else
810  tol_eta = cs%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
811  endif
812  tol_vel = cs%tol_vel
813 
814  do i=ish-1,ieh
815  if (uh_err(i) > 0.0) then ; du_max(i) = du(i)
816  elseif (uh_err(i) < 0.0) then ; du_min(i) = du(i)
817  else ; do_i(i) = .false. ; endif
818  enddo
819  domore = .false.
820  do i=ish-1,ieh ; if (do_i(i)) then
821  if ((dt*min(g%IareaT(i,j),g%IareaT(i+1,j))*abs(uh_err(i)) > tol_eta) .or.&
822  (cs%better_iter .and. ((abs(uh_err(i)) > tol_vel * duhdu_tot(i)) .or.&
823  (abs(uh_err(i)) > uh_err_best(i))) )) then
824  ! Use Newton's method, provided it stays bounded. Otherwise bisect
825  ! the value with the appropriate bound.
826  if (full_prec) then
827  ddu = -uh_err(i) / duhdu_tot(i)
828  du_prev = du(i)
829  du(i) = du(i) + ddu
830  if (abs(ddu) < 1.0e-15*abs(du(i))) then
831  do_i(i) = .false. ! ddu is small enough to quit.
832  elseif (ddu > 0.0) then
833  if (du(i) >= du_max(i)) then
834  du(i) = 0.5*(du_prev + du_max(i))
835  if (du_max(i) - du_prev < 1.0e-15*abs(du(i))) do_i(i) = .false.
836  endif
837  else ! ddu < 0.0
838  if (du(i) <= du_min(i)) then
839  du(i) = 0.5*(du_prev + du_min(i))
840  if (du_prev - du_min(i) < 1.0e-15*abs(du(i))) do_i(i) = .false.
841  endif
842  endif
843  else
844  ! Use Newton's method, provided it stays bounded, just like above.
845  du(i) = du(i) - uh_err(i) / duhdu_tot(i)
846  if ((du(i) >= du_max(i)) .or. (du(i) <= du_min(i))) &
847  du(i) = 0.5*(du_max(i) + du_min(i))
848  endif
849  if (do_i(i)) domore = .true.
850  else
851  do_i(i) = .false.
852  endif
853  endif ; enddo
854  if (.not.domore) exit
855 
856  if ((itt < max_itts) .or. present(uh_3d)) then ; do k=1,nz
857  do i=ish-1,ieh ; u_new(i) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
858  call zonal_flux_layer(u_new, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
859  uh_aux(:,k), duhdu(:,k), visc_rem(:,k), &
860  dt, g, j, ish, ieh, do_i, cs%vol_CFL, obc)
861  enddo ; endif
862 
863  if (itt < max_itts) then
864  do i=ish-1,ieh
865  uh_err(i) = -uhbt(i) ; duhdu_tot(i) = 0.0
866  enddo
867  do k=1,nz ; do i=ish-1,ieh
868  uh_err(i) = uh_err(i) + uh_aux(i,k)
869  duhdu_tot(i) = duhdu_tot(i) + duhdu(i,k)
870  enddo ; enddo
871  do i=ish-1,ieh
872  uh_err_best(i) = min(uh_err_best(i), abs(uh_err(i)))
873  enddo
874  endif
875  enddo ! itt-loop
876  ! If there are any faces which have not converged to within the tolerance,
877  ! so-be-it, or else use a final upwind correction?
878  ! This never seems to happen with 20 iterations as max_itt.
879 
880  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
881  uh_3d(i,j,k) = uh_aux(i,k)
882  enddo ; enddo ; endif
883 
884 end subroutine zonal_flux_adjust
885 
886 !> Sets a structure that describes the zonal barotropic volume or mass fluxes as a
887 !! function of barotropic flow to agree closely with the sum of the layer's transports.
888 subroutine set_zonal_bt_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, &
889  du_max_CFL, du_min_CFL, dt, G, US, CS, visc_rem, &
890  visc_rem_max, j, ish, ieh, do_I)
891  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
892  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [m s-1].
893  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to
894  !! calculate fluxes [H ~> m or kg m-2].
895  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the
896  !! reconstruction [H ~> m or kg m-2].
897  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the
898  !! reconstruction [H ~> m or kg m-2].
899  type(bt_cont_type), intent(inout) :: BT_cont !< A structure with elements
900  !! that describe the effective open face areas as a function of barotropic flow.
901  real, dimension(SZIB_(G)), intent(in) :: uh_tot_0 !< The summed transport
902  !! with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
903  real, dimension(SZIB_(G)), intent(in) :: duhdu_tot_0 !< The partial derivative
904  !! of du_err with du at 0 adjustment [H m ~> m2 or kg m-1].
905  real, dimension(SZIB_(G)), intent(in) :: du_max_CFL !< Maximum acceptable
906  !! value of du [m s-1].
907  real, dimension(SZIB_(G)), intent(in) :: du_min_CFL !< Minimum acceptable
908  !! value of du [m s-1].
909  real, intent(in) :: dt !< Time increment [s].
910  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
911  type(continuity_ppm_cs), pointer :: CS !< This module's control structure.
912  real, dimension(SZIB_(G),SZK_(G)), intent(in) :: visc_rem !< Both the fraction of the
913  !! momentum originally in a layer that remains after a time-step of viscosity, and
914  !! the fraction of a time-step's worth of a barotropic acceleration that a layer
915  !! experiences after viscosity is applied.
916  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
917  real, dimension(SZIB_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem.
918  integer, intent(in) :: j !< Spatial index.
919  integer, intent(in) :: ish !< Start of index range.
920  integer, intent(in) :: ieh !< End of index range.
921  logical, dimension(SZIB_(G)), intent(in) :: do_I !< A logical flag indicating
922  !! which I values to work on.
923  ! Local variables
924  real, dimension(SZIB_(G)) :: &
925  du0, & ! The barotropic velocity increment that gives 0 transport [m s-1].
926  duL, duR, & ! The barotropic velocity increments that give the westerly
927  ! (duL) and easterly (duR) test velocities.
928  zeros, & ! An array of full of 0's.
929  du_cfl, & ! The velocity increment that corresponds to CFL_min [m s-1].
930  u_l, u_r, & ! The westerly (u_L), easterly (u_R), and zero-barotropic
931  u_0, & ! transport (u_0) layer test velocities [m s-1].
932  fa_marg_l, & ! The effective layer marginal face areas with the westerly
933  fa_marg_r, & ! (_L), easterly (_R), and zero-barotropic (_0) test
934  fa_marg_0, & ! velocities [H m ~> m2 or kg m-1].
935  uh_l, uh_r, & ! The layer transports with the westerly (_L), easterly (_R),
936  uh_0, & ! and zero-barotropic (_0) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
937  famt_l, famt_r, & ! The summed effective marginal face areas for the 3
938  famt_0, & ! test velocities [H m ~> m2 or kg m-1].
939  uhtot_l, & ! The summed transport with the westerly (uhtot_L) and
940  uhtot_r ! and easterly (uhtot_R) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
941  real :: FA_0 ! The effective face area with 0 barotropic transport [m H ~> m2 or kg m].
942  real :: FA_avg ! The average effective face area [m H ~> m2 or kg m], nominally given by
943  ! the realized transport divided by the barotropic velocity.
944  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim] This
945  ! limiting is necessary to keep the inverse of visc_rem
946  ! from leading to large CFL numbers.
947  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
948  ! in finding the barotropic velocity that changes the
949  ! flow direction. This is necessary to keep the inverse
950  ! of visc_rem from leading to large CFL numbers.
951  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
952  ! flow is truly upwind [nondim]
953  real :: Idt ! The inverse of the time step [s-1].
954  logical :: domore
955  integer :: i, k, nz
956 
957  nz = g%ke ; idt = 1.0/dt
958  min_visc_rem = 0.1 ; cfl_min = 1e-6
959 
960  ! Diagnose the zero-transport correction, du0.
961  do i=ish-1,ieh ; zeros(i) = 0.0 ; enddo
962  call zonal_flux_adjust(u, h_in, h_l, h_r, zeros, uh_tot_0, &
963  duhdu_tot_0, du0, du_max_cfl, du_min_cfl, dt, g, &
964  cs, visc_rem, j, ish, ieh, do_i, .true.)
965 
966  ! Determine the westerly- and easterly- fluxes. Choose a sufficiently
967  ! negative velocity correction for the easterly-flux, and a sufficiently
968  ! positive correction for the westerly-flux.
969  domore = .false.
970  do i=ish-1,ieh
971  if (do_i(i)) domore = .true.
972  du_cfl(i) = (cfl_min * idt) * g%dxCu(i,j)
973  dur(i) = min(0.0,du0(i) - du_cfl(i))
974  dul(i) = max(0.0,du0(i) + du_cfl(i))
975  famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
976  uhtot_l(i) = 0.0 ; uhtot_r(i) = 0.0
977  enddo
978 
979  if (.not.domore) then
980  do k=1,nz ; do i=ish-1,ieh
981  bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
982  bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
983  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
984  enddo ; enddo
985  return
986  endif
987 
988  do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
989  visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
990  if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
991  if (u(i,j,k) + dur(i)*visc_rem_lim > -du_cfl(i)*visc_rem(i,k)) &
992  dur(i) = -(u(i,j,k) + du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
993  if (u(i,j,k) + dul(i)*visc_rem_lim < du_cfl(i)*visc_rem(i,k)) &
994  dul(i) = -(u(i,j,k) - du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
995  endif
996  endif ; enddo ; enddo
997 
998  do k=1,nz
999  do i=ish-1,ieh ; if (do_i(i)) then
1000  u_l(i) = u(i,j,k) + dul(i) * visc_rem(i,k)
1001  u_r(i) = u(i,j,k) + dur(i) * visc_rem(i,k)
1002  u_0(i) = u(i,j,k) + du0(i) * visc_rem(i,k)
1003  endif ; enddo
1004  call zonal_flux_layer(u_0, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_0, &
1005  fa_marg_0, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
1006  cs%vol_CFL)
1007  call zonal_flux_layer(u_l, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_l, &
1008  fa_marg_l, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
1009  cs%vol_CFL)
1010  call zonal_flux_layer(u_r, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_r, &
1011  fa_marg_r, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
1012  cs%vol_CFL)
1013  do i=ish-1,ieh ; if (do_i(i)) then
1014  famt_0(i) = famt_0(i) + fa_marg_0(i)
1015  famt_l(i) = famt_l(i) + fa_marg_l(i)
1016  famt_r(i) = famt_r(i) + fa_marg_r(i)
1017  uhtot_l(i) = uhtot_l(i) + uh_l(i)
1018  uhtot_r(i) = uhtot_r(i) + uh_r(i)
1019  endif ; enddo
1020  enddo
1021  do i=ish-1,ieh ; if (do_i(i)) then
1022  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1023  if ((dul(i) - du0(i)) /= 0.0) &
1024  fa_avg = uhtot_l(i) / (dul(i) - du0(i))
1025  if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))
1026  elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif
1027 
1028  bt_cont%FA_u_W0(i,j) = us%m_to_L*fa_0 ; bt_cont%FA_u_WW(i,j) = us%m_to_L*famt_l(i)
1029  if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%uBT_WW(i,j) = 0.0 ; else
1030  bt_cont%uBT_WW(i,j) = us%m_s_to_L_T*(1.5 * (dul(i) - du0(i))) * &
1031  ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
1032  endif
1033 
1034  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1035  if ((dur(i) - du0(i)) /= 0.0) &
1036  fa_avg = uhtot_r(i) / (dur(i) - du0(i))
1037  if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))
1038  elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif
1039 
1040  bt_cont%FA_u_E0(i,j) = us%m_to_L*fa_0 ; bt_cont%FA_u_EE(i,j) = us%m_to_L*famt_r(i)
1041  if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%uBT_EE(i,j) = 0.0 ; else
1042  bt_cont%uBT_EE(i,j) = us%m_s_to_L_T*(1.5 * (dur(i) - du0(i))) * &
1043  ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
1044  endif
1045  else
1046  bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
1047  bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
1048  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
1049  endif ; enddo
1050 
1051 end subroutine set_zonal_bt_cont
1052 
1053 !> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.
1054 subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, &
1055  visc_rem_v, v_cor, vhbt_aux, v_cor_aux, BT_cont)
1056  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1057  type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1058  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [m s-1].
1059  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to
1060  !! calculate fluxes [H ~> m or kg m-2].
1061  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh !< Volume flux through meridional
1062  !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
1063  real, intent(in) :: dt !< Time increment [s].
1064  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1065  type(continuity_ppm_cs), pointer :: CS !< This module's control structure.
1066  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1067  type(ocean_obc_type), optional, pointer :: OBC !< Open boundary condition type
1068  !! specifies whether, where, and what open boundary conditions are used.
1069  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1070  optional, intent(in) :: visc_rem_v !< Both the fraction of the momentum
1071  !! originally in a layer that remains after a time-step of viscosity,
1072  !! and the fraction of a time-step's worth of a barotropic acceleration
1073  !! that a layer experiences after viscosity is applied. Nondimensional between
1074  !! 0 (at the bottom) and 1 (far above the bottom).
1075  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt !< The summed volume flux through
1076  !< meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
1077  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt_aux !< A second set of summed volume fluxes
1078  !! through meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
1079  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1080  optional, intent(out) :: v_cor
1081  !< The meridional velocitiess (v with a barotropic correction)
1082  !! that give vhbt as the depth-integrated transport [m s-1].
1083  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1084  optional, intent(out) :: v_cor_aux
1085  !< The meridional velocities (v with a barotropic correction)
1086  !! that give vhbt_aux as the depth-integrated transports [m s-1].
1087  type(bt_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
1088  !! the effective open face areas as a function of barotropic flow.
1089  ! Local variables
1090  real, dimension(SZI_(G),SZK_(G)) :: &
1091  dvhdv ! Partial derivative of vh with v [H m ~> m2 or kg m-1].
1092  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1093  h_L, h_R ! Left and right face thicknesses [H ~> m or kg m-2].
1094  real, dimension(SZI_(G)) :: &
1095  dv, & ! Corrective barotropic change in the velocity [m s-1].
1096  dv_min_CFL, & ! Min/max limits on dv correction
1097  dv_max_CFL, & ! to avoid CFL violations
1098  dvhdv_tot_0, & ! Summed partial derivative of vh with v [H m ~> m2 or kg m-1].
1099  vh_tot_0, & ! Summed transport with no barotropic correction [H m2 s-1 ~> m3 s-1 or kg s-1].
1100  visc_rem_max ! The column maximum of visc_rem.
1101  logical, dimension(SZI_(G)) :: do_I
1102  real, dimension(SZI_(G)) :: FAvi ! A list of sums of meridional face areas [H m ~> m2 or kg m-1].
1103  real :: FA_v ! A sum of meridional face areas [H m ~> m2 or kg m-1].
1104  real, dimension(SZI_(G),SZK_(G)) :: &
1105  visc_rem ! A 2-D copy of visc_rem_v or an array of 1's.
1106  real :: I_vrm ! 1.0 / visc_rem_max, nondim.
1107  real :: CFL_dt ! The maximum CFL ratio of the adjusted velocities divided by
1108  ! the time step [s-1].
1109  real :: I_dt ! 1.0 / dt [s-1].
1110  real :: dv_lim ! The velocity change that give a relative CFL of 1 [m s-1].
1111  real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [m].
1112  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1113  logical :: do_aux, local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
1114  logical :: local_Flather_OBC, is_simple, local_open_BC
1115  type(obc_segment_type), pointer :: segment => null()
1116 
1117  do_aux = (present(vhbt_aux) .and. present(v_cor_aux))
1118  use_visc_rem = present(visc_rem_v)
1119  local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
1120  local_open_bc = .false.
1121  if (present(bt_cont)) set_bt_cont = (associated(bt_cont))
1122  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
1123  local_specified_bc = obc%specified_v_BCs_exist_globally
1124  local_flather_obc = obc%Flather_v_BCs_exist_globally
1125  local_open_bc = obc%open_v_BCs_exist_globally
1126  endif ; endif ; endif
1127  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1128 
1129  cfl_dt = cs%CFL_limit_adjust / dt
1130  i_dt = 1.0 / dt
1131  if (cs%aggress_adjust) cfl_dt = i_dt
1132 
1133  call cpu_clock_begin(id_clock_update)
1134 !$OMP parallel do default(none) shared(nz,ish,ieh,jsh,jeh,h_in,h_L,h_R,G,GV,LB,CS,visc_rem,OBC)
1135  do k=1,nz
1136  ! This sets h_L and h_R.
1137  if (cs%upwind_1st) then
1138  do j=jsh-1,jeh+1 ; do i=ish,ieh
1139  h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
1140  enddo ; enddo
1141  else
1142  call ppm_reconstruction_y(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
1143  2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
1144  endif
1145  do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo
1146  enddo
1147  call cpu_clock_end(id_clock_update)
1148 
1149  call cpu_clock_begin(id_clock_correct)
1150 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,v,h_in,h_L,h_R,vh,use_visc_rem, &
1151 !$OMP visc_rem_v,dt,G,GV,CS,local_specified_BC,OBC,vhbt,do_aux, &
1152 !$OMP set_BT_cont,CFL_dt,I_dt,v_cor,vhbt_aux, &
1153 !$OMP v_cor_aux,BT_cont, local_Flather_OBC ) &
1154 !$OMP private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0, &
1155 !$OMP dvhdv_tot_0,visc_rem_max,I_vrm,dv_lim,dy_N, &
1156 !$OMP is_simple,FAvi,dy_S,any_simple_OBC ) &
1157 !$OMP firstprivate(visc_rem)
1158  do j=jsh-1,jeh
1159  do i=ish,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ; enddo
1160  ! This sets vh and dvhdv.
1161  do k=1,nz
1162  if (use_visc_rem) then ; do i=ish,ieh
1163  visc_rem(i,k) = visc_rem_v(i,j,k)
1164  visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
1165  enddo ; endif
1166  call merid_flux_layer(v(:,j,k), h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
1167  vh(:,j,k), dvhdv(:,k), visc_rem(:,k), &
1168  dt, g, j, ish, ieh, do_i, cs%vol_CFL, obc)
1169  if (local_specified_bc) then
1170  do i=ish,ieh
1171  if (obc%segment(obc%segnum_v(i,j))%specified) &
1172  vh(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k)
1173  enddo
1174  endif
1175  enddo ! k-loop
1176  if ((.not.use_visc_rem) .or. (.not.cs%use_visc_rem_max)) then ; do i=ish,ieh
1177  visc_rem_max(i) = 1.0
1178  enddo ; endif
1179 
1180  if (present(vhbt) .or. do_aux .or. set_bt_cont) then
1181  ! Set limits on dv that will keep the CFL number between -1 and 1.
1182  ! This should be adequate to keep the root bracketed in all cases.
1183  do i=ish,ieh
1184  i_vrm = 0.0
1185  if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
1186  if (cs%vol_CFL) then
1187  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1188  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1189  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1190  dv_max_cfl(i) = 2.0 * (cfl_dt * dy_s) * i_vrm
1191  dv_min_cfl(i) = -2.0 * (cfl_dt * dy_n) * i_vrm
1192  vh_tot_0(i) = 0.0 ; dvhdv_tot_0(i) = 0.0
1193  enddo
1194  do k=1,nz ; do i=ish,ieh
1195  dvhdv_tot_0(i) = dvhdv_tot_0(i) + dvhdv(i,k)
1196  vh_tot_0(i) = vh_tot_0(i) + vh(i,j,k)
1197  enddo ; enddo
1198 
1199  if (use_visc_rem) then
1200  if (cs%aggress_adjust) then
1201  do k=1,nz ; do i=ish,ieh
1202  if (cs%vol_CFL) then
1203  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1204  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1205  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1206  dv_lim = 0.499*((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k)))
1207  if (dv_max_cfl(i) * visc_rem(i,k) > dv_lim) &
1208  dv_max_cfl(i) = dv_lim / visc_rem(i,k)
1209 
1210  dv_lim = 0.499*((-dy_n*cfl_dt - v(i,j,k)) + max(0.0,v(i,j+1,k)))
1211  if (dv_min_cfl(i) * visc_rem(i,k) < dv_lim) &
1212  dv_min_cfl(i) = dv_lim / visc_rem(i,k)
1213  enddo ; enddo
1214  else
1215  do k=1,nz ; do i=ish,ieh
1216  if (cs%vol_CFL) then
1217  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1218  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1219  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1220  if (dv_max_cfl(i) * visc_rem(i,k) > dy_s*cfl_dt - v(i,j,k)) &
1221  dv_max_cfl(i) = (dy_s*cfl_dt - v(i,j,k)) / visc_rem(i,k)
1222  if (dv_min_cfl(i) * visc_rem(i,k) < -dy_n*cfl_dt - v(i,j,k)) &
1223  dv_min_cfl(i) = -(dy_n*cfl_dt + v(i,j,k)) / visc_rem(i,k)
1224  enddo ; enddo
1225  endif
1226  else
1227  if (cs%aggress_adjust) then
1228  do k=1,nz ; do i=ish,ieh
1229  if (cs%vol_CFL) then
1230  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1231  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1232  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1233  dv_max_cfl(i) = min(dv_max_cfl(i), 0.499 * &
1234  ((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k))) )
1235  dv_min_cfl(i) = max(dv_min_cfl(i), 0.499 * &
1236  ((-dy_n*i_dt - v(i,j,k)) + max(0.0,v(i,j+1,k))) )
1237  enddo ; enddo
1238  else
1239  do k=1,nz ; do i=ish,ieh
1240  if (cs%vol_CFL) then
1241  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1242  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1243  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1244  dv_max_cfl(i) = min(dv_max_cfl(i), dy_s*cfl_dt - v(i,j,k))
1245  dv_min_cfl(i) = max(dv_min_cfl(i), -(dy_n*cfl_dt + v(i,j,k)))
1246  enddo ; enddo
1247  endif
1248  endif
1249  do i=ish,ieh
1250  dv_max_cfl(i) = max(dv_max_cfl(i),0.0)
1251  dv_min_cfl(i) = min(dv_min_cfl(i),0.0)
1252  enddo
1253 
1254  ! Up to this point, everything is shared between vhbt and vhbt_aux.
1255 
1256  any_simple_obc = .false.
1257  if (present(vhbt) .or. do_aux .or. set_bt_cont) then
1258  if (local_specified_bc .or. local_flather_obc) then ; do i=ish,ieh
1259  ! Avoid reconciling barotropic/baroclinic transports if transport is specified
1260  is_simple = obc%segment(obc%segnum_v(i,j))%specified
1261  do_i(i) = .not.(obc%segnum_v(i,j) /= obc_none .and. is_simple)
1262  any_simple_obc = any_simple_obc .or. is_simple
1263  enddo ; else ; do i=ish,ieh
1264  do_i(i) = .true.
1265  enddo ; endif
1266  endif
1267 
1268  if (present(vhbt)) then
1269  call meridional_flux_adjust(v, h_in, h_l, h_r, vhbt(:,j), vh_tot_0, &
1270  dvhdv_tot_0, dv, dv_max_cfl, dv_min_cfl, dt, g, &
1271  cs, visc_rem, j, ish, ieh, do_i, .true., vh, obc=obc)
1272 
1273  if (present(v_cor)) then ; do k=1,nz
1274  do i=ish,ieh ; v_cor(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
1275  if (local_specified_bc) then ; do i=ish,ieh
1276  if (obc%segment(obc%segnum_v(i,j))%specified) &
1277  v_cor(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1278  enddo ; endif
1279  enddo ; endif ! v-corrected
1280  endif
1281 
1282  if (do_aux) then
1283  call meridional_flux_adjust(v, h_in, h_l, h_r, vhbt_aux(:,j), vh_tot_0, &
1284  dvhdv_tot_0, dv, dv_max_cfl, dv_min_cfl, dt, g, &
1285  cs, visc_rem, j, ish, ieh, do_i, .false., obc=obc)
1286 
1287  do k=1,nz
1288  do i=ish,ieh ; v_cor_aux(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
1289  if (local_specified_bc) then ; do i=ish,ieh
1290  if (obc%segment(obc%segnum_v(i,j))%specified) &
1291  v_cor_aux(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1292  enddo ; endif
1293  enddo
1294  endif ! do_aux
1295 
1296  if (set_bt_cont) then
1297  call set_merid_bt_cont(v, h_in, h_l, h_r, bt_cont, vh_tot_0, dvhdv_tot_0,&
1298  dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1299  visc_rem_max, j, ish, ieh, do_i)
1300  if (any_simple_obc) then
1301  do i=ish,ieh
1302  do_i(i) = (obc%segment(obc%segnum_v(i,j))%specified)
1303  if (do_i(i)) favi(i) = gv%H_subroundoff*g%dx_Cv(i,j)
1304  enddo
1305  do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then
1306  if ((abs(obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
1307  (obc%segment(obc%segnum_v(i,j))%specified)) &
1308  favi(i) = favi(i) + &
1309  obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k) / &
1310  obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1311  endif ; enddo ; enddo
1312  do i=ish,ieh ; if (do_i(i)) then
1313  bt_cont%FA_v_S0(i,j) = us%m_to_L*favi(i) ; bt_cont%FA_v_N0(i,j) = us%m_to_L*favi(i)
1314  bt_cont%FA_v_SS(i,j) = us%m_to_L*favi(i) ; bt_cont%FA_v_NN(i,j) = us%m_to_L*favi(i)
1315  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1316  endif ; enddo
1317  endif
1318  endif ! set_BT_cont
1319 
1320  endif ! present(vhbt) or do_aux or set_BT_cont
1321  enddo ! j-loop
1322 
1323  if (local_open_bc .and. set_bt_cont) then
1324  do n = 1, obc%number_of_segments
1325  if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
1326  j = obc%segment(n)%HI%JsdB
1327  if (obc%segment(n)%direction == obc_direction_n) then
1328  do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
1329  fa_v = 0.0
1330  do k=1,nz ; fa_v = fa_v + h_in(i,j,k)*us%m_to_L*g%dx_Cv(i,j) ; enddo
1331  bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
1332  bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
1333  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1334  enddo
1335  else
1336  do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
1337  fa_v = 0.0
1338  do k=1,nz ; fa_v = fa_v + h_in(i,j+1,k)*us%m_to_L*g%dx_Cv(i,j) ; enddo
1339  bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
1340  bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
1341  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1342  enddo
1343  endif
1344  endif
1345  enddo
1346  endif
1347  call cpu_clock_end(id_clock_correct)
1348 
1349  if (set_bt_cont) then ; if (allocated(bt_cont%h_v)) then
1350  if (present(v_cor)) then
1351  call merid_face_thickness(v_cor, h_in, h_l, h_r, bt_cont%h_v, dt, g, lb, &
1352  cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
1353  else
1354  call merid_face_thickness(v, h_in, h_l, h_r, bt_cont%h_v, dt, g, lb, &
1355  cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
1356  endif
1357  endif ; endif
1358 
1359 end subroutine meridional_mass_flux
1360 
1361 !> Evaluates the meridional mass or volume fluxes in a layer.
1362 subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, J, &
1363  ish, ieh, do_I, vol_CFL, OBC)
1364  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1365  real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity [m s-1].
1366  real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the
1367  !! momentum originally in a layer that remains after a time-step
1368  !! of viscosity, and the fraction of a time-step's worth of a barotropic
1369  !! acceleration that a layer experiences after viscosity is applied.
1370  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1371  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,
1372  !! [H ~> m or kg m-2].
1373  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_L !< Left thickness in the reconstruction
1374  !! [H ~> m or kg m-2].
1375  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_R !< Right thickness in the reconstruction
1376  !! [H ~> m or kg m-2].
1377  real, dimension(SZI_(G)), intent(inout) :: vh !< Meridional mass or volume transport
1378  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
1379  real, dimension(SZI_(G)), intent(inout) :: dvhdv !< Partial derivative of vh with v
1380  !! [H m ~> m2 or kg m-1].
1381  real, intent(in) :: dt !< Time increment [s].
1382  integer, intent(in) :: j !< Spatial index.
1383  integer, intent(in) :: ish !< Start of index range.
1384  integer, intent(in) :: ieh !< End of index range.
1385  logical, dimension(SZI_(G)), intent(in) :: do_I !< Which i values to work on.
1386  logical, intent(in) :: vol_CFL !< If true, rescale the
1387  !! ratio of face areas to the cell areas when estimating the CFL number.
1388  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
1389  ! Local variables
1390  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1391  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1392  ! with the same units as h, i.e. [H ~> m or kg m-2].
1393  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1394  integer :: i
1395  logical :: local_open_BC
1396 
1397  local_open_bc = .false.
1398  if (present(obc)) then ; if (associated(obc)) then
1399  local_open_bc = obc%open_u_BCs_exist_globally
1400  endif ; endif
1401 
1402  do i=ish,ieh ; if (do_i(i)) then
1403  if (v(i) > 0.0) then
1404  if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1405  else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
1406  curv_3 = h_l(i,j) + h_r(i,j) - 2.0*h(i,j)
1407  vh(i) = g%dx_Cv(i,j) * v(i) * ( h_r(i,j) + cfl * &
1408  (0.5*(h_l(i,j) - h_r(i,j)) + curv_3*(cfl - 1.5)) )
1409  h_marg = h_r(i,j) + cfl * ((h_l(i,j) - h_r(i,j)) + &
1410  3.0*curv_3*(cfl - 1.0))
1411  elseif (v(i) < 0.0) then
1412  if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1413  else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
1414  curv_3 = h_l(i,j+1) + h_r(i,j+1) - 2.0*h(i,j+1)
1415  vh(i) = g%dx_Cv(i,j) * v(i) * ( h_l(i,j+1) + cfl * &
1416  (0.5*(h_r(i,j+1)-h_l(i,j+1)) + curv_3*(cfl - 1.5)) )
1417  h_marg = h_l(i,j+1) + cfl * ((h_r(i,j+1)-h_l(i,j+1)) + &
1418  3.0*curv_3*(cfl - 1.0))
1419  else
1420  vh(i) = 0.0
1421  h_marg = 0.5 * (h_l(i,j+1) + h_r(i,j))
1422  endif
1423  dvhdv(i) = g%dx_Cv(i,j) * h_marg * visc_rem(i)
1424  endif ; enddo
1425 
1426  if (local_open_bc) then
1427  do i=ish,ieh ; if (do_i(i)) then
1428  if (obc%segment(obc%segnum_v(i,j))%open) then
1429  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
1430  vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j)
1431  dvhdv(i) = g%dx_Cv(i,j) * h(i,j) * visc_rem(i)
1432  else
1433  vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j+1)
1434  dvhdv(i) = g%dx_Cv(i,j) * h(i,j+1) * visc_rem(i)
1435  endif
1436  endif
1437  endif ; enddo
1438  endif
1439 end subroutine merid_flux_layer
1440 
1441 !> Sets the effective interface thickness at each meridional velocity point.
1442 subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, LB, vol_CFL, &
1443  marginal, visc_rem_v, OBC)
1444  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1445  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [m s-1].
1446  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,
1447  !! [H ~> m or kg m-2].
1448  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the reconstruction,
1449  !! [H ~> m or kg m-2].
1450  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the reconstruction,
1451  !! [H ~> m or kg m-2].
1452  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h_v !< Thickness at meridional faces,
1453  !! [H ~> m or kg m-2].
1454  real, intent(in) :: dt !< Time increment [s].
1455  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1456  logical, intent(in) :: vol_CFL !< If true, rescale the ratio
1457  !! of face areas to the cell areas when estimating the CFL number.
1458  logical, intent(in) :: marginal !< If true, report the marginal
1459  !! face thicknesses; otherwise report transport-averaged thicknesses.
1460  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(in) :: visc_rem_v !< Both the fraction
1461  !! of the momentum originally in a layer that remains after a time-step of
1462  !! viscosity, and the fraction of a time-step's worth of a barotropic
1463  !! acceleration that a layer experiences after viscosity is applied.
1464  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1465  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
1466 
1467  ! Local variables
1468  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1469  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1470  ! with the same units as h_in.
1471  real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].
1472  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1473  logical :: local_open_BC
1474  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1475  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1476 
1477  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
1478  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
1479  if (v(i,j,k) > 0.0) then
1480  if (vol_cfl) then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1481  else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ; endif
1482  curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
1483  h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
1484  h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + &
1485  3.0*curv_3*(cfl - 1.0))
1486  elseif (v(i,j,k) < 0.0) then
1487  if (vol_cfl) then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1488  else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ; endif
1489  curv_3 = h_l(i,j+1,k) + h_r(i,j+1,k) - 2.0*h(i,j+1,k)
1490  h_avg = h_l(i,j+1,k) + cfl * (0.5*(h_r(i,j+1,k)-h_l(i,j+1,k)) + curv_3*(cfl - 1.5))
1491  h_marg = h_l(i,j+1,k) + cfl * ((h_r(i,j+1,k)-h_l(i,j+1,k)) + &
1492  3.0*curv_3*(cfl - 1.0))
1493  else
1494  h_avg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1495  ! The choice to use the arithmetic mean here is somewhat arbitrariy, but
1496  ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
1497  h_marg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1498  ! h_marg = (2.0 * h_L(i,j+1,k) * h_R(i,j,k)) / &
1499  ! (h_L(i,j+1,k) + h_R(i,j,k) + GV%H_subroundoff)
1500  endif
1501 
1502  if (marginal) then ; h_v(i,j,k) = h_marg
1503  else ; h_v(i,j,k) = h_avg ; endif
1504  enddo ; enddo ; enddo
1505 
1506  if (present(visc_rem_v)) then
1507  !$OMP parallel do default(shared)
1508  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
1509  h_v(i,j,k) = h_v(i,j,k) * visc_rem_v(i,j,k)
1510  enddo ; enddo ; enddo
1511  endif
1512 
1513  local_open_bc = .false.
1514  if (present(obc)) then ; if (associated(obc)) then
1515  local_open_bc = obc%open_u_BCs_exist_globally
1516  endif ; endif
1517  if (local_open_bc) then
1518  do n = 1, obc%number_of_segments
1519  if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
1520  j = obc%segment(n)%HI%JsdB
1521  if (obc%segment(n)%direction == obc_direction_n) then
1522  if (present(visc_rem_v)) then ; do k=1,nz
1523  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1524  h_v(i,j,k) = h(i,j,k) * visc_rem_v(i,j,k)
1525  enddo
1526  enddo ; else ; do k=1,nz
1527  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1528  h_v(i,j,k) = h(i,j,k)
1529  enddo
1530  enddo ; endif
1531  else
1532  if (present(visc_rem_v)) then ; do k=1,nz
1533  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1534  h_v(i,j,k) = h(i,j+1,k) * visc_rem_v(i,j,k)
1535  enddo
1536  enddo ; else ; do k=1,nz
1537  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1538  h_v(i,j,k) = h(i,j+1,k)
1539  enddo
1540  enddo ; endif
1541  endif
1542  endif
1543  enddo
1544  endif
1545 
1546 end subroutine merid_face_thickness
1547 
1548 !> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.
1549 subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0, &
1550  dv, dv_max_CFL, dv_min_CFL, dt, G, CS, visc_rem, &
1551  j, ish, ieh, do_I_in, full_precision, vh_3d, OBC)
1552  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1553  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1554  intent(in) :: v !< Meridional velocity [m s-1].
1555  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1556  intent(in) :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
1557  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),&
1558  intent(in) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
1559  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1560  intent(in) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
1561  real, dimension(SZI_(G),SZK_(G)), intent(in) :: visc_rem
1562  !< Both the fraction of the momentum originally
1563  !! in a layer that remains after a time-step of viscosity, and the
1564  !! fraction of a time-step's worth of a barotropic acceleration that
1565  !! a layer experiences after viscosity is applied. Non-dimensional
1566  !! between 0 (at the bottom) and 1 (far above the bottom).
1567  real, dimension(SZI_(G)), &
1568  optional, intent(in) :: vhbt !< The summed volume flux through meridional faces
1569  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
1570  real, dimension(SZI_(G)), intent(in) :: dv_max_CFL !< Maximum acceptable value of dv [m s-1].
1571  real, dimension(SZI_(G)), intent(in) :: dv_min_CFL !< Minimum acceptable value of dv [m s-1].
1572  real, dimension(SZI_(G)), intent(in) :: vh_tot_0 !< The summed transport with 0 adjustment
1573  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
1574  real, dimension(SZI_(G)), intent(in) :: dvhdv_tot_0 !< The partial derivative of dv_err with
1575  !! dv at 0 adjustment [H m ~> m2 or kg m-1].
1576  real, dimension(SZI_(G)), intent(out) :: dv !< The barotropic velocity adjustment [m s-1].
1577  real, intent(in) :: dt !< Time increment [s].
1578  type(continuity_ppm_cs), pointer :: CS !< This module's control structure.
1579  integer, intent(in) :: j !< Spatial index.
1580  integer, intent(in) :: ish !< Start of index range.
1581  integer, intent(in) :: ieh !< End of index range.
1582  logical, dimension(SZI_(G)), &
1583  intent(in) :: do_I_in !< A flag indicating which I values to work on.
1584  logical, optional, intent(in) :: full_precision !< A flag indicating how carefully to
1585  !! iterate. The default is .true. (more accurate).
1586  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1587  optional, intent(inout) :: vh_3d !< Volume flux through meridional
1588  !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
1589  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
1590  ! Local variables
1591  real, dimension(SZI_(G),SZK_(G)) :: &
1592  vh_aux, & ! An auxiliary meridional volume flux [H m s-1 ~> m2 s-1 or kg m-1 s-1].
1593  dvhdv ! Partial derivative of vh with v [H m ~> m2 or kg m-1].
1594  real, dimension(SZI_(G)) :: &
1595  vh_err, & ! Difference between vhbt and the summed vh [H m2 s-1 ~> m3 s-1 or kg s-1].
1596  vh_err_best, & ! The smallest value of vh_err found so far [H m2 s-1 ~> m3 s-1 or kg s-1].
1597  v_new, & ! The velocity with the correction added [m s-1].
1598  dvhdv_tot,&! Summed partial derivative of vh with u [H m ~> m2 or kg m-1].
1599  dv_min, & ! Min/max limits on dv correction based on CFL limits
1600  dv_max ! and previous iterations [m s-1].
1601  real :: dv_prev ! The previous value of dv [m s-1].
1602  real :: ddv ! The change in dv from the previous iteration [m s-1].
1603  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
1604  real :: tol_vel ! The tolerance for velocity in the current iteration [m s-1].
1605  integer :: i, k, nz, itt, max_itts = 20
1606  logical :: full_prec, domore, do_I(SZI_(G))
1607 
1608  nz = g%ke
1609  full_prec = .true. ; if (present(full_precision)) full_prec = full_precision
1610 
1611  vh_aux(:,:) = 0.0 ; dvhdv(:,:) = 0.0
1612 
1613  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
1614  vh_aux(i,k) = vh_3d(i,j,k)
1615  enddo ; enddo ; endif
1616 
1617  do i=ish,ieh
1618  dv(i) = 0.0 ; do_i(i) = do_i_in(i)
1619  dv_max(i) = dv_max_cfl(i) ; dv_min(i) = dv_min_cfl(i)
1620  vh_err(i) = vh_tot_0(i) - vhbt(i) ; dvhdv_tot(i) = dvhdv_tot_0(i)
1621  vh_err_best(i) = abs(vh_err(i))
1622  enddo
1623 
1624  do itt=1,max_itts
1625  if (full_prec) then
1626  select case (itt)
1627  case (:1) ; tol_eta = 1e-6 * cs%tol_eta
1628  case (2) ; tol_eta = 1e-4 * cs%tol_eta
1629  case (3) ; tol_eta = 1e-2 * cs%tol_eta
1630  case default ; tol_eta = cs%tol_eta
1631  end select
1632  else
1633  tol_eta = cs%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
1634  endif
1635  tol_vel = cs%tol_vel
1636 
1637  do i=ish,ieh
1638  if (vh_err(i) > 0.0) then ; dv_max(i) = dv(i)
1639  elseif (vh_err(i) < 0.0) then ; dv_min(i) = dv(i)
1640  else ; do_i(i) = .false. ; endif
1641  enddo
1642  domore = .false.
1643  do i=ish,ieh ; if (do_i(i)) then
1644  if ((dt*min(g%IareaT(i,j),g%IareaT(i,j+1))*abs(vh_err(i)) > tol_eta) .or.&
1645  (cs%better_iter .and. ((abs(vh_err(i)) > tol_vel * dvhdv_tot(i)) .or.&
1646  (abs(vh_err(i)) > vh_err_best(i))) )) then
1647  ! Use Newton's method, provided it stays bounded. Otherwise bisect
1648  ! the value with the appropriate bound.
1649  if (full_prec) then
1650  ddv = -vh_err(i) / dvhdv_tot(i)
1651  dv_prev = dv(i)
1652  dv(i) = dv(i) + ddv
1653  if (abs(ddv) < 1.0e-15*abs(dv(i))) then
1654  do_i(i) = .false. ! ddv is small enough to quit.
1655  elseif (ddv > 0.0) then
1656  if (dv(i) >= dv_max(i)) then
1657  dv(i) = 0.5*(dv_prev + dv_max(i))
1658  if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_i(i) = .false.
1659  endif
1660  else ! dvv(i) < 0.0
1661  if (dv(i) <= dv_min(i)) then
1662  dv(i) = 0.5*(dv_prev + dv_min(i))
1663  if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_i(i) = .false.
1664  endif
1665  endif
1666  else
1667  ! Use Newton's method, provided it stays bounded, just like above.
1668  dv(i) = dv(i) - vh_err(i) / dvhdv_tot(i)
1669  if ((dv(i) >= dv_max(i)) .or. (dv(i) <= dv_min(i))) &
1670  dv(i) = 0.5*(dv_max(i) + dv_min(i))
1671  endif
1672  if (do_i(i)) domore = .true.
1673  else
1674  do_i(i) = .false.
1675  endif
1676  endif ; enddo
1677  if (.not.domore) exit
1678 
1679  if ((itt < max_itts) .or. present(vh_3d)) then ; do k=1,nz
1680  do i=ish,ieh ; v_new(i) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
1681  call merid_flux_layer(v_new, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
1682  vh_aux(:,k), dvhdv(:,k), visc_rem(:,k), &
1683  dt, g, j, ish, ieh, do_i, cs%vol_CFL, obc)
1684  enddo ; endif
1685 
1686  if (itt < max_itts) then
1687  do i=ish,ieh
1688  vh_err(i) = -vhbt(i) ; dvhdv_tot(i) = 0.0
1689  enddo
1690  do k=1,nz ; do i=ish,ieh
1691  vh_err(i) = vh_err(i) + vh_aux(i,k)
1692  dvhdv_tot(i) = dvhdv_tot(i) + dvhdv(i,k)
1693  enddo ; enddo
1694  do i=ish,ieh
1695  vh_err_best(i) = min(vh_err_best(i), abs(vh_err(i)))
1696  enddo
1697  endif
1698  enddo ! itt-loop
1699  ! If there are any faces which have not converged to within the tolerance,
1700  ! so-be-it, or else use a final upwind correction?
1701  ! This never seems to happen with 20 iterations as max_itt.
1702 
1703  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
1704  vh_3d(i,j,k) = vh_aux(i,k)
1705  enddo ; enddo ; endif
1706 
1707 end subroutine meridional_flux_adjust
1708 
1709 !> Sets of a structure that describes the meridional barotropic volume or mass fluxes as a
1710 !! function of barotropic flow to agree closely with the sum of the layer's transports.
1711 subroutine set_merid_bt_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, &
1712  dv_max_CFL, dv_min_CFL, dt, G, US, CS, visc_rem, &
1713  visc_rem_max, j, ish, ieh, do_I)
1714  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1715  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [m s-1].
1716  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to calculate fluxes,
1717  !! [H ~> m or kg m-2].
1718  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the reconstruction,
1719  !! [H ~> m or kg m-2].
1720  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the reconstruction,
1721  !! [H ~> m or kg m-2].
1722  type(bt_cont_type), intent(inout) :: BT_cont !< A structure with elements
1723  !! that describe the effective open face areas as a function of barotropic flow.
1724  real, dimension(SZI_(G)), intent(in) :: vh_tot_0 !< The summed transport
1725  !! with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
1726  real, dimension(SZI_(G)), intent(in) :: dvhdv_tot_0 !< The partial derivative
1727  !! of du_err with dv at 0 adjustment [H m ~> m2 or kg m-1].
1728  real, dimension(SZI_(G)), intent(in) :: dv_max_CFL !< Maximum acceptable value of dv [m s-1].
1729  real, dimension(SZI_(G)), intent(in) :: dv_min_CFL !< Minimum acceptable value of dv [m s-1].
1730  real, intent(in) :: dt !< Time increment [s].
1731  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1732  type(continuity_ppm_cs), pointer :: CS !< This module's control structure.
1733  real, dimension(SZI_(G),SZK_(G)), intent(in) :: visc_rem !< Both the fraction of the
1734  !! momentum originally in a layer that remains after a time-step
1735  !! of viscosity, and the fraction of a time-step's worth of a barotropic
1736  !! acceleration that a layer experiences after viscosity is applied.
1737  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1738  real, dimension(SZI_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem.
1739  integer, intent(in) :: j !< Spatial index.
1740  integer, intent(in) :: ish !< Start of index range.
1741  integer, intent(in) :: ieh !< End of index range.
1742  logical, dimension(SZI_(G)), intent(in) :: do_I !< A logical flag indicating
1743  !! which I values to work on.
1744  ! Local variables
1745  real, dimension(SZI_(G)) :: &
1746  dv0, & ! The barotropic velocity increment that gives 0 transport [m s-1].
1747  dvL, dvR, & ! The barotropic velocity increments that give the southerly
1748  ! (dvL) and northerly (dvR) test velocities.
1749  zeros, & ! An array of full of 0's.
1750  dv_cfl, & ! The velocity increment that corresponds to CFL_min [m s-1].
1751  v_l, v_r, & ! The southerly (v_L), northerly (v_R), and zero-barotropic
1752  v_0, & ! transport (v_0) layer test velocities [m s-1].
1753  fa_marg_l, & ! The effective layer marginal face areas with the southerly
1754  fa_marg_r, & ! (_L), northerly (_R), and zero-barotropic (_0) test
1755  fa_marg_0, & ! velocities [H m ~> m2 or kg m-1].
1756  vh_l, vh_r, & ! The layer transports with the southerly (_L), northerly (_R)
1757  vh_0, & ! and zero-barotropic (_0) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
1758  famt_l, famt_r, & ! The summed effective marginal face areas for the 3
1759  famt_0, & ! test velocities [H m ~> m2 or kg m-1].
1760  vhtot_l, & ! The summed transport with the southerly (vhtot_L) and
1761  vhtot_r ! and northerly (vhtot_R) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
1762  real :: FA_0 ! The effective face area with 0 barotropic transport [H m ~> m2 or kg m-1].
1763  real :: FA_avg ! The average effective face area [H m ~> m2 or kg m-1], nominally given by
1764  ! the realized transport divided by the barotropic velocity.
1765  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim] This
1766  ! limiting is necessary to keep the inverse of visc_rem
1767  ! from leading to large CFL numbers.
1768  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
1769  ! in finding the barotropic velocity that changes the
1770  ! flow direction. This is necessary to keep the inverse
1771  ! of visc_rem from leading to large CFL numbers.
1772  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
1773  ! flow is truly upwind [nondim]
1774  real :: Idt ! The inverse of the time step [s-1].
1775  logical :: domore
1776  integer :: i, k, nz
1777 
1778  nz = g%ke ; idt = 1.0/dt
1779  min_visc_rem = 0.1 ; cfl_min = 1e-6
1780 
1781  ! Diagnose the zero-transport correction, dv0.
1782  do i=ish,ieh ; zeros(i) = 0.0 ; enddo
1783  call meridional_flux_adjust(v, h_in, h_l, h_r, zeros, vh_tot_0, &
1784  dvhdv_tot_0, dv0, dv_max_cfl, dv_min_cfl, dt, g, &
1785  cs, visc_rem, j, ish, ieh, do_i, .true.)
1786 
1787  ! Determine the southerly- and northerly- fluxes. Choose a sufficiently
1788  ! negative velocity correction for the northerly-flux, and a sufficiently
1789  ! positive correction for the southerly-flux.
1790  domore = .false.
1791  do i=ish,ieh ; if (do_i(i)) then
1792  domore = .true.
1793  dv_cfl(i) = (cfl_min * idt) * g%dyCv(i,j)
1794  dvr(i) = min(0.0,dv0(i) - dv_cfl(i))
1795  dvl(i) = max(0.0,dv0(i) + dv_cfl(i))
1796  famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
1797  vhtot_l(i) = 0.0 ; vhtot_r(i) = 0.0
1798  endif ; enddo
1799 
1800  if (.not.domore) then
1801  do k=1,nz ; do i=ish,ieh
1802  bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
1803  bt_cont%vBT_SS(i,j) = 0.0
1804  bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
1805  bt_cont%vBT_NN(i,j) = 0.0
1806  enddo ; enddo
1807  return
1808  endif
1809 
1810  do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then
1811  visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
1812  if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
1813  if (v(i,j,k) + dvr(i)*visc_rem_lim > -dv_cfl(i)*visc_rem(i,k)) &
1814  dvr(i) = -(v(i,j,k) + dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
1815  if (v(i,j,k) + dvl(i)*visc_rem_lim < dv_cfl(i)*visc_rem(i,k)) &
1816  dvl(i) = -(v(i,j,k) - dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
1817  endif
1818  endif ; enddo ; enddo
1819  do k=1,nz
1820  do i=ish,ieh ; if (do_i(i)) then
1821  v_l(i) = v(i,j,k) + dvl(i) * visc_rem(i,k)
1822  v_r(i) = v(i,j,k) + dvr(i) * visc_rem(i,k)
1823  v_0(i) = v(i,j,k) + dv0(i) * visc_rem(i,k)
1824  endif ; enddo
1825  call merid_flux_layer(v_0, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_0, &
1826  fa_marg_0, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
1827  cs%vol_CFL)
1828  call merid_flux_layer(v_l, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_l, &
1829  fa_marg_l, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
1830  cs%vol_CFL)
1831  call merid_flux_layer(v_r, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_r, &
1832  fa_marg_r, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
1833  cs%vol_CFL)
1834  do i=ish,ieh ; if (do_i(i)) then
1835  famt_0(i) = famt_0(i) + fa_marg_0(i)
1836  famt_l(i) = famt_l(i) + fa_marg_l(i)
1837  famt_r(i) = famt_r(i) + fa_marg_r(i)
1838  vhtot_l(i) = vhtot_l(i) + vh_l(i)
1839  vhtot_r(i) = vhtot_r(i) + vh_r(i)
1840  endif ; enddo
1841  enddo
1842  do i=ish,ieh ; if (do_i(i)) then
1843  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1844  if ((dvl(i) - dv0(i)) /= 0.0) &
1845  fa_avg = vhtot_l(i) / (dvl(i) - dv0(i))
1846  if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))
1847  elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif
1848  bt_cont%FA_v_S0(i,j) = us%m_to_L*fa_0 ; bt_cont%FA_v_SS(i,j) = us%m_to_L*famt_l(i)
1849  if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%vBT_SS(i,j) = 0.0 ; else
1850  bt_cont%vBT_SS(i,j) = us%m_s_to_L_T*(1.5 * (dvl(i) - dv0(i))) * &
1851  ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
1852  endif
1853 
1854  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1855  if ((dvr(i) - dv0(i)) /= 0.0) &
1856  fa_avg = vhtot_r(i) / (dvr(i) - dv0(i))
1857  if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))
1858  elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif
1859  bt_cont%FA_v_N0(i,j) = us%m_to_L*fa_0 ; bt_cont%FA_v_NN(i,j) = us%m_to_L*famt_r(i)
1860  if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%vBT_NN(i,j) = 0.0 ; else
1861  bt_cont%vBT_NN(i,j) = us%m_s_to_L_T*(1.5 * (dvr(i) - dv0(i))) * &
1862  ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
1863  endif
1864  else
1865  bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
1866  bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
1867  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1868  endif ; enddo
1869 
1870 end subroutine set_merid_bt_cont
1871 
1872 !> Calculates left/right edge values for PPM reconstruction.
1873 subroutine ppm_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
1874  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1875  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
1876  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_L !< Left thickness in the reconstruction,
1877  !! [H ~> m or kg m-2].
1878  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_R !< Right thickness in the reconstruction,
1879  !! [H ~> m or kg m-2].
1880  type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
1881  real, intent(in) :: h_min !< The minimum thickness
1882  !! that can be obtained by a concave parabolic fit.
1883  logical, optional, intent(in) :: monotonic !< If true, use the
1884  !! Colella & Woodward monotonic limiter.
1885  !! Otherwise use a simple positive-definite limiter.
1886  logical, optional, intent(in) :: simple_2nd !< If true, use the
1887  !! arithmetic mean thicknesses as the default edge values
1888  !! for a simple 2nd order scheme.
1889  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
1890 
1891  ! Local variables with useful mnemonic names.
1892  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1893  real, parameter :: oneSixth = 1./6.
1894  real :: h_ip1, h_im1
1895  real :: dMx, dMn
1896  logical :: use_CW84, use_2nd
1897  character(len=256) :: mesg
1898  integer :: i, j, isl, iel, jsl, jel, n, stencil
1899  logical :: local_open_BC
1900  type(obc_segment_type), pointer :: segment => null()
1901 
1902  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
1903  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1904 
1905  local_open_bc = .false.
1906  if (present(obc)) then ; if (associated(obc)) then
1907  local_open_bc = obc%open_u_BCs_exist_globally
1908  endif ; endif
1909 
1910  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1911 
1912  ! This is the stencil of the reconstruction, not the scheme overall.
1913  stencil = 2 ; if (use_2nd) stencil = 1
1914 
1915  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
1916  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1917  & "x-halo that needs to be increased by ",i2,".")') &
1918  stencil + max(g%isd-isl,iel-g%ied)
1919  call mom_error(fatal,mesg)
1920  endif
1921  if ((jsl < g%jsd) .or. (jel > g%jed)) then
1922  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1923  & "y-halo that needs to be increased by ",i2,".")') &
1924  max(g%jsd-jsl,jel-g%jed)
1925  call mom_error(fatal,mesg)
1926  endif
1927 
1928  if (use_2nd) then
1929  do j=jsl,jel ; do i=isl,iel
1930  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1931  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1932  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1933  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1934  enddo ; enddo
1935  else
1936  do j=jsl,jel ; do i=isl-1,iel+1
1937  if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
1938  slp(i,j) = 0.0
1939  else
1940  ! This uses a simple 2nd order slope.
1941  slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1942  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1943  dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1944  dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1945  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1946  ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
1947  endif
1948  enddo ; enddo
1949 
1950  if (local_open_bc) then
1951  do n=1, obc%number_of_segments
1952  segment => obc%segment(n)
1953  if (.not. segment%on_pe) cycle
1954  if (segment%direction == obc_direction_e .or. &
1955  segment%direction == obc_direction_w) then
1956  i=segment%HI%IsdB
1957  do j=segment%HI%jsd,segment%HI%jed
1958  slp(i+1,j) = 0.0
1959  slp(i,j) = 0.0
1960  enddo
1961  endif
1962  enddo
1963  endif
1964 
1965  do j=jsl,jel ; do i=isl,iel
1966  ! Neighboring values should take into account any boundaries. The 3
1967  ! following sets of expressions are equivalent.
1968  ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
1969  ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
1970  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1971  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1972  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
1973  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1974  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1975  enddo ; enddo
1976  endif
1977 
1978  if (local_open_bc) then
1979  do n=1, obc%number_of_segments
1980  segment => obc%segment(n)
1981  if (.not. segment%on_pe) cycle
1982  if (segment%direction == obc_direction_e) then
1983  i=segment%HI%IsdB
1984  do j=segment%HI%jsd,segment%HI%jed
1985  h_l(i+1,j) = h_in(i,j)
1986  h_r(i+1,j) = h_in(i,j)
1987  h_l(i,j) = h_in(i,j)
1988  h_r(i,j) = h_in(i,j)
1989  enddo
1990  elseif (segment%direction == obc_direction_w) then
1991  i=segment%HI%IsdB
1992  do j=segment%HI%jsd,segment%HI%jed
1993  h_l(i,j) = h_in(i+1,j)
1994  h_r(i,j) = h_in(i+1,j)
1995  h_l(i+1,j) = h_in(i+1,j)
1996  h_r(i+1,j) = h_in(i+1,j)
1997  enddo
1998  endif
1999  enddo
2000  endif
2001 
2002  if (use_cw84) then
2003  call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
2004  else
2005  call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
2006  endif
2007 
2008  return
2009 end subroutine ppm_reconstruction_x
2010 
2011 !> Calculates left/right edge values for PPM reconstruction.
2012 subroutine ppm_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
2013  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2014  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2015  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_L !< Left thickness in the reconstruction,
2016  !! [H ~> m or kg m-2].
2017  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_R !< Right thickness in the reconstruction,
2018  !! [H ~> m or kg m-2].
2019  type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
2020  real, intent(in) :: h_min !< The minimum thickness
2021  !! that can be obtained by a concave parabolic fit.
2022  logical, optional, intent(in) :: monotonic !< If true, use the
2023  !! Colella & Woodward monotonic limiter.
2024  !! Otherwise use a simple positive-definite limiter.
2025  logical, optional, intent(in) :: simple_2nd !< If true, use the
2026  !! arithmetic mean thicknesses as the default edge values
2027  !! for a simple 2nd order scheme.
2028  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
2029 
2030  ! Local variables with useful mnemonic names.
2031  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
2032  real, parameter :: oneSixth = 1./6.
2033  real :: h_jp1, h_jm1
2034  real :: dMx, dMn
2035  logical :: use_CW84, use_2nd
2036  character(len=256) :: mesg
2037  integer :: i, j, isl, iel, jsl, jel, n, stencil
2038  logical :: local_open_BC
2039  type(obc_segment_type), pointer :: segment => null()
2040 
2041  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
2042  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
2043 
2044  local_open_bc = .false.
2045  if (present(obc)) then ; if (associated(obc)) then
2046  local_open_bc = obc%open_u_BCs_exist_globally
2047  endif ; endif
2048 
2049  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
2050 
2051  ! This is the stencil of the reconstruction, not the scheme overall.
2052  stencil = 2 ; if (use_2nd) stencil = 1
2053 
2054  if ((isl < g%isd) .or. (iel > g%ied)) then
2055  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
2056  & "x-halo that needs to be increased by ",i2,".")') &
2057  max(g%isd-isl,iel-g%ied)
2058  call mom_error(fatal,mesg)
2059  endif
2060  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
2061  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
2062  & "y-halo that needs to be increased by ",i2,".")') &
2063  stencil + max(g%jsd-jsl,jel-g%jed)
2064  call mom_error(fatal,mesg)
2065  endif
2066 
2067  if (use_2nd) then
2068  do j=jsl,jel ; do i=isl,iel
2069  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2070  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2071  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
2072  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
2073  enddo ; enddo
2074  else
2075  do j=jsl-1,jel+1 ; do i=isl,iel
2076  if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
2077  slp(i,j) = 0.0
2078  else
2079  ! This uses a simple 2nd order slope.
2080  slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
2081  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
2082  dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
2083  dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
2084  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
2085  ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
2086  endif
2087  enddo ; enddo
2088 
2089  if (local_open_bc) then
2090  do n=1, obc%number_of_segments
2091  segment => obc%segment(n)
2092  if (.not. segment%on_pe) cycle
2093  if (segment%direction == obc_direction_s .or. &
2094  segment%direction == obc_direction_n) then
2095  j=segment%HI%JsdB
2096  do i=segment%HI%isd,segment%HI%ied
2097  slp(i,j+1) = 0.0
2098  slp(i,j) = 0.0
2099  enddo
2100  endif
2101  enddo
2102  endif
2103 
2104  do j=jsl,jel ; do i=isl,iel
2105  ! Neighboring values should take into account any boundaries. The 3
2106  ! following sets of expressions are equivalent.
2107  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2108  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2109  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2110  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2111  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2112  enddo ; enddo
2113  endif
2114 
2115  if (local_open_bc) then
2116  do n=1, obc%number_of_segments
2117  segment => obc%segment(n)
2118  if (.not. segment%on_pe) cycle
2119  if (segment%direction == obc_direction_n) then
2120  j=segment%HI%JsdB
2121  do i=segment%HI%isd,segment%HI%ied
2122  h_l(i,j+1) = h_in(i,j)
2123  h_r(i,j+1) = h_in(i,j)
2124  h_l(i,j) = h_in(i,j)
2125  h_r(i,j) = h_in(i,j)
2126  enddo
2127  elseif (segment%direction == obc_direction_s) then
2128  j=segment%HI%JsdB
2129  do i=segment%HI%isd,segment%HI%ied
2130  h_l(i,j) = h_in(i,j+1)
2131  h_r(i,j) = h_in(i,j+1)
2132  h_l(i,j+1) = h_in(i,j+1)
2133  h_r(i,j+1) = h_in(i,j+1)
2134  enddo
2135  endif
2136  enddo
2137  endif
2138 
2139  if (use_cw84) then
2140  call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
2141  else
2142  call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
2143  endif
2144 
2145  return
2146 end subroutine ppm_reconstruction_y
2147 
2148 !> This subroutine limits the left/right edge values of the PPM reconstruction
2149 !! to give a reconstruction that is positive-definite. Here this is
2150 !! reinterpreted as giving a constant thickness if the mean thickness is less
2151 !! than h_min, with a minimum of h_min otherwise.
2152 subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
2153  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2154  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2155  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
2156  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
2157  real, intent(in) :: h_min !< The minimum thickness
2158  !! that can be obtained by a concave parabolic fit.
2159  integer, intent(in) :: iis !< Start of i index range.
2160  integer, intent(in) :: iie !< End of i index range.
2161  integer, intent(in) :: jis !< Start of j index range.
2162  integer, intent(in) :: jie !< End of j index range.
2163 
2164 ! Local variables
2165  real :: curv, dh, scale
2166  character(len=256) :: mesg
2167  integer :: i,j
2168 
2169  do j=jis,jie ; do i=iis,iie
2170  ! This limiter prevents undershooting minima within the domain with
2171  ! values less than h_min.
2172  curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2173  if (curv > 0.0) then ! Only minima are limited.
2174  dh = h_r(i,j) - h_l(i,j)
2175  if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
2176  if (h_in(i,j) <= h_min) then
2177  h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2178  elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
2179  ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
2180  ! be limited in this case. 0 < scale < 1.
2181  scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2182  h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2183  h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2184  endif
2185  endif
2186  endif
2187  enddo ; enddo
2188 
2189 end subroutine ppm_limit_pos
2190 
2191 !> This subroutine limits the left/right edge values of the PPM reconstruction
2192 !! according to the monotonic prescription of Colella and Woodward, 1984.
2193 subroutine ppm_limit_cw84(h_in, h_L, h_R, G, iis, iie, jis, jie)
2194  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2195  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2196  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction,
2197  !! [H ~> m or kg m-2].
2198  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction,
2199  !! [H ~> m or kg m-2].
2200  integer, intent(in) :: iis !< Start of i index range.
2201  integer, intent(in) :: iie !< End of i index range.
2202  integer, intent(in) :: jis !< Start of j index range.
2203  integer, intent(in) :: jie !< End of j index range.
2204 
2205  ! Local variables
2206  real :: h_i, RLdiff, RLdiff2, RLmean, FunFac
2207  character(len=256) :: mesg
2208  integer :: i,j
2209 
2210  do j=jis,jie ; do i=iis,iie
2211  ! This limiter monotonizes the parabola following
2212  ! Colella and Woodward, 1984, Eq. 1.10
2213  h_i = h_in(i,j)
2214  if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. ) then
2215  h_l(i,j) = h_i ; h_r(i,j) = h_i
2216  else
2217  rldiff = h_r(i,j) - h_l(i,j) ! Difference of edge values
2218  rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) ) ! Mean of edge values
2219  funfac = 6. * rldiff * ( h_i - rlmean ) ! Some funny factor
2220  rldiff2 = rldiff * rldiff ! Square of difference
2221  if ( funfac > rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)
2222  if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)
2223  endif
2224  enddo ; enddo
2225 
2226  return
2227 end subroutine ppm_limit_cw84
2228 
2229 !> Return the maximum ratio of a/b or maxrat.
2230 function ratio_max(a, b, maxrat) result(ratio)
2231  real, intent(in) :: a !< Numerator
2232  real, intent(in) :: b !< Denominator
2233  real, intent(in) :: maxrat !< Maximum value of ratio.
2234  real :: ratio !< Return value.
2235 
2236  if (abs(a) > abs(maxrat*b)) then
2237  ratio = maxrat
2238  else
2239  ratio = a / b
2240  endif
2241 end function ratio_max
2242 
2243 !> Initializes continuity_ppm_cs
2244 subroutine continuity_ppm_init(Time, G, GV, param_file, diag, CS)
2245  type(time_type), target, intent(in) :: time !< Time increment [s].
2246  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2247  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
2248  type(param_file_type), intent(in) :: param_file !< A structure indicating
2249  !! the open file to parse for model parameter values.
2250  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
2251  !! regulate diagnostic output.
2252  type(continuity_ppm_cs), pointer :: cs !< Module's control structure.
2253 !> This include declares and sets the variable "version".
2254 #include "version_variable.h"
2255  real :: tol_eta_m ! An unscaled version of tol_eta [m].
2256  character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name.
2257 
2258  if (associated(cs)) then
2259  call mom_error(warning, "continuity_PPM_init called with associated control structure.")
2260  return
2261  endif
2262  allocate(cs)
2263 
2264 ! Read all relevant parameters and write them to the model log.
2265  call log_version(param_file, mdl, version, "")
2266  call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", cs%monotonic, &
2267  "If true, CONTINUITY_PPM uses the Colella and Woodward "//&
2268  "monotonic limiter. The default (false) is to use a "//&
2269  "simple positive definite limiter.", default=.false.)
2270  call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &
2271  "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2272  "(arithmetic mean) interpolation of the edge values. "//&
2273  "This may give better PV conservation properties. While "//&
2274  "it formally reduces the accuracy of the continuity "//&
2275  "solver itself in the strongly advective limit, it does "//&
2276  "not reduce the overall order of accuracy of the dynamic "//&
2277  "core.", default=.false.)
2278  call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
2279  "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
2280  "continuity solver. This scheme is highly diffusive "//&
2281  "but may be useful for debugging or in single-column "//&
2282  "mode where its minimal stencil is useful.", default=.false.)
2283  call get_param(param_file, mdl, "ETA_TOLERANCE", cs%tol_eta, &
2284  "The tolerance for the differences between the "//&
2285  "barotropic and baroclinic estimates of the sea surface "//&
2286  "height due to the fluxes through each face. The total "//&
2287  "tolerance for SSH is 4 times this value. The default "//&
2288  "is 0.5*NK*ANGSTROM, and this should not be set less "//&
2289  "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=gv%m_to_H, &
2290  default=0.5*g%ke*gv%Angstrom_m, unscaled=tol_eta_m)
2291 
2292  call get_param(param_file, mdl, "ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
2293  "The tolerance for free-surface height discrepancies "//&
2294  "between the barotropic solution and the sum of the "//&
2295  "layer thicknesses when calculating the auxiliary "//&
2296  "corrected velocities. By default, this is the same as "//&
2297  "ETA_TOLERANCE, but can be made larger for efficiency.", &
2298  units="m", default=tol_eta_m, scale=gv%m_to_H)
2299  call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &
2300  "The tolerance for barotropic velocity discrepancies "//&
2301  "between the barotropic solution and the sum of the "//&
2302  "layer thicknesses.", units="m s-1", default=3.0e8) ! The speed of light is the default.
2303 
2304  call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&
2305  "If true, allow the adjusted velocities to have a "//&
2306  "relative CFL change up to 0.5.", default=.false.)
2307  cs%vol_CFL = cs%aggress_adjust
2308  call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &
2309  "If true, use the ratio of the open face lengths to the "//&
2310  "tracer cell areas when estimating CFL numbers. The "//&
2311  "default is set by CONT_PPM_AGGRESS_ADJUST.", &
2312  default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)
2313  call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &
2314  "The maximum CFL of the adjusted velocities.", units="nondim", &
2315  default=0.5)
2316  call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", cs%better_iter, &
2317  "If true, stop corrective iterations using a velocity "//&
2318  "based criterion and only stop if the iteration is "//&
2319  "better than all predecessors.", default=.true.)
2320  call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", &
2321  cs%use_visc_rem_max, &
2322  "If true, use more appropriate limiting bounds for "//&
2323  "corrections in strongly viscous columns.", default=.true.)
2324  call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &
2325  "If true, use the marginal face areas from the continuity "//&
2326  "solver for use as the weights in the barotropic solver. "//&
2327  "Otherwise use the transport averaged areas.", default=.true.)
2328 
2329  cs%diag => diag
2330 
2331  id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=clock_routine)
2332  id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=clock_routine)
2333 
2334 end subroutine continuity_ppm_init
2335 
2336 !> continuity_PPM_stencil returns the continuity solver stencil size
2337 function continuity_ppm_stencil(CS) result(stencil)
2338  type(continuity_ppm_cs), pointer :: cs !< Module's control structure.
2339  integer :: stencil !< The continuity solver stencil size with the current settings.
2340 
2341  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1
2342 
2343 end function continuity_ppm_stencil
2344 
2345 !> Destructor for continuity_ppm_cs
2346 subroutine continuity_ppm_end(CS)
2347  type(continuity_ppm_cs), pointer :: cs !< Module's control structure.
2348  deallocate(cs)
2349 end subroutine continuity_ppm_end
2350 
2351 !> \namespace mom_continuity_ppm
2352 !!
2353 !! This module contains the subroutines that advect layer
2354 !! thickness. The scheme here uses a Piecewise-Parabolic method with
2355 !! a positive definite limiter.
2356 
2357 end module mom_continuity_ppm
mom_continuity_ppm::loop_bounds_type
A container for loop bounds.
Definition: MOM_continuity_PPM.F90:66
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_variables::bt_cont_type
Container for information about the summed layer transports and how they will vary as the barotropic ...
Definition: MOM_variables.F90:266
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_continuity_ppm
Solve the layer continuity equation using the PPM method for layer fluxes.
Definition: MOM_continuity_PPM.F90:2
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_continuity_ppm::continuity_ppm_cs
Control structure for mom_continuity_ppm.
Definition: MOM_continuity_PPM.F90:28
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:186
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_open_boundary::obc_segment_type
Open boundary segment data structure.
Definition: MOM_open_boundary.F90:107
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:25
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239