MOM6
Data Types | Functions/Subroutines
mom_continuity_ppm Module Reference

Detailed Description

Solve the layer continuity equation using the PPM method for layer fluxes.

This module contains the subroutines that advect layer thickness. The scheme here uses a Piecewise-Parabolic method with a positive definite limiter.

Data Types

type  continuity_ppm_cs
 Control structure for mom_continuity_ppm. More...
 
type  loop_bounds_type
 A container for loop bounds. More...
 

Functions/Subroutines

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, based on Lin (1994). More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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 barotropic flow to agree closely with the sum of the layer's transports. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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 barotropic flow to agree closely with the sum of the layer's transports. More...
 
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. More...
 
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. More...
 
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 that is positive-definite. Here this is reinterpreted as giving a constant thickness if the mean thickness is less than h_min, with a minimum of h_min otherwise. More...
 
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 monotonic prescription of Colella and Woodward, 1984. More...
 
real function ratio_max (a, b, maxrat)
 Return the maximum ratio of a/b or maxrat. More...
 
subroutine, public continuity_ppm_init (Time, G, GV, US, param_file, diag, CS)
 Initializes continuity_ppm_cs. More...
 
integer function, public continuity_ppm_stencil (CS)
 continuity_PPM_stencil returns the continuity solver stencil size More...
 
subroutine, public continuity_ppm_end (CS)
 Destructor for continuity_ppm_cs. More...
 

Variables

integer id_clock_update
 CPU time clock IDs.
 
integer id_clock_correct
 CPU time clock IDs.
 

Function/Subroutine Documentation

◆ continuity_ppm()

subroutine, public mom_continuity_ppm::continuity_ppm ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  hin,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  uh,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)  vh,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional  uhbt,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in), optional  vhbt,
type(ocean_obc_type), optional, pointer  OBC,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in), optional  visc_rem_u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in), optional  visc_rem_v,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional  u_cor,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out), optional  v_cor,
type(bt_cont_type), optional, pointer  BT_cont 
)

Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme, based on Lin (1994).

Parameters
[in,out]gThe ocean's grid structure.
csModule's control structure.
[in]uZonal velocity [L T-1 ~> m s-1].
[in]vMeridional velocity [L T-1 ~> m s-1].
[in]hinInitial layer thickness [H ~> m or kg m-2].
[in,out]hFinal layer thickness [H ~> m or kg m-2].
[out]uhZonal volume flux, u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
[out]vhMeridional volume flux, v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].
[in]dtTime increment [T ~> s].
[in]gvVertical grid structure.
[in]usA dimensional unit scaling type
[in]uhbtThe summed volume flux through zonal faces
[in]vhbtThe summed volume flux through meridional faces
obcOpen boundaries control structure.
[in]visc_rem_uThe fraction of zonal momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]visc_rem_vThe fraction of meridional momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[out]u_corThe zonal velocities that give uhbt as the depth-integrated transport [L T-1 ~> m s-1].
[out]v_corThe meridional velocities that give vhbt as the depth-integrated transport [L T-1 ~> m s-1].
bt_contA structure with elements that describe the effective open face areas as a function of barotropic flow.

Definition at line 78 of file MOM_continuity_PPM.F90.

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 

◆ continuity_ppm_end()

subroutine, public mom_continuity_ppm::continuity_ppm_end ( type(continuity_ppm_cs), pointer  CS)

Destructor for continuity_ppm_cs.

Parameters
csModule's control structure.

Definition at line 2335 of file MOM_continuity_PPM.F90.

2335  type(continuity_PPM_CS), pointer :: CS !< Module's control structure.
2336  deallocate(cs)

◆ continuity_ppm_init()

subroutine, public mom_continuity_ppm::continuity_ppm_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(continuity_ppm_cs), pointer  CS 
)

Initializes continuity_ppm_cs.

Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure.
[in]gvVertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure indicating the open file to parse for model parameter values.
[in,out]diagA structure that is used to regulate diagnostic output.
csModule's control structure.

This include declares and sets the variable "version".

Definition at line 2230 of file MOM_continuity_PPM.F90.

2230  type(time_type), target, intent(in) :: Time !< The current model time.
2231  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2232  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
2233  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2234  type(param_file_type), intent(in) :: param_file !< A structure indicating
2235  !! the open file to parse for model parameter values.
2236  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
2237  !! regulate diagnostic output.
2238  type(continuity_PPM_CS), pointer :: CS !< Module's control structure.
2239 !> This include declares and sets the variable "version".
2240 #include "version_variable.h"
2241  real :: tol_eta_m ! An unscaled version of tol_eta [m].
2242  character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name.
2243 
2244  if (associated(cs)) then
2245  call mom_error(warning, "continuity_PPM_init called with associated control structure.")
2246  return
2247  endif
2248  allocate(cs)
2249 
2250 ! Read all relevant parameters and write them to the model log.
2251  call log_version(param_file, mdl, version, "")
2252  call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", cs%monotonic, &
2253  "If true, CONTINUITY_PPM uses the Colella and Woodward "//&
2254  "monotonic limiter. The default (false) is to use a "//&
2255  "simple positive definite limiter.", default=.false.)
2256  call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &
2257  "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2258  "(arithmetic mean) interpolation of the edge values. "//&
2259  "This may give better PV conservation properties. While "//&
2260  "it formally reduces the accuracy of the continuity "//&
2261  "solver itself in the strongly advective limit, it does "//&
2262  "not reduce the overall order of accuracy of the dynamic "//&
2263  "core.", default=.false.)
2264  call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
2265  "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
2266  "continuity solver. This scheme is highly diffusive "//&
2267  "but may be useful for debugging or in single-column "//&
2268  "mode where its minimal stencil is useful.", default=.false.)
2269  call get_param(param_file, mdl, "ETA_TOLERANCE", cs%tol_eta, &
2270  "The tolerance for the differences between the "//&
2271  "barotropic and baroclinic estimates of the sea surface "//&
2272  "height due to the fluxes through each face. The total "//&
2273  "tolerance for SSH is 4 times this value. The default "//&
2274  "is 0.5*NK*ANGSTROM, and this should not be set less "//&
2275  "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=gv%m_to_H, &
2276  default=0.5*g%ke*gv%Angstrom_m, unscaled=tol_eta_m)
2277 
2278  !### ETA_TOLERANCE_AUX can be obsoleted.
2279  call get_param(param_file, mdl, "ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
2280  "The tolerance for free-surface height discrepancies "//&
2281  "between the barotropic solution and the sum of the "//&
2282  "layer thicknesses when calculating the auxiliary "//&
2283  "corrected velocities. By default, this is the same as "//&
2284  "ETA_TOLERANCE, but can be made larger for efficiency.", &
2285  units="m", default=tol_eta_m, scale=gv%m_to_H)
2286  call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &
2287  "The tolerance for barotropic velocity discrepancies "//&
2288  "between the barotropic solution and the sum of the "//&
2289  "layer thicknesses.", units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)
2290  ! The speed of light is the default.
2291 
2292  call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&
2293  "If true, allow the adjusted velocities to have a "//&
2294  "relative CFL change up to 0.5.", default=.false.)
2295  cs%vol_CFL = cs%aggress_adjust
2296  call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &
2297  "If true, use the ratio of the open face lengths to the "//&
2298  "tracer cell areas when estimating CFL numbers. The "//&
2299  "default is set by CONT_PPM_AGGRESS_ADJUST.", &
2300  default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)
2301  call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &
2302  "The maximum CFL of the adjusted velocities.", units="nondim", &
2303  default=0.5)
2304  call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", cs%better_iter, &
2305  "If true, stop corrective iterations using a velocity "//&
2306  "based criterion and only stop if the iteration is "//&
2307  "better than all predecessors.", default=.true.)
2308  call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", &
2309  cs%use_visc_rem_max, &
2310  "If true, use more appropriate limiting bounds for "//&
2311  "corrections in strongly viscous columns.", default=.true.)
2312  call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &
2313  "If true, use the marginal face areas from the continuity "//&
2314  "solver for use as the weights in the barotropic solver. "//&
2315  "Otherwise use the transport averaged areas.", default=.true.)
2316 
2317  cs%diag => diag
2318 
2319  id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=clock_routine)
2320  id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=clock_routine)
2321 

◆ continuity_ppm_stencil()

integer function, public mom_continuity_ppm::continuity_ppm_stencil ( type(continuity_ppm_cs), pointer  CS)

continuity_PPM_stencil returns the continuity solver stencil size

Parameters
csModule's control structure.
Returns
The continuity solver stencil size with the current settings.

Definition at line 2326 of file MOM_continuity_PPM.F90.

2326  type(continuity_PPM_CS), pointer :: CS !< Module's control structure.
2327  integer :: stencil !< The continuity solver stencil size with the current settings.
2328 
2329  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1
2330 

◆ merid_face_thickness()

subroutine mom_continuity_ppm::merid_face_thickness ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_L,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_R,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)  h_v,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(loop_bounds_type), intent(in)  LB,
logical, intent(in)  vol_CFL,
logical, intent(in)  marginal,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in), optional  visc_rem_v,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Sets the effective interface thickness at each meridional velocity point.

Parameters
[in,out]gOcean's grid structure.
[in]vMeridional velocity [L T-1 ~> m s-1].
[in]hLayer thickness used to calculate fluxes, [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction, [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction, [H ~> m or kg m-2].
[in,out]h_vThickness at meridional faces, [H ~> m or kg m-2].
[in]dtTime increment [T ~> s].
[in]lbLoop bounds structure.
[in]usA dimensional unit scaling type
[in]vol_cflIf true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
[in]marginalIf true, report the marginal face thicknesses; otherwise report transport-averaged thicknesses.
[in]visc_rem_vBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
obcOpen boundaries control structure.

Definition at line 1428 of file MOM_continuity_PPM.F90.

1428  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1429  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1430  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,
1431  !! [H ~> m or kg m-2].
1432  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the reconstruction,
1433  !! [H ~> m or kg m-2].
1434  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the reconstruction,
1435  !! [H ~> m or kg m-2].
1436  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h_v !< Thickness at meridional faces,
1437  !! [H ~> m or kg m-2].
1438  real, intent(in) :: dt !< Time increment [T ~> s].
1439  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1440  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1441  logical, intent(in) :: vol_CFL !< If true, rescale the ratio
1442  !! of face areas to the cell areas when estimating the CFL number.
1443  logical, intent(in) :: marginal !< If true, report the marginal
1444  !! face thicknesses; otherwise report transport-averaged thicknesses.
1445  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(in) :: visc_rem_v !< Both the fraction
1446  !! of the momentum originally in a layer that remains after a time-step of
1447  !! viscosity, and the fraction of a time-step's worth of a barotropic
1448  !! acceleration that a layer experiences after viscosity is applied.
1449  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1450  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1451 
1452  ! Local variables
1453  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1454  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1455  ! with the same units as h [H ~> m or kg m-2] .
1456  real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].
1457  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1458  logical :: local_open_BC
1459  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1460  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1461 
1462  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
1463  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
1464  if (v(i,j,k) > 0.0) then
1465  if (vol_cfl) then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1466  else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ; endif
1467  curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
1468  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))
1469  h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + &
1470  3.0*curv_3*(cfl - 1.0))
1471  elseif (v(i,j,k) < 0.0) then
1472  if (vol_cfl) then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1473  else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ; endif
1474  curv_3 = h_l(i,j+1,k) + h_r(i,j+1,k) - 2.0*h(i,j+1,k)
1475  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))
1476  h_marg = h_l(i,j+1,k) + cfl * ((h_r(i,j+1,k)-h_l(i,j+1,k)) + &
1477  3.0*curv_3*(cfl - 1.0))
1478  else
1479  h_avg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1480  ! The choice to use the arithmetic mean here is somewhat arbitrariy, but
1481  ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
1482  h_marg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1483  ! h_marg = (2.0 * h_L(i,j+1,k) * h_R(i,j,k)) / &
1484  ! (h_L(i,j+1,k) + h_R(i,j,k) + GV%H_subroundoff)
1485  endif
1486 
1487  if (marginal) then ; h_v(i,j,k) = h_marg
1488  else ; h_v(i,j,k) = h_avg ; endif
1489  enddo ; enddo ; enddo
1490 
1491  if (present(visc_rem_v)) then
1492  !$OMP parallel do default(shared)
1493  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
1494  h_v(i,j,k) = h_v(i,j,k) * visc_rem_v(i,j,k)
1495  enddo ; enddo ; enddo
1496  endif
1497 
1498  local_open_bc = .false.
1499  if (present(obc)) then ; if (associated(obc)) then
1500  local_open_bc = obc%open_v_BCs_exist_globally
1501  endif ; endif
1502  if (local_open_bc) then
1503  do n = 1, obc%number_of_segments
1504  if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
1505  j = obc%segment(n)%HI%JsdB
1506  if (obc%segment(n)%direction == obc_direction_n) then
1507  if (present(visc_rem_v)) then ; do k=1,nz
1508  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1509  h_v(i,j,k) = h(i,j,k) * visc_rem_v(i,j,k)
1510  enddo
1511  enddo ; else ; do k=1,nz
1512  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1513  h_v(i,j,k) = h(i,j,k)
1514  enddo
1515  enddo ; endif
1516  else
1517  if (present(visc_rem_v)) then ; do k=1,nz
1518  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1519  h_v(i,j,k) = h(i,j+1,k) * visc_rem_v(i,j,k)
1520  enddo
1521  enddo ; else ; do k=1,nz
1522  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1523  h_v(i,j,k) = h(i,j+1,k)
1524  enddo
1525  enddo ; endif
1526  endif
1527  endif
1528  enddo
1529  endif
1530 

◆ merid_flux_layer()

