MOM6
Data Types
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...
 
integer id_clock_update
 CPU time clock IDs. More...
 
integer id_clock_correct
 CPU time clock IDs. More...
 
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...
 

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 

References id_clock_update, meridional_mass_flux(), mom_error_handler::mom_error(), and zonal_mass_flux().

Here is the call graph for this function:

◆ 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 2286 of file MOM_continuity_PPM.F90.

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

Referenced by mom_continuity::continuity_end().

Here is the caller graph for this function:

◆ 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 2181 of file MOM_continuity_PPM.F90.

2181  type(time_type), target, intent(in) :: Time !< The current model time.
2182  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2183  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
2184  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2185  type(param_file_type), intent(in) :: param_file !< A structure indicating
2186  !! the open file to parse for model parameter values.
2187  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
2188  !! regulate diagnostic output.
2189  type(continuity_PPM_CS), pointer :: CS !< Module's control structure.
2190 !> This include declares and sets the variable "version".
2191 #include "version_variable.h"
2192  real :: tol_eta_m ! An unscaled version of tol_eta [m].
2193  character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name.
2194 
2195  if (associated(cs)) then
2196  call mom_error(warning, "continuity_PPM_init called with associated control structure.")
2197  return
2198  endif
2199  allocate(cs)
2200 
2201 ! Read all relevant parameters and write them to the model log.
2202  call log_version(param_file, mdl, version, "")
2203  call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", cs%monotonic, &
2204  "If true, CONTINUITY_PPM uses the Colella and Woodward "//&
2205  "monotonic limiter. The default (false) is to use a "//&
2206  "simple positive definite limiter.", default=.false.)
2207  call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &
2208  "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2209  "(arithmetic mean) interpolation of the edge values. "//&
2210  "This may give better PV conservation properties. While "//&
2211  "it formally reduces the accuracy of the continuity "//&
2212  "solver itself in the strongly advective limit, it does "//&
2213  "not reduce the overall order of accuracy of the dynamic "//&
2214  "core.", default=.false.)
2215  call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
2216  "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
2217  "continuity solver. This scheme is highly diffusive "//&
2218  "but may be useful for debugging or in single-column "//&
2219  "mode where its minimal stencil is useful.", default=.false.)
2220  call get_param(param_file, mdl, "ETA_TOLERANCE", cs%tol_eta, &
2221  "The tolerance for the differences between the "//&
2222  "barotropic and baroclinic estimates of the sea surface "//&
2223  "height due to the fluxes through each face. The total "//&
2224  "tolerance for SSH is 4 times this value. The default "//&
2225  "is 0.5*NK*ANGSTROM, and this should not be set less "//&
2226  "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=gv%m_to_H, &
2227  default=0.5*g%ke*gv%Angstrom_m, unscaled=tol_eta_m)
2228 
2229  !### ETA_TOLERANCE_AUX can be obsoleted.
2230  call get_param(param_file, mdl, "ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
2231  "The tolerance for free-surface height discrepancies "//&
2232  "between the barotropic solution and the sum of the "//&
2233  "layer thicknesses when calculating the auxiliary "//&
2234  "corrected velocities. By default, this is the same as "//&
2235  "ETA_TOLERANCE, but can be made larger for efficiency.", &
2236  units="m", default=tol_eta_m, scale=gv%m_to_H)
2237  call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &
2238  "The tolerance for barotropic velocity discrepancies "//&
2239  "between the barotropic solution and the sum of the "//&
2240  "layer thicknesses.", units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)
2241  ! The speed of light is the default.
2242 
2243  call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&
2244  "If true, allow the adjusted velocities to have a "//&
2245  "relative CFL change up to 0.5.", default=.false.)
2246  cs%vol_CFL = cs%aggress_adjust
2247  call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &
2248  "If true, use the ratio of the open face lengths to the "//&
2249  "tracer cell areas when estimating CFL numbers. The "//&
2250  "default is set by CONT_PPM_AGGRESS_ADJUST.", &
2251  default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)
2252  call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &
2253  "The maximum CFL of the adjusted velocities.", units="nondim", &
2254  default=0.5)
2255  call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", cs%better_iter, &
2256  "If true, stop corrective iterations using a velocity "//&
2257  "based criterion and only stop if the iteration is "//&
2258  "better than all predecessors.", default=.true.)
2259  call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", &
2260  cs%use_visc_rem_max, &
2261  "If true, use more appropriate limiting bounds for "//&
2262  "corrections in strongly viscous columns.", default=.true.)
2263  call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &
2264  "If true, use the marginal face areas from the continuity "//&
2265  "solver for use as the weights in the barotropic solver. "//&
2266  "Otherwise use the transport averaged areas.", default=.true.)
2267 
2268  cs%diag => diag
2269 
2270  id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=clock_routine)
2271  id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=clock_routine)
2272 

