MOM6
Data Types
mom_coriolisadv Module Reference

Detailed Description

Accelerations due to the Coriolis force and momentum advection.

This file contains the subroutine that calculates the time derivatives of the velocities due to Coriolis acceleration and momentum advection. This subroutine uses either a vorticity advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or Sadourny's (JAS 1975) energy conserving scheme. Both have been modified to use general orthogonal coordinates as described in Arakawa and Lamb, Mon. Wea. Rev. 1981. Both schemes are second order accurate, and allow for vanishingly small layer thicknesses. The Arakawa and Hsu scheme globally conserves both total energy and potential enstrophy in the limit of nondivergent flow. Sadourny's energy conserving scheme conserves energy if the flow is nondivergent or centered difference thickness fluxes are used.

Two sets of boundary conditions have been coded in the definition of relative vorticity. These are written as: NOSLIP defined (in spherical coordinates): relvort = dv/dx (east & west), with v = 0. relvort = -sec(Q) * d(u cos(Q))/dy (north & south), with u = 0.

NOSLIP not defined (free slip): relvort = 0 (all boundaries)

with Q temporarily defined as latitude. The free slip boundary condition is much more natural on a C-grid.

A small fragment of the grid is shown below: 

    j+1  x ^ x ^ x   At x:  q, CoriolisBu
    j+1  > o > o >   At ^:  v, CAv, vh
    j    x ^ x ^ x   At >:  u, CAu, uh, a, b, c, d
    j    > o > o >   At o:  h, KE
    j-1  x ^ x ^ x
        i-1  i  i+1  At x & ^:
           i  i+1    At > & o:

The boundaries always run through q grid points (x).

Data Types

type  coriolisadv_cs
 Control structure for mom_coriolisadv. More...
 
integer, parameter sadourny75_energy = 1
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter arakawa_hsu90 = 2
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter robust_enstro = 3
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter sadourny75_enstro = 4
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter arakawa_lamb81 = 5
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter al_blend = 6
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter sadourny75_energy_string = "SADOURNY75_ENERGY"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter arakawa_hsu_string = "ARAKAWA_HSU90"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter robust_enstro_string = "ROBUST_ENSTRO"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter sadourny75_enstro_string = "SADOURNY75_ENSTRO"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter arakawa_lamb_string = "ARAKAWA_LAMB81"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter al_blend_string = "ARAKAWA_LAMB_BLEND"
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter ke_arakawa = 10
 Enumeration values for KE_Scheme. More...
 
integer, parameter ke_simple_gudonov = 11
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter ke_gudonov = 12
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter ke_arakawa_string = "KE_ARAKAWA"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter ke_simple_gudonov_string = "KE_SIMPLE_GUDONOV"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter ke_gudonov_string = "KE_GUDONOV"
 Enumeration values for Coriolis_Scheme. More...
 
integer, parameter pv_adv_centered = 21
 Enumeration values for PV_Adv_Scheme. More...
 
integer, parameter pv_adv_upwind1 = 22
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter pv_adv_centered_string = "PV_ADV_CENTERED"
 Enumeration values for Coriolis_Scheme. More...
 
character *(20), parameter pv_adv_upwind1_string = "PV_ADV_UPWIND1"
 Enumeration values for Coriolis_Scheme. More...
 
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. More...
 
subroutine gradke (u, v, h, KE, KEx, KEy, k, OBC, G, US, CS)
 Calculates the acceleration due to the gradient of kinetic energy. More...
 
subroutine, public coriolisadv_init (Time, G, GV, US, param_file, diag, AD, CS)
 Initializes the control structure for coriolisadv_cs. More...
 
subroutine, public coriolisadv_end (CS)
 Destructor for coriolisadv_cs. More...
 

Function/Subroutine Documentation

◆ coradcalc()

subroutine, public mom_coriolisadv::coradcalc ( real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  uh,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  vh,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)  CAu,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)  CAv,
type(ocean_obc_type), pointer  OBC,
type(accel_diag_ptrs), intent(inout)  AD,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(coriolisadv_cs), pointer  CS 
)

Calculates the Coriolis and momentum advection contributions to the acceleration.

