src/core/MOM_CoriolisAdv.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		#include "do_concurrent_compat.h"
6
7		!> Accelerations due to the Coriolis force and momentum advection
8		module MOM_CoriolisAdv
9
10		!> \author Robert Hallberg, April 1994 - June 2002
11
12		use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
13		use MOM_diag_mediator, only : post_product_u, post_product_sum_u
14		use MOM_diag_mediator, only : post_product_v, post_product_sum_v
15		use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
16		use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING
17		use MOM_file_parser, only : get_param, log_version, param_file_type
18		use MOM_grid, only : ocean_grid_type
19		use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
20		use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S
21		use MOM_open_boundary, only : OBC_VORTICITY_ZERO, OBC_VORTICITY_FREESLIP
22		use MOM_open_boundary, only : OBC_VORTICITY_COMPUTED, OBC_VORTICITY_SPECIFIED
23		use MOM_string_functions, only : uppercase
24		use MOM_unit_scaling, only : unit_scale_type
25		use MOM_variables, only : accel_diag_ptrs, porous_barrier_type
26		use MOM_verticalGrid, only : verticalGrid_type
27		use MOM_wave_interface, only : wave_parameters_CS
28
29		implicit none ; private
30
31		public CorAdCalc, CoriolisAdv_init, CoriolisAdv_end, CoriolisAdv_stencil
32
33		#include <MOM_memory.h>
34
35		!> Control structure for mom_coriolisadv
36		type, public :: CoriolisAdv_CS ; private
37		logical :: initialized = .false. !< True if this control structure has been initialized.
38		integer :: Coriolis_Scheme !< Selects the discretization for the Coriolis terms.
39		!! Valid values are:
40		!! - SADOURNY75_ENERGY - Sadourny, 1975
41		!! - ARAKAWA_HSU90 - Arakawa & Hsu, 1990, Energy & non-div. Enstrophy
42		!! - ROBUST_ENSTRO - Pseudo-enstrophy scheme
43		!! - SADOURNY75_ENSTRO - Sadourny, JAS 1975, Enstrophy
44		!! - ARAKAWA_LAMB81 - Arakawa & Lamb, MWR 1981, Energy & Enstrophy
45		!! - ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with Arakawa & Hsu and Sadourny energy.
46		!! - WENOVI3RD_PV_ENSTRO - 3rd-order WENO scheme for PV reconstruction
47		!! - WENOVI5TH_PV_ENSTRO - 5th-order WENO scheme for PV reconstruction
48		!! - WENOVI7TH_PV_ENSTRO - 7th-order WENO scheme for PV reconstruction
49		!! The default, SADOURNY75_ENERGY, is the safest choice then the
50		!! deformation radius is poorly resolved.
51		integer :: KE_Scheme !< KE_SCHEME selects the discretization for
52		!! the kinetic energy. Valid values are:
53		!! KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV
54		logical :: KE_use_limiter !< If true, use the Koren limiter for KE_UP3 scheme
55		integer :: PV_Adv_Scheme !< PV_ADV_SCHEME selects the discretization for PV advection
56		!! Valid values are:
57		!! - PV_ADV_CENTERED - centered (aka Sadourny, 75)
58		!! - PV_ADV_UPWIND1 - upwind, first order
59		integer :: nkblock !< The k block size used in Coriolis/advection calculations [nondim].
60		real :: F_eff_max_blend !< The factor by which the maximum effective Coriolis
61		!! acceleration from any point can be increased when
62		!! blending different discretizations with the
63		!! ARAKAWA_LAMB_BLEND Coriolis scheme [nondim].
64		!! This must be greater than 2.0, and is 4.0 by default.
65		real :: wt_lin_blend !< A weighting value beyond which the blending between
66		!! Sadourny and Arakawa & Hsu goes linearly to 0 [nondim].
67		!! This must be between 1 and 1e-15, often 1/8.
68		logical :: no_slip !< If true, no slip boundary conditions are used.
69		!! Otherwise free slip boundary conditions are assumed.
70		!! The implementation of the free slip boundary
71		!! conditions on a C-grid is much cleaner than the
72		!! no slip boundary conditions. The use of free slip
73		!! b.c.s is strongly encouraged. The no slip b.c.s
74		!! are not implemented with the biharmonic viscosity.
75		logical :: bound_Coriolis !< If true, the Coriolis terms at u points are
76		!! bounded by the four estimates of (f+rv)v from the
77		!! four neighboring v points, and similarly at v
78		!! points. This option would have no effect on the
79		!! SADOURNY75_ENERGY scheme if it were possible to
80		!! use centered difference thickness fluxes.
81		logical :: Coriolis_En_Dis !< If CORIOLIS_EN_DIS is defined, two estimates of
82		!! the thickness fluxes are used to estimate the
83		!! Coriolis term, and the one that dissipates energy
84		!! relative to the other one is used. This is only
85		!! available at present if Coriolis scheme is
86		!! SADOURNY75_ENERGY.
87		logical :: weno_velocity_smooth !< If true, use velocity to compute the smoothness indicator for WENO
88		type(time_type), pointer :: Time !< A pointer to the ocean model's clock.
89		type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output.
90		!>@{ Diagnostic IDs
91		integer :: id_rv = -1, id_PV = -1, id_gKEu = -1, id_gKEv = -1
92		integer :: id_rvxu = -1, id_rvxv = -1
93		! integer :: id_hf_gKEu = -1, id_hf_gKEv = -1
94		integer :: id_hf_gKEu_2d = -1, id_hf_gKEv_2d = -1
95		integer :: id_intz_gKEu_2d = -1, id_intz_gKEv_2d = -1
96		! integer :: id_hf_rvxu = -1, id_hf_rvxv = -1
97		integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1
98		integer :: id_h_gKEu = -1, id_h_gKEv = -1
99		integer :: id_h_rvxu = -1, id_h_rvxv = -1
100		integer :: id_intz_rvxu_2d = -1, id_intz_rvxv_2d = -1
101		integer :: id_CAuS = -1, id_CAvS = -1
102		!>@}
103		end type CoriolisAdv_CS
104
105		!>@{ Enumeration values for Coriolis_Scheme
106		integer, parameter :: SADOURNY75_ENERGY = 1
107		integer, parameter :: ARAKAWA_HSU90 = 2
108		integer, parameter :: ROBUST_ENSTRO = 3
109		integer, parameter :: SADOURNY75_ENSTRO = 4
110		integer, parameter :: ARAKAWA_LAMB81 = 5
111		integer, parameter :: AL_BLEND = 6
112		integer, parameter :: wenovi7th_PV_ENSTRO = 7
113		integer, parameter :: wenovi5th_PV_ENSTRO = 8
114		integer, parameter :: wenovi3rd_PV_ENSTRO = 9
115		character*(20), parameter :: SADOURNY75_ENERGY_STRING = "SADOURNY75_ENERGY"
116		character*(20), parameter :: ARAKAWA_HSU_STRING = "ARAKAWA_HSU90"
117		character*(20), parameter :: ROBUST_ENSTRO_STRING = "ROBUST_ENSTRO"
118		character*(20), parameter :: SADOURNY75_ENSTRO_STRING = "SADOURNY75_ENSTRO"
119		character*(20), parameter :: ARAKAWA_LAMB_STRING = "ARAKAWA_LAMB81"
120		character*(20), parameter :: AL_BLEND_STRING = "ARAKAWA_LAMB_BLEND"
121		character*(20), parameter :: WENOVI7TH_PV_ENSTRO_STRING = "WENOVI7TH_PV_ENSTRO"
122		character*(20), parameter :: WENOVI5TH_PV_ENSTRO_STRING = "WENOVI5TH_PV_ENSTRO"
123		character*(20), parameter :: WENOVI3RD_PV_ENSTRO_STRING = "WENOVI3RD_PV_ENSTRO"
124		!>@}
125		!>@{ Enumeration values for KE_Scheme
126		integer, parameter :: KE_ARAKAWA = 10
127		integer, parameter :: KE_SIMPLE_GUDONOV = 11
128		integer, parameter :: KE_GUDONOV = 12
129		integer, parameter :: KE_UP3 = 13
130		character*(20), parameter :: KE_ARAKAWA_STRING = "KE_ARAKAWA"
131		character*(20), parameter :: KE_SIMPLE_GUDONOV_STRING = "KE_SIMPLE_GUDONOV"
132		character*(20), parameter :: KE_GUDONOV_STRING = "KE_GUDONOV"
133		character*(20), parameter :: KE_UP3_STRING = "KE_UP3"
134		!>@}
135		!>@{ Enumeration values for PV_Adv_Scheme
136		integer, parameter :: PV_ADV_CENTERED = 21
137		integer, parameter :: PV_ADV_UPWIND1 = 22
138		character*(20), parameter :: PV_ADV_CENTERED_STRING = "PV_ADV_CENTERED"
139		character*(20), parameter :: PV_ADV_UPWIND1_STRING = "PV_ADV_UPWIND1"
140		!>@}
141
142		contains
143
144		!> Calculates the Coriolis and momentum advection contributions to the acceleration.
145	49	subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv, Waves)
146		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
147		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
148		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
149		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
150		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
151		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: uh !< Zonal transport uhdy
152		!! [H L2 T-1 ~> m3 s-1 or kg s-1]
153		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vh !< Meridional transport vhdx
154		!! [H L2 T-1 ~> m3 s-1 or kg s-1]
155		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: CAu !< Zonal acceleration due to Coriolis
156		!! and momentum advection [L T-2 ~> m s-2].
157		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: CAv !< Meridional acceleration due to Coriolis
158		!! and momentum advection [L T-2 ~> m s-2].
159		type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
160		type(accel_diag_ptrs), intent(inout) :: AD !< Storage for acceleration diagnostics
161		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
162		type(CoriolisAdv_CS), intent(in) :: CS !< Control structure for MOM_CoriolisAdv
163		type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
164		type(Wave_parameters_CS), optional, pointer :: Waves !< An optional pointer to Stokes drift CS
165
166		! Local variables
167		real, dimension(SZIB_(G),SZJB_(G)) :: &
168	98	Area_q ! The sum of the ocean areas at the 4 adjacent thickness points [L2 ~> m2].
169		real, dimension(SZIB_(G),SZJB_(G),merge(GV%ke,CS%nkblock,CS%nkblock==0)) :: &
170	49	q, & ! Layer potential vorticity [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
171	49	qS, & ! Layer Stokes vorticity [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
172	49	Ih_q, & ! The inverse of thickness interpolated to q points [H-1 ~> m-1 or m2 kg-1].
173	49	h_q, & ! The thickness interpolated to q points [H-1 ~> m-1 or m2 kg-1].
174	98	dvdx, dudy, & ! Contributions to the circulation around q-points [L2 T-1 ~> m2 s-1]
175	98	dvSdx, duSdy, & ! idem. for Stokes drift [L2 T-1 ~> m2 s-1]
176	49	rel_vort, & ! Relative vorticity at q-points [T-1 ~> s-1].
177	49	abs_vort, & ! Absolute vorticity at q-points [T-1 ~> s-1].
178	49	stk_vort, & ! Stokes vorticity at q-points [T-1 ~> s-1].
179	49	q2 ! Relative vorticity over thickness [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
180
181		real, dimension(SZIB_(G),SZJ_(G),merge(GV%ke,CS%nkblock,CS%nkblock==0)) :: &
182	147	a, b, c, d, & ! a, b, c, & d are combinations of the potential vorticities
183		! surrounding an h grid point. At small scales, a = q/4,
184		! b = q/4, etc. All are in [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1],
185		! and use the indexing of the corresponding u point.
186	49	hArea_u, & ! The cell area weighted thickness interpolated to u points
187		! times the effective areas [H L2 ~> m3 or kg].
188	49	KEx, & ! The zonal gradient of Kinetic energy per unit mass [L T-2 ~> m s-2],
189		! KEx = d/dx KE.
190	49	uh_center ! Transport based on arithmetic mean h at u-points [H L2 T-1 ~> m3 s-1 or kg s-1]
191
192		real, dimension(SZI_(G),SZJ_(G)) :: &
193	49	Area_h ! The ocean area at h points [L2 ~> m2]. Area_h is used to find the
194		! average thickness in the denominator of q. 0 for land points.
195		real, dimension(SZI_(G),SZJ_(G),merge(GV%ke,CS%nkblock,CS%nkblock==0)) :: &
196	49	KE, & ! Kinetic energy per unit mass [L2 T-2 ~> m2 s-2], KE = (u^2 + v^2)/2.
197	49	uh_min, uh_max, & ! The smallest and largest estimates of the zonal volume fluxes through
198		! the faces (i.e. uhdy) [H L2 T-1 ~> m3 s-1 or kg s-1]
199	49	vh_min, vh_max, & ! The smallest and largest estimates of the meridional volume fluxes through
200		! the faces (i.e. vhdx) [H L2 T-1 ~> m3 s-1 or kg s-1]
201	49	ep_u, ep_v ! Additional pseudo-Coriolis terms in the Arakawa and Lamb
202		! discretization [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
203		real, dimension(SZI_(G),SZJB_(G),merge(GV%ke,CS%nkblock,CS%nkblock==0)) :: &
204	49	hArea_v, & ! The cell area weighted thickness interpolated to v points
205		! times the effective areas [H L2 ~> m3 or kg].
206	49	KEy, & ! The meridional gradient of Kinetic energy per unit mass [L T-2 ~> m s-2],
207		! KEy = d/dy KE.
208	49	vh_center ! Transport based on arithmetic mean h at v-points [H L2 T-1 ~> m3 s-1 or kg s-1]
209		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
210	49	PV, & ! A diagnostic array of the potential vorticities [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
211	49	RV ! A diagnostic array of the relative vorticities [T-1 ~> s-1].
212	98	real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: CAuS ! Stokes contribution to CAu [L T-2 ~> m s-2]
213	49	real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: CAvS ! Stokes contribution to CAv [L T-2 ~> m s-2]
214		real :: fv1, fv2, fv3, fv4 ! (f+rv)*v at the 4 points surrounding a u points[L T-2 ~> m s-2]
215		real :: fu1, fu2, fu3, fu4 ! -(f+rv)*u at the 4 points surrounding a v point [L T-2 ~> m s-2]
216		real :: max_fv, max_fu ! The maximum of the neighboring Coriolis accelerations [L T-2 ~> m s-2]
217		real :: min_fv, min_fu ! The minimum of the neighboring Coriolis accelerations [L T-2 ~> m s-2]
218
219		real, parameter :: C1_12 = 1.0 / 12.0 ! C1_12 = 1/12 [nondim]
220		real, parameter :: C1_24 = 1.0 / 24.0 ! C1_24 = 1/24 [nondim]
221		real :: max_Ihq, min_Ihq ! The maximum and minimum of the nearby Ihq [H-1 ~> m-1 or m2 kg-1].
222		real :: hArea_q ! The sum of area times thickness of the cells
223		! surrounding a q point [H L2 ~> m3 or kg].
224		real :: vol_neglect ! A volume so small that is expected to be
225		! lost in roundoff [H L2 ~> m3 or kg].
226		real :: area_neglect ! An area so small that is expected to be
227		! lost in roundoff [L2 ~> m2].
228		real :: temp1, temp2 ! Temporary variables [L2 T-2 ~> m2 s-2].
229		real :: eps_vel ! A tiny, positive velocity [L T-1 ~> m s-1].
230
231		real :: uhc, vhc ! Centered estimates of uh and vh [H L2 T-1 ~> m3 s-1 or kg s-1].
232		real :: uhm, vhm ! The input estimates of uh and vh [H L2 T-1 ~> m3 s-1 or kg s-1].
233		real :: c1, c2, c3, slope ! Nondimensional parameters for the Coriolis limiter scheme [nondim]
234
235		real :: Fe_m2 ! Temporary variable associated with the ARAKAWA_LAMB_BLEND scheme [nondim]
236		real :: rat_lin ! Temporary variable associated with the ARAKAWA_LAMB_BLEND scheme [nondim]
237		real :: rat_m1 ! The ratio of the maximum neighboring inverse thickness
238		! to the minimum inverse thickness minus 1 [nondim]. rat_m1 >= 0.
239		real :: AL_wt ! The relative weight of the Arakawa & Lamb scheme to the
240		! Arakawa & Hsu scheme [nondim], between 0 and 1.
241		real :: Sad_wt ! The relative weight of the Sadourny energy scheme to
242		! the other two with the ARAKAWA_LAMB_BLEND scheme [nondim],
243		! between 0 and 1.
244
245		real :: Heff1, Heff2 ! Temporary effective H at U or V points [H ~> m or kg m-2].
246		real :: Heff3, Heff4 ! Temporary effective H at U or V points [H ~> m or kg m-2].
247		real :: h_tiny ! A very small thickness [H ~> m or kg m-2].
248		real :: UHeff, VHeff ! More temporary variables [H L2 T-1 ~> m3 s-1 or kg s-1].
249		real :: QUHeff,QVHeff ! More temporary variables [H L2 T-2 ~> m3 s-2 or kg s-2].
250		integer :: i, j, k, n, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkblock, k_start, k_end, kmax, kk
251		integer :: Is_q, Ie_q, Js_q, Je_q ! The scheme-dependent range of values at which vorticity is set.
252		logical :: Stokes_VF
253		real :: u_v, v_u ! u_v is the u velocity at v point, v_u is the v velocity at u point [L T-1 ~> m s-1]
254		real :: q_v, q_u ! PV at the u and v points [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1]
255		logical :: use_weno ! True if using one of the WENO schemes
256		integer :: seventh_order, fifth_order, third_order ! Order of accuracy for the WENO calculations
257		real :: u_q8(8) ! Eight-point zonal velocity at WENO stencils [L T-1 ~> m s-1]
258		real :: u_q6(6) ! Six-point zonal velocity at WENO stencils [L T-1 ~> m s-1]
259		real :: u_q4(4) ! Four-point zonal velocity at WENO stencils [L T-1 ~> m s-1]
260		real :: v_q8(8) ! Eight-point meridional velocity at WENO stencils [L T-1 ~> m s-1]
261		real :: v_q6(6) ! Six-point meridional velocity at WENO stencils [L T-1 ~> m s-1]
262		real :: v_q4(4) ! Four-point meridional velocity at WENO stencils [L T-1 ~> m s-1]
263		integer :: stencil ! Stencil size of WENO scheme
264
265		! To work, the following fields must be set outside of the usual
266		! is to ie range before this subroutine is called:
267		! v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2),
268		! uh(is-1,ie,js:je+1) and vh(is:ie+1,js-1:je).
269
270	49	if (.not.CS%initialized) call MOM_error(FATAL, &
271	0	"MOM_CoriolisAdv: Module must be initialized before it is used.")
272
273	49	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
274	49	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
275	49	nkblock = merge(GV%ke, CS%nkblock, CS%nkblock==0)
276	49	vol_neglect = GV%H_subroundoff * (1e-4 * US%m_to_L)**2
277	49	area_neglect = (1e-4 * US%m_to_L)**2
278	49	eps_vel = 1.0e-10*US%m_s_to_L_T
279	49	h_tiny = GV%Angstrom_H ! Perhaps this should be set to h_neglect instead.
280
281	49	stencil = CoriolisAdv_stencil(CS)
282
283		use_weno = CS%Coriolis_Scheme == wenovi7th_PV_ENSTRO &
284		.or. CS%Coriolis_Scheme == wenovi5th_PV_ENSTRO &
285	49	.or. CS%Coriolis_Scheme == wenovi3rd_PV_ENSTRO
286
287	49	if (use_weno) then
288	0	Is_q = is - stencil ; Ie_q = ie + stencil - 1 ; Js_q = js - stencil ; Je_q = je + stencil - 1
289		else
290	49	Is_q = G%IscB - 1 ; Ie_q = G%IecB + 1 ; Js_q = G%JscB - 1 ; Je_q = G%JecB + 1
291		endif
292
293		!$omp target enter data map(alloc: Area_h, Area_q)
294
295	6125	do concurrent (j=Js_q:Je_q+1, I=Is_q:Ie_q+1)
296	394989	Area_h(i,j) = G%mask2dT(i,j) * G%areaT(i,j)
297		enddo
298
299	49	if (associated(OBC)) then
300		!$omp target update from(Area_h)
301
302	0	do n=1,OBC%number_of_segments
303	0	if (.not. OBC%segment(n)%on_pe) cycle
304	0	I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
305	0	if (OBC%segment(n)%is_N_or_S .and. (J >= Js_q) .and. (J <= Je_q)) then
306	0	do i = max(Is_q,OBC%segment(n)%HI%isd), min(Ie_q+1,OBC%segment(n)%HI%ied)
307	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
308	0	Area_h(i,j+1) = Area_h(i,j)
309		else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
310	0	Area_h(i,j) = Area_h(i,j+1)
311		endif
312		enddo
313	0	elseif (OBC%segment(n)%is_E_or_W .and. (I >= Is_q) .and. (I <= Ie_q)) then
314	0	do j = max(Js_q,OBC%segment(n)%HI%jsd), min(Je_q+1,OBC%segment(n)%HI%jed)
315	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
316	0	Area_h(i+1,j) = Area_h(i,j)
317		else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
318	0	Area_h(i,j) = Area_h(i+1,j)
319		endif
320		enddo
321		endif
322		enddo
323
324		!$omp target update to(Area_h)
325		endif
326
327	49	do concurrent (J=Js_q:Je_q, I=Is_q:Ie_q)
328		Area_q(i,j) = (Area_h(i,j) + Area_h(i+1,j+1)) + &
329	385777	(Area_h(i+1,j) + Area_h(i,j+1))
330		enddo
331
332	49	Stokes_VF = .false.
333	49	if (present(Waves)) then ; if (associated(Waves)) then
334	48	Stokes_VF = Waves%Stokes_VF
335		endif ; endif
336
337		!$omp target enter data map(alloc: dvdx, dudy)
338		!$omp target enter data map(alloc: hArea_u, hArea_v)
339		!$omp target enter data map(alloc: rel_vort, abs_vort, q, Ih_q)
340		!$omp target enter data map(alloc: h_q) if (use_weno)
341		!$omp target enter data map(alloc: a, b, c, d, ep_u, ep_v)
342		!$omp target enter data map(alloc: KE, KEx, KEy)
343		! TODO: These Stokes_VF fields seem associated with diagnostics
344		!$omp target enter data map(alloc: dvSdx, duSdy, stk_vort, qS) if (Stokes_VF)
345		!$omp target enter data map(alloc: CAuS, CAvS) if (Stokes_VF)
346		!$omp target enter data map(alloc: uh_center, vh_center) if (CS%Coriolis_En_Dis)
347		! TODO: May also need SADOURNEY75_ENERGY
348		!$omp target enter data map(alloc: uh_min, vh_min) if (CS%Coriolis_En_Dis)
349		!$omp target enter data map(alloc: uh_max, vh_max) if (CS%Coriolis_En_Dis)
350
351		! Diagnostics
352		!$omp target enter data map(alloc: RV) if (CS%id_RV > 0)
353		!$omp target enter data map(alloc: PV) if (CS%id_PV > 0)
354		!$omp target enter data map(alloc: q2) &
355		!$omp if(associated(AD%rv_x_u) .or. associated(AD%rv_x_v))
356		!$omp target enter data map(alloc: AD%gradKEu) if (associated(AD%gradKEu))
357		!$omp target enter data map(alloc: AD%gradKEv) if (associated(AD%gradKEv))
358		!$omp target enter data map(alloc: AD%rv_x_u) if (associated(AD%rv_x_u))
359		!$omp target enter data map(alloc: AD%rv_x_v) if (associated(AD%rv_x_v))
360
361		! TODO: Do this outside of the function
362		!$omp target enter data map(to: pbv, pbv%por_face_areaU, pbv%por_face_areaV) &
363		!$omp if (CS%Coriolis_En_Dis)
364
365		! Hoist AL_BLEND k-independent scalars out of the block loop.
366	49	Fe_m2 = 0.0 ; rat_lin = 0.0
367	49	if (CS%Coriolis_Scheme == AL_BLEND) then
368	0	Fe_m2 = CS%F_eff_max_blend - 2.0
369	0	rat_lin = 1.5 * Fe_m2 / max(CS%wt_lin_blend, 1.0e-16)
370	0	if (CS%F_eff_max_blend <= 2.0) then ; Fe_m2 = -1. ; rat_lin = -1.0 ; endif
371		endif
372
373	3724	do k_start=1,nz,nkblock
374	3675	k_end = min(k_start+nkblock-1, nz)
375	3675	kmax = k_end - k_start + 1
376		! Here the second order accurate layer potential vorticities, q,
377		! are calculated. hq is second order accurate in space. Relative
378		! vorticity is second order accurate everywhere with free slip b.c.s,
379		! but only first order accurate at boundaries with no slip b.c.s.
380		! First calculate the contributions to the circulation around the q-point.
381	3675	if (Stokes_VF) then
382	0	if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then
383	0	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q) DO_LOCALITY(local(k))
384	0	k = k_start + kk - 1
385		dvSdx(I,J,kk) = (-Waves%us_y(i+1,J,k)*G%dyCv(i+1,J)) - &
386	0	(-Waves%us_y(i,J,k)*G%dyCv(i,J))
387		duSdy(I,J,kk) = (-Waves%us_x(I,j+1,k)*G%dxCu(I,j+1)) - &
388	0	(-Waves%us_x(I,j,k)*G%dxCu(I,j))
389		enddo
390		endif
391	0	if (.not. Waves%Passive_Stokes_VF) then
392	0	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q) DO_LOCALITY(local(k))
393	0	k = k_start + kk - 1
394		dvdx(I,J,kk) = ((v(i+1,J,k)-Waves%us_y(i+1,J,k))*G%dyCv(i+1,J)) - &
395	0	((v(i,J,k)-Waves%us_y(i,J,k))*G%dyCv(i,J))
396		dudy(I,J,kk) = ((u(I,j+1,k)-Waves%us_x(I,j+1,k))*G%dxCu(I,j+1)) - &
397	0	((u(I,j,k)-Waves%us_x(I,j,k))*G%dxCu(I,j))
398		enddo
399		else
400	0	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q) DO_LOCALITY(local(k))
401	0	k = k_start + kk - 1
402	0	dvdx(I,J,kk) = (v(i+1,J,k)G%dyCv(i+1,J)) - (v(i,J,k)G%dyCv(i,J))
403	0	dudy(I,J,kk) = (u(I,j+1,k)G%dxCu(I,j+1)) - (u(I,j,k)G%dxCu(I,j))
404		enddo
405		endif
406		else
407	3675	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q) DO_LOCALITY(local(k))
408	28477575	k = k_start + kk - 1
409	28477575	dvdx(I,J,kk) = (v(i+1,J,k)G%dyCv(i+1,J)) - (v(i,J,k)G%dyCv(i,J))
410	57410850	dudy(I,J,kk) = (u(I,j+1,k)G%dxCu(I,j+1)) - (u(I,j,k)G%dxCu(I,j))
411		enddo
412		endif
413
414	3675	do concurrent (kk=1:kmax, J=Js_q:Je_q, i=Is_q:Ie_q+1) DO_LOCALITY(local(k))
415	28709100	k = k_start + kk - 1
416	57877575	hArea_v(i,J,kk) = 0.5((Area_h(i,j) h(i,j,k)) + (Area_h(i,j+1) * h(i,j+1,k)))
417		enddo
418
419	455700	do concurrent (kk=1:kmax, j=Js_q:Je_q+1, I=Is_q:Ie_q) DO_LOCALITY(local(k))
420	28929600	k = k_start + kk - 1
421	58314900	hArea_u(I,j,kk) = 0.5((Area_h(i,j) h(i,j,k)) + (Area_h(i+1,j) * h(i+1,j,k)))
422		enddo
423
424	3675	if (CS%Coriolis_En_Dis) then
425	0	do concurrent (kk=1:kmax, J=Jsq:Jeq+1, I=is-1:ie) DO_LOCALITY(local(k))
426	0	k = k_start + kk - 1
427	0	uh_center(I,j,kk) = 0.5 * ((G%dy_Cu(I,j)pbv%por_face_areaU(I,j,k)) u(I,j,k)) * (h(i,j,k) + h(i+1,j,k))
428		enddo
429
430	0	do concurrent (kk=1:kmax, J=js-1:je, i=Isq:Ieq+1) DO_LOCALITY(local(k))
431	0	k = k_start + kk - 1
432	0	vh_center(i,J,kk) = 0.5 * ((G%dx_Cv(i,J)pbv%por_face_areaV(i,J,k)) v(i,J,k)) * (h(i,j,k) + h(i,j+1,k))
433		enddo
434		endif
435
436		! Adjust circulation components to relative vorticity and thickness projected onto
437		! velocity points on open boundaries.
438	3675	if (associated(OBC)) then
439		!$omp target update from(Area_h)
440	0	do k=k_start,k_end ! TODO: port (OBC GPU path not yet implemented)
441	0	kk = k - k_start + 1
442		!$omp target update from(dvdx(:,:,kk), dudy(:,:,kk))
443		!$omp target update from(hArea_u(:,:,kk), hArea_v(:,:,kk))
444		!$omp target update from(uh_center(:,:,kk), vh_center(:,:,kk))
445
446	0	do n=1,OBC%number_of_segments
447	0	if (.not. OBC%segment(n)%on_pe) cycle
448
449	0	I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
450
451	0	if (OBC%segment(n)%is_N_or_S .and. (J >= Js_q) .and. (J <= Je_q)) then
452	0	select case (OBC%vorticity_config)
453		case (OBC_VORTICITY_ZERO)
454	0	do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB
455	0	dvdx(I,J,kk) = 0. ; dudy(I,J,kk) = 0.
456		enddo
457		case (OBC_VORTICITY_FREESLIP)
458	0	do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB
459	0	dudy(I,J,kk) = 0.
460		enddo
461		case (OBC_VORTICITY_COMPUTED)
462	0	do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB
463	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
464	0	dudy(I,J,kk) = 2.0(OBC%segment(n)%tangential_vel(I,J,k) - u(I,j,k))G%dxCu(I,j)
465		else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
466	0	dudy(I,J,kk) = 2.0(u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))G%dxCu(I,j+1)
467		endif
468		enddo
469		case (OBC_VORTICITY_SPECIFIED)
470	0	do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB
471	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
472	0	dudy(I,J,kk) = OBC%segment(n)%tangential_grad(I,J,k)G%dxCu(I,j)G%dyBu(I,J)
473		else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
474	0	dudy(I,J,kk) = OBC%segment(n)%tangential_grad(I,J,k)G%dxCu(I,j+1)G%dyBu(I,J)
475		endif
476		enddo
477		end select
478
479		! Project thicknesses across OBC points with a no-gradient condition.
480	0	do i=max(Is_q,OBC%segment(n)%HI%isd), min(Ie_q+1,OBC%segment(n)%HI%ied)
481	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
482	0	hArea_v(i,J,kk) = 0.5 * (Area_h(i,j) + Area_h(i,j+1)) * h(i,j,k)
483		else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
484	0	hArea_v(i,J,kk) = 0.5 * (Area_h(i,j) + Area_h(i,j+1)) * h(i,j+1,k)
485		endif
486		enddo
487
488	0	if (CS%Coriolis_En_Dis) then
489	0	do i=max(Isq,OBC%segment(n)%HI%isd), min(Ieq+1,OBC%segment(n)%HI%ied)
490	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
491	0	vh_center(i,J,kk) = (G%dx_Cv(i,J)pbv%por_face_areaV(i,J,k)) v(i,J,k) * h(i,j,k)
492		else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
493	0	vh_center(i,J,kk) = (G%dx_Cv(i,J)pbv%por_face_areaV(i,J,k)) v(i,J,k) * h(i,j+1,k)
494		endif
495		enddo
496		endif
497	0	elseif (OBC%segment(n)%is_E_or_W .and. (I >= Is_q) .and. (I <= Ie_q)) then
498	0	select case (OBC%vorticity_config)
499		case (OBC_VORTICITY_ZERO)
500	0	do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB
501	0	dvdx(I,J,kk) = 0. ; dudy(I,J,kk) = 0.
502		enddo
503		case (OBC_VORTICITY_FREESLIP)
504	0	do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB
505	0	dvdx(I,J,kk) = 0.
506		enddo
507		case (OBC_VORTICITY_COMPUTED)
508	0	do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB
509	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
510	0	dvdx(I,J,kk) = 2.0(OBC%segment(n)%tangential_vel(I,J,k) - v(i,J,k))G%dyCv(i,J)
511		else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
512	0	dvdx(I,J,kk) = 2.0(v(i+1,J,k) - OBC%segment(n)%tangential_vel(I,J,k))G%dyCv(i+1,J)
513		endif
514		enddo
515		case (OBC_VORTICITY_SPECIFIED)
516	0	do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB
517	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
518	0	dvdx(I,J,kk) = OBC%segment(n)%tangential_grad(I,J,k)G%dyCv(i,J)G%dxBu(I,J)
519		else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
520	0	dvdx(I,J,kk) = OBC%segment(n)%tangential_grad(I,J,k)G%dyCv(i+1,J)G%dxBu(I,J)
521		endif
522		enddo
523		end select
524
525		! Project thicknesses across OBC points with a no-gradient condition.
526	0	do j=max(Js_q,OBC%segment(n)%HI%jsd), min(Je_q+1,OBC%segment(n)%HI%jed)
527	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
528	0	hArea_u(I,j,kk) = 0.5(Area_h(i,j) + Area_h(i+1,j)) h(i,j,k)
529		else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
530	0	hArea_u(I,j,kk) = 0.5(Area_h(i,j) + Area_h(i+1,j)) h(i+1,j,k)
531		endif
532		enddo
533
534	0	if (CS%Coriolis_En_Dis) then
535	0	do j=max(Jsq,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
536	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
537	0	uh_center(I,j,kk) = (G%dy_Cu(I,j)pbv%por_face_areaU(I,j,k)) u(I,j,k) * h(i,j,k)
538		else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
539	0	uh_center(I,j,kk) = (G%dy_Cu(I,j)pbv%por_face_areaU(I,j,k)) u(I,j,k) * h(i+1,j,k)
540		endif
541		enddo
542		endif
543		endif
544		enddo
545
546		! Now project thicknesses across cell-corner points in the OBCs. The two
547		! projections have to occur in sequence and can not be combined easily.
548	0	do n=1,OBC%number_of_segments
549	0	if (.not. OBC%segment(n)%on_pe) cycle
550	0	I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
551	0	if (OBC%segment(n)%is_N_or_S .and. (J >= Js_q) .and. (J <= Je_q)) then
552	0	do I = max(Is_q,OBC%segment(n)%HI%IsdB), min(Ie_q,OBC%segment(n)%HI%IedB)
553	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
554	0	if (Area_h(i,j) + Area_h(i+1,j) > 0.0) then
555		hArea_u(I,j+1,kk) = hArea_u(I,j,kk) * ((Area_h(i,j+1) + Area_h(i+1,j+1)) / &
556	0	(Area_h(i,j) + Area_h(i+1,j)))
557	0	else ; hArea_u(I,j+1,kk) = 0.0 ; endif
558		else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
559	0	if (Area_h(i,j+1) + Area_h(i+1,j+1) > 0.0) then
560		hArea_u(I,j,kk) = hArea_u(I,j+1,kk) * ((Area_h(i,j) + Area_h(i+1,j)) / &
561	0	(Area_h(i,j+1) + Area_h(i+1,j+1)))
562	0	else ; hArea_u(I,j,kk) = 0.0 ; endif
563		endif
564		enddo
565	0	elseif (OBC%segment(n)%is_E_or_W .and. (I >= Is_q) .and. (I <= Ie_q)) then
566	0	do J = max(Js_q,OBC%segment(n)%HI%JsdB), min(Je_q,OBC%segment(n)%HI%JedB)
567	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
568	0	if (Area_h(i,j) + Area_h(i,j+1) > 0.0) then
569		hArea_v(i+1,J,kk) = hArea_v(i,J,kk) * ((Area_h(i+1,j) + Area_h(i+1,j+1)) / &
570	0	(Area_h(i,j) + Area_h(i,j+1)))
571	0	else ; hArea_v(i+1,J,kk) = 0.0 ; endif
572		else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
573	0	hArea_v(i,J,kk) = 0.5 * (Area_h(i,j) + Area_h(i,j+1)) * h(i,j+1,k)
574	0	if (Area_h(i+1,j) + Area_h(i+1,j+1) > 0.0) then
575		hArea_v(i,J,kk) = hArea_v(i+1,J,kk) * ((Area_h(i,j) + Area_h(i,j+1)) / &
576	0	(Area_h(i+1,j) + Area_h(i+1,j+1)))
577	0	else ; hArea_v(i,J,kk) = 0.0 ; endif
578		endif
579		enddo
580		endif
581		enddo
582
583		!$omp target update to(dvdx(:,:,kk), dudy(:,:,kk))
584		!$omp target update to(hArea_u(:,:,kk), hArea_v(:,:,kk))
585		!$omp target update to(uh_center(:,:,kk), vh_center(:,:,kk))
586		enddo ! k=k_start,k_end OBC
587		endif
588
589	3675	if (CS%no_slip) then
590	0	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q)
591	0	rel_vort(I,J,kk) = (2.0 - G%mask2dBu(I,J)) * (dvdx(I,J,kk) - dudy(I,J,kk)) * G%IareaBu(I,J)
592		enddo
593
594	0	if (Stokes_VF) then
595	0	if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then
596	0	do concurrent (kk=1:kmax, J=Jsq-1:Jeq+1, I=Isq-1:Ieq+1)
597	0	stk_vort(I,J,kk) = (2.0 - G%mask2dBu(I,J)) * (dvSdx(I,J,kk) - duSdy(I,J,kk)) * G%IareaBu(I,J)
598		enddo
599		endif
600		endif
601		else
602	3675	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q)
603	57410850	rel_vort(I,J,kk) = G%mask2dBu(I,J) * (dvdx(I,J,kk) - dudy(I,J,kk)) * G%IareaBu(I,J)
604		enddo
605
606	3675	if (Stokes_VF) then
607	0	if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then
608	0	do concurrent (kk=1:kmax, J=Jsq-1:Jeq+1, I=Isq-1:Ieq+1)
609	0	stk_vort(I,J,kk) = (2.0 - G%mask2dBu(I,J)) * (dvSdx(I,J,kk) - duSdy(I,J,kk)) * G%IareaBu(I,J)
610		enddo
611		endif
612		endif
613		endif
614
615	3675	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q)
616	57410850	abs_vort(I,J,kk) = G%CoriolisBu(I,J) + rel_vort(I,J,kk)
617		enddo
618
619	3675	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q) DO_LOCALITY(local(hArea_q))
620	28477575	hArea_q = (hArea_u(I,j,kk) + hArea_u(I,j+1,kk)) + (hArea_v(i,J,kk) + hArea_v(i+1,J,kk))
621	28477575	Ih_q(I,J,kk) = Area_q(I,J) / (hArea_q + vol_neglect)
622	57410850	q(I,J,kk) = abs_vort(I,J,kk) * Ih_q(I,J,kk)
623		enddo
624
625		! NOTE: `h_q` is only used by WENO and was pulled out of the above loop to
626		! improve GPU performance, but it may need to be moved back.
627	3675	if (use_weno) then
628	0	do concurrent (kk=1:kmax, J=Js_q:Je_q, I=Is_q:Ie_q) DO_LOCALITY(local(hArea_q))
629	0	hArea_q = (hArea_u(I,j,kk) + hArea_u(I,j+1,kk)) + (hArea_v(i,J,kk) + hArea_v(i+1,J,kk))
630	0	h_q(I,J,kk) = hArea_q / max(Area_q(I,J), area_neglect)
631		enddo
632		endif
633
634	3675	if (Stokes_VF) then
635	0	if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then
636	0	do concurrent (kk=1:kmax, J=js-1:Jeq, I=is-1:Ieq)
637	0	qS(I,J,kk) = stk_vort(I,J,kk) * Ih_q(I,J,kk)
638		enddo
639		endif
640		endif
641
642	3675	if (CS%id_rv > 0) then
643	0	do concurrent (kk=1:kmax, J=Jsq-1:Jeq+1, I=Isq-1:Ieq+1) DO_LOCALITY(local(k))
644	0	k = k_start + kk - 1
645	0	RV(I,J,k) = rel_vort(I,J,kk)
646		enddo
647		endif
648
649	3675	if (CS%id_PV > 0) then
650	0	do concurrent (kk=1:kmax, J=Jsq-1:Jeq+1, I=Isq-1:Ieq+1) DO_LOCALITY(local(k))
651	0	k = k_start + kk - 1
652	0	PV(I,J,k) = q(I,J,kk)
653		enddo
654		endif
655
656	3675	if (associated(AD%rv_x_v) .or. associated(AD%rv_x_u)) then
657	0	do concurrent (kk=1:kmax, J=Jsq-1:Jeq+1, I=Isq-1:Ieq+1)
658	0	q2(I,J,kk) = rel_vort(I,J,kk) * Ih_q(I,J,kk)
659		enddo
660		endif
661
662		! a, b, c, and d are combinations of neighboring potential
663		! vorticities which form the Arakawa and Hsu vorticity advection
664		! scheme. All are defined at u grid points.
665
666	3675	if (CS%Coriolis_Scheme == ARAKAWA_HSU90) then
667	0	do concurrent (kk=1:kmax, j=Jsq:Jeq+1, I=is-1:Ieq)
668	0	a(I,j,kk) = (q(I,J,kk) + (q(I+1,J,kk) + q(I,J-1,kk))) * C1_12
669	0	d(I,j,kk) = ((q(I,J,kk) + q(I+1,J-1,kk)) + q(I,J-1,kk)) * C1_12
670		enddo
671
672	0	do concurrent (kk=1:kmax, j=Jsq:Jeq+1, I=Isq:Ieq)
673	0	b(I,j,kk) = (q(I,J,kk) + (q(I-1,J,kk) + q(I,J-1,kk))) * C1_12
674	0	c(I,j,kk) = ((q(I,J,kk) + q(I-1,J-1,kk)) + q(I,J-1,kk)) * C1_12
675		enddo
676	3675	elseif (CS%Coriolis_Scheme == ARAKAWA_LAMB81) then
677	0	do concurrent (kk=1:kmax, j=Jsq:Jeq+1, I=Isq:Ieq+1)
678	0	a(I-1,j,kk) = (2.0(q(I,J,kk) + q(I-1,J-1,kk)) + (q(I-1,J,kk) + q(I,J-1,kk))) C1_24
679	0	d(I-1,j,kk) = ((q(I,j,kk) + q(I-1,J-1,kk)) + 2.0(q(I-1,J,kk) + q(I,J-1,kk))) C1_24
680	0	b(I,j,kk) = ((q(I,J,kk) + q(I-1,J-1,kk)) + 2.0(q(I-1,J,kk) + q(I,J-1,kk))) C1_24
681	0	c(I,j,kk) = (2.0(q(I,J,kk) + q(I-1,J-1,kk)) + (q(I-1,J,kk) + q(I,J-1,kk))) C1_24
682	0	ep_u(i,j,kk) = ((q(I,J,kk) - q(I-1,J-1,kk)) + (q(I-1,J,kk) - q(I,J-1,kk))) * C1_24
683	0	ep_v(i,j,kk) = (-(q(I,J,kk) - q(I-1,J-1,kk)) + (q(I-1,J,kk) - q(I,J-1,kk))) * C1_24
684		enddo
685	3675	elseif (CS%Coriolis_Scheme == AL_BLEND) then
686		! Fe_m2 and rat_lin are k-independent; computed before the block loop.
687		do concurrent (kk=1:kmax, j=Jsq:Jeq+1, I=Isq:Ieq+1) &
688	0	DO_LOCALITY(local(min_Ihq, max_Ihq, rat_m1, AL_wt, Sad_wt))
689	0	min_Ihq = MIN(Ih_q(I-1,J-1,kk), Ih_q(I,J-1,kk), Ih_q(I-1,J,kk), Ih_q(I,J,kk))
690	0	max_Ihq = MAX(Ih_q(I-1,J-1,kk), Ih_q(I,J-1,kk), Ih_q(I-1,J,kk), Ih_q(I,J,kk))
691	0	rat_m1 = 1.0e15
692	0	if (max_Ihq < 1.0e15*min_Ihq) rat_m1 = max_Ihq / min_Ihq - 1.0
693		! The weights used here are designed to keep the effective Coriolis
694		! acceleration from any one point on its neighbors within a factor
695		! of F_eff_max. The minimum permitted value is 2 (the factor for
696		! Sadourny's energy conserving scheme).
697
698		! Determine the relative weights of Arakawa & Lamb vs. Arakawa and Hsu.
699	0	if (rat_m1 <= Fe_m2) then ; AL_wt = 1.0
700	0	elseif (rat_m1 < 1.5Fe_m2) then ; AL_wt = 3.0Fe_m2 / rat_m1 - 2.0
701	0	else ; AL_wt = 0.0 ; endif
702
703		! Determine the relative weights of Sadourny Energy vs. the other two.
704	0	if (rat_m1 <= 1.5*Fe_m2) then ; Sad_wt = 0.0
705	0	elseif (rat_m1 <= rat_lin) then
706	0	Sad_wt = 1.0 - (1.5*Fe_m2) / rat_m1
707	0	elseif (rat_m1 < 2.0*rat_lin) then
708	0	Sad_wt = 1.0 - (CS%wt_lin_blend / rat_lin) * (rat_m1 - 2.0*rat_lin)
709	0	else ; Sad_wt = 1.0 ; endif
710
711		a(I-1,j,kk) = Sad_wt * 0.25 * q(I-1,J,kk) + (1.0 - Sad_wt) * &
712		( ((2.0-AL_wt)* q(I-1,J,kk) + AL_wt*q(I,J-1,kk)) + &
713	0	2.0 * (q(I,J,kk) + q(I-1,J-1,kk)) ) * C1_24
714		d(I-1,j,kk) = Sad_wt * 0.25 * q(I-1,J-1,kk) + (1.0 - Sad_wt) * &
715		( ((2.0-AL_wt)* q(I-1,J-1,kk) + AL_wt*q(I,J,kk)) + &
716	0	2.0 * (q(I-1,J,kk) + q(I,J-1,kk)) ) * C1_24
717		b(I,j,kk) = Sad_wt * 0.25 * q(I,J,kk) + (1.0 - Sad_wt) * &
718		( ((2.0-AL_wt)* q(I,J,kk) + AL_wt*q(I-1,J-1,kk)) + &
719	0	2.0 * (q(I-1,J,kk) + q(I,J-1,kk)) ) * C1_24
720		c(I,j,kk) = Sad_wt * 0.25 * q(I,J-1,kk) + (1.0 - Sad_wt) * &
721		( ((2.0-AL_wt)* q(I,J-1,kk) + AL_wt*q(I-1,J,kk)) + &
722	0	2.0 * (q(I,J,kk) + q(I-1,J-1,kk)) ) * C1_24
723	0	ep_u(i,j,kk) = AL_wt * ((q(I,J,kk) - q(I-1,J-1,kk)) + (q(I-1,J,kk) - q(I,J-1,kk))) * C1_24
724	0	ep_v(i,j,kk) = AL_wt * (-(q(I,J,kk) - q(I-1,J-1,kk)) + (q(I-1,J,kk) - q(I,J-1,kk))) * C1_24
725		enddo
726		endif
727
728		! .and. SADOURNEY75_ENERGY ??
729	3675	if (CS%Coriolis_En_Dis) then
730		! c1 = 1.0-1.5*RANGE ; c2 = 1.0-RANGE ; c3 = 2.0 ; slope = 0.5
731	0	c1 = 1.0-1.5*0.5 ; c2 = 1.0-0.5 ; c3 = 2.0 ; slope = 0.5
732
733	0	do concurrent (kk=1:kmax, j=Jsq:Jeq+1, I=is-1:ie) DO_LOCALITY(local(k, uhc, uhm))
734	0	k = k_start + kk - 1
735	0	uhc = uh_center(I,j,kk)
736	0	uhm = uh(I,j,k)
737		! This sometimes matters with some types of open boundary conditions.
738	0	if (G%dy_Cu(I,j) == 0.0) uhc = uhm
739
740	0	if (abs(uhc) < 0.1*abs(uhm)) then
741	0	uhm = 10.0*uhc
742	0	elseif (abs(uhc) > c1*abs(uhm)) then
743	0	if (abs(uhc) < c2abs(uhm)) then ; uhc = (3.0uhc+(1.0-c23.0)uhm)
744	0	elseif (abs(uhc) <= c3*abs(uhm)) then ; uhc = uhm
745	0	else ; uhc = slopeuhc+(1.0-c3slope)*uhm
746		endif
747		endif
748
749	0	if (uhc > uhm) then
750	0	uh_min(I,j,kk) = uhm ; uh_max(I,j,kk) = uhc
751		else
752	0	uh_max(I,j,kk) = uhm ; uh_min(I,j,kk) = uhc
753		endif
754		enddo
755
756	0	do concurrent (kk=1:kmax, J=js-1:je, i=Isq:Ieq+1) DO_LOCALITY(local(k, vhc, vhm))
757	0	k = k_start + kk - 1
758	0	vhc = vh_center(i,J,kk)
759	0	vhm = vh(i,J,k)
760		! This sometimes matters with some types of open boundary conditions.
761	0	if (G%dx_Cv(i,J) == 0.0) vhc = vhm
762
763	0	if (abs(vhc) < 0.1*abs(vhm)) then
764	0	vhm = 10.0*vhc
765	0	elseif (abs(vhc) > c1*abs(vhm)) then
766	0	if (abs(vhc) < c2abs(vhm)) then ; vhc = (3.0vhc+(1.0-c23.0)vhm)
767	0	elseif (abs(vhc) <= c3*abs(vhm)) then ; vhc = vhm
768	0	else ; vhc = slopevhc+(1.0-c3slope)*vhm
769		endif
770		endif
771
772	0	if (vhc > vhm) then
773	0	vh_min(i,J,kk) = vhm ; vh_max(i,J,kk) = vhc
774		else
775	0	vh_max(i,J,kk) = vhm ; vh_min(i,J,kk) = vhc
776		endif
777		enddo
778		endif
779
780		! Calculate KE and the gradient of KE
781	3675	call gradKE(u, v, h, KE, KEx, KEy, k_start, k_end, nkblock, G, GV, US, CS)
782		! TODO: Can KE be removed from this function?
783
784		! Calculate the tendencies of zonal velocity due to the Coriolis
785		! force and momentum advection. On a Cartesian grid, this is
786		! CAu = q * vh - d(KE)/dx.
787	3675	if (CS%Coriolis_Scheme == SADOURNY75_ENERGY) then
788	0	if (CS%Coriolis_En_Dis) then
789		! Energy dissipating biased scheme, Hallberg 200x
790	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k, temp1, temp2))
791	0	k = k_start + kk - 1
792	0	if (q(I,J,kk)*u(I,j,k) == 0.0) then
793		temp1 = q(I,J,kk) * ( (vh_max(i,j,kk)+vh_max(i+1,j,kk)) &
794	0	+ (vh_min(i,j,kk)+vh_min(i+1,j,kk)) )*0.5
795	0	elseif (q(I,J,kk)*u(I,j,k) < 0.0) then
796	0	temp1 = q(I,J,kk) * (vh_max(i,j,kk)+vh_max(i+1,j,kk))
797		else
798	0	temp1 = q(I,J,kk) * (vh_min(i,j,kk)+vh_min(i+1,j,kk))
799		endif
800	0	if (q(I,J-1,kk)*u(I,j,k) == 0.0) then
801		temp2 = q(I,J-1,kk) * ( (vh_max(i,j-1,kk)+vh_max(i+1,j-1,kk)) &
802	0	+ (vh_min(i,j-1,kk)+vh_min(i+1,j-1,kk)) )*0.5
803	0	elseif (q(I,J-1,kk)*u(I,j,k) < 0.0) then
804	0	temp2 = q(I,J-1,kk) * (vh_max(i,j-1,kk)+vh_max(i+1,j-1,kk))
805		else
806	0	temp2 = q(I,J-1,kk) * (vh_min(i,j-1,kk)+vh_min(i+1,j-1,kk))
807		endif
808	0	CAu(I,j,k) = 0.25 * G%IdxCu(I,j) * (temp1 + temp2)
809		enddo
810		else
811		! Energy conserving scheme, Sadourny 1975
812	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
813	0	k = k_start + kk - 1
814		CAu(I,j,k) = 0.25 * &
815		((q(I,J,kk) * (vh(i+1,J,k) + vh(i,J,k))) + &
816	0	(q(I,J-1,kk) * (vh(i,J-1,k) + vh(i+1,J-1,k)))) * G%IdxCu(I,j)
817		enddo
818		endif
819	3675	elseif (CS%Coriolis_Scheme == SADOURNY75_ENSTRO) then
820	3675	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
821	26680500	k = k_start + kk - 1
822		CAu(I,j,k) = 0.125 * (G%IdxCu(I,j) * (q(I,J,kk) + q(I,J-1,kk))) * &
823	53809350	((vh(i+1,J,k) + vh(i,J,k)) + (vh(i,J-1,k) + vh(i+1,J-1,k)))
824		enddo
825		elseif ((CS%Coriolis_Scheme == ARAKAWA_HSU90) .or. &
826	0	(CS%Coriolis_Scheme == ARAKAWA_LAMB81) .or. &
827		(CS%Coriolis_Scheme == AL_BLEND)) then
828		! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
829	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
830	0	k = k_start + kk - 1
831		CAu(I,j,k) = (((a(I,j,kk) * vh(i+1,J,k)) + (c(I,j,kk) * vh(i,J-1,k))) + &
832	0	((b(I,j,kk) * vh(i,J,k)) + (d(I,j,kk) * vh(i+1,J-1,k)))) * G%IdxCu(I,j)
833		enddo
834	0	elseif (CS%Coriolis_Scheme == ROBUST_ENSTRO) then
835		! An enstrophy conserving scheme robust to vanishing layers
836		! Note: Heffs are in lieu of h_at_v that should be returned by the
837		! continuity solver. AJA
838		do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) &
839	0	DO_LOCALITY(local(k, Heff1, Heff2, Heff3, Heff4, VHeff, QVHeff))
840	0	k = k_start + kk - 1
841	0	Heff1 = abs(vh(i,J,k) * G%IdxCv(i,J)) / (eps_vel+abs(v(i,J,k)))
842	0	Heff1 = max(Heff1, min(h(i,j,k),h(i,j+1,k)))
843	0	Heff1 = min(Heff1, max(h(i,j,k),h(i,j+1,k)))
844	0	Heff2 = abs(vh(i,J-1,k) * G%IdxCv(i,J-1)) / (eps_vel+abs(v(i,J-1,k)))
845	0	Heff2 = max(Heff2, min(h(i,j-1,k),h(i,j,k)))
846	0	Heff2 = min(Heff2, max(h(i,j-1,k),h(i,j,k)))
847	0	Heff3 = abs(vh(i+1,J,k) * G%IdxCv(i+1,J)) / (eps_vel+abs(v(i+1,J,k)))
848	0	Heff3 = max(Heff3, min(h(i+1,j,k),h(i+1,j+1,k)))
849	0	Heff3 = min(Heff3, max(h(i+1,j,k),h(i+1,j+1,k)))
850	0	Heff4 = abs(vh(i+1,J-1,k) * G%IdxCv(i+1,J-1)) / (eps_vel+abs(v(i+1,J-1,k)))
851	0	Heff4 = max(Heff4, min(h(i+1,j-1,k),h(i+1,j,k)))
852	0	Heff4 = min(Heff4, max(h(i+1,j-1,k),h(i+1,j,k)))
853	0	if (CS%PV_Adv_Scheme == PV_ADV_CENTERED) then
854		CAu(I,j,k) = 0.5(abs_vort(I,J,kk)+abs_vort(I,J-1,kk)) &
855		((vh(i,J,k) + vh(i+1,J-1,k)) + (vh(i,J-1,k) + vh(i+1,J,k)) ) / &
856	0	(h_tiny + ((Heff1+Heff4) + (Heff2+Heff3)) ) * G%IdxCu(I,j)
857	0	elseif (CS%PV_Adv_Scheme == PV_ADV_UPWIND1) then
858	0	VHeff = ((vh(i,J,k) + vh(i+1,J-1,k)) + (vh(i,J-1,k) + vh(i+1,J,k)) )
859		QVHeff = 0.5( ((abs_vort(I,J,kk)+abs_vort(I,J-1,kk))VHeff) &
860	0	- ((abs_vort(I,J,kk)-abs_vort(I,J-1,kk))*abs(VHeff)) )
861	0	CAu(I,j,k) = (QVHeff / ( h_tiny + ((Heff1+Heff4) + (Heff2+Heff3)) ) ) * G%IdxCu(I,j)
862		endif
863		enddo
864	0	elseif (CS%Coriolis_Scheme == wenovi7th_PV_ENSTRO) then
865	0	do k=k_start,k_end ! TODO: port
866	0	kk = k - k_start + 1
867		!$omp target update from(u(:,:,k), vh(:,:,k), abs_vort(:,:,kk), h_q(:,:,kk), q(:,:,kk))
868	0	do j=js,je ; do I=Isq,Ieq
869	0	v_u = 0.25G%IdxCu(I,j)((vh(i+1,J,k) + vh(i,J,k)) + (vh(i,J-1,k) + vh(i+1,J-1,k)))
870		! check whether there is masked land points in the stencil
871		third_order = (G%mask2dCu(I,j-2) * G%mask2dCu(I,j-1) * G%mask2dCu(I,j) * &
872	0	G%mask2dCu(I,j+1) * G%mask2dCu(I,j+2))
873
874	0	fifth_order = third_order * G%mask2dCu(I,j-3) * G%mask2dCu(I,j+3)
875	0	seventh_order = fifth_order * G%mask2dCu(I,j-4) * G%mask2dCu(I,j+4)
876
877
878		! compute the masking to make sure that inland values are not used
879	0	if (seventh_order == 1) then
880		! all values are valid, we use seventh order reconstruction
881	0	u_q8(:) = (u(I,j-4:j+3,k) + u(I,j-3:j+4,k)) * 0.5
882		call weno_seven_h_weight_reconstruction(abs_vort(I,J-4:J+3,kk), &
883		h_q(I,J-4:J+3,kk), &
884		u_q8, &
885	0	GV%H_subroundoff, v_u, q_u, cs%weno_velocity_smooth)
886	0	CAu(I,j,k) = (q_u * v_u)
887
888	0	elseif (fifth_order == 1) then
889		! all values are valid, we use fifth order reconstruction
890	0	u_q6(:) = (u(I,j-3:j+2,k) + u(I,j-2:j+3,k)) * 0.5
891		call weno_five_h_weight_reconstruction(abs_vort(I,J-3:J+2,kk), &
892		h_q(I,J-3:J+2,kk), &
893		u_q6, &
894	0	GV%H_subroundoff, v_u, q_u, CS%weno_velocity_smooth)
895	0	CAu(I,j,k) = (q_u * v_u)
896
897	0	elseif (third_order == 1) then
898		! only the middle values are valid, we use third order reconstruction
899	0	u_q4(:) = (u(I,j-2:j+1,k) + u(I,j-1:j+2,k)) * 0.5
900		call weno_three_h_weight_reconstruction(abs_vort(I,J-2:J+1,kk), &
901		h_q(I,J-2:J+1,kk), &
902		u_q4, &
903	0	GV%H_subroundoff, v_u, q_u, CS%weno_velocity_smooth)
904	0	CAu(I,j,k) = (q_u * v_u)
905		else ! Upwind first order
906	0	if (v_u>0.) then
907	0	q_u = q(I,J-1,kk)
908		else
909	0	q_u = q(I,J,kk)
910		endif
911	0	CAu(I,j,k) = (q_u * v_u)
912
913		endif
914		enddo ; enddo
915		!$omp target update to(CAu(:,:,k))
916		enddo
917	0	elseif (CS%Coriolis_Scheme == wenovi5th_PV_ENSTRO) then
918	0	do k=k_start,k_end ! TODO: port
919	0	kk = k - k_start + 1
920		!$omp target update from(u(:,:,k), vh(:,:,k), abs_vort(:,:,kk), h_q(:,:,kk), q(:,:,kk))
921	0	do j=js,je ; do I=Isq,Ieq
922	0	v_u = 0.25G%IdxCu(I,j)((vh(i+1,J,k) + vh(i,J,k)) + (vh(i,J-1,k) + vh(i+1,J-1,k)))
923		third_order = (G%mask2dCu(I,j-2) * G%mask2dCu(I,j-1) * G%mask2dCu(I,j) * &
924	0	G%mask2dCu(I,j+1) * G%mask2dCu(I,j+2))
925
926	0	fifth_order = third_order * G%mask2dCu(I,j-3) * G%mask2dCu(I,j+3)
927
928	0	if (fifth_order == 1) then
929		! all values are valid, we use fifth order reconstruction
930	0	u_q6(:) = (u(I,j-3:j+2,k) + u(I,j-2:j+3,k)) * 0.5
931		call weno_five_h_weight_reconstruction(abs_vort(I,J-3:J+2,kk), &
932		h_q(I,J-3:J+2,kk), &
933		u_q6, &
934	0	GV%H_subroundoff, v_u, q_u, CS%weno_velocity_smooth)
935	0	CAu(I,j,k) = (q_u * v_u)
936
937	0	elseif (third_order == 1) then
938		! only the middle values are valid, we use third order reconstruction
939	0	u_q4(:) = (u(I,j-2:j+1,k) + u(I,j-1:j+2,k)) * 0.5
940		call weno_three_h_weight_reconstruction(abs_vort(I,J-2:J+1,kk), &
941		h_q(I,J-2:J+1,kk), &
942		u_q4, &
943	0	GV%H_subroundoff, v_u, q_u, CS%weno_velocity_smooth)
944	0	CAu(I,j,k) = (q_u * v_u)
945
946		else ! Upwind first order
947	0	if (v_u>0.) then
948	0	q_u = q(I,J-1,kk)
949		else
950	0	q_u = q(I,J,kk)
951		endif
952	0	CAu(I,j,k) = (q_u * v_u)
953		endif
954		enddo ; enddo
955		!$omp target update to(CAu(:,:,k))
956		enddo
957	0	elseif (CS%Coriolis_Scheme == wenovi3rd_PV_ENSTRO) then
958	0	do k=k_start,k_end ! TODO: port
959	0	kk = k - k_start + 1
960		!$omp target update from(u(:,:,k), vh(:,:,k), abs_vort(:,:,kk), h_q(:,:,kk), q(:,:,kk))
961	0	do j=js,je ; do I=Isq,Ieq
962	0	v_u = 0.25G%IdxCu(I,j)((vh(i+1,J,k) + vh(i,J,k)) + (vh(i,J-1,k) + vh(i+1,J-1,k)))
963		third_order = (G%mask2dCu(I,j-2) * G%mask2dCu(I,j-1) * G%mask2dCu(I,j) * &
964	0	G%mask2dCu(I,j+1) * G%mask2dCu(I,j+2))
965
966
967	0	if (third_order == 1) then
968		! only the middle values are valid, we use third order reconstruction
969	0	u_q4(:) = (u(I,j-2:j+1,k) + u(I,j-1:j+2,k)) * 0.5
970		call weno_three_h_weight_reconstruction(abs_vort(I,J-2:J+1,kk), &
971		h_q(I,J-2:J+1,kk), &
972		u_q4, &
973	0	GV%H_subroundoff, v_u, q_u, CS%weno_velocity_smooth)
974	0	CAu(I,j,k) = (q_u * v_u)
975
976		else ! Upwind first order
977	0	if (v_u>0.) then
978	0	q_u = q(I,J-1,kk)
979		else
980	0	q_u = q(I,J,kk)
981		endif
982	0	CAu(I,j,k) = (q_u * v_u)
983		endif
984		enddo ; enddo
985		!$omp target update to(CAu(:,:,k))
986		enddo
987		endif
988
989		! Add in the additional terms with Arakawa & Lamb.
990	3675	if ((CS%Coriolis_Scheme == ARAKAWA_LAMB81) .or. &
991		(CS%Coriolis_Scheme == AL_BLEND)) then
992	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
993	0	k = k_start + kk - 1
994		CAu(I,j,k) = CAu(I,j,k) + &
995	0	((ep_u(i,j,kk)uh(I-1,j,k)) - (ep_u(i+1,j,kk)uh(I+1,j,k))) * G%IdxCu(I,j)
996		enddo
997		endif
998
999	3675	if (Stokes_VF) then
1000	0	if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then
1001		! Computing the diagnostic Stokes contribution to CAu
1002	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
1003	0	k = k_start + kk - 1
1004		CAuS(I,j,k) = 0.25 * &
1005		((qS(I,J,kk) * (vh(i+1,J,k) + vh(i,J,k))) + &
1006	0	(qS(I,J-1,kk) * (vh(i,J-1,k) + vh(i+1,J-1,k)))) * G%IdxCu(I,j)
1007		enddo
1008		endif
1009		endif
1010
1011	3675	if (CS%bound_Coriolis) then
1012	3675	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k, fv1, fv2, fv3, fv4, max_fv, min_fv))
1013	26680500	k = k_start + kk - 1
1014	26680500	fv1 = abs_vort(I,J,kk) * v(i+1,J,k)
1015	26680500	fv2 = abs_vort(I,J,kk) * v(i,J,k)
1016	26680500	fv3 = abs_vort(I,J-1,kk) * v(i+1,J-1,k)
1017	26680500	fv4 = abs_vort(I,J-1,kk) * v(i,J-1,k)
1018
1019	26680500	max_fv = max(fv1, fv2, fv3, fv4)
1020	26680500	min_fv = min(fv1, fv2, fv3, fv4)
1021
1022	26680500	CAu(I,j,k) = min(CAu(I,j,k), max_fv)
1023	53809350	CAu(I,j,k) = max(CAu(I,j,k), min_fv)
1024		enddo
1025		endif
1026
1027		! Term - d(KE)/dx.
1028	3675	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
1029	26680500	k = k_start + kk - 1
1030	53809350	CAu(I,j,k) = CAu(I,j,k) - KEx(I,j,kk)
1031		enddo
1032
1033	3675	if (associated(AD%gradKEu)) then
1034	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
1035	0	k = k_start + kk - 1
1036	0	AD%gradKEu(I,j,k) = -KEx(I,j,kk)
1037		enddo
1038		endif
1039
1040		! Calculate the tendencies of meridional velocity due to the Coriolis
1041		! force and momentum advection. On a Cartesian grid, this is
1042		! CAv = - q * uh - d(KE)/dy.
1043	3675	if (CS%Coriolis_Scheme == SADOURNY75_ENERGY) then
1044	0	if (CS%Coriolis_En_Dis) then
1045		! Energy dissipating biased scheme, Hallberg 200x
1046	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k, temp1, temp2))
1047	0	k = k_start + kk - 1
1048	0	if (q(I-1,J,kk)*v(i,J,k) == 0.0) then
1049		temp1 = q(I-1,J,kk) * ( (uh_max(i-1,j,kk)+uh_max(i-1,j+1,kk)) &
1050	0	+ (uh_min(i-1,j,kk)+uh_min(i-1,j+1,kk)) )*0.5
1051	0	elseif (q(I-1,J,kk)*v(i,J,k) > 0.0) then
1052	0	temp1 = q(I-1,J,kk) * (uh_max(i-1,j,kk)+uh_max(i-1,j+1,kk))
1053		else
1054	0	temp1 = q(I-1,J,kk) * (uh_min(i-1,j,kk)+uh_min(i-1,j+1,kk))
1055		endif
1056	0	if (q(I,J,kk)*v(i,J,k) == 0.0) then
1057		temp2 = q(I,J,kk) * ( (uh_max(i,j,kk)+uh_max(i,j+1,kk)) &
1058	0	+ (uh_min(i,j,kk)+uh_min(i,j+1,kk)) )*0.5
1059	0	elseif (q(I,J,kk)*v(i,J,k) > 0.0) then
1060	0	temp2 = q(I,J,kk) * (uh_max(i,j,kk)+uh_max(i,j+1,kk))
1061		else
1062	0	temp2 = q(I,J,kk) * (uh_min(i,j,kk)+uh_min(i,j+1,kk))
1063		endif
1064	0	CAv(i,J,k) = -0.25 * G%IdyCv(i,J) * (temp1 + temp2)
1065		enddo
1066		else
1067		! Energy conserving scheme, Sadourny 1975
1068	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1069	0	k = k_start + kk - 1
1070		CAv(i,J,k) = - 0.25* &
1071		((q(I-1,J,kk)*(uh(I-1,j,k) + uh(I-1,j+1,k))) + &
1072	0	(q(I,J,kk)(uh(I,j,k) + uh(I,j+1,k)))) G%IdyCv(i,J)
1073		enddo
1074		endif
1075	3675	elseif (CS%Coriolis_Scheme == SADOURNY75_ENSTRO) then
1076	3675	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1077	26901000	k = k_start + kk - 1
1078		CAv(i,J,k) = -0.125 * (G%IdyCv(i,J) * (q(I-1,J,kk) + q(I,J,kk))) * &
1079	54246675	((uh(I-1,j,k) + uh(I-1,j+1,k)) + (uh(I,j,k) + uh(I,j+1,k)))
1080		enddo
1081		elseif ((CS%Coriolis_Scheme == ARAKAWA_HSU90) .or. &
1082	0	(CS%Coriolis_Scheme == ARAKAWA_LAMB81) .or. &
1083		(CS%Coriolis_Scheme == AL_BLEND)) then
1084		! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
1085	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1086	0	k = k_start + kk - 1
1087		CAv(i,J,k) = - (((a(I-1,j,kk) * uh(I-1,j,k)) + &
1088		(c(I,j+1,kk) * uh(I,j+1,k))) &
1089		+ ((b(I,j,kk) * uh(I,j,k)) + &
1090	0	(d(I-1,j+1,kk) * uh(I-1,j+1,k)))) * G%IdyCv(i,J)
1091		enddo
1092	0	elseif (CS%Coriolis_Scheme == ROBUST_ENSTRO) then
1093		! An enstrophy conserving scheme robust to vanishing layers
1094		! Note: Heffs are in lieu of h_at_u that should be returned by the
1095		! continuity solver. AJA
1096		do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) &
1097	0	DO_LOCALITY(local(k, Heff1, Heff2, Heff3, Heff4, UHeff, QUHeff))
1098	0	k = k_start + kk - 1
1099	0	Heff1 = abs(uh(I,j,k) * G%IdyCu(I,j)) / (eps_vel+abs(u(I,j,k)))
1100	0	Heff1 = max(Heff1, min(h(i,j,k),h(i+1,j,k)))
1101	0	Heff1 = min(Heff1, max(h(i,j,k),h(i+1,j,k)))
1102	0	Heff2 = abs(uh(I-1,j,k) * G%IdyCu(I-1,j)) / (eps_vel+abs(u(I-1,j,k)))
1103	0	Heff2 = max(Heff2, min(h(i-1,j,k),h(i,j,k)))
1104	0	Heff2 = min(Heff2, max(h(i-1,j,k),h(i,j,k)))
1105	0	Heff3 = abs(uh(I,j+1,k) * G%IdyCu(I,j+1)) / (eps_vel+abs(u(I,j+1,k)))
1106	0	Heff3 = max(Heff3, min(h(i,j+1,k),h(i+1,j+1,k)))
1107	0	Heff3 = min(Heff3, max(h(i,j+1,k),h(i+1,j+1,k)))
1108	0	Heff4 = abs(uh(I-1,j+1,k) * G%IdyCu(I-1,j+1)) / (eps_vel+abs(u(I-1,j+1,k)))
1109	0	Heff4 = max(Heff4, min(h(i-1,j+1,k),h(i,j+1,k)))
1110	0	Heff4 = min(Heff4, max(h(i-1,j+1,k),h(i,j+1,k)))
1111	0	if (CS%PV_Adv_Scheme == PV_ADV_CENTERED) then
1112		CAv(i,J,k) = - 0.5(abs_vort(I,J,kk)+abs_vort(I-1,J,kk)) &
1113		((uh(I ,j ,k)+uh(I-1,j+1,k)) + &
1114		(uh(I-1,j ,k)+uh(I ,j+1,k)) ) / &
1115	0	(h_tiny + ((Heff1+Heff4) +(Heff2+Heff3)) ) * G%IdyCv(i,J)
1116	0	elseif (CS%PV_Adv_Scheme == PV_ADV_UPWIND1) then
1117		UHeff = ((uh(I ,j ,k)+uh(I-1,j+1,k)) + &
1118	0	(uh(I-1,j ,k)+uh(I ,j+1,k)) )
1119		QUHeff = 0.5( ((abs_vort(I,J,kk)+abs_vort(I-1,J,kk))UHeff) &
1120	0	- ((abs_vort(I,J,kk)-abs_vort(I-1,J,kk))*abs(UHeff)) )
1121		CAv(i,J,k) = - QUHeff / &
1122	0	(h_tiny + ((Heff1+Heff4) +(Heff2+Heff3)) ) * G%IdyCv(i,J)
1123		endif
1124		enddo
1125	0	elseif (CS%Coriolis_Scheme == wenovi7th_PV_ENSTRO) then
1126	0	do k=k_start,k_end ! TODO: port
1127	0	kk = k - k_start + 1
1128		!$omp target update from(v(:,:,k), uh(:,:,k), abs_vort(:,:,kk), h_q(:,:,kk), q(:,:,kk))
1129	0	do J=Jsq,Jeq ; do i=is,ie
1130	0	u_v = 0.25G%IdyCv(i,J)((uh(I-1,j,k) + uh(I-1,j+1,k)) + (uh(I,j,k) + uh(I,j+1,k)))
1131
1132		! check whether there is any masked land values within the stencils
1133		third_order = (G%mask2dCv(i-2,J) * G%mask2dCv(i-1,J) * G%mask2dCv(i,J) * G%mask2dCv(i+1,J) * &
1134	0	G%mask2dCv(i+2,J))
1135	0	fifth_order = third_order * G%mask2dCv(i-3,J) * G%mask2dCv(i+3,J)
1136	0	seventh_order = fifth_order * G%mask2dCv(i-4,J) * G%mask2dCv(i+4,J)
1137
1138
1139
1140		! compute the masking to make sure that inland values are not used
1141	0	if (seventh_order == 1) then
1142	0	v_q8(:) = (v(i-4:i+3,J,k) + v(i-3:i+4,J,k)) * 0.5
1143		! all values are valid, we use seventh order reconstruction
1144		call weno_seven_h_weight_reconstruction(abs_vort(I-4:I+3,J,kk), &
1145		h_q(I-4:I+3,J,kk), &
1146		v_q8, &
1147	0	GV%H_subroundoff, u_v, q_v, CS%weno_velocity_smooth)
1148	0	CAv(i,J,k) = - (q_v * u_v)
1149
1150	0	elseif (fifth_order == 1) then
1151	0	v_q6(:) = (v(i-3:i+2,J,k) + v(i-2:i+3,J,k)) * 0.5
1152		! all values are valid, we use fifth order reconstruction
1153		call weno_five_h_weight_reconstruction(abs_vort(I-3:I+2,J,kk), &
1154		h_q(I-3:I+2,J,kk), &
1155		v_q6, &
1156	0	GV%H_subroundoff, u_v, q_v, CS%weno_velocity_smooth)
1157	0	CAv(i,J,k) = - (q_v * u_v)
1158
1159	0	elseif (third_order == 1) then
1160	0	v_q4(:) = (v(i-2:i+1,J,k) + v(i-1:i+2,J,k)) * 0.5
1161		! ! only the middle values are valid, we use third order reconstruction
1162		call weno_three_h_weight_reconstruction(abs_vort(I-2:I+1,J,kk), &
1163		h_q(I-2:I+1,J,kk), &
1164		v_q4, &
1165	0	GV%H_subroundoff, u_v, q_v, CS%weno_velocity_smooth)
1166	0	CAv(i,J,k) = - (q_v * u_v)
1167		else ! Upwind first order!
1168	0	if (u_v>0.) then
1169	0	q_v = q(I-1,J,kk)
1170		else
1171	0	q_v = q(I,J,kk)
1172		endif
1173	0	CAv(i,J,k) = - (q_v * u_v)
1174		endif
1175
1176		enddo ; enddo
1177		!$omp target update to(CAv(:,:,k))
1178		enddo
1179	0	elseif (CS%Coriolis_Scheme == wenovi5th_PV_ENSTRO) then
1180	0	do k=k_start,k_end ! TODO: port
1181	0	kk = k - k_start + 1
1182		!$omp target update from(v(:,:,k), uh(:,:,k), abs_vort(:,:,kk), h_q(:,:,kk), q(:,:,kk))
1183	0	do J=Jsq,Jeq ; do i=is,ie
1184	0	u_v = 0.25G%IdyCv(i,J)((uh(I-1,j,k) + uh(I-1,j+1,k)) + (uh(I,j,k) + uh(I,j+1,k)))
1185
1186		third_order = (G%mask2dCv(i-2,J) * G%mask2dCv(i-1,J) * G%mask2dCv(i,J) * G%mask2dCv(i+1,J) * &
1187	0	G%mask2dCv(i+2,J))
1188	0	fifth_order = third_order * G%mask2dCv(i-3,J) * G%mask2dCv(i+3,J)
1189
1190
1191		! compute the masking to make sure that inland values are not used
1192	0	if (fifth_order == 1) then
1193	0	v_q6(:) = (v(i-3:i+2,J,k) + v(i-2:i+3,J,k)) * 0.5
1194		! all values are valid, we use fifth order reconstruction
1195		call weno_five_h_weight_reconstruction(abs_vort(I-3:I+2,J,kk), &
1196		h_q(I-3:I+2,J,kk), &
1197		v_q6, &
1198	0	GV%H_subroundoff, u_v, q_v, CS%weno_velocity_smooth)
1199	0	CAv(i,J,k) = - (q_v * u_v)
1200
1201	0	elseif (third_order == 1) then
1202	0	v_q4(:) = (v(i-2:i+1,J,k) + v(i-1:i+2,J,k)) * 0.5
1203		! ! only the middle values are valid, we use third order reconstruction
1204		call weno_three_h_weight_reconstruction(abs_vort(I-2:I+1,J,kk), &
1205		h_q(I-2:I+1,J,kk), &
1206		v_q4, &
1207	0	GV%H_subroundoff, u_v, q_v, CS%weno_velocity_smooth)
1208	0	CAv(i,J,k) = - (q_v * u_v)
1209
1210		else
1211	0	if (u_v>0.) then
1212	0	q_v = q(I-1,J,kk)
1213		else
1214	0	q_v = q(I,J,kk)
1215		endif
1216	0	CAv(i,J,k) = - (q_v * u_v)
1217		endif
1218
1219		enddo ; enddo
1220		!$omp target update to(CAv(:,:,k))
1221		enddo
1222	0	elseif (CS%Coriolis_Scheme == wenovi3rd_PV_ENSTRO) then
1223	0	do k=k_start,k_end ! TODO: port
1224	0	kk = k - k_start + 1
1225		!$omp target update from(v(:,:,k), uh(:,:,k), abs_vort(:,:,kk), h_q(:,:,kk), q(:,:,kk))
1226	0	do J=Jsq,Jeq ; do i=is,ie
1227	0	u_v = 0.25G%IdyCv(i,J)((uh(I-1,j,k) + uh(I-1,j+1,k)) + (uh(I,j,k) + uh(I,j+1,k)))
1228
1229		third_order = (G%mask2dCv(i-2,J) * G%mask2dCv(i-1,J) * G%mask2dCv(i,J) * G%mask2dCv(i+1,J) * &
1230	0	G%mask2dCv(i+2,J))
1231
1232
1233		! compute the masking to make sure that inland values are not used
1234	0	if (third_order == 1) then
1235	0	v_q4(:) = (v(i-2:i+1,J,k) + v(i-1:i+2,J,k)) * 0.5
1236		! ! only the middle values are valid, we use third order reconstruction
1237		call weno_three_h_weight_reconstruction(abs_vort(I-2:I+1,J,kk), &
1238		h_q(I-2:I+1,J,kk), &
1239		v_q4, &
1240	0	GV%H_subroundoff, u_v, q_v, CS%weno_velocity_smooth)
1241	0	CAv(i,J,k) = - (q_v * u_v)
1242
1243		else
1244	0	if (u_v>0.) then
1245	0	q_v = q(I-1,J,kk)
1246		else
1247	0	q_v = q(I,J,kk)
1248		endif
1249	0	CAv(i,J,k) = - (q_v * u_v)
1250		endif
1251
1252		enddo ; enddo
1253		!$omp target update to(CAv(:,:,k))
1254		enddo
1255		endif
1256		! Add in the additonal terms with Arakawa & Lamb.
1257	3675	if ((CS%Coriolis_Scheme == ARAKAWA_LAMB81) .or. &
1258		(CS%Coriolis_Scheme == AL_BLEND)) then
1259	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1260	0	k = k_start + kk - 1
1261		CAv(i,J,k) = CAv(i,J,k) + &
1262	0	((ep_v(i,j,kk)vh(i,J-1,k)) - (ep_v(i,j+1,kk)vh(i,J+1,k))) * G%IdyCv(i,J)
1263		enddo
1264		endif
1265
1266	3675	if (Stokes_VF) then
1267	0	if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then
1268		! Computing the diagnostic Stokes contribution to CAv
1269	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1270	0	k = k_start + kk - 1
1271		CAvS(i,J,k) = 0.25 * &
1272		((qS(I,J,kk) * (uh(I,j+1,k) + uh(I,j,k))) + &
1273	0	(qS(I-1,J,kk) * (uh(I-1,j,k) + uh(I-1,j+1,k)))) * G%IdyCv(i,J)
1274		enddo
1275		endif
1276		endif
1277
1278	3675	if (CS%bound_Coriolis) then
1279	3675	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k, fu1, fu2, fu3, fu4, max_fu, min_fu))
1280	26901000	k = k_start + kk - 1
1281	26901000	fu1 = -abs_vort(I,J,kk) * u(I,j+1,k)
1282	26901000	fu2 = -abs_vort(I,J,kk) * u(I,j,k)
1283	26901000	fu3 = -abs_vort(I-1,J,kk) * u(I-1,j+1,k)
1284	26901000	fu4 = -abs_vort(I-1,J,kk) * u(I-1,j,k)
1285
1286	26901000	max_fu = max(fu1, fu2, fu3, fu4)
1287	26901000	min_fu = min(fu1, fu2, fu3, fu4)
1288
1289	26901000	CAv(I,j,k) = min(CAv(I,j,k), max_fu)
1290	54246675	CAv(I,j,k) = max(CAv(I,j,k), min_fu)
1291		enddo
1292		endif
1293
1294		! Term - d(KE)/dy.
1295	3675	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1296	26901000	k = k_start + kk - 1
1297	54246675	CAv(i,J,k) = CAv(i,J,k) - KEy(i,J,kk)
1298		enddo
1299	3675	if (associated(AD%gradKEv)) then
1300	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1301	0	k = k_start + kk - 1
1302	0	AD%gradKEv(i,J,k) = -KEy(i,J,kk)
1303		enddo
1304		endif
1305
1306	3724	if (associated(AD%rv_x_u) .or. associated(AD%rv_x_v)) then
1307		! Calculate the Coriolis-like acceleration due to relative vorticity.
1308	0	if (CS%Coriolis_Scheme == SADOURNY75_ENERGY) then
1309	0	if (associated(AD%rv_x_u)) then
1310	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1311	0	k = k_start + kk - 1
1312		AD%rv_x_u(i,J,k) = - 0.25* &
1313		((q2(I-1,j,kk)*(uh(I-1,j,k) + uh(I-1,j+1,k))) + &
1314	0	(q2(I,j,kk)(uh(I,j,k) + uh(I,j+1,k)))) G%IdyCv(i,J)
1315		enddo
1316		endif
1317
1318	0	if (associated(AD%rv_x_v)) then
1319	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
1320	0	k = k_start + kk - 1
1321		AD%rv_x_v(I,j,k) = 0.25 * &
1322		((q2(I,j,kk) * (vh(i+1,J,k) + vh(i,J,k))) + &
1323	0	(q2(I,j-1,kk) * (vh(i,J-1,k) + vh(i+1,J-1,k)))) * G%IdxCu(I,j)
1324		enddo
1325		endif
1326		else
1327	0	if (associated(AD%rv_x_u)) then
1328	0	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(local(k))
1329	0	k = k_start + kk - 1
1330		AD%rv_x_u(i,J,k) = -G%IdyCv(i,J) * C1_12 * &
1331		(((((q2(I,J,kk) + q2(I-1,J-1,kk)) + q2(I-1,J,kk)) * uh(I-1,j,k)) + &
1332		(((q2(I-1,J,kk) + q2(I,J+1,kk)) + q2(I,J,kk)) * uh(I,j+1,k))) + &
1333		((((q2(I-1,J,kk) + q2(I,J-1,kk)) + q2(I,J,kk)) * uh(I,j,k))+ &
1334	0	(((q2(I,J,kk) + q2(I-1,J+1,kk)) + q2(I-1,J,kk)) * uh(I-1,j+1,k))))
1335		enddo
1336		endif
1337
1338	0	if (associated(AD%rv_x_v)) then
1339	0	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq) DO_LOCALITY(local(k))
1340	0	k = k_start + kk - 1
1341		AD%rv_x_v(I,j,k) = G%IdxCu(I,j) * C1_12 * &
1342		(((((q2(I+1,J,kk) + q2(I,J-1,kk)) + q2(I,J,kk)) * vh(i+1,J,k)) + &
1343		(((q2(I-1,J-1,kk) + q2(I,J,kk)) + q2(I,J-1,kk)) * vh(i,J-1,k))) + &
1344		((((q2(I-1,J,kk) + q2(I,J-1,kk)) + q2(I,J,kk)) * vh(i,J,k)) + &
1345	0	(((q2(I+1,J-1,kk) + q2(I,J,kk)) + q2(I,J-1,kk)) * vh(i+1,J-1,k))))
1346		enddo
1347		endif
1348		endif
1349		endif
1350		enddo ! end of k_start block loop.
1351
1352		!$omp target exit data map(delete: Area_h, Area_q)
1353		!$omp target exit data map(delete: dvdx, dudy)
1354		!$omp target exit data map(delete: hArea_u, hArea_v)
1355		!$omp target exit data map(delete: rel_vort, abs_vort, q, Ih_q)
1356		!$omp target exit data map(delete: h_q) if (use_weno)
1357		!$omp target exit data map(delete: a, b, c, d, ep_u, ep_v)
1358		!$omp target exit data map(delete: KE, KEx, KEy)
1359		!$omp target exit data map(delete: dvSdx, duSdy, stk_vort, qS) if (Stokes_VF)
1360		!$omp target exit data map(delete: uh_center, vh_center) if (CS%Coriolis_En_Dis)
1361		!$omp target exit data map(delete: uh_min, vh_min) if (CS%Coriolis_En_Dis)
1362		!$omp target exit data map(delete: uh_max, vh_max) if (CS%Coriolis_En_Dis)
1363		!$omp target exit data map(delete: q2) &
1364		!$omp if(associated(AD%rv_x_u) .or. associated(AD%rv_x_v))
1365
1366		! TODO: Move outside function
1367		!$omp target exit data map(delete: pbv, pbv%por_face_areaU, pbv%por_face_areaV) &
1368		!$omp if (CS%Coriolis_En_Dis)
1369
1370		! Diagnostics
1371		!$omp target exit data map(from: RV) if (CS%id_RV > 0)
1372		!$omp target exit data map(from: PV) if (CS%id_RV > 0)
1373		!$omp target exit data map(from: AD%gradKEu) if (associated(AD%gradKEu))
1374		!$omp target exit data map(from: AD%gradKEv) if (associated(AD%gradKEv))
1375		!$omp target exit data map(from: AD%rv_x_u) if (associated(AD%rv_x_u))
1376		!$omp target exit data map(from: AD%rv_x_u) if (associated(AD%rv_x_u))
1377		!$omp target exit data map(from: CAuS, CAvS) if (Stokes_VF)
1378
1379		! Here the various Coriolis-related derived quantities are offered for averaging.
1380	49	if (query_averaging_enabled(CS%diag)) then
1381	24	if (CS%id_rv > 0) call post_data(CS%id_rv, RV, CS%diag)
1382	24	if (CS%id_PV > 0) call post_data(CS%id_PV, PV, CS%diag)
1383	24	if (CS%id_gKEu>0) call post_data(CS%id_gKEu, AD%gradKEu, CS%diag)
1384	24	if (CS%id_gKEv>0) call post_data(CS%id_gKEv, AD%gradKEv, CS%diag)
1385	24	if (CS%id_rvxu > 0) call post_data(CS%id_rvxu, AD%rv_x_u, CS%diag)
1386	24	if (CS%id_rvxv > 0) call post_data(CS%id_rvxv, AD%rv_x_v, CS%diag)
1387	24	if (Stokes_VF) then
1388	0	if (CS%id_CAuS > 0) call post_data(CS%id_CAuS, CAuS, CS%diag)
1389	0	if (CS%id_CAvS > 0) call post_data(CS%id_CAvS, CAvS, CS%diag)
1390		endif
1391
1392		! Diagnostics for terms multiplied by fractional thicknesses
1393
1394		! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option.
1395		! The code is retained for debugging purposes in the future.
1396		! if (CS%id_hf_gKEu > 0) call post_product_u(CS%id_hf_gKEu, AD%gradKEu, AD%diag_hfrac_u, G, nz, CS%diag)
1397		! if (CS%id_hf_gKEv > 0) call post_product_v(CS%id_hf_gKEv, AD%gradKEv, AD%diag_hfrac_v, G, nz, CS%diag)
1398		! if (CS%id_hf_rvxv > 0) call post_product_u(CS%id_hf_rvxv, AD%rv_x_v, AD%diag_hfrac_u, G, nz, CS%diag)
1399		! if (CS%id_hf_rvxu > 0) call post_product_v(CS%id_hf_rvxu, AD%rv_x_u, AD%diag_hfrac_v, G, nz, CS%diag)
1400
1401	24	if (CS%id_hf_gKEu_2d > 0) call post_product_sum_u(CS%id_hf_gKEu_2d, AD%gradKEu, AD%diag_hfrac_u, G, nz, CS%diag)
1402	24	if (CS%id_hf_gKEv_2d > 0) call post_product_sum_v(CS%id_hf_gKEv_2d, AD%gradKEv, AD%diag_hfrac_v, G, nz, CS%diag)
1403	24	if (CS%id_intz_gKEu_2d > 0) call post_product_sum_u(CS%id_intz_gKEu_2d, AD%gradKEu, AD%diag_hu, G, nz, CS%diag)
1404	24	if (CS%id_intz_gKEv_2d > 0) call post_product_sum_v(CS%id_intz_gKEv_2d, AD%gradKEv, AD%diag_hv, G, nz, CS%diag)
1405
1406	24	if (CS%id_hf_rvxv_2d > 0) call post_product_sum_u(CS%id_hf_rvxv_2d, AD%rv_x_v, AD%diag_hfrac_u, G, nz, CS%diag)
1407	24	if (CS%id_hf_rvxu_2d > 0) call post_product_sum_v(CS%id_hf_rvxu_2d, AD%rv_x_u, AD%diag_hfrac_v, G, nz, CS%diag)
1408
1409	24	if (CS%id_h_gKEu > 0) call post_product_u(CS%id_h_gKEu, AD%gradKEu, AD%diag_hu, G, nz, CS%diag)
1410	24	if (CS%id_h_gKEv > 0) call post_product_v(CS%id_h_gKEv, AD%gradKEv, AD%diag_hv, G, nz, CS%diag)
1411	24	if (CS%id_h_rvxv > 0) call post_product_u(CS%id_h_rvxv, AD%rv_x_v, AD%diag_hu, G, nz, CS%diag)
1412	24	if (CS%id_h_rvxu > 0) call post_product_v(CS%id_h_rvxu, AD%rv_x_u, AD%diag_hv, G, nz, CS%diag)
1413
1414	24	if (CS%id_intz_rvxv_2d > 0) call post_product_sum_u(CS%id_intz_rvxv_2d, AD%rv_x_v, AD%diag_hu, G, nz, CS%diag)
1415	24	if (CS%id_intz_rvxu_2d > 0) call post_product_sum_v(CS%id_intz_rvxu_2d, AD%rv_x_u, AD%diag_hv, G, nz, CS%diag)
1416		endif
1417	49	end subroutine CorAdCalc
1418
1419
1420		!> Calculates the acceleration due to the gradient of kinetic energy for a k-block.
1421	3675	subroutine gradKE(u, v, h, KE, KEx, KEy, k_start, k_end, nkblock, G, GV, US, CS)
1422		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1423		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1424		integer, intent(in) :: k_start !< First layer in the k-block [nondim].
1425		integer, intent(in) :: k_end !< Last layer in the k-block [nondim].
1426		integer, intent(in) :: nkblock !< Size of the k-block [nondim].
1427		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
1428		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
1429		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1430		real, dimension(SZI_(G),SZJ_(G),nkblock), intent(out) :: KE !< Kinetic energy per unit mass [L2 T-2 ~> m2 s-2]
1431		real, dimension(SZIB_(G),SZJ_(G),nkblock), intent(out) :: KEx !< Zonal acceleration due to kinetic
1432		!! energy gradient [L T-2 ~> m s-2]
1433		real, dimension(SZI_(G),SZJB_(G),nkblock), intent(out) :: KEy !< Meridional acceleration due to kinetic
1434		!! energy gradient [L T-2 ~> m s-2]
1435		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1436		type(CoriolisAdv_CS), intent(in) :: CS !< Control structure for MOM_CoriolisAdv
1437		! Local variables
1438		real :: um, up, vm, vp ! Temporary variables [L T-1 ~> m s-1].
1439		real :: um2, up2, vm2, vp2 ! Temporary variables [L2 T-2 ~> m2 s-2].
1440		real :: um2a, up2a, vm2a, vp2a ! Temporary variables [L4 T-2 ~> m4 s-2].
1441		real :: third_order_u, third_order_v ! Product of mask values to determine the boundary
1442		integer :: i, j, k, kk, kmax, is, ie, js, je, Isq, Ieq, Jsq, Jeq
1443		real, parameter :: C1_12 = 1.0/12.0 ! The ratio of 1/12 [nondim]
1444
1445	3675	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1446	3675	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1447	3675	kmax = k_end - k_start + 1
1448
1449		! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term).
1450	3675	if (CS%KE_Scheme == KE_ARAKAWA) then
1451		! The following calculation of Kinetic energy includes the metric terms
1452		! identified in Arakawa & Lamb 1982 as important for KE conservation. It
1453		! also includes the possibility of partially-blocked tracer cell faces.
1454	452025	do concurrent (kk=1:kmax, j=Jsq:Jeq+1, i=Isq:Ieq+1) DO_LOCALITY(local(k))
1455	27797700	k = k_start + kk - 1
1456		KE(i,j,kk) = ( ( (G%areaCu( I ,j)(u( I ,j,k)u( I ,j,k))) + &
1457		(G%areaCu(I-1,j)(u(I-1,j,k)u(I-1,j,k))) ) + &
1458		( (G%areaCv(i, J )(v(i, J ,k)v(i, J ,k))) + &
1459	56047425	(G%areaCv(i,J-1)(v(i,J-1,k)v(i,J-1,k))) ) )0.25G%IareaT(i,j)
1460		enddo
1461	0	elseif (CS%KE_Scheme == KE_SIMPLE_GUDONOV) then
1462		! The following discretization of KE is based on the one-dimensional Gudonov
1463		! scheme which does not take into account any geometric factors
1464		do concurrent (kk=1:kmax, j=Jsq:Jeq+1, i=Isq:Ieq+1) &
1465	0	DO_LOCALITY(local(k, up, um, vp, vm, up2, um2, vp2, vm2))
1466	0	k = k_start + kk - 1
1467	0	up = 0.5( u(I-1,j,k) + ABS( u(I-1,j,k) ) ) ; up2 = upup
1468	0	um = 0.5( u( I ,j,k) - ABS( u( I ,j,k) ) ) ; um2 = umum
1469	0	vp = 0.5( v(i,J-1,k) + ABS( v(i,J-1,k) ) ) ; vp2 = vpvp
1470	0	vm = 0.5( v(i, J ,k) - ABS( v(i, J ,k) ) ) ; vm2 = vmvm
1471	0	KE(i,j,kk) = ( max(up2,um2) + max(vp2,vm2) ) *0.5
1472		enddo
1473	0	elseif (CS%KE_Scheme == KE_GUDONOV) then
1474		! The following discretization of KE is based on the one-dimensional Gudonov
1475		! scheme but has been adapted to take horizontal grid factors into account
1476		do concurrent (kk=1:kmax, j=Jsq:Jeq+1, i=Isq:Ieq+1) &
1477	0	DO_LOCALITY(local(k, up, um, vp, vm, up2a, um2a, vp2a, vm2a))
1478	0	k = k_start + kk - 1
1479	0	up = 0.5( u(I-1,j,k) + ABS( u(I-1,j,k) ) ) ; up2a = upup*G%areaCu(I-1,j)
1480	0	um = 0.5( u( I ,j,k) - ABS( u( I ,j,k) ) ) ; um2a = umum*G%areaCu( I ,j)
1481	0	vp = 0.5( v(i,J-1,k) + ABS( v(i,J-1,k) ) ) ; vp2a = vpvp*G%areaCv(i,J-1)
1482	0	vm = 0.5( v(i, J ,k) - ABS( v(i, J ,k) ) ) ; vm2a = vmvm*G%areaCv(i, J )
1483	0	KE(i,j,kk) = ( max(um2a,up2a) + max(vm2a,vp2a) )0.5G%IareaT(i,j)
1484		enddo
1485	0	elseif (CS%KE_Scheme == KE_UP3) then
1486		! The following discretization of KE is based on the one-dimensional third-order
1487		! upwind scheme which does not take horizontal grid factors into account
1488
1489	0	do k=k_start,k_end ! TODO: port
1490	0	kk = k - k_start + 1
1491		! TODO: GPU data tansfers?
1492	0	if (CS%KE_use_limiter) then
1493	0	do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1494		! compute the masking to make sure that inland values are not used
1495		third_order_u = (G%mask2dCu(I-2,j) * G%mask2dCu(I-1,j)* &
1496	0	G%mask2dCu(I,j) * G%mask2dCu(I+1,j))
1497
1498	0	if (third_order_u == 1) then
1499	0	up = (7.0 * (u(I-1,j,k) + u(I,j,k)) - (u(I-2,j,k) + u(I+1,j,k))) * C1_12
1500	0	call UP3_Koren_limiter_reconstruction(u(I-2:I+1,j,k), up, um)
1501		else
1502	0	up = (u(I-1,j,k) + u(I,j,k))*0.5
1503	0	if (up>0.) then
1504	0	um = u(I-1,j,k)
1505	0	elseif (up<0.) then
1506	0	um = u(I,j,k)
1507		else
1508	0	um = up
1509		endif
1510		endif
1511
1512		third_order_v = (G%mask2dCv(i,J-2) * G%mask2dCv(i,J-1)* &
1513	0	G%mask2dCv(i,J) * G%mask2dCv(i,J+1))
1514	0	if (third_order_v ==1) then
1515	0	vp = (7.0 * (v(i,J-1,k) + v(i,J,k)) - (v(i,J-2,k) + v(i,J+1,k))) * C1_12
1516	0	call UP3_Koren_limiter_reconstruction(v(i,J-2:J+1,k), vp, vm)
1517		else
1518	0	vp = (v(i,J-1,k) + v(i,J,k))*0.5
1519	0	if (vp>0.) then
1520	0	vm = v(i,J-1,k)
1521	0	elseif (vp<0.) then
1522	0	vm = v(i,J,k)
1523		else
1524	0	vm = vp
1525		endif
1526		endif
1527
1528	0	KE(i,j,kk) = ( (umum) + (vmvm) )*0.5
1529		enddo ; enddo
1530		else
1531	0	do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1532		! compute the masking to make sure that inland values are not used
1533		third_order_u = (G%mask2dCu(I-2,j) * G%mask2dCu(I-1,j)* &
1534	0	G%mask2dCu(I,j) * G%mask2dCu(I+1,j))
1535
1536	0	if (third_order_u == 1) then
1537	0	up = (7.0 * (u(I-1,j,k) + u(I,j,k)) - (u(I-2,j,k) + u(I+1,j,k))) * C1_12
1538	0	call UP3_reconstruction(u(I-2:I+1,j,k), up, um)
1539		else
1540	0	up = (u(I-1,j,k) + u(I,j,k))*0.5
1541	0	if (up>0.) then
1542	0	um = u(I-1,j,k)
1543	0	elseif (up<0.) then
1544	0	um = u(I,j,k)
1545		else
1546	0	um = up
1547		endif
1548		endif
1549
1550		third_order_v = (G%mask2dCv(i,J-2) * G%mask2dCv(i,J-1)* &
1551	0	G%mask2dCv(i,J) * G%mask2dCv(i,J+1))
1552	0	if (third_order_v ==1) then
1553	0	vp = (7.0 * (v(i,J-1,k) + v(i,J,k)) - (v(i,J-2,k) + v(i,J+1,k))) * C1_12
1554	0	call UP3_reconstruction(v(i,J-2:J+1,k), vp, vm)
1555		else
1556	0	vp = (v(i,J-1,k) + v(i,J,k))*0.5
1557	0	if (vp>0.) then
1558	0	vm = v(i,J-1,k)
1559	0	elseif (vp<0.) then
1560	0	vm = v(i,J,k)
1561		else
1562	0	vm = vp
1563		endif
1564		endif
1565
1566	0	KE(i,j,kk) = ( (umum) + (vmvm) )*0.5
1567		enddo ; enddo
1568		endif
1569		enddo
1570		endif
1571
1572		! Term - d(KE)/dx.
1573	3675	do concurrent (kk=1:kmax, j=js:je, I=Isq:Ieq)
1574	53809350	KEx(I,j,kk) = (KE(i+1,j,kk) - KE(i,j,kk)) * G%IdxCu_OBCmask(I,j)
1575		enddo
1576
1577		! Term - d(KE)/dy.
1578	3675	do concurrent (kk=1:kmax, J=Jsq:Jeq, i=is:ie)
1579	54246675	KEy(i,J,kk) = (KE(i,j+1,kk) - KE(i,j,kk)) * G%IdyCv_OBCmask(i,J)
1580		enddo
1581	3675	end subroutine gradKE
1582
1583		!> Reconstruct the scalar (e.g., pv, vorticity) onto point i-1/2 using a third-order upwind scheme
1584	0	subroutine UP3_reconstruction(q4,u,qr)
1585		real, intent(in) :: q4(4) !< Tracer values on points i-2, i-1, i, i+1 [A ~> a]
1586		real, intent(in) :: u !< Velocity or thickness flux on point i-1/2
1587		!! [l t-1 ~> m s-1] or [l2 t-1 ~> m2 s-1]
1588		real, intent(inout) :: qr !< Reconstruction of tracer q at point i-1/2 [A ~> a]
1589		real, parameter :: C1_6 = 1.0/6.0 ! The ratio of 1/6 [nondim]
1590
1591	0	if (u>0.) then
1592	0	qr = ((2.q4(3) + 5.q4(2)) - q4(1)) * C1_6
1593		else
1594	0	qr = ((2.q4(2) + 5.q4(3)) - q4(4)) * C1_6
1595		endif
1596
1597	0	end subroutine UP3_reconstruction
1598
1599
1600		!> Reconstruct the scalar (e.g., PV, vorticity) onto point i-1/2
1601		!! using a third-order upwind scheme with the Koren flux limiter
1602	0	subroutine UP3_Koren_limiter_reconstruction(q4,u,qr)
1603		real, intent(in) :: q4(4) !< Tracer values on points i-2, i-1, i, i+1 [A ~> a]
1604		real, intent(in) :: u !< Velocity or thickness flux on point i-1/2
1605		!! [L T-1 ~> m s-1] or [L2 T-1 ~> m2 s-1]
1606		real, intent(inout) :: qr !< Reconstruction of tracer q on point i-1/2 [A ~> a]
1607		real :: theta ! Ratio of gradient [nondim]
1608		real :: psi ! Limiter function [nondim]
1609		real, parameter :: C1_3 = 1.0/3.0 ! The ratio of 1/3 [nondim]
1610		real, parameter :: C1_6 = 1.0/6.0 ! The ratio of 1/6 [nondim]
1611
1612	0	if (u>0.) then
1613	0	if (q4(3) == q4(2)) then
1614	0	qr = q4(2)
1615		else
1616	0	theta = (q4(2) - q4(1))/(q4(3) - q4(2))
1617	0	psi = max(0., min(1., C1_3 + C1_6*theta, theta)) ! limiter introduced by Koren (1993)
1618	0	qr = q4(2) + psi*(q4(3) - q4(2))
1619		endif
1620		else
1621	0	if (q4(3) == q4(2)) then
1622	0	qr = q4(3)
1623		else
1624	0	theta = (q4(4) - q4(3))/(q4(3) - q4(2))
1625	0	psi = max(0., min(1., C1_3 + C1_6*theta, theta))
1626	0	qr = q4(3) + psi*(q4(2) - q4(3))
1627		endif
1628		endif
1629
1630	0	end subroutine UP3_Koren_limiter_reconstruction
1631
1632		!> Compute the factor for the WENO weights
1633	0	function fac_fn(tau, b) result(fac)
1634		real, intent(in) :: tau !< Difference of the smoothness indicator [A ~> a]
1635		real, intent(in) :: b !< The smoothness indicator [A ~> a]
1636		real :: fac !< The factor for the weight [nondim]
1637
1638	0	fac = 1.0e40 ; if (abs(b) > 1.0e-20tau) fac = (1 + tau / b)*2
1639
1640	0	end function fac_fn
1641
1642
1643		!> Reconstruct the tracer (e.g., PV, vorticity) onto the point i-1/2 using a third-order WENO scheme
1644		!! This reconstruction is thickness-weighted
1645	0	subroutine weno_three_h_weight_reconstruction(q4, h4, u4, &
1646		h_tiny, u, qr, velocity_smoothing)
1647		real, intent(in) :: q4(4) !< Tracer value times thickness on points i-2, i-1, i, i+1 [A ~> a]
1648		real, intent(in) :: h4(4) !< Thickness values on points i-2, i-1, i, i+1 [L ~> m]
1649		real, optional, intent(in) :: u4(4) !< Velocity values on points i-2, i-1, i, i+1
1650		!![L T-1 ~> m s-1]
1651		real, intent(in) :: h_tiny !< A tiny thickness to prevent division by zero [L ~> m]
1652		real, intent(in) :: u !< Velocity or thickness flux on point i-1/2
1653		!! [L T-1 ~> m s-1] or [L2 T-1 ~> m2 s-1]
1654		real, intent(inout) :: qr !< Reconstruction of tracer q on point i-1/2 [A ~> a]
1655		logical, intent(in) :: velocity_smoothing !< If true, use velocity to compute smoothness indicator
1656		real :: vr ! Reconstruction of hq [A ~> a]
1657		real :: hr ! Reconstruction of h [L ~> m]
1658		real :: c0, c1 ! Intermediate reconstruction of q [A ~> a]
1659		real :: d0, d1 ! Intermediate reconstruction of h [L ~> m]
1660		real :: b0, b1 ! Smoothness indicator [A ~> a]
1661		real :: tau ! Difference of smoothness indicator [A ~> a]
1662		real :: w0, w1 ! Weights [nondim]
1663		real :: s ! Temporary variables [nondim]
1664		real, parameter :: C2_3 = 2.0/3.0 ! The ratio of 2/3 [nondim]
1665		real, parameter :: C1_3 = 1.0/3.0 ! The ratio of 1/3 [nondim]
1666
1667	0	if (u>0.) then
1668	0	call weno_three_reconstruction_0(q4(2:3), c0) ! Reconstruction in the second upwind stencil
1669	0	call weno_three_reconstruction_1(q4(1:2), c1) ! Reconstruction in the first upwind stencil
1670
1671	0	call weno_three_reconstruction_0(h4(2:3), d0)
1672	0	call weno_three_reconstruction_1(h4(1:2), d1)
1673	0	if (velocity_smoothing) then
1674	0	call weno_three_weight(u4(2:3), b0) ! Smoothness indicator the second upwind stencil
1675	0	call weno_three_weight(u4(1:2), b1) ! Smoothness indicator the first upwind stencil
1676		else
1677	0	call weno_three_weight(q4(2:3), b0) ! Smoothness indicator the second upwind stencil
1678	0	call weno_three_weight(q4(1:2), b1) ! Smoothness indicator the first upwind stencil
1679		endif
1680		else
1681	0	call weno_three_reconstruction_0(q4(3:2:-1), c0) ! Reconstruction in the second upwind stencil
1682	0	call weno_three_reconstruction_1(q4(4:3:-1), c1) ! Reconstruction in the first upwind stencil
1683
1684	0	call weno_three_reconstruction_0(h4(3:2:-1), d0)
1685	0	call weno_three_reconstruction_1(h4(4:3:-1), d1)
1686	0	if (velocity_smoothing) then
1687	0	call weno_three_weight(u4(3:2:-1), b0) ! Smoothness indicator the second upwind stencil
1688	0	call weno_three_weight(u4(4:3:-1), b1) ! Smoothness indicator the first upwind stencil
1689		else
1690	0	call weno_three_weight(q4(3:2:-1), b0) ! Smoothness indicator the second upwind stencil
1691	0	call weno_three_weight(q4(4:3:-1), b1) ! Smoothness indicator the first upwind stencil
1692		endif
1693		endif
1694
1695	0	tau = abs(b0-b1)
1696	0	w0 = C2_3 * fac_fn(tau, b0)
1697	0	w1 = C1_3 * fac_fn(tau, b1)
1698
1699	0	s = 1. / (w0 + w1)
1700	0	w0 = w0 * s ! Weights of stencils
1701	0	w1 = w1 * s
1702
1703	0	vr = (w0 * c0) + (w1 * c1)
1704	0	hr = (w0 * d0) + (w1 * d1)
1705		! vr = min(max(q4(3), q4(2)), vr) ; vr = max(min(q4(3), q4(2)), vr) !Impose a monotonicity limiter
1706	0	hr = min(max(h4(3), h4(2)), hr) ; hr = max(min(h4(3), h4(2)), hr) ! A monotonicity limiter
1707
1708	0	qr = vr / max(hr, h_tiny)
1709
1710	0	end subroutine weno_three_h_weight_reconstruction
1711
1712		!> Compute the smoothness indicator for the two-point stencil of the third-order WENO scheme
1713	0	subroutine weno_three_weight(q2, w0)
1714		real, intent(in) :: q2(2) !< Tracer values on the two-point stencil [A ~> a]
1715		real, intent(inout) :: w0 !< Smoothness indicator for this stencil [A2 ~> a2]
1716
1717	0	w0 = (q2(1) - q2(2))**2
1718
1719	0	end subroutine weno_three_weight
1720
1721		!> Reconstruction in the second upwind stencil of the third-order WENO scheme
1722	0	subroutine weno_three_reconstruction_0(q2, w0)
1723		real, intent(in) :: q2(2) !< Tracer values on the two-point stencil [A ~> a]
1724		real, intent(inout) :: w0 !< Reconstruction of the quantity [A2 ~> a2]
1725
1726	0	w0 = (q2(1) + q2(2)) * 0.5
1727
1728	0	end subroutine weno_three_reconstruction_0
1729
1730		!> Reconstruction in the first upwind stencil for third-order WENO scheme
1731	0	subroutine weno_three_reconstruction_1(q2, w0)
1732		real, intent(in) :: q2(2) !< Tracer values on the two-point stencil [A ~> a]
1733		real, intent(inout) :: w0 !< Reconstruction of the quantity [A ~> a]
1734
1735	0	w0 = (- q2(1) + 3 * q2(2)) * 0.5
1736
1737	0	end subroutine weno_three_reconstruction_1
1738
1739
1740		!> Reconstruct the tracer (e.g., PV, vorticity) onto point i-1/2 using a fifth-order WENO scheme
1741		!! The reconstruction is weighted by the thickness
1742	0	subroutine weno_five_h_weight_reconstruction(q6, h6, u6, &
1743		h_tiny, u, qr, velocity_smoothing)
1744		real, intent(in) :: q6(6)
1745		!< Tracer values on points i-3, i-2, i-1, i, i+1, i+2 [A ~> a]
1746		real, intent(in) :: h6(6)
1747		!< Thickness values on points i-3, i-2, i-1, i, i+1, i+2 [L ~> m]
1748		real, optional, intent(in) :: u6(6)
1749		!< Velocity values on points i-3, i-2, i-1, i, i+1, i+2 [L T-1 ~> m s-1]
1750		real, intent(in) :: h_tiny !< A tiny thickness to prevent division by zero [L ~> m]
1751		real, intent(in) :: u !< Velocity or thickness flux on point i-1/2
1752		!! [L T-1 ~> m s-1] or [L2 T-1 ~> m2 s-1]
1753		logical, intent(in) :: velocity_smoothing !< If ture, use velocity to compute the smoothness indicator
1754		real, intent(inout) :: qr !< Reconstruction of tracer q on point i-1/2 [A ~> a]
1755		real :: vr ! Reconstruction of hq [A ~> a]
1756		real :: hr ! Reconstruction of h [L ~> m]
1757		real :: c0, c1, c2 ! Intermediate reconstruction of hq[A ~> a]
1758		real :: d0, d1, d2 ! Intermediate reconstruction of h [L ~> m]
1759		real :: b0, b1, b2 ! Smoothness indicator [A ~> a]
1760		real :: tau ! Difference of smoothness indicators [A ~> a]
1761		real :: w0, w1, w2 ! Weights [nondim]
1762		real :: s ! Temporary variables [nondim]
1763		real, parameter :: C3_10 = 3.0/10.0 ! The ratio of 3/10 [nondim]
1764		real, parameter :: C3_5 = 3.0/5.0 ! The ratio of 3/5 [nondim]
1765		real, parameter :: C1_10 = 1.0/10.0 ! The ratio of 1/10 [nondim]
1766
1767	0	if (u>0.) then
1768	0	call weno_five_reconstruction_0(q6(3:5), c0) ! Reconstruction in the third upwind stencil
1769	0	call weno_five_reconstruction_1(q6(2:4), c1) ! Reconstruction in the second upwind stencil
1770	0	call weno_five_reconstruction_2(q6(1:3), c2) ! Reconstruction in the first upwind stencil
1771
1772	0	call weno_five_reconstruction_0(h6(3:5), d0)
1773	0	call weno_five_reconstruction_1(h6(2:4), d1)
1774	0	call weno_five_reconstruction_2(h6(1:3), d2)
1775	0	if (velocity_smoothing) then
1776	0	call weno_five_weight_0(u6(3:5), b0) ! Smoothness indicator of the third upwind stencil
1777	0	call weno_five_weight_1(u6(2:4), b1) ! Smoothness indicator of the second upwind stencil
1778	0	call weno_five_weight_2(u6(1:3), b2) ! Smoothness indicator of the first upwind stencil
1779		else
1780	0	call weno_five_weight_0(q6(3:5), b0)
1781	0	call weno_five_weight_1(q6(2:4), b1)
1782	0	call weno_five_weight_2(q6(1:3), b2)
1783		endif
1784		else
1785	0	call weno_five_reconstruction_0(q6(4:2:-1), c0) ! Reconstruction in the third upwind stencil
1786	0	call weno_five_reconstruction_1(q6(5:3:-1), c1) ! Reconstruction in the second upwind stencil
1787	0	call weno_five_reconstruction_2(q6(6:4:-1), c2) ! Reconstruction in the first upwind stencil
1788
1789	0	call weno_five_reconstruction_0(h6(4:2:-1), d0)
1790	0	call weno_five_reconstruction_1(h6(5:3:-1), d1)
1791	0	call weno_five_reconstruction_2(h6(6:4:-1), d2)
1792	0	if (velocity_smoothing) then
1793	0	call weno_five_weight_0(u6(4:2:-1), b0) ! Smoothness indicator of the third upwind stencil
1794	0	call weno_five_weight_1(u6(5:3:-1), b1) ! Smoothness indicator of the second upwind stencil
1795	0	call weno_five_weight_2(u6(6:4:-1), b2) ! Smoothness indicator of the first upwind stencil
1796		else
1797	0	call weno_five_weight_0(q6(4:2:-1), b0)
1798	0	call weno_five_weight_1(q6(5:3:-1), b1)
1799	0	call weno_five_weight_2(q6(6:4:-1), b2)
1800		endif
1801		endif
1802
1803	0	tau = abs(b0 - b2)
1804	0	w0 = C3_10 * fac_fn(tau, b0)
1805	0	w1 = C3_5 * fac_fn(tau, b1)
1806	0	w2 = C1_10 * fac_fn(tau, b2)
1807
1808	0	s = 1. / ((w0 + w1) + w2)
1809	0	w0 = w0 * s ! Weights of stencils
1810	0	w1 = w1 * s
1811	0	w2 = w2 * s
1812
1813	0	vr = ((w0 * c0) + (w1 * c1)) + (w2 * c2)
1814	0	hr = ((w0 * d0) + (w1 * d1)) + (w2 * d2)
1815		! vr = min(max(q6(3), q6(4)), vr) ; vr = max(min(q6(3), q6(4)), vr) !Impose a monotonicity limiter
1816	0	hr = min(max(h6(3), h6(4)), hr) ; hr = max(min(h6(3), h6(4)), hr) !Impose a monotonicity limiter
1817
1818	0	qr = vr / max(hr, h_tiny)
1819
1820	0	end subroutine weno_five_h_weight_reconstruction
1821
1822		!> Compute the smoothness indicator for the third upwind stencil of the fifth-order WENO scheme
1823	0	subroutine weno_five_weight_0(q3, w0)
1824		real, intent(in) :: q3(3) !< Tracer values on the three-point stencil [A ~> a]
1825		real, intent(inout) :: w0 !< Smoothness indicator for this stencil [A2 ~> a2]
1826
1827		w0 = (q3(1) * ((10 * q3(1) - 31 * q3(2)) + 11 * q3(3))) + &
1828	0	((q3(2) * (25 * q3(2) - 19 * q3(3))) + 4 * (q3(3) * q3(3)))
1829
1830	0	end subroutine weno_five_weight_0
1831
1832		!> Compute the smoothness indicator for the second upwind stencil of the fifth-order WENO scheme
1833	0	subroutine weno_five_weight_1(q3, w1)
1834		real, intent(in) :: q3(3) !< Tracer values on the three-point stencil [A ~> a]
1835		real, intent(inout) :: w1 !< Smoothness indicator for this stencil [A2 ~> a2]
1836
1837		w1 = (q3(1) * ((4 * q3(1) - 13 * q3(2)) + 5 * q3(3))) + &
1838	0	((q3(2) * (13 * q3(2) - 13 * q3(3))) + 4 * (q3(3) * q3(3)))
1839
1840	0	end subroutine weno_five_weight_1
1841
1842		!> Compute the smoothness indicator for the first upwind stencil of the fifth-order WENO scheme
1843	0	subroutine weno_five_weight_2(q3, w2)
1844		real, intent(in) :: q3(3) !< Tracer values on the three-point stencil [A ~> a]
1845		real, intent(inout) :: w2 !< Smoothness indicator for this stencil [A2 ~> a2]
1846
1847		w2 = (q3(1) * ((4 * q3(1) - 19 * q3(2)) + 11 * q3(3))) + &
1848	0	((q3(2) * (25 * q3(2) - 31 * q3(3))) + 10 * (q3(3) * q3(3)))
1849
1850	0	end subroutine weno_five_weight_2
1851
1852		!> Reconstruction in the third upwind stencil of the fifth-order WENO scheme
1853	0	subroutine weno_five_reconstruction_0(q3, p0)
1854		real, intent(in) :: q3(3) !< Tracer values on three points [A ~> a]
1855		real, intent(inout) :: p0 !< Reconstruction of the quantity [A ~> a]
1856		real, parameter :: C1_6 = 1.0/6.0 ! One sixth [nondim]
1857
1858	0	p0 = ((2q3(1) + 5q3(2)) - q3(3)) * C1_6
1859
1860	0	end subroutine weno_five_reconstruction_0
1861
1862		!> Reconstruction in the second upwind stencil of the fifth-order WENO scheme
1863	0	subroutine weno_five_reconstruction_1(q3, p1)
1864		real, intent(in) :: q3(3) !< Tracer values on the three-point stencil [A ~> a]
1865		real, intent(inout) :: p1 !< Reconstruction of the quantity [A ~> a]
1866		real, parameter :: C1_6 = 1.0/6.0 ! One sixth [nondim]
1867
1868	0	p1 = ((-q3(1) + 5q3(2)) + 2q3(3)) * C1_6
1869
1870	0	end subroutine weno_five_reconstruction_1
1871
1872		!> Reconstruction in the first upwind stencil of the fifth-order WENO scheme
1873	0	subroutine weno_five_reconstruction_2(q3, p2)
1874		real, intent(in) :: q3(3) !< Tracer values on the three-point stencil [A ~> a]
1875		real, intent(inout) :: p2 !< Reconstruction of the quantity [A ~> a]
1876		real, parameter :: C1_6 = 1.0/6.0 ! One sixth [nondim]
1877
1878	0	p2 = ((2q3(1) - 7q3(2)) + 11q3(3)) C1_6
1879
1880	0	end subroutine weno_five_reconstruction_2
1881
1882
1883		!> Reconstruct the tracer (e.g., PV, vorticity) onto point i-1/2 using a seventh-order WENO scheme
1884		!! This reconstruction computes a thickness weighted average of PV
1885	0	subroutine weno_seven_h_weight_reconstruction(q8, h8, u8, &
1886		h_tiny, u, qr, velocity_smoothing)
1887		real, intent(in) :: q8(8)
1888		!< Tracer values on points i-4, i-3, i-2, i-1, i, i+1, i+2, i+3
1889		real, intent(in) :: h8(8)
1890		!< Thickness on the same tracer points i-4, i-3, i-2, i-1, i, i+1, i+2, i+3 [L ~> m]
1891		real, optional, intent(in) :: u8(8)
1892		!< Velocity values on points i-4, i-3, i-2, i-1, i, i+1, i+2, i+3 [L T-1 ~> m s-1]
1893		real, intent(in) :: h_tiny !< A tiny thickness to prevent division by zero [L ~> m]
1894		real, intent(in) :: u !< Velocity or thickness flux on point i-1/2
1895		!! [L T-1 ~> m s-1] or [L2 T-1 ~> m2 s-1]
1896		logical, intent(in) :: velocity_smoothing !< If true, use velocity to compute the smoothness indicator
1897		real, intent(inout) :: qr !< Reconstruction of tracer q on point i-1/2 [A ~> a]
1898		real :: vr ! Reconstruction of hq [A ~> a]
1899		real :: hr ! Reconstruction of h [L ~> m]
1900		real :: c0, c1, c2, c3 ! Intermediate reconstruction of hq [A ~> a]
1901		real :: d0, d1, d2, d3 ! Intermediate reconstruction of h [L ~> m]
1902		real :: b0, b1, b2, b3 ! Smoothness indicator [A ~> a]
1903		real :: tau ! Difference of smoothness indicators [A ~> a]
1904		real :: w0, w1, w2, w3 ! Weights [nondim]
1905		real :: s ! Temporary variables [nondim]
1906		real, parameter :: C4_35 = 4.0/35.0 ! The ratio of 4/35 [nondim]
1907		real, parameter :: C18_35 = 18.0/35.0 ! The ratio of 18/35 [nondim]
1908		real, parameter :: C12_35 = 12.0/35.0 ! The ratio of 12/35 [nondim]
1909		real, parameter :: C1_35 = 1.0/35.0 ! The ratio of 1/35 [nondim]
1910
1911	0	if (u>0.) then
1912	0	call weno_seven_reconstruction_0(q8(4:7), c0) ! Reconstruction in the fourth upwind stencil
1913	0	call weno_seven_reconstruction_1(q8(3:6), c1) ! Reconstruction in the third upwind stencil
1914	0	call weno_seven_reconstruction_2(q8(2:5), c2) ! Reconstruction in the second upwind stencil
1915	0	call weno_seven_reconstruction_3(q8(1:4), c3) ! Reconstruction in the first upwind stencil
1916
1917	0	call weno_seven_reconstruction_0(h8(4:7), d0) ! Reconstruction in the fourth upwind stencil
1918	0	call weno_seven_reconstruction_1(h8(3:6), d1) ! Reconstruction in the third upwind stencil
1919	0	call weno_seven_reconstruction_2(h8(2:5), d2) ! Reconstruction in the second upwind stencil
1920	0	call weno_seven_reconstruction_3(h8(1:4), d3) ! Reconstruction in the first upwind stencil
1921	0	if (velocity_smoothing) then
1922	0	call weno_seven_weight_0(u8(4:7), b0) ! Smoothness indicator of the fourth upwind stencil
1923	0	call weno_seven_weight_1(u8(3:6), b1) ! Smoothness indicator of the third upwind stencil
1924	0	call weno_seven_weight_2(u8(2:5), b2) ! Smoothness indicator of the second upwind stencil
1925	0	call weno_seven_weight_3(u8(1:4), b3) ! Smoothness indicator of the first upwind stencil
1926		else
1927	0	call weno_seven_weight_0(q8(4:7), b0)
1928	0	call weno_seven_weight_1(q8(3:6), b1)
1929	0	call weno_seven_weight_2(q8(2:5), b2)
1930	0	call weno_seven_weight_3(q8(1:4), b3)
1931		endif
1932		else
1933	0	call weno_seven_reconstruction_0(q8(5:2:-1), c0) ! Reconstruction in the fourth upwind stencil
1934	0	call weno_seven_reconstruction_1(q8(6:3:-1), c1) ! Reconstruction in the third upwind stencil
1935	0	call weno_seven_reconstruction_2(q8(7:4:-1), c2) ! Reconstruction in the second upwind stencil
1936	0	call weno_seven_reconstruction_3(q8(8:5:-1), c3) ! Reconstruction in the first upwind stencil
1937
1938	0	call weno_seven_reconstruction_0(h8(5:2:-1), d0)
1939	0	call weno_seven_reconstruction_1(h8(6:3:-1), d1)
1940	0	call weno_seven_reconstruction_2(h8(7:4:-1), d2)
1941	0	call weno_seven_reconstruction_3(h8(8:5:-1), d3)
1942	0	if (velocity_smoothing) then
1943	0	call weno_seven_weight_0(u8(5:2:-1), b0) ! Smoothness indicator of the fourth upwind stencil
1944	0	call weno_seven_weight_1(u8(6:3:-1), b1) ! Smoothness indicator of the third upwind stencil
1945	0	call weno_seven_weight_2(u8(7:4:-1), b2) ! Smoothness indicator of the second upwind stencil
1946	0	call weno_seven_weight_3(u8(8:5:-1), b3) ! Smoothness indicator of the first upwind stencil
1947		else
1948	0	call weno_seven_weight_0(q8(5:2:-1), b0)
1949	0	call weno_seven_weight_1(q8(6:3:-1), b1)
1950	0	call weno_seven_weight_2(q8(7:4:-1), b2)
1951	0	call weno_seven_weight_3(q8(8:5:-1), b3)
1952		endif
1953		endif
1954
1955	0	tau = abs((b0 - b3) + 3 * (b1 - b2))
1956	0	w0 = C4_35 * fac_fn(tau, b0)
1957	0	w1 = C18_35 * fac_fn(tau, b1)
1958	0	w2 = C12_35 * fac_fn(tau, b2)
1959	0	w3 = C1_35 * fac_fn(tau, b3)
1960
1961	0	s = 1. / ((w0 + w1) + (w2 + w3))
1962	0	w0 = w0 * s ! Weights of the stencils
1963	0	w1 = w1 * s
1964	0	w2 = w2 * s
1965	0	w3 = w3 * s
1966
1967	0	vr = ((w0 * c0) + (w1 * c1)) + ((w2 * c2) + (w3 * c3))
1968	0	hr = ((w0 * d0) + (w1 * d1)) + ((w2 * d2) + (w3 * d3))
1969
1970		! vr = min(max(q4, q5), vr) ; vr = max(min(q4, q5), vr)
1971	0	hr = min(max(h8(4), h8(5)), hr) ; hr = max(min(h8(4), h8(5)), hr) ! Impose a monotonicity limiter
1972
1973	0	qr = vr / max(hr, h_tiny)
1974
1975	0	end subroutine weno_seven_h_weight_reconstruction
1976
1977		!> Compute the smoothness indicator for the fourth upwind stencil of the seventh-order WENO scheme
1978	0	subroutine weno_seven_weight_0(q4, w0)
1979		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
1980		real, intent(inout) :: w0 !< Smoothness indicator for this stencil [A2 ~> a2]
1981
1982		! Coefficients from Balsara and Shu (2000). The division by 1000 will be normalized out by fac_fn
1983		w0 = ((q4(1) * ((2.107 * q4(1) - 9.402 * q4(2)) + (7.042 * q4(3) - 1.854 * q4(4)))) + &
1984		(q4(2) * ((11.003 * q4(2) - 17.246 * q4(3)) + 4.642 * q4(4)))) + &
1985	0	((q4(3) * (7.043 * q4(3) - 3.882 * q4(4))) + 0.547 * (q4(4) * q4(4)))
1986
1987	0	end subroutine weno_seven_weight_0
1988
1989		!> Compute the smoothness indicator for the third upwind stencil of the seventh-order WENO scheme
1990	0	subroutine weno_seven_weight_1(q4, w1)
1991		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
1992		real, intent(inout) :: w1 !< Smoothness indicator for this stencil [A2 ~> a2]
1993
1994		! Coefficients from Balsara and Shu (2000). The division by 1000 will be normalized out by fac_fn
1995		w1 = ((q4(1) * ((0.547 * q4(1) - 2.522 * q4(2)) + (1.922 * q4(3) - 0.494 * q4(4)))) + &
1996		(q4(2) * ((3.443 * q4(2) - 5.966 * q4(3)) + 1.602 * q4(4)))) + &
1997	0	((q4(3) * (2.843 * q4(3) - 1.642 * q4(4))) + 0.267 * (q4(4) * q4(4)))
1998
1999	0	end subroutine weno_seven_weight_1
2000
2001		!> Compute the smoothness indicator for the second upwind stencil of the seventh-order WENO scheme
2002	0	subroutine weno_seven_weight_2(q4, w2)
2003		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
2004		real, intent(inout) :: w2 !< Smoothness indicator for this stencil [A2 ~> a2]
2005
2006		! Coefficients from Balsara and Shu (2000). The division by 1000 will be normalized out by fac_fn
2007		w2 = ((q4(1) * ((0.267 * q4(1) - 1.642 * q4(2)) + (1.602 * q4(3) - 0.494 * q4(4)))) + &
2008		(q4(2) * ((2.843 * q4(2) - 5.966 * q4(3)) + 1.922 * q4(4)))) + &
2009	0	((q4(3) * (3.443 * q4(3) - 2.522 * q4(4))) + 0.547 * (q4(4) * q4(4)))
2010
2011	0	end subroutine weno_seven_weight_2
2012
2013		!> Compute smoothness indicator for the first upwind stencil of the seventh-order WENO scheme
2014	0	subroutine weno_seven_weight_3(q4, w3)
2015		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
2016		real, intent(inout) :: w3 !< Smoothness indicator for this stencil [A2 ~> a2]
2017
2018		! Coefficients from Balsara and Shu (2000). The division by 1000 will be normalized out by fac_fn
2019		w3 = ((q4(1) * ((0.547 * q4(1) - 3.882 * q4(2)) + (4.642 * q4(3) - 1.854 * q4(4)))) + &
2020		(q4(2) * ((7.043 * q4(2) - 17.246 * q4(3)) + 7.042 * q4(4)))) + &
2021	0	((q4(3) * (11.003 * q4(3) - 9.402 * q4(4))) + 2.107 * (q4(4) * q4(4)))
2022
2023	0	end subroutine weno_seven_weight_3
2024
2025		!> Reconstruction in the fourth upwind stencil for seventh-order WENO scheme
2026	0	subroutine weno_seven_reconstruction_0(q4, p0)
2027		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
2028		real, intent(inout) :: p0 !< Reconstruction of the quantity [A ~> a]
2029		real, parameter :: C1_24 = 1.0/24.0 ! One twenty fourth [nondim]
2030
2031	0	p0 = (((6 * q4(1) + 26 * q4(2)) - 10 * q4(3)) + 2 * q4(4)) * C1_24
2032
2033	0	end subroutine weno_seven_reconstruction_0
2034
2035		!> Reconstruction in the third upwind stencil for seventh-order WENO scheme
2036	0	subroutine weno_seven_reconstruction_1(q4, p1)
2037		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
2038		real, intent(inout) :: p1 !< Reconstruction of the quantity [A ~> a]
2039		real, parameter :: C1_24 = 1.0/24.0 ! One twenty fourth [nondim]
2040
2041	0	p1 = (14 * (q4(2) + q4(3)) - 2 * (q4(1) + q4(4))) * C1_24
2042
2043	0	end subroutine weno_seven_reconstruction_1
2044
2045		!> Reconstruction in the second upwind stencil for seventh-order WENO scheme
2046	0	subroutine weno_seven_reconstruction_2(q4, p2)
2047		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
2048		real, intent(inout) :: p2 !< Reconstruction of the quantity [A ~> a]
2049		real, parameter :: C1_24 = 1.0/24.0 ! One twenty fourth [nondim]
2050
2051	0	p2 = (((2 * q4(1) - 10 * q4(2)) + 26 * q4(3)) + 6 * q4(4)) * C1_24
2052
2053	0	end subroutine weno_seven_reconstruction_2
2054
2055		!> Reconstruction in the first upwind stencil for seventh-order WENO scheme
2056	0	subroutine weno_seven_reconstruction_3(q4, p3)
2057		real, intent(in) :: q4(4) !< Tracer values on the four-point stencil [A ~> a]
2058		real, intent(inout) :: p3 !< Reconstruction of the quantity [A ~> a]
2059		real, parameter :: C1_24 = 1.0/24.0 ! One twenty fourth [nondim]
2060
2061	0	p3 = (((-6 * q4(1) + 26 * q4(2)) - 46 * q4(3)) + 50 * q4(4)) * C1_24
2062
2063	0	end subroutine weno_seven_reconstruction_3
2064
2065	75	function CoriolisAdv_stencil(CS) result(stencil)
2066		type(CoriolisAdv_CS), intent(in) :: CS !< Control structure for MOM_CoriolisAdv
2067		integer :: stencil !< The halo stencil size for the Coriolis advection scheme
2068
2069	75	stencil = 2
2070	75	if (CS%Coriolis_Scheme == wenovi7th_PV_ENSTRO) stencil = 4
2071	75	if (CS%Coriolis_Scheme == wenovi5th_PV_ENSTRO) stencil = 3
2072
2073	75	end function CoriolisAdv_stencil
2074
2075
2076		!> Initializes the control structure for MOM_CoriolisAdv
2077	1	subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
2078		type(time_type), target, intent(in) :: Time !< Current model time
2079		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2080		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
2081		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2082		type(param_file_type), intent(in) :: param_file !< Runtime parameter handles
2083		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
2084		type(accel_diag_ptrs), target, intent(inout) :: AD !< Storage for acceleration diagnostics
2085		type(CoriolisAdv_CS), intent(inout) :: CS !< Control structure for MOM_CoriolisAdv
2086		! Local variables
2087		! This include declares and sets the variable "version".
2088		#include "version_variable.h"
2089		character(len=40) :: mdl = "MOM_CoriolisAdv" ! This module's name.
2090		character(len=20) :: tmpstr
2091		character(len=400) :: mesg
2092		logical :: use_weno
2093		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
2094		#ifdef __NVCOMPILER_OPENMP_GPU
2095		integer, parameter :: default_nkblock = 0
2096		#else
2097		integer, parameter :: default_nkblock = 1
2098		#endif
2099
2100	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
2101	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
2102
2103	1	CS%initialized = .true.
2104	1	CS%diag => diag ; CS%Time => Time
2105
2106		! Read all relevant parameters and write them to the model log.
2107	1	call log_version(param_file, mdl, version, "")
2108		call get_param(param_file, mdl, "CORIOLIS_ADV_NKBLOCK", CS%nkblock, &
2109		"The k-direction block size used in Coriolis and momentum advection "//&
2110		"calculations. The default 0 setting dynamically uses the full vertical column.", &
2111	1	default=default_nkblock, layoutParam=.true.)
2112	1	if (CS%nkblock < 0) call MOM_error(FATAL, "CORIOLIS_ADV_NKBLOCK must be >= 0.")
2113		call get_param(param_file, mdl, "NOSLIP", CS%no_slip, &
2114		"If true, no slip boundary conditions are used; otherwise "//&
2115		"free slip boundary conditions are assumed. The "//&
2116		"implementation of the free slip BCs on a C-grid is much "//&
2117		"cleaner than the no slip BCs. The use of free slip BCs "//&
2118		"is strongly encouraged, and no slip BCs are not used with "//&
2119	1	"the biharmonic viscosity.", default=.false.)
2120
2121		call get_param(param_file, mdl, "CORIOLIS_EN_DIS", CS%Coriolis_En_Dis, &
2122		"If true, two estimates of the thickness fluxes are used "//&
2123		"to estimate the Coriolis term, and the one that "//&
2124		"dissipates energy relative to the other one is used.", &
2125	1	default=.false.)
2126
2127		! Set %Coriolis_Scheme
2128		! (Select the baseline discretization for the Coriolis term)
2129		call get_param(param_file, mdl, "CORIOLIS_SCHEME", tmpstr, &
2130		"CORIOLIS_SCHEME selects the discretization for the "//&
2131		"Coriolis terms. Valid values are: \n"//&
2132		"\t SADOURNY75_ENERGY - Sadourny, 1975; energy cons. \n"//&
2133		"\t ARAKAWA_HSU90 - Arakawa & Hsu, 1990 \n"//&
2134		"\t SADOURNY75_ENSTRO - Sadourny, 1975; enstrophy cons. \n"//&
2135		"\t ARAKAWA_LAMB81 - Arakawa & Lamb, 1981; En. + Enst.\n"//&
2136		"\t ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with \n"//&
2137		"\t Arakawa & Hsu and Sadourny energy \n"//&
2138		"\t WENOVI5TH_PV_ENSTRO - 5th-order WENO PV enstrophy \n"//&
2139		"\t WENOVI3RD_PV_ENSTRO - 3rd-order WENO PV enstrophy \n"//&
2140		"\t WENOVI7TH_PV_ENSTRO - 7th-order WENO PV enstrophy \n", &
2141	1	default=SADOURNY75_ENERGY_STRING)
2142	1	tmpstr = uppercase(tmpstr)
2143	0	select case (tmpstr)
2144		case (SADOURNY75_ENERGY_STRING)
2145	0	CS%Coriolis_Scheme = SADOURNY75_ENERGY
2146		case (ARAKAWA_HSU_STRING)
2147	0	CS%Coriolis_Scheme = ARAKAWA_HSU90
2148		case (SADOURNY75_ENSTRO_STRING)
2149	1	CS%Coriolis_Scheme = SADOURNY75_ENSTRO
2150		case (ARAKAWA_LAMB_STRING)
2151	0	CS%Coriolis_Scheme = ARAKAWA_LAMB81
2152		case (AL_BLEND_STRING)
2153	0	CS%Coriolis_Scheme = AL_BLEND
2154		case (ROBUST_ENSTRO_STRING)
2155	0	CS%Coriolis_Scheme = ROBUST_ENSTRO
2156	0	CS%Coriolis_En_Dis = .false.
2157		case (WENOVI7TH_PV_ENSTRO_STRING)
2158	0	CS%Coriolis_Scheme = wenovi7th_PV_ENSTRO
2159		case (WENOVI5TH_PV_ENSTRO_STRING)
2160	0	CS%Coriolis_Scheme = wenovi5th_PV_ENSTRO
2161		case (WENOVI3RD_PV_ENSTRO_STRING)
2162	0	CS%Coriolis_Scheme = wenovi3rd_PV_ENSTRO
2163		case default
2164	0	call MOM_mesg('CoriolisAdv_init: Coriolis_Scheme ="'//trim(tmpstr)//'"', 0)
2165		call MOM_error(FATAL, "CoriolisAdv_init: Unrecognized setting "// &
2166	1	"#define CORIOLIS_SCHEME "//trim(tmpstr)//" found in input file.")
2167		end select
2168
2169		use_weno = CS%Coriolis_Scheme == wenovi7th_PV_ENSTRO &
2170		.or. CS%Coriolis_Scheme == wenovi5th_PV_ENSTRO &
2171	1	.or. CS%Coriolis_Scheme == wenovi3rd_PV_ENSTRO
2172
2173	1	if (use_weno) then
2174		call get_param(param_file, mdl, "WENO_VELOCITY_SMOOTH", CS%weno_velocity_smooth, &
2175		"If true, use velocity to compute weighting for WENO. ", &
2176	0	default=.false.)
2177		endif
2178
2179	1	if (CS%Coriolis_Scheme == AL_BLEND) then
2180		call get_param(param_file, mdl, "CORIOLIS_BLEND_WT_LIN", CS%wt_lin_blend, &
2181		"A weighting value for the ratio of inverse thicknesses, "//&
2182		"beyond which the blending between Sadourny Energy and "//&
2183		"Arakawa & Hsu goes linearly to 0 when CORIOLIS_SCHEME "//&
2184		"is ARAWAKA_LAMB_BLEND. This must be between 1 and 1e-16.", &
2185	0	units="nondim", default=0.125)
2186		call get_param(param_file, mdl, "CORIOLIS_BLEND_F_EFF_MAX", CS%F_eff_max_blend, &
2187		"The factor by which the maximum effective Coriolis "//&
2188		"acceleration from any point can be increased when "//&
2189		"blending different discretizations with the "//&
2190		"ARAKAWA_LAMB_BLEND Coriolis scheme. This must be "//&
2191		"greater than 2.0 (the max value for Sadourny energy).", &
2192	0	units="nondim", default=4.0)
2193	0	CS%wt_lin_blend = min(1.0, max(CS%wt_lin_blend,1e-16))
2194	0	if (CS%F_eff_max_blend < 2.0) call MOM_error(WARNING, "CoriolisAdv_init: "//&
2195	0	"CORIOLIS_BLEND_F_EFF_MAX should be at least 2.")
2196		endif
2197
2198		mesg = "If true, the Coriolis terms at u-points are bounded by "//&
2199		"the four estimates of (f+rv)v from the four neighboring "//&
2200	1	"v-points, and similarly at v-points."
2201	1	if (CS%Coriolis_En_Dis .and. (CS%Coriolis_Scheme == SADOURNY75_ENERGY)) then
2202		mesg = trim(mesg)//" This option is "//&
2203		"always effectively false with CORIOLIS_EN_DIS defined and "//&
2204	0	"CORIOLIS_SCHEME set to "//trim(SADOURNY75_ENERGY_STRING)//"."
2205		else
2206		mesg = trim(mesg)//" This option would "//&
2207		"have no effect on the SADOURNY Coriolis scheme if it "//&
2208	1	"were possible to use centered difference thickness fluxes."
2209		endif
2210		call get_param(param_file, mdl, "BOUND_CORIOLIS", CS%bound_Coriolis, mesg, &
2211	1	default=.false.)
2212	1	if ((CS%Coriolis_En_Dis .and. (CS%Coriolis_Scheme == SADOURNY75_ENERGY)) .or. &
2213	0	(CS%Coriolis_Scheme == ROBUST_ENSTRO)) CS%bound_Coriolis = .false.
2214
2215		! Set KE_Scheme (selects discretization of KE)
2216		call get_param(param_file, mdl, "KE_SCHEME", tmpstr, &
2217		"KE_SCHEME selects the discretization for acceleration "//&
2218		"due to the kinetic energy gradient. Valid values are: \n"//&
2219		"\t KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV, KE_UP3", &
2220	1	default=KE_ARAKAWA_STRING)
2221	1	tmpstr = uppercase(tmpstr)
2222	1	select case (tmpstr)
2223	1	case (KE_ARAKAWA_STRING); CS%KE_Scheme = KE_ARAKAWA
2224	0	case (KE_SIMPLE_GUDONOV_STRING); CS%KE_Scheme = KE_SIMPLE_GUDONOV
2225	0	case (KE_GUDONOV_STRING); CS%KE_Scheme = KE_GUDONOV
2226	0	case (KE_UP3_STRING); CS%KE_Scheme = KE_UP3
2227		case default
2228	0	call MOM_mesg('CoriolisAdv_init: KE_Scheme ="'//trim(tmpstr)//'"', 0)
2229		call MOM_error(FATAL, "CoriolisAdv_init: "// &
2230	1	"#define KE_SCHEME "//trim(tmpstr)//" in input file is invalid.")
2231		end select
2232
2233	1	if (CS%KE_Scheme == KE_UP3) then
2234		call get_param(param_file, mdl, "KE_USE_LIMITER", CS%KE_use_limiter, &
2235		"If true, use Koren limiter for KE_UP3 scheme", &
2236	0	default=.True.)
2237		endif
2238
2239		! Set PV_Adv_Scheme (selects discretization of PV advection)
2240		call get_param(param_file, mdl, "PV_ADV_SCHEME", tmpstr, &
2241		"PV_ADV_SCHEME selects the discretization for PV "//&
2242		"advection. Valid values are: \n"//&
2243		"\t PV_ADV_CENTERED - centered (aka Sadourny, 75) \n"//&
2244		"\t PV_ADV_UPWIND1 - upwind, first order", &
2245	1	default=PV_ADV_CENTERED_STRING)
2246	2	select case (uppercase(tmpstr))
2247		case (PV_ADV_CENTERED_STRING)
2248	1	CS%PV_Adv_Scheme = PV_ADV_CENTERED
2249		case (PV_ADV_UPWIND1_STRING)
2250	0	CS%PV_Adv_Scheme = PV_ADV_UPWIND1
2251		case default
2252	0	call MOM_mesg('CoriolisAdv_init: PV_Adv_Scheme ="'//trim(tmpstr)//'"', 0)
2253		call MOM_error(FATAL, "CoriolisAdv_init: "// &
2254	1	"#DEFINE PV_ADV_SCHEME in input file is invalid.")
2255		end select
2256
2257		CS%id_rv = register_diag_field('ocean_model', 'RV', diag%axesBL, Time, &
2258	1	'Relative Vorticity', 's-1', conversion=US%s_to_T)
2259
2260		CS%id_PV = register_diag_field('ocean_model', 'PV', diag%axesBL, Time, &
2261	1	'Potential Vorticity', 'm-1 s-1', conversion=GV%m_to_H*US%s_to_T)
2262
2263		CS%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, Time, &
2264	1	'Zonal Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2)
2265
2266		CS%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, Time, &
2267	1	'Meridional Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2)
2268
2269		CS%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, Time, &
2270	1	'Meridional Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
2271
2272		CS%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, Time, &
2273	1	'Zonal Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
2274
2275		CS%id_CAuS = register_diag_field('ocean_model', 'CAu_Stokes', diag%axesCuL, Time, &
2276	1	'Zonal Acceleration from Stokes Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
2277		! add to AD
2278
2279		CS%id_CAvS = register_diag_field('ocean_model', 'CAv_Stokes', diag%axesCvL, Time, &
2280	1	'Meridional Acceleration from Stokes Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
2281		! add to AD
2282
2283		!CS%id_hf_gKEu = register_diag_field('ocean_model', 'hf_gKEu', diag%axesCuL, Time, &
2284		! 'Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
2285		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
2286		CS%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCu1, Time, &
2287		'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
2288	1	'm s-2', conversion=US%L_T2_to_m_s2)
2289
2290		!CS%id_hf_gKEv = register_diag_field('ocean_model', 'hf_gKEv', diag%axesCvL, Time, &
2291		! 'Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
2292		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
2293		CS%id_hf_gKEv_2d = register_diag_field('ocean_model', 'hf_gKEv_2d', diag%axesCv1, Time, &
2294		'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
2295	1	'm s-2', conversion=US%L_T2_to_m_s2)
2296
2297		CS%id_h_gKEu = register_diag_field('ocean_model', 'h_gKEu', diag%axesCuL, Time, &
2298		'Thickness Multiplied Zonal Acceleration from Grad. Kinetic Energy', &
2299	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2300		CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
2301		'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', &
2302	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2303
2304		CS%id_h_gKEv = register_diag_field('ocean_model', 'h_gKEv', diag%axesCvL, Time, &
2305		'Thickness Multiplied Meridional Acceleration from Grad. Kinetic Energy', &
2306	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2307		CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, &
2308		'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', &
2309	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2310
2311		!CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, &
2312		! 'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
2313		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
2314		CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCv1, Time, &
2315		'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
2316	1	'm s-2', conversion=US%L_T2_to_m_s2)
2317
2318		!CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, &
2319		! 'Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
2320		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
2321		CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCu1, Time, &
2322		'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
2323	1	'm s-2', conversion=US%L_T2_to_m_s2)
2324
2325		CS%id_h_rvxu = register_diag_field('ocean_model', 'h_rvxu', diag%axesCvL, Time, &
2326		'Thickness Multiplied Meridional Acceleration from Relative Vorticity', &
2327	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2328		CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
2329		'Depth-integral of Meridional Acceleration from Relative Vorticity', &
2330	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2331
2332		CS%id_h_rvxv = register_diag_field('ocean_model', 'h_rvxv', diag%axesCuL, Time, &
2333		'Thickness Multiplied Zonal Acceleration from Relative Vorticity', &
2334	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2335		CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
2336		'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
2337	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2338
2339		! Allocate memory needed for the diagnostics that have been enabled.
2340		if ((CS%id_gKEu > 0) .or. (CS%id_hf_gKEu_2d > 0) .or. &
2341		! (CS%id_hf_gKEu > 0) .or. &
2342	1	(CS%id_h_gKEu > 0) .or. (CS%id_intz_gKEu_2d > 0)) then
2343	0	call safe_alloc_ptr(AD%gradKEu, IsdB, IedB, jsd, jed, nz)
2344		endif
2345		if ((CS%id_gKEv > 0) .or. (CS%id_hf_gKEv_2d > 0) .or. &
2346		! (CS%id_hf_gKEv > 0) .or. &
2347	1	(CS%id_h_gKEv > 0) .or. (CS%id_intz_gKEv_2d > 0)) then
2348	0	call safe_alloc_ptr(AD%gradKEv, isd, ied, JsdB, JedB, nz)
2349		endif
2350		if ((CS%id_rvxu > 0) .or. (CS%id_hf_rvxu_2d > 0) .or. &
2351		! (CS%id_hf_rvxu > 0) .or. &
2352	1	(CS%id_h_rvxu > 0) .or. (CS%id_intz_rvxu_2d > 0)) then
2353	0	call safe_alloc_ptr(AD%rv_x_u, isd, ied, JsdB, JedB, nz)
2354		endif
2355		if ((CS%id_rvxv > 0) .or. (CS%id_hf_rvxv_2d > 0) .or. &
2356		! (CS%id_hf_rvxv > 0) .or. &
2357	1	(CS%id_h_rvxv > 0) .or. (CS%id_intz_rvxv_2d > 0)) then
2358	0	call safe_alloc_ptr(AD%rv_x_v, IsdB, IedB, jsd, jed, nz)
2359		endif
2360
2361	1	if ((CS%id_hf_gKEv_2d > 0) .or. &
2362		! (CS%id_hf_gKEv > 0) .or. (CS%id_hf_rvxu > 0) .or. &
2363		(CS%id_hf_rvxu_2d > 0)) then
2364	0	call safe_alloc_ptr(AD%diag_hfrac_v, isd, ied, JsdB, JedB, nz)
2365		endif
2366	1	if ((CS%id_hf_gKEu_2d > 0) .or. &
2367		! (CS%id_hf_gKEu > 0) .or. (CS%id_hf_rvxv > 0) .or. &
2368		(CS%id_hf_rvxv_2d > 0)) then
2369	0	call safe_alloc_ptr(AD%diag_hfrac_u, IsdB, IedB, jsd, jed, nz)
2370		endif
2371		if ((CS%id_h_gKEu > 0) .or. (CS%id_intz_gKEu_2d > 0) .or. &
2372	1	(CS%id_h_rvxv > 0) .or. (CS%id_intz_rvxv_2d > 0)) then
2373	0	call safe_alloc_ptr(AD%diag_hu, IsdB, IedB, jsd, jed, nz)
2374		endif
2375		if ((CS%id_h_gKEv > 0) .or. (CS%id_intz_gKEv_2d > 0) .or. &
2376	1	(CS%id_h_rvxu > 0) .or. (CS%id_intz_rvxu_2d > 0)) then
2377	0	call safe_alloc_ptr(AD%diag_hv, isd, ied, JsdB, JedB, nz)
2378		endif
2379
2380	1	end subroutine CoriolisAdv_init
2381
2382		!> Destructor for coriolisadv_cs
2383	1	subroutine CoriolisAdv_end(CS)
2384		type(CoriolisAdv_CS), intent(inout) :: CS !< Control structure for MOM_CoriolisAdv
2385	1	end subroutine CoriolisAdv_end
2386
2387		!> \namespace mom_coriolisadv
2388		!!
2389		!! This file contains the subroutine that calculates the time
2390		!! derivatives of the velocities due to Coriolis acceleration and
2391		!! momentum advection. This subroutine uses either a vorticity
2392		!! advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or
2393		!! Sadourny's (JAS 1975) energy conserving scheme. Both have been
2394		!! modified to use general orthogonal coordinates as described in
2395		!! Arakawa and Lamb, Mon. Wea. Rev. 1981. Both schemes are second
2396		!! order accurate, and allow for vanishingly small layer thicknesses.
2397		!! The Arakawa and Hsu scheme globally conserves both total energy
2398		!! and potential enstrophy in the limit of nondivergent flow.
2399		!! Sadourny's energy conserving scheme conserves energy if the flow
2400		!! is nondivergent or centered difference thickness fluxes are used.
2401		!!
2402		!! A small fragment of the grid is shown below:
2403		!! \verbatim
2404		!!
2405		!! j+1 x ^ x ^ x At x: q, CoriolisBu
2406		!! j+1 > o > o > At ^: v, CAv, vh
2407		!! j x ^ x ^ x At >: u, CAu, uh, a, b, c, d
2408		!! j > o > o > At o: h, KE
2409		!! j-1 x ^ x ^ x
2410		!! i-1 i i+1 At x & ^:
2411		!! i i+1 At > & o:
2412		!! \endverbatim
2413		!!
2414		!! The boundaries always run through q grid points (x).
2415
2416	0	end module MOM_CoriolisAdv