MOM6
MOM_CoriolisAdv.F90
Go to the documentation of this file.
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 [L T-1 ~> m s-1]
115  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> 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 L2 T-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 L2 T-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 [L T-2 ~> m s-2].
123  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: cav !< Meridional acceleration due to Coriolis
124  !! and momentum advection [L T-2 ~> 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 [H-1 T-1 ~> m-1 s-1 or m2 kg-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 [L2 ~> 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 T-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 [L2 ~> 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 [L2 T-2 ~> 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 L2 ~> m3 or kg].
149  kex, & ! The zonal gradient of Kinetic energy per unit mass [L T-2 ~> m s-2],
150  ! KEx = d/dx KE.
151  uh_center ! Transport based on arithmetic mean h at u-points [H L2 T-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 L2 ~> m3 or kg].
155  key, & ! The meridonal gradient of Kinetic energy per unit mass [L T-2 ~> m s-2],
156  ! KEy = d/dy KE.
157  vh_center ! Transport based on arithmetic mean h at v-points [H L2 T-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 L2 T-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 [L2 T-1 ~> m2 s-1]
166  abs_vort, & ! Absolute vorticity at q-points [T-1 ~> s-1].
167  q2, & ! Relative vorticity over thickness [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
168  max_fvq, & ! The maximum of the adjacent values of (-u) times absolute vorticity [L T-2 ~> m s-2].
169  min_fvq, & ! The minimum of the adjacent values of (-u) times absolute vorticity [L T-2 ~> m s-2].
170  max_fuq, & ! The maximum of the adjacent values of u times absolute vorticity [L T-2 ~> m s-2].
171  min_fuq ! The minimum of the adjacent values of u times absolute vorticity [L T-2 ~> m s-2].
172  real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
173  pv, & ! A diagnostic array of the potential vorticities [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
174  rv ! A diagnostic array of the relative vorticities [T-1 ~> s-1].
175  real :: fv1, fv2, fu1, fu2 ! (f+rv)*v or (f+rv)*u [L T-2 ~> m s-2].
176  real :: max_fv, max_fu ! The maximum or minimum of the neighboring Coriolis
177  real :: min_fv, min_fu ! accelerations [L T-2 ~> 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 [T-1 ~> s-1].
182  real :: relative_vorticity ! Relative vorticity [T-1 ~> 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 L2 ~> 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 [L2 T-2 ~> m2 s-2].
190  real :: eps_vel ! A tiny, positive velocity [L T-1 ~> m s-1].
191 
192  real :: uhc, vhc ! Centered estimates of uh and vh [H L2 T-1 ~> m3 s-1 or kg s-1].
193  real :: uhm, vhm ! The input estimates of uh and vh [H L2 T-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 L2 T-1 ~> m3 s-1 or kg s-1].
210  real :: quheff,qvheff ! More temporary variables [H L2 T-1 s-1 ~> m3 s-2 or kg s-2].
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  eps_vel = 1.0e-10*us%m_s_to_L_T
224  h_tiny = gv%Angstrom_H ! Perhaps this should be set to h_neglect instead.
225 
226  !$OMP parallel do default(private) shared(Isq,Ieq,Jsq,Jeq,G,Area_h)
227  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
228  area_h(i,j) = g%mask2dT(i,j) * g%areaT(i,j)
229  enddo ; enddo
230  if (associated(obc)) then ; do n=1,obc%number_of_segments
231  if (.not. obc%segment(n)%on_pe) cycle
232  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
233  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
234  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
235  if (obc%segment(n)%direction == obc_direction_n) then
236  area_h(i,j+1) = area_h(i,j)
237  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
238  area_h(i,j) = area_h(i,j+1)
239  endif
240  enddo
241  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
242  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
243  if (obc%segment(n)%direction == obc_direction_e) then
244  area_h(i+1,j) = area_h(i,j)
245  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
246  area_h(i,j) = area_h(i+1,j)
247  endif
248  enddo
249  endif
250  enddo ; endif
251  !$OMP parallel do default(private) shared(Isq,Ieq,Jsq,Jeq,G,Area_h,Area_q)
252  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
253  area_q(i,j) = (area_h(i,j) + area_h(i+1,j+1)) + &
254  (area_h(i+1,j) + area_h(i,j+1))
255  enddo ; enddo
256 
257  !$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,&
258  !$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,h_neglect,h_tiny,OBC)
259  do k=1,nz
260 
261  ! Here the second order accurate layer potential vorticities, q,
262  ! are calculated. hq is second order accurate in space. Relative
263  ! vorticity is second order accurate everywhere with free slip b.c.s,
264  ! but only first order accurate at boundaries with no slip b.c.s.
265  ! First calculate the contributions to the circulation around the q-point.
266  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
267  dvdx(i,j) = (v(i+1,j,k)*g%dyCv(i+1,j) - v(i,j,k)*g%dyCv(i,j))
268  dudy(i,j) = (u(i,j+1,k)*g%dxCu(i,j+1) - u(i,j,k)*g%dxCu(i,j))
269  enddo ; enddo
270  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+2
271  harea_v(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i,j+1) * h(i,j+1,k))
272  enddo ; enddo
273  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+1
274  harea_u(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i+1,j) * h(i+1,j,k))
275  enddo ; enddo
276  if (cs%Coriolis_En_Dis) then
277  do j=jsq,jeq+1 ; do i=is-1,ie
278  uh_center(i,j) = 0.5 * (g%dy_Cu(i,j) * u(i,j,k)) * (h(i,j,k) + h(i+1,j,k))
279  enddo ; enddo
280  do j=js-1,je ; do i=isq,ieq+1
281  vh_center(i,j) = 0.5 * (g%dx_Cv(i,j) * v(i,j,k)) * (h(i,j,k) + h(i,j+1,k))
282  enddo ; enddo
283  endif
284 
285  ! Adjust circulation components to relative vorticity and thickness projected onto
286  ! velocity points on open boundaries.
287  if (associated(obc)) then ; do n=1,obc%number_of_segments
288  if (.not. obc%segment(n)%on_pe) cycle
289  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
290  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
291  if (obc%zero_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
292  dvdx(i,j) = 0. ; dudy(i,j) = 0.
293  enddo ; endif
294  if (obc%freeslip_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
295  dudy(i,j) = 0.
296  enddo ; endif
297  if (obc%computed_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
298  if (obc%segment(n)%direction == obc_direction_n) then
299  dudy(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%dxCu(i,j)
300  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
301  dudy(i,j) = 2.0*(u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dxCu(i,j+1)
302  endif
303  enddo ; endif
304  if (obc%specified_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
305  if (obc%segment(n)%direction == obc_direction_n) then
306  dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j)*g%dyBu(i,j)
307  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
308  dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j+1)*g%dyBu(i,j)
309  endif
310  enddo ; endif
311 
312  ! Project thicknesses across OBC points with a no-gradient condition.
313  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
314  if (obc%segment(n)%direction == obc_direction_n) then
315  harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j,k)
316  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
317  harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
318  endif
319  enddo
320 
321  if (cs%Coriolis_En_Dis) then
322  do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
323  if (obc%segment(n)%direction == obc_direction_n) then
324  vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j,k)
325  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
326  vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j+1,k)
327  endif
328  enddo
329  endif
330  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
331  if (obc%zero_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
332  dvdx(i,j) = 0. ; dudy(i,j) = 0.
333  enddo ; endif
334  if (obc%freeslip_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
335  dvdx(i,j) = 0.
336  enddo ; endif
337  if (obc%computed_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
338  if (obc%segment(n)%direction == obc_direction_e) then
339  dvdx(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%dyCv(i,j)
340  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
341  dvdx(i,j) = 2.0*(v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dyCv(i+1,j)
342  endif
343  enddo ; endif
344  if (obc%specified_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
345  if (obc%segment(n)%direction == obc_direction_e) then
346  dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i,j)*g%dxBu(i,j)
347  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
348  dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i+1,j)*g%dxBu(i,j)
349  endif
350  enddo ; endif
351 
352  ! Project thicknesses across OBC points with a no-gradient condition.
353  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
354  if (obc%segment(n)%direction == obc_direction_e) then
355  harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i,j,k)
356  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
357  harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i+1,j,k)
358  endif
359  enddo
360  if (cs%Coriolis_En_Dis) then
361  do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
362  if (obc%segment(n)%direction == obc_direction_e) then
363  uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i,j,k)
364  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
365  uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i+1,j,k)
366  endif
367  enddo
368  endif
369  endif
370  enddo ; endif
371 
372  if (associated(obc)) then ; do n=1,obc%number_of_segments
373  if (.not. obc%segment(n)%on_pe) cycle
374  ! Now project thicknesses across cell-corner points in the OBCs. The two
375  ! projections have to occur in sequence and can not be combined easily.
376  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
377  if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
378  do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
379  if (obc%segment(n)%direction == obc_direction_n) then
380  if (area_h(i,j) + area_h(i+1,j) > 0.0) then
381  harea_u(i,j+1) = harea_u(i,j) * ((area_h(i,j+1) + area_h(i+1,j+1)) / &
382  (area_h(i,j) + area_h(i+1,j)))
383  else ; harea_u(i,j+1) = 0.0 ; endif
384  else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
385  if (area_h(i,j+1) + area_h(i+1,j+1) > 0.0) then
386  harea_u(i,j) = harea_u(i,j+1) * ((area_h(i,j) + area_h(i+1,j)) / &
387  (area_h(i,j+1) + area_h(i+1,j+1)))
388  else ; harea_u(i,j) = 0.0 ; endif
389  endif
390  enddo
391  elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
392  do j = max(jsq-1,obc%segment(n)%HI%JsdB), min(jeq+1,obc%segment(n)%HI%JedB)
393  if (obc%segment(n)%direction == obc_direction_e) then
394  if (area_h(i,j) + area_h(i,j+1) > 0.0) then
395  harea_v(i+1,j) = harea_v(i,j) * ((area_h(i+1,j) + area_h(i+1,j+1)) / &
396  (area_h(i,j) + area_h(i,j+1)))
397  else ; harea_v(i+1,j) = 0.0 ; endif
398  else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
399  harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
400  if (area_h(i+1,j) + area_h(i+1,j+1) > 0.0) then
401  harea_v(i,j) = harea_v(i+1,j) * ((area_h(i,j) + area_h(i,j+1)) / &
402  (area_h(i+1,j) + area_h(i+1,j+1)))
403  else ; harea_v(i,j) = 0.0 ; endif
404  endif
405  enddo
406  endif
407  enddo ; endif
408 
409  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
410  if (cs%no_slip ) then
411  relative_vorticity = (2.0-g%mask2dBu(i,j)) * (dvdx(i,j) - dudy(i,j)) * g%IareaBu(i,j)
412  else
413  relative_vorticity = g%mask2dBu(i,j) * (dvdx(i,j) - dudy(i,j)) * g%IareaBu(i,j)
414  endif
415  absolute_vorticity = 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, us, 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) + c(i,j) * vh(i,j-1,k)) + &
614  (b(i,j) * vh(i,j,k) + d(i,j) * vh(i+1,j-1,k))) * g%IdxCu(i,j)
615  enddo ; enddo
616  elseif (cs%Coriolis_Scheme == robust_enstro) then
617  ! An enstrophy conserving scheme robust to vanishing layers
618  ! Note: Heffs are in lieu of h_at_v that should be returned by the
619  ! continuity solver. AJA
620  do j=js,je ; do i=isq,ieq
621  heff1 = abs(vh(i,j,k) * g%IdxCv(i,j)) / (eps_vel+abs(v(i,j,k)))
622  heff1 = max(heff1, min(h(i,j,k),h(i,j+1,k)))
623  heff1 = min(heff1, max(h(i,j,k),h(i,j+1,k)))
624  heff2 = abs(vh(i,j-1,k) * g%IdxCv(i,j-1)) / (eps_vel+abs(v(i,j-1,k)))
625  heff2 = max(heff2, min(h(i,j-1,k),h(i,j,k)))
626  heff2 = min(heff2, max(h(i,j-1,k),h(i,j,k)))
627  heff3 = abs(vh(i+1,j,k) * g%IdxCv(i+1,j)) / (eps_vel+abs(v(i+1,j,k)))
628  heff3 = max(heff3, min(h(i+1,j,k),h(i+1,j+1,k)))
629  heff3 = min(heff3, max(h(i+1,j,k),h(i+1,j+1,k)))
630  heff4 = abs(vh(i+1,j-1,k) * g%IdxCv(i+1,j-1)) / (eps_vel+abs(v(i+1,j-1,k)))
631  heff4 = max(heff4, min(h(i+1,j-1,k),h(i+1,j,k)))
632  heff4 = min(heff4, max(h(i+1,j-1,k),h(i+1,j,k)))
633  if (cs%PV_Adv_Scheme == pv_adv_centered) then
634  cau(i,j,k) = 0.5*(abs_vort(i,j)+abs_vort(i,j-1)) * &
635  ((vh(i,j,k) + vh(i+1,j-1,k)) + (vh(i,j-1,k) + vh(i+1,j,k)) ) / &
636  (h_tiny + ((heff1+heff4) + (heff2+heff3)) ) * g%IdxCu(i,j)
637  elseif (cs%PV_Adv_Scheme == pv_adv_upwind1) then
638  vheff = ((vh(i,j,k) + vh(i+1,j-1,k)) + (vh(i,j-1,k) + vh(i+1,j,k)) )
639  qvheff = 0.5*( (abs_vort(i,j)+abs_vort(i,j-1))*vheff &
640  -(abs_vort(i,j)-abs_vort(i,j-1))*abs(vheff) )
641  cau(i,j,k) = (qvheff / ( h_tiny + ((heff1+heff4) + (heff2+heff3)) ) ) * g%IdxCu(i,j)
642  endif
643  enddo ; enddo
644  endif
645  ! Add in the additonal terms with Arakawa & Lamb.
646  if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
647  (cs%Coriolis_Scheme == al_blend)) then ; do j=js,je ; do i=isq,ieq
648  cau(i,j,k) = cau(i,j,k) + &
649  (ep_u(i,j)*uh(i-1,j,k) - ep_u(i+1,j)*uh(i+1,j,k)) * g%IdxCu(i,j)
650  enddo ; enddo ; endif
651 
652 
653  if (cs%bound_Coriolis) then
654  do j=js,je ; do i=isq,ieq
655  max_fv = max(max_fvq(i,j), max_fvq(i,j-1))
656  min_fv = min(min_fvq(i,j), min_fvq(i,j-1))
657  ! CAu(I,j,k) = min( CAu(I,j,k), max_fv )
658  ! CAu(I,j,k) = max( CAu(I,j,k), min_fv )
659  if (cau(i,j,k) > max_fv) then
660  cau(i,j,k) = max_fv
661  else
662  if (cau(i,j,k) < min_fv) cau(i,j,k) = min_fv
663  endif
664  enddo ; enddo
665  endif
666 
667  ! Term - d(KE)/dx.
668  do j=js,je ; do i=isq,ieq
669  cau(i,j,k) = cau(i,j,k) - kex(i,j)
670  if (associated(ad%gradKEu)) ad%gradKEu(i,j,k) = -kex(i,j)
671  enddo ; enddo
672 
673 
674  ! Calculate the tendencies of meridional velocity due to the Coriolis
675  ! force and momentum advection. On a Cartesian grid, this is
676  ! CAv = - q * uh - d(KE)/dy.
677  if (cs%Coriolis_Scheme == sadourny75_energy) then
678  if (cs%Coriolis_En_Dis) then
679  ! Energy dissipating biased scheme, Hallberg 200x
680  do j=jsq,jeq ; do i=is,ie
681  if (q(i-1,j)*v(i,j,k) == 0.0) then
682  temp1 = q(i-1,j) * ( (uh_max(i-1,j)+uh_max(i-1,j+1)) &
683  + (uh_min(i-1,j)+uh_min(i-1,j+1)) )*0.5
684  elseif (q(i-1,j)*v(i,j,k) > 0.0) then
685  temp1 = q(i-1,j) * (uh_max(i-1,j)+uh_max(i-1,j+1))
686  else
687  temp1 = q(i-1,j) * (uh_min(i-1,j)+uh_min(i-1,j+1))
688  endif
689  if (q(i,j)*v(i,j,k) == 0.0) then
690  temp2 = q(i,j) * ( (uh_max(i,j)+uh_max(i,j+1)) &
691  + (uh_min(i,j)+uh_min(i,j+1)) )*0.5
692  elseif (q(i,j)*v(i,j,k) > 0.0) then
693  temp2 = q(i,j) * (uh_max(i,j)+uh_max(i,j+1))
694  else
695  temp2 = q(i,j) * (uh_min(i,j)+uh_min(i,j+1))
696  endif
697  cav(i,j,k) = -0.25 * g%IdyCv(i,j) * (temp1 + temp2)
698  enddo ; enddo
699  else
700  ! Energy conserving scheme, Sadourny 1975
701  do j=jsq,jeq ; do i=is,ie
702  cav(i,j,k) = - 0.25* &
703  (q(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
704  q(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
705  enddo ; enddo
706  endif
707  elseif (cs%Coriolis_Scheme == sadourny75_enstro) then
708  do j=jsq,jeq ; do i=is,ie
709  cav(i,j,k) = -0.125 * (g%IdyCv(i,j) * (q(i-1,j) + q(i,j))) * &
710  ((uh(i-1,j,k) + uh(i-1,j+1,k)) + (uh(i,j,k) + uh(i,j+1,k)))
711  enddo ; enddo
712  elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
713  (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
714  (cs%Coriolis_Scheme == al_blend)) then
715  ! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
716  do j=jsq,jeq ; do i=is,ie
717  cav(i,j,k) = - ((a(i-1,j) * uh(i-1,j,k) + &
718  c(i,j+1) * uh(i,j+1,k)) &
719  + (b(i,j) * uh(i,j,k) + &
720  d(i-1,j+1) * uh(i-1,j+1,k))) * g%IdyCv(i,j)
721  enddo ; enddo
722  elseif (cs%Coriolis_Scheme == robust_enstro) then
723  ! An enstrophy conserving scheme robust to vanishing layers
724  ! Note: Heffs are in lieu of h_at_u that should be returned by the
725  ! continuity solver. AJA
726  do j=jsq,jeq ; do i=is,ie
727  heff1 = abs(uh(i,j,k) * g%IdyCu(i,j)) / (eps_vel+abs(u(i,j,k)))
728  heff1 = max(heff1, min(h(i,j,k),h(i+1,j,k)))
729  heff1 = min(heff1, max(h(i,j,k),h(i+1,j,k)))
730  heff2 = abs(uh(i-1,j,k) * g%IdyCu(i-1,j)) / (eps_vel+abs(u(i-1,j,k)))
731  heff2 = max(heff2, min(h(i-1,j,k),h(i,j,k)))
732  heff2 = min(heff2, max(h(i-1,j,k),h(i,j,k)))
733  heff3 = abs(uh(i,j+1,k) * g%IdyCu(i,j+1)) / (eps_vel+abs(u(i,j+1,k)))
734  heff3 = max(heff3, min(h(i,j+1,k),h(i+1,j+1,k)))
735  heff3 = min(heff3, max(h(i,j+1,k),h(i+1,j+1,k)))
736  heff4 = abs(uh(i-1,j+1,k) * g%IdyCu(i-1,j+1)) / (eps_vel+abs(u(i-1,j+1,k)))
737  heff4 = max(heff4, min(h(i-1,j+1,k),h(i,j+1,k)))
738  heff4 = min(heff4, max(h(i-1,j+1,k),h(i,j+1,k)))
739  if (cs%PV_Adv_Scheme == pv_adv_centered) then
740  cav(i,j,k) = - 0.5*(abs_vort(i,j)+abs_vort(i-1,j)) * &
741  ((uh(i ,j ,k)+uh(i-1,j+1,k)) + &
742  (uh(i-1,j ,k)+uh(i ,j+1,k)) ) / &
743  (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
744  elseif (cs%PV_Adv_Scheme == pv_adv_upwind1) then
745  uheff = ((uh(i ,j ,k)+uh(i-1,j+1,k)) + &
746  (uh(i-1,j ,k)+uh(i ,j+1,k)) )
747  quheff = 0.5*( (abs_vort(i,j)+abs_vort(i-1,j))*uheff &
748  -(abs_vort(i,j)-abs_vort(i-1,j))*abs(uheff) )
749  cav(i,j,k) = - quheff / &
750  (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
751  endif
752  enddo ; enddo
753  endif
754  ! Add in the additonal terms with Arakawa & Lamb.
755  if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
756  (cs%Coriolis_Scheme == al_blend)) then ; do j=jsq,jeq ; do i=is,ie
757  cav(i,j,k) = cav(i,j,k) + &
758  (ep_v(i,j)*vh(i,j-1,k) - ep_v(i,j+1)*vh(i,j+1,k)) * g%IdyCv(i,j)
759  enddo ; enddo ; endif
760 
761  if (cs%bound_Coriolis) then
762  do j=jsq,jeq ; do i=is,ie
763  max_fu = max(max_fuq(i,j),max_fuq(i-1,j))
764  min_fu = min(min_fuq(i,j),min_fuq(i-1,j))
765  if (cav(i,j,k) > max_fu) then
766  cav(i,j,k) = max_fu
767  else
768  if (cav(i,j,k) < min_fu) cav(i,j,k) = min_fu
769  endif
770  enddo ; enddo
771  endif
772 
773  ! Term - d(KE)/dy.
774  do j=jsq,jeq ; do i=is,ie
775  cav(i,j,k) = cav(i,j,k) - key(i,j)
776  if (associated(ad%gradKEv)) ad%gradKEv(i,j,k) = -key(i,j)
777  enddo ; enddo
778 
779  if (associated(ad%rv_x_u) .or. associated(ad%rv_x_v)) then
780  ! Calculate the Coriolis-like acceleration due to relative vorticity.
781  if (cs%Coriolis_Scheme == sadourny75_energy) then
782  if (associated(ad%rv_x_u)) then
783  do j=jsq,jeq ; do i=is,ie
784  ad%rv_x_u(i,j,k) = - 0.25* &
785  (q2(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
786  q2(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
787  enddo ; enddo
788  endif
789 
790  if (associated(ad%rv_x_v)) then
791  do j=js,je ; do i=isq,ieq
792  ad%rv_x_v(i,j,k) = 0.25 * &
793  (q2(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
794  q2(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
795  enddo ; enddo
796  endif
797  else
798  if (associated(ad%rv_x_u)) then
799  do j=jsq,jeq ; do i=is,ie
800  ad%rv_x_u(i,j,k) = -g%IdyCv(i,j) * c1_12 * &
801  ((q2(i,j) + q2(i-1,j) + q2(i-1,j-1)) * uh(i-1,j,k) + &
802  (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * uh(i,j,k) + &
803  (q2(i-1,j) + q2(i,j+1) + q2(i,j)) * uh(i,j+1,k) + &
804  (q2(i,j) + q2(i-1,j+1) + q2(i-1,j)) * uh(i-1,j+1,k))
805  enddo ; enddo
806  endif
807 
808  if (associated(ad%rv_x_v)) then
809  do j=js,je ; do i=isq,ieq
810  ad%rv_x_v(i,j,k) = g%IdxCu(i,j) * c1_12 * &
811  ((q2(i+1,j) + q2(i,j) + q2(i,j-1)) * vh(i+1,j,k) + &
812  (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * vh(i,j,k) + &
813  (q2(i-1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i,j-1,k) + &
814  (q2(i+1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i+1,j-1,k))
815  enddo ; enddo
816  endif
817  endif
818  endif
819 
820  enddo ! k-loop.
821 
822  ! Here the various Coriolis-related derived quantities are offered for averaging.
823  if (query_averaging_enabled(cs%diag)) then
824  if (cs%id_rv > 0) call post_data(cs%id_rv, rv, cs%diag)
825  if (cs%id_PV > 0) call post_data(cs%id_PV, pv, cs%diag)
826  if (cs%id_gKEu>0) call post_data(cs%id_gKEu, ad%gradKEu, cs%diag)
827  if (cs%id_gKEv>0) call post_data(cs%id_gKEv, ad%gradKEv, cs%diag)
828  if (cs%id_rvxu > 0) call post_data(cs%id_rvxu, ad%rv_x_u, cs%diag)
829  if (cs%id_rvxv > 0) call post_data(cs%id_rvxv, ad%rv_x_v, cs%diag)
830  endif
831 
832 end subroutine coradcalc
833 
834 
835 !> Calculates the acceleration due to the gradient of kinetic energy.
836 subroutine gradke(u, v, h, KE, KEx, KEy, k, OBC, G, US, CS)
837  type(ocean_grid_type), intent(in) :: G !< Ocen grid structure
838  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
839  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
840  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
841  real, dimension(SZI_(G) ,SZJ_(G) ), intent(out) :: KE !< Kinetic energy per unit mass [L2 T-2 ~> m2 s-2]
842  real, dimension(SZIB_(G),SZJ_(G) ), intent(out) :: KEx !< Zonal acceleration due to kinetic
843  !! energy gradient [L T-2 ~> m s-2]
844  real, dimension(SZI_(G) ,SZJB_(G)), intent(out) :: KEy !< Meridional acceleration due to kinetic
845  !! energy gradient [L T-2 ~> m s-2]
846  integer, intent(in) :: k !< Layer number to calculate for
847  type(ocean_obc_type), pointer :: OBC !< Open boundary control structure
848  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
849  type(coriolisadv_cs), pointer :: CS !< Control structure for MOM_CoriolisAdv
850  ! Local variables
851  real :: um, up, vm, vp ! Temporary variables [L T-1 ~> m s-1].
852  real :: um2, up2, vm2, vp2 ! Temporary variables [L2 T-2 ~> m2 s-2].
853  real :: um2a, up2a, vm2a, vp2a ! Temporary variables [L4 T-2 ~> m4 s-2].
854  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
855 
856  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
857  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
858 
859 
860  ! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term).
861  if (cs%KE_Scheme == ke_arakawa) then
862  ! The following calculation of Kinetic energy includes the metric terms
863  ! identified in Arakawa & Lamb 1982 as important for KE conservation. It
864  ! also includes the possibility of partially-blocked tracer cell faces.
865  do j=jsq,jeq+1 ; do i=isq,ieq+1
866  ke(i,j) = ( ( g%areaCu( i ,j)*(u( i ,j,k)*u( i ,j,k)) + &
867  g%areaCu(i-1,j)*(u(i-1,j,k)*u(i-1,j,k)) ) + &
868  ( g%areaCv(i, j )*(v(i, j ,k)*v(i, j ,k)) + &
869  g%areaCv(i,j-1)*(v(i,j-1,k)*v(i,j-1,k)) ) )*0.25*g%IareaT(i,j)
870  enddo ; enddo
871  elseif (cs%KE_Scheme == ke_simple_gudonov) then
872  ! The following discretization of KE is based on the one-dimensinal Gudonov
873  ! scheme which does not take into account any geometric factors
874  do j=jsq,jeq+1 ; do i=isq,ieq+1
875  up = 0.5*( u(i-1,j,k) + abs( u(i-1,j,k) ) ) ; up2 = up*up
876  um = 0.5*( u( i ,j,k) - abs( u( i ,j,k) ) ) ; um2 = um*um
877  vp = 0.5*( v(i,j-1,k) + abs( v(i,j-1,k) ) ) ; vp2 = vp*vp
878  vm = 0.5*( v(i, j ,k) - abs( v(i, j ,k) ) ) ; vm2 = vm*vm
879  ke(i,j) = ( max(up2,um2) + max(vp2,vm2) ) *0.5
880  enddo ; enddo
881  elseif (cs%KE_Scheme == ke_gudonov) then
882  ! The following discretization of KE is based on the one-dimensinal Gudonov
883  ! scheme but has been adapted to take horizontal grid factors into account
884  do j=jsq,jeq+1 ; do i=isq,ieq+1
885  up = 0.5*( u(i-1,j,k) + abs( u(i-1,j,k) ) ) ; up2a = up*up*g%areaCu(i-1,j)
886  um = 0.5*( u( i ,j,k) - abs( u( i ,j,k) ) ) ; um2a = um*um*g%areaCu( i ,j)
887  vp = 0.5*( v(i,j-1,k) + abs( v(i,j-1,k) ) ) ; vp2a = vp*vp*g%areaCv(i,j-1)
888  vm = 0.5*( v(i, j ,k) - abs( v(i, j ,k) ) ) ; vm2a = vm*vm*g%areaCv(i, j )
889  ke(i,j) = ( max(um2a,up2a) + max(vm2a,vp2a) )*0.5*g%IareaT(i,j)
890  enddo ; enddo
891  endif
892 
893  ! Term - d(KE)/dx.
894  do j=js,je ; do i=isq,ieq
895  kex(i,j) = (ke(i+1,j) - ke(i,j)) * g%IdxCu(i,j)
896  enddo ; enddo
897 
898  ! Term - d(KE)/dy.
899  do j=jsq,jeq ; do i=is,ie
900  key(i,j) = (ke(i,j+1) - ke(i,j)) * g%IdyCv(i,j)
901  enddo ; enddo
902 
903  if (associated(obc)) then
904  do n=1,obc%number_of_segments
905  if (obc%segment(n)%is_N_or_S) then
906  do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
907  key(i,obc%segment(n)%HI%JsdB) = 0.
908  enddo
909  elseif (obc%segment(n)%is_E_or_W) then
910  do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
911  kex(obc%segment(n)%HI%IsdB,j) = 0.
912  enddo
913  endif
914  enddo
915  endif
916 
917 end subroutine gradke
918 
919 !> Initializes the control structure for coriolisadv_cs
920 subroutine coriolisadv_init(Time, G, GV, US, param_file, diag, AD, CS)
921  type(time_type), target, intent(in) :: time !< Current model time
922  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
923  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
924  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
925  type(param_file_type), intent(in) :: param_file !< Runtime parameter handles
926  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
927  type(accel_diag_ptrs), target, intent(inout) :: ad !< Strorage for acceleration diagnostics
928  type(coriolisadv_cs), pointer :: cs !< Control structure fro MOM_CoriolisAdv
929  ! Local variables
930 ! This include declares and sets the variable "version".
931 #include "version_variable.h"
932  character(len=40) :: mdl = "MOM_CoriolisAdv" ! This module's name.
933  character(len=20) :: tmpstr
934  character(len=400) :: mesg
935  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
936 
937  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
938  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
939 
940  if (associated(cs)) then
941  call mom_error(warning, "CoriolisAdv_init called with associated control structure.")
942  return
943  endif
944  allocate(cs)
945 
946  cs%diag => diag ; cs%Time => time
947 
948  ! Read all relevant parameters and write them to the model log.
949  call log_version(param_file, mdl, version, "")
950  call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &
951  "If true, no slip boundary conditions are used; otherwise "//&
952  "free slip boundary conditions are assumed. The "//&
953  "implementation of the free slip BCs on a C-grid is much "//&
954  "cleaner than the no slip BCs. The use of free slip BCs "//&
955  "is strongly encouraged, and no slip BCs are not used with "//&
956  "the biharmonic viscosity.", default=.false.)
957 
958  call get_param(param_file, mdl, "CORIOLIS_EN_DIS", cs%Coriolis_En_Dis, &
959  "If true, two estimates of the thickness fluxes are used "//&
960  "to estimate the Coriolis term, and the one that "//&
961  "dissipates energy relative to the other one is used.", &
962  default=.false.)
963 
964  ! Set %Coriolis_Scheme
965  ! (Select the baseline discretization for the Coriolis term)
966  call get_param(param_file, mdl, "CORIOLIS_SCHEME", tmpstr, &
967  "CORIOLIS_SCHEME selects the discretization for the "//&
968  "Coriolis terms. Valid values are: \n"//&
969  "\t SADOURNY75_ENERGY - Sadourny, 1975; energy cons. \n"//&
970  "\t ARAKAWA_HSU90 - Arakawa & Hsu, 1990 \n"//&
971  "\t SADOURNY75_ENSTRO - Sadourny, 1975; enstrophy cons. \n"//&
972  "\t ARAKAWA_LAMB81 - Arakawa & Lamb, 1981; En. + Enst.\n"//&
973  "\t ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with \n"//&
974  "\t Arakawa & Hsu and Sadourny energy", &
975  default=sadourny75_energy_string)
976  tmpstr = uppercase(tmpstr)
977  select case (tmpstr)
978  case (sadourny75_energy_string)
979  cs%Coriolis_Scheme = sadourny75_energy
980  case (arakawa_hsu_string)
981  cs%Coriolis_Scheme = arakawa_hsu90
982  case (sadourny75_enstro_string)
983  cs%Coriolis_Scheme = sadourny75_enstro
984  case (arakawa_lamb_string)
985  cs%Coriolis_Scheme = arakawa_lamb81
986  case (al_blend_string)
987  cs%Coriolis_Scheme = al_blend
988  case (robust_enstro_string)
989  cs%Coriolis_Scheme = robust_enstro
990  cs%Coriolis_En_Dis = .false.
991  case default
992  call mom_mesg('CoriolisAdv_init: Coriolis_Scheme ="'//trim(tmpstr)//'"', 0)
993  call mom_error(fatal, "CoriolisAdv_init: Unrecognized setting "// &
994  "#define CORIOLIS_SCHEME "//trim(tmpstr)//" found in input file.")
995  end select
996  if (cs%Coriolis_Scheme == al_blend) then
997  call get_param(param_file, mdl, "CORIOLIS_BLEND_WT_LIN", cs%wt_lin_blend, &
998  "A weighting value for the ratio of inverse thicknesses, "//&
999  "beyond which the blending between Sadourny Energy and "//&
1000  "Arakawa & Hsu goes linearly to 0 when CORIOLIS_SCHEME "//&
1001  "is ARAWAKA_LAMB_BLEND. This must be between 1 and 1e-16.", &
1002  units="nondim", default=0.125)
1003  call get_param(param_file, mdl, "CORIOLIS_BLEND_F_EFF_MAX", cs%F_eff_max_blend, &
1004  "The factor by which the maximum effective Coriolis "//&
1005  "acceleration from any point can be increased when "//&
1006  "blending different discretizations with the "//&
1007  "ARAKAWA_LAMB_BLEND Coriolis scheme. This must be "//&
1008  "greater than 2.0 (the max value for Sadourny energy).", &
1009  units="nondim", default=4.0)
1010  cs%wt_lin_blend = min(1.0, max(cs%wt_lin_blend,1e-16))
1011  if (cs%F_eff_max_blend < 2.0) call mom_error(warning, "CoriolisAdv_init: "//&
1012  "CORIOLIS_BLEND_F_EFF_MAX should be at least 2.")
1013  endif
1014 
1015  mesg = "If true, the Coriolis terms at u-points are bounded by "//&
1016  "the four estimates of (f+rv)v from the four neighboring "//&
1017  "v-points, and similarly at v-points."
1018  if (cs%Coriolis_En_Dis .and. (cs%Coriolis_Scheme == sadourny75_energy)) then
1019  mesg = trim(mesg)//" This option is "//&
1020  "always effectively false with CORIOLIS_EN_DIS defined and "//&
1021  "CORIOLIS_SCHEME set to "//trim(sadourny75_energy_string)//"."
1022  else
1023  mesg = trim(mesg)//" This option would "//&
1024  "have no effect on the SADOURNY Coriolis scheme if it "//&
1025  "were possible to use centered difference thickness fluxes."
1026  endif
1027  call get_param(param_file, mdl, "BOUND_CORIOLIS", cs%bound_Coriolis, mesg, &
1028  default=.false.)
1029  if ((cs%Coriolis_En_Dis .and. (cs%Coriolis_Scheme == sadourny75_energy)) .or. &
1030  (cs%Coriolis_Scheme == robust_enstro)) cs%bound_Coriolis = .false.
1031 
1032  ! Set KE_Scheme (selects discretization of KE)
1033  call get_param(param_file, mdl, "KE_SCHEME", tmpstr, &
1034  "KE_SCHEME selects the discretization for acceleration "//&
1035  "due to the kinetic energy gradient. Valid values are: \n"//&
1036  "\t KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV", &
1037  default=ke_arakawa_string)
1038  tmpstr = uppercase(tmpstr)
1039  select case (tmpstr)
1040  case (ke_arakawa_string); cs%KE_Scheme = ke_arakawa
1041  case (ke_simple_gudonov_string); cs%KE_Scheme = ke_simple_gudonov
1042  case (ke_gudonov_string); cs%KE_Scheme = ke_gudonov
1043  case default
1044  call mom_mesg('CoriolisAdv_init: KE_Scheme ="'//trim(tmpstr)//'"', 0)
1045  call mom_error(fatal, "CoriolisAdv_init: "// &
1046  "#define KE_SCHEME "//trim(tmpstr)//" in input file is invalid.")
1047  end select
1048 
1049  ! Set PV_Adv_Scheme (selects discretization of PV advection)
1050  call get_param(param_file, mdl, "PV_ADV_SCHEME", tmpstr, &
1051  "PV_ADV_SCHEME selects the discretization for PV "//&
1052  "advection. Valid values are: \n"//&
1053  "\t PV_ADV_CENTERED - centered (aka Sadourny, 75) \n"//&
1054  "\t PV_ADV_UPWIND1 - upwind, first order", &
1055  default=pv_adv_centered_string)
1056  select case (uppercase(tmpstr))
1057  case (pv_adv_centered_string)
1058  cs%PV_Adv_Scheme = pv_adv_centered
1059  case (pv_adv_upwind1_string)
1060  cs%PV_Adv_Scheme = pv_adv_upwind1
1061  case default
1062  call mom_mesg('CoriolisAdv_init: PV_Adv_Scheme ="'//trim(tmpstr)//'"', 0)
1063  call mom_error(fatal, "CoriolisAdv_init: "// &
1064  "#DEFINE PV_ADV_SCHEME in input file is invalid.")
1065  end select
1066 
1067  cs%id_rv = register_diag_field('ocean_model', 'RV', diag%axesBL, time, &
1068  'Relative Vorticity', 's-1', conversion=us%s_to_T)
1069 
1070  cs%id_PV = register_diag_field('ocean_model', 'PV', diag%axesBL, time, &
1071  'Potential Vorticity', 'm-1 s-1', conversion=gv%m_to_H*us%s_to_T)
1072 
1073  cs%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, time, &
1074  'Zonal Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
1075  if (cs%id_gKEu > 0) call safe_alloc_ptr(ad%gradKEu,isdb,iedb,jsd,jed,nz)
1076 
1077  cs%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, time, &
1078  'Meridional Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
1079  if (cs%id_gKEv > 0) call safe_alloc_ptr(ad%gradKEv,isd,ied,jsdb,jedb,nz)
1080 
1081  cs%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, time, &
1082  'Meridional Acceleration from Relative Vorticity', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
1083  if (cs%id_rvxu > 0) call safe_alloc_ptr(ad%rv_x_u,isd,ied,jsdb,jedb,nz)
1084 
1085  cs%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, time, &
1086  'Zonal Acceleration from Relative Vorticity', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
1087  if (cs%id_rvxv > 0) call safe_alloc_ptr(ad%rv_x_v,isdb,iedb,jsd,jed,nz)
1088 
1089 end subroutine coriolisadv_init
1090 
1091 !> Destructor for coriolisadv_cs
1092 subroutine coriolisadv_end(CS)
1093  type(coriolisadv_cs), pointer :: cs !< Control structure fro MOM_CoriolisAdv
1094  deallocate(cs)
1095 end subroutine coriolisadv_end
1096 
1097 !> \namespace mom_coriolisadv
1098 !!
1099 !! This file contains the subroutine that calculates the time
1100 !! derivatives of the velocities due to Coriolis acceleration and
1101 !! momentum advection. This subroutine uses either a vorticity
1102 !! advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or
1103 !! Sadourny's (JAS 1975) energy conserving scheme. Both have been
1104 !! modified to use general orthogonal coordinates as described in
1105 !! Arakawa and Lamb, Mon. Wea. Rev. 1981. Both schemes are second
1106 !! order accurate, and allow for vanishingly small layer thicknesses.
1107 !! The Arakawa and Hsu scheme globally conserves both total energy
1108 !! and potential enstrophy in the limit of nondivergent flow.
1109 !! Sadourny's energy conserving scheme conserves energy if the flow
1110 !! is nondivergent or centered difference thickness fluxes are used.
1111 !!
1112 !! Two sets of boundary conditions have been coded in the
1113 !! definition of relative vorticity. These are written as:
1114 !! NOSLIP defined (in spherical coordinates):
1115 !! relvort = dv/dx (east & west), with v = 0.
1116 !! relvort = -sec(Q) * d(u cos(Q))/dy (north & south), with u = 0.
1117 !!
1118 !! NOSLIP not defined (free slip):
1119 !! relvort = 0 (all boundaries)
1120 !!
1121 !! with Q temporarily defined as latitude. The free slip boundary
1122 !! condition is much more natural on a C-grid.
1123 !!
1124 !! A small fragment of the grid is shown below:
1125 !! \verbatim
1126 !!
1127 !! j+1 x ^ x ^ x At x: q, CoriolisBu
1128 !! j+1 > o > o > At ^: v, CAv, vh
1129 !! j x ^ x ^ x At >: u, CAu, uh, a, b, c, d
1130 !! j > o > o > At o: h, KE
1131 !! j-1 x ^ x ^ x
1132 !! i-1 i i+1 At x & ^:
1133 !! i i+1 At > & o:
1134 !! \endverbatim
1135 !!
1136 !! The boundaries always run through q grid points (x).
1137 
1138 end module mom_coriolisadv
mom_coriolisadv::sadourny75_energy_string
character *(20), parameter sadourny75_energy_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:86
mom_open_boundary::obc_direction_n
integer, parameter, public obc_direction_n
Indicates the boundary is an effective northern boundary.
Definition: MOM_open_boundary.F90:63
mom_coriolisadv::arakawa_hsu_string
character *(20), parameter arakawa_hsu_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:87
mom_open_boundary::obc_direction_s
integer, parameter, public obc_direction_s
Indicates the boundary is an effective southern boundary.
Definition: MOM_open_boundary.F90:64
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
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_coriolisadv::pv_adv_upwind1
integer, parameter pv_adv_upwind1
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:103
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_open_boundary::obc_direction_e
integer, parameter, public obc_direction_e
Indicates the boundary is an effective eastern boundary.
Definition: MOM_open_boundary.F90:65
mom_coriolisadv::al_blend
integer, parameter al_blend
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:85
mom_open_boundary::obc_direction_w
integer, parameter, public obc_direction_w
Indicates the boundary is an effective western boundary.
Definition: MOM_open_boundary.F90:66
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_coriolisadv::robust_enstro_string
character *(20), parameter robust_enstro_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:88
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
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_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:42
mom_coriolisadv::ke_arakawa
integer, parameter ke_arakawa
Enumeration values for KE_Scheme.
Definition: MOM_CoriolisAdv.F90:94
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_coriolisadv::ke_simple_gudonov
integer, parameter ke_simple_gudonov
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:95
mom_coriolisadv::al_blend_string
character *(20), parameter al_blend_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:91
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_coriolisadv::ke_arakawa_string
character *(20), parameter ke_arakawa_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:97
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_coriolisadv::arakawa_lamb_string
character *(20), parameter arakawa_lamb_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:90
mom_coriolisadv::sadourny75_energy
integer, parameter sadourny75_energy
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:80
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_coriolisadv::arakawa_lamb81
integer, parameter arakawa_lamb81
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:84
mom_coriolisadv::pv_adv_centered_string
character *(20), parameter pv_adv_centered_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:104
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:151
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:195
mom_coriolisadv::ke_simple_gudonov_string
character *(20), parameter ke_simple_gudonov_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:98
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_coriolisadv::robust_enstro
integer, parameter robust_enstro
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:82
mom_coriolisadv::pv_adv_centered
integer, parameter pv_adv_centered
Enumeration values for PV_Adv_Scheme.
Definition: MOM_CoriolisAdv.F90:102
mom_coriolisadv::ke_gudonov
integer, parameter ke_gudonov
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:96
mom_coriolisadv::coriolisadv_end
subroutine, public coriolisadv_end(CS)
Destructor for coriolisadv_cs.
Definition: MOM_CoriolisAdv.F90:1093
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_coriolisadv::coradcalc
subroutine, public coradcalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
Calculates the Coriolis and momentum advection contributions to the acceleration.
Definition: MOM_CoriolisAdv.F90:112
mom_coriolisadv::gradke
subroutine gradke(u, v, h, KE, KEx, KEy, k, OBC, G, US, CS)
Calculates the acceleration due to the gradient of kinetic energy.
Definition: MOM_CoriolisAdv.F90:837
mom_coriolisadv::sadourny75_enstro
integer, parameter sadourny75_enstro
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:83
mom_coriolisadv
Accelerations due to the Coriolis force and momentum advection.
Definition: MOM_CoriolisAdv.F90:2
mom_coriolisadv::pv_adv_upwind1_string
character *(20), parameter pv_adv_upwind1_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:105
mom_coriolisadv::coriolisadv_init
subroutine, public coriolisadv_init(Time, G, GV, US, param_file, diag, AD, CS)
Initializes the control structure for coriolisadv_cs.
Definition: MOM_CoriolisAdv.F90:921
mom_coriolisadv::sadourny75_enstro_string
character *(20), parameter sadourny75_enstro_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:89
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_coriolisadv::arakawa_hsu90
integer, parameter arakawa_hsu90
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:81
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_coriolisadv::ke_gudonov_string
character *(20), parameter ke_gudonov_string
Enumeration values for Coriolis_Scheme.
Definition: MOM_CoriolisAdv.F90:99
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