subroutine mom_continuity_ppm::merid_flux_layer ( real, dimension(szi_(g)), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_L,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_R,
real, dimension( g %isd: g %ied), intent(inout)  vh,
real, dimension( g %isd: g %ied), intent(inout)  dvhdv,
real, dimension(szi_(g)), intent(in)  visc_rem,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  J,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, dimension( g %isd: g %ied), intent(in)  do_I,
logical, intent(in)  vol_CFL,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Evaluates the meridional mass or volume fluxes in a layer.

Parameters
[in,out]gOcean's grid structure.
[in]vMeridional velocity [L T-1 ~> m s-1].
[in]visc_remBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]hLayer thickness used to calculate fluxes, [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction [H ~> m or kg m-2].
[in,out]vhMeridional mass or volume transport [H L2 T-1 ~> m3 s-1 or kg s-1].
[in,out]dvhdvPartial derivative of vh with v [H L ~> m2 or kg m-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
[in]jSpatial index.
[in]ishStart of index range.
[in]iehEnd of index range.
[in]do_iWhich i values to work on.
[in]vol_cflIf true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
obcOpen boundaries control structure.

Definition at line 1342 of file MOM_continuity_PPM.F90.

1342  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1343  real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1344  real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the
1345  !! momentum originally in a layer that remains after a time-step
1346  !! of viscosity, and the fraction of a time-step's worth of a barotropic
1347  !! acceleration that a layer experiences after viscosity is applied.
1348  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1349  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,
1350  !! [H ~> m or kg m-2].
1351  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_L !< Left thickness in the reconstruction
1352  !! [H ~> m or kg m-2].
1353  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_R !< Right thickness in the reconstruction
1354  !! [H ~> m or kg m-2].
1355  real, dimension(SZI_(G)), intent(inout) :: vh !< Meridional mass or volume transport
1356  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
1357  real, dimension(SZI_(G)), intent(inout) :: dvhdv !< Partial derivative of vh with v
1358  !! [H L ~> m2 or kg m-1].
1359  real, intent(in) :: dt !< Time increment [T ~> s].
1360  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1361  integer, intent(in) :: j !< Spatial index.
1362  integer, intent(in) :: ish !< Start of index range.
1363  integer, intent(in) :: ieh !< End of index range.
1364  logical, dimension(SZI_(G)), intent(in) :: do_I !< Which i values to work on.
1365  logical, intent(in) :: vol_CFL !< If true, rescale the
1366  !! ratio of face areas to the cell areas when estimating the CFL number.
1367  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1368  ! Local variables
1369  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1370  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1371  ! with the same units as h, i.e. [H ~> m or kg m-2].
1372  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1373  integer :: i
1374  integer :: l_seg
1375  logical :: local_open_BC
1376 
1377  local_open_bc = .false.
1378  if (present(obc)) then ; if (associated(obc)) then
1379  local_open_bc = obc%open_v_BCs_exist_globally
1380  endif ; endif
1381 
1382  do i=ish,ieh ; if (do_i(i)) then
1383  if (v(i) > 0.0) then
1384  if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1385  else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
1386  curv_3 = h_l(i,j) + h_r(i,j) - 2.0*h(i,j)
1387  vh(i) = g%dx_Cv(i,j) * v(i) * ( h_r(i,j) + cfl * &
1388  (0.5*(h_l(i,j) - h_r(i,j)) + curv_3*(cfl - 1.5)) )
1389  h_marg = h_r(i,j) + cfl * ((h_l(i,j) - h_r(i,j)) + &
1390  3.0*curv_3*(cfl - 1.0))
1391  elseif (v(i) < 0.0) then
1392  if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1393  else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
1394  curv_3 = h_l(i,j+1) + h_r(i,j+1) - 2.0*h(i,j+1)
1395  vh(i) = g%dx_Cv(i,j) * v(i) * ( h_l(i,j+1) + cfl * &
1396  (0.5*(h_r(i,j+1)-h_l(i,j+1)) + curv_3*(cfl - 1.5)) )
1397  h_marg = h_l(i,j+1) + cfl * ((h_r(i,j+1)-h_l(i,j+1)) + &
1398  3.0*curv_3*(cfl - 1.0))
1399  else
1400  vh(i) = 0.0
1401  h_marg = 0.5 * (h_l(i,j+1) + h_r(i,j))
1402  endif
1403  dvhdv(i) = g%dx_Cv(i,j) * h_marg * visc_rem(i)
1404  endif ; enddo
1405 
1406  if (local_open_bc) then
1407  do i=ish,ieh ; if (do_i(i)) then
1408  l_seg = obc%segnum_v(i,j)
1409 
1410  if (l_seg /= obc_none) then
1411  if (obc%segment(l_seg)%open) then
1412  if (obc%segment(l_seg)%direction == obc_direction_n) then
1413  vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j)
1414  dvhdv(i) = g%dx_Cv(i,j) * h(i,j) * visc_rem(i)
1415  else
1416  vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j+1)
1417  dvhdv(i) = g%dx_Cv(i,j) * h(i,j+1) * visc_rem(i)
1418  endif
1419  endif
1420  endif
1421  endif ; enddo
1422  endif

◆ meridional_flux_adjust()

subroutine mom_continuity_ppm::meridional_flux_adjust ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_in,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_L,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_R,
real, dimension( g %isd: g %ied), intent(in), optional  vhbt,
real, dimension( g %isd: g %ied), intent(in)  vh_tot_0,
real, dimension( g %isd: g %ied), intent(in)  dvhdv_tot_0,
real, dimension( g %isd: g %ied), intent(out)  dv,
real, dimension( g %isd: g %ied), intent(in)  dv_max_CFL,
real, dimension( g %isd: g %ied), intent(in)  dv_min_CFL,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
real, dimension(szi_(g),szk_(g)), intent(in)  visc_rem,
integer, intent(in)  j,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, dimension(szi_(g)), intent(in)  do_I_in,
logical, intent(in), optional  full_precision,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout), optional  vh_3d,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.

Parameters
[in,out]gOcean's grid structure.
[in]vMeridional velocity [L T-1 ~> m s-1].
[in]h_inLayer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction [H ~> m or kg m-2].
[in]visc_remBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]vhbtThe summed volume flux through meridional faces
[in]dv_max_cflMaximum acceptable value of dv [L T-1 ~> m s-1].
[in]dv_min_cflMinimum acceptable value of dv [L T-1 ~> m s-1].
[in]vh_tot_0The summed transport with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
[in]dvhdv_tot_0The partial derivative of dv_err with dv at 0 adjustment [H L ~> m2 or kg m-1].
[out]dvThe barotropic velocity adjustment [L T-1 ~> m s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThis module's control structure.
[in]jSpatial index.
[in]ishStart of index range.
[in]iehEnd of index range.
[in]do_i_inA flag indicating which I values to work on.
[in]full_precisionA flag indicating how carefully to iterate. The default is .true. (more accurate).
[in,out]vh_3dVolume flux through meridional
obcOpen boundaries control structure.

Definition at line 1537 of file MOM_continuity_PPM.F90.

1537  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1538  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1539  intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1540  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1541  intent(in) :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
1542  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),&
1543  intent(in) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
1544  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1545  intent(in) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
1546  real, dimension(SZI_(G),SZK_(G)), intent(in) :: visc_rem
1547  !< Both the fraction of the momentum originally
1548  !! in a layer that remains after a time-step of viscosity, and the
1549  !! fraction of a time-step's worth of a barotropic acceleration that
1550  !! a layer experiences after viscosity is applied. Non-dimensional
1551  !! between 0 (at the bottom) and 1 (far above the bottom).
1552  real, dimension(SZI_(G)), &
1553  optional, intent(in) :: vhbt !< The summed volume flux through meridional faces
1554  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
1555  real, dimension(SZI_(G)), intent(in) :: dv_max_CFL !< Maximum acceptable value of dv [L T-1 ~> m s-1].
1556  real, dimension(SZI_(G)), intent(in) :: dv_min_CFL !< Minimum acceptable value of dv [L T-1 ~> m s-1].
1557  real, dimension(SZI_(G)), intent(in) :: vh_tot_0 !< The summed transport with 0 adjustment
1558  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
1559  real, dimension(SZI_(G)), intent(in) :: dvhdv_tot_0 !< The partial derivative of dv_err with
1560  !! dv at 0 adjustment [H L ~> m2 or kg m-1].
1561  real, dimension(SZI_(G)), intent(out) :: dv !< The barotropic velocity adjustment [L T-1 ~> m s-1].
1562  real, intent(in) :: dt !< Time increment [T ~> s].
1563  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1564  type(continuity_PPM_CS), pointer :: CS !< This module's control structure.
1565  integer, intent(in) :: j !< Spatial index.
1566  integer, intent(in) :: ish !< Start of index range.
1567  integer, intent(in) :: ieh !< End of index range.
1568  logical, dimension(SZI_(G)), &
1569  intent(in) :: do_I_in !< A flag indicating which I values to work on.
1570  logical, optional, intent(in) :: full_precision !< A flag indicating how carefully to
1571  !! iterate. The default is .true. (more accurate).
1572  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1573  optional, intent(inout) :: vh_3d !< Volume flux through meridional
1574  !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].
1575  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1576  ! Local variables
1577  real, dimension(SZI_(G),SZK_(G)) :: &
1578  vh_aux, & ! An auxiliary meridional volume flux [H L2 s-1 ~> m3 s-1 or kg s-1].
1579  dvhdv ! Partial derivative of vh with v [H m ~> m2 or kg m-1].
1580  real, dimension(SZI_(G)) :: &
1581  vh_err, & ! Difference between vhbt and the summed vh [H L2 T-1 ~> m3 s-1 or kg s-1].
1582  vh_err_best, & ! The smallest value of vh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
1583  v_new, & ! The velocity with the correction added [L T-1 ~> m s-1].
1584  dvhdv_tot,&! Summed partial derivative of vh with u [H L ~> m2 or kg m-1].
1585  dv_min, & ! Min/max limits on dv correction based on CFL limits
1586  dv_max ! and previous iterations [L T-1 ~> m s-1].
1587  real :: dv_prev ! The previous value of dv [L T-1 ~> m s-1].
1588  real :: ddv ! The change in dv from the previous iteration [L T-1 ~> m s-1].
1589  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
1590  real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
1591  integer :: i, k, nz, itt, max_itts = 20
1592  logical :: full_prec, domore, do_I(SZI_(G))
1593 
1594  nz = g%ke
1595  full_prec = .true. ; if (present(full_precision)) full_prec = full_precision
1596 
1597  vh_aux(:,:) = 0.0 ; dvhdv(:,:) = 0.0
1598 
1599  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
1600  vh_aux(i,k) = vh_3d(i,j,k)
1601  enddo ; enddo ; endif
1602 
1603  do i=ish,ieh
1604  dv(i) = 0.0 ; do_i(i) = do_i_in(i)
1605  dv_max(i) = dv_max_cfl(i) ; dv_min(i) = dv_min_cfl(i)
1606  vh_err(i) = vh_tot_0(i) - vhbt(i) ; dvhdv_tot(i) = dvhdv_tot_0(i)
1607  vh_err_best(i) = abs(vh_err(i))
1608  enddo
1609 
1610  do itt=1,max_itts
1611  if (full_prec) then
1612  select case (itt)
1613  case (:1) ; tol_eta = 1e-6 * cs%tol_eta
1614  case (2) ; tol_eta = 1e-4 * cs%tol_eta
1615  case (3) ; tol_eta = 1e-2 * cs%tol_eta
1616  case default ; tol_eta = cs%tol_eta
1617  end select
1618  else
1619  tol_eta = cs%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
1620  endif
1621  tol_vel = cs%tol_vel
1622 
1623  do i=ish,ieh
1624  if (vh_err(i) > 0.0) then ; dv_max(i) = dv(i)
1625  elseif (vh_err(i) < 0.0) then ; dv_min(i) = dv(i)
1626  else ; do_i(i) = .false. ; endif
1627  enddo
1628  domore = .false.
1629  do i=ish,ieh ; if (do_i(i)) then
1630  if ((dt * min(g%IareaT(i,j),g%IareaT(i,j+1))*abs(vh_err(i)) > tol_eta) .or. &
1631  (cs%better_iter .and. ((abs(vh_err(i)) > tol_vel * dvhdv_tot(i)) .or. &
1632  (abs(vh_err(i)) > vh_err_best(i))) )) then
1633  ! Use Newton's method, provided it stays bounded. Otherwise bisect
1634  ! the value with the appropriate bound.
1635  if (full_prec) then
1636  ddv = -vh_err(i) / dvhdv_tot(i)
1637  dv_prev = dv(i)
1638  dv(i) = dv(i) + ddv
1639  if (abs(ddv) < 1.0e-15*abs(dv(i))) then
1640  do_i(i) = .false. ! ddv is small enough to quit.
1641  elseif (ddv > 0.0) then
1642  if (dv(i) >= dv_max(i)) then
1643  dv(i) = 0.5*(dv_prev + dv_max(i))
1644  if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_i(i) = .false.
1645  endif
1646  else ! dvv(i) < 0.0
1647  if (dv(i) <= dv_min(i)) then
1648  dv(i) = 0.5*(dv_prev + dv_min(i))
1649  if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_i(i) = .false.
1650  endif
1651  endif
1652  else
1653  ! Use Newton's method, provided it stays bounded, just like above.
1654  dv(i) = dv(i) - vh_err(i) / dvhdv_tot(i)
1655  if ((dv(i) >= dv_max(i)) .or. (dv(i) <= dv_min(i))) &
1656  dv(i) = 0.5*(dv_max(i) + dv_min(i))
1657  endif
1658  if (do_i(i)) domore = .true.
1659  else
1660  do_i(i) = .false.
1661  endif
1662  endif ; enddo
1663  if (.not.domore) exit
1664 
1665  if ((itt < max_itts) .or. present(vh_3d)) then ; do k=1,nz
1666  do i=ish,ieh ; v_new(i) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
1667  call merid_flux_layer(v_new, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
1668  vh_aux(:,k), dvhdv(:,k), visc_rem(:,k), &
1669  dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
1670  enddo ; endif
1671 
1672  if (itt < max_itts) then
1673  do i=ish,ieh
1674  vh_err(i) = -vhbt(i) ; dvhdv_tot(i) = 0.0
1675  enddo
1676  do k=1,nz ; do i=ish,ieh
1677  vh_err(i) = vh_err(i) + vh_aux(i,k)
1678  dvhdv_tot(i) = dvhdv_tot(i) + dvhdv(i,k)
1679  enddo ; enddo
1680  do i=ish,ieh
1681  vh_err_best(i) = min(vh_err_best(i), abs(vh_err(i)))
1682  enddo
1683  endif
1684  enddo ! itt-loop
1685  ! If there are any faces which have not converged to within the tolerance,
1686  ! so-be-it, or else use a final upwind correction?
1687  ! This never seems to happen with 20 iterations as max_itt.
1688 
1689  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
1690  vh_3d(i,j,k) = vh_aux(i,k)
1691  enddo ; enddo ; endif
1692 

◆ meridional_mass_flux()

subroutine mom_continuity_ppm::meridional_mass_flux ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_in,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)  vh,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
type(loop_bounds_type), intent(in)  LB,
real, dimension(szi_(g),szjb_(g)), intent(in), optional  vhbt,
type(ocean_obc_type), optional, pointer  OBC,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in), optional  visc_rem_v,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out), optional  v_cor,
type(bt_cont_type), optional, pointer  BT_cont 
)
private