References id_clock_correct, id_clock_update, and mom_error_handler::mom_error().

Here is the call graph for this function:

◆ 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 2277 of file MOM_continuity_PPM.F90.

2277  type(continuity_PPM_CS), pointer :: CS !< Module's control structure.
2278  integer :: stencil !< The continuity solver stencil size with the current settings.
2279 
2280  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1
2281 

◆ 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 1379 of file MOM_continuity_PPM.F90.

1379  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1380  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1381  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,
1382  !! [H ~> m or kg m-2].
1383  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_L !< Left thickness in the reconstruction,
1384  !! [H ~> m or kg m-2].
1385  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_R !< Right thickness in the reconstruction,
1386  !! [H ~> m or kg m-2].
1387  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h_v !< Thickness at meridional faces,
1388  !! [H ~> m or kg m-2].
1389  real, intent(in) :: dt !< Time increment [T ~> s].
1390  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1391  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1392  logical, intent(in) :: vol_CFL !< If true, rescale the ratio
1393  !! of face areas to the cell areas when estimating the CFL number.
1394  logical, intent(in) :: marginal !< If true, report the marginal
1395  !! face thicknesses; otherwise report transport-averaged thicknesses.
1396  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(in) :: visc_rem_v !< Both the fraction
1397  !! of the momentum originally in a layer that remains after a time-step of
1398  !! viscosity, and the fraction of a time-step's worth of a barotropic
1399  !! acceleration that a layer experiences after viscosity is applied.
1400  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1401  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1402 
1403  ! Local variables
1404  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1405  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1406  ! with the same units as h [H ~> m or kg m-2] .
1407  real :: h_avg ! The average thickness of a flux [H ~> m or kg m-2].
1408  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1409  logical :: local_open_BC
1410  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1411  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1412 
1413  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
1414  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
1415  if (v(i,j,k) > 0.0) then
1416  if (vol_cfl) then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1417  else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ; endif
1418  curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
1419  h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
1420  h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + &
1421  3.0*curv_3*(cfl - 1.0))
1422  elseif (v(i,j,k) < 0.0) then
1423  if (vol_cfl) then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1424  else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ; endif
1425  curv_3 = h_l(i,j+1,k) + h_r(i,j+1,k) - 2.0*h(i,j+1,k)
1426  h_avg = h_l(i,j+1,k) + cfl * (0.5*(h_r(i,j+1,k)-h_l(i,j+1,k)) + curv_3*(cfl - 1.5))
1427  h_marg = h_l(i,j+1,k) + cfl * ((h_r(i,j+1,k)-h_l(i,j+1,k)) + &
1428  3.0*curv_3*(cfl - 1.0))
1429  else
1430  h_avg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1431  ! The choice to use the arithmetic mean here is somewhat arbitrariy, but
1432  ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
1433  h_marg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1434  ! h_marg = (2.0 * h_L(i,j+1,k) * h_R(i,j,k)) / &
1435  ! (h_L(i,j+1,k) + h_R(i,j,k) + GV%H_subroundoff)
1436  endif
1437 
1438  if (marginal) then ; h_v(i,j,k) = h_marg
1439  else ; h_v(i,j,k) = h_avg ; endif
1440  enddo ; enddo ; enddo
1441 
1442  if (present(visc_rem_v)) then
1443  !$OMP parallel do default(shared)
1444  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
1445  h_v(i,j,k) = h_v(i,j,k) * visc_rem_v(i,j,k)
1446  enddo ; enddo ; enddo
1447  endif
1448 
1449  local_open_bc = .false.
1450  if (present(obc)) then ; if (associated(obc)) then
1451  local_open_bc = obc%open_u_BCs_exist_globally
1452  endif ; endif
1453  if (local_open_bc) then
1454  do n = 1, obc%number_of_segments
1455  if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
1456  j = obc%segment(n)%HI%JsdB
1457  if (obc%segment(n)%direction == obc_direction_n) then
1458  if (present(visc_rem_v)) then ; do k=1,nz
1459  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1460  h_v(i,j,k) = h(i,j,k) * visc_rem_v(i,j,k)
1461  enddo
1462  enddo ; else ; do k=1,nz
1463  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1464  h_v(i,j,k) = h(i,j,k)
1465  enddo
1466  enddo ; endif
1467  else
1468  if (present(visc_rem_v)) then ; do k=1,nz
1469  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1470  h_v(i,j,k) = h(i,j+1,k) * visc_rem_v(i,j,k)
1471  enddo
1472  enddo ; else ; do k=1,nz
1473  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1474  h_v(i,j,k) = h(i,j+1,k)
1475  enddo
1476  enddo ; endif
1477  endif
1478  endif
1479  enddo
1480  endif
1481 