Parameters
[in]gOcen grid structure
[in]gvVertical grid structure
[in]uZonal velocity [L T-1 ~> m s-1]
[in]vMeridional velocity [L T-1 ~> m s-1]
[in]hLayer thickness [H ~> m or kg m-2]
[in]uhZonal transport u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]
[in]vhMeridional transport v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]
[out]cauZonal acceleration due to Coriolis and momentum advection [L T-2 ~> m s-2].
[out]cavMeridional acceleration due to Coriolis and momentum advection [L T-2 ~> m s-2].
obcOpen boundary control structure
[in,out]adStorage for acceleration diagnostics
[in]usA dimensional unit scaling type
csControl structure for MOM_CoriolisAdv

Definition at line 112 of file MOM_CoriolisAdv.F90.

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 

References al_blend, arakawa_hsu90, arakawa_lamb81, gradke(), mom_error_handler::mom_error(), mom_open_boundary::obc_direction_n, pv_adv_centered, pv_adv_upwind1, robust_enstro, sadourny75_energy, and sadourny75_enstro.

Referenced by mom_dynamics_split_rk2::step_mom_dyn_split_rk2(), mom_dynamics_unsplit::step_mom_dyn_unsplit(), and mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ coriolisadv_end()

subroutine, public mom_coriolisadv::coriolisadv_end ( type(coriolisadv_cs), pointer  CS)

Destructor for coriolisadv_cs.

Parameters
csControl structure fro MOM_CoriolisAdv

Definition at line 1093 of file MOM_CoriolisAdv.F90.

1093  type(CoriolisAdv_CS), pointer :: CS !< Control structure fro MOM_CoriolisAdv
1094  deallocate(cs)

◆ coriolisadv_init()

subroutine, public mom_coriolisadv::coriolisadv_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(accel_diag_ptrs), intent(inout), target  AD,
type(coriolisadv_cs), pointer  CS 
)

Initializes the control structure for coriolisadv_cs.

Parameters
[in]timeCurrent model time
[in]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]param_fileRuntime parameter handles
[in,out]diagDiagnostics control structure
[in,out]adStrorage for acceleration diagnostics
csControl structure fro MOM_CoriolisAdv

Definition at line 921 of file MOM_CoriolisAdv.F90.

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 

References al_blend, al_blend_string, arakawa_hsu90, arakawa_hsu_string, arakawa_lamb81, arakawa_lamb_string, ke_arakawa, ke_arakawa_string, ke_gudonov, ke_gudonov_string, ke_simple_gudonov, ke_simple_gudonov_string, mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), pv_adv_centered, pv_adv_centered_string, pv_adv_upwind1, pv_adv_upwind1_string, mom_diag_mediator::register_diag_field(), robust_enstro, robust_enstro_string, sadourny75_energy, sadourny75_energy_string, sadourny75_enstro, sadourny75_enstro_string, and mom_string_functions::uppercase().

Referenced by mom_dynamics_split_rk2::initialize_dyn_split_rk2(), mom_dynamics_unsplit::initialize_dyn_unsplit(), and mom_dynamics_unsplit_rk2::initialize_dyn_unsplit_rk2().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ gradke()

subroutine mom_coriolisadv::gradke ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied , g %jsd: g %jed ), intent(out)  KE,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed ), intent(out)  KEx,
real, dimension( g %isd: g %ied , g %jsdb: g %jedb), intent(out)  KEy,
integer, intent(in)  k,
type(ocean_obc_type), pointer  OBC,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(coriolisadv_cs), pointer  CS 
)
private

Calculates the acceleration due to the gradient of kinetic energy.

Parameters
[in]gOcen grid structure
[in]uZonal velocity [L T-1 ~> m s-1]
[in]vMeridional velocity [L T-1 ~> m s-1]
[in]hLayer thickness [H ~> m or kg m-2]
[out]keKinetic energy per unit mass [L2 T-2 ~> m2 s-2]
[out]kexZonal acceleration due to kinetic energy gradient [L T-2 ~> m s-2]
[out]keyMeridional acceleration due to kinetic energy gradient [L T-2 ~> m s-2]
[in]kLayer number to calculate for
obcOpen boundary control structure
[in]usA dimensional unit scaling type
csControl structure for MOM_CoriolisAdv

Definition at line 837 of file MOM_CoriolisAdv.F90.

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 