Calculates the mass or volume fluxes through the meridional faces, and other related quantities.

Parameters
[in,out]gOcean's grid structure.
[in]gvOcean's vertical grid structure.
[in]vMeridional velocity [L T-1 ~> m s-1].
[in]h_inLayer thickness used to calculate fluxes [H ~> m or kg m-2].
[out]vhVolume flux through meridional faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThis module's control structure.G
[in]lbLoop bounds structure.
obcOpen boundary condition type specifies whether, where, and what open boundary conditions are used.
[in]visc_rem_vBoth the fraction of the momentum
[in]vhbtThe summed volume flux through meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
[out]v_corThe meridional velocitiess (v with a barotropic correction) that give vhbt as the depth-integrated transport [L T-1 ~> m s-1].
bt_contA structure with elements that describe the effective open face areas as a function of barotropic flow.

Definition at line 1039 of file MOM_continuity_PPM.F90.

1039  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1040  type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
1041  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1042  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to
1043  !! calculate fluxes [H ~> m or kg m-2].
1044  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh !< Volume flux through meridional
1045  !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
1046  real, intent(in) :: dt !< Time increment [T ~> s].
1047  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1048  type(continuity_PPM_CS), pointer :: CS !< This module's control structure.G
1049  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1050  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundary condition type
1051  !! specifies whether, where, and what open boundary conditions are used.
1052  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1053  optional, intent(in) :: visc_rem_v !< Both the fraction of the momentum
1054  !! originally in a layer that remains after a time-step of viscosity,
1055  !! and the fraction of a time-step's worth of a barotropic acceleration
1056  !! that a layer experiences after viscosity is applied. Nondimensional between
1057  !! 0 (at the bottom) and 1 (far above the bottom).
1058  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt !< The summed volume flux through
1059  !< meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
1060  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1061  optional, intent(out) :: v_cor
1062  !< The meridional velocitiess (v with a barotropic correction)
1063  !! that give vhbt as the depth-integrated transport [L T-1 ~> m s-1].
1064  type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
1065  !! the effective open face areas as a function of barotropic flow.
1066  ! Local variables
1067  real, dimension(SZI_(G),SZK_(G)) :: &
1068  dvhdv ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
1069  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1070  h_L, h_R ! Left and right face thicknesses [H ~> m or kg m-2].
1071  real, dimension(SZI_(G)) :: &
1072  dv, & ! Corrective barotropic change in the velocity [L T-1 ~> m s-1].
1073  dv_min_CFL, & ! Min/max limits on dv correction
1074  dv_max_CFL, & ! to avoid CFL violations
1075  dvhdv_tot_0, & ! Summed partial derivative of vh with v [H L ~> m2 or kg m-1].
1076  vh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
1077  visc_rem_max ! The column maximum of visc_rem.
1078  logical, dimension(SZI_(G)) :: do_I
1079  real, dimension(SZI_(G)) :: FAvi ! A list of sums of meridional face areas [H L ~> m2 or kg m-1].
1080  real :: FA_v ! A sum of meridional face areas [H m ~> m2 or kg m-1].
1081  real, dimension(SZI_(G),SZK_(G)) :: &
1082  visc_rem ! A 2-D copy of visc_rem_v or an array of 1's.
1083  real :: I_vrm ! 1.0 / visc_rem_max, nondim.
1084  real :: CFL_dt ! The maximum CFL ratio of the adjusted velocities divided by
1085  ! the time step [T-1 ~> s-1].
1086  real :: I_dt ! 1.0 / dt [T-1 ~> s-1].
1087  real :: dv_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].
1088  real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [L ~> m].
1089  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1090  integer :: l_seg
1091  logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
1092  logical :: local_Flather_OBC, is_simple, local_open_BC
1093  type(OBC_segment_type), pointer :: segment => null()
1094 
1095  use_visc_rem = present(visc_rem_v)
1096  local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
1097  local_open_bc = .false.
1098  if (present(bt_cont)) set_bt_cont = (associated(bt_cont))
1099  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
1100  local_specified_bc = obc%specified_v_BCs_exist_globally
1101  local_flather_obc = obc%Flather_v_BCs_exist_globally
1102  local_open_bc = obc%open_v_BCs_exist_globally
1103  endif ; endif ; endif
1104  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1105 
1106  cfl_dt = cs%CFL_limit_adjust / dt
1107  i_dt = 1.0 / dt
1108  if (cs%aggress_adjust) cfl_dt = i_dt
1109 
1110  call cpu_clock_begin(id_clock_update)
1111 !$OMP parallel do default(none) shared(nz,ish,ieh,jsh,jeh,h_in,h_L,h_R,G,GV,LB,CS,visc_rem,OBC)
1112  do k=1,nz
1113  ! This sets h_L and h_R.
1114  if (cs%upwind_1st) then
1115  do j=jsh-1,jeh+1 ; do i=ish,ieh
1116  h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
1117  enddo ; enddo
1118  else
1119  call ppm_reconstruction_y(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
1120  2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
1121  endif
1122  do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo
1123  enddo
1124  call cpu_clock_end(id_clock_update)
1125 
1126  call cpu_clock_begin(id_clock_correct)
1127 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,v,h_in,h_L,h_R,vh,use_visc_rem, &
1128 !$OMP visc_rem_v,dt,US,G,GV,CS,local_specified_BC,OBC,vhbt, &
1129 !$OMP set_BT_cont,CFL_dt,I_dt,v_cor,BT_cont, local_Flather_OBC ) &
1130 !$OMP private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0, &
1131 !$OMP dvhdv_tot_0,visc_rem_max,I_vrm,dv_lim,dy_N, &
1132 !$OMP is_simple,FAvi,dy_S,any_simple_OBC,l_seg) &
1133 !$OMP firstprivate(visc_rem)
1134  do j=jsh-1,jeh
1135  do i=ish,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ; enddo
1136  ! This sets vh and dvhdv.
1137  do k=1,nz
1138  if (use_visc_rem) then ; do i=ish,ieh
1139  visc_rem(i,k) = visc_rem_v(i,j,k)
1140  visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
1141  enddo ; endif
1142  call merid_flux_layer(v(:,j,k), h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
1143  vh(:,j,k), dvhdv(:,k), visc_rem(:,k), &
1144  dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
1145  if (local_specified_bc) then
1146  do i=ish,ieh
1147  l_seg = obc%segnum_v(i,j)
1148 
1149  if (l_seg /= obc_none) then
1150  if (obc%segment(l_seg)%specified) &
1151  vh(i,j,k) = obc%segment(l_seg)%normal_trans(i,j,k)
1152  endif
1153  enddo
1154  endif
1155  enddo ! k-loop
1156  if ((.not.use_visc_rem) .or. (.not.cs%use_visc_rem_max)) then ; do i=ish,ieh
1157  visc_rem_max(i) = 1.0
1158  enddo ; endif
1159 
1160  if (present(vhbt) .or. set_bt_cont) then
1161  ! Set limits on dv that will keep the CFL number between -1 and 1.
1162  ! This should be adequate to keep the root bracketed in all cases.
1163  do i=ish,ieh
1164  i_vrm = 0.0
1165  if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
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  dv_max_cfl(i) = 2.0 * (cfl_dt * dy_s) * i_vrm
1171  dv_min_cfl(i) = -2.0 * (cfl_dt * dy_n) * i_vrm
1172  vh_tot_0(i) = 0.0 ; dvhdv_tot_0(i) = 0.0
1173  enddo
1174  do k=1,nz ; do i=ish,ieh
1175  dvhdv_tot_0(i) = dvhdv_tot_0(i) + dvhdv(i,k)
1176  vh_tot_0(i) = vh_tot_0(i) + vh(i,j,k)
1177  enddo ; enddo
1178 
1179  if (use_visc_rem) then
1180  if (cs%aggress_adjust) then
1181  do k=1,nz ; do i=ish,ieh
1182  if (cs%vol_CFL) then
1183  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1184  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1185  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1186  dv_lim = 0.499*((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k)))
1187  if (dv_max_cfl(i) * visc_rem(i,k) > dv_lim) &
1188  dv_max_cfl(i) = dv_lim / visc_rem(i,k)
1189 
1190  dv_lim = 0.499*((-dy_n*cfl_dt - v(i,j,k)) + max(0.0,v(i,j+1,k)))
1191  if (dv_min_cfl(i) * visc_rem(i,k) < dv_lim) &
1192  dv_min_cfl(i) = dv_lim / visc_rem(i,k)
1193  enddo ; enddo
1194  else
1195  do k=1,nz ; do i=ish,ieh
1196  if (cs%vol_CFL) then
1197  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1198  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1199  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1200  if (dv_max_cfl(i) * visc_rem(i,k) > dy_s*cfl_dt - v(i,j,k)) &
1201  dv_max_cfl(i) = (dy_s*cfl_dt - v(i,j,k)) / visc_rem(i,k)
1202  if (dv_min_cfl(i) * visc_rem(i,k) < -dy_n*cfl_dt - v(i,j,k)) &
1203  dv_min_cfl(i) = -(dy_n*cfl_dt + v(i,j,k)) / visc_rem(i,k)
1204  enddo ; enddo
1205  endif
1206  else
1207  if (cs%aggress_adjust) then
1208  do k=1,nz ; do i=ish,ieh
1209  if (cs%vol_CFL) then
1210  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1211  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1212  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1213  dv_max_cfl(i) = min(dv_max_cfl(i), 0.499 * &
1214  ((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k))) )
1215  dv_min_cfl(i) = max(dv_min_cfl(i), 0.499 * &
1216  ((-dy_n*i_dt - v(i,j,k)) + max(0.0,v(i,j+1,k))) )
1217  enddo ; enddo
1218  else
1219  do k=1,nz ; do i=ish,ieh
1220  if (cs%vol_CFL) then
1221  dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1222  dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1223  else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
1224  dv_max_cfl(i) = min(dv_max_cfl(i), dy_s*cfl_dt - v(i,j,k))
1225  dv_min_cfl(i) = max(dv_min_cfl(i), -(dy_n*cfl_dt + v(i,j,k)))
1226  enddo ; enddo
1227  endif
1228  endif
1229  do i=ish,ieh
1230  dv_max_cfl(i) = max(dv_max_cfl(i),0.0)
1231  dv_min_cfl(i) = min(dv_min_cfl(i),0.0)
1232  enddo
1233 
1234  any_simple_obc = .false.
1235  if (present(vhbt) .or. set_bt_cont) then
1236  if (local_specified_bc .or. local_flather_obc) then ; do i=ish,ieh
1237  l_seg = obc%segnum_v(i,j)
1238 
1239  ! Avoid reconciling barotropic/baroclinic transports if transport is specified
1240  is_simple = .false.
1241  if (l_seg /= obc_none) &
1242  is_simple = obc%segment(l_seg)%specified
1243  do_i(i) = .not.(l_seg /= obc_none .and. is_simple)
1244  any_simple_obc = any_simple_obc .or. is_simple
1245  enddo ; else ; do i=ish,ieh
1246  do_i(i) = .true.
1247  enddo ; endif
1248  endif
1249 
1250  if (present(vhbt)) then
1251  call meridional_flux_adjust(v, h_in, h_l, h_r, vhbt(:,j), vh_tot_0, dvhdv_tot_0, dv, &
1252  dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1253  j, ish, ieh, do_i, .true., vh, obc=obc)
1254 
1255  if (present(v_cor)) then ; do k=1,nz
1256  do i=ish,ieh ; v_cor(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
1257  if (local_specified_bc) then ; do i=ish,ieh
1258  l_seg = obc%segnum_v(i,j)
1259 
1260  if (l_seg /= obc_none) then
1261  if (obc%segment(obc%segnum_v(i,j))%specified) &
1262  v_cor(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1263  endif
1264  enddo ; endif
1265  enddo ; endif ! v-corrected
1266  endif
1267 
1268  if (set_bt_cont) then
1269  call set_merid_bt_cont(v, h_in, h_l, h_r, bt_cont, vh_tot_0, dvhdv_tot_0,&
1270  dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1271  visc_rem_max, j, ish, ieh, do_i)
1272  if (any_simple_obc) then
1273  do i=ish,ieh
1274  l_seg = obc%segnum_v(i,j)
1275 
1276  do_i(i) = .false.
1277  if(l_seg /= obc_none) &
1278  do_i(i) = (obc%segment(l_seg)%specified)
1279 
1280  if (do_i(i)) favi(i) = gv%H_subroundoff*g%dx_Cv(i,j)
1281  enddo
1282  ! NOTE: do_I(I) should prevent access to segment OBC_NONE
1283  do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then
1284  if ((abs(obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
1285  (obc%segment(obc%segnum_v(i,j))%specified)) &
1286  favi(i) = favi(i) + obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k) / &
1287  obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1288  endif ; enddo ; enddo
1289  do i=ish,ieh ; if (do_i(i)) then
1290  bt_cont%FA_v_S0(i,j) = favi(i) ; bt_cont%FA_v_N0(i,j) = favi(i)
1291  bt_cont%FA_v_SS(i,j) = favi(i) ; bt_cont%FA_v_NN(i,j) = favi(i)
1292  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1293  endif ; enddo
1294  endif
1295  endif ! set_BT_cont
1296 
1297  endif ! present(vhbt) or set_BT_cont
1298 
1299  enddo ! j-loop
1300 
1301  if (local_open_bc .and. set_bt_cont) then
1302  do n = 1, obc%number_of_segments
1303  if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
1304  j = obc%segment(n)%HI%JsdB
1305  if (obc%segment(n)%direction == obc_direction_n) then
1306  do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
1307  fa_v = 0.0
1308  do k=1,nz ; fa_v = fa_v + h_in(i,j,k)*g%dx_Cv(i,j) ; enddo
1309  bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
1310  bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
1311  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1312  enddo
1313  else
1314  do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
1315  fa_v = 0.0
1316  do k=1,nz ; fa_v = fa_v + h_in(i,j+1,k)*g%dx_Cv(i,j) ; enddo
1317  bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
1318  bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
1319  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1320  enddo
1321  endif
1322  endif
1323  enddo
1324  endif
1325  call cpu_clock_end(id_clock_correct)
1326 
1327  if (set_bt_cont) then ; if (allocated(bt_cont%h_v)) then
1328  if (present(v_cor)) then
1329  call merid_face_thickness(v_cor, h_in, h_l, h_r, bt_cont%h_v, dt, g, us, lb, &
1330  cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
1331  else
1332  call merid_face_thickness(v, h_in, h_l, h_r, bt_cont%h_v, dt, g, us, lb, &
1333  cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
1334  endif
1335  endif ; endif
1336 

◆ ppm_limit_cw84()

subroutine mom_continuity_ppm::ppm_limit_cw84 ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in,
real, dimension(szi_(g),szj_(g)), intent(inout)  h_L,
real, dimension(szi_(g),szj_(g)), intent(inout)  h_R,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  iis,
integer, intent(in)  iie,
integer, intent(in)  jis,
integer, intent(in)  jie 
)
private

This subroutine limits the left/right edge values of the PPM reconstruction according to the monotonic prescription of Colella and Woodward, 1984.

Parameters
[in]gOcean's grid structure.
[in]h_inLayer thickness [H ~> m or kg m-2].
[in,out]h_lLeft thickness in the reconstruction, [H ~> m or kg m-2].
[in,out]h_rRight thickness in the reconstruction, [H ~> m or kg m-2].
[in]iisStart of i index range.
[in]iieEnd of i index range.
[in]jisStart of j index range.
[in]jieEnd of j index range.

Definition at line 2179 of file MOM_continuity_PPM.F90.

2179  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2180  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2181  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction,
2182  !! [H ~> m or kg m-2].
2183  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction,
2184  !! [H ~> m or kg m-2].
2185  integer, intent(in) :: iis !< Start of i index range.
2186  integer, intent(in) :: iie !< End of i index range.
2187  integer, intent(in) :: jis !< Start of j index range.
2188  integer, intent(in) :: jie !< End of j index range.
2189 
2190  ! Local variables
2191  real :: h_i, RLdiff, RLdiff2, RLmean, FunFac
2192  character(len=256) :: mesg
2193  integer :: i,j
2194 
2195  do j=jis,jie ; do i=iis,iie
2196  ! This limiter monotonizes the parabola following
2197  ! Colella and Woodward, 1984, Eq. 1.10
2198  h_i = h_in(i,j)
2199  if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. ) then
2200  h_l(i,j) = h_i ; h_r(i,j) = h_i
2201  else
2202  rldiff = h_r(i,j) - h_l(i,j) ! Difference of edge values
2203  rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) ) ! Mean of edge values
2204  funfac = 6. * rldiff * ( h_i - rlmean ) ! Some funny factor
2205  rldiff2 = rldiff * rldiff ! Square of difference
2206  if ( funfac > rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)
2207  if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)
2208  endif
2209  enddo ; enddo
2210 
2211  return