References mom_open_boundary::obc_direction_n.

Referenced by meridional_mass_flux().

Here is the caller graph for this function:

◆ 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 1298 of file MOM_continuity_PPM.F90.

1298  type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure.
1299  real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
1300  real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the
1301  !! momentum originally in a layer that remains after a time-step
1302  !! of viscosity, and the fraction of a time-step's worth of a barotropic
1303  !! acceleration that a layer experiences after viscosity is applied.
1304  !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
1305  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Layer thickness used to calculate fluxes,
1306  !! [H ~> m or kg m-2].
1307  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_L !< Left thickness in the reconstruction
1308  !! [H ~> m or kg m-2].
1309  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_R !< Right thickness in the reconstruction
1310  !! [H ~> m or kg m-2].
1311  real, dimension(SZI_(G)), intent(inout) :: vh !< Meridional mass or volume transport
1312  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
1313  real, dimension(SZI_(G)), intent(inout) :: dvhdv !< Partial derivative of vh with v
1314  !! [H L ~> m2 or kg m-1].
1315  real, intent(in) :: dt !< Time increment [T ~> s].
1316  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1317  integer, intent(in) :: j !< Spatial index.
1318  integer, intent(in) :: ish !< Start of index range.
1319  integer, intent(in) :: ieh !< End of index range.
1320  logical, dimension(SZI_(G)), intent(in) :: do_I !< Which i values to work on.
1321  logical, intent(in) :: vol_CFL !< If true, rescale the
1322  !! ratio of face areas to the cell areas when estimating the CFL number.
1323  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1324  ! Local variables
1325  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
1326  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1327  ! with the same units as h, i.e. [H ~> m or kg m-2].
1328  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
1329  integer :: i
1330  logical :: local_open_BC
1331 
1332  local_open_bc = .false.
1333  if (present(obc)) then ; if (associated(obc)) then
1334  local_open_bc = obc%open_u_BCs_exist_globally
1335  endif ; endif
1336 
1337  do i=ish,ieh ; if (do_i(i)) then
1338  if (v(i) > 0.0) then
1339  if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1340  else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
1341  curv_3 = h_l(i,j) + h_r(i,j) - 2.0*h(i,j)
1342  vh(i) = g%dx_Cv(i,j) * v(i) * ( h_r(i,j) + cfl * &
1343  (0.5*(h_l(i,j) - h_r(i,j)) + curv_3*(cfl - 1.5)) )
1344  h_marg = h_r(i,j) + cfl * ((h_l(i,j) - h_r(i,j)) + &
1345  3.0*curv_3*(cfl - 1.0))
1346  elseif (v(i) < 0.0) then
1347  if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1348  else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
1349  curv_3 = h_l(i,j+1) + h_r(i,j+1) - 2.0*h(i,j+1)
1350  vh(i) = g%dx_Cv(i,j) * v(i) * ( h_l(i,j+1) + cfl * &
1351  (0.5*(h_r(i,j+1)-h_l(i,j+1)) + curv_3*(cfl - 1.5)) )
1352  h_marg = h_l(i,j+1) + cfl * ((h_r(i,j+1)-h_l(i,j+1)) + &
1353  3.0*curv_3*(cfl - 1.0))
1354  else
1355  vh(i) = 0.0
1356  h_marg = 0.5 * (h_l(i,j+1) + h_r(i,j))
1357  endif
1358  dvhdv(i) = g%dx_Cv(i,j) * h_marg * visc_rem(i)
1359  endif ; enddo
1360 
1361  if (local_open_bc) then
1362  do i=ish,ieh ; if (do_i(i)) then
1363  if (obc%segment(obc%segnum_v(i,j))%open) then
1364  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
1365  vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j)
1366  dvhdv(i) = g%dx_Cv(i,j) * h(i,j) * visc_rem(i)
1367  else
1368  vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j+1)
1369  dvhdv(i) = g%dx_Cv(i,j) * h(i,j+1) * visc_rem(i)
1370  endif
1371  endif
1372  endif ; enddo
1373  endif