References ke_arakawa, ke_gudonov, and ke_simple_gudonov.

Referenced by coradcalc().

Here is the caller graph for this function:

Variable Documentation

◆ al_blend

integer, parameter mom_coriolisadv::al_blend = 6
private

Enumeration values for Coriolis_Scheme.

Definition at line 85 of file MOM_CoriolisAdv.F90.

85 integer, parameter :: AL_BLEND = 6

Referenced by coradcalc(), and coriolisadv_init().

◆ al_blend_string

character*(20), parameter mom_coriolisadv::al_blend_string = "ARAKAWA_LAMB_BLEND"
private

Enumeration values for Coriolis_Scheme.

Definition at line 91 of file MOM_CoriolisAdv.F90.

91 character*(20), parameter :: AL_BLEND_STRING = "ARAKAWA_LAMB_BLEND"

Referenced by coriolisadv_init().

◆ arakawa_hsu90

integer, parameter mom_coriolisadv::arakawa_hsu90 = 2
private

Enumeration values for Coriolis_Scheme.

Definition at line 81 of file MOM_CoriolisAdv.F90.

81 integer, parameter :: ARAKAWA_HSU90 = 2

Referenced by coradcalc(), and coriolisadv_init().

◆ arakawa_hsu_string

character*(20), parameter mom_coriolisadv::arakawa_hsu_string = "ARAKAWA_HSU90"
private

Enumeration values for Coriolis_Scheme.

Definition at line 87 of file MOM_CoriolisAdv.F90.

87 character*(20), parameter :: ARAKAWA_HSU_STRING = "ARAKAWA_HSU90"

Referenced by coriolisadv_init().

◆ arakawa_lamb81

integer, parameter mom_coriolisadv::arakawa_lamb81 = 5
private

Enumeration values for Coriolis_Scheme.

Definition at line 84 of file MOM_CoriolisAdv.F90.

84 integer, parameter :: ARAKAWA_LAMB81 = 5

Referenced by coradcalc(), and coriolisadv_init().

◆ arakawa_lamb_string

character*(20), parameter mom_coriolisadv::arakawa_lamb_string = "ARAKAWA_LAMB81"
private

Enumeration values for Coriolis_Scheme.

Definition at line 90 of file MOM_CoriolisAdv.F90.

90 character*(20), parameter :: ARAKAWA_LAMB_STRING = "ARAKAWA_LAMB81"

Referenced by coriolisadv_init().

◆ ke_arakawa

integer, parameter mom_coriolisadv::ke_arakawa = 10
private

Enumeration values for KE_Scheme.

Definition at line 94 of file MOM_CoriolisAdv.F90.

94 integer, parameter :: KE_ARAKAWA = 10

Referenced by coriolisadv_init(), and gradke().

◆ ke_arakawa_string

character*(20), parameter mom_coriolisadv::ke_arakawa_string = "KE_ARAKAWA"
private

Enumeration values for Coriolis_Scheme.

Definition at line 97 of file MOM_CoriolisAdv.F90.

97 character*(20), parameter :: KE_ARAKAWA_STRING = "KE_ARAKAWA"

Referenced by coriolisadv_init().

◆ ke_gudonov

integer, parameter mom_coriolisadv::ke_gudonov = 12
private

Enumeration values for Coriolis_Scheme.

Definition at line 96 of file MOM_CoriolisAdv.F90.

96 integer, parameter :: KE_GUDONOV = 12

Referenced by coriolisadv_init(), and gradke().

◆ ke_gudonov_string

character*(20), parameter mom_coriolisadv::ke_gudonov_string = "KE_GUDONOV"
private

Enumeration values for Coriolis_Scheme.

Definition at line 99 of file MOM_CoriolisAdv.F90.

99 character*(20), parameter :: KE_GUDONOV_STRING = "KE_GUDONOV"

Referenced by coriolisadv_init().

◆ ke_simple_gudonov

integer, parameter mom_coriolisadv::ke_simple_gudonov = 11
private

Enumeration values for Coriolis_Scheme.

Definition at line 95 of file MOM_CoriolisAdv.F90.

95 integer, parameter :: KE_SIMPLE_GUDONOV = 11

Referenced by coriolisadv_init(), and gradke().

◆ ke_simple_gudonov_string