◆ ppm_limit_pos()

subroutine mom_continuity_ppm::ppm_limit_pos ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in,
real, dimension(szi_(g),szj_(g)), intent(inout)  h_L,
real, dimension(szi_(g),szj_(g)), intent(inout)  h_R,
real, intent(in)  h_min,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  iis,
integer, intent(in)  iie,
integer, intent(in)  jis,
integer, intent(in)  jie 
)
private

This subroutine limits the left/right edge values of the PPM reconstruction to give a reconstruction that is positive-definite. Here this is reinterpreted as giving a constant thickness if the mean thickness is less than h_min, with a minimum of h_min otherwise.

Parameters
[in]gOcean's grid structure.
[in]h_inLayer thickness [H ~> m or kg m-2].
[in,out]h_lLeft thickness in the reconstruction [H ~> m or kg m-2].
[in,out]h_rRight thickness in the reconstruction [H ~> m or kg m-2].
[in]h_minThe minimum thickness that can be obtained by a concave parabolic fit.
[in]iisStart of i index range.
[in]iieEnd of i index range.
[in]jisStart of j index range.
[in]jieEnd of j index range.

Definition at line 2138 of file MOM_continuity_PPM.F90.

2138  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2139  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2140  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
2141  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
2142  real, intent(in) :: h_min !< The minimum thickness
2143  !! that can be obtained by a concave parabolic fit.
2144  integer, intent(in) :: iis !< Start of i index range.
2145  integer, intent(in) :: iie !< End of i index range.
2146  integer, intent(in) :: jis !< Start of j index range.
2147  integer, intent(in) :: jie !< End of j index range.
2148 
2149 ! Local variables
2150  real :: curv, dh, scale
2151  character(len=256) :: mesg
2152  integer :: i,j
2153 
2154  do j=jis,jie ; do i=iis,iie
2155  ! This limiter prevents undershooting minima within the domain with
2156  ! values less than h_min.
2157  curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2158  if (curv > 0.0) then ! Only minima are limited.
2159  dh = h_r(i,j) - h_l(i,j)
2160  if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
2161  if (h_in(i,j) <= h_min) then
2162  h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2163  elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
2164  ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
2165  ! be limited in this case. 0 < scale < 1.
2166  scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2167  h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2168  h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2169  endif
2170  endif
2171  endif
2172  enddo ; enddo
2173 

◆ ppm_reconstruction_x()

subroutine mom_continuity_ppm::ppm_reconstruction_x ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in,
real, dimension(szi_(g),szj_(g)), intent(out)  h_L,
real, dimension(szi_(g),szj_(g)), intent(out)  h_R,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB,
real, intent(in)  h_min,
logical, intent(in), optional  monotonic,
logical, intent(in), optional  simple_2nd,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Calculates left/right edge values for PPM reconstruction.

Parameters
[in]gOcean's grid structure.
[in]h_inLayer thickness [H ~> m or kg m-2].
[out]h_lLeft thickness in the reconstruction, [H ~> m or kg m-2].
[out]h_rRight thickness in the reconstruction, [H ~> m or kg m-2].
[in]lbActive loop bounds structure.
[in]h_minThe minimum thickness that can be obtained by a concave parabolic fit.
[in]monotonicIf true, use the Colella & Woodward monotonic limiter. Otherwise use a simple positive-definite limiter.
[in]simple_2ndIf true, use the arithmetic mean thicknesses as the default edge values for a simple 2nd order scheme.
obcOpen boundaries control structure.

Definition at line 1859 of file MOM_continuity_PPM.F90.

1859  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1860  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
1861  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_L !< Left thickness in the reconstruction,
1862  !! [H ~> m or kg m-2].
1863  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_R !< Right thickness in the reconstruction,
1864  !! [H ~> m or kg m-2].
1865  type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
1866  real, intent(in) :: h_min !< The minimum thickness
1867  !! that can be obtained by a concave parabolic fit.
1868  logical, optional, intent(in) :: monotonic !< If true, use the
1869  !! Colella & Woodward monotonic limiter.
1870  !! Otherwise use a simple positive-definite limiter.
1871  logical, optional, intent(in) :: simple_2nd !< If true, use the
1872  !! arithmetic mean thicknesses as the default edge values
1873  !! for a simple 2nd order scheme.
1874  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1875 
1876  ! Local variables with useful mnemonic names.
1877  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1878  real, parameter :: oneSixth = 1./6.
1879  real :: h_ip1, h_im1
1880  real :: dMx, dMn
1881  logical :: use_CW84, use_2nd
1882  character(len=256) :: mesg
1883  integer :: i, j, isl, iel, jsl, jel, n, stencil
1884  logical :: local_open_BC
1885  type(OBC_segment_type), pointer :: segment => null()
1886 
1887  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
1888  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1889 
1890  local_open_bc = .false.
1891  if (present(obc)) then ; if (associated(obc)) then
1892  local_open_bc = obc%open_u_BCs_exist_globally
1893  endif ; endif
1894 
1895  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1896 
1897  ! This is the stencil of the reconstruction, not the scheme overall.
1898  stencil = 2 ; if (use_2nd) stencil = 1
1899 
1900  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
1901  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1902  & "x-halo that needs to be increased by ",i2,".")') &
1903  stencil + max(g%isd-isl,iel-g%ied)
1904  call mom_error(fatal,mesg)
1905  endif
1906  if ((jsl < g%jsd) .or. (jel > g%jed)) then
1907  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1908  & "y-halo that needs to be increased by ",i2,".")') &
1909  max(g%jsd-jsl,jel-g%jed)
1910  call mom_error(fatal,mesg)
1911  endif
1912 
1913  if (use_2nd) then
1914  do j=jsl,jel ; do i=isl,iel
1915  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1916  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1917  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1918  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1919  enddo ; enddo
1920  else
1921  do j=jsl,jel ; do i=isl-1,iel+1
1922  if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
1923  slp(i,j) = 0.0
1924  else
1925  ! This uses a simple 2nd order slope.
1926  slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1927  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1928  dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1929  dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1930  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1931  ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
1932  endif
1933  enddo ; enddo
1934 
1935  if (local_open_bc) then
1936  do n=1, obc%number_of_segments
1937  segment => obc%segment(n)
1938  if (.not. segment%on_pe) cycle
1939  if (segment%direction == obc_direction_e .or. &
1940  segment%direction == obc_direction_w) then
1941  i=segment%HI%IsdB
1942  do j=segment%HI%jsd,segment%HI%jed
1943  slp(i+1,j) = 0.0
1944  slp(i,j) = 0.0
1945  enddo
1946  endif
1947  enddo
1948  endif
1949 
1950  do j=jsl,jel ; do i=isl,iel
1951  ! Neighboring values should take into account any boundaries. The 3
1952  ! following sets of expressions are equivalent.
1953  ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
1954  ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
1955  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1956  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1957  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
1958  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1959  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1960  enddo ; enddo
1961  endif
1962 
1963  if (local_open_bc) then
1964  do n=1, obc%number_of_segments
1965  segment => obc%segment(n)
1966  if (.not. segment%on_pe) cycle
1967  if (segment%direction == obc_direction_e) then
1968  i=segment%HI%IsdB
1969  do j=segment%HI%jsd,segment%HI%jed
1970  h_l(i+1,j) = h_in(i,j)
1971  h_r(i+1,j) = h_in(i,j)
1972  h_l(i,j) = h_in(i,j)
1973  h_r(i,j) = h_in(i,j)
1974  enddo
1975  elseif (segment%direction == obc_direction_w) then
1976  i=segment%HI%IsdB
1977  do j=segment%HI%jsd,segment%HI%jed
1978  h_l(i,j) = h_in(i+1,j)
1979  h_r(i,j) = h_in(i+1,j)
1980  h_l(i+1,j) = h_in(i+1,j)
1981  h_r(i+1,j) = h_in(i+1,j)
1982  enddo
1983  endif
1984  enddo
1985  endif
1986 
1987  if (use_cw84) then
1988  call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
1989  else
1990  call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
1991  endif
1992 
1993  return

◆ ppm_reconstruction_y()