References mom_open_boundary::obc_direction_n.

Referenced by meridional_flux_adjust(), meridional_mass_flux(), and set_merid_bt_cont().

Here is the caller graph for this function:

◆ 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 1488 of file MOM_continuity_PPM.F90.

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

References merid_flux_layer().

Referenced by meridional_mass_flux(), and set_merid_bt_cont().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 1014 of file MOM_continuity_PPM.F90.

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

References id_clock_correct, id_clock_update, merid_face_thickness(), merid_flux_layer(), meridional_flux_adjust(), mom_open_boundary::obc_direction_n, ppm_reconstruction_y(), ratio_max(), and set_merid_bt_cont().

Referenced by continuity_ppm().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 2130 of file MOM_continuity_PPM.F90.

2130  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2131  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2132  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction,
2133  !! [H ~> m or kg m-2].
2134  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction,
2135  !! [H ~> m or kg m-2].
2136  integer, intent(in) :: iis !< Start of i index range.
2137  integer, intent(in) :: iie !< End of i index range.
2138  integer, intent(in) :: jis !< Start of j index range.
2139  integer, intent(in) :: jie !< End of j index range.
2140 
2141  ! Local variables
2142  real :: h_i, RLdiff, RLdiff2, RLmean, FunFac
2143  character(len=256) :: mesg
2144  integer :: i,j
2145 
2146  do j=jis,jie ; do i=iis,iie
2147  ! This limiter monotonizes the parabola following
2148  ! Colella and Woodward, 1984, Eq. 1.10
2149  h_i = h_in(i,j)
2150  if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. ) then
2151  h_l(i,j) = h_i ; h_r(i,j) = h_i
2152  else
2153  rldiff = h_r(i,j) - h_l(i,j) ! Difference of edge values
2154  rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) ) ! Mean of edge values
2155  funfac = 6. * rldiff * ( h_i - rlmean ) ! Some funny factor
2156  rldiff2 = rldiff * rldiff ! Square of difference
2157  if ( funfac > rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)
2158  if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)
2159  endif
2160  enddo ; enddo
2161 
2162  return

Referenced by ppm_reconstruction_x(), and ppm_reconstruction_y().

Here is the caller graph for this function:

◆ 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 2089 of file MOM_continuity_PPM.F90.

2089  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
2090  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
2091  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
2092  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
2093  real, intent(in) :: h_min !< The minimum thickness
2094  !! that can be obtained by a concave parabolic fit.
2095  integer, intent(in) :: iis !< Start of i index range.
2096  integer, intent(in) :: iie !< End of i index range.
2097  integer, intent(in) :: jis !< Start of j index range.
2098  integer, intent(in) :: jie !< End of j index range.
2099 
2100 ! Local variables
2101  real :: curv, dh, scale
2102  character(len=256) :: mesg
2103  integer :: i,j
2104 
2105  do j=jis,jie ; do i=iis,iie
2106  ! This limiter prevents undershooting minima within the domain with
2107  ! values less than h_min.
2108  curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2109  if (curv > 0.0) then ! Only minima are limited.
2110  dh = h_r(i,j) - h_l(i,j)
2111  if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
2112  if (h_in(i,j) <= h_min) then
2113  h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2114  elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
2115  ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
2116  ! be limited in this case. 0 < scale < 1.
2117  scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2118  h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2119  h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2120  endif
2121  endif
2122  endif
2123  enddo ; enddo
2124 