character*(20), parameter mom_coriolisadv::ke_simple_gudonov_string = "KE_SIMPLE_GUDONOV"
private

Enumeration values for Coriolis_Scheme.

Definition at line 98 of file MOM_CoriolisAdv.F90.

98 character*(20), parameter :: KE_SIMPLE_GUDONOV_STRING = "KE_SIMPLE_GUDONOV"

Referenced by coriolisadv_init().

◆ pv_adv_centered

integer, parameter mom_coriolisadv::pv_adv_centered = 21
private

Enumeration values for PV_Adv_Scheme.

Definition at line 102 of file MOM_CoriolisAdv.F90.

102 integer, parameter :: PV_ADV_CENTERED = 21

Referenced by coradcalc(), and coriolisadv_init().

◆ pv_adv_centered_string

character*(20), parameter mom_coriolisadv::pv_adv_centered_string = "PV_ADV_CENTERED"
private

Enumeration values for Coriolis_Scheme.

Definition at line 104 of file MOM_CoriolisAdv.F90.

104 character*(20), parameter :: PV_ADV_CENTERED_STRING = "PV_ADV_CENTERED"

Referenced by coriolisadv_init().

◆ pv_adv_upwind1

integer, parameter mom_coriolisadv::pv_adv_upwind1 = 22
private

Enumeration values for Coriolis_Scheme.

Definition at line 103 of file MOM_CoriolisAdv.F90.

103 integer, parameter :: PV_ADV_UPWIND1 = 22

Referenced by coradcalc(), and coriolisadv_init().

◆ pv_adv_upwind1_string

character*(20), parameter mom_coriolisadv::pv_adv_upwind1_string = "PV_ADV_UPWIND1"
private

Enumeration values for Coriolis_Scheme.

Definition at line 105 of file MOM_CoriolisAdv.F90.

105 character*(20), parameter :: PV_ADV_UPWIND1_STRING = "PV_ADV_UPWIND1"

Referenced by coriolisadv_init().

◆ robust_enstro

integer, parameter mom_coriolisadv::robust_enstro = 3
private

Enumeration values for Coriolis_Scheme.

Definition at line 82 of file MOM_CoriolisAdv.F90.

82 integer, parameter :: ROBUST_ENSTRO = 3

Referenced by coradcalc(), and coriolisadv_init().

◆ robust_enstro_string

character*(20), parameter mom_coriolisadv::robust_enstro_string = "ROBUST_ENSTRO"
private

Enumeration values for Coriolis_Scheme.

Definition at line 88 of file MOM_CoriolisAdv.F90.

88 character*(20), parameter :: ROBUST_ENSTRO_STRING = "ROBUST_ENSTRO"

Referenced by coriolisadv_init().

◆ sadourny75_energy

integer, parameter mom_coriolisadv::sadourny75_energy = 1

Enumeration values for Coriolis_Scheme.

Definition at line 80 of file MOM_CoriolisAdv.F90.

80 integer, parameter :: SADOURNY75_ENERGY = 1

Referenced by coradcalc(), and coriolisadv_init().

◆ sadourny75_energy_string

character*(20), parameter mom_coriolisadv::sadourny75_energy_string = "SADOURNY75_ENERGY"
private

Enumeration values for Coriolis_Scheme.

Definition at line 86 of file MOM_CoriolisAdv.F90.

86 character*(20), parameter :: SADOURNY75_ENERGY_STRING = "SADOURNY75_ENERGY"

Referenced by coriolisadv_init().

◆ sadourny75_enstro

integer, parameter mom_coriolisadv::sadourny75_enstro = 4
private

Enumeration values for Coriolis_Scheme.

Definition at line 83 of file MOM_CoriolisAdv.F90.

83 integer, parameter :: SADOURNY75_ENSTRO = 4

Referenced by coradcalc(), and coriolisadv_init().

◆ sadourny75_enstro_string

character*(20), parameter mom_coriolisadv::sadourny75_enstro_string = "SADOURNY75_ENSTRO"
private

Enumeration values for Coriolis_Scheme.

Definition at line 89 of file MOM_CoriolisAdv.F90.

89 character*(20), parameter :: SADOURNY75_ENSTRO_STRING = "SADOURNY75_ENSTRO"

Referenced by coriolisadv_init().