subroutine mom_continuity_ppm::ppm_reconstruction_y ( real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_in,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  h_L,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  h_R,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB,
real, intent(in)  h_min,
logical, intent(in), optional  monotonic,
logical, intent(in), optional  simple_2nd,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Calculates left/right edge values for PPM reconstruction.

Parameters
[in]gOcean's grid structure.
[in]h_inLayer thickness [H ~> m or kg m-2].
[out]h_lLeft thickness in the reconstruction, [H ~> m or kg m-2].
[out]h_rRight thickness in the reconstruction, [H ~> m or kg m-2].
[in]lbActive loop bounds structure.
[in]h_minThe minimum thickness that can be obtained by a concave parabolic fit.
[in]monotonicIf true, use the Colella & Woodward monotonic limiter. Otherwise use a simple positive-definite limiter.
[in]simple_2ndIf true, use the arithmetic mean thicknesses as the default edge values for a simple 2nd order scheme.
obcOpen boundaries control structure.

Definition at line 1998 of file MOM_continuity_PPM.F90.

1998  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1999  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2000  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_L !< Left thickness in the reconstruction,
2001  !! [H ~> m or kg m-2].
2002  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_R !< Right thickness in the reconstruction,
2003  !! [H ~> m or kg m-2].
2004  type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
2005  real, intent(in) :: h_min !< The minimum thickness
2006  !! that can be obtained by a concave parabolic fit.
2007  logical, optional, intent(in) :: monotonic !< If true, use the
2008  !! Colella & Woodward monotonic limiter.
2009  !! Otherwise use a simple positive-definite limiter.
2010  logical, optional, intent(in) :: simple_2nd !< If true, use the
2011  !! arithmetic mean thicknesses as the default edge values
2012  !! for a simple 2nd order scheme.
2013  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
2014 
2015  ! Local variables with useful mnemonic names.
2016  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
2017  real, parameter :: oneSixth = 1./6.
2018  real :: h_jp1, h_jm1
2019  real :: dMx, dMn
2020  logical :: use_CW84, use_2nd
2021  character(len=256) :: mesg
2022  integer :: i, j, isl, iel, jsl, jel, n, stencil
2023  logical :: local_open_BC
2024  type(OBC_segment_type), pointer :: segment => null()
2025 
2026  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
2027  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
2028 
2029  local_open_bc = .false.
2030  if (present(obc)) then ; if (associated(obc)) then
2031  local_open_bc = obc%open_v_BCs_exist_globally
2032  endif ; endif
2033 
2034  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
2035 
2036  ! This is the stencil of the reconstruction, not the scheme overall.
2037  stencil = 2 ; if (use_2nd) stencil = 1
2038 
2039  if ((isl < g%isd) .or. (iel > g%ied)) then
2040  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
2041  & "x-halo that needs to be increased by ",i2,".")') &
2042  max(g%isd-isl,iel-g%ied)
2043  call mom_error(fatal,mesg)
2044  endif
2045  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
2046  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
2047  & "y-halo that needs to be increased by ",i2,".")') &
2048  stencil + max(g%jsd-jsl,jel-g%jed)
2049  call mom_error(fatal,mesg)
2050  endif
2051 
2052  if (use_2nd) then
2053  do j=jsl,jel ; do i=isl,iel
2054  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2055  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2056  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
2057  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
2058  enddo ; enddo
2059  else
2060  do j=jsl-1,jel+1 ; do i=isl,iel
2061  if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
2062  slp(i,j) = 0.0
2063  else
2064  ! This uses a simple 2nd order slope.
2065  slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
2066  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
2067  dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
2068  dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
2069  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
2070  ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
2071  endif
2072  enddo ; enddo
2073 
2074  if (local_open_bc) then
2075  do n=1, obc%number_of_segments
2076  segment => obc%segment(n)
2077  if (.not. segment%on_pe) cycle
2078  if (segment%direction == obc_direction_s .or. &
2079  segment%direction == obc_direction_n) then
2080  j=segment%HI%JsdB
2081  do i=segment%HI%isd,segment%HI%ied
2082  slp(i,j+1) = 0.0
2083  slp(i,j) = 0.0
2084  enddo
2085  endif
2086  enddo
2087  endif
2088 
2089  do j=jsl,jel ; do i=isl,iel
2090  ! Neighboring values should take into account any boundaries. The 3
2091  ! following sets of expressions are equivalent.
2092  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2093  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2094  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2095  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2096  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2097  enddo ; enddo
2098  endif
2099 
2100  if (local_open_bc) then
2101  do n=1, obc%number_of_segments
2102  segment => obc%segment(n)
2103  if (.not. segment%on_pe) cycle
2104  if (segment%direction == obc_direction_n) then
2105  j=segment%HI%JsdB
2106  do i=segment%HI%isd,segment%HI%ied
2107  h_l(i,j+1) = h_in(i,j)
2108  h_r(i,j+1) = h_in(i,j)
2109  h_l(i,j) = h_in(i,j)
2110  h_r(i,j) = h_in(i,j)
2111  enddo
2112  elseif (segment%direction == obc_direction_s) then
2113  j=segment%HI%JsdB
2114  do i=segment%HI%isd,segment%HI%ied
2115  h_l(i,j) = h_in(i,j+1)
2116  h_r(i,j) = h_in(i,j+1)
2117  h_l(i,j+1) = h_in(i,j+1)
2118  h_r(i,j+1) = h_in(i,j+1)
2119  enddo
2120  endif
2121  enddo
2122  endif
2123 
2124  if (use_cw84) then
2125  call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
2126  else
2127  call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
2128  endif
2129 
2130  return

◆ ratio_max()

real function mom_continuity_ppm::ratio_max ( real, intent(in)  a,
real, intent(in)  b,
real, intent(in)  maxrat 
)
private

Return the maximum ratio of a/b or maxrat.

Parameters
[in]aNumerator
[in]bDenominator
[in]maxratMaximum value of ratio.
Returns
Return value.

Definition at line 2216 of file MOM_continuity_PPM.F90.

2216  real, intent(in) :: a !< Numerator
2217  real, intent(in) :: b !< Denominator
2218  real, intent(in) :: maxrat !< Maximum value of ratio.
2219  real :: ratio !< Return value.
2220 
2221  if (abs(a) > abs(maxrat*b)) then
2222  ratio = maxrat
2223  else
2224  ratio = a / b
2225  endif

◆ set_merid_bt_cont()

subroutine mom_continuity_ppm::set_merid_bt_cont ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_in,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_L,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_R,
type(bt_cont_type), intent(inout)  BT_cont,
real, dimension(szi_(g)), intent(in)  vh_tot_0,
real, dimension(szi_(g)), intent(in)  dvhdv_tot_0,
real, dimension(szi_(g)), intent(in)  dv_max_CFL,
real, dimension(szi_(g)), intent(in)  dv_min_CFL,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %ke), intent(in)  visc_rem,
real, dimension(szi_(g)), intent(in)  visc_rem_max,
integer, intent(in)  j,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, dimension(szi_(g)), intent(in)  do_I 
)
private

Sets of a structure that describes the meridional barotropic volume or mass fluxes as a function of barotropic flow to agree closely with the sum of the layer's transports.