Referenced by ppm_reconstruction_x(), and ppm_reconstruction_y().

Here is the caller graph for this function:

◆ 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 1810 of file MOM_continuity_PPM.F90.

1810  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1811  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
1812  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_L !< Left thickness in the reconstruction,
1813  !! [H ~> m or kg m-2].
1814  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_R !< Right thickness in the reconstruction,
1815  !! [H ~> m or kg m-2].
1816  type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
1817  real, intent(in) :: h_min !< The minimum thickness
1818  !! that can be obtained by a concave parabolic fit.
1819  logical, optional, intent(in) :: monotonic !< If true, use the
1820  !! Colella & Woodward monotonic limiter.
1821  !! Otherwise use a simple positive-definite limiter.
1822  logical, optional, intent(in) :: simple_2nd !< If true, use the
1823  !! arithmetic mean thicknesses as the default edge values
1824  !! for a simple 2nd order scheme.
1825  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1826 
1827  ! Local variables with useful mnemonic names.
1828  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1829  real, parameter :: oneSixth = 1./6.
1830  real :: h_ip1, h_im1
1831  real :: dMx, dMn
1832  logical :: use_CW84, use_2nd
1833  character(len=256) :: mesg
1834  integer :: i, j, isl, iel, jsl, jel, n, stencil
1835  logical :: local_open_BC
1836  type(OBC_segment_type), pointer :: segment => null()
1837 
1838  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
1839  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1840 
1841  local_open_bc = .false.
1842  if (present(obc)) then ; if (associated(obc)) then
1843  local_open_bc = obc%open_u_BCs_exist_globally
1844  endif ; endif
1845 
1846  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1847 
1848  ! This is the stencil of the reconstruction, not the scheme overall.
1849  stencil = 2 ; if (use_2nd) stencil = 1
1850 
1851  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
1852  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1853  & "x-halo that needs to be increased by ",i2,".")') &
1854  stencil + max(g%isd-isl,iel-g%ied)
1855  call mom_error(fatal,mesg)
1856  endif
1857  if ((jsl < g%jsd) .or. (jel > g%jed)) then
1858  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1859  & "y-halo that needs to be increased by ",i2,".")') &
1860  max(g%jsd-jsl,jel-g%jed)
1861  call mom_error(fatal,mesg)
1862  endif
1863 
1864  if (use_2nd) then
1865  do j=jsl,jel ; do i=isl,iel
1866  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1867  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1868  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1869  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1870  enddo ; enddo
1871  else
1872  do j=jsl,jel ; do i=isl-1,iel+1
1873  if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
1874  slp(i,j) = 0.0
1875  else
1876  ! This uses a simple 2nd order slope.
1877  slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1878  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1879  dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1880  dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1881  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1882  ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
1883  endif
1884  enddo ; enddo
1885 
1886  if (local_open_bc) then
1887  do n=1, obc%number_of_segments
1888  segment => obc%segment(n)
1889  if (.not. segment%on_pe) cycle
1890  if (segment%direction == obc_direction_e .or. &
1891  segment%direction == obc_direction_w) then
1892  i=segment%HI%IsdB
1893  do j=segment%HI%jsd,segment%HI%jed
1894  slp(i+1,j) = 0.0
1895  slp(i,j) = 0.0
1896  enddo
1897  endif
1898  enddo
1899  endif
1900 
1901  do j=jsl,jel ; do i=isl,iel
1902  ! Neighboring values should take into account any boundaries. The 3
1903  ! following sets of expressions are equivalent.
1904  ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
1905  ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
1906  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1907  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1908  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
1909  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1910  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1911  enddo ; enddo
1912  endif
1913 
1914  if (local_open_bc) then
1915  do n=1, obc%number_of_segments
1916  segment => obc%segment(n)
1917  if (.not. segment%on_pe) cycle
1918  if (segment%direction == obc_direction_e) then
1919  i=segment%HI%IsdB
1920  do j=segment%HI%jsd,segment%HI%jed
1921  h_l(i+1,j) = h_in(i,j)
1922  h_r(i+1,j) = h_in(i,j)
1923  h_l(i,j) = h_in(i,j)
1924  h_r(i,j) = h_in(i,j)
1925  enddo
1926  elseif (segment%direction == obc_direction_w) then
1927  i=segment%HI%IsdB
1928  do j=segment%HI%jsd,segment%HI%jed
1929  h_l(i,j) = h_in(i+1,j)
1930  h_r(i,j) = h_in(i+1,j)
1931  h_l(i+1,j) = h_in(i+1,j)
1932  h_r(i+1,j) = h_in(i+1,j)
1933  enddo
1934  endif
1935  enddo
1936  endif
1937 
1938  if (use_cw84) then
1939  call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
1940  else
1941  call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
1942  endif
1943 
1944  return

