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