Parameters
[in,out]gOcean's grid structure.
[in]vMeridional velocity [L T-1 ~> m s-1].
[in]h_inLayer thickness used to calculate fluxes, [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction, [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction, [H ~> m or kg m-2].
[in,out]bt_contA structure with elements that describe the effective open face areas as a function of barotropic flow.
[in]vh_tot_0The summed transport with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
[in]dvhdv_tot_0The partial derivative of du_err with dv at 0 adjustment [H L ~> m2 or kg m-1].
[in]dv_max_cflMaximum acceptable value of dv [L T-1 ~> m s-1].
[in]dv_min_cflMinimum acceptable value of dv [L T-1 ~> m s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThis module's control structure.
[in]visc_remBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]visc_rem_maxMaximum allowable visc_rem.
[in]jSpatial index.
[in]ishStart of index range.
[in]iehEnd of index range.
[in]do_iA logical flag indicating which I values to work on.

Definition at line 1700 of file MOM_continuity_PPM.F90.

1700  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1701  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1702  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to calculate fluxes,
1703  !! [H ~> m or kg m-2].
1704  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the reconstruction,
1705  !! [H ~> m or kg m-2].
1706  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the reconstruction,
1707  !! [H ~> m or kg m-2].
1708  type(BT_cont_type), intent(inout) :: BT_cont !< A structure with elements
1709  !! that describe the effective open face areas as a function of barotropic flow.
1710  real, dimension(SZI_(G)), intent(in) :: vh_tot_0 !< The summed transport
1711  !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
1712  real, dimension(SZI_(G)), intent(in) :: dvhdv_tot_0 !< The partial derivative
1713  !! of du_err with dv at 0 adjustment [H L ~> m2 or kg m-1].
1714  real, dimension(SZI_(G)), intent(in) :: dv_max_CFL !< Maximum acceptable value
1715  !! of dv [L T-1 ~> m s-1].
1716  real, dimension(SZI_(G)), intent(in) :: dv_min_CFL !< Minimum acceptable value
1717  !! of dv [L T-1 ~> m s-1].
1718  real, intent(in) :: dt !< Time increment [T ~> s].
1719  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1720  type(continuity_PPM_CS), pointer :: CS !< This module's control structure.
1721  real, dimension(SZI_(G),SZK_(G)), intent(in) :: visc_rem !< Both the fraction of the
1722  !! momentum originally in a layer that remains after a time-step
1723  !! of viscosity, and the fraction of a time-step's worth of a barotropic
1724  !! acceleration that a layer experiences after viscosity is applied.
1725  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1726  real, dimension(SZI_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem.
1727  integer, intent(in) :: j !< Spatial index.
1728  integer, intent(in) :: ish !< Start of index range.
1729  integer, intent(in) :: ieh !< End of index range.
1730  logical, dimension(SZI_(G)), intent(in) :: do_I !< A logical flag indicating
1731  !! which I values to work on.
1732  ! Local variables
1733  real, dimension(SZI_(G)) :: &
1734  dv0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
1735  dvL, dvR, & ! The barotropic velocity increments that give the southerly
1736  ! (dvL) and northerly (dvR) test velocities [L T-1 ~> m s-1].
1737  zeros, & ! An array of full of 0's.
1738  dv_cfl, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
1739  v_l, v_r, & ! The southerly (v_L), northerly (v_R), and zero-barotropic
1740  v_0, & ! transport (v_0) layer test velocities [L T-1 ~> m s-1].
1741  dvhdv_l, & ! The effective layer marginal face areas with the southerly
1742  dvhdv_r, & ! (_L), northerly (_R), and zero-barotropic (_0) test
1743  dvhdv_0, & ! velocities [H L ~> m2 or kg m-1].
1744  vh_l, vh_r, & ! The layer transports with the southerly (_L), northerly (_R)
1745  vh_0, & ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
1746  famt_l, famt_r, & ! The summed effective marginal face areas for the 3
1747  famt_0, & ! test velocities [H m ~> m2 or kg m-1].
1748  vhtot_l, & ! The summed transport with the southerly (vhtot_L) and
1749  vhtot_r ! and northerly (vhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
1750  real :: FA_0 ! The effective face area with 0 barotropic transport [H L ~> m2 or kg m-1].
1751  real :: FA_avg ! The average effective face area [H L ~> m2 or kg m-1], nominally given by
1752  ! the realized transport divided by the barotropic velocity.
1753  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim] This
1754  ! limiting is necessary to keep the inverse of visc_rem
1755  ! from leading to large CFL numbers.
1756  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
1757  ! in finding the barotropic velocity that changes the
1758  ! flow direction. This is necessary to keep the inverse
1759  ! of visc_rem from leading to large CFL numbers.
1760  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
1761  ! flow is truly upwind [nondim]
1762  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
1763  logical :: domore
1764  integer :: i, k, nz
1765 
1766  nz = g%ke ; idt = 1.0 / dt
1767  min_visc_rem = 0.1 ; cfl_min = 1e-6
1768 
1769  ! Diagnose the zero-transport correction, dv0.
1770  do i=ish,ieh ; zeros(i) = 0.0 ; enddo
1771  call meridional_flux_adjust(v, h_in, h_l, h_r, zeros, vh_tot_0, dvhdv_tot_0, dv0, &
1772  dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1773  j, ish, ieh, do_i, .true.)
1774 
1775  ! Determine the southerly- and northerly- fluxes. Choose a sufficiently
1776  ! negative velocity correction for the northerly-flux, and a sufficiently
1777  ! positive correction for the southerly-flux.
1778  domore = .false.
1779  do i=ish,ieh ; if (do_i(i)) then
1780  domore = .true.
1781  dv_cfl(i) = (cfl_min * idt) * g%dyCv(i,j)
1782  dvr(i) = min(0.0,dv0(i) - dv_cfl(i))
1783  dvl(i) = max(0.0,dv0(i) + dv_cfl(i))
1784  famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
1785  vhtot_l(i) = 0.0 ; vhtot_r(i) = 0.0
1786  endif ; enddo
1787 
1788  if (.not.domore) then
1789  do k=1,nz ; do i=ish,ieh
1790  bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
1791  bt_cont%vBT_SS(i,j) = 0.0
1792  bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
1793  bt_cont%vBT_NN(i,j) = 0.0
1794  enddo ; enddo
1795  return
1796  endif
1797 
1798  do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then
1799  visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
1800  if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
1801  if (v(i,j,k) + dvr(i)*visc_rem_lim > -dv_cfl(i)*visc_rem(i,k)) &
1802  dvr(i) = -(v(i,j,k) + dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
1803  if (v(i,j,k) + dvl(i)*visc_rem_lim < dv_cfl(i)*visc_rem(i,k)) &
1804  dvl(i) = -(v(i,j,k) - dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
1805  endif
1806  endif ; enddo ; enddo
1807  do k=1,nz
1808  do i=ish,ieh ; if (do_i(i)) then
1809  v_l(i) = v(i,j,k) + dvl(i) * visc_rem(i,k)
1810  v_r(i) = v(i,j,k) + dvr(i) * visc_rem(i,k)
1811  v_0(i) = v(i,j,k) + dv0(i) * visc_rem(i,k)
1812  endif ; enddo
1813  call merid_flux_layer(v_0, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_0, dvhdv_0, &
1814  visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
1815  call merid_flux_layer(v_l, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_l, dvhdv_l, &
1816  visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
1817  call merid_flux_layer(v_r, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_r, dvhdv_r, &
1818  visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
1819  do i=ish,ieh ; if (do_i(i)) then
1820  famt_0(i) = famt_0(i) + dvhdv_0(i)
1821  famt_l(i) = famt_l(i) + dvhdv_l(i)
1822  famt_r(i) = famt_r(i) + dvhdv_r(i)
1823  vhtot_l(i) = vhtot_l(i) + vh_l(i)
1824  vhtot_r(i) = vhtot_r(i) + vh_r(i)
1825  endif ; enddo
1826  enddo
1827  do i=ish,ieh ; if (do_i(i)) then
1828  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1829  if ((dvl(i) - dv0(i)) /= 0.0) &
1830  fa_avg = vhtot_l(i) / (dvl(i) - dv0(i))
1831  if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))
1832  elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif
1833  bt_cont%FA_v_S0(i,j) = fa_0 ; bt_cont%FA_v_SS(i,j) = famt_l(i)
1834  if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%vBT_SS(i,j) = 0.0 ; else
1835  bt_cont%vBT_SS(i,j) = (1.5 * (dvl(i) - dv0(i))) * &
1836  ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
1837  endif
1838 
1839  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1840  if ((dvr(i) - dv0(i)) /= 0.0) &
1841  fa_avg = vhtot_r(i) / (dvr(i) - dv0(i))
1842  if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))
1843  elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif
1844  bt_cont%FA_v_N0(i,j) = fa_0 ; bt_cont%FA_v_NN(i,j) = famt_r(i)
1845  if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%vBT_NN(i,j) = 0.0 ; else
1846  bt_cont%vBT_NN(i,j) = (1.5 * (dvr(i) - dv0(i))) * &
1847  ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
1848  endif
1849  else
1850  bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
1851  bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
1852  bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1853  endif ; enddo
1854 

◆ set_zonal_bt_cont()

subroutine mom_continuity_ppm::set_zonal_bt_cont ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_in,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_L,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_R,
type(bt_cont_type), intent(inout)  BT_cont,
real, dimension( g %isdb: g %iedb), intent(in)  uh_tot_0,
real, dimension( g %isdb: g %iedb), intent(in)  duhdu_tot_0,
real, dimension( g %isdb: g %iedb), intent(in)  du_max_CFL,
real, dimension( g %isdb: g %iedb), intent(in)  du_min_CFL,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
real, dimension(szib_(g),szk_(g)), intent(in)  visc_rem,
real, dimension( g %isdb: g %iedb), intent(in)  visc_rem_max,
integer, intent(in)  j,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, dimension( g %isdb: g %iedb), intent(in)  do_I 
)
private

Sets a structure that describes the zonal barotropic volume or mass fluxes as a function of barotropic flow to agree closely with the sum of the layer's transports.

Parameters
[in,out]gOcean's grid structure.
[in]uZonal velocity [L T-1 ~> m s-1].
[in]h_inLayer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction [H ~> m or kg m-2].
[in,out]bt_contA structure with elements that describe the effective open face areas as a function of barotropic flow.
[in]uh_tot_0The summed transport with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
[in]duhdu_tot_0The partial derivative of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
[in]du_max_cflMaximum acceptable value of du [L T-1 ~> m s-1].
[in]du_min_cflMinimum acceptable value of du [L T-1 ~> m s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThis module's control structure.
[in]visc_remBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]visc_rem_maxMaximum allowable visc_rem.
[in]jSpatial index.
[in]ishStart of index range.
[in]iehEnd of index range.
[in]do_iA logical flag indicating which I values to work on.

Definition at line 877 of file MOM_continuity_PPM.F90.

877  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
878  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
879  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to
880  !! calculate fluxes [H ~> m or kg m-2].
881  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the
882  !! reconstruction [H ~> m or kg m-2].
883  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the
884  !! reconstruction [H ~> m or kg m-2].
885  type(BT_cont_type), intent(inout) :: BT_cont !< A structure with elements
886  !! that describe the effective open face areas as a function of barotropic flow.
887  real, dimension(SZIB_(G)), intent(in) :: uh_tot_0 !< The summed transport
888  !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
889  real, dimension(SZIB_(G)), intent(in) :: duhdu_tot_0 !< The partial derivative
890  !! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
891  real, dimension(SZIB_(G)), intent(in) :: du_max_CFL !< Maximum acceptable
892  !! value of du [L T-1 ~> m s-1].
893  real, dimension(SZIB_(G)), intent(in) :: du_min_CFL !< Minimum acceptable
894  !! value of du [L T-1 ~> m s-1].
895  real, intent(in) :: dt !< Time increment [T ~> s].
896  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
897  type(continuity_PPM_CS), pointer :: CS !< This module's control structure.
898  real, dimension(SZIB_(G),SZK_(G)), intent(in) :: visc_rem !< Both the fraction of the
899  !! momentum originally in a layer that remains after a time-step of viscosity, and
900  !! the fraction of a time-step's worth of a barotropic acceleration that a layer
901  !! experiences after viscosity is applied.
902  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
903  real, dimension(SZIB_(G)), intent(in) :: visc_rem_max !< Maximum allowable visc_rem.
904  integer, intent(in) :: j !< Spatial index.
905  integer, intent(in) :: ish !< Start of index range.
906  integer, intent(in) :: ieh !< End of index range.
907  logical, dimension(SZIB_(G)), intent(in) :: do_I !< A logical flag indicating
908  !! which I values to work on.
909  ! Local variables
910  real, dimension(SZIB_(G)) :: &
911  du0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
912  duL, duR, & ! The barotropic velocity increments that give the westerly
913  ! (duL) and easterly (duR) test velocities [L T-1 ~> m s-1].
914  zeros, & ! An array of full of 0's.
915  du_cfl, & ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
916  u_l, u_r, & ! The westerly (u_L), easterly (u_R), and zero-barotropic
917  u_0, & ! transport (u_0) layer test velocities [L T-1 ~> m s-1].
918  duhdu_l, & ! The effective layer marginal face areas with the westerly
919  duhdu_r, & ! (_L), easterly (_R), and zero-barotropic (_0) test
920  duhdu_0, & ! velocities [H L ~> m2 or kg m-1].
921  uh_l, uh_r, & ! The layer transports with the westerly (_L), easterly (_R),
922  uh_0, & ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
923  famt_l, famt_r, & ! The summed effective marginal face areas for the 3
924  famt_0, & ! test velocities [H L ~> m2 or kg m-1].
925  uhtot_l, & ! The summed transport with the westerly (uhtot_L) and
926  uhtot_r ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
927  real :: FA_0 ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m].
928  real :: FA_avg ! The average effective face area [L H ~> m2 or kg m], nominally given by
929  ! the realized transport divided by the barotropic velocity.
930  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim] This
931  ! limiting is necessary to keep the inverse of visc_rem
932  ! from leading to large CFL numbers.
933  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
934  ! in finding the barotropic velocity that changes the
935  ! flow direction. This is necessary to keep the inverse
936  ! of visc_rem from leading to large CFL numbers.
937  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
938  ! flow is truly upwind [nondim]
939  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
940  logical :: domore
941  integer :: i, k, nz
942 
943  nz = g%ke ; idt = 1.0 / dt
944  min_visc_rem = 0.1 ; cfl_min = 1e-6
945 
946  ! Diagnose the zero-transport correction, du0.
947  do i=ish-1,ieh ; zeros(i) = 0.0 ; enddo
948  call zonal_flux_adjust(u, h_in, h_l, h_r, zeros, uh_tot_0, duhdu_tot_0, du0, &
949  du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
950  j, ish, ieh, do_i, .true.)
951 
952  ! Determine the westerly- and easterly- fluxes. Choose a sufficiently
953  ! negative velocity correction for the easterly-flux, and a sufficiently
954  ! positive correction for the westerly-flux.
955  domore = .false.
956  do i=ish-1,ieh
957  if (do_i(i)) domore = .true.
958  du_cfl(i) = (cfl_min * idt) * g%dxCu(i,j)
959  dur(i) = min(0.0,du0(i) - du_cfl(i))
960  dul(i) = max(0.0,du0(i) + du_cfl(i))
961  famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
962  uhtot_l(i) = 0.0 ; uhtot_r(i) = 0.0
963  enddo
964 
965  if (.not.domore) then
966  do k=1,nz ; do i=ish-1,ieh
967  bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
968  bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
969  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
970  enddo ; enddo
971  return
972  endif
973 
974  do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
975  visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
976  if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
977  if (u(i,j,k) + dur(i)*visc_rem_lim > -du_cfl(i)*visc_rem(i,k)) &
978  dur(i) = -(u(i,j,k) + du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
979  if (u(i,j,k) + dul(i)*visc_rem_lim < du_cfl(i)*visc_rem(i,k)) &
980  dul(i) = -(u(i,j,k) - du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
981  endif
982  endif ; enddo ; enddo
983 
984  do k=1,nz
985  do i=ish-1,ieh ; if (do_i(i)) then
986  u_l(i) = u(i,j,k) + dul(i) * visc_rem(i,k)
987  u_r(i) = u(i,j,k) + dur(i) * visc_rem(i,k)
988  u_0(i) = u(i,j,k) + du0(i) * visc_rem(i,k)
989  endif ; enddo
990  call zonal_flux_layer(u_0, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_0, duhdu_0, &
991  visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
992  call zonal_flux_layer(u_l, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_l, duhdu_l, &
993  visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
994  call zonal_flux_layer(u_r, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_r, duhdu_r, &
995  visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
996  do i=ish-1,ieh ; if (do_i(i)) then
997  famt_0(i) = famt_0(i) + duhdu_0(i)
998  famt_l(i) = famt_l(i) + duhdu_l(i)
999  famt_r(i) = famt_r(i) + duhdu_r(i)
1000  uhtot_l(i) = uhtot_l(i) + uh_l(i)
1001  uhtot_r(i) = uhtot_r(i) + uh_r(i)
1002  endif ; enddo
1003  enddo
1004  do i=ish-1,ieh ; if (do_i(i)) then
1005  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1006  if ((dul(i) - du0(i)) /= 0.0) &
1007  fa_avg = uhtot_l(i) / (dul(i) - du0(i))
1008  if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))
1009  elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif
1010 
1011  bt_cont%FA_u_W0(i,j) = fa_0 ; bt_cont%FA_u_WW(i,j) = famt_l(i)
1012  if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%uBT_WW(i,j) = 0.0 ; else
1013  bt_cont%uBT_WW(i,j) = (1.5 * (dul(i) - du0(i))) * &
1014  ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
1015  endif
1016 
1017  fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1018  if ((dur(i) - du0(i)) /= 0.0) &
1019  fa_avg = uhtot_r(i) / (dur(i) - du0(i))
1020  if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))
1021  elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif
1022 
1023  bt_cont%FA_u_E0(i,j) = fa_0 ; bt_cont%FA_u_EE(i,j) = famt_r(i)
1024  if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%uBT_EE(i,j) = 0.0 ; else
1025  bt_cont%uBT_EE(i,j) = (1.5 * (dur(i) - du0(i))) * &
1026  ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
1027  endif
1028  else
1029  bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
1030  bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
1031  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
1032  endif ; enddo
1033 

◆ zonal_face_thickness()

subroutine mom_continuity_ppm::zonal_face_thickness ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_L,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_R,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  h_u,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(loop_bounds_type), intent(in)  LB,
logical, intent(in)  vol_CFL,
logical, intent(in)  marginal,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in), optional  visc_rem_u,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Sets the effective interface thickness at each zonal velocity point.

Parameters
[in,out]gOcean's grid structure.
[in]uZonal velocity [L T-1 ~> m s-1].
[in]hLayer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction [H ~> m or kg m-2].
[in,out]h_uThickness at zonal faces [H ~> m or kg m-2].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
[in]lbLoop bounds structure.
[in]vol_cflIf true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
[in]marginalIf true, report the marginal face thicknesses; otherwise report transport-averaged thicknesses.
[in]visc_rem_uBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
obcOpen boundaries control structure.

Definition at line 605 of file MOM_continuity_PPM.F90.

605  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
606  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
607  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness used to
608  !! calculate fluxes [H ~> m or kg m-2].
609  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the
610  !! reconstruction [H ~> m or kg m-2].
611  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the
612  !! reconstruction [H ~> m or kg m-2].
613  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_u !< Thickness at zonal faces [H ~> m or kg m-2].
614  real, intent(in) :: dt !< Time increment [T ~> s].
615  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
616  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
617  logical, intent(in) :: vol_CFL !< If true, rescale the ratio
618  !! of face areas to the cell areas when estimating the CFL number.
619  logical, intent(in) :: marginal !< If true, report the
620  !! marginal face thicknesses; otherwise report transport-averaged thicknesses.
621  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
622  optional, intent(in) :: visc_rem_u
623  !< Both the fraction of the momentum originally in a layer that remains after
624  !! a time-step of viscosity, and the fraction of a time-step's worth of a
625  !! barotropic acceleration that a layer experiences after viscosity is applied.
626  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
627  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
628 
629  ! Local variables
630  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
631  real :: curv_3 ! A measure of the thickness curvature over a grid length,
632  ! with the same units as h_in.
633  real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].
634  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
635  logical :: local_open_BC
636  integer :: i, j, k, ish, ieh, jsh, jeh, nz, n
637  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
638 
639  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
640  do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
641  if (u(i,j,k) > 0.0) then
642  if (vol_cfl) then ; cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
643  else ; cfl = u(i,j,k) * dt * g%IdxT(i,j) ; endif
644  curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
645  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))
646  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))
647  elseif (u(i,j,k) < 0.0) then
648  if (vol_cfl) then ; cfl = (-u(i,j,k)*dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
649  else ; cfl = -u(i,j,k) * dt * g%IdxT(i+1,j) ; endif
650  curv_3 = h_l(i+1,j,k) + h_r(i+1,j,k) - 2.0*h(i+1,j,k)
651  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))
652  h_marg = h_l(i+1,j,k) + cfl * ((h_r(i+1,j,k)-h_l(i+1,j,k)) + &
653  3.0*curv_3*(cfl - 1.0))
654  else
655  h_avg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
656  ! The choice to use the arithmetic mean here is somewhat arbitrariy, but
657  ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
658  h_marg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
659  ! h_marg = (2.0 * h_L(i+1,j,k) * h_R(i,j,k)) / &
660  ! (h_L(i+1,j,k) + h_R(i,j,k) + GV%H_subroundoff)
661  endif
662 
663  if (marginal) then ; h_u(i,j,k) = h_marg
664  else ; h_u(i,j,k) = h_avg ; endif
665  enddo ; enddo ; enddo
666  if (present(visc_rem_u)) then
667  !$OMP parallel do default(shared)
668  do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
669  h_u(i,j,k) = h_u(i,j,k) * visc_rem_u(i,j,k)
670  enddo ; enddo ; enddo
671  endif
672 
673  local_open_bc = .false.
674  if (present(obc)) then ; if (associated(obc)) then
675  local_open_bc = obc%open_u_BCs_exist_globally
676  endif ; endif
677  if (local_open_bc) then
678  do n = 1, obc%number_of_segments
679  if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
680  i = obc%segment(n)%HI%IsdB
681  if (obc%segment(n)%direction == obc_direction_e) then
682  if (present(visc_rem_u)) then ; do k=1,nz
683  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
684  h_u(i,j,k) = h(i,j,k) * visc_rem_u(i,j,k)
685  enddo
686  enddo ; else ; do k=1,nz
687  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
688  h_u(i,j,k) = h(i,j,k)
689  enddo
690  enddo ; endif
691  else
692  if (present(visc_rem_u)) then ; do k=1,nz
693  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
694  h_u(i,j,k) = h(i+1,j,k) * visc_rem_u(i,j,k)
695  enddo
696  enddo ; else ; do k=1,nz
697  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
698  h_u(i,j,k) = h(i+1,j,k)
699  enddo
700  enddo ; endif
701  endif
702  endif
703  enddo
704  endif
705 