References mom_error_handler::mom_error(), mom_open_boundary::obc_direction_e, mom_open_boundary::obc_direction_w, ppm_limit_cw84(), and ppm_limit_pos().

Referenced by zonal_mass_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 1949 of file MOM_continuity_PPM.F90.

1949  type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
1950  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2].
1951  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_L !< Left thickness in the reconstruction,
1952  !! [H ~> m or kg m-2].
1953  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_R !< Right thickness in the reconstruction,
1954  !! [H ~> m or kg m-2].
1955  type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure.
1956  real, intent(in) :: h_min !< The minimum thickness
1957  !! that can be obtained by a concave parabolic fit.
1958  logical, optional, intent(in) :: monotonic !< If true, use the
1959  !! Colella & Woodward monotonic limiter.
1960  !! Otherwise use a simple positive-definite limiter.
1961  logical, optional, intent(in) :: simple_2nd !< If true, use the
1962  !! arithmetic mean thicknesses as the default edge values
1963  !! for a simple 2nd order scheme.
1964  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
1965 
1966  ! Local variables with useful mnemonic names.
1967  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1968  real, parameter :: oneSixth = 1./6.
1969  real :: h_jp1, h_jm1
1970  real :: dMx, dMn
1971  logical :: use_CW84, use_2nd
1972  character(len=256) :: mesg
1973  integer :: i, j, isl, iel, jsl, jel, n, stencil
1974  logical :: local_open_BC
1975  type(OBC_segment_type), pointer :: segment => null()
1976 
1977  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
1978  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1979 
1980  local_open_bc = .false.
1981  if (present(obc)) then ; if (associated(obc)) then
1982  local_open_bc = obc%open_u_BCs_exist_globally
1983  endif ; endif
1984 
1985  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
1986 
1987  ! This is the stencil of the reconstruction, not the scheme overall.
1988  stencil = 2 ; if (use_2nd) stencil = 1
1989 
1990  if ((isl < g%isd) .or. (iel > g%ied)) then
1991  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
1992  & "x-halo that needs to be increased by ",i2,".")') &
1993  max(g%isd-isl,iel-g%ied)
1994  call mom_error(fatal,mesg)
1995  endif
1996  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
1997  write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
1998  & "y-halo that needs to be increased by ",i2,".")') &
1999  stencil + max(g%jsd-jsl,jel-g%jed)
2000  call mom_error(fatal,mesg)
2001  endif
2002 
2003  if (use_2nd) then
2004  do j=jsl,jel ; do i=isl,iel
2005  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2006  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2007  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
2008  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
2009  enddo ; enddo
2010  else
2011  do j=jsl-1,jel+1 ; do i=isl,iel
2012  if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
2013  slp(i,j) = 0.0
2014  else
2015  ! This uses a simple 2nd order slope.
2016  slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
2017  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
2018  dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
2019  dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
2020  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
2021  ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
2022  endif
2023  enddo ; enddo
2024 
2025  if (local_open_bc) then
2026  do n=1, obc%number_of_segments
2027  segment => obc%segment(n)
2028  if (.not. segment%on_pe) cycle
2029  if (segment%direction == obc_direction_s .or. &
2030  segment%direction == obc_direction_n) then
2031  j=segment%HI%JsdB
2032  do i=segment%HI%isd,segment%HI%ied
2033  slp(i,j+1) = 0.0
2034  slp(i,j) = 0.0
2035  enddo
2036  endif
2037  enddo
2038  endif
2039 
2040  do j=jsl,jel ; do i=isl,iel
2041  ! Neighboring values should take into account any boundaries. The 3
2042  ! following sets of expressions are equivalent.
2043  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2044  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2045  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2046  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2047  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2048  enddo ; enddo
2049  endif
2050 
2051  if (local_open_bc) then
2052  do n=1, obc%number_of_segments
2053  segment => obc%segment(n)
2054  if (.not. segment%on_pe) cycle
2055  if (segment%direction == obc_direction_n) then
2056  j=segment%HI%JsdB
2057  do i=segment%HI%isd,segment%HI%ied
2058  h_l(i,j+1) = h_in(i,j)
2059  h_r(i,j+1) = h_in(i,j)
2060  h_l(i,j) = h_in(i,j)
2061  h_r(i,j) = h_in(i,j)
2062  enddo
2063  elseif (segment%direction == obc_direction_s) then
2064  j=segment%HI%JsdB
2065  do i=segment%HI%isd,segment%HI%ied
2066  h_l(i,j) = h_in(i,j+1)
2067  h_r(i,j) = h_in(i,j+1)
2068  h_l(i,j+1) = h_in(i,j+1)
2069  h_r(i,j+1) = h_in(i,j+1)
2070  enddo
2071  endif
2072  enddo
2073  endif
2074 
2075  if (use_cw84) then
2076  call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
2077  else
2078  call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
2079  endif
2080 
2081  return

