MOM_CoriolisAdv.F90

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