◆ zonal_flux_adjust()

subroutine mom_continuity_ppm::zonal_flux_adjust ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_in,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_L,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_R,
real, dimension(szib_(g)), intent(in), optional  uhbt,
real, dimension(szib_(g)), intent(in)  uh_tot_0,
real, dimension(szib_(g)), intent(in)  duhdu_tot_0,
real, dimension(szib_(g)), intent(out)  du,
real, dimension(szib_(g)), intent(in)  du_max_CFL,
real, dimension(szib_(g)), intent(in)  du_min_CFL,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
real, dimension( g %isdb: g %iedb, g %ke), intent(in)  visc_rem,
integer, intent(in)  j,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, dimension( g %isdb: g %iedb), intent(in)  do_I_in,
logical, intent(in), optional  full_precision,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout), optional  uh_3d,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.

Parameters
[in,out]gOcean's grid structure.
[in]uZonal velocity [L T-1 ~> m s-1].
[in]h_inLayer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]h_lLeft thickness in the reconstruction [H ~> m or kg m-2].
[in]h_rRight thickness in the reconstruction [H ~> m or kg m-2].
[in]visc_remBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]uhbtThe summed volume flux through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
[in]du_max_cflMaximum acceptable value of du [L T-1 ~> m s-1].
[in]du_min_cflMinimum acceptable value of du [L T-1 ~> m s-1].
[in]uh_tot_0The summed transport with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
[in]duhdu_tot_0The partial derivative of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
[out]duThe barotropic velocity adjustment [L T-1 ~> m s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThis module's control structure.
[in]jSpatial index.
[in]ishStart of index range.
[in]iehEnd of index range.
[in]do_i_inA logical flag indicating which I values to work on.
[in]full_precisionA flag indicating how carefully to iterate. The default is .true. (more accurate).
[in,out]uh_3dVolume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
obcOpen boundaries control structure.

Definition at line 713 of file MOM_continuity_PPM.F90.

713  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
714  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
715  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to
716  !! calculate fluxes [H ~> m or kg m-2].
717  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the
718  !! reconstruction [H ~> m or kg m-2].
719  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the
720  !! reconstruction [H ~> m or kg m-2].
721  real, dimension(SZIB_(G),SZK_(G)), intent(in) :: visc_rem !< Both the fraction of the
722  !! momentum originally in a layer that remains after a time-step of viscosity, and
723  !! the fraction of a time-step's worth of a barotropic acceleration that a layer
724  !! experiences after viscosity is applied.
725  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
726  real, dimension(SZIB_(G)), optional, intent(in) :: uhbt !< The summed volume flux
727  !! through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
728 
729  real, dimension(SZIB_(G)), intent(in) :: du_max_CFL !< Maximum acceptable
730  !! value of du [L T-1 ~> m s-1].
731  real, dimension(SZIB_(G)), intent(in) :: du_min_CFL !< Minimum acceptable
732  !! value of du [L T-1 ~> m s-1].
733  real, dimension(SZIB_(G)), intent(in) :: uh_tot_0 !< The summed transport
734  !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
735  real, dimension(SZIB_(G)), intent(in) :: duhdu_tot_0 !< The partial derivative
736  !! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
737  real, dimension(SZIB_(G)), intent(out) :: du !<
738  !! The barotropic velocity adjustment [L T-1 ~> m s-1].
739  real, intent(in) :: dt !< Time increment [T ~> s].
740  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
741  type(continuity_PPM_CS), pointer :: CS !< This module's control structure.
742  integer, intent(in) :: j !< Spatial index.
743  integer, intent(in) :: ish !< Start of index range.
744  integer, intent(in) :: ieh !< End of index range.
745  logical, dimension(SZIB_(G)), intent(in) :: do_I_in !<
746  !! A logical flag indicating which I values to work on.
747  logical, optional, intent(in) :: full_precision !<
748  !! A flag indicating how carefully to iterate. The
749  !! default is .true. (more accurate).
750  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(inout) :: uh_3d !<
751  !! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
752  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
753  ! Local variables
754  real, dimension(SZIB_(G),SZK_(G)) :: &
755  uh_aux, & ! An auxiliary zonal volume flux [H L2 s-1 ~> m3 s-1 or kg s-1].
756  duhdu ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
757  real, dimension(SZIB_(G)) :: &
758  uh_err, & ! Difference between uhbt and the summed uh [H L2 T-1 ~> m3 s-1 or kg s-1].
759  uh_err_best, & ! The smallest value of uh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
760  u_new, & ! The velocity with the correction added [L T-1 ~> m s-1].
761  duhdu_tot,&! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
762  du_min, & ! Min/max limits on du correction based on CFL limits
763  du_max ! and previous iterations [L T-1 ~> m s-1].
764  real :: du_prev ! The previous value of du [L T-1 ~> m s-1].
765  real :: ddu ! The change in du from the previous iteration [L T-1 ~> m s-1].
766  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
767  real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
768  integer :: i, k, nz, itt, max_itts = 20
769  logical :: full_prec, domore, do_I(SZIB_(G))
770 
771  nz = g%ke
772  full_prec = .true. ; if (present(full_precision)) full_prec = full_precision
773 
774  uh_aux(:,:) = 0.0 ; duhdu(:,:) = 0.0
775 
776  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
777  uh_aux(i,k) = uh_3d(i,j,k)
778  enddo ; enddo ; endif
779 
780  do i=ish-1,ieh
781  du(i) = 0.0 ; do_i(i) = do_i_in(i)
782  du_max(i) = du_max_cfl(i) ; du_min(i) = du_min_cfl(i)
783  uh_err(i) = uh_tot_0(i) - uhbt(i) ; duhdu_tot(i) = duhdu_tot_0(i)
784  uh_err_best(i) = abs(uh_err(i))
785  enddo
786 
787  do itt=1,max_itts
788  if (full_prec) then
789  select case (itt)
790  case (:1) ; tol_eta = 1e-6 * cs%tol_eta
791  case (2) ; tol_eta = 1e-4 * cs%tol_eta
792  case (3) ; tol_eta = 1e-2 * cs%tol_eta
793  case default ; tol_eta = cs%tol_eta
794  end select
795  else
796  tol_eta = cs%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
797  endif
798  tol_vel = cs%tol_vel
799 
800  do i=ish-1,ieh
801  if (uh_err(i) > 0.0) then ; du_max(i) = du(i)
802  elseif (uh_err(i) < 0.0) then ; du_min(i) = du(i)
803  else ; do_i(i) = .false. ; endif
804  enddo
805  domore = .false.
806  do i=ish-1,ieh ; if (do_i(i)) then
807  if ((dt * min(g%IareaT(i,j),g%IareaT(i+1,j))*abs(uh_err(i)) > tol_eta) .or. &
808  (cs%better_iter .and. ((abs(uh_err(i)) > tol_vel * duhdu_tot(i)) .or. &
809  (abs(uh_err(i)) > uh_err_best(i))) )) then
810  ! Use Newton's method, provided it stays bounded. Otherwise bisect
811  ! the value with the appropriate bound.
812  if (full_prec) then
813  ddu = -uh_err(i) / duhdu_tot(i)
814  du_prev = du(i)
815  du(i) = du(i) + ddu
816  if (abs(ddu) < 1.0e-15*abs(du(i))) then
817  do_i(i) = .false. ! ddu is small enough to quit.
818  elseif (ddu > 0.0) then
819  if (du(i) >= du_max(i)) then
820  du(i) = 0.5*(du_prev + du_max(i))
821  if (du_max(i) - du_prev < 1.0e-15*abs(du(i))) do_i(i) = .false.
822  endif
823  else ! ddu < 0.0
824  if (du(i) <= du_min(i)) then
825  du(i) = 0.5*(du_prev + du_min(i))
826  if (du_prev - du_min(i) < 1.0e-15*abs(du(i))) do_i(i) = .false.
827  endif
828  endif
829  else
830  ! Use Newton's method, provided it stays bounded, just like above.
831  du(i) = du(i) - uh_err(i) / duhdu_tot(i)
832  if ((du(i) >= du_max(i)) .or. (du(i) <= du_min(i))) &
833  du(i) = 0.5*(du_max(i) + du_min(i))
834  endif
835  if (do_i(i)) domore = .true.
836  else
837  do_i(i) = .false.
838  endif
839  endif ; enddo
840  if (.not.domore) exit
841 
842  if ((itt < max_itts) .or. present(uh_3d)) then ; do k=1,nz
843  do i=ish-1,ieh ; u_new(i) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
844  call zonal_flux_layer(u_new, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
845  uh_aux(:,k), duhdu(:,k), visc_rem(:,k), &
846  dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
847  enddo ; endif
848 
849  if (itt < max_itts) then
850  do i=ish-1,ieh
851  uh_err(i) = -uhbt(i) ; duhdu_tot(i) = 0.0
852  enddo
853  do k=1,nz ; do i=ish-1,ieh
854  uh_err(i) = uh_err(i) + uh_aux(i,k)
855  duhdu_tot(i) = duhdu_tot(i) + duhdu(i,k)
856  enddo ; enddo
857  do i=ish-1,ieh
858  uh_err_best(i) = min(uh_err_best(i), abs(uh_err(i)))
859  enddo
860  endif
861  enddo ! itt-loop
862  ! If there are any faces which have not converged to within the tolerance,
863  ! so-be-it, or else use a final upwind correction?
864  ! This never seems to happen with 20 iterations as max_itt.
865 
866  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
867  uh_3d(i,j,k) = uh_aux(i,k)
868  enddo ; enddo ; endif
869 

◆ zonal_flux_layer()

subroutine mom_continuity_ppm::zonal_flux_layer ( real, dimension( g %isdb: g %iedb), intent(in)  u,
real, dimension(szi_(g)), intent(in)  h,
real, dimension(szi_(g)), intent(in)  h_L,
real, dimension(szi_(g)), intent(in)  h_R,
real, dimension(szib_(g)), intent(inout)  uh,
real, dimension(szib_(g)), intent(inout)  duhdu,
real, dimension( g %isdb: g %iedb), intent(in)  visc_rem,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  j,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, dimension(szib_(g)), intent(in)  do_I,
logical, intent(in)  vol_CFL,
type(ocean_obc_type), optional, pointer  OBC 
)
private

Evaluates the zonal mass or volume fluxes in a layer.

Parameters
[in,out]gOcean's grid structure.
[in]uZonal velocity [L T-1 ~> m s-1].
[in]visc_remBoth the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]hLayer thickness [H ~> m or kg m-2].
[in]h_lLeft thickness [H ~> m or kg m-2].
[in]h_rRight thickness [H ~> m or kg m-2].
[in,out]uhZonal mass or volume transport [H L2 T-1 ~> m3 s-1 or kg s-1].
[in,out]duhduPartial derivative of uh with u [H L ~> m2 or kg m-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
[in]jSpatial index.
[in]ishStart of index range.
[in]iehEnd of index range.
[in]do_iWhich i values to work on.
[in]vol_cflIf true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
obcOpen boundaries control structure.

Definition at line 523 of file MOM_continuity_PPM.F90.

523  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
524  real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
525  real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the
526  !! momentum originally in a layer that remains after a time-step
527  !! of viscosity, and the fraction of a time-step's worth of a barotropic
528  !! acceleration that a layer experiences after viscosity is applied.
529  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
530  real, dimension(SZI_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
531  real, dimension(SZI_(G)), intent(in) :: h_L !< Left thickness [H ~> m or kg m-2].
532  real, dimension(SZI_(G)), intent(in) :: h_R !< Right thickness [H ~> m or kg m-2].
533  real, dimension(SZIB_(G)), intent(inout) :: uh !< Zonal mass or volume
534  !! transport [H L2 T-1 ~> m3 s-1 or kg s-1].
535  real, dimension(SZIB_(G)), intent(inout) :: duhdu !< Partial derivative of uh
536  !! with u [H L ~> m2 or kg m-1].
537  real, intent(in) :: dt !< Time increment [T ~> s].
538  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
539  integer, intent(in) :: j !< Spatial index.
540  integer, intent(in) :: ish !< Start of index range.
541  integer, intent(in) :: ieh !< End of index range.
542  logical, dimension(SZIB_(G)), intent(in) :: do_I !< Which i values to work on.
543  logical, intent(in) :: vol_CFL !< If true, rescale the
544  !! ratio of face areas to the cell areas when estimating the CFL number.
545  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
546  ! Local variables
547  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
548  real :: curv_3 ! A measure of the thickness curvature over a grid length,
549  ! with the same units as h_in.
550  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
551  integer :: i
552  integer :: l_seg
553  logical :: local_open_BC
554 
555  local_open_bc = .false.
556  if (present(obc)) then ; if (associated(obc)) then
557  local_open_bc = obc%open_u_BCs_exist_globally
558  endif ; endif
559 
560  do i=ish-1,ieh ; if (do_i(i)) then
561  ! Set new values of uh and duhdu.
562  if (u(i) > 0.0) then
563  if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
564  else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
565  curv_3 = h_l(i) + h_r(i) - 2.0*h(i)
566  uh(i) = g%dy_Cu(i,j) * u(i) * &
567  (h_r(i) + cfl * (0.5*(h_l(i) - h_r(i)) + curv_3*(cfl - 1.5)))
568  h_marg = h_r(i) + cfl * ((h_l(i) - h_r(i)) + 3.0*curv_3*(cfl - 1.0))
569  elseif (u(i) < 0.0) then
570  if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
571  else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
572  curv_3 = h_l(i+1) + h_r(i+1) - 2.0*h(i+1)
573  uh(i) = g%dy_Cu(i,j) * u(i) * &
574  (h_l(i+1) + cfl * (0.5*(h_r(i+1)-h_l(i+1)) + curv_3*(cfl - 1.5)))
575  h_marg = h_l(i+1) + cfl * ((h_r(i+1)-h_l(i+1)) + 3.0*curv_3*(cfl - 1.0))
576  else
577  uh(i) = 0.0
578  h_marg = 0.5 * (h_l(i+1) + h_r(i))
579  endif
580  duhdu(i) = g%dy_Cu(i,j) * h_marg * visc_rem(i)
581  endif ; enddo
582 
583  if (local_open_bc) then
584  do i=ish-1,ieh ; if (do_i(i)) then
585  l_seg = obc%segnum_u(i,j)
586 
587  if (l_seg /= obc_none) then
588  if (obc%segment(l_seg)%open) then
589  if (obc%segment(l_seg)%direction == obc_direction_e) then
590  uh(i) = g%dy_Cu(i,j) * u(i) * h(i)
591  duhdu(i) = g%dy_Cu(i,j) * h(i) * visc_rem(i)
592  else
593  uh(i) = g%dy_Cu(i,j) * u(i) * h(i+1)
594  duhdu(i) = g%dy_Cu(i,j) * h(i+1) * visc_rem(i)
595  endif
596  endif
597  endif
598  endif ; enddo
599  endif

◆ zonal_mass_flux()

subroutine mom_continuity_ppm::zonal_mass_flux ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_in,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  uh,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(continuity_ppm_cs), pointer  CS,
type(loop_bounds_type), intent(in)  LB,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional  uhbt,
type(ocean_obc_type), optional, pointer  OBC,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in), optional  visc_rem_u,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional  u_cor,
type(bt_cont_type), optional, pointer  BT_cont 
)
private

Calculates the mass or volume fluxes through the zonal faces, and other related quantities.

Parameters
[in,out]gOcean's grid structure.
[in]gvOcean's vertical grid structure.
[in]uZonal velocity [L T-1 ~> m s-1].
[in]h_inLayer thickness used to calculate fluxes [H ~> m or kg m-2].
[out]uhVolume flux through zonal faces = u*h*dy
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThis module's control structure.
[in]lbLoop bounds structure.
obcOpen boundaries control structure.
[in]visc_rem_uThe fraction of zonal momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]uhbtThe summed volume flux through zonal faces
[out]u_corThe zonal velocitiess (u with a barotropic correction) that give uhbt as the depth-integrated transport, m s-1.
bt_contA structure with elements that describe the effective open face areas as a function of barotropic flow.

