MOM6
MOM_CoriolisAdv.F90
1 !> Accelerations due to the Coriolis force and momentum advection
2 module mom_coriolisadv
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 !> \author Robert Hallberg, April 1994 - June 2002
7 
8 use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
9 use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
10 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
11 use mom_file_parser, only : get_param, log_version, param_file_type
12 use mom_grid, only : ocean_grid_type
13 use mom_open_boundary, only : ocean_obc_type, obc_direction_e, obc_direction_w
14 use mom_open_boundary, only : obc_direction_n, obc_direction_s
15 use mom_string_functions, only : uppercase
16 use mom_unit_scaling, only : unit_scale_type
17 use mom_variables, only : accel_diag_ptrs
18 use mom_verticalgrid, only : verticalgrid_type
19 
20 implicit none ; private
21 
22 public coradcalc, coriolisadv_init, coriolisadv_end
23 
24 #include <MOM_memory.h>
25 
26 !> Control structure for mom_coriolisadv
27 type, public :: coriolisadv_cs ; private
28  integer :: coriolis_scheme !< Selects the discretization for the Coriolis terms.
29  !! Valid values are:
30  !! - SADOURNY75_ENERGY - Sadourny, 1975
31  !! - ARAKAWA_HSU90 - Arakawa & Hsu, 1990, Energy & non-div. Enstrophy
32  !! - ROBUST_ENSTRO - Pseudo-enstrophy scheme
33  !! - SADOURNY75_ENSTRO - Sadourny, JAS 1975, Enstrophy
34  !! - ARAKAWA_LAMB81 - Arakawa & Lamb, MWR 1981, Energy & Enstrophy
35  !! - ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with Arakawa & Hsu and Sadourny energy.
36  !! The default, SADOURNY75_ENERGY, is the safest choice then the
37  !! deformation radius is poorly resolved.
38  integer :: ke_scheme !< KE_SCHEME selects the discretization for
39  !! the kinetic energy. Valid values are:
40  !! KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV
41  integer :: pv_adv_scheme !< PV_ADV_SCHEME selects the discretization for PV advection
42  !! Valid values are:
43  !! - PV_ADV_CENTERED - centered (aka Sadourny, 75)
44  !! - PV_ADV_UPWIND1 - upwind, first order
45  real :: f_eff_max_blend !< The factor by which the maximum effective Coriolis
46  !! acceleration from any point can be increased when
47  !! blending different discretizations with the
48  !! ARAKAWA_LAMB_BLEND Coriolis scheme. This must be
49  !! greater than 2.0, and is 4.0 by default.
50  real :: wt_lin_blend !< A weighting value beyond which the blending between
51  !! Sadourny and Arakawa & Hsu goes linearly to 0.
52  !! This must be between 1 and 1e-15, often 1/8.
53  logical :: no_slip !< If true, no slip boundary conditions are used.
54  !! Otherwise free slip boundary conditions are assumed.
55  !! The implementation of the free slip boundary
56  !! conditions on a C-grid is much cleaner than the
57  !! no slip boundary conditions. The use of free slip
58  !! b.c.s is strongly encouraged. The no slip b.c.s
59  !! are not implemented with the biharmonic viscosity.
60  logical :: bound_coriolis !< If true, the Coriolis terms at u points are
61  !! bounded by the four estimates of (f+rv)v from the
62  !! four neighboring v points, and similarly at v
63  !! points. This option would have no effect on the
64  !! SADOURNY75_ENERGY scheme if it were possible to
65  !! use centered difference thickness fluxes.
66  logical :: coriolis_en_dis !< If CORIOLIS_EN_DIS is defined, two estimates of
67  !! the thickness fluxes are used to estimate the
68  !! Coriolis term, and the one that dissipates energy
69  !! relative to the other one is used. This is only
70  !! available at present if Coriolis scheme is
71  !! SADOURNY75_ENERGY.
72  type(time_type), pointer :: time !< A pointer to the ocean model's clock.
73  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output.
74  !>@{ Diagnostic IDs
75  integer :: id_rv = -1, id_pv = -1, id_gkeu = -1, id_gkev = -1
76  integer :: id_rvxu = -1, id_rvxv = -1 !!@}
77 end type coriolisadv_cs
78 
79 !>@{ Enumeration values for Coriolis_Scheme
80 integer, parameter :: sadourny75_energy = 1
81 integer, parameter :: arakawa_hsu90 = 2
82 integer, parameter :: robust_enstro = 3
83 integer, parameter :: sadourny75_enstro = 4
84 integer, parameter :: arakawa_lamb81 = 5
85 integer, parameter :: al_blend = 6
86 character*(20), parameter :: sadourny75_energy_string = "SADOURNY75_ENERGY"
87 character*(20), parameter :: arakawa_hsu_string = "ARAKAWA_HSU90"
88 character*(20), parameter :: robust_enstro_string = "ROBUST_ENSTRO"
89 character*(20), parameter :: sadourny75_enstro_string = "SADOURNY75_ENSTRO"
90 character*(20), parameter :: arakawa_lamb_string = "ARAKAWA_LAMB81"
91 character*(20), parameter :: al_blend_string = "ARAKAWA_LAMB_BLEND"
92 !!@}
93 !>@{ Enumeration values for KE_Scheme
94 integer, parameter :: ke_arakawa = 10
95 integer, parameter :: ke_simple_gudonov = 11
96 integer, parameter :: ke_gudonov = 12
97 character*(20), parameter :: ke_arakawa_string = "KE_ARAKAWA"
98 character*(20), parameter :: ke_simple_gudonov_string = "KE_SIMPLE_GUDONOV"
99 character*(20), parameter :: ke_gudonov_string = "KE_GUDONOV"
100 !!@}
101 !>@{ Enumeration values for PV_Adv_Scheme
102 integer, parameter :: pv_adv_centered = 21
103 integer, parameter :: pv_adv_upwind1 = 22
104 character*(20), parameter :: pv_adv_centered_string = "PV_ADV_CENTERED"
105 character*(20), parameter :: pv_adv_upwind1_string = "PV_ADV_UPWIND1"
106 !!@}
107 
108 contains
109 
110 !> Calculates the Coriolis and momentum advection contributions to the acceleration.
111 subroutine coradcalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
112  type(ocean_grid_type), intent(in) :: g !< Ocen grid structure
113  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
114  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [m s-1]
115  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [m s-1]
116  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
117  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uh !< Zonal transport u*h*dy
118  !! [H m2 s-1 ~> m3 s-1 or kg s-1]
119  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vh !< Meridional transport v*h*dx
120  !! [H m2 s-1 ~> m3 s-1 or kg s-1]
121  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: cau !< Zonal acceleration due to Coriolis
122  !! and momentum advection [m s-2].
123  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: cav !< Meridional acceleration due to Coriolis
124  !! and momentum advection [m s-2].
125  type(ocean_obc_type), pointer :: obc !< Open boundary control structure
126  type(accel_diag_ptrs), intent(inout) :: ad !< Storage for acceleration diagnostics
127  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
128  type(coriolisadv_cs), pointer :: cs !< Control structure for MOM_CoriolisAdv
129 
130  ! Local variables
131  real, dimension(SZIB_(G),SZJB_(G)) :: &
132  q, & ! Layer potential vorticity [m-1 s-1].
133  ih_q, & ! The inverse of thickness interpolated to q points [H-1 ~> m-1 or m2 kg-1].
134  area_q ! The sum of the ocean areas at the 4 adjacent thickness points [m2].
135 
136  real, dimension(SZIB_(G),SZJ_(G)) :: &
137  a, b, c, d ! a, b, c, & d are combinations of the potential vorticities
138  ! surrounding an h grid point. At small scales, a = q/4,
139  ! b = q/4, etc. All are in [H-1 s-1 ~> m-1 s-1 or m2 kg-1 s-1],
140  ! and use the indexing of the corresponding u point.
141 
142  real, dimension(SZI_(G),SZJ_(G)) :: &
143  area_h, & ! The ocean area at h points [m2]. Area_h is used to find the
144  ! average thickness in the denominator of q. 0 for land points.
145  ke ! Kinetic energy per unit mass [m2 s-2], KE = (u^2 + v^2)/2.
146  real, dimension(SZIB_(G),SZJ_(G)) :: &
147  harea_u, & ! The cell area weighted thickness interpolated to u points
148  ! times the effective areas [H m2 ~> m3 or kg].
149  kex, & ! The zonal gradient of Kinetic energy per unit mass [m s-2],
150  ! KEx = d/dx KE.
151  uh_center ! Transport based on arithmetic mean h at u-points [H m2 s-1 ~> m3 s-1 or kg s-1]
152  real, dimension(SZI_(G),SZJB_(G)) :: &
153  harea_v, & ! The cell area weighted thickness interpolated to v points
154  ! times the effective areas [H m2 ~> m3 or kg].
155  key, & ! The meridonal gradient of Kinetic energy per unit mass [m s-2],
156  ! KEy = d/dy KE.
157  vh_center ! Transport based on arithmetic mean h at v-points [H m2 s-1 ~> m3 s-1 or kg s-1]
158  real, dimension(SZI_(G),SZJ_(G)) :: &
159  uh_min, uh_max, & ! The smallest and largest estimates of the volume
160  vh_min, vh_max, & ! fluxes through the faces (i.e. u*h*dy & v*h*dx)
161  ! [H m2 s-1 ~> m3 s-1 or kg s-1].
162  ep_u, ep_v ! Additional pseudo-Coriolis terms in the Arakawa and Lamb
163  ! discretization [H-1 s-1 ~> m-1 s-1 or m2 kg-1 s-1].
164  real, dimension(SZIB_(G),SZJB_(G)) :: &
165  dvdx,dudy, &! Contributions to the circulation around q-points [m2 s-1]
166  abs_vort, & ! Absolute vorticity at q-points [s-1].
167  q2, & ! Relative vorticity over thickness [H-1 s-1 ~> m-1 s-1 or m2 kg-1 s-1].
168  max_fvq, & ! The maximum or minimum of the
169  min_fvq, & ! adjacent values of (-u) or v times
170  max_fuq, & ! the absolute vorticity [m s-2].
171  min_fuq ! All are defined at q points.
172  real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
173  pv, & ! A diagnostic array of the potential vorticities [m-1 s-1].
174  rv ! A diagnostic array of the relative vorticities [s-1].
175  real :: fv1, fv2, fu1, fu2 ! (f+rv)*v or (f+rv)*u [m s-2].
176  real :: max_fv, max_fu ! The maximum or minimum of the neighboring Coriolis
177  real :: min_fv, min_fu ! accelerations [m s-2], i.e. max(min)_fu(v)q.
178 
179  real, parameter :: c1_12=1.0/12.0 ! C1_12 = 1/12
180  real, parameter :: c1_24=1.0/24.0 ! C1_24 = 1/24
181  real :: absolute_vorticity ! Absolute vorticity [s-1].
182  real :: relative_vorticity ! Relative vorticity [s-1].
183  real :: ih ! Inverse of thickness [H-1 ~> m-1 or m2 kg-1].
184  real :: max_ihq, min_ihq ! The maximum and minimum of the nearby Ihq [H-1 ~> m-1 or m2 kg-1].
185  real :: harea_q ! The sum of area times thickness of the cells
186  ! surrounding a q point [H m2 ~> m3 or kg].
187  real :: h_neglect ! A thickness that is so small it is usually
188  ! lost in roundoff and can be neglected [H ~> m or kg m-2].
189  real :: temp1, temp2 ! Temporary variables [m2 s-2].
190  real, parameter :: eps_vel=1.0e-10 ! A tiny, positive velocity [m s-1].
191 
192  real :: uhc, vhc ! Centered estimates of uh and vh [H m2 s-1 ~> m3 s-1 or kg s-1].
193  real :: uhm, vhm ! The input estimates of uh and vh [H m2 s-1 ~> m3 s-1 or kg s-1].
194  real :: c1, c2, c3, slope ! Nondimensional parameters for the Coriolis limiter scheme.
195 
196  real :: fe_m2 ! Nondimensional temporary variables asssociated with
197  real :: rat_lin ! the ARAKAWA_LAMB_BLEND scheme.
198  real :: rat_m1 ! The ratio of the maximum neighboring inverse thickness
199  ! to the minimum inverse thickness minus 1. rat_m1 >= 0.
200  real :: al_wt ! The relative weight of the Arakawa & Lamb scheme to the
201  ! Arakawa & Hsu scheme, nondimensional between 0 and 1.
202  real :: sad_wt ! The relative weight of the Sadourny energy scheme to
203  ! the other two with the ARAKAWA_LAMB_BLEND scheme,
204  ! nondimensional between 0 and 1.
205 
206  real :: heff1, heff2 ! Temporary effective H at U or V points [H ~> m or kg m-2].
207  real :: heff3, heff4 ! Temporary effective H at U or V points [H ~> m or kg m-2].
208  real :: h_tiny ! A very small thickness [H ~> m or kg m-2].
209  real :: uheff, vheff ! More temporary variables [H m2 s-1 ~> m3 s-1 or kg s-1].
210  real :: quheff,qvheff ! More temporary variables [H m2 s-1 ~> m3 s-1 or kg s-1].
211  integer :: i, j, k, n, is, ie, js, je, isq, ieq, jsq, jeq, nz
212 
213 ! To work, the following fields must be set outside of the usual
214 ! is to ie range before this subroutine is called:
215 ! 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),
216 ! uh(is-1,ie,js:je+1) and vh(is:ie+1,js-1:je).
217 
218  if (.not.associated(cs)) call mom_error(fatal, &
219  "MOM_CoriolisAdv: Module must be initialized before it is used.")
220  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
221  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
222  h_neglect = gv%H_subroundoff
223  h_tiny = gv%Angstrom_H ! Perhaps this should be set to h_neglect instead.
224 
225  !$OMP parallel do default(private) shared(Isq,Ieq,Jsq,Jeq,G,Area_h)
226  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
227  area_h(i,j) = g%mask2dT(i,j) * g%areaT(i,j)
228  enddo ; enddo
229  if (associated(obc)) then ; do n=1,obc%number_of_segments
230  if (.not. obc%segment(n)%on_pe) cycle
231  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
232  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
233  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
234  if (obc%segment(n)%direction == obc_direction_n) then
235  area_h(i,j+1) = area_h(i,j)
236  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
237  area_h(i,j) = area_h(i,j+1)
238  endif
239  enddo
240  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
241  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
242  if (obc%segment(n)%direction == obc_direction_e) then
243  area_h(i+1,j) = area_h(i,j)
244  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
245  area_h(i,j) = area_h(i+1,j)
246  endif
247  enddo
248  endif
249  enddo ; endif
250  !$OMP parallel do default(private) shared(Isq,Ieq,Jsq,Jeq,G,Area_h,Area_q)
251  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
252  area_q(i,j) = (area_h(i,j) + area_h(i+1,j+1)) + &
253  (area_h(i+1,j) + area_h(i,j+1))
254  enddo ; enddo
255 
256  !$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,&
257  !$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,h_neglect,h_tiny,OBC)
258  do k=1,nz
259  ! Here the second order accurate layer potential vorticities, q,
260  ! are calculated. hq is second order accurate in space. Relative
261  ! vorticity is second order accurate everywhere with free slip b.c.s,
262  ! but only first order accurate at boundaries with no slip b.c.s.
263  ! First calculate the contributions to the circulation around the q-point.
264  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
265  dvdx(i,j) = v(i+1,j,k)*g%dyCv(i+1,j) - v(i,j,k)*g%dyCv(i,j)
266  dudy(i,j) = u(i,j+1,k)*g%dxCu(i,j+1) - u(i,j,k)*g%dxCu(i,j)
267  enddo ; enddo
268  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+2
269  harea_v(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i,j+1) * h(i,j+1,k))
270  enddo ; enddo
271  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+1
272  harea_u(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i+1,j) * h(i+1,j,k))
273  enddo ; enddo
274  if (cs%Coriolis_En_Dis) then
275  do j=jsq,jeq+1 ; do i=is-1,ie
276  uh_center(i,j) = 0.5 * (g%dy_Cu(i,j) * u(i,j,k)) * (h(i,j,k) + h(i+1,j,k))
277  enddo ; enddo
278  do j=js-1,je ; do i=isq,ieq+1
279  vh_center(i,j) = 0.5 * (g%dx_Cv(i,j) * v(i,j,k)) * (h(i,j,k) + h(i,j+1,k))
280  enddo ; enddo
281  endif
282 
283  ! Adjust circulation components to relative vorticity and thickness projected onto
284  ! velocity points on open boundaries.
285  if (associated(obc)) then ; do n=1,obc%number_of_segments
286  if (.not. obc%segment(n)%on_pe) cycle
287  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
288  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
289  if (obc%zero_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
290  dvdx(i,j) = 0. ; dudy(i,j) = 0.
291  enddo ; endif
292  if (obc%freeslip_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
293  dudy(i,j) = 0.
294  enddo ; endif
295  if (obc%computed_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
296  if (obc%segment(n)%direction == obc_direction_n) then
297  dudy(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%dxCu(i,j)
298  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
299  dudy(i,j) = 2.0*(u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dxCu(i,j+1)
300  endif
301  enddo ; endif
302  if (obc%specified_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
303  if (obc%segment(n)%direction == obc_direction_n) then
304  dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j)*g%dyBu(i,j)
305  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
306  dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j+1)*g%dyBu(i,j)
307  endif
308  enddo ; endif
309 
310  ! Project thicknesses across OBC points with a no-gradient condition.
311  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
312  if (obc%segment(n)%direction == obc_direction_n) then
313  harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j,k)
314  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
315  harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
316  endif
317  enddo
318 
319  if (cs%Coriolis_En_Dis) then
320  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
321  if (obc%segment(n)%direction == obc_direction_n) then
322  vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j,k)
323  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
324  vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j+1,k)
325  endif
326  enddo
327  endif
328  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
329  if (obc%zero_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
330  dvdx(i,j) = 0. ; dudy(i,j) = 0.
331  enddo ; endif
332  if (obc%freeslip_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
333  dvdx(i,j) = 0.
334  enddo ; endif
335  if (obc%computed_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
336  if (obc%segment(n)%direction == obc_direction_e) then
337  dvdx(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%dyCv(i,j)
338  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
339  dvdx(i,j) = 2.0*(v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dyCv(i+1,j)
340  endif
341  enddo ; endif
342  if (obc%specified_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
343  if (obc%segment(n)%direction == obc_direction_e) then
344  dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i,j)*g%dxBu(i,j)
345  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
346  dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i+1,j)*g%dxBu(i,j)
347  endif
348  enddo ; endif
349 
350  ! Project thicknesses across OBC points with a no-gradient condition.
351  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
352  if (obc%segment(n)%direction == obc_direction_e) then
353  harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i,j,k)
354  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
355  harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i+1,j,k)
356  endif
357  enddo
358  if (cs%Coriolis_En_Dis) then
359  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
360  if (obc%segment(n)%direction == obc_direction_e) then
361  uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i,j,k)
362  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
363  uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i+1,j,k)
364  endif
365  enddo
366  endif
367  endif
368  enddo ; endif
369 
370  if (associated(obc)) then ; do n=1,obc%number_of_segments
371  if (.not. obc%segment(n)%on_pe) cycle
372  ! Now project thicknesses across cell-corner points in the OBCs. The two
373  ! projections have to occur in sequence and can not be combined easily.
374  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
375  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
376  do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
377  if (obc%segment(n)%direction == obc_direction_n) then
378  if (area_h(i,j) + area_h(i+1,j) > 0.0) then
379  harea_u(i,j+1) = harea_u(i,j) * ((area_h(i,j+1) + area_h(i+1,j+1)) / &
380  (area_h(i,j) + area_h(i+1,j)))
381  else ; harea_u(i,j+1) = 0.0 ; endif
382  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
383  if (area_h(i,j+1) + area_h(i+1,j+1) > 0.0) then
384  harea_u(i,j) = harea_u(i,j+1) * ((area_h(i,j) + area_h(i+1,j)) / &
385  (area_h(i,j+1) + area_h(i+1,j+1)))
386  else ; harea_u(i,j) = 0.0 ; endif
387  endif
388  enddo
389  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
390  do j = max(jsq-1,obc%segment(n)%HI%JsdB), min(jeq+1,obc%segment(n)%HI%JedB)
391  if (obc%segment(n)%direction == obc_direction_e) then
392  if (area_h(i,j) + area_h(i,j+1) > 0.0) then
393  harea_v(i+1,j) = harea_v(i,j) * ((area_h(i+1,j) + area_h(i+1,j+1)) / &
394  (area_h(i,j) + area_h(i,j+1)))
395  else ; harea_v(i+1,j) = 0.0 ; endif
396  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
397  harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
398  if (area_h(i+1,j) + area_h(i+1,j+1) > 0.0) then
399  harea_v(i,j) = harea_v(i+1,j) * ((area_h(i,j) + area_h(i,j+1)) / &
400  (area_h(i+1,j) + area_h(i+1,j+1)))
401  else ; harea_v(i,j) = 0.0 ; endif
402  endif
403  enddo
404  endif
405  enddo ; endif
406 
407  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
408  if (cs%no_slip ) then
409  relative_vorticity = (2.0-g%mask2dBu(i,j)) * (dvdx(i,j) - dudy(i,j)) * &
410  g%IareaBu(i,j)
411  else
412  relative_vorticity = g%mask2dBu(i,j) * (dvdx(i,j) - dudy(i,j)) * &
413  g%IareaBu(i,j)
414  endif
415  absolute_vorticity = us%s_to_T*g%CoriolisBu(i,j) + relative_vorticity
416  ih = 0.0
417  if (area_q(i,j) > 0.0) then
418  harea_q = (harea_u(i,j) + harea_u(i,j+1)) + (harea_v(i,j) + harea_v(i+1,j))
419  ih = area_q(i,j) / (harea_q + h_neglect*area_q(i,j))
420  endif
421  q(i,j) = absolute_vorticity * ih
422  abs_vort(i,j) = absolute_vorticity
423  ih_q(i,j) = ih
424 
425  if (cs%bound_Coriolis) then
426  fv1 = absolute_vorticity*v(i+1,j,k)
427  fv2 = absolute_vorticity*v(i,j,k)
428  fu1 = -absolute_vorticity*u(i,j+1,k)
429  fu2 = -absolute_vorticity*u(i,j,k)
430  if (fv1 > fv2) then
431  max_fvq(i,j) = fv1 ; min_fvq(i,j) = fv2
432  else
433  max_fvq(i,j) = fv2 ; min_fvq(i,j) = fv1
434  endif
435  if (fu1 > fu2) then
436  max_fuq(i,j) = fu1 ; min_fuq(i,j) = fu2
437  else
438  max_fuq(i,j) = fu2 ; min_fuq(i,j) = fu1
439  endif
440  endif
441 
442  if (cs%id_rv > 0) rv(i,j,k) = relative_vorticity
443  if (cs%id_PV > 0) pv(i,j,k) = q(i,j)
444  if (associated(ad%rv_x_v) .or. associated(ad%rv_x_u)) &
445  q2(i,j) = relative_vorticity * ih
446  enddo ; enddo
447 
448  ! a, b, c, and d are combinations of neighboring potential
449  ! vorticities which form the Arakawa and Hsu vorticity advection
450  ! scheme. All are defined at u grid points.
451 
452  if (cs%Coriolis_Scheme == arakawa_hsu90) then
453  do j=jsq,jeq+1
454  do i=is-1,ieq
455  a(i,j) = (q(i,j) + (q(i+1,j) + q(i,j-1))) * c1_12
456  d(i,j) = ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) * c1_12
457  enddo
458  do i=isq,ieq
459  b(i,j) = (q(i,j) + (q(i-1,j) + q(i,j-1))) * c1_12
460  c(i,j) = ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) * c1_12
461  enddo
462  enddo
463  elseif (cs%Coriolis_Scheme == arakawa_lamb81) then
464  do j=jsq,jeq+1 ; do i=isq,ieq+1
465  a(i-1,j) = (2.0*(q(i,j) + q(i-1,j-1)) + (q(i-1,j) + q(i,j-1))) * c1_24
466  d(i-1,j) = ((q(i,j) + q(i-1,j-1)) + 2.0*(q(i-1,j) + q(i,j-1))) * c1_24
467  b(i,j) = ((q(i,j) + q(i-1,j-1)) + 2.0*(q(i-1,j) + q(i,j-1))) * c1_24
468  c(i,j) = (2.0*(q(i,j) + q(i-1,j-1)) + (q(i-1,j) + q(i,j-1))) * c1_24
469  ep_u(i,j) = ((q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
470  ep_v(i,j) = (-(q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
471  enddo ; enddo
472  elseif (cs%Coriolis_Scheme == al_blend) then
473  fe_m2 = cs%F_eff_max_blend - 2.0
474  rat_lin = 1.5 * fe_m2 / max(cs%wt_lin_blend, 1.0e-16)
475 
476  ! This allows the code to always give Sadourny Energy
477  if (cs%F_eff_max_blend <= 2.0) then ; fe_m2 = -1. ; rat_lin = -1.0 ; endif
478 
479  do j=jsq,jeq+1 ; do i=isq,ieq+1
480  min_ihq = min(ih_q(i-1,j-1), ih_q(i,j-1), ih_q(i-1,j), ih_q(i,j))
481  max_ihq = max(ih_q(i-1,j-1), ih_q(i,j-1), ih_q(i-1,j), ih_q(i,j))
482  rat_m1 = 1.0e15
483  if (max_ihq < 1.0e15*min_ihq) rat_m1 = max_ihq / min_ihq - 1.0
484  ! The weights used here are designed to keep the effective Coriolis
485  ! acceleration from any one point on its neighbors within a factor
486  ! of F_eff_max. The minimum permitted value is 2 (the factor for
487  ! Sadourny's energy conserving scheme).
488 
489  ! Determine the relative weights of Arakawa & Lamb vs. Arakawa and Hsu.
490  if (rat_m1 <= fe_m2) then ; al_wt = 1.0
491  elseif (rat_m1 < 1.5*fe_m2) then ; al_wt = 3.0*fe_m2 / rat_m1 - 2.0
492  else ; al_wt = 0.0 ; endif
493 
494  ! Determine the relative weights of Sadourny Energy vs. the other two.
495  if (rat_m1 <= 1.5*fe_m2) then ; sad_wt = 0.0
496  elseif (rat_m1 <= rat_lin) then
497  sad_wt = 1.0 - (1.5*fe_m2) / rat_m1
498  elseif (rat_m1 < 2.0*rat_lin) then
499  sad_wt = 1.0 - (cs%wt_lin_blend / rat_lin) * (rat_m1 - 2.0*rat_lin)
500  else ; sad_wt = 1.0 ; endif
501 
502  a(i-1,j) = sad_wt * 0.25 * q(i-1,j) + (1.0 - sad_wt) * &
503  ( ((2.0-al_wt)* q(i-1,j) + al_wt*q(i,j-1)) + &
504  2.0 * (q(i,j) + q(i-1,j-1)) ) * c1_24
505  d(i-1,j) = sad_wt * 0.25 * q(i-1,j-1) + (1.0 - sad_wt) * &
506  ( ((2.0-al_wt)* q(i-1,j-1) + al_wt*q(i,j)) + &
507  2.0 * (q(i-1,j) + q(i,j-1)) ) * c1_24
508  b(i,j) = sad_wt * 0.25 * q(i,j) + (1.0 - sad_wt) * &
509  ( ((2.0-al_wt)* q(i,j) + al_wt*q(i-1,j-1)) + &
510  2.0 * (q(i-1,j) + q(i,j-1)) ) * c1_24
511  c(i,j) = sad_wt * 0.25 * q(i,j-1) + (1.0 - sad_wt) * &
512  ( ((2.0-al_wt)* q(i,j-1) + al_wt*q(i-1,j)) + &
513  2.0 * (q(i,j) + q(i-1,j-1)) ) * c1_24
514  ep_u(i,j) = al_wt * ((q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
515  ep_v(i,j) = al_wt * (-(q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
516  enddo ; enddo
517  endif
518 
519  if (cs%Coriolis_En_Dis) then
520  ! c1 = 1.0-1.5*RANGE ; c2 = 1.0-RANGE ; c3 = 2.0 ; slope = 0.5
521  c1 = 1.0-1.5*0.5 ; c2 = 1.0-0.5 ; c3 = 2.0 ; slope = 0.5
522 
523  do j=jsq,jeq+1 ; do i=is-1,ie
524  uhc = uh_center(i,j)
525  uhm = uh(i,j,k)
526  ! This sometimes matters with some types of open boundary conditions.
527  if (g%dy_Cu(i,j) == 0.0) uhc = uhm
528 
529  if (abs(uhc) < 0.1*abs(uhm)) then
530  uhm = 10.0*uhc
531  elseif (abs(uhc) > c1*abs(uhm)) then
532  if (abs(uhc) < c2*abs(uhm)) then ; uhc = (3.0*uhc+(1.0-c2*3.0)*uhm)
533  elseif (abs(uhc) <= c3*abs(uhm)) then ; uhc = uhm
534  else ; uhc = slope*uhc+(1.0-c3*slope)*uhm
535  endif
536  endif
537 
538  if (uhc > uhm) then
539  uh_min(i,j) = uhm ; uh_max(i,j) = uhc
540  else
541  uh_max(i,j) = uhm ; uh_min(i,j) = uhc
542  endif
543  enddo ; enddo
544  do j=js-1,je ; do i=isq,ieq+1
545  vhc = vh_center(i,j)
546  vhm = vh(i,j,k)
547  ! This sometimes matters with some types of open boundary conditions.
548  if (g%dx_Cv(i,j) == 0.0) vhc = vhm
549 
550  if (abs(vhc) < 0.1*abs(vhm)) then
551  vhm = 10.0*vhc
552  elseif (abs(vhc) > c1*abs(vhm)) then
553  if (abs(vhc) < c2*abs(vhm)) then ; vhc = (3.0*vhc+(1.0-c2*3.0)*vhm)
554  elseif (abs(vhc) <= c3*abs(vhm)) then ; vhc = vhm
555  else ; vhc = slope*vhc+(1.0-c3*slope)*vhm
556  endif
557  endif
558 
559  if (vhc > vhm) then
560  vh_min(i,j) = vhm ; vh_max(i,j) = vhc
561  else
562  vh_max(i,j) = vhm ; vh_min(i,j) = vhc
563  endif
564  enddo ; enddo
565  endif
566 
567  ! Calculate KE and the gradient of KE
568  call gradke(u, v, h, ke, kex, key, k, obc, g, cs)
569 
570  ! Calculate the tendencies of zonal velocity due to the Coriolis
571  ! force and momentum advection. On a Cartesian grid, this is
572  ! CAu = q * vh - d(KE)/dx.
573  if (cs%Coriolis_Scheme == sadourny75_energy) then
574  if (cs%Coriolis_En_Dis) then
575  ! Energy dissipating biased scheme, Hallberg 200x
576  do j=js,je ; do i=isq,ieq
577  if (q(i,j)*u(i,j,k) == 0.0) then
578  temp1 = q(i,j) * ( (vh_max(i,j)+vh_max(i+1,j)) &
579  + (vh_min(i,j)+vh_min(i+1,j)) )*0.5
580  elseif (q(i,j)*u(i,j,k) < 0.0) then
581  temp1 = q(i,j) * (vh_max(i,j)+vh_max(i+1,j))
582  else
583  temp1 = q(i,j) * (vh_min(i,j)+vh_min(i+1,j))
584  endif
585  if (q(i,j-1)*u(i,j,k) == 0.0) then
586  temp2 = q(i,j-1) * ( (vh_max(i,j-1)+vh_max(i+1,j-1)) &
587  + (vh_min(i,j-1)+vh_min(i+1,j-1)) )*0.5
588  elseif (q(i,j-1)*u(i,j,k) < 0.0) then
589  temp2 = q(i,j-1) * (vh_max(i,j-1)+vh_max(i+1,j-1))
590  else
591  temp2 = q(i,j-1) * (vh_min(i,j-1)+vh_min(i+1,j-1))
592  endif
593  cau(i,j,k) = 0.25 * g%IdxCu(i,j) * (temp1 + temp2)
594  enddo ; enddo
595  else
596  ! Energy conserving scheme, Sadourny 1975
597  do j=js,je ; do i=isq,ieq
598  cau(i,j,k) = 0.25 * &
599  (q(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
600  q(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
601  enddo ; enddo
602  endif
603  elseif (cs%Coriolis_Scheme == sadourny75_enstro) then
604  do j=js,je ; do i=isq,ieq
605  cau(i,j,k) = 0.125 * (g%IdxCu(i,j) * (q(i,j) + q(i,j-1))) * &
606  ((vh(i+1,j,k) + vh(i,j,k)) + (vh(i,j-1,k) + vh(i+1,j-1,k)))
607  enddo ; enddo
608  elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
609  (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
610  (cs%Coriolis_Scheme == al_blend)) then
611  ! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
612  do j=js,je ; do i=isq,ieq
613  cau(i,j,k) = ((a(i,j) * vh(i+1,j,k) + &
614  c(i,j) * vh(i,j-1,k)) &
615  + (b(i,j) * vh(i,j,k) + &
616  d(i,j) * vh(i+1,j-1,k))) * g%IdxCu(i,j)
617  enddo ; enddo
618  elseif (cs%Coriolis_Scheme == robust_enstro) then
619  ! An enstrophy conserving scheme robust to vanishing layers
620  ! Note: Heffs are in lieu of h_at_v that should be returned by the
621  ! continuity solver. AJA
622  do j=js,je ; do i=isq,ieq
623  heff1 = abs(vh(i,j,k)*g%IdxCv(i,j))/(eps_vel+abs(v(i,j,k)))
624  heff1 = max(heff1,min(h(i,j,k),h(i,j+1,k)))
625  heff1 = min(heff1,max(h(i,j,k),h(i,j+1,k)))
626  heff2 = abs(vh(i,j-1,k)*g%IdxCv(i,j-1))/(eps_vel+abs(v(i,j-1,k)))
627  heff2 = max(heff2,min(h(i,j-1,k),h(i,j,k)))
628  heff2 = min(heff2,max(h(i,j-1,k),h(i,j,k)))
629  heff3 = abs(vh(i+1,j,k)*g%IdxCv(i+1,j))/(eps_vel+abs(v(i+1,j,k)))
630  heff3 = max(heff3,min(h(i+1,j,k),h(i+1,j+1,k)))
631  heff3 = min(heff3,max(h(i+1,j,k),h(i+1,j+1,k)))
632  heff4 = abs(vh(i+1,j-1,k)*g%IdxCv(i+1,j-1))/(eps_vel+abs(v(i+1,j-1,k)))
633  heff4 = max(heff4,min(h(i+1,j-1,k),h(i+1,j,k)))
634  heff4 = min(heff4,max(h(i+1,j-1,k),h(i+1,j,k)))
635  if (cs%PV_Adv_Scheme == pv_adv_centered) then
636  cau(i,j,k) = 0.5*(abs_vort(i,j)+abs_vort(i,j-1)) * &
637  ((vh(i ,j ,k)+vh(i+1,j-1,k)) + &
638  (vh(i ,j-1,k)+vh(i+1,j ,k)) ) / &
639  (h_tiny +((heff1+heff4) +(heff2+heff3)) ) * g%IdxCu(i,j)
640  elseif (cs%PV_Adv_Scheme == pv_adv_upwind1) then
641  vheff = ((vh(i ,j ,k)+vh(i+1,j-1,k)) + &
642  (vh(i ,j-1,k)+vh(i+1,j ,k)) )
643  qvheff = 0.5*( (abs_vort(i,j)+abs_vort(i,j-1))*vheff &
644  -(abs_vort(i,j)-abs_vort(i,j-1))*abs(vheff) )
645  cau(i,j,k) = qvheff / &
646  (h_tiny +((heff1+heff4) +(heff2+heff3)) ) * g%IdxCu(i,j)
647  endif
648  enddo ; enddo
649  endif
650  ! Add in the additonal terms with Arakawa & Lamb.
651  if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
652  (cs%Coriolis_Scheme == al_blend)) then ; do j=js,je ; do i=isq,ieq
653  cau(i,j,k) = cau(i,j,k) + &
654  (ep_u(i,j)*uh(i-1,j,k) - ep_u(i+1,j)*uh(i+1,j,k)) * g%IdxCu(i,j)
655  enddo ; enddo ; endif
656 
657 
658  if (cs%bound_Coriolis) then
659  do j=js,je ; do i=isq,ieq
660  max_fv = max(max_fvq(i,j),max_fvq(i,j-1))
661  min_fv = min(min_fvq(i,j),min_fvq(i,j-1))
662  ! CAu(I,j,k) = min( CAu(I,j,k), max_fv )
663  ! CAu(I,j,k) = max( CAu(I,j,k), min_fv )
664  if (cau(i,j,k) > max_fv) then
665  cau(i,j,k) = max_fv
666  else
667  if (cau(i,j,k) < min_fv) cau(i,j,k) = min_fv
668  endif
669  enddo ; enddo
670  endif
671 
672  ! Term - d(KE)/dx.
673  do j=js,je ; do i=isq,ieq
674  cau(i,j,k) = cau(i,j,k) - kex(i,j)
675  if (associated(ad%gradKEu)) ad%gradKEu(i,j,k) = -kex(i,j)
676  enddo ; enddo
677 
678 
679  ! Calculate the tendencies of meridional velocity due to the Coriolis
680  ! force and momentum advection. On a Cartesian grid, this is
681  ! CAv = - q * uh - d(KE)/dy.
682  if (cs%Coriolis_Scheme == sadourny75_energy) then
683  if (cs%Coriolis_En_Dis) then
684  ! Energy dissipating biased scheme, Hallberg 200x
685  do j=jsq,jeq ; do i=is,ie
686  if (q(i-1,j)*v(i,j,k) == 0.0) then
687  temp1 = q(i-1,j) * ( (uh_max(i-1,j)+uh_max(i-1,j+1)) &
688  + (uh_min(i-1,j)+uh_min(i-1,j+1)) )*0.5
689  elseif (q(i-1,j)*v(i,j,k) > 0.0) then
690  temp1 = q(i-1,j) * (uh_max(i-1,j)+uh_max(i-1,j+1))
691  else
692  temp1 = q(i-1,j) * (uh_min(i-1,j)+uh_min(i-1,j+1))
693  endif
694  if (q(i,j)*v(i,j,k) == 0.0) then
695  temp2 = q(i,j) * ( (uh_max(i,j)+uh_max(i,j+1)) &
696  + (uh_min(i,j)+uh_min(i,j+1)) )*0.5
697  elseif (q(i,j)*v(i,j,k) > 0.0) then
698  temp2 = q(i,j) * (uh_max(i,j)+uh_max(i,j+1))
699  else
700  temp2 = q(i,j) * (uh_min(i,j)+uh_min(i,j+1))
701  endif
702  cav(i,j,k) = - 0.25 * g%IdyCv(i,j) * (temp1 + temp2)
703  enddo ; enddo
704  else
705  ! Energy conserving scheme, Sadourny 1975
706  do j=jsq,jeq ; do i=is,ie
707  cav(i,j,k) = - 0.25* &
708  (q(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
709  q(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
710  enddo ; enddo
711  endif
712  elseif (cs%Coriolis_Scheme == sadourny75_enstro) then
713  do j=jsq,jeq ; do i=is,ie
714  cav(i,j,k) = -0.125 * (g%IdyCv(i,j) * (q(i-1,j) + q(i,j))) * &
715  ((uh(i-1,j,k) + uh(i-1,j+1,k)) + (uh(i,j,k) + uh(i,j+1,k)))
716  enddo ; enddo
717  elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
718  (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
719  (cs%Coriolis_Scheme == al_blend)) then
720  ! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
721  do j=jsq,jeq ; do i=is,ie
722  cav(i,j,k) = - ((a(i-1,j) * uh(i-1,j,k) + &
723  c(i,j+1) * uh(i,j+1,k)) &
724  + (b(i,j) * uh(i,j,k) + &
725  d(i-1,j+1) * uh(i-1,j+1,k))) * g%IdyCv(i,j)
726  enddo ; enddo
727  elseif (cs%Coriolis_Scheme == robust_enstro) then
728  ! An enstrophy conserving scheme robust to vanishing layers
729  ! Note: Heffs are in lieu of h_at_u that should be returned by the
730  ! continuity solver. AJA
731  do j=jsq,jeq ; do i=is,ie
732  heff1 = abs(uh(i,j,k)*g%IdyCu(i,j))/(eps_vel+abs(u(i,j,k)))
733  heff1 = max(heff1,min(h(i,j,k),h(i+1,j,k)))
734  heff1 = min(heff1,max(h(i,j,k),h(i+1,j,k)))
735  heff2 = abs(uh(i-1,j,k)*g%IdyCu(i-1,j))/(eps_vel+abs(u(i-1,j,k)))
736  heff2 = max(heff2,min(h(i-1,j,k),h(i,j,k)))
737  heff2 = min(heff2,max(h(i-1,j,k),h(i,j,k)))
738  heff3 = abs(uh(i,j+1,k)*g%IdyCu(i,j+1))/(eps_vel+abs(u(i,j+1,k)))
739  heff3 = max(heff3,min(h(i,j+1,k),h(i+1,j+1,k)))
740  heff3 = min(heff3,max(h(i,j+1,k),h(i+1,j+1,k)))
741  heff4 = abs(uh(i-1,j+1,k)*g%IdyCu(i-1,j+1))/(eps_vel+abs(u(i-1,j+1,k)))
742  heff4 = max(heff4,min(h(i-1,j+1,k),h(i,j+1,k)))
743  heff4 = min(heff4,max(h(i-1,j+1,k),h(i,j+1,k)))
744  if (cs%PV_Adv_Scheme == pv_adv_centered) then
745  cav(i,j,k) = - 0.5*(abs_vort(i,j)+abs_vort(i-1,j)) * &
746  ((uh(i ,j ,k)+uh(i-1,j+1,k)) + &
747  (uh(i-1,j ,k)+uh(i ,j+1,k)) ) / &
748  (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
749  elseif (cs%PV_Adv_Scheme == pv_adv_upwind1) then
750  uheff = ((uh(i ,j ,k)+uh(i-1,j+1,k)) + &
751  (uh(i-1,j ,k)+uh(i ,j+1,k)) )
752  quheff = 0.5*( (abs_vort(i,j)+abs_vort(i-1,j))*uheff &
753  -(abs_vort(i,j)-abs_vort(i-1,j))*abs(uheff) )
754  cav(i,j,k) = - quheff / &
755  (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
756  endif
757  enddo ; enddo
758  endif
759  ! Add in the additonal terms with Arakawa & Lamb.
760  if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
761  (cs%Coriolis_Scheme == al_blend)) then ; do j=jsq,jeq ; do i=is,ie
762  cav(i,j,k) = cav(i,j,k) + &
763  (ep_v(i,j)*vh(i,j-1,k) - ep_v(i,j+1)*vh(i,j+1,k)) * g%IdyCv(i,j)
764  enddo ; enddo ; endif
765 
766  if (cs%bound_Coriolis) then
767  do j=jsq,jeq ; do i=is,ie
768  max_fu = max(max_fuq(i,j),max_fuq(i-1,j))
769  min_fu = min(min_fuq(i,j),min_fuq(i-1,j))
770  if (cav(i,j,k) > max_fu) then
771  cav(i,j,k) = max_fu
772  else
773  if (cav(i,j,k) < min_fu) cav(i,j,k) = min_fu
774  endif
775  enddo ; enddo
776  endif
777 
778  ! Term - d(KE)/dy.
779  do j=jsq,jeq ; do i=is,ie
780  cav(i,j,k) = cav(i,j,k) - key(i,j)
781  if (associated(ad%gradKEv)) ad%gradKEv(i,j,k) = -key(i,j)
782  enddo ; enddo
783 
784  if (associated(ad%rv_x_u) .or. associated(ad%rv_x_v)) then
785  ! Calculate the Coriolis-like acceleration due to relative vorticity.
786  if (cs%Coriolis_Scheme == sadourny75_energy) then
787  if (associated(ad%rv_x_u)) then
788  do j=jsq,jeq ; do i=is,ie
789  ad%rv_x_u(i,j,k) = - 0.25* &
790  (q2(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
791  q2(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
792  enddo ; enddo
793  endif
794 
795  if (associated(ad%rv_x_v)) then
796  do j=js,je ; do i=isq,ieq
797  ad%rv_x_v(i,j,k) = 0.25 * &
798  (q2(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
799  q2(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
800  enddo ; enddo
801  endif
802  else
803  if (associated(ad%rv_x_u)) then
804  do j=jsq,jeq ; do i=is,ie
805  ad%rv_x_u(i,j,k) = -g%IdyCv(i,j) * c1_12 * &
806  ((q2(i,j) + q2(i-1,j) + q2(i-1,j-1)) * uh(i-1,j,k) + &
807  (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * uh(i,j,k) + &
808  (q2(i-1,j) + q2(i,j+1) + q2(i,j)) * uh(i,j+1,k) + &
809  (q2(i,j) + q2(i-1,j+1) + q2(i-1,j)) * uh(i-1,j+1,k))
810  enddo ; enddo
811  endif
812 
813  if (associated(ad%rv_x_v)) then
814  do j=js,je ; do i=isq,ieq
815  ad%rv_x_v(i,j,k) = g%IdxCu(i,j) * c1_12 * &
816  ((q2(i+1,j) + q2(i,j) + q2(i,j-1)) * vh(i+1,j,k) + &
817  (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * vh(i,j,k) + &
818  (q2(i-1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i,j-1,k) + &
819  (q2(i+1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i+1,j-1,k))
820  enddo ; enddo
821  endif
822  endif
823  endif
824 
825  enddo ! k-loop.
826 
827  ! Here the various Coriolis-related derived quantities are offered for averaging.
828  if (query_averaging_enabled(cs%diag)) then
829  if (cs%id_rv > 0) call post_data(cs%id_rv, rv, cs%diag)
830  if (cs%id_PV > 0) call post_data(cs%id_PV, pv, cs%diag)
831  if (cs%id_gKEu>0) call post_data(cs%id_gKEu, ad%gradKEu, cs%diag)
832  if (cs%id_gKEv>0) call post_data(cs%id_gKEv, ad%gradKEv, cs%diag)
833  if (cs%id_rvxu > 0) call post_data(cs%id_rvxu, ad%rv_x_u, cs%diag)
834  if (cs%id_rvxv > 0) call post_data(cs%id_rvxv, ad%rv_x_v, cs%diag)
835  endif
836 
837 end subroutine coradcalc
838 
839 
840 !> Calculates the acceleration due to the gradient of kinetic energy.
841 subroutine gradke(u, v, h, KE, KEx, KEy, k, OBC, G, CS)
842  type(ocean_grid_type), intent(in) :: G !< Ocen grid structure
843  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [m s-1]
844  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [m s-1]
845  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
846  real, dimension(SZI_(G) ,SZJ_(G) ), intent(out) :: KE !< Kinetic energy [m2 s-2]
847  real, dimension(SZIB_(G),SZJ_(G) ), intent(out) :: KEx !< Zonal acceleration due to kinetic
848  !! energy gradient [m s-2]
849  real, dimension(SZI_(G) ,SZJB_(G)), intent(out) :: KEy !< Meridional acceleration due to kinetic
850  !! energy gradient [m s-2]
851  integer, intent(in) :: k !< Layer number to calculate for
852  type(ocean_obc_type), pointer :: OBC !< Open boundary control structure
853  type(coriolisadv_cs), pointer :: CS !< Control structure for MOM_CoriolisAdv
854  ! Local variables
855  real :: um, up, vm, vp ! Temporary variables [m s-1].
856  real :: um2, up2, vm2, vp2 ! Temporary variables [m2 s-2].
857  real :: um2a, up2a, vm2a, vp2a ! Temporary variables [m4 s-2].
858  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
859 
860  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
861  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
862 
863 
864  ! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term).
865  if (cs%KE_Scheme == ke_arakawa) then
866  ! The following calculation of Kinetic energy includes the metric terms
867  ! identified in Arakawa & Lamb 1982 as important for KE conservation. It
868  ! also includes the possibility of partially-blocked tracer cell faces.
869  do j=jsq,jeq+1 ; do i=isq,ieq+1
870  ke(i,j) = ( ( g%areaCu( i ,j)*(u( i ,j,k)*u( i ,j,k)) &
871  +g%areaCu(i-1,j)*(u(i-1,j,k)*u(i-1,j,k)) ) &
872  +( g%areaCv(i, j )*(v(i, j ,k)*v(i, j ,k)) &
873  +g%areaCv(i,j-1)*(v(i,j-1,k)*v(i,j-1,k)) ) &
874  )*0.25*g%IareaT(i,j)
875  enddo ; enddo
876  elseif (cs%KE_Scheme == ke_simple_gudonov) then
877  ! The following discretization of KE is based on the one-dimensinal Gudonov
878  ! scheme which does not take into account any geometric factors
879  do j=jsq,jeq+1 ; do i=isq,ieq+1
880  up = 0.5*( u(i-1,j,k) + abs( u(i-1,j,k) ) ) ; up2 = up*up
881  um = 0.5*( u( i ,j,k) - abs( u( i ,j,k) ) ) ; um2 = um*um
882  vp = 0.5*( v(i,j-1,k) + abs( v(i,j-1,k) ) ) ; vp2 = vp*vp
883  vm = 0.5*( v(i, j ,k) - abs( v(i, j ,k) ) ) ; vm2 = vm*vm
884  ke(i,j) = ( max(up2,um2) + max(vp2,vm2) ) *0.5
885  enddo ; enddo
886  elseif (cs%KE_Scheme == ke_gudonov) then
887  ! The following discretization of KE is based on the one-dimensinal Gudonov
888  ! scheme but has been adapted to take horizontal grid factors into account
889  do j=jsq,jeq+1 ; do i=isq,ieq+1
890  up = 0.5*( u(i-1,j,k) + abs( u(i-1,j,k) ) ) ; up2a = up*up*g%areaCu(i-1,j)
891  um = 0.5*( u( i ,j,k) - abs( u( i ,j,k) ) ) ; um2a = um*um*g%areaCu( i ,j)
892  vp = 0.5*( v(i,j-1,k) + abs( v(i,j-1,k) ) ) ; vp2a = vp*vp*g%areaCv(i,j-1)
893  vm = 0.5*( v(i, j ,k) - abs( v(i, j ,k) ) ) ; vm2a = vm*vm*g%areaCv(i, j )
894  ke(i,j) = ( max(um2a,up2a) + max(vm2a,vp2a) )*0.5*g%IareaT(i,j)
895  enddo ; enddo
896  endif
897 
898  ! Term - d(KE)/dx.
899  do j=js,je ; do i=isq,ieq
900  kex(i,j) = (ke(i+1,j) - ke(i,j)) * g%IdxCu(i,j)
901  enddo ; enddo
902 
903  ! Term - d(KE)/dy.
904  do j=jsq,jeq ; do i=is,ie
905  key(i,j) = (ke(i,j+1) - ke(i,j)) * g%IdyCv(i,j)
906  enddo ; enddo
907 
908  if (associated(obc)) then
909  do n=1,obc%number_of_segments
910  if (obc%segment(n)%is_N_or_S) then
911  do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
912  key(i,obc%segment(n)%HI%JsdB) = 0.
913  enddo
914  elseif (obc%segment(n)%is_E_or_W) then
915  do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
916  kex(obc%segment(n)%HI%IsdB,j) = 0.
917  enddo
918  endif
919  enddo
920  endif
921 
922 end subroutine gradke
923 
924 !> Initializes the control structure for coriolisadv_cs
925 subroutine coriolisadv_init(Time, G, param_file, diag, AD, CS)
926  type(time_type), target, intent(in) :: time !< Current model time
927  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
928  type(param_file_type), intent(in) :: param_file !< Runtime parameter handles
929  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
930  type(accel_diag_ptrs), target, intent(inout) :: ad !< Strorage for acceleration diagnostics
931  type(coriolisadv_cs), pointer :: cs !< Control structure fro MOM_CoriolisAdv
932  ! Local variables
933 ! This include declares and sets the variable "version".
934 #include "version_variable.h"
935  character(len=40) :: mdl = "MOM_CoriolisAdv" ! This module's name.
936  character(len=20) :: tmpstr
937  character(len=400) :: mesg
938  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
939 
940  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
941  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
942 
943  if (associated(cs)) then
944  call mom_error(warning, "CoriolisAdv_init called with associated control structure.")
945  return
946  endif
947  allocate(cs)
948 
949  cs%diag => diag ; cs%Time => time
950 
951  ! Read all relevant parameters and write them to the model log.
952  call log_version(param_file, mdl, version, "")
953  call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &
954  "If true, no slip boundary conditions are used; otherwise "//&
955  "free slip boundary conditions are assumed. The "//&
956  "implementation of the free slip BCs on a C-grid is much "//&
957  "cleaner than the no slip BCs. The use of free slip BCs "//&
958  "is strongly encouraged, and no slip BCs are not used with "//&
959  "the biharmonic viscosity.", default=.false.)
960 
961  call get_param(param_file, mdl, "CORIOLIS_EN_DIS", cs%Coriolis_En_Dis, &
962  "If true, two estimates of the thickness fluxes are used "//&
963  "to estimate the Coriolis term, and the one that "//&
964  "dissipates energy relative to the other one is used.", &
965  default=.false.)
966 
967  ! Set %Coriolis_Scheme
968  ! (Select the baseline discretization for the Coriolis term)
969  call get_param(param_file, mdl, "CORIOLIS_SCHEME", tmpstr, &
970  "CORIOLIS_SCHEME selects the discretization for the "//&
971  "Coriolis terms. Valid values are: \n"//&
972  "\t SADOURNY75_ENERGY - Sadourny, 1975; energy cons. \n"//&
973  "\t ARAKAWA_HSU90 - Arakawa & Hsu, 1990 \n"//&
974  "\t SADOURNY75_ENSTRO - Sadourny, 1975; enstrophy cons. \n"//&
975  "\t ARAKAWA_LAMB81 - Arakawa & Lamb, 1981; En. + Enst.\n"//&
976  "\t ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with \n"//&
977  "\t Arakawa & Hsu and Sadourny energy", &
978  default=sadourny75_energy_string)
979  tmpstr = uppercase(tmpstr)
980  select case (tmpstr)
981  case (sadourny75_energy_string)
982  cs%Coriolis_Scheme = sadourny75_energy
983  case (arakawa_hsu_string)
984  cs%Coriolis_Scheme = arakawa_hsu90
985  case (sadourny75_enstro_string)
986  cs%Coriolis_Scheme = sadourny75_enstro
987  case (arakawa_lamb_string)
988  cs%Coriolis_Scheme = arakawa_lamb81
989  case (al_blend_string)
990  cs%Coriolis_Scheme = al_blend
991  case (robust_enstro_string)
992  cs%Coriolis_Scheme = robust_enstro
993  cs%Coriolis_En_Dis = .false.
994  case default
995  call mom_mesg('CoriolisAdv_init: Coriolis_Scheme ="'//trim(tmpstr)//'"', 0)
996  call mom_error(fatal, "CoriolisAdv_init: Unrecognized setting "// &
997  "#define CORIOLIS_SCHEME "//trim(tmpstr)//" found in input file.")
998  end select
999  if (cs%Coriolis_Scheme == al_blend) then
1000  call get_param(param_file, mdl, "CORIOLIS_BLEND_WT_LIN", cs%wt_lin_blend, &
1001  "A weighting value for the ratio of inverse thicknesses, "//&
1002  "beyond which the blending between Sadourny Energy and "//&
1003  "Arakawa & Hsu goes linearly to 0 when CORIOLIS_SCHEME "//&
1004  "is ARAWAKA_LAMB_BLEND. This must be between 1 and 1e-16.", &
1005  units="nondim", default=0.125)
1006  call get_param(param_file, mdl, "CORIOLIS_BLEND_F_EFF_MAX", cs%F_eff_max_blend, &
1007  "The factor by which the maximum effective Coriolis "//&
1008  "acceleration from any point can be increased when "//&
1009  "blending different discretizations with the "//&
1010  "ARAKAWA_LAMB_BLEND Coriolis scheme. This must be "//&
1011  "greater than 2.0 (the max value for Sadourny energy).", &
1012  units="nondim", default=4.0)
1013  cs%wt_lin_blend = min(1.0, max(cs%wt_lin_blend,1e-16))
1014  if (cs%F_eff_max_blend < 2.0) call mom_error(warning, "CoriolisAdv_init: "//&
1015  "CORIOLIS_BLEND_F_EFF_MAX should be at least 2.")
1016  endif
1017 
1018  mesg = "If true, the Coriolis terms at u-points are bounded by "//&
1019  "the four estimates of (f+rv)v from the four neighboring "//&
1020  "v-points, and similarly at v-points."
1021  if (cs%Coriolis_En_Dis .and. (cs%Coriolis_Scheme == sadourny75_energy)) then
1022  mesg = trim(mesg)//" This option is "//&
1023  "always effectively false with CORIOLIS_EN_DIS defined and "//&
1024  "CORIOLIS_SCHEME set to "//trim(sadourny75_energy_string)//"."
1025  else
1026  mesg = trim(mesg)//" This option would "//&
1027  "have no effect on the SADOURNY Coriolis scheme if it "//&
1028  "were possible to use centered difference thickness fluxes."
1029  endif
1030  call get_param(param_file, mdl, "BOUND_CORIOLIS", cs%bound_Coriolis, mesg, &
1031  default=.false.)
1032  if ((cs%Coriolis_En_Dis .and. (cs%Coriolis_Scheme == sadourny75_energy)) .or. &
1033  (cs%Coriolis_Scheme == robust_enstro)) cs%bound_Coriolis = .false.
1034 
1035  ! Set KE_Scheme (selects discretization of KE)
1036  call get_param(param_file, mdl, "KE_SCHEME", tmpstr, &
1037  "KE_SCHEME selects the discretization for acceleration "//&
1038  "due to the kinetic energy gradient. Valid values are: \n"//&
1039  "\t KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV", &
1040  default=ke_arakawa_string)
1041  tmpstr = uppercase(tmpstr)
1042  select case (tmpstr)
1043  case (ke_arakawa_string); cs%KE_Scheme = ke_arakawa
1044  case (ke_simple_gudonov_string); cs%KE_Scheme = ke_simple_gudonov
1045  case (ke_gudonov_string); cs%KE_Scheme = ke_gudonov
1046  case default
1047  call mom_mesg('CoriolisAdv_init: KE_Scheme ="'//trim(tmpstr)//'"', 0)
1048  call mom_error(fatal, "CoriolisAdv_init: "// &
1049  "#define KE_SCHEME "//trim(tmpstr)//" in input file is invalid.")
1050  end select
1051 
1052  ! Set PV_Adv_Scheme (selects discretization of PV advection)
1053  call get_param(param_file, mdl, "PV_ADV_SCHEME", tmpstr, &
1054  "PV_ADV_SCHEME selects the discretization for PV "//&
1055  "advection. Valid values are: \n"//&
1056  "\t PV_ADV_CENTERED - centered (aka Sadourny, 75) \n"//&
1057  "\t PV_ADV_UPWIND1 - upwind, first order", &
1058  default=pv_adv_centered_string)
1059  select case (uppercase(tmpstr))
1060  case (pv_adv_centered_string)
1061  cs%PV_Adv_Scheme = pv_adv_centered
1062  case (pv_adv_upwind1_string)
1063  cs%PV_Adv_Scheme = pv_adv_upwind1
1064  case default
1065  call mom_mesg('CoriolisAdv_init: PV_Adv_Scheme ="'//trim(tmpstr)//'"', 0)
1066  call mom_error(fatal, "CoriolisAdv_init: "// &
1067  "#DEFINE PV_ADV_SCHEME in input file is invalid.")
1068  end select
1069 
1070  cs%id_rv = register_diag_field('ocean_model', 'RV', diag%axesBL, time, &
1071  'Relative Vorticity', 's-1')
1072 
1073  cs%id_PV = register_diag_field('ocean_model', 'PV', diag%axesBL, time, &
1074  'Potential Vorticity', 'm-1 s-1')
1075 
1076  cs%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, time, &
1077  'Zonal Acceleration from Grad. Kinetic Energy', 'm-1 s-2')
1078  if (cs%id_gKEu > 0) call safe_alloc_ptr(ad%gradKEu,isdb,iedb,jsd,jed,nz)
1079 
1080  cs%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, time, &
1081  'Meridional Acceleration from Grad. Kinetic Energy', 'm-1 s-2')
1082  if (cs%id_gKEv > 0) call safe_alloc_ptr(ad%gradKEv,isd,ied,jsdb,jedb,nz)
1083 
1084  cs%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, time, &
1085  'Meridional Acceleration from Relative Vorticity', 'm-1 s-2')
1086  if (cs%id_rvxu > 0) call safe_alloc_ptr(ad%rv_x_u,isd,ied,jsdb,jedb,nz)
1087 
1088  cs%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, time, &
1089  'Zonal Acceleration from Relative Vorticity', 'm-1 s-2')
1090  if (cs%id_rvxv > 0) call safe_alloc_ptr(ad%rv_x_v,isdb,iedb,jsd,jed,nz)
1091 
1092 end subroutine coriolisadv_init
1093 
1094 !> Destructor for coriolisadv_cs
1095 subroutine coriolisadv_end(CS)
1096  type(coriolisadv_cs), pointer :: cs !< Control structure fro MOM_CoriolisAdv
1097  deallocate(cs)
1098 end subroutine coriolisadv_end
1099 
1100 !> \namespace mom_coriolisadv
1101 !!
1102 !! This file contains the subroutine that calculates the time
1103 !! derivatives of the velocities due to Coriolis acceleration and
1104 !! momentum advection. This subroutine uses either a vorticity
1105 !! advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or
1106 !! Sadourny's (JAS 1975) energy conserving scheme. Both have been
1107 !! modified to use general orthogonal coordinates as described in
1108 !! Arakawa and Lamb, Mon. Wea. Rev. 1981. Both schemes are second
1109 !! order accurate, and allow for vanishingly small layer thicknesses.
1110 !! The Arakawa and Hsu scheme globally conserves both total energy
1111 !! and potential enstrophy in the limit of nondivergent flow.
1112 !! Sadourny's energy conserving scheme conserves energy if the flow
1113 !! is nondivergent or centered difference thickness fluxes are used.
1114 !!
1115 !! Two sets of boundary conditions have been coded in the
1116 !! definition of relative vorticity. These are written as:
1117 !! NOSLIP defined (in spherical coordinates):
1118 !! relvort = dv/dx (east & west), with v = 0.
1119 !! relvort = -sec(Q) * d(u cos(Q))/dy (north & south), with u = 0.
1120 !!
1121 !! NOSLIP not defined (free slip):
1122 !! relvort = 0 (all boundaries)
1123 !!
1124 !! with Q temporarily defined as latitude. The free slip boundary
1125 !! condition is much more natural on a C-grid.
1126 !!
1127 !! A small fragment of the grid is shown below:
1128 !! \verbatim
1129 !!
1130 !! j+1 x ^ x ^ x At x: q, CoriolisBu
1131 !! j+1 > o > o > At ^: v, CAv, vh
1132 !! j x ^ x ^ x At >: u, CAu, uh, a, b, c, d
1133 !! j > o > o > At o: h, KE
1134 !! j-1 x ^ x ^ x
1135 !! i-1 i i+1 At x & ^:
1136 !! i i+1 At > & o:
1137 !! \endverbatim
1138 !!
1139 !! The boundaries always run through q grid points (x).
1140 
1141 end module mom_coriolisadv
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_coriolisadv::coriolisadv_cs
Control structure for mom_coriolisadv.
Definition: MOM_CoriolisAdv.F90:27
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_variables::accel_diag_ptrs
Pointers to arrays with accelerations, which can later be used for derived diagnostics,...
Definition: MOM_variables.F90:155
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:186
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_coriolisadv
Accelerations due to the Coriolis force and momentum advection.
Definition: MOM_CoriolisAdv.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:25
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239