References mom_error_handler::mom_error(), mom_open_boundary::obc_direction_n, mom_open_boundary::obc_direction_s, ppm_limit_cw84(), and ppm_limit_pos().

Referenced by meridional_mass_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 2167 of file MOM_continuity_PPM.F90.

2167  real, intent(in) :: a !< Numerator
2168  real, intent(in) :: b !< Denominator
2169  real, intent(in) :: maxrat !< Maximum value of ratio.
2170  real :: ratio !< Return value.
2171 
2172  if (abs(a) > abs(maxrat*b)) then
2173  ratio = maxrat
2174  else
2175  ratio = a / b
2176  endif

Referenced by meridional_mass_flux(), and zonal_mass_flux().

Here is the caller graph for this function:

◆ 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 1651 of file MOM_continuity_PPM.F90.

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

References merid_flux_layer(), and meridional_flux_adjust().

Referenced by meridional_mass_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 852 of file MOM_continuity_PPM.F90.

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

References zonal_flux_adjust(), and zonal_flux_layer().

Referenced by zonal_mass_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 580 of file MOM_continuity_PPM.F90.

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

References mom_open_boundary::obc_direction_e.

Referenced by zonal_mass_flux().

Here is the caller graph for this function:

◆ 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 688 of file MOM_continuity_PPM.F90.

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

References zonal_flux_layer().

Referenced by set_zonal_bt_cont(), and zonal_mass_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 503 of file MOM_continuity_PPM.F90.

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

References mom_open_boundary::obc_direction_e.

Referenced by set_zonal_bt_cont(), zonal_flux_adjust(), and zonal_mass_flux().

Here is the caller graph for this function:

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

References id_clock_correct, id_clock_update, mom_open_boundary::obc_direction_e, ppm_reconstruction_x(), ratio_max(), set_zonal_bt_cont(), zonal_face_thickness(), zonal_flux_adjust(), and zonal_flux_layer().

Referenced by continuity_ppm().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ id_clock_correct

integer mom_continuity_ppm::id_clock_correct
private

CPU time clock IDs.

Definition at line 24 of file MOM_continuity_PPM.F90.

Referenced by continuity_ppm_init(), meridional_mass_flux(), and zonal_mass_flux().

◆ id_clock_update

integer mom_continuity_ppm::id_clock_update
private

CPU time clock IDs.

Definition at line 24 of file MOM_continuity_PPM.F90.

24 integer :: id_clock_update, id_clock_correct

Referenced by continuity_ppm(), continuity_ppm_init(), meridional_mass_flux(), and zonal_mass_flux().