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