Definition at line 213 of file MOM_continuity_PPM.F90.

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  integer :: l_seg
267  logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
268  logical :: local_Flather_OBC, local_open_BC, is_simple
269  type(OBC_segment_type), pointer :: segment => null()
270 
271  use_visc_rem = present(visc_rem_u)
272  local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
273  local_open_bc = .false.
274  if (present(bt_cont)) set_bt_cont = (associated(bt_cont))
275  if (present(obc)) then ; if (associated(obc)) then
276  local_specified_bc = obc%specified_u_BCs_exist_globally
277  local_flather_obc = obc%Flather_u_BCs_exist_globally
278  local_open_bc = obc%open_u_BCs_exist_globally
279  endif ; endif
280  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
281 
282  cfl_dt = cs%CFL_limit_adjust / dt
283  i_dt = 1.0 / dt
284  if (cs%aggress_adjust) cfl_dt = i_dt
285 
286  call cpu_clock_begin(id_clock_update)
287 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,CS,h_L,h_in,h_R,G,GV,LB,visc_rem,OBC)
288  do k=1,nz
289  ! This sets h_L and h_R.
290  if (cs%upwind_1st) then
291  do j=jsh,jeh ; do i=ish-1,ieh+1
292  h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
293  enddo ; enddo
294  else
295  call ppm_reconstruction_x(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
296  2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
297  endif
298  do i=ish-1,ieh ; visc_rem(i,k) = 1.0 ; enddo
299  enddo
300  call cpu_clock_end(id_clock_update)
301 
302  call cpu_clock_begin(id_clock_correct)
303 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,u,h_in,h_L,h_R,use_visc_rem,visc_rem_u, &
304 !$OMP uh,dt,US,G,GV,CS,local_specified_BC,OBC,uhbt,set_BT_cont, &
305 !$OMP CFL_dt,I_dt,u_cor,BT_cont, local_Flather_OBC) &
306 !$OMP private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0,duhdu_tot_0, &
307 !$OMP is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W, &
308 !$OMP any_simple_OBC,l_seg) &
309 !$OMP firstprivate(visc_rem)
310  do j=jsh,jeh
311  do i=ish-1,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ; enddo
312  ! Set uh and duhdu.
313  do k=1,nz
314  if (use_visc_rem) then ; do i=ish-1,ieh
315  visc_rem(i,k) = visc_rem_u(i,j,k)
316  visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
317  enddo ; endif
318  call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
319  uh(:,j,k), duhdu(:,k), visc_rem(:,k), &
320  dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
321  if (local_specified_bc) then
322  do i=ish-1,ieh
323  l_seg = obc%segnum_u(i,j)
324 
325  if (l_seg /= obc_none) then
326  if (obc%segment(l_seg)%specified) &
327  uh(i,j,k) = obc%segment(l_seg)%normal_trans(i,j,k)
328  endif
329  enddo
330  endif
331  enddo
332 
333  if ((.not.use_visc_rem).or.(.not.cs%use_visc_rem_max)) then ; do i=ish-1,ieh
334  visc_rem_max(i) = 1.0
335  enddo ; endif
336 
337  if (present(uhbt) .or. set_bt_cont) then
338  ! Set limits on du that will keep the CFL number between -1 and 1.
339  ! This should be adequate to keep the root bracketed in all cases.
340  do i=ish-1,ieh
341  i_vrm = 0.0
342  if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
343  if (cs%vol_CFL) then
344  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
345  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
346  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
347  du_max_cfl(i) = 2.0* (cfl_dt * dx_w) * i_vrm
348  du_min_cfl(i) = -2.0 * (cfl_dt * dx_e) * i_vrm
349  uh_tot_0(i) = 0.0 ; duhdu_tot_0(i) = 0.0
350  enddo
351  do k=1,nz ; do i=ish-1,ieh
352  duhdu_tot_0(i) = duhdu_tot_0(i) + duhdu(i,k)
353  uh_tot_0(i) = uh_tot_0(i) + uh(i,j,k)
354  enddo ; enddo
355  if (use_visc_rem) then
356  if (cs%aggress_adjust) then
357  do k=1,nz ; do i=ish-1,ieh
358  if (cs%vol_CFL) then
359  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
360  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
361  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
362 
363  du_lim = 0.499*((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k)))
364  if (du_max_cfl(i) * visc_rem(i,k) > du_lim) &
365  du_max_cfl(i) = du_lim / visc_rem(i,k)
366 
367  du_lim = 0.499*((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k)))
368  if (du_min_cfl(i) * visc_rem(i,k) < du_lim) &
369  du_min_cfl(i) = du_lim / visc_rem(i,k)
370  enddo ; enddo
371  else
372  do k=1,nz ; do i=ish-1,ieh
373  if (cs%vol_CFL) then
374  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
375  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
376  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
377 
378  if (du_max_cfl(i) * visc_rem(i,k) > dx_w*cfl_dt - u(i,j,k)) &
379  du_max_cfl(i) = (dx_w*cfl_dt - u(i,j,k)) / visc_rem(i,k)
380  if (du_min_cfl(i) * visc_rem(i,k) < -dx_e*cfl_dt - u(i,j,k)) &
381  du_min_cfl(i) = -(dx_e*cfl_dt + u(i,j,k)) / visc_rem(i,k)
382  enddo ; enddo
383  endif
384  else
385  if (cs%aggress_adjust) then
386  do k=1,nz ; do i=ish-1,ieh
387  if (cs%vol_CFL) then
388  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
389  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
390  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
391 
392  du_max_cfl(i) = min(du_max_cfl(i), 0.499 * &
393  ((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k))) )
394  du_min_cfl(i) = max(du_min_cfl(i), 0.499 * &
395  ((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k))) )
396  enddo ; enddo
397  else
398  do k=1,nz ; do i=ish-1,ieh
399  if (cs%vol_CFL) then
400  dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
401  dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
402  else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
403 
404  du_max_cfl(i) = min(du_max_cfl(i), dx_w*cfl_dt - u(i,j,k))
405  du_min_cfl(i) = max(du_min_cfl(i), -(dx_e*cfl_dt + u(i,j,k)))
406  enddo ; enddo
407  endif
408  endif
409  do i=ish-1,ieh
410  du_max_cfl(i) = max(du_max_cfl(i),0.0)
411  du_min_cfl(i) = min(du_min_cfl(i),0.0)
412  enddo
413 
414  any_simple_obc = .false.
415  if (present(uhbt) .or. set_bt_cont) then
416  if (local_specified_bc .or. local_flather_obc) then ; do i=ish-1,ieh
417  l_seg = obc%segnum_u(i,j)
418 
419  ! Avoid reconciling barotropic/baroclinic transports if transport is specified
420  is_simple = .false.
421  if (l_seg /= obc_none) &
422  is_simple = obc%segment(l_seg)%specified
423  do_i(i) = .not. (l_seg /= obc_none .and. is_simple)
424  any_simple_obc = any_simple_obc .or. is_simple
425  enddo ; else ; do i=ish-1,ieh
426  do_i(i) = .true.
427  enddo ; endif
428  endif
429 
430  if (present(uhbt)) then
431  call zonal_flux_adjust(u, h_in, h_l, h_r, uhbt(:,j), uh_tot_0, duhdu_tot_0, du, &
432  du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
433  j, ish, ieh, do_i, .true., uh, obc=obc)
434 
435  if (present(u_cor)) then ; do k=1,nz
436  do i=ish-1,ieh ; u_cor(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
437  if (local_specified_bc) then ; do i=ish-1,ieh
438  l_seg = obc%segnum_u(i,j)
439 
440  if (l_seg /= obc_none) then
441  if (obc%segment(l_seg)%specified) &
442  u_cor(i,j,k) = obc%segment(l_seg)%normal_vel(i,j,k)
443  endif
444  enddo ; endif
445  enddo ; endif ! u-corrected
446 
447  endif
448 
449  if (set_bt_cont) then
450  call set_zonal_bt_cont(u, h_in, h_l, h_r, bt_cont, uh_tot_0, duhdu_tot_0,&
451  du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
452  visc_rem_max, j, ish, ieh, do_i)
453  if (any_simple_obc) then
454  do i=ish-1,ieh
455  l_seg = obc%segnum_u(i,j)
456 
457  do_i(i) = .false.
458  if (l_seg /= obc_none) &
459  do_i(i) = obc%segment(l_seg)%specified
460 
461  if (do_i(i)) faui(i) = gv%H_subroundoff*g%dy_Cu(i,j)
462  enddo
463  ! NOTE: do_I(I) should prevent access to segment OBC_NONE
464  do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
465  if ((abs(obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
466  (obc%segment(obc%segnum_u(i,j))%specified)) &
467  faui(i) = faui(i) + obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k) / &
468  obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
469  endif ; enddo ; enddo
470  do i=ish-1,ieh ; if (do_i(i)) then
471  bt_cont%FA_u_W0(i,j) = faui(i) ; bt_cont%FA_u_E0(i,j) = faui(i)
472  bt_cont%FA_u_WW(i,j) = faui(i) ; bt_cont%FA_u_EE(i,j) = faui(i)
473  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
474  endif ; enddo
475  endif
476  endif ! set_BT_cont
477 
478  endif ! present(uhbt) or set_BT_cont
479 
480  enddo ! j-loop
481 
482  if (local_open_bc .and. set_bt_cont) then
483  do n = 1, obc%number_of_segments
484  if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
485  i = obc%segment(n)%HI%IsdB
486  if (obc%segment(n)%direction == obc_direction_e) then
487  do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
488  fa_u = 0.0
489  do k=1,nz ; fa_u = fa_u + h_in(i,j,k)*g%dy_Cu(i,j) ; enddo
490  bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
491  bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
492  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
493  enddo
494  else
495  do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
496  fa_u = 0.0
497  do k=1,nz ; fa_u = fa_u + h_in(i+1,j,k)*g%dy_Cu(i,j) ; enddo
498  bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
499  bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
500  bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
501  enddo
502  endif
503  endif
504  enddo
505  endif
506  call cpu_clock_end(id_clock_correct)
507 
508  if (set_bt_cont) then ; if (allocated(bt_cont%h_u)) then
509  if (present(u_cor)) then
510  call zonal_face_thickness(u_cor, h_in, h_l, h_r, bt_cont%h_u, dt, g, us, lb, &
511  cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
512  else
513  call zonal_face_thickness(u, h_in, h_l, h_r, bt_cont%h_u, dt, g, us, lb, &
514  cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
515  endif
516  endif ; endif
517