FV3DYCORE  Version 2.0.0
sw_core.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the FV3 dynamical core.
5 !*
6 !* The FV3 dynamical core is free software: you can redistribute it
7 !* and/or modify it under the terms of the
8 !* GNU Lesser General Public License as published by the
9 !* Free Software Foundation, either version 3 of the License, or
10 !* (at your option) any later version.
11 !*
12 !* The FV3 dynamical core is distributed in the hope that it will be
13 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
14 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15 !* See the GNU General Public License for more details.
16 !*
17 !* You should have received a copy of the GNU Lesser General Public
18 !* License along with the FV3 dynamical core.
19 !* If not, see <http://www.gnu.org/licenses/>.
20 !***********************************************************************
21 
25 
26  module sw_core_mod
27 
28 ! Modules Included:
29 ! <table>
30 ! <tr>
31 ! <th>Module Name</th>
32 ! <th>Functions Included</th>
33 ! </tr>
34 ! <tr>
35 ! <td>a2b_edge_mod</td>
36 ! <td>a2b_ord4</td>
37 ! <tr>
38 ! <td>fv_arrays_mod</td>
39 ! <td>ffv_grid_type, fv_grid_bounds_type, fv_flags_type</td>
40 ! </tr>
41 ! <tr>
42 ! <td>fv_mp_mod</td>
43 ! <td>ng,fill_corners, XDir, YDir</td>
44 ! </tr>
45 ! </table>
46 
47  use fv_mp_mod, only: ng
48  use tp_core_mod, only: fv_tp_2d, pert_ppm, copy_corners
49  use fv_mp_mod, only: fill_corners, xdir, ydir
50  use fv_arrays_mod, only: fv_grid_type, fv_grid_bounds_type, fv_flags_type
51  use a2b_edge_mod, only: a2b_ord4
52 
53 #ifdef SW_DYNAMICS
54  use test_cases_mod, only: test_case
55 #endif
56 
57  implicit none
58 
59  real, parameter:: r3 = 1./3.
60  real, parameter:: t11=27./28., t12=-13./28., t13=3./7., t14=6./7., t15=3./28.
61  real, parameter:: s11=11./14., s13=-13./14., s14=4./7., s15=3./14.
62  real, parameter:: near_zero = 1.e-9
63 #ifdef OVERLOAD_R4
64  real, parameter:: big_number = 1.e8
65 #else
66  real, parameter:: big_number = 1.e30
67 #endif
68 !----------------------
69 ! PPM volume mean form:
70 !----------------------
71  real, parameter:: p1 = 7./12.
72  real, parameter:: p2 = -1./12.
73 !----------------------------
74 ! 4-pt Lagrange interpolation
75 !----------------------------
76  real, parameter:: a1 = 0.5625
77  real, parameter:: a2 = -0.0625
78 !----------------------------------------------
79 ! volume-conserving cubic with 2nd drv=0 at end point:
80  real, parameter:: c1 = -2./14.
81  real, parameter:: c2 = 11./14.
82  real, parameter:: c3 = 5./14.
83 ! 3-pt off-center intp formular:
84 ! real, parameter:: c1 = -0.125
85 ! real, parameter:: c2 = 0.75
86 ! real, parameter:: c3 = 0.375
87 !----------------------------------------------
88 ! scheme 2.1: perturbation form
89  real, parameter:: b1 = 1./30.
90  real, parameter:: b2 = -13./60.
91  real, parameter:: b3 = -13./60.
92  real, parameter:: b4 = 0.45
93  real, parameter:: b5 = -0.05
94 
95 
96  private
97  public :: c_sw, d_sw, fill_4corners, del6_vt_flux, divergence_corner, divergence_corner_nest
98 
99  contains
100 
102  subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
103  ut, vt, divg_d, nord, dt2, hydrostatic, dord4, &
104  bd, gridstruct, flagstruct)
105 
106  type(fv_grid_bounds_type), intent(IN) :: bd
107  real, intent(INOUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed+1) :: u, vc
108  real, intent(INOUT), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ) :: v, uc
109  real, intent(INOUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed ) :: delp, pt, ua, va, ut, vt
110  real, intent(INOUT), dimension(bd%isd: , bd%jsd: ) :: w
111  real, intent(OUT ), dimension(bd%isd:bd%ied, bd%jsd:bd%jed ) :: delpc, ptc, wc
112  real, intent(OUT ), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed+1) :: divg_d
113  integer, intent(IN) :: nord
114  real, intent(IN) :: dt2
115  logical, intent(IN) :: hydrostatic
116  logical, intent(IN) :: dord4
117  type(fv_grid_type), intent(IN), target :: gridstruct
118  type(fv_flags_type), intent(IN), target :: flagstruct
119 
120 ! Local:
121  logical:: sw_corner, se_corner, ne_corner, nw_corner
122  real, dimension(bd%is-1:bd%ie+1,bd%js-1:bd%je+1):: vort, ke
123  real, dimension(bd%is-1:bd%ie+2,bd%js-1:bd%je+1):: fx, fx1, fx2
124  real, dimension(bd%is-1:bd%ie+1,bd%js-1:bd%je+2):: fy, fy1, fy2
125  real :: dt4
126  integer :: i,j
127  integer iep1, jep1
128 
129  integer :: is, ie, js, je
130  integer :: isd, ied, jsd, jed
131  integer :: npx, npy
132  logical :: bounded_domain
133 
134  real, pointer, dimension(:,:,:) :: sin_sg, cos_sg
135  real, pointer, dimension(:,:) :: cosa_u, cosa_v
136  real, pointer, dimension(:,:) :: sina_u, sina_v
137 
138  real, pointer, dimension(:,:) :: dx, dy, dxc, dyc
139 
140  is = bd%is
141  ie = bd%ie
142  js = bd%js
143  je = bd%je
144  isd = bd%isd
145  ied = bd%ied
146  jsd = bd%jsd
147  jed = bd%jed
148 
149  npx = flagstruct%npx
150  npy = flagstruct%npy
151  bounded_domain = gridstruct%bounded_domain
152 
153  sin_sg => gridstruct%sin_sg
154  cos_sg => gridstruct%cos_sg
155  cosa_u => gridstruct%cosa_u
156  cosa_v => gridstruct%cosa_v
157  sina_u => gridstruct%sina_u
158  sina_v => gridstruct%sina_v
159  dx => gridstruct%dx
160  dy => gridstruct%dy
161  dxc => gridstruct%dxc
162  dyc => gridstruct%dyc
163 
164  sw_corner = gridstruct%sw_corner
165  se_corner = gridstruct%se_corner
166  nw_corner = gridstruct%nw_corner
167  ne_corner = gridstruct%ne_corner
168 
169  iep1 = ie+1; jep1 = je+1
170 
171  call d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, bd, &
172  npx, npy, bounded_domain, flagstruct%grid_type)
173 
174  if( nord > 0 ) then
175  if (bounded_domain) then
176  call divergence_corner_nest(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
177  else
178  call divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
179  endif
180  endif
181 
182  do j=js-1,jep1
183  do i=is-1,iep1+1
184  if (ut(i,j) > 0.) then
185  ut(i,j) = dt2*ut(i,j)*dy(i,j)*sin_sg(i-1,j,3)
186  else
187  ut(i,j) = dt2*ut(i,j)*dy(i,j)*sin_sg(i,j,1)
188  end if
189  enddo
190  enddo
191  do j=js-1,je+2
192  do i=is-1,iep1
193  if (vt(i,j) > 0.) then
194  vt(i,j) = dt2*vt(i,j)*dx(i,j)*sin_sg(i,j-1,4)
195  else
196  vt(i,j) = dt2*vt(i,j)*dx(i,j)*sin_sg(i,j, 2)
197  end if
198  enddo
199  enddo
200 
201 !----------------
202 ! Transport delp:
203 !----------------
204 ! Xdir:
205  if (flagstruct%grid_type < 3 .and. .not. bounded_domain) call fill2_4corners(delp, pt, 1, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
206 
207  if ( hydrostatic ) then
208 #ifdef SW_DYNAMICS
209  do j=js-1,jep1
210  do i=is-1,ie+2
211  if ( ut(i,j) > 0. ) then
212  fx1(i,j) = delp(i-1,j)
213  else
214  fx1(i,j) = delp(i,j)
215  endif
216  fx1(i,j) = ut(i,j)*fx1(i,j)
217  enddo
218  enddo
219 #else
220  do j=js-1,jep1
221  do i=is-1,ie+2
222  if ( ut(i,j) > 0. ) then
223  fx1(i,j) = delp(i-1,j)
224  fx(i,j) = pt(i-1,j)
225  else
226  fx1(i,j) = delp(i,j)
227  fx(i,j) = pt(i,j)
228  endif
229  fx1(i,j) = ut(i,j)*fx1(i,j)
230  fx(i,j) = fx1(i,j)* fx(i,j)
231  enddo
232  enddo
233 #endif
234  else
235  if (flagstruct%grid_type < 3) &
236  call fill_4corners(w, 1, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
237  do j=js-1,je+1
238  do i=is-1,ie+2
239  if ( ut(i,j) > 0. ) then
240  fx1(i,j) = delp(i-1,j)
241  fx(i,j) = pt(i-1,j)
242  fx2(i,j) = w(i-1,j)
243  else
244  fx1(i,j) = delp(i,j)
245  fx(i,j) = pt(i,j)
246  fx2(i,j) = w(i,j)
247  endif
248  fx1(i,j) = ut(i,j)*fx1(i,j)
249  fx(i,j) = fx1(i,j)* fx(i,j)
250  fx2(i,j) = fx1(i,j)*fx2(i,j)
251  enddo
252  enddo
253  endif
254 
255 ! Ydir:
256  if (flagstruct%grid_type < 3 .and. .not. bounded_domain) call fill2_4corners(delp, pt, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
257  if ( hydrostatic ) then
258  do j=js-1,jep1+1
259  do i=is-1,iep1
260  if ( vt(i,j) > 0. ) then
261  fy1(i,j) = delp(i,j-1)
262  fy(i,j) = pt(i,j-1)
263  else
264  fy1(i,j) = delp(i,j)
265  fy(i,j) = pt(i,j)
266  endif
267  fy1(i,j) = vt(i,j)*fy1(i,j)
268  fy(i,j) = fy1(i,j)* fy(i,j)
269  enddo
270  enddo
271  do j=js-1,jep1
272  do i=is-1,iep1
273  delpc(i,j) = delp(i,j) + (fx1(i,j)-fx1(i+1,j)+fy1(i,j)-fy1(i,j+1))*gridstruct%rarea(i,j)
274 #ifdef SW_DYNAMICS
275  ptc(i,j) = pt(i,j)
276 #else
277  ptc(i,j) = (pt(i,j)*delp(i,j) + &
278  (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*gridstruct%rarea(i,j))/delpc(i,j)
279 #endif
280  enddo
281  enddo
282  else
283  if (flagstruct%grid_type < 3) call fill_4corners(w, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
284  do j=js-1,je+2
285  do i=is-1,ie+1
286  if ( vt(i,j) > 0. ) then
287  fy1(i,j) = delp(i,j-1)
288  fy(i,j) = pt(i,j-1)
289  fy2(i,j) = w(i,j-1)
290  else
291  fy1(i,j) = delp(i,j)
292  fy(i,j) = pt(i,j)
293  fy2(i,j) = w(i,j)
294  endif
295  fy1(i,j) = vt(i,j)*fy1(i,j)
296  fy(i,j) = fy1(i,j)* fy(i,j)
297  fy2(i,j) = fy1(i,j)*fy2(i,j)
298  enddo
299  enddo
300  do j=js-1,je+1
301  do i=is-1,ie+1
302  delpc(i,j) = delp(i,j) + (fx1(i,j)-fx1(i+1,j)+fy1(i,j)-fy1(i,j+1))*gridstruct%rarea(i,j)
303  ptc(i,j) = (pt(i,j)*delp(i,j) + &
304  (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*gridstruct%rarea(i,j))/delpc(i,j)
305  wc(i,j) = (w(i,j)*delp(i,j) + (fx2(i,j)-fx2(i+1,j) + &
306  fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j))/delpc(i,j)
307  enddo
308  enddo
309  endif
310 
311 !------------
312 ! Compute KE:
313 !------------
314 
315 !Since uc = u*, i.e. the covariant wind perpendicular to the face edge, if we want to compute kinetic energy we will need the true coordinate-parallel covariant wind, computed through u = uc*sina + v*cosa.
316 !Use the alpha for the cell KE is being computed in.
317 !!! TO DO:
318 !!! Need separate versions for nesting/single-tile
319 !!! and for cubed-sphere
320  if (bounded_domain .or. flagstruct%grid_type >=3 ) then
321  do j=js-1,jep1
322  do i=is-1,iep1
323  if ( ua(i,j) > 0. ) then
324  ke(i,j) = uc(i,j)
325  else
326  ke(i,j) = uc(i+1,j)
327  endif
328  enddo
329  enddo
330  do j=js-1,jep1
331  do i=is-1,iep1
332  if ( va(i,j) > 0. ) then
333  vort(i,j) = vc(i,j)
334  else
335  vort(i,j) = vc(i,j+1)
336  endif
337  enddo
338  enddo
339  else
340  do j=js-1,jep1
341  do i=is-1,iep1
342  if ( ua(i,j) > 0. ) then
343  if ( i==1 ) then
344  ke(1,j) = uc(1, j)*sin_sg(1,j,1)+v(1,j)*cos_sg(1,j,1)
345  elseif ( i==npx ) then
346  ke(i,j) = uc(npx,j)*sin_sg(npx,j,1)+v(npx,j)*cos_sg(npx,j,1)
347  else
348  ke(i,j) = uc(i,j)
349  endif
350  else
351  if ( i==0 ) then
352  ke(0,j) = uc(1, j)*sin_sg(0,j,3)+v(1,j)*cos_sg(0,j,3)
353  elseif ( i==(npx-1) ) then
354  ke(i,j) = uc(npx,j)*sin_sg(npx-1,j,3)+v(npx,j)*cos_sg(npx-1,j,3)
355  else
356  ke(i,j) = uc(i+1,j)
357  endif
358  endif
359  enddo
360  enddo
361  do j=js-1,jep1
362  do i=is-1,iep1
363  if ( va(i,j) > 0. ) then
364  if ( j==1 ) then
365  vort(i,1) = vc(i, 1)*sin_sg(i,1,2)+u(i, 1)*cos_sg(i,1,2)
366  elseif ( j==npy ) then
367  vort(i,j) = vc(i,npy)*sin_sg(i,npy,2)+u(i,npy)*cos_sg(i,npy,2)
368  else
369  vort(i,j) = vc(i,j)
370  endif
371  else
372  if ( j==0 ) then
373  vort(i,0) = vc(i, 1)*sin_sg(i,0,4)+u(i, 1)*cos_sg(i,0,4)
374  elseif ( j==(npy-1) ) then
375  vort(i,j) = vc(i,npy)*sin_sg(i,npy-1,4)+u(i,npy)*cos_sg(i,npy-1,4)
376  else
377  vort(i,j) = vc(i,j+1)
378  endif
379  endif
380  enddo
381  enddo
382  endif
383 
384  dt4 = 0.5*dt2
385  do j=js-1,jep1
386  do i=is-1,iep1
387  ke(i,j) = dt4*(ua(i,j)*ke(i,j) + va(i,j)*vort(i,j))
388  enddo
389  enddo
390 
391 !------------------------------
392 ! Compute circulation on C grid
393 !------------------------------
394 ! To consider using true co-variant winds at face edges?
395  do j=js-1,je+1
396  do i=is,ie+1
397  fx(i,j) = uc(i,j) * dxc(i,j)
398  enddo
399  enddo
400 
401  do j=js,je+1
402  do i=is-1,ie+1
403  fy(i,j) = vc(i,j) * dyc(i,j)
404  enddo
405  enddo
406 
407  do j=js,je+1
408  do i=is,ie+1
409  vort(i,j) = fx(i,j-1) - fx(i,j) - fy(i-1,j) + fy(i,j)
410  enddo
411  enddo
412 
413 ! Remove the extra term at the corners:
414  if ( sw_corner ) vort(1, 1) = vort(1, 1) + fy(0, 1)
415  if ( se_corner ) vort(npx ,1) = vort(npx, 1) - fy(npx, 1)
416  if ( ne_corner ) vort(npx,npy) = vort(npx,npy) - fy(npx,npy)
417  if ( nw_corner ) vort(1, npy) = vort(1, npy) + fy(0, npy)
418 
419 !----------------------------
420 ! Compute absolute vorticity
421 !----------------------------
422  do j=js,je+1
423  do i=is,ie+1
424  vort(i,j) = gridstruct%fC(i,j) + gridstruct%rarea_c(i,j) * vort(i,j)
425  enddo
426  enddo
427 
428 !----------------------------------
429 ! Transport absolute vorticity:
430 !----------------------------------
431 !To go from v to contravariant v at the edges, we divide by sin_sg;
432 ! but we then must multiply by sin_sg to get the proper flux.
433 ! These cancel, leaving us with fy1 = dt2*v at the edges.
434 ! (For the same reason we only divide by sin instead of sin**2 in the interior)
435 
436 !! TO DO: separate versions for nesting/single-tile and cubed-sphere
437  if (bounded_domain .or. flagstruct%grid_type >= 3) then
438  do j=js,je
439  do i=is,iep1
440  fy1(i,j) = dt2*(v(i,j)-uc(i,j)*cosa_u(i,j))/sina_u(i,j)
441  if ( fy1(i,j) > 0. ) then
442  fy(i,j) = vort(i,j)
443  else
444  fy(i,j) = vort(i,j+1)
445  endif
446  enddo
447  enddo
448  do j=js,jep1
449  do i=is,ie
450  fx1(i,j) = dt2*(u(i,j)-vc(i,j)*cosa_v(i,j))/sina_v(i,j)
451  if ( fx1(i,j) > 0. ) then
452  fx(i,j) = vort(i,j)
453  else
454  fx(i,j) = vort(i+1,j)
455  endif
456  enddo
457  enddo
458  else
459  do j=js,je
460 !DEC$ VECTOR ALWAYS
461  do i=is,iep1
462  if ( ( i==1 .or. i==npx ) ) then
463  fy1(i,j) = dt2*v(i,j)
464  else
465  fy1(i,j) = dt2*(v(i,j)-uc(i,j)*cosa_u(i,j))/sina_u(i,j)
466  endif
467  if ( fy1(i,j) > 0. ) then
468  fy(i,j) = vort(i,j)
469  else
470  fy(i,j) = vort(i,j+1)
471  endif
472  enddo
473  enddo
474  do j=js,jep1
475  if ( ( j==1 .or. j==npy ) ) then
476 !DEC$ VECTOR ALWAYS
477  do i=is,ie
478  fx1(i,j) = dt2*u(i,j)
479  if ( fx1(i,j) > 0. ) then
480  fx(i,j) = vort(i,j)
481  else
482  fx(i,j) = vort(i+1,j)
483  endif
484  enddo
485  else
486 !DEC$ VECTOR ALWAYS
487  do i=is,ie
488  fx1(i,j) = dt2*(u(i,j)-vc(i,j)*cosa_v(i,j))/sina_v(i,j)
489  if ( fx1(i,j) > 0. ) then
490  fx(i,j) = vort(i,j)
491  else
492  fx(i,j) = vort(i+1,j)
493  endif
494  enddo
495  endif
496  enddo
497  endif
498 
499 ! Update time-centered winds on the C-Grid
500  do j=js,je
501  do i=is,iep1
502  uc(i,j) = uc(i,j) + fy1(i,j)*fy(i,j) + gridstruct%rdxc(i,j)*(ke(i-1,j)-ke(i,j))
503  enddo
504  enddo
505  do j=js,jep1
506  do i=is,ie
507  vc(i,j) = vc(i,j) - fx1(i,j)*fx(i,j) + gridstruct%rdyc(i,j)*(ke(i,j-1)-ke(i,j))
508  enddo
509  enddo
510 
511  end subroutine c_sw
512 
513 
514 
515 ! d_sw :: D-Grid Shallow Water Routine
516 
519  subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
520  ua, va, divg_d, xflux, yflux, cx, cy, &
521  crx_adv, cry_adv, xfx_adv, yfx_adv, q_con, z_rat, kgb, heat_source,diss_est, &
522  zvir, sphum, nq, q, k, km, inline_q, &
523  dt, hord_tr, hord_mt, hord_vt, hord_tm, hord_dp, nord, &
524  nord_v, nord_w, nord_t, dddmp, d2_bg, d4_bg, damp_v, damp_w, &
525  damp_t, d_con, hydrostatic, gridstruct, flagstruct, bd)
526 
527  integer, intent(IN):: hord_tr, hord_mt, hord_vt, hord_tm, hord_dp
528  integer, intent(IN):: nord
529  integer, intent(IN):: nord_v
530  integer, intent(IN):: nord_w
531  integer, intent(IN):: nord_t
532  integer, intent(IN):: sphum, nq, k, km
533  real , intent(IN):: dt, dddmp, d2_bg, d4_bg, d_con
534  real , intent(IN):: zvir
535  real, intent(in):: damp_v, damp_w, damp_t, kgb
536  type(fv_grid_bounds_type), intent(IN) :: bd
537  real, intent(inout):: divg_d(bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
538  real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: z_rat
539  real, intent(INOUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: delp, pt, ua, va
540  real, intent(INOUT), dimension(bd%isd: , bd%jsd: ):: w, q_con
541  real, intent(INOUT), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed+1):: u, vc
542  real, intent(INOUT), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ):: v, uc
543  real, intent(INOUT):: q(bd%isd:bd%ied,bd%jsd:bd%jed,km,nq)
544  real, intent(OUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed) :: delpc, ptc
545  real, intent(OUT), dimension(bd%is:bd%ie,bd%js:bd%je):: heat_source
546  real, intent(OUT), dimension(bd%is:bd%ie,bd%js:bd%je):: diss_est
547 ! The flux capacitors:
548  real, intent(INOUT):: xflux(bd%is:bd%ie+1,bd%js:bd%je )
549  real, intent(INOUT):: yflux(bd%is:bd%ie ,bd%js:bd%je+1)
550 !------------------------
551  real, intent(INOUT):: cx(bd%is:bd%ie+1,bd%jsd:bd%jed )
552  real, intent(INOUT):: cy(bd%isd:bd%ied,bd%js:bd%je+1)
553  logical, intent(IN):: hydrostatic
554  logical, intent(IN):: inline_q
555  real, intent(OUT), dimension(bd%is:bd%ie+1,bd%jsd:bd%jed):: crx_adv, xfx_adv
556  real, intent(OUT), dimension(bd%isd:bd%ied,bd%js:bd%je+1):: cry_adv, yfx_adv
557  type(fv_grid_type), intent(IN), target :: gridstruct
558  type(fv_flags_type), intent(IN), target :: flagstruct
559 ! Local:
560  logical:: sw_corner, se_corner, ne_corner, nw_corner
561  real :: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed)
562  real :: vt(bd%isd:bd%ied, bd%jsd:bd%jed+1)
563 !---
564  real :: fx2(bd%isd:bd%ied+1,bd%jsd:bd%jed)
565  real :: fy2(bd%isd:bd%ied, bd%jsd:bd%jed+1)
566  real :: dw(bd%is:bd%ie,bd%js:bd%je)
567 !---
568  real, dimension(bd%is:bd%ie+1,bd%js:bd%je+1):: ub, vb
569  real :: wk(bd%isd:bd%ied,bd%jsd:bd%jed)
570  real :: ke(bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
571  real :: vort(bd%isd:bd%ied,bd%jsd:bd%jed)
572  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
573  real :: fy(bd%is:bd%ie ,bd%js:bd%je+1)
574  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
575  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
576  real :: gx(bd%is:bd%ie+1,bd%js:bd%je )
577  real :: gy(bd%is:bd%ie ,bd%js:bd%je+1)
578  logical :: fill_c
579 
580  real :: dt2, dt4, dt5, dt6
581  real :: damp, damp2, damp4, dd8, u2, v2, du2, dv2
582  real :: u_lon
583  integer :: i,j, is2, ie1, js2, je1, n, nt, n2, iq
584 
585  real, pointer, dimension(:,:) :: area, area_c, rarea
586 
587  real, pointer, dimension(:,:,:) :: sin_sg
588  real, pointer, dimension(:,:) :: cosa_u, cosa_v, cosa_s
589  real, pointer, dimension(:,:) :: sina_u, sina_v
590  real, pointer, dimension(:,:) :: rsin_u, rsin_v, rsina
591  real, pointer, dimension(:,:) :: f0, rsin2, divg_u, divg_v
592 
593  real, pointer, dimension(:,:) :: cosa, dx, dy, dxc, dyc, rdxa, rdya, rdx, rdy
594 
595  integer :: is, ie, js, je
596  integer :: isd, ied, jsd, jed
597  integer :: npx, npy, ng
598  logical :: bounded_domain
599 
600  is = bd%is
601  ie = bd%ie
602  js = bd%js
603  je = bd%je
604  isd = bd%isd
605  ied = bd%ied
606  jsd = bd%jsd
607  jed = bd%jed
608  ng = bd%ng
609 
610  npx = flagstruct%npx
611  npy = flagstruct%npy
612  bounded_domain = gridstruct%bounded_domain
613 
614  area => gridstruct%area
615  rarea => gridstruct%rarea
616  sin_sg => gridstruct%sin_sg
617  cosa_u => gridstruct%cosa_u
618  cosa_v => gridstruct%cosa_v
619  cosa_s => gridstruct%cosa_s
620  sina_u => gridstruct%sina_u
621  sina_v => gridstruct%sina_v
622  rsin_u => gridstruct%rsin_u
623  rsin_v => gridstruct%rsin_v
624  rsina => gridstruct%rsina
625  f0 => gridstruct%f0
626  rsin2 => gridstruct%rsin2
627  divg_u => gridstruct%divg_u
628  divg_v => gridstruct%divg_v
629  cosa => gridstruct%cosa
630  dx => gridstruct%dx
631  dy => gridstruct%dy
632  dxc => gridstruct%dxc
633  dyc => gridstruct%dyc
634  rdxa => gridstruct%rdxa
635  rdya => gridstruct%rdya
636  rdx => gridstruct%rdx
637  rdy => gridstruct%rdy
638 
639  sw_corner = gridstruct%sw_corner
640  se_corner = gridstruct%se_corner
641  nw_corner = gridstruct%nw_corner
642  ne_corner = gridstruct%ne_corner
643 
644 #ifdef SW_DYNAMICS
645  if ( test_case == 1 ) then
646  do j=jsd,jed
647  do i=is,ie+1
648  xfx_adv(i,j) = dt * uc(i,j) / sina_u(i,j)
649  if (xfx_adv(i,j) > 0.) then
650  crx_adv(i,j) = xfx_adv(i,j) * rdxa(i-1,j)
651  else
652  crx_adv(i,j) = xfx_adv(i,j) * rdxa(i,j)
653  endif
654  xfx_adv(i,j) = dy(i,j)*xfx_adv(i,j)*sina_u(i,j)
655  enddo
656  enddo
657 
658  do j=js,je+1
659  do i=isd,ied
660  yfx_adv(i,j) = dt * vc(i,j) / sina_v(i,j)
661  if (yfx_adv(i,j) > 0.) then
662  cry_adv(i,j) = yfx_adv(i,j) * rdya(i,j-1)
663  else
664  cry_adv(i,j) = yfx_adv(i,j) * rdya(i,j)
665  endif
666  yfx_adv(i,j) = dx(i,j)*yfx_adv(i,j)*sina_v(i,j)
667  enddo
668  enddo
669  else
670 #endif
671  if ( flagstruct%grid_type < 3 ) then
672 
673 !!! TO DO: separate versions for nesting and for cubed-sphere
674  if (bounded_domain) then
675  do j=jsd,jed
676  do i=is-1,ie+2
677  ut(i,j) = ( uc(i,j) - 0.25 * cosa_u(i,j) * &
678  (vc(i-1,j)+vc(i,j)+vc(i-1,j+1)+vc(i,j+1)))*rsin_u(i,j)
679  enddo
680  enddo
681  do j=js-1,je+2
682  do i=isd,ied
683  vt(i,j) = ( vc(i,j) - 0.25 * cosa_v(i,j) * &
684  (uc(i,j-1)+uc(i+1,j-1)+uc(i,j)+uc(i+1,j)))*rsin_v(i,j)
685  enddo
686  enddo
687  else
688  do j=jsd,jed
689  if( j/=0 .and. j/=1 .and. j/=(npy-1) .and. j/=npy) then
690  do i=is-1,ie+2
691  ut(i,j) = ( uc(i,j) - 0.25 * cosa_u(i,j) * &
692  (vc(i-1,j)+vc(i,j)+vc(i-1,j+1)+vc(i,j+1)))*rsin_u(i,j)
693  enddo
694  endif
695  enddo
696  do j=js-1,je+2
697 
698  if( j/=1 .and. j/=npy ) then
699 
700  do i=isd,ied
701  vt(i,j) = ( vc(i,j) - 0.25 * cosa_v(i,j) * &
702  (uc(i,j-1)+uc(i+1,j-1)+uc(i,j)+uc(i+1,j)))*rsin_v(i,j)
703  enddo
704  endif
705  enddo
706  endif
707 
708  if (.not. bounded_domain) then
709 ! West edge:
710  if ( is==1 ) then
711  do j=jsd,jed
712  if ( uc(1,j)*dt > 0. ) then
713  ut(1,j) = uc(1,j) / sin_sg(0,j,3)
714  else
715  ut(1,j) = uc(1,j) / sin_sg(1,j,1)
716  endif
717  enddo
718  do j=max(3,js), min(npy-2,je+1)
719  vt(0,j) = vc(0,j) - 0.25*cosa_v(0,j)* &
720  (ut(0,j-1)+ut(1,j-1)+ut(0,j)+ut(1,j))
721  vt(1,j) = vc(1,j) - 0.25*cosa_v(1,j)* &
722  (ut(1,j-1)+ut(2,j-1)+ut(1,j)+ut(2,j))
723  enddo
724  endif ! West face
725 
726 ! East edge:
727  if ( (ie+1)==npx ) then
728  do j=jsd,jed
729  if ( uc(npx,j)*dt > 0. ) then
730  ut(npx,j) = uc(npx,j) / sin_sg(npx-1,j,3)
731  else
732  ut(npx,j) = uc(npx,j) / sin_sg(npx,j,1)
733  endif
734  enddo
735 
736  do j=max(3,js), min(npy-2,je+1)
737  vt(npx-1,j) = vc(npx-1,j) - 0.25*cosa_v(npx-1,j)* &
738  (ut(npx-1,j-1)+ut(npx,j-1)+ut(npx-1,j)+ut(npx,j))
739  vt(npx,j) = vc(npx,j) - 0.25*cosa_v(npx,j)* &
740  (ut(npx,j-1)+ut(npx+1,j-1)+ut(npx,j)+ut(npx+1,j))
741  enddo
742  endif
743 
744 ! South (Bottom) edge:
745  if ( js==1 ) then
746 
747  do i=isd,ied
748  if ( vc(i,1)*dt > 0. ) then
749  vt(i,1) = vc(i,1) / sin_sg(i,0,4)
750  else
751  vt(i,1) = vc(i,1) / sin_sg(i,1,2)
752  endif
753  enddo
754 
755  do i=max(3,is),min(npx-2,ie+1)
756  ut(i,0) = uc(i,0) - 0.25*cosa_u(i,0)* &
757  (vt(i-1,0)+vt(i,0)+vt(i-1,1)+vt(i,1))
758  ut(i,1) = uc(i,1) - 0.25*cosa_u(i,1)* &
759  (vt(i-1,1)+vt(i,1)+vt(i-1,2)+vt(i,2))
760  enddo
761  endif
762 
763 ! North edge:
764  if ( (je+1)==npy ) then
765  do i=isd,ied
766  if ( vc(i,npy)*dt > 0. ) then
767  vt(i,npy) = vc(i,npy) / sin_sg(i,npy-1,4)
768  else
769  vt(i,npy) = vc(i,npy) / sin_sg(i,npy,2)
770  endif
771  enddo
772  do i=max(3,is),min(npx-2,ie+1)
773  ut(i,npy-1) = uc(i,npy-1) - 0.25*cosa_u(i,npy-1)* &
774  (vt(i-1,npy-1)+vt(i,npy-1)+vt(i-1,npy)+vt(i,npy))
775  ut(i,npy) = uc(i,npy) - 0.25*cosa_u(i,npy)* &
776  (vt(i-1,npy)+vt(i,npy)+vt(i-1,npy+1)+vt(i,npy+1))
777  enddo
778  endif
779 
780 ! The following code solves a 2x2 system to get the interior parallel-to-edge uc,vc values
781 ! near the corners (ex: for the sw corner ut(2,1) and vt(1,2) are solved for simultaneously).
782 ! It then computes the halo uc, vc values so as to be consistent with the computations on
783 ! the facing panel.
784 
785  !The system solved is:
786  ! ut(2,1) = uc(2,1) - avg(vt)*cosa_u(2,1)
787  ! vt(1,2) = vc(1,2) - avg(ut)*cosa_v(1,2)
788  ! in which avg(vt) includes vt(1,2) and avg(ut) includes ut(2,1)
789 
790  if( sw_corner ) then
791  damp = 1. / (1.-0.0625*cosa_u(2,0)*cosa_v(1,0))
792  ut(2,0) = (uc(2,0)-0.25*cosa_u(2,0)*(vt(1,1)+vt(2,1)+vt(2,0) +vc(1,0) - &
793  0.25*cosa_v(1,0)*(ut(1,0)+ut(1,-1)+ut(2,-1))) ) * damp
794  damp = 1. / (1.-0.0625*cosa_u(0,1)*cosa_v(0,2))
795  vt(0,2) = (vc(0,2)-0.25*cosa_v(0,2)*(ut(1,1)+ut(1,2)+ut(0,2)+uc(0,1) - &
796  0.25*cosa_u(0,1)*(vt(0,1)+vt(-1,1)+vt(-1,2))) ) * damp
797 
798  damp = 1. / (1.-0.0625*cosa_u(2,1)*cosa_v(1,2))
799  ut(2,1) = (uc(2,1)-0.25*cosa_u(2,1)*(vt(1,1)+vt(2,1)+vt(2,2)+vc(1,2) - &
800  0.25*cosa_v(1,2)*(ut(1,1)+ut(1,2)+ut(2,2))) ) * damp
801 
802  vt(1,2) = (vc(1,2)-0.25*cosa_v(1,2)*(ut(1,1)+ut(1,2)+ut(2,2)+uc(2,1) - &
803  0.25*cosa_u(2,1)*(vt(1,1)+vt(2,1)+vt(2,2))) ) * damp
804  endif
805 
806  if( se_corner ) then
807  damp = 1. / (1. - 0.0625*cosa_u(npx-1,0)*cosa_v(npx-1,0))
808  ut(npx-1,0) = ( uc(npx-1,0)-0.25*cosa_u(npx-1,0)*( &
809  vt(npx-1,1)+vt(npx-2,1)+vt(npx-2,0)+vc(npx-1,0) - &
810  0.25*cosa_v(npx-1,0)*(ut(npx,0)+ut(npx,-1)+ut(npx-1,-1))) ) * damp
811  damp = 1. / (1. - 0.0625*cosa_u(npx+1,1)*cosa_v(npx,2))
812  vt(npx, 2) = ( vc(npx,2)-0.25*cosa_v(npx,2)*( &
813  ut(npx,1)+ut(npx,2)+ut(npx+1,2)+uc(npx+1,1) - &
814  0.25*cosa_u(npx+1,1)*(vt(npx,1)+vt(npx+1,1)+vt(npx+1,2))) ) * damp
815 
816  damp = 1. / (1. - 0.0625*cosa_u(npx-1,1)*cosa_v(npx-1,2))
817  ut(npx-1,1) = ( uc(npx-1,1)-0.25*cosa_u(npx-1,1)*( &
818  vt(npx-1,1)+vt(npx-2,1)+vt(npx-2,2)+vc(npx-1,2) - &
819  0.25*cosa_v(npx-1,2)*(ut(npx,1)+ut(npx,2)+ut(npx-1,2))) ) * damp
820  vt(npx-1,2) = ( vc(npx-1,2)-0.25*cosa_v(npx-1,2)*( &
821  ut(npx,1)+ut(npx,2)+ut(npx-1,2)+uc(npx-1,1) - &
822  0.25*cosa_u(npx-1,1)*(vt(npx-1,1)+vt(npx-2,1)+vt(npx-2,2))) ) * damp
823  endif
824 
825  if( ne_corner ) then
826  damp = 1. / (1. - 0.0625*cosa_u(npx-1,npy)*cosa_v(npx-1,npy+1))
827  ut(npx-1,npy) = ( uc(npx-1,npy)-0.25*cosa_u(npx-1,npy)*( &
828  vt(npx-1,npy)+vt(npx-2,npy)+vt(npx-2,npy+1)+vc(npx-1,npy+1) - &
829  0.25*cosa_v(npx-1,npy+1)*(ut(npx,npy)+ut(npx,npy+1)+ut(npx-1,npy+1))) ) * damp
830  damp = 1. / (1. - 0.0625*cosa_u(npx+1,npy-1)*cosa_v(npx,npy-1))
831  vt(npx, npy-1) = ( vc(npx,npy-1)-0.25*cosa_v(npx,npy-1)*( &
832  ut(npx,npy-1)+ut(npx,npy-2)+ut(npx+1,npy-2)+uc(npx+1,npy-1) - &
833  0.25*cosa_u(npx+1,npy-1)*(vt(npx,npy)+vt(npx+1,npy)+vt(npx+1,npy-1))) ) * damp
834 
835  damp = 1. / (1. - 0.0625*cosa_u(npx-1,npy-1)*cosa_v(npx-1,npy-1))
836  ut(npx-1,npy-1) = ( uc(npx-1,npy-1)-0.25*cosa_u(npx-1,npy-1)*( &
837  vt(npx-1,npy)+vt(npx-2,npy)+vt(npx-2,npy-1)+vc(npx-1,npy-1) - &
838  0.25*cosa_v(npx-1,npy-1)*(ut(npx,npy-1)+ut(npx,npy-2)+ut(npx-1,npy-2))) ) * damp
839  vt(npx-1,npy-1) = ( vc(npx-1,npy-1)-0.25*cosa_v(npx-1,npy-1)*( &
840  ut(npx,npy-1)+ut(npx,npy-2)+ut(npx-1,npy-2)+uc(npx-1,npy-1) - &
841  0.25*cosa_u(npx-1,npy-1)*(vt(npx-1,npy)+vt(npx-2,npy)+vt(npx-2,npy-1))) ) * damp
842  endif
843 
844  if( nw_corner ) then
845  damp = 1. / (1. - 0.0625*cosa_u(2,npy)*cosa_v(1,npy+1))
846  ut(2,npy) = ( uc(2,npy)-0.25*cosa_u(2,npy)*( &
847  vt(1,npy)+vt(2,npy)+vt(2,npy+1)+vc(1,npy+1) - &
848  0.25*cosa_v(1,npy+1)*(ut(1,npy)+ut(1,npy+1)+ut(2,npy+1))) ) * damp
849  damp = 1. / (1. - 0.0625*cosa_u(0,npy-1)*cosa_v(0,npy-1))
850  vt(0,npy-1) = ( vc(0,npy-1)-0.25*cosa_v(0,npy-1)*( &
851  ut(1,npy-1)+ut(1,npy-2)+ut(0,npy-2)+uc(0,npy-1) - &
852  0.25*cosa_u(0,npy-1)*(vt(0,npy)+vt(-1,npy)+vt(-1,npy-1))) ) * damp
853 
854  damp = 1. / (1. - 0.0625*cosa_u(2,npy-1)*cosa_v(1,npy-1))
855  ut(2,npy-1) = ( uc(2,npy-1)-0.25*cosa_u(2,npy-1)*( &
856  vt(1,npy)+vt(2,npy)+vt(2,npy-1)+vc(1,npy-1) - &
857  0.25*cosa_v(1,npy-1)*(ut(1,npy-1)+ut(1,npy-2)+ut(2,npy-2))) ) * damp
858 
859  vt(1,npy-1) = ( vc(1,npy-1)-0.25*cosa_v(1,npy-1)*( &
860  ut(1,npy-1)+ut(1,npy-2)+ut(2,npy-2)+uc(2,npy-1) - &
861  0.25*cosa_u(2,npy-1)*(vt(1,npy)+vt(2,npy)+vt(2,npy-1))) ) * damp
862  endif
863 
864  end if !.not. bounded_domain
865 
866  else
867 ! flagstruct%grid_type >= 3
868  do j=jsd,jed
869  do i=is,ie+1
870  ut(i,j) = uc(i,j)
871  enddo
872  enddo
873 
874  do j=js,je+1
875  do i=isd,ied
876  vt(i,j) = vc(i,j)
877  enddo
878  enddo
879  endif ! end grid_type choices
880 
881  do j=jsd,jed
882  do i=is,ie+1
883  xfx_adv(i,j) = dt*ut(i,j)
884  enddo
885  enddo
886 
887  do j=js,je+1
888  do i=isd,ied
889  yfx_adv(i,j) = dt*vt(i,j)
890  enddo
891  enddo
892 
893 ! Explanation of the following code:
894 ! xfx_adv = dt*ut*dy
895 ! crx_adv = dt*ut/dx
896 
897  do j=jsd,jed
898 !DEC$ VECTOR ALWAYS
899  do i=is,ie+1
900  if ( xfx_adv(i,j) > 0. ) then
901  crx_adv(i,j) = xfx_adv(i,j) * rdxa(i-1,j)
902  xfx_adv(i,j) = dy(i,j)*xfx_adv(i,j)*sin_sg(i-1,j,3)
903  else
904  crx_adv(i,j) = xfx_adv(i,j) * rdxa(i,j)
905  xfx_adv(i,j) = dy(i,j)*xfx_adv(i,j)*sin_sg(i,j,1)
906  end if
907  enddo
908  enddo
909  do j=js,je+1
910 !DEC$ VECTOR ALWAYS
911  do i=isd,ied
912  if ( yfx_adv(i,j) > 0. ) then
913  cry_adv(i,j) = yfx_adv(i,j) * rdya(i,j-1)
914  yfx_adv(i,j) = dx(i,j)*yfx_adv(i,j)*sin_sg(i,j-1,4)
915  else
916  cry_adv(i,j) = yfx_adv(i,j) * rdya(i,j)
917  yfx_adv(i,j) = dx(i,j)*yfx_adv(i,j)*sin_sg(i,j,2)
918  endif
919  enddo
920  enddo
921 
922 #ifdef SW_DYNAMICS
923  endif
924 #endif
925 
926  do j=jsd,jed
927  do i=is,ie
928  ra_x(i,j) = area(i,j) + xfx_adv(i,j) - xfx_adv(i+1,j)
929  enddo
930  enddo
931  do j=js,je
932  do i=isd,ied
933  ra_y(i,j) = area(i,j) + yfx_adv(i,j) - yfx_adv(i,j+1)
934  enddo
935  enddo
936 
937 
938  call fv_tp_2d(delp, crx_adv, cry_adv, npx, npy, hord_dp, fx, fy, &
939  xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac, &
940  nord=nord_v, damp_c=damp_v)
941 
942 ! <<< Save the mass fluxes to the "Flux Capacitor" for tracer transport >>>
943  do j=jsd,jed
944  do i=is,ie+1
945  cx(i,j) = cx(i,j) + crx_adv(i,j)
946  enddo
947  enddo
948  do j=js,je
949  do i=is,ie+1
950  xflux(i,j) = xflux(i,j) + fx(i,j)
951  enddo
952  enddo
953  do j=js,je+1
954  do i=isd,ied
955  cy(i,j) = cy(i,j) + cry_adv(i,j)
956  enddo
957  do i=is,ie
958  yflux(i,j) = yflux(i,j) + fy(i,j)
959  enddo
960  enddo
961 
962 #ifndef SW_DYNAMICS
963  do j=js,je
964  do i=is,ie
965  heat_source(i,j) = 0.
966  diss_est(i,j) = 0.
967  enddo
968  enddo
969 
970  if ( .not. hydrostatic ) then
971  if ( damp_w>1.e-5 ) then
972  dd8 = kgb*abs(dt)
973  damp4 = (damp_w*gridstruct%da_min_c)**(nord_w+1)
974  call del6_vt_flux(nord_w, npx, npy, damp4, w, wk, fx2, fy2, gridstruct, bd)
975  do j=js,je
976  do i=is,ie
977  dw(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*rarea(i,j)
978 ! 0.5 * [ (w+dw)**2 - w**2 ] = w*dw + 0.5*dw*dw
979 ! heat_source(i,j) = -d_con*dw(i,j)*(w(i,j)+0.5*dw(i,j))
980  heat_source(i,j) = dd8 - dw(i,j)*(w(i,j)+0.5*dw(i,j))
981  diss_est(i,j) = heat_source(i,j)
982  enddo
983  enddo
984  endif
985  call fv_tp_2d(w, crx_adv,cry_adv, npx, npy, hord_vt, gx, gy, xfx_adv, yfx_adv, &
986  gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac, &
987  mfx=fx, mfy=fy)
988  do j=js,je
989  do i=is,ie
990  w(i,j) = delp(i,j)*w(i,j) + (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j)
991  enddo
992  enddo
993  endif
994 
995 #ifdef USE_COND
996  call fv_tp_2d(q_con, crx_adv,cry_adv, npx, npy, hord_dp, gx, gy, &
997  xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac,&
998  mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t)
999  do j=js,je
1000  do i=is,ie
1001  q_con(i,j) = delp(i,j)*q_con(i,j) + (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j)
1002  enddo
1003  enddo
1004 #endif
1005 
1006 ! if ( inline_q .and. zvir>0.01 ) then
1007 ! do j=jsd,jed
1008 ! do i=isd,ied
1009 ! pt(i,j) = pt(i,j)/(1.+zvir*q(i,j,k,sphum))
1010 ! enddo
1011 ! enddo
1012 ! endif
1013  call fv_tp_2d(pt, crx_adv,cry_adv, npx, npy, hord_tm, gx, gy, &
1014  xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac, &
1015  mfx=fx, mfy=fy, mass=delp, nord=nord_v, damp_c=damp_v)
1016 ! mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t)
1017 #endif
1018 
1019  if ( inline_q ) then
1020  do j=js,je
1021  do i=is,ie
1022  wk(i,j) = delp(i,j)
1023  delp(i,j) = wk(i,j) + (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j)
1024 #ifdef SW_DYNAMICS
1025  ptc(i,j) = pt(i,j)
1026 #else
1027  pt(i,j) = (pt(i,j)*wk(i,j) + &
1028  (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j))/delp(i,j)
1029 #endif
1030  enddo
1031  enddo
1032  do iq=1,nq
1033  call fv_tp_2d(q(isd,jsd,k,iq), crx_adv,cry_adv, npx, npy, hord_tr, gx, gy, &
1034  xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac, &
1035  mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t)
1036  do j=js,je
1037  do i=is,ie
1038  q(i,j,k,iq) = (q(i,j,k,iq)*wk(i,j) + &
1039  (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j))/delp(i,j)
1040  enddo
1041  enddo
1042  enddo
1043 ! if ( zvir>0.01 ) then
1044 ! do j=js,je
1045 ! do i=is,ie
1046 ! pt(i,j) = pt(i,j)*(1.+zvir*q(i,j,k,sphum))
1047 ! enddo
1048 ! enddo
1049 ! endif
1050 
1051  else
1052  do j=js,je
1053  do i=is,ie
1054 #ifndef SW_DYNAMICS
1055  pt(i,j) = pt(i,j)*delp(i,j) + &
1056  (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j)
1057 #endif
1058  delp(i,j) = delp(i,j) + &
1059  (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j)
1060 #ifndef SW_DYNAMICS
1061  pt(i,j) = pt(i,j) / delp(i,j)
1062 
1063 #endif
1064  enddo
1065  enddo
1066  endif
1067 
1068 #ifdef SW_DYNAMICS
1069  if (test_case > 1) then
1070 #endif
1071 
1072 !----------------------
1073 ! Kinetic Energy Fluxes
1074 !----------------------
1075 ! Compute B grid contra-variant components for KE:
1076 
1077  dt5 = 0.5 *dt
1078  dt4 = 0.25*dt
1079 
1080  if (bounded_domain) then
1081  is2 = is; ie1 = ie+1
1082  js2 = js; je1 = je+1
1083  else
1084  is2 = max(2,is); ie1 = min(npx-1,ie+1)
1085  js2 = max(2,js); je1 = min(npy-1,je+1)
1086  end if
1087 
1088  if (flagstruct%grid_type < 3) then
1089 
1090  if (bounded_domain) then
1091  do j=js2,je1
1092  do i=is2,ie1
1093  vb(i,j) = dt5*(vc(i-1,j)+vc(i,j)-(uc(i,j-1)+uc(i,j))*cosa(i,j))*rsina(i,j)
1094  enddo
1095  enddo
1096  else
1097  if ( js==1 ) then
1098  do i=is,ie+1
1099  vb(i,1) = dt5*(vt(i-1,1)+vt(i,1)) ! corner values are incorrect
1100  enddo
1101  endif
1102 
1103  do j=js2,je1
1104  do i=is2,ie1
1105  vb(i,j) = dt5*(vc(i-1,j)+vc(i,j)-(uc(i,j-1)+uc(i,j))*cosa(i,j))*rsina(i,j)
1106  enddo
1107 
1108  if ( is==1 ) then
1109  ! 2-pt extrapolation from both sides:
1110  vb(1,j) = dt4*(-vt(-1,j) + 3.*(vt(0,j)+vt(1,j)) - vt(2,j))
1111  endif
1112  if ( (ie+1)==npx ) then
1113  ! 2-pt extrapolation from both sides:
1114  vb(npx,j) = dt4*(-vt(npx-2,j) + 3.*(vt(npx-1,j)+vt(npx,j)) - vt(npx+1,j))
1115  endif
1116  enddo
1117 
1118  if ( (je+1)==npy ) then
1119  do i=is,ie+1
1120  vb(i,npy) = dt5*(vt(i-1,npy)+vt(i,npy)) ! corner values are incorrect
1121  enddo
1122  endif
1123  endif
1124 
1125  else
1126  do j=js,je+1
1127  do i=is,ie+1
1128  vb(i,j) = dt5*(vc(i-1,j)+vc(i,j))
1129  enddo
1130  enddo
1131  endif
1132 
1133  call ytp_v(is,ie,js,je,isd,ied,jsd,jed, vb, u, v, ub, hord_mt, gridstruct%dy, gridstruct%rdy, &
1134  npx, npy, flagstruct%grid_type, flagstruct%lim_fac, bounded_domain)
1135 
1136  do j=js,je+1
1137  do i=is,ie+1
1138  ke(i,j) = vb(i,j)*ub(i,j)
1139  enddo
1140  enddo
1141 
1142  if (flagstruct%grid_type < 3) then
1143 
1144  if (bounded_domain) then
1145 
1146  do j=js,je+1
1147 
1148  do i=is2,ie1
1149  ub(i,j) = dt5*(uc(i,j-1)+uc(i,j)-(vc(i-1,j)+vc(i,j))*cosa(i,j))*rsina(i,j)
1150  enddo
1151 
1152  enddo
1153 
1154 
1155  else
1156  if ( is==1 ) then
1157  do j=js,je+1
1158  ub(1,j) = dt5*(ut(1,j-1)+ut(1,j)) ! corner values are incorrect
1159  enddo
1160  endif
1161 
1162  do j=js,je+1
1163  if ( (j==1 .or. j==npy) ) then
1164  do i=is2,ie1
1165  ! 2-pt extrapolation from both sides:
1166  ub(i,j) = dt4*(-ut(i,j-2) + 3.*(ut(i,j-1)+ut(i,j)) - ut(i,j+1))
1167  enddo
1168  else
1169  do i=is2,ie1
1170  ub(i,j) = dt5*(uc(i,j-1)+uc(i,j)-(vc(i-1,j)+vc(i,j))*cosa(i,j))*rsina(i,j)
1171  enddo
1172  endif
1173  enddo
1174 
1175  if ( (ie+1)==npx ) then
1176  do j=js,je+1
1177  ub(npx,j) = dt5*(ut(npx,j-1)+ut(npx,j)) ! corner values are incorrect
1178  enddo
1179  endif
1180  endif
1181 
1182  else
1183  do j=js,je+1
1184  do i=is,ie+1
1185  ub(i,j) = dt5*(uc(i,j-1)+uc(i,j))
1186  enddo
1187  enddo
1188  endif
1189 
1190  call xtp_u(is,ie,js,je, isd,ied,jsd,jed, ub, u, v, vb, hord_mt, gridstruct%dx, gridstruct%rdx, &
1191  npx, npy, flagstruct%grid_type, flagstruct%lim_fac, bounded_domain)
1192 
1193  do j=js,je+1
1194  do i=is,ie+1
1195  ke(i,j) = 0.5*(ke(i,j) + ub(i,j)*vb(i,j))
1196  enddo
1197  enddo
1198 
1199 !-----------------------------------------
1200 ! Fix KE at the 4 corners of the face:
1201 !-----------------------------------------
1202  if (.not. bounded_domain) then
1203  dt6 = dt / 6.
1204  if ( sw_corner ) then
1205  ke(1,1) = dt6*( (ut(1,1) + ut(1,0)) * u(1,1) + &
1206  (vt(1,1) + vt(0,1)) * v(1,1) + &
1207  (ut(1,1) + vt(1,1)) * u(0,1) )
1208  endif
1209  if ( se_corner ) then
1210  i = npx
1211  ke(i,1) = dt6*( (ut(i,1) + ut(i, 0)) * u(i-1,1) + &
1212  (vt(i,1) + vt(i-1,1)) * v(i, 1) + &
1213  (ut(i,1) - vt(i-1,1)) * u(i, 1) )
1214  endif
1215  if ( ne_corner ) then
1216  i = npx; j = npy
1217  ke(i,j) = dt6*( (ut(i,j ) + ut(i,j-1)) * u(i-1,j) + &
1218  (vt(i,j ) + vt(i-1,j)) * v(i,j-1) + &
1219  (ut(i,j-1) + vt(i-1,j)) * u(i,j ) )
1220  endif
1221  if ( nw_corner ) then
1222  j = npy
1223  ke(1,j) = dt6*( (ut(1, j) + ut(1,j-1)) * u(1,j ) + &
1224  (vt(1, j) + vt(0, j)) * v(1,j-1) + &
1225  (ut(1,j-1) - vt(1, j)) * u(0,j ) )
1226  endif
1227  end if
1228 
1229 ! Compute vorticity:
1230  do j=jsd,jed+1
1231  do i=isd,ied
1232  vt(i,j) = u(i,j)*dx(i,j)
1233  enddo
1234  enddo
1235  do j=jsd,jed
1236  do i=isd,ied+1
1237  ut(i,j) = v(i,j)*dy(i,j)
1238  enddo
1239  enddo
1240 
1241 ! wk is "volume-mean" relative vorticity
1242  do j=jsd,jed
1243  do i=isd,ied
1244  wk(i,j) = rarea(i,j)*(vt(i,j)-vt(i,j+1)-ut(i,j)+ut(i+1,j))
1245  enddo
1246  enddo
1247 
1248  if ( .not. hydrostatic ) then
1249  if( flagstruct%do_f3d ) then
1250 #ifdef ROT3
1251  dt2 = 2.*dt
1252  do j=js,je
1253  do i=is,ie
1254  w(i,j) = w(i,j)/delp(i,j) + dt2*gridstruct%w00(i,j) * &
1255  ( gridstruct%a11(i,j)*(u(i,j)+u(i,j+1)) + &
1256  gridstruct%a12(i,j)*(v(i,j)+v(i+1,j)) )
1257  enddo
1258  enddo
1259 #endif
1260  else
1261  do j=js,je
1262  do i=is,ie
1263  w(i,j) = w(i,j)/delp(i,j)
1264  enddo
1265  enddo
1266  endif
1267  if ( damp_w>1.e-5 ) then
1268  do j=js,je
1269  do i=is,ie
1270  w(i,j) = w(i,j) + dw(i,j)
1271  enddo
1272  enddo
1273  endif
1274 
1275  endif
1276 #ifdef USE_COND
1277  do j=js,je
1278  do i=is,ie
1279  q_con(i,j) = q_con(i,j)/delp(i,j)
1280  enddo
1281  enddo
1282 #endif
1283 
1284 !-----------------------------
1285 ! Compute divergence damping
1286 !-----------------------------
1287 ! damp = dddmp * da_min_c
1288 
1289  if ( nord==0 ) then
1290 ! area ~ dxb*dyb*sin(alpha)
1291 
1292  if (bounded_domain) then
1293 
1294  do j=js,je+1
1295  do i=is-1,ie+1
1296  ptc(i,j) = (u(i,j)-0.5*(va(i,j-1)+va(i,j))*cosa_v(i,j)) &
1297  *dyc(i,j)*sina_v(i,j)
1298  enddo
1299  enddo
1300 
1301  do j=js-1,je+1
1302  do i=is2,ie1
1303  vort(i,j) = (v(i,j) - 0.5*(ua(i-1,j)+ua(i,j))*cosa_u(i,j)) &
1304  *dxc(i,j)*sina_u(i,j)
1305  enddo
1306  enddo
1307 
1308  else
1309  do j=js,je+1
1310 
1311  if ( (j==1 .or. j==npy) ) then
1312  do i=is-1,ie+1
1313  if (vc(i,j) > 0) then
1314  ptc(i,j) = u(i,j)*dyc(i,j)*sin_sg(i,j-1,4)
1315  else
1316  ptc(i,j) = u(i,j)*dyc(i,j)*sin_sg(i,j,2)
1317  end if
1318  enddo
1319  else
1320  do i=is-1,ie+1
1321  ptc(i,j) = (u(i,j)-0.5*(va(i,j-1)+va(i,j))*cosa_v(i,j)) &
1322  *dyc(i,j)*sina_v(i,j)
1323  enddo
1324  endif
1325  enddo
1326 
1327  do j=js-1,je+1
1328  do i=is2,ie1
1329  vort(i,j) = (v(i,j) - 0.5*(ua(i-1,j)+ua(i,j))*cosa_u(i,j)) &
1330  *dxc(i,j)*sina_u(i,j)
1331  enddo
1332  if ( is == 1 ) then
1333  if (uc(1,j) > 0) then
1334  vort(1, j) = v(1, j)*dxc(1, j)*sin_sg(0,j,3)
1335  else
1336  vort(1, j) = v(1, j)*dxc(1, j)*sin_sg(1,j,1)
1337  end if
1338  end if
1339  if ( (ie+1)==npx ) then
1340  if (uc(npx,j) > 0) then
1341  vort(npx,j) = v(npx,j)*dxc(npx,j)* &
1342  sin_sg(npx-1,j,3)
1343  else
1344  vort(npx,j) = v(npx,j)*dxc(npx,j)* &
1345  sin_sg(npx,j,1)
1346  end if
1347  end if
1348  enddo
1349  endif
1350 
1351  do j=js,je+1
1352  do i=is,ie+1
1353  delpc(i,j) = vort(i,j-1) - vort(i,j) + ptc(i-1,j) - ptc(i,j)
1354  enddo
1355  enddo
1356 
1357 ! Remove the extra term at the corners:
1358  if (sw_corner) delpc(1, 1) = delpc(1, 1) - vort(1, 0)
1359  if (se_corner) delpc(npx, 1) = delpc(npx, 1) - vort(npx, 0)
1360  if (ne_corner) delpc(npx,npy) = delpc(npx,npy) + vort(npx,npy)
1361  if (nw_corner) delpc(1, npy) = delpc(1, npy) + vort(1, npy)
1362 
1363  do j=js,je+1
1364  do i=is,ie+1
1365  delpc(i,j) = gridstruct%rarea_c(i,j)*delpc(i,j)
1366  damp = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*abs(delpc(i,j)*dt)))
1367  vort(i,j) = damp*delpc(i,j)
1368  ke(i,j) = ke(i,j) + vort(i,j)
1369  enddo
1370  enddo
1371  else
1372 !--------------------------
1373 ! Higher order divg damping
1374 !--------------------------
1375  do j=js,je+1
1376  do i=is,ie+1
1377 ! Save divergence for external mode filter
1378  delpc(i,j) = divg_d(i,j)
1379  enddo
1380  enddo
1381 
1382  n2 = nord + 1 ! N > 1
1383  do n=1,nord
1384  nt = nord-n
1385 
1386  fill_c = (nt/=0) .and. (flagstruct%grid_type<3) .and. &
1387  ( sw_corner .or. se_corner .or. ne_corner .or. nw_corner ) &
1388  .and. .not. bounded_domain
1389 
1390  if ( fill_c ) call fill_corners(divg_d, npx, npy, fill=xdir, bgrid=.true.)
1391  do j=js-nt,je+1+nt
1392  do i=is-1-nt,ie+1+nt
1393  vc(i,j) = (divg_d(i+1,j)-divg_d(i,j))*divg_u(i,j)
1394  enddo
1395  enddo
1396 
1397  if ( fill_c ) call fill_corners(divg_d, npx, npy, fill=ydir, bgrid=.true.)
1398  do j=js-1-nt,je+1+nt
1399  do i=is-nt,ie+1+nt
1400  uc(i,j) = (divg_d(i,j+1)-divg_d(i,j))*divg_v(i,j)
1401  enddo
1402  enddo
1403 
1404  if ( fill_c ) call fill_corners(vc, uc, npx, npy, vector=.true., dgrid=.true.)
1405  do j=js-nt,je+1+nt
1406  do i=is-nt,ie+1+nt
1407  divg_d(i,j) = uc(i,j-1) - uc(i,j) + vc(i-1,j) - vc(i,j)
1408  enddo
1409  enddo
1410 
1411 ! Remove the extra term at the corners:
1412  if (sw_corner) divg_d(1, 1) = divg_d(1, 1) - uc(1, 0)
1413  if (se_corner) divg_d(npx, 1) = divg_d(npx, 1) - uc(npx, 0)
1414  if (ne_corner) divg_d(npx,npy) = divg_d(npx,npy) + uc(npx,npy)
1415  if (nw_corner) divg_d(1, npy) = divg_d(1, npy) + uc(1, npy)
1416 
1417  if ( .not. gridstruct%stretched_grid ) then
1418  do j=js-nt,je+1+nt
1419  do i=is-nt,ie+1+nt
1420  divg_d(i,j) = divg_d(i,j)*gridstruct%rarea_c(i,j)
1421  enddo
1422  enddo
1423  endif
1424 
1425  enddo ! n-loop
1426 
1427  if ( dddmp<1.e-5) then
1428  vort(:,:) = 0.
1429  else
1430  if ( flagstruct%grid_type < 3 ) then
1431 ! Interpolate relative vort to cell corners
1432  call a2b_ord4(wk, vort, gridstruct, npx, npy, is, ie, js, je, ng, .false.)
1433  do j=js,je+1
1434  do i=is,ie+1
1435 ! The following is an approxi form of Smagorinsky diffusion
1436  vort(i,j) = abs(dt)*sqrt(delpc(i,j)**2 + vort(i,j)**2)
1437  enddo
1438  enddo
1439  else ! Correct form: works only for doubly preiodic domain
1440  call smag_corner(abs(dt), u, v, ua, va, vort, bd, npx, npy, gridstruct, ng)
1441  endif
1442  endif
1443 
1444  if (gridstruct%stretched_grid ) then
1445 ! Stretched grid with variable damping ~ area
1446  dd8 = gridstruct%da_min * d4_bg**n2
1447  else
1448  dd8 = ( gridstruct%da_min_c*d4_bg )**n2
1449  endif
1450 
1451  do j=js,je+1
1452  do i=is,ie+1
1453  damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2
1454  vort(i,j) = damp2*delpc(i,j) + dd8*divg_d(i,j)
1455  ke(i,j) = ke(i,j) + vort(i,j)
1456  enddo
1457  enddo
1458 
1459  endif
1460 
1461  if ( d_con > 1.e-5 ) then
1462  do j=js,je+1
1463  do i=is,ie
1464  ub(i,j) = vort(i,j) - vort(i+1,j)
1465  enddo
1466  enddo
1467  do j=js,je
1468  do i=is,ie+1
1469  vb(i,j) = vort(i,j) - vort(i,j+1)
1470  enddo
1471  enddo
1472  endif
1473 
1474 ! Vorticity transport
1475  if ( hydrostatic ) then
1476  do j=jsd,jed
1477  do i=isd,ied
1478  vort(i,j) = wk(i,j) + f0(i,j)
1479  enddo
1480  enddo
1481  else
1482  if ( flagstruct%do_f3d ) then
1483  do j=jsd,jed
1484  do i=isd,ied
1485  vort(i,j) = wk(i,j) + f0(i,j)*z_rat(i,j)
1486  enddo
1487  enddo
1488  else
1489  do j=jsd,jed
1490  do i=isd,ied
1491  vort(i,j) = wk(i,j) + f0(i,j)
1492  enddo
1493  enddo
1494  endif
1495  endif
1496 
1497  call fv_tp_2d(vort, crx_adv, cry_adv, npx, npy, hord_vt, fx, fy, &
1498  xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac)
1499  do j=js,je+1
1500  do i=is,ie
1501  u(i,j) = vt(i,j) + ke(i,j) - ke(i+1,j) + fy(i,j)
1502  enddo
1503  enddo
1504  do j=js,je
1505  do i=is,ie+1
1506  v(i,j) = ut(i,j) + ke(i,j) - ke(i,j+1) - fx(i,j)
1507  enddo
1508  enddo
1509 
1510 !--------------------------------------------------------
1511 ! damping applied to relative vorticity (wk):
1512  if ( damp_v>1.e-5 ) then
1513  damp4 = (damp_v*gridstruct%da_min_c)**(nord_v+1)
1514  call del6_vt_flux(nord_v, npx, npy, damp4, wk, vort, ut, vt, gridstruct, bd)
1515  endif
1516 
1517  if ( d_con > 1.e-5 .or. flagstruct%do_skeb ) then
1518  do j=js,je+1
1519  do i=is,ie
1520  ub(i,j) = (ub(i,j) + vt(i,j))*rdx(i,j)
1521  fy(i,j) = u(i,j)*rdx(i,j)
1522  gy(i,j) = fy(i,j)*ub(i,j)
1523  enddo
1524  enddo
1525  do j=js,je
1526  do i=is,ie+1
1527  vb(i,j) = (vb(i,j) - ut(i,j))*rdy(i,j)
1528  fx(i,j) = v(i,j)*rdy(i,j)
1529  gx(i,j) = fx(i,j)*vb(i,j)
1530  enddo
1531  enddo
1532 !----------------------------------
1533 ! Heating due to damping:
1534 !----------------------------------
1535  damp = 0.25*d_con
1536  do j=js,je
1537  do i=is,ie
1538  u2 = fy(i,j) + fy(i,j+1)
1539  du2 = ub(i,j) + ub(i,j+1)
1540  v2 = fx(i,j) + fx(i+1,j)
1541  dv2 = vb(i,j) + vb(i+1,j)
1542 ! Total energy conserving:
1543 ! Convert lost KE due to divergence damping to "heat"
1544  heat_source(i,j) = delp(i,j)*(heat_source(i,j) - damp*rsin2(i,j)*( &
1545  (ub(i,j)**2 + ub(i,j+1)**2 + vb(i,j)**2 + vb(i+1,j)**2) &
1546  + 2.*(gy(i,j)+gy(i,j+1)+gx(i,j)+gx(i+1,j)) &
1547  - cosa_s(i,j)*(u2*dv2 + v2*du2 + du2*dv2)) )
1548  if (flagstruct%do_skeb) then
1549  diss_est(i,j) = diss_est(i,j)-rsin2(i,j)*( &
1550  (ub(i,j)**2 + ub(i,j+1)**2 + vb(i,j)**2 + vb(i+1,j)**2) &
1551  + 2.*(gy(i,j)+gy(i,j+1)+gx(i,j)+gx(i+1,j)) &
1552  - cosa_s(i,j)*(u2*dv2 + v2*du2 + du2*dv2))
1553  endif
1554  enddo
1555  enddo
1556  endif
1557 
1558 ! Add diffusive fluxes to the momentum equation:
1559  if ( damp_v>1.e-5 ) then
1560  do j=js,je+1
1561  do i=is,ie
1562  u(i,j) = u(i,j) + vt(i,j)
1563  enddo
1564  enddo
1565  do j=js,je
1566  do i=is,ie+1
1567  v(i,j) = v(i,j) - ut(i,j)
1568  enddo
1569  enddo
1570  endif
1571 
1572 #ifdef SW_DYNAMICS
1573  endif ! test_case
1574 #endif
1575 
1576  end subroutine d_sw
1577 
1580  subroutine del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd)
1581 ! Del-nord damping for the relative vorticity
1582 ! nord must be <= 2
1583 !------------------
1584 ! nord = 0: del-2
1585 ! nord = 1: del-4
1586 ! nord = 2: del-6
1587 !------------------
1588  integer, intent(in):: nord, npx, npy
1589  real, intent(in):: damp
1590  type(fv_grid_bounds_type), intent(IN) :: bd
1591  real, intent(inout):: q(bd%isd:bd%ied, bd%jsd:bd%jed) ! rel. vorticity ghosted on input
1592  type(fv_grid_type), intent(IN), target :: gridstruct
1593 ! Work arrays:
1594  real, intent(out):: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
1595  real, intent(out):: fx2(bd%isd:bd%ied+1,bd%jsd:bd%jed), fy2(bd%isd:bd%ied,bd%jsd:bd%jed+1)
1596  integer i,j, nt, n, i1, i2, j1, j2
1597 
1598  logical :: bounded_domain
1599 
1600 #ifdef USE_SG
1601  real, pointer, dimension(:,:,:) :: sin_sg
1602  real, pointer, dimension(:,:) :: rdxc, rdyc, dx,dy
1603 #endif
1604 
1605  integer :: is, ie, js, je
1606 
1607 #ifdef USE_SG
1608  sin_sg => gridstruct%sin_sg
1609  rdxc => gridstruct%rdxc
1610  rdyc => gridstruct%rdyc
1611  dx => gridstruct%dx
1612  dy => gridstruct%dy
1613 #endif
1614  bounded_domain = gridstruct%bounded_domain
1615 
1616  is = bd%is
1617  ie = bd%ie
1618  js = bd%js
1619  je = bd%je
1620 
1621  i1 = is-1-nord; i2 = ie+1+nord
1622  j1 = js-1-nord; j2 = je+1+nord
1623 
1624  do j=j1, j2
1625  do i=i1, i2
1626  d2(i,j) = damp*q(i,j)
1627  enddo
1628  enddo
1629 
1630  if( nord>0 .and. .not. bounded_domain) call copy_corners(d2, npx, npy, 1, bounded_domain, bd, gridstruct%sw_corner, &
1631  gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1632  do j=js-nord,je+nord
1633  do i=is-nord,ie+nord+1
1634 #ifdef USE_SG
1635  fx2(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*(d2(i-1,j)-d2(i,j))*rdxc(i,j)
1636 #else
1637  fx2(i,j) = gridstruct%del6_v(i,j)*(d2(i-1,j)-d2(i,j))
1638 #endif
1639  enddo
1640  enddo
1641 
1642  if( nord>0 .and. .not. bounded_domain) call copy_corners(d2, npx, npy, 2, bounded_domain, bd, gridstruct%sw_corner, &
1643  gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1644  do j=js-nord,je+nord+1
1645  do i=is-nord,ie+nord
1646 #ifdef USE_SG
1647  fy2(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*(d2(i,j-1)-d2(i,j))*rdyc(i,j)
1648 #else
1649  fy2(i,j) = gridstruct%del6_u(i,j)*(d2(i,j-1)-d2(i,j))
1650 #endif
1651  enddo
1652  enddo
1653 
1654  if ( nord>0 ) then
1655  do n=1, nord
1656  nt = nord-n
1657  do j=js-nt-1,je+nt+1
1658  do i=is-nt-1,ie+nt+1
1659  d2(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j)
1660  enddo
1661  enddo
1662 
1663  if (.not. bounded_domain) call copy_corners(d2, npx, npy, 1, bounded_domain, bd, gridstruct%sw_corner, &
1664  gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1665 
1666  do j=js-nt,je+nt
1667  do i=is-nt,ie+nt+1
1668 #ifdef USE_SG
1669  fx2(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))*dy(i,j)*(d2(i,j)-d2(i-1,j))*rdxc(i,j)
1670 #else
1671  fx2(i,j) = gridstruct%del6_v(i,j)*(d2(i,j)-d2(i-1,j))
1672 #endif
1673  enddo
1674  enddo
1675 
1676  if (.not. bounded_domain) call copy_corners(d2, npx, npy, 2, bounded_domain, bd, &
1677  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1678 
1679  do j=js-nt,je+nt+1
1680  do i=is-nt,ie+nt
1681 #ifdef USE_SG
1682  fy2(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))*dx(i,j)*(d2(i,j)-d2(i,j-1))*rdyc(i,j)
1683 #else
1684  fy2(i,j) = gridstruct%del6_u(i,j)*(d2(i,j)-d2(i,j-1))
1685 #endif
1686  enddo
1687  enddo
1688  enddo
1689  endif
1690 
1691  end subroutine del6_vt_flux
1692 
1695  subroutine divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
1696  type(fv_grid_bounds_type), intent(IN) :: bd
1697  real, intent(in), dimension(bd%isd:bd%ied, bd%jsd:bd%jed+1):: u
1698  real, intent(in), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ):: v
1699  real, intent(in), dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: ua, va
1700  real, intent(out), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed+1):: divg_d
1701  type(fv_grid_type), intent(IN), target :: gridstruct
1702  type(fv_flags_type), intent(IN), target :: flagstruct
1703 ! local
1704  real uf(bd%is-2:bd%ie+2,bd%js-1:bd%je+2)
1705  real vf(bd%is-1:bd%ie+2,bd%js-2:bd%je+2)
1706  integer i,j
1707  integer is2, ie1
1708 
1709  real, pointer, dimension(:,:,:) :: sin_sg, cos_sg
1710  real, pointer, dimension(:,:) :: dxc,dyc
1711 
1712  integer :: is, ie, js, je
1713  integer :: npx, npy
1714  logical :: bounded_domain
1715 
1716  is = bd%is
1717  ie = bd%ie
1718  js = bd%js
1719  je = bd%je
1720 
1721  npx = flagstruct%npx
1722  npy = flagstruct%npy
1723  bounded_domain = gridstruct%bounded_domain
1724 
1725  sin_sg => gridstruct%sin_sg
1726  cos_sg => gridstruct%cos_sg
1727  dxc => gridstruct%dxc
1728  dyc => gridstruct%dyc
1729 
1730  if (bounded_domain) then
1731  is2 = is; ie1 = ie+1
1732  else
1733  is2 = max(2,is); ie1 = min(npx-1,ie+1)
1734  end if
1735 
1736  if (flagstruct%grid_type==4) then
1737  do j=js-1,je+2
1738  do i=is-2,ie+2
1739  uf(i,j) = u(i,j)*dyc(i,j)
1740  enddo
1741  enddo
1742  do j=js-2,je+2
1743  do i=is-1,ie+2
1744  vf(i,j) = v(i,j)*dxc(i,j)
1745  enddo
1746  enddo
1747  do j=js-1,je+2
1748  do i=is-1,ie+2
1749  divg_d(i,j) = gridstruct%rarea_c(i,j)*(vf(i,j-1)-vf(i,j)+uf(i-1,j)-uf(i,j))
1750  enddo
1751  enddo
1752  else
1753 ! 9---4---8
1754 ! | |
1755 ! 1 5 3
1756 ! | |
1757 ! 6---2---7
1758  do j=js,je+1
1759  if ( j==1 .or. j==npy ) then
1760  do i=is-1,ie+1
1761  uf(i,j) = u(i,j)*dyc(i,j)*0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))
1762  enddo
1763  else
1764  do i=is-1,ie+1
1765  uf(i,j) = (u(i,j)-0.25*(va(i,j-1)+va(i,j))*(cos_sg(i,j-1,4)+cos_sg(i,j,2))) &
1766  * dyc(i,j)*0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))
1767  enddo
1768  endif
1769  enddo
1770 
1771  do j=js-1,je+1
1772  do i=is2,ie1
1773  vf(i,j) = (v(i,j) - 0.25*(ua(i-1,j)+ua(i,j))*(cos_sg(i-1,j,3)+cos_sg(i,j,1))) &
1774  *dxc(i,j)*0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))
1775  enddo
1776  if ( is == 1 ) vf(1, j) = v(1, j)*dxc(1, j)*0.5*(sin_sg(0,j,3)+sin_sg(1,j,1))
1777  if ( (ie+1)==npx ) vf(npx,j) = v(npx,j)*dxc(npx,j)*0.5*(sin_sg(npx-1,j,3)+sin_sg(npx,j,1))
1778  enddo
1779 
1780  do j=js,je+1
1781  do i=is,ie+1
1782  divg_d(i,j) = vf(i,j-1) - vf(i,j) + uf(i-1,j) - uf(i,j)
1783  enddo
1784  enddo
1785 
1786 ! Remove the extra term at the corners:
1787  if (gridstruct%sw_corner) divg_d(1, 1) = divg_d(1, 1) - vf(1, 0)
1788  if (gridstruct%se_corner) divg_d(npx, 1) = divg_d(npx, 1) - vf(npx, 0)
1789  if (gridstruct%ne_corner) divg_d(npx,npy) = divg_d(npx,npy) + vf(npx,npy)
1790  if (gridstruct%nw_corner) divg_d(1, npy) = divg_d(1, npy) + vf(1, npy)
1791 
1792  do j=js,je+1
1793  do i=is,ie+1
1794  divg_d(i,j) = gridstruct%rarea_c(i,j)*divg_d(i,j)
1795  enddo
1796  enddo
1797 
1798  endif
1799 
1800  end subroutine divergence_corner
1801 
1802  subroutine divergence_corner_nest(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
1803  type(fv_grid_bounds_type), intent(IN) :: bd
1804  real, intent(in), dimension(bd%isd:bd%ied, bd%jsd:bd%jed+1):: u
1805  real, intent(in), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed):: v
1806  real, intent(in), dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: ua, va
1807  real, intent(out), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed+1):: divg_d
1808  type(fv_grid_type), intent(IN), target :: gridstruct
1809  type(fv_flags_type), intent(IN), target :: flagstruct
1810 
1811 ! local
1812  real uf(bd%isd:bd%ied,bd%jsd:bd%jed+1)
1813  real vf(bd%isd:bd%ied+1,bd%jsd:bd%jed)
1814  integer i,j
1815 
1816 
1817  real, pointer, dimension(:,:) :: rarea_c
1818 
1819  real, pointer, dimension(:,:,:) :: sin_sg, cos_sg
1820  real, pointer, dimension(:,:) :: cosa_u, cosa_v
1821  real, pointer, dimension(:,:) :: sina_u, sina_v
1822  real, pointer, dimension(:,:) :: dxc,dyc
1823 
1824  integer :: isd, ied, jsd, jed
1825  integer :: npx, npy
1826 
1827  isd = bd%isd
1828  ied = bd%ied
1829  jsd = bd%jsd
1830  jed = bd%jed
1831 
1832  npx = flagstruct%npx
1833  npy = flagstruct%npy
1834 
1835  rarea_c => gridstruct%rarea_c
1836  sin_sg => gridstruct%sin_sg
1837  cos_sg => gridstruct%cos_sg
1838  cosa_u => gridstruct%cosa_u
1839  cosa_v => gridstruct%cosa_v
1840  sina_u => gridstruct%sina_u
1841  sina_v => gridstruct%sina_v
1842  dxc => gridstruct%dxc
1843  dyc => gridstruct%dyc
1844 
1845  divg_d = 1.e25
1846 
1847  if (flagstruct%grid_type==4) then
1848  do j=jsd,jed
1849  do i=isd,ied
1850  uf(i,j) = u(i,j)*dyc(i,j)
1851  enddo
1852  enddo
1853  do j=jsd,jed
1854  do i=isd,ied
1855  vf(i,j) = v(i,j)*dxc(i,j)
1856  enddo
1857  enddo
1858  do j=jsd+1,jed
1859  do i=isd+1,ied
1860  divg_d(i,j) = rarea_c(i,j)*(vf(i,j-1)-vf(i,j)+uf(i-1,j)-uf(i,j))
1861  enddo
1862  enddo
1863  else
1864 
1865  do j=jsd+1,jed
1866  do i=isd,ied
1867  uf(i,j) = (u(i,j)-0.25*(va(i,j-1)+va(i,j))*(cos_sg(i,j-1,4)+cos_sg(i,j,2))) &
1868  * dyc(i,j)*0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))
1869  enddo
1870  enddo
1871 
1872  do j=jsd,jed
1873  do i=isd+1,ied
1874  vf(i,j) = (v(i,j) - 0.25*(ua(i-1,j)+ua(i,j))*(cos_sg(i-1,j,3)+cos_sg(i,j,1))) &
1875  *dxc(i,j)*0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))
1876  enddo
1877  enddo
1878 
1879  do j=jsd+1,jed
1880  do i=isd+1,ied
1881  divg_d(i,j) = (vf(i,j-1) - vf(i,j) + uf(i-1,j) - uf(i,j))*rarea_c(i,j)
1882  enddo
1883  enddo
1884 
1885 !!$ !Edges
1886 !!$
1887 !!$ !West, East
1888 !!$ do j=jsd+1,jed
1889 !!$ divg_d(isd ,j) = (vf(isd,j-1) - vf(isd,j) + uf(isd,j) - uf(isd+1,j))*rarea_c(isd,j)
1890 !!$ divg_d(ied+1,j) = (vf(ied+1,j-1) - vf(ied+1,j) + uf(ied-1,j) - uf(ied,j))*rarea_c(ied,j)
1891 !!$ end do
1892 !!$
1893 !!$ !North, South
1894 !!$ do i=isd+1,ied
1895 !!$ divg_d(i,jsd ) = (vf(i,jsd) - vf(i,jsd+1) + uf(i-1,jsd) - uf(i,jsd))*rarea_c(i,jsd)
1896 !!$ divg_d(i,jed+1) = (vf(i,jed-1) - vf(i,jed) + uf(i-1,jed+1) - uf(i,jed+1))*rarea_c(i,jed)
1897 !!$ end do
1898 !!$
1899 !!$ !Corners (just use next corner value)
1900 !!$ divg_d(isd,jsd) = divg_d(isd+1,jsd+1)
1901 !!$ divg_d(isd,jed+1) = divg_d(isd+1,jed)
1902 !!$ divg_d(ied+1,jsd) = divg_d(ied,jsd+1)
1903 !!$ divg_d(ied+1,jed+1) = divg_d(ied,jed)
1904 
1905  endif
1906 
1907 
1908 end subroutine divergence_corner_nest
1909 
1911  subroutine smag_corner(dt, u, v, ua, va, smag_c, bd, npx, npy, gridstruct, ng)
1912 
1914  type(fv_grid_bounds_type), intent(IN) :: bd
1915  real, intent(in):: dt
1916  integer, intent(IN) :: npx, npy, ng
1917  real, intent(in), dimension(bd%isd:bd%ied, bd%jsd:bd%jed+1):: u
1918  real, intent(in), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ):: v
1919  real, intent(in), dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: ua, va
1920  real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: smag_c
1921  type(fv_grid_type), intent(IN), target :: gridstruct
1922 ! local
1923  real:: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed)
1924  real:: vt(bd%isd:bd%ied, bd%jsd:bd%jed+1)
1925  real:: wk(bd%isd:bd%ied,bd%jsd:bd%jed)
1926  real:: sh(bd%isd:bd%ied,bd%jsd:bd%jed)
1927  integer i,j
1928  integer is2, ie1
1929 
1930  real, pointer, dimension(:,:) :: dxc, dyc, dx, dy, rarea, rarea_c
1931 
1932  integer :: is, ie, js, je
1933  integer :: isd, ied, jsd, jed
1934 
1935  is = bd%is
1936  ie = bd%ie
1937  js = bd%js
1938  je = bd%je
1939 
1940  isd = bd%isd
1941  ied = bd%ied
1942  jsd = bd%jsd
1943  jed = bd%jed
1944 
1945  dxc => gridstruct%dxc
1946  dyc => gridstruct%dyc
1947  dx => gridstruct%dx
1948  dy => gridstruct%dy
1949  rarea => gridstruct%rarea
1950  rarea_c => gridstruct%rarea_c
1951 
1952  is2 = max(2,is); ie1 = min(npx-1,ie+1)
1953 
1954 ! Smag = sqrt [ T**2 + S**2 ]: unit = 1/s
1955 ! where T = du/dx - dv/dy; S = du/dy + dv/dx
1956 ! Compute tension strain at corners:
1957  do j=js,je+1
1958  do i=is-1,ie+1
1959  ut(i,j) = u(i,j)*dyc(i,j)
1960  enddo
1961  enddo
1962  do j=js-1,je+1
1963  do i=is,ie+1
1964  vt(i,j) = v(i,j)*dxc(i,j)
1965  enddo
1966  enddo
1967  do j=js,je+1
1968  do i=is,ie+1
1969  smag_c(i,j) = rarea_c(i,j)*(vt(i,j-1)-vt(i,j)-ut(i-1,j)+ut(i,j))
1970  enddo
1971  enddo
1972 ! Fix the corners?? if grid_type /= 4
1973 
1974 ! Compute shear strain:
1975  do j=jsd,jed+1
1976  do i=isd,ied
1977  vt(i,j) = u(i,j)*dx(i,j)
1978  enddo
1979  enddo
1980  do j=jsd,jed
1981  do i=isd,ied+1
1982  ut(i,j) = v(i,j)*dy(i,j)
1983  enddo
1984  enddo
1985 
1986  do j=jsd,jed
1987  do i=isd,ied
1988  wk(i,j) = rarea(i,j)*(vt(i,j)-vt(i,j+1)+ut(i,j)-ut(i+1,j))
1989  enddo
1990  enddo
1991  call a2b_ord4(wk, sh, gridstruct, npx, npy, is, ie, js, je, ng, .false.)
1992  do j=js,je+1
1993  do i=is,ie+1
1994  smag_c(i,j) = dt*sqrt( sh(i,j)**2 + smag_c(i,j)**2 )
1995  enddo
1996  enddo
1997 
1998  end subroutine smag_corner
1999 
2000 
2001  subroutine xtp_u(is,ie,js,je,isd,ied,jsd,jed,c, u, v, flux, iord, dx, rdx, npx, npy, grid_type, lim_fac,bounded_domain)
2002 
2003  integer, intent(in):: is,ie,js,je, isd,ied,jsd,jed
2004  real, INTENT(IN):: u(isd:ied,jsd:jed+1)
2005  real, INTENT(IN):: v(isd:ied+1,jsd:jed)
2006  real, INTENT(IN):: c(is:ie+1,js:je+1)
2007  real, INTENT(out):: flux(is:ie+1,js:je+1)
2008  real, INTENT(IN) :: dx(isd:ied, jsd:jed+1)
2009  real, INTENT(IN) :: rdx(isd:ied, jsd:jed+1)
2010  integer, INTENT(IN) :: iord, npx, npy, grid_type
2011  logical, INTENT(IN) :: bounded_domain
2012  real, INTENT(IN) :: lim_fac
2013 ! Local
2014  real, dimension(is-1:ie+1):: bl, br, b0
2015  logical, dimension(is-1:ie+1):: smt5, smt6
2016  logical, dimension(is:ie+1):: hi5, hi6
2017  real:: fx0(is:ie+1)
2018  real al(is-1:ie+2), dm(is-2:ie+2)
2019  real dq(is-3:ie+2)
2020  real dl, dr, xt, pmp, lac, cfl
2021  real pmp_1, lac_1, pmp_2, lac_2
2022  real x0, x1, x0L, x0R
2023  integer i, j
2024  integer is3, ie3
2025  integer is2, ie2
2026 
2027  if ( bounded_domain .or. grid_type>3 ) then
2028  is3 = is-1 ; ie3 = ie+1
2029  else
2030  is3 = max(3,is-1) ; ie3 = min(npx-3,ie+1)
2031  end if
2032 
2033 
2034  if ( iord < 8 ) then
2035 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
2036 
2037  do j=js,je+1
2038 
2039  do i=is3,ie3+1
2040  al(i) = p1*(u(i-1,j)+u(i,j)) + p2*(u(i-2,j)+u(i+1,j))
2041  enddo
2042  do i=is3,ie3
2043  bl(i) = al(i ) - u(i,j)
2044  br(i) = al(i+1) - u(i,j)
2045  enddo
2046 
2047  if ( (.not.bounded_domain) .and. grid_type < 3) then
2048  if ( is==1 ) then
2049  xt = c3*u(1,j) + c2*u(2,j) + c1*u(3,j)
2050  br(1) = xt - u(1,j)
2051  bl(2) = xt - u(2,j)
2052  br(2) = al(3) - u(2,j)
2053  if( j==1 .or. j==npy ) then
2054  bl(0) = 0. ! out
2055  br(0) = 0. ! edge
2056  bl(1) = 0. ! edge
2057  br(1) = 0. ! in
2058  else
2059  bl(0) = c1*u(-2,j) + c2*u(-1,j) + c3*u(0,j) - u(0,j)
2060  xt = 0.5*( ((2.*dx(0,j)+dx(-1,j))*(u(0,j))-dx(0,j)*u(-1,j))/(dx(0,j)+dx(-1,j)) &
2061  + ((2.*dx(1,j)+dx( 2,j))*(u(1,j))-dx(1,j)*u( 2,j))/(dx(1,j)+dx( 2,j)) )
2062  br(0) = xt - u(0,j)
2063  bl(1) = xt - u(1,j)
2064  endif
2065 ! call pert_ppm(1, u(2,j), bl(2), br(2), -1)
2066  endif
2067  if ( (ie+1)==npx ) then
2068  bl(npx-2) = al(npx-2) - u(npx-2,j)
2069  xt = c1*u(npx-3,j) + c2*u(npx-2,j) + c3*u(npx-1,j)
2070  br(npx-2) = xt - u(npx-2,j)
2071  bl(npx-1) = xt - u(npx-1,j)
2072  if( j==1 .or. j==npy ) then
2073  bl(npx-1) = 0. ! in
2074  br(npx-1) = 0. ! edge
2075  bl(npx ) = 0. ! edge
2076  br(npx ) = 0. ! out
2077  else
2078  xt = 0.5*( ( (2.*dx(npx-1,j)+dx(npx-2,j))*u(npx-1,j)-dx(npx-1,j)*u(npx-2,j))/(dx(npx-1,j)+dx(npx-2,j)) &
2079  + ( (2.*dx(npx ,j)+dx(npx+1,j))*u(npx ,j)-dx(npx ,j)*u(npx+1,j))/(dx(npx ,j)+dx(npx+1,j)) )
2080  br(npx-1) = xt - u(npx-1,j)
2081  bl(npx ) = xt - u(npx ,j)
2082  br(npx) = c3*u(npx,j) + c2*u(npx+1,j) + c1*u(npx+2,j) - u(npx,j)
2083  endif
2084 ! call pert_ppm(1, u(npx-2,j), bl(npx-2), br(npx-2), -1)
2085  endif
2086  endif
2087 
2088  do i=is-1,ie+1
2089  b0(i) = bl(i) + br(i)
2090  enddo
2091 
2092  if ( iord==1 ) then
2093 
2094  do i=is-1, ie+1
2095  smt5(i) = abs(lim_fac*b0(i)) < abs(bl(i)-br(i))
2096  enddo
2097 !DEC$ VECTOR ALWAYS
2098  do i=is,ie+1
2099  if( c(i,j)>0. ) then
2100  cfl = c(i,j)*rdx(i-1,j)
2101  fx0(i) = (1.-cfl)*(br(i-1)-cfl*b0(i-1))
2102  flux(i,j) = u(i-1,j)
2103  else
2104  cfl = c(i,j)*rdx(i,j)
2105  fx0(i) = (1.+cfl)*(bl(i)+cfl*b0(i))
2106  flux(i,j) = u(i,j)
2107  endif
2108  if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx0(i)
2109  enddo
2110 
2111  elseif ( iord==2 ) then ! Perfectly linear
2112 
2113 !DEC$ VECTOR ALWAYS
2114  do i=is,ie+1
2115  if( c(i,j)>0. ) then
2116  cfl = c(i,j)*rdx(i-1,j)
2117  flux(i,j) = u(i-1,j) + (1.-cfl)*(br(i-1)-cfl*b0(i-1))
2118  else
2119  cfl = c(i,j)*rdx(i,j)
2120  flux(i,j) = u(i,j) + (1.+cfl)*(bl(i)+cfl*b0(i))
2121  endif
2122  enddo
2123 
2124  elseif ( iord==3 ) then
2125 
2126  do i=is-1, ie+1
2127  x0 = abs(b0(i))
2128  x1 = abs(bl(i)-br(i))
2129  smt5(i) = x0 < x1
2130  smt6(i) = 3.*x0 < x1
2131  enddo
2132  do i=is, ie+1
2133  fx0(i) = 0.
2134  hi5(i) = smt5(i-1) .and. smt5(i)
2135  hi6(i) = smt6(i-1) .or. smt6(i)
2136  enddo
2137  do i=is, ie+1
2138  if( c(i,j)>0. ) then
2139  cfl = c(i,j)*rdx(i-1,j)
2140  if ( hi6(i) ) then
2141  fx0(i) = br(i-1) - cfl*b0(i-1)
2142  elseif( hi5(i) ) then
2143  fx0(i) = sign(min(abs(bl(i-1)),abs(br(i-1))), br(i-1))
2144  endif
2145  flux(i,j) = u(i-1,j) + (1.-cfl)*fx0(i)
2146  else
2147  cfl = c(i,j)*rdx(i,j)
2148  if ( hi6(i) ) then
2149  fx0(i) = bl(i) + cfl*b0(i)
2150  elseif( hi5(i) ) then
2151  fx0(i) = sign(min(abs(bl(i)),abs(br(i))), bl(i))
2152  endif
2153  flux(i,j) = u(i,j) + (1.+cfl)*fx0(i)
2154  endif
2155  enddo
2156 
2157  elseif ( iord==4 ) then
2158 
2159  do i=is-1, ie+1
2160  x0 = abs(b0(i))
2161  x1 = abs(bl(i)-br(i))
2162  smt5(i) = x0 < x1
2163  smt6(i) = 3.*x0 < x1
2164  enddo
2165  do i=is, ie+1
2166  hi5(i) = smt5(i-1) .and. smt5(i)
2167  hi6(i) = smt6(i-1) .or. smt6(i)
2168  hi5(i) = hi5(i) .or. hi6(i)
2169  enddo
2170 !DEC$ VECTOR ALWAYS
2171  do i=is,ie+1
2172  if( c(i,j)>0. ) then
2173  cfl = c(i,j)*rdx(i-1,j)
2174  fx0(i) = (1.-cfl)*(br(i-1)-cfl*b0(i-1))
2175  flux(i,j) = u(i-1,j)
2176  else
2177  cfl = c(i,j)*rdx(i,j)
2178  fx0(i) = (1.+cfl)*(bl(i)+cfl*b0(i))
2179  flux(i,j) = u(i,j)
2180  endif
2181  if ( hi5(i) ) flux(i,j) = flux(i,j) + fx0(i)
2182 
2183  enddo
2184 
2185  else ! iord=5,6,7
2186 
2187  if ( iord==5 ) then
2188  do i=is-1, ie+1
2189  smt5(i) = bl(i)*br(i) < 0.
2190  enddo
2191  else
2192  do i=is-1, ie+1
2193  smt5(i) = 3.*abs(b0(i)) < abs(bl(i)-br(i))
2194  enddo
2195  endif
2196 
2197 !DEC$ VECTOR ALWAYS
2198  do i=is,ie+1
2199  if( c(i,j)>0. ) then
2200  cfl = c(i,j)*rdx(i-1,j)
2201  fx0(i) = (1.-cfl)*(br(i-1)-cfl*b0(i-1))
2202  flux(i,j) = u(i-1,j)
2203  else
2204  cfl = c(i,j)*rdx(i,j)
2205  fx0(i) = (1.+cfl)*(bl(i)+cfl*b0(i))
2206  flux(i,j) = u(i,j)
2207  endif
2208  if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx0(i)
2209  enddo
2210 
2211  endif
2212  enddo
2213 
2214  else
2215  ! iord = 8, 9, 10, 11
2216 
2217  do j=js,je+1
2218  do i=is-2,ie+2
2219  xt = 0.25*(u(i+1,j) - u(i-1,j))
2220  dm(i) = sign(min(abs(xt), max(u(i-1,j), u(i,j), u(i+1,j)) - u(i,j), &
2221  u(i,j) - min(u(i-1,j), u(i,j), u(i+1,j))), xt)
2222  enddo
2223  do i=is-3,ie+2
2224  dq(i) = u(i+1,j) - u(i,j)
2225  enddo
2226 
2227  if (grid_type < 3) then
2228 
2229  do i=is3,ie3+1
2230  al(i) = 0.5*(u(i-1,j)+u(i,j)) + r3*(dm(i-1) - dm(i))
2231  enddo
2232 
2233 ! Perturbation form:
2234  if( iord==8 ) then
2235  do i=is3,ie3
2236  xt = 2.*dm(i)
2237  bl(i) = -sign(min(abs(xt), abs(al(i )-u(i,j))), xt)
2238  br(i) = sign(min(abs(xt), abs(al(i+1)-u(i,j))), xt)
2239  enddo
2240  elseif( iord==9 ) then
2241  do i=is3,ie3
2242  pmp_1 = -2.*dq(i)
2243  lac_1 = pmp_1 + 1.5*dq(i+1)
2244  bl(i) = min(max(0., pmp_1, lac_1), max(al(i )-u(i,j), min(0., pmp_1, lac_1)))
2245  pmp_2 = 2.*dq(i-1)
2246  lac_2 = pmp_2 - 1.5*dq(i-2)
2247  br(i) = min(max(0., pmp_2, lac_2), max(al(i+1)-u(i,j), min(0., pmp_2, lac_2)))
2248  enddo
2249  elseif( iord==10 ) then
2250  do i=is3,ie3
2251  bl(i) = al(i ) - u(i,j)
2252  br(i) = al(i+1) - u(i,j)
2253 ! if ( abs(dm(i-1))+abs(dm(i))+abs(dm(i+1)) < near_zero ) then
2254  if ( abs(dm(i)) < near_zero ) then
2255  if ( abs(dm(i-1))+abs(dm(i+1)) < near_zero ) then
2256 ! 2-delta-x structure detected within 3 cells
2257  bl(i) = 0.
2258  br(i) = 0.
2259  endif
2260  elseif( abs(3.*(bl(i)+br(i))) > abs(bl(i)-br(i)) ) then
2261  pmp_1 = -2.*dq(i)
2262  lac_1 = pmp_1 + 1.5*dq(i+1)
2263  bl(i) = min(max(0., pmp_1, lac_1), max(bl(i), min(0., pmp_1, lac_1)))
2264  pmp_2 = 2.*dq(i-1)
2265  lac_2 = pmp_2 - 1.5*dq(i-2)
2266  br(i) = min(max(0., pmp_2, lac_2), max(br(i), min(0., pmp_2, lac_2)))
2267  endif
2268  enddo
2269  else
2270 ! un-limited: 11
2271  do i=is3,ie3
2272  bl(i) = al(i ) - u(i,j)
2273  br(i) = al(i+1) - u(i,j)
2274  enddo
2275  endif
2276 
2277 !--------------
2278 ! fix the edges
2279 !--------------
2280 !!! TO DO: separate versions for bounded_domain and for cubed-sphere
2281  if ( is==1 .and. .not. bounded_domain) then
2282  br(2) = al(3) - u(2,j)
2283  xt = s15*u(1,j) + s11*u(2,j) - s14*dm(2)
2284  bl(2) = xt - u(2,j)
2285  br(1) = xt - u(1,j)
2286  if( j==1 .or. j==npy ) then
2287  bl(0) = 0. ! out
2288  br(0) = 0. ! edge
2289  bl(1) = 0. ! edge
2290  br(1) = 0. ! in
2291  else
2292  bl(0) = s14*dm(-1) - s11*dq(-1)
2293  x0l = 0.5*((2.*dx(0,j)+dx(-1,j))*(u(0,j)) &
2294  - dx(0,j)*(u(-1,j)))/(dx(0,j)+dx(-1,j))
2295  x0r = 0.5*((2.*dx(1,j)+dx(2,j))*(u(1,j)) &
2296  - dx(1,j)*(u(2,j)))/(dx(1,j)+dx(2,j))
2297  xt = x0l + x0r
2298  br(0) = xt - u(0,j)
2299  bl(1) = xt - u(1,j)
2300  endif
2301  call pert_ppm(1, u(2,j), bl(2), br(2), -1)
2302  endif
2303 
2304  if ( (ie+1)==npx .and. .not. bounded_domain) then
2305  bl(npx-2) = al(npx-2) - u(npx-2,j)
2306  xt = s15*u(npx-1,j) + s11*u(npx-2,j) + s14*dm(npx-2)
2307  br(npx-2) = xt - u(npx-2,j)
2308  bl(npx-1) = xt - u(npx-1,j)
2309  if( j==1 .or. j==npy ) then
2310  bl(npx-1) = 0. ! in
2311  br(npx-1) = 0. ! edge
2312  bl(npx ) = 0. ! edge
2313  br(npx ) = 0. ! out
2314  else
2315  br(npx) = s11*dq(npx) - s14*dm(npx+1)
2316  x0l = 0.5*( (2.*dx(npx-1,j)+dx(npx-2,j))*(u(npx-1,j)) &
2317  - dx(npx-1,j)*(u(npx-2,j)))/(dx(npx-1,j)+dx(npx-2,j))
2318  x0r = 0.5*( (2.*dx(npx,j)+dx(npx+1,j))*(u(npx,j)) &
2319  - dx(npx,j)*(u(npx+1,j)))/(dx(npx,j)+dx(npx+1,j))
2320  xt = x0l + x0r
2321  br(npx-1) = xt - u(npx-1,j)
2322  bl(npx ) = xt - u(npx ,j)
2323  endif
2324  call pert_ppm(1, u(npx-2,j), bl(npx-2), br(npx-2), -1)
2325  endif
2326 
2327  else
2328 ! Other grids:
2329  do i=is-1,ie+2
2330  al(i) = 0.5*(u(i-1,j)+u(i,j)) + r3*(dm(i-1) - dm(i))
2331  enddo
2332 
2333  do i=is-1,ie+1
2334  pmp = -2.*dq(i)
2335  lac = pmp + 1.5*dq(i+1)
2336  bl(i) = min(max(0., pmp, lac), max(al(i )-u(i,j), min(0.,pmp, lac)))
2337  pmp = 2.*dq(i-1)
2338  lac = pmp - 1.5*dq(i-2)
2339  br(i) = min(max(0., pmp, lac), max(al(i+1)-u(i,j), min(0.,pmp, lac)))
2340  enddo
2341  endif
2342 
2343  do i=is,ie+1
2344  if( c(i,j)>0. ) then
2345  cfl = c(i,j)*rdx(i-1,j)
2346  flux(i,j) = u(i-1,j) + (1.-cfl)*(br(i-1)-cfl*(bl(i-1)+br(i-1)))
2347  else
2348  cfl = c(i,j)*rdx(i,j)
2349  flux(i,j) = u(i, j) + (1.+cfl)*(bl(i )+cfl*(bl(i )+br(i )))
2350  endif
2351  enddo
2352  enddo
2353 
2354  endif
2355 
2356  end subroutine xtp_u
2357 
2358 
2359  subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx, npy, grid_type, lim_fac, bounded_domain)
2360  integer, intent(in):: is,ie,js,je, isd,ied,jsd,jed
2361  integer, intent(IN):: jord
2362  real, INTENT(IN) :: u(isd:ied,jsd:jed+1)
2363  real, INTENT(IN) :: v(isd:ied+1,jsd:jed)
2364  real, INTENT(IN) :: c(is:ie+1,js:je+1)
2365  real, INTENT(OUT):: flux(is:ie+1,js:je+1)
2366  real, INTENT(IN) :: dy(isd:ied+1,jsd:jed)
2367  real, INTENT(IN) :: rdy(isd:ied+1,jsd:jed)
2368  integer, INTENT(IN) :: npx, npy, grid_type
2369  logical, INTENT(IN) :: bounded_domain
2370  real, INTENT(IN) :: lim_fac
2371 ! Local:
2372  logical, dimension(is:ie+1,js-1:je+1):: smt5, smt6
2373  logical, dimension(is:ie+1):: hi5, hi6
2374  real:: fx0(is:ie+1)
2375  real dm(is:ie+1,js-2:je+2)
2376  real al(is:ie+1,js-1:je+2)
2377  real, dimension(is:ie+1,js-1:je+1):: bl, br, b0
2378  real dq(is:ie+1,js-3:je+2)
2379  real xt, dl, dr, pmp, lac, cfl
2380  real pmp_1, lac_1, pmp_2, lac_2
2381  real x0, x1, x0R, x0L
2382  integer i, j, is1, ie1, js3, je3
2383 
2384  if ( bounded_domain .or. grid_type>3 ) then
2385  js3 = js-1; je3 = je+1
2386  else
2387  js3 = max(3,js-1); je3 = min(npy-3,je+1)
2388  end if
2389 
2390  if ( jord<8 ) then
2391 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
2392 
2393  do j=js3,je3+1
2394  do i=is,ie+1
2395  al(i,j) = p1*(v(i,j-1)+v(i,j)) + p2*(v(i,j-2)+v(i,j+1))
2396  enddo
2397  enddo
2398  do j=js3,je3
2399  do i=is,ie+1
2400  bl(i,j) = al(i,j ) - v(i,j)
2401  br(i,j) = al(i,j+1) - v(i,j)
2402  enddo
2403  enddo
2404 
2405  if ( (.not.bounded_domain) .and. grid_type < 3) then
2406  if( js==1 ) then
2407  do i=is,ie+1
2408  bl(i,0) = c1*v(i,-2) + c2*v(i,-1) + c3*v(i,0) - v(i,0)
2409  xt = 0.5*( ((2.*dy(i,0)+dy(i,-1))*v(i,0)-dy(i,0)*v(i,-1))/(dy(i,0)+dy(i,-1)) &
2410  + ((2.*dy(i,1)+dy(i, 2))*v(i,1)-dy(i,1)*v(i, 2))/(dy(i,1)+dy(i, 2)) )
2411  br(i,0) = xt - v(i,0)
2412  bl(i,1) = xt - v(i,1)
2413  xt = c3*v(i,1) + c2*v(i,2) + c1*v(i,3)
2414  br(i,1) = xt - v(i,1)
2415  bl(i,2) = xt - v(i,2)
2416  br(i,2) = al(i,3) - v(i,2)
2417  enddo
2418  if ( is==1 ) then
2419  bl(1,0) = 0. ! out
2420  br(1,0) = 0. ! edge
2421  bl(1,1) = 0. ! edge
2422  br(1,1) = 0. ! in
2423  endif
2424  if ( (ie+1)==npx ) then
2425  bl(npx,0) = 0. ! out
2426  br(npx,0) = 0. ! edge
2427  bl(npx,1) = 0. ! edge
2428  br(npx,1) = 0. ! in
2429  endif
2430 ! j=2
2431 ! call pert_ppm(ie-is+2, v(is,j), bl(is,j), br(is,j), -1)
2432  endif
2433  if( (je+1)==npy ) then
2434  do i=is,ie+1
2435  bl(i,npy-2) = al(i,npy-2) - v(i,npy-2)
2436  xt = c1*v(i,npy-3) + c2*v(i,npy-2) + c3*v(i,npy-1)
2437  br(i,npy-2) = xt - v(i,npy-2)
2438  bl(i,npy-1) = xt - v(i,npy-1)
2439  xt = 0.5*( ((2.*dy(i,npy-1)+dy(i,npy-2))*v(i,npy-1)-dy(i,npy-1)*v(i,npy-2))/(dy(i,npy-1)+dy(i,npy-2)) &
2440  + ((2.*dy(i,npy )+dy(i,npy+1))*v(i,npy )-dy(i,npy )*v(i,npy+1))/(dy(i,npy )+dy(i,npy+1)) )
2441  br(i,npy-1) = xt - v(i,npy-1)
2442  bl(i,npy ) = xt - v(i,npy)
2443  br(i,npy) = c3*v(i,npy)+ c2*v(i,npy+1) + c1*v(i,npy+2) - v(i,npy)
2444  enddo
2445  if ( is==1 ) then
2446  bl(1,npy-1) = 0. ! in
2447  br(1,npy-1) = 0. ! edge
2448  bl(1,npy ) = 0. ! edge
2449  br(1,npy ) = 0. ! out
2450  endif
2451  if ( (ie+1)==npx ) then
2452  bl(npx,npy-1) = 0. ! in
2453  br(npx,npy-1) = 0. ! edge
2454  bl(npx,npy ) = 0. ! edge
2455  br(npx,npy ) = 0. ! out
2456  endif
2457 ! j=npy-2
2458 ! call pert_ppm(ie-is+2, v(is,j), bl(is,j), br(is,j), -1)
2459  endif
2460  endif
2461 
2462  do j=js-1,je+1
2463  do i=is,ie+1
2464  b0(i,j) = bl(i,j) + br(i,j)
2465  enddo
2466  enddo
2467 
2468 
2469  if ( jord==1 ) then ! Perfectly linear
2470 
2471  do j=js-1,je+1
2472  do i=is,ie+1
2473  smt5(i,j) = abs(lim_fac*b0(i,j)) < abs(bl(i,j)-br(i,j))
2474  enddo
2475  enddo
2476  do j=js,je+1
2477 !DEC$ VECTOR ALWAYS
2478  do i=is,ie+1
2479  if( c(i,j)>0. ) then
2480  cfl = c(i,j)*rdy(i,j-1)
2481  fx0(i) = (1.-cfl)*(br(i,j-1)-cfl*b0(i,j-1))
2482  flux(i,j) = v(i,j-1)
2483  else
2484  cfl = c(i,j)*rdy(i,j)
2485  fx0(i) = (1.+cfl)*(bl(i,j)+cfl*b0(i,j))
2486  flux(i,j) = v(i,j)
2487  endif
2488  if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx0(i)
2489  enddo
2490  enddo
2491 
2492  elseif ( jord==2 ) then ! Perfectly linear
2493  do j=js,je+1
2494 !DEC$ VECTOR ALWAYS
2495  do i=is,ie+1
2496  if( c(i,j)>0. ) then
2497  cfl = c(i,j)*rdy(i,j-1)
2498  flux(i,j) = v(i,j-1) + (1.-cfl)*(br(i,j-1)-cfl*b0(i,j-1))
2499  else
2500  cfl = c(i,j)*rdy(i,j)
2501  flux(i,j) = v(i,j) + (1.+cfl)*(bl(i,j)+cfl*b0(i,j))
2502  endif
2503  enddo
2504  enddo
2505 
2506  elseif ( jord==3 ) then
2507 
2508  do j=js-1,je+1
2509  do i=is,ie+1
2510  x0 = abs(b0(i,j))
2511  x1 = abs(bl(i,j)-br(i,j))
2512  smt5(i,j) = x0 < x1
2513  smt6(i,j) = 3.*x0 < x1
2514  enddo
2515  enddo
2516  do j=js,je+1
2517  do i=is,ie+1
2518  fx0(i) = 0.
2519  hi5(i) = smt5(i,j-1) .and. smt5(i,j)
2520  hi6(i) = smt6(i,j-1) .or. smt6(i,j)
2521  enddo
2522  do i=is,ie+1
2523  if( c(i,j)>0. ) then
2524  cfl = c(i,j)*rdy(i,j-1)
2525  if ( hi6(i) ) then
2526  fx0(i) = br(i,j-1) - cfl*b0(i,j-1)
2527  elseif ( hi5(i) ) then ! piece-wise linear
2528  fx0(i) = sign(min(abs(bl(i,j-1)),abs(br(i,j-1))), br(i,j-1))
2529  endif
2530  flux(i,j) = v(i,j-1) + (1.-cfl)*fx0(i)
2531  else
2532  cfl = c(i,j)*rdy(i,j)
2533  if ( hi6(i) ) then
2534  fx0(i) = bl(i,j) + cfl*b0(i,j)
2535  elseif ( hi5(i) ) then ! piece-wise linear
2536  fx0(i) = sign(min(abs(bl(i,j)),abs(br(i,j))), bl(i,j))
2537  endif
2538  flux(i,j) = v(i,j) + (1.+cfl)*fx0(i)
2539  endif
2540  enddo
2541  enddo
2542 
2543  elseif ( jord==4 ) then
2544 
2545  do j=js-1,je+1
2546  do i=is,ie+1
2547  x0 = abs(b0(i,j))
2548  x1 = abs(bl(i,j)-br(i,j))
2549  smt5(i,j) = x0 < x1
2550  smt6(i,j) = 3.*x0 < x1
2551  enddo
2552  enddo
2553  do j=js,je+1
2554  do i=is,ie+1
2555  fx0(i) = 0.
2556  hi5(i) = smt5(i,j-1) .and. smt5(i,j)
2557  hi6(i) = smt6(i,j-1) .or. smt6(i,j)
2558  hi5(i) = hi5(i) .or. hi6(i)
2559  enddo
2560 !DEC$ VECTOR ALWAYS
2561  do i=is,ie+1
2562  if( c(i,j)>0. ) then
2563  cfl = c(i,j)*rdy(i,j-1)
2564  fx0(i) = (1.-cfl)*(br(i,j-1)-cfl*b0(i,j-1))
2565  flux(i,j) = v(i,j-1)
2566  else
2567  cfl = c(i,j)*rdy(i,j)
2568  fx0(i) = (1.+cfl)*(bl(i,j)+cfl*b0(i,j))
2569  flux(i,j) = v(i,j)
2570  endif
2571  if ( hi5(i) ) flux(i,j) = flux(i,j) + fx0(i)
2572  enddo
2573  enddo
2574 
2575  else ! jord = 5,6,7
2576  if ( jord==5 ) then
2577  do j=js-1,je+1
2578  do i=is,ie+1
2579  smt5(i,j) = bl(i,j)*br(i,j) < 0.
2580  enddo
2581  enddo
2582  do j=js,je+1
2583 !DEC$ VECTOR ALWAYS
2584  do i=is,ie+1
2585  if( c(i,j)>0. ) then
2586  cfl = c(i,j)*rdy(i,j-1)
2587  fx0(i) = (1.-cfl)*(br(i,j-1)-cfl*b0(i,j-1))
2588  flux(i,j) = v(i,j-1)
2589  else
2590  cfl = c(i,j)*rdy(i,j)
2591  fx0(i) = (1.+cfl)*(bl(i,j)+cfl*b0(i,j))
2592  flux(i,j) = v(i,j)
2593  endif
2594  if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx0(i)
2595  enddo
2596  enddo
2597  else
2598 ! hord=6
2599  do j=js-1,je+1
2600  do i=is,ie+1
2601  smt6(i,j) = 3.*abs(b0(i,j)) < abs(bl(i,j)-br(i,j))
2602  enddo
2603  enddo
2604  do j=js,je+1
2605 !DEC$ VECTOR ALWAYS
2606  do i=is,ie+1
2607  if( c(i,j)>0. ) then
2608  cfl = c(i,j)*rdy(i,j-1)
2609  fx0(i) = (1.-cfl)*(br(i,j-1)-cfl*b0(i,j-1))
2610  flux(i,j) = v(i,j-1)
2611  else
2612  cfl = c(i,j)*rdy(i,j)
2613  fx0(i) = (1.+cfl)*(bl(i,j)+cfl*b0(i,j))
2614  flux(i,j) = v(i,j)
2615  endif
2616  if (smt6(i,j-1).or.smt6(i,j)) flux(i,j) = flux(i,j) + fx0(i)
2617  enddo
2618  enddo
2619  endif
2620 
2621  endif
2622 
2623  else
2624 ! jord= 8, 9, 10
2625 
2626  do j=js-2,je+2
2627  do i=is,ie+1
2628  xt = 0.25*(v(i,j+1) - v(i,j-1))
2629  dm(i,j) = sign(min(abs(xt), max(v(i,j-1), v(i,j), v(i,j+1)) - v(i,j), &
2630  v(i,j) - min(v(i,j-1), v(i,j), v(i,j+1))), xt)
2631  enddo
2632  enddo
2633 
2634  do j=js-3,je+2
2635  do i=is,ie+1
2636  dq(i,j) = v(i,j+1) - v(i,j)
2637  enddo
2638  enddo
2639 
2640  if (grid_type < 3) then
2641  do j=js3,je3+1
2642  do i=is,ie+1
2643  al(i,j) = 0.5*(v(i,j-1)+v(i,j)) + r3*(dm(i,j-1)-dm(i,j))
2644  enddo
2645  enddo
2646 
2647  if ( jord==8 ) then
2648  do j=js3,je3
2649  do i=is,ie+1
2650  xt = 2.*dm(i,j)
2651  bl(i,j) = -sign(min(abs(xt), abs(al(i,j)-v(i,j))), xt)
2652  br(i,j) = sign(min(abs(xt), abs(al(i,j+1)-v(i,j))), xt)
2653  enddo
2654  enddo
2655  elseif ( jord==9 ) then
2656  do j=js3,je3
2657  do i=is,ie+1
2658  pmp_1 = -2.*dq(i,j)
2659  lac_1 = pmp_1 + 1.5*dq(i,j+1)
2660  bl(i,j) = min(max(0., pmp_1, lac_1), max(al(i,j)-v(i,j), min(0., pmp_1, lac_1)))
2661  pmp_2 = 2.*dq(i,j-1)
2662  lac_2 = pmp_2 - 1.5*dq(i,j-2)
2663  br(i,j) = min(max(0., pmp_2, lac_2), max(al(i,j+1)-v(i,j), min(0., pmp_2, lac_2)))
2664  enddo
2665  enddo
2666  elseif ( jord==10 ) then
2667  do j=js3,je3
2668  do i=is,ie+1
2669  bl(i,j) = al(i,j ) - v(i,j)
2670  br(i,j) = al(i,j+1) - v(i,j)
2671 ! if ( abs(dm(i,j-1))+abs(dm(i,j))+abs(dm(i,j+1)) < near_zero ) then
2672  if ( abs(dm(i,j)) < near_zero ) then
2673  if ( abs(dm(i,j-1))+abs(dm(i,j+1)) < near_zero ) then
2674  bl(i,j) = 0.
2675  br(i,j) = 0.
2676  endif
2677  elseif( abs(3.*(bl(i,j)+br(i,j))) > abs(bl(i,j)-br(i,j)) ) then
2678  pmp_1 = -2.*dq(i,j)
2679  lac_1 = pmp_1 + 1.5*dq(i,j+1)
2680  bl(i,j) = min(max(0., pmp_1, lac_1), max(bl(i,j), min(0., pmp_1, lac_1)))
2681  pmp_2 = 2.*dq(i,j-1)
2682  lac_2 = pmp_2 - 1.5*dq(i,j-2)
2683  br(i,j) = min(max(0., pmp_2, lac_2), max(br(i,j), min(0., pmp_2, lac_2)))
2684  endif
2685  enddo
2686  enddo
2687  else
2688 ! Unlimited:
2689  do j=js3,je3
2690  do i=is,ie+1
2691  bl(i,j) = al(i,j ) - v(i,j)
2692  br(i,j) = al(i,j+1) - v(i,j)
2693  enddo
2694  enddo
2695  endif
2696 
2697 !--------------
2698 ! fix the edges
2699 !--------------
2700  if( js==1 .and. .not. bounded_domain) then
2701  do i=is,ie+1
2702  br(i,2) = al(i,3) - v(i,2)
2703  xt = s15*v(i,1) + s11*v(i,2) - s14*dm(i,2)
2704  br(i,1) = xt - v(i,1)
2705  bl(i,2) = xt - v(i,2)
2706 
2707  bl(i,0) = s14*dm(i,-1) - s11*dq(i,-1)
2708 
2709 #ifdef ONE_SIDE
2710  xt = t14*v(i,1) + t12*v(i,2) + t15*v(i,3)
2711  bl(i,1) = 2.*xt - v(i,1)
2712  xt = t14*v(i,0) + t12*v(i,-1) + t15*v(i,-2)
2713  br(i,0) = 2.*xt - v(i,0)
2714 #else
2715  x0l = 0.5*( (2.*dy(i,0)+dy(i,-1))*(v(i,0)) &
2716  - dy(i,0)*(v(i,-1)))/(dy(i,0)+dy(i,-1))
2717  x0r = 0.5*( (2.*dy(i,1)+dy(i,2))*(v(i,1)) &
2718  - dy(i,1)*(v(i,2)))/(dy(i,1)+dy(i,2))
2719  xt = x0l + x0r
2720 
2721  bl(i,1) = xt - v(i,1)
2722  br(i,0) = xt - v(i,0)
2723 #endif
2724  enddo
2725  if ( is==1 ) then
2726  bl(1,0) = 0. ! out
2727  br(1,0) = 0. ! edge
2728  bl(1,1) = 0. ! edge
2729  br(1,1) = 0. ! in
2730  endif
2731  if ( (ie+1)==npx ) then
2732  bl(npx,0) = 0. ! out
2733  br(npx,0) = 0. ! edge
2734  bl(npx,1) = 0. ! edge
2735  br(npx,1) = 0. ! in
2736  endif
2737  j=2
2738  call pert_ppm(ie-is+2, v(is,j), bl(is,j), br(is,j), -1)
2739  endif
2740  if( (je+1)==npy .and. .not. bounded_domain) then
2741  do i=is,ie+1
2742  bl(i,npy-2) = al(i,npy-2) - v(i,npy-2)
2743  xt = s15*v(i,npy-1) + s11*v(i,npy-2) + s14*dm(i,npy-2)
2744  br(i,npy-2) = xt - v(i,npy-2)
2745  bl(i,npy-1) = xt - v(i,npy-1)
2746  br(i,npy) = s11*dq(i,npy) - s14*dm(i,npy+1)
2747 #ifdef ONE_SIDE
2748  xt = t14*v(i,npy-1) + t12*v(i,npy-2) + t15*v(i,npy-3)
2749  br(i,npy-1) = 2.*xt - v(i,npy-1)
2750  xt = t14*v(i,npy) + t12*v(i,npy+1) + t15*v(i,npy+2)
2751  bl(i,npy ) = 2.*xt - v(i,npy)
2752 #else
2753  x0l= 0.5*((2.*dy(i,npy-1)+dy(i,npy-2))*(v(i,npy-1)) - &
2754  dy(i,npy-1)*(v(i,npy-2)))/(dy(i,npy-1)+dy(i,npy-2))
2755  x0r= 0.5*((2.*dy(i,npy)+dy(i,npy+1))*(v(i,npy)) - &
2756  dy(i,npy)*(v(i,npy+1)))/(dy(i,npy)+dy(i,npy+1))
2757  xt = x0l + x0r
2758 
2759  br(i,npy-1) = xt - v(i,npy-1)
2760  bl(i,npy ) = xt - v(i,npy)
2761 #endif
2762  enddo
2763  if ( is==1 ) then
2764  bl(1,npy-1) = 0. ! in
2765  br(1,npy-1) = 0. ! edge
2766  bl(1,npy ) = 0. ! edge
2767  br(1,npy ) = 0. ! out
2768  endif
2769  if ( (ie+1)==npx ) then
2770  bl(npx,npy-1) = 0. ! in
2771  br(npx,npy-1) = 0. ! edge
2772  bl(npx,npy ) = 0. ! edge
2773  br(npx,npy ) = 0. ! out
2774  endif
2775  j=npy-2
2776  call pert_ppm(ie-is+2, v(is,j), bl(is,j), br(is,j), -1)
2777  endif
2778 
2779  else
2780 
2781  do j=js-1,je+2
2782  do i=is,ie+1
2783  al(i,j) = 0.5*(v(i,j-1)+v(i,j)) + r3*(dm(i,j-1)-dm(i,j))
2784  enddo
2785  enddo
2786 
2787  do j=js-1,je+1
2788  do i=is,ie+1
2789  pmp = 2.*dq(i,j-1)
2790  lac = pmp - 1.5*dq(i,j-2)
2791  br(i,j) = min(max(0.,pmp,lac), max(al(i,j+1)-v(i,j), min(0.,pmp,lac)))
2792  pmp = -2.*dq(i,j)
2793  lac = pmp + 1.5*dq(i,j+1)
2794  bl(i,j) = min(max(0.,pmp,lac), max(al(i,j)-v(i,j), min(0.,pmp,lac)))
2795  enddo
2796  enddo
2797 
2798  endif
2799 
2800  do j=js,je+1
2801  do i=is,ie+1
2802  if(c(i,j)>0.) then
2803  cfl = c(i,j)*rdy(i,j-1)
2804  flux(i,j) = v(i,j-1) + (1.-cfl)*(br(i,j-1)-cfl*(bl(i,j-1)+br(i,j-1)))
2805  else
2806  cfl = c(i,j)*rdy(i,j)
2807  flux(i,j) = v(i,j ) + (1.+cfl)*(bl(i,j )+cfl*(bl(i,j )+br(i,j )))
2808  endif
2809  enddo
2810  enddo
2811 
2812  endif
2813 
2814 end subroutine ytp_v
2815 
2816 
2817 !There is a limit to how far this routine can fill uc and vc in the
2818 ! halo, and so either mpp_update_domains or some sort of boundary
2819 ! routine (extrapolation, outflow, interpolation from a bounded_domain grid)
2820 
2821 ! is needed after c_sw is completed if these variables are needed
2822 ! in the halo
2823  subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
2824  bd, npx, npy, bounded_domain, grid_type)
2825  type(fv_grid_bounds_type), intent(IN) :: bd
2826  logical, intent(in):: dord4
2827  real, intent(in) :: u(bd%isd:bd%ied,bd%jsd:bd%jed+1)
2828  real, intent(in) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed)
2829  real, intent(out), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed ):: uc
2830  real, intent(out), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed+1):: vc
2831  real, intent(out), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ):: ua, va, ut, vt
2832  integer, intent(IN) :: npx, npy, grid_type
2833  logical, intent(IN) :: bounded_domain
2834  type(fv_grid_type), intent(IN), target :: gridstruct
2835 ! Local
2836  real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: utmp, vtmp
2837  integer npt, i, j, ifirst, ilast, id
2838  integer :: is, ie, js, je
2839  integer :: isd, ied, jsd, jed
2840 
2841  real, pointer, dimension(:,:,:) :: sin_sg
2842  real, pointer, dimension(:,:) :: cosa_u, cosa_v, cosa_s
2843  real, pointer, dimension(:,:) :: rsin_u, rsin_v, rsin2
2844  real, pointer, dimension(:,:) :: dxa,dya
2845 
2846  is = bd%is
2847  ie = bd%ie
2848  js = bd%js
2849  je = bd%je
2850  isd = bd%isd
2851  ied = bd%ied
2852  jsd = bd%jsd
2853  jed = bd%jed
2854 
2855  sin_sg => gridstruct%sin_sg
2856  cosa_u => gridstruct%cosa_u
2857  cosa_v => gridstruct%cosa_v
2858  cosa_s => gridstruct%cosa_s
2859  rsin_u => gridstruct%rsin_u
2860  rsin_v => gridstruct%rsin_v
2861  rsin2 => gridstruct%rsin2
2862  dxa => gridstruct%dxa
2863  dya => gridstruct%dya
2864 
2865  if ( dord4 ) then
2866  id = 1
2867  else
2868  id = 0
2869  endif
2870 
2871  if (grid_type < 3 .and. .not. bounded_domain) then
2872  npt = 4
2873  else
2874  npt = -2
2875  endif
2876 
2877 ! Initialize the non-existing corner regions
2878  utmp(:,:) = big_number
2879  vtmp(:,:) = big_number
2880 
2881  if ( bounded_domain ) then
2882 
2883  do j=jsd+1,jed-1
2884  do i=isd,ied
2885  utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
2886  enddo
2887  enddo
2888  do i=isd,ied
2889  j = jsd
2890  utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
2891  j = jed
2892  utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
2893  end do
2894 
2895  do j=jsd,jed
2896  do i=isd+1,ied-1
2897  vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
2898  enddo
2899  i = isd
2900  vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
2901  i = ied
2902  vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
2903  enddo
2904 
2905  do j=jsd,jed
2906  do i=isd,ied
2907  ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
2908  va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
2909  enddo
2910  enddo
2911 
2912  else
2913  !----------
2914  ! Interior:
2915  !----------
2916  do j=max(npt,js-1),min(npy-npt,je+1)
2917  do i=max(npt,isd),min(npx-npt,ied)
2918  utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
2919  enddo
2920  enddo
2921  do j=max(npt,jsd),min(npy-npt,jed)
2922  do i=max(npt,is-1),min(npx-npt,ie+1)
2923  vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
2924  enddo
2925  enddo
2926 
2927  !----------
2928  ! edges:
2929  !----------
2930  if (grid_type < 3) then
2931 
2932  if ( js==1 .or. jsd<npt) then
2933  do j=jsd,npt-1
2934  do i=isd,ied
2935  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
2936  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
2937  enddo
2938  enddo
2939  endif
2940  if ( (je+1)==npy .or. jed>=(npy-npt)) then
2941  do j=npy-npt+1,jed
2942  do i=isd,ied
2943  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
2944  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
2945  enddo
2946  enddo
2947  endif
2948 
2949  if ( is==1 .or. isd<npt ) then
2950  do j=max(npt,jsd),min(npy-npt,jed)
2951  do i=isd,npt-1
2952  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
2953  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
2954  enddo
2955  enddo
2956  endif
2957  if ( (ie+1)==npx .or. ied>=(npx-npt)) then
2958  do j=max(npt,jsd),min(npy-npt,jed)
2959  do i=npx-npt+1,ied
2960  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
2961  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
2962  enddo
2963  enddo
2964  endif
2965 
2966  endif
2967 
2968 ! Contra-variant components at cell center:
2969  do j=js-1-id,je+1+id
2970  do i=is-1-id,ie+1+id
2971  ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
2972  va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
2973  enddo
2974  enddo
2975 
2976  end if
2977 
2978 ! A -> C
2979 !--------------
2980 ! Fix the edges
2981 !--------------
2982 ! Xdir:
2983  if( gridstruct%sw_corner ) then
2984  do i=-2,0
2985  utmp(i,0) = -vtmp(0,1-i)
2986  enddo
2987  endif
2988  if( gridstruct%se_corner ) then
2989  do i=0,2
2990  utmp(npx+i,0) = vtmp(npx,i+1)
2991  enddo
2992  endif
2993  if( gridstruct%ne_corner ) then
2994  do i=0,2
2995  utmp(npx+i,npy) = -vtmp(npx,je-i)
2996  enddo
2997  endif
2998  if( gridstruct%nw_corner ) then
2999  do i=-2,0
3000  utmp(i,npy) = vtmp(0,je+i)
3001  enddo
3002  endif
3003 
3004  if (grid_type < 3 .and. .not. bounded_domain) then
3005  ifirst = max(3, is-1)
3006  ilast = min(npx-2,ie+2)
3007  else
3008  ifirst = is-1
3009  ilast = ie+2
3010  endif
3011 !---------------------------------------------
3012 ! 4th order interpolation for interior points:
3013 !---------------------------------------------
3014  do j=js-1,je+1
3015  do i=ifirst,ilast
3016  uc(i,j) = a2*(utmp(i-2,j)+utmp(i+1,j)) + a1*(utmp(i-1,j)+utmp(i,j))
3017  ut(i,j) = (uc(i,j) - v(i,j)*cosa_u(i,j))*rsin_u(i,j)
3018  enddo
3019  enddo
3020 
3021  if (grid_type < 3) then
3022 ! Xdir:
3023  if( gridstruct%sw_corner ) then
3024  ua(-1,0) = -va(0,2)
3025  ua( 0,0) = -va(0,1)
3026  endif
3027  if( gridstruct%se_corner ) then
3028  ua(npx, 0) = va(npx,1)
3029  ua(npx+1,0) = va(npx,2)
3030  endif
3031  if( gridstruct%ne_corner ) then
3032  ua(npx, npy) = -va(npx,npy-1)
3033  ua(npx+1,npy) = -va(npx,npy-2)
3034  endif
3035  if( gridstruct%nw_corner ) then
3036  ua(-1,npy) = va(0,npy-2)
3037  ua( 0,npy) = va(0,npy-1)
3038  endif
3039 
3040  if( is==1 .and. .not. bounded_domain ) then
3041  do j=js-1,je+1
3042  uc(0,j) = c1*utmp(-2,j) + c2*utmp(-1,j) + c3*utmp(0,j)
3043  ut(1,j) = edge_interpolate4(ua(-1:2,j), dxa(-1:2,j))
3044  !Want to use the UPSTREAM value
3045  if (ut(1,j) > 0.) then
3046  uc(1,j) = ut(1,j)*sin_sg(0,j,3)
3047  else
3048  uc(1,j) = ut(1,j)*sin_sg(1,j,1)
3049  end if
3050  uc(2,j) = c1*utmp(3,j) + c2*utmp(2,j) + c3*utmp(1,j)
3051  ut(0,j) = (uc(0,j) - v(0,j)*cosa_u(0,j))*rsin_u(0,j)
3052  ut(2,j) = (uc(2,j) - v(2,j)*cosa_u(2,j))*rsin_u(2,j)
3053  enddo
3054  endif
3055 
3056  if( (ie+1)==npx .and. .not. bounded_domain ) then
3057  do j=js-1,je+1
3058  uc(npx-1,j) = c1*utmp(npx-3,j)+c2*utmp(npx-2,j)+c3*utmp(npx-1,j)
3059  ut(npx, j) = edge_interpolate4(ua(npx-2:npx+1,j), dxa(npx-2:npx+1,j))
3060  if (ut(npx,j) > 0.) then
3061  uc(npx,j) = ut(npx,j)*sin_sg(npx-1,j,3)
3062  else
3063  uc(npx,j) = ut(npx,j)*sin_sg(npx,j,1)
3064  end if
3065  uc(npx+1,j) = c3*utmp(npx,j) + c2*utmp(npx+1,j) + c1*utmp(npx+2,j)
3066  ut(npx-1,j) = (uc(npx-1,j)-v(npx-1,j)*cosa_u(npx-1,j))*rsin_u(npx-1,j)
3067  ut(npx+1,j) = (uc(npx+1,j)-v(npx+1,j)*cosa_u(npx+1,j))*rsin_u(npx+1,j)
3068  enddo
3069  endif
3070 
3071  endif
3072 
3073 !------
3074 ! Ydir:
3075 !------
3076  if( gridstruct%sw_corner ) then
3077  do j=-2,0
3078  vtmp(0,j) = -utmp(1-j,0)
3079  enddo
3080  endif
3081  if( gridstruct%nw_corner ) then
3082  do j=0,2
3083  vtmp(0,npy+j) = utmp(j+1,npy)
3084  enddo
3085  endif
3086  if( gridstruct%se_corner ) then
3087  do j=-2,0
3088  vtmp(npx,j) = utmp(ie+j,0)
3089  enddo
3090  endif
3091  if( gridstruct%ne_corner ) then
3092  do j=0,2
3093  vtmp(npx,npy+j) = -utmp(ie-j,npy)
3094  enddo
3095  endif
3096  if( gridstruct%sw_corner ) then
3097  va(0,-1) = -ua(2,0)
3098  va(0, 0) = -ua(1,0)
3099  endif
3100  if( gridstruct%se_corner ) then
3101  va(npx, 0) = ua(npx-1,0)
3102  va(npx,-1) = ua(npx-2,0)
3103  endif
3104  if( gridstruct%ne_corner ) then
3105  va(npx,npy ) = -ua(npx-1,npy)
3106  va(npx,npy+1) = -ua(npx-2,npy)
3107  endif
3108  if( gridstruct%nw_corner ) then
3109  va(0,npy) = ua(1,npy)
3110  va(0,npy+1) = ua(2,npy)
3111  endif
3112 
3113  if (grid_type < 3) then
3114 
3115  do j=js-1,je+2
3116  if ( j==1 .and. .not. bounded_domain ) then
3117  do i=is-1,ie+1
3118  vt(i,j) = edge_interpolate4(va(i,-1:2), dya(i,-1:2))
3119  if (vt(i,j) > 0.) then
3120  vc(i,j) = vt(i,j)*sin_sg(i,j-1,4)
3121  else
3122  vc(i,j) = vt(i,j)*sin_sg(i,j,2)
3123  end if
3124  enddo
3125  elseif ( j==0 .or. j==(npy-1) .and. .not. bounded_domain ) then
3126  do i=is-1,ie+1
3127  vc(i,j) = c1*vtmp(i,j-2) + c2*vtmp(i,j-1) + c3*vtmp(i,j)
3128  vt(i,j) = (vc(i,j) - u(i,j)*cosa_v(i,j))*rsin_v(i,j)
3129  enddo
3130  elseif ( j==2 .or. j==(npy+1) .and. .not. bounded_domain ) then
3131  do i=is-1,ie+1
3132  vc(i,j) = c1*vtmp(i,j+1) + c2*vtmp(i,j) + c3*vtmp(i,j-1)
3133  vt(i,j) = (vc(i,j) - u(i,j)*cosa_v(i,j))*rsin_v(i,j)
3134  enddo
3135  elseif ( j==npy .and. .not. bounded_domain ) then
3136  do i=is-1,ie+1
3137  vt(i,j) = edge_interpolate4(va(i,j-2:j+1), dya(i,j-2:j+1))
3138  if (vt(i,j) > 0.) then
3139  vc(i,j) = vt(i,j)*sin_sg(i,j-1,4)
3140  else
3141  vc(i,j) = vt(i,j)*sin_sg(i,j,2)
3142  end if
3143  enddo
3144  else
3145 ! 4th order interpolation for interior points:
3146  do i=is-1,ie+1
3147  vc(i,j) = a2*(vtmp(i,j-2)+vtmp(i,j+1)) + a1*(vtmp(i,j-1)+vtmp(i,j))
3148  vt(i,j) = (vc(i,j) - u(i,j)*cosa_v(i,j))*rsin_v(i,j)
3149  enddo
3150  endif
3151  enddo
3152  else
3153 ! 4th order interpolation:
3154  do j=js-1,je+2
3155  do i=is-1,ie+1
3156  vc(i,j) = a2*(vtmp(i,j-2)+vtmp(i,j+1)) + a1*(vtmp(i,j-1)+vtmp(i,j))
3157  vt(i,j) = vc(i,j)
3158  enddo
3159  enddo
3160  endif
3161 
3162  end subroutine d2a2c_vect
3163 
3164 
3165  real function edge_interpolate4(ua, dxa)
3166 
3167  real, intent(in) :: ua(4)
3168  real, intent(in) :: dxa(4)
3169  real:: t1, t2
3170 
3171  t1 = dxa(1) + dxa(2)
3172  t2 = dxa(3) + dxa(4)
3173  edge_interpolate4 = 0.5*( ((t1+dxa(2))*ua(2)-dxa(2)*ua(1)) / t1 + &
3174  ((t2+dxa(3))*ua(3)-dxa(3)*ua(4)) / t2 )
3175 
3176  end function edge_interpolate4
3177 
3179  subroutine fill3_4corners(q1, q2, q3, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
3180  type(fv_grid_bounds_type), intent(IN) :: bd
3181  integer, intent(in):: dir
3182  real, intent(inout):: q1(bd%isd:bd%ied,bd%jsd:bd%jed)
3183  real, intent(inout):: q2(bd%isd:bd%ied,bd%jsd:bd%jed)
3184  real, intent(inout):: q3(bd%isd:bd%ied,bd%jsd:bd%jed)
3185  logical, intent(IN) :: sw_corner, se_corner, ne_corner, nw_corner
3186  integer, intent(IN) :: npx, npy
3187  integer i,j
3188 
3189  integer :: is, ie, js, je
3190  integer :: isd, ied, jsd, jed
3191 
3192  is = bd%is
3193  ie = bd%ie
3194  js = bd%js
3195  je = bd%je
3196  isd = bd%isd
3197  ied = bd%ied
3198  jsd = bd%jsd
3199  jed = bd%jed
3200 
3201  select case(dir)
3202  case(1)
3203  if ( sw_corner ) then
3204  q1(-1,0) = q1(0,2); q1(0,0) = q1(0,1); q1(0,-1) = q1(-1,1)
3205  q2(-1,0) = q2(0,2); q2(0,0) = q2(0,1); q2(0,-1) = q2(-1,1)
3206  q3(-1,0) = q3(0,2); q3(0,0) = q3(0,1); q3(0,-1) = q3(-1,1)
3207  endif
3208  if ( se_corner ) then
3209  q1(npx+1,0) = q1(npx,2); q1(npx,0) = q1(npx,1); q1(npx,-1) = q1(npx+1,1)
3210  q2(npx+1,0) = q2(npx,2); q2(npx,0) = q2(npx,1); q2(npx,-1) = q2(npx+1,1)
3211  q3(npx+1,0) = q3(npx,2); q3(npx,0) = q3(npx,1); q3(npx,-1) = q3(npx+1,1)
3212  endif
3213  if ( ne_corner ) then
3214  q1(npx,npy) = q1(npx,npy-1); q1(npx+1,npy) = q1(npx,npy-2); q1(npx,npy+1) = q1(npx+1,npy-1)
3215  q2(npx,npy) = q2(npx,npy-1); q2(npx+1,npy) = q2(npx,npy-2); q2(npx,npy+1) = q2(npx+1,npy-1)
3216  q3(npx,npy) = q3(npx,npy-1); q3(npx+1,npy) = q3(npx,npy-2); q3(npx,npy+1) = q3(npx+1,npy-1)
3217  endif
3218  if ( nw_corner ) then
3219  q1(0,npy) = q1(0,npy-1); q1(-1,npy) = q1(0,npy-2); q1(0,npy+1) = q1(-1,npy-1)
3220  q2(0,npy) = q2(0,npy-1); q2(-1,npy) = q2(0,npy-2); q2(0,npy+1) = q2(-1,npy-1)
3221  q3(0,npy) = q3(0,npy-1); q3(-1,npy) = q3(0,npy-2); q3(0,npy+1) = q3(-1,npy-1)
3222  endif
3223 
3224  case(2)
3225  if ( sw_corner ) then
3226  q1(0,0) = q1(1,0); q1(0,-1) = q1(2,0); q1(-1,0) = q1(1,-1)
3227  q2(0,0) = q2(1,0); q2(0,-1) = q2(2,0); q2(-1,0) = q2(1,-1)
3228  q3(0,0) = q3(1,0); q3(0,-1) = q3(2,0); q3(-1,0) = q3(1,-1)
3229  endif
3230  if ( se_corner ) then
3231  q1(npx,0) = q1(npx-1,0); q1(npx,-1) = q1(npx-2,0); q1(npx+1,0) = q1(npx-1,-1)
3232  q2(npx,0) = q2(npx-1,0); q2(npx,-1) = q2(npx-2,0); q2(npx+1,0) = q2(npx-1,-1)
3233  q3(npx,0) = q3(npx-1,0); q3(npx,-1) = q3(npx-2,0); q3(npx+1,0) = q3(npx-1,-1)
3234  endif
3235  if ( ne_corner ) then
3236  q1(npx,npy) = q1(npx-1,npy); q1(npx,npy+1) = q1(npx-2,npy); q1(npx+1,npy) = q1(npx-1,npy+1)
3237  q2(npx,npy) = q2(npx-1,npy); q2(npx,npy+1) = q2(npx-2,npy); q2(npx+1,npy) = q2(npx-1,npy+1)
3238  q3(npx,npy) = q3(npx-1,npy); q3(npx,npy+1) = q3(npx-2,npy); q3(npx+1,npy) = q3(npx-1,npy+1)
3239  endif
3240  if ( nw_corner ) then
3241  q1(0,npy) = q1(1,npy); q1(0,npy+1) = q1(2,npy); q1(-1,npy) = q1(1,npy+1)
3242  q2(0,npy) = q2(1,npy); q2(0,npy+1) = q2(2,npy); q2(-1,npy) = q2(1,npy+1)
3243  q3(0,npy) = q3(1,npy); q3(0,npy+1) = q3(2,npy); q3(-1,npy) = q3(1,npy+1)
3244  endif
3245 
3246  end select
3247  end subroutine fill3_4corners
3248 
3250  subroutine fill2_4corners(q1, q2, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
3251  type(fv_grid_bounds_type), intent(IN) :: bd
3252  integer, intent(in):: dir
3253  real, intent(inout):: q1(bd%isd:bd%ied,bd%jsd:bd%jed)
3254  real, intent(inout):: q2(bd%isd:bd%ied,bd%jsd:bd%jed)
3255  logical, intent(IN) :: sw_corner, se_corner, ne_corner, nw_corner
3256  integer, intent(IN) :: npx, npy
3257 
3258  integer :: is, ie, js, je
3259  integer :: isd, ied, jsd, jed
3260 
3261  is = bd%is
3262  ie = bd%ie
3263  js = bd%js
3264  je = bd%je
3265  isd = bd%isd
3266  ied = bd%ied
3267  jsd = bd%jsd
3268  jed = bd%jed
3269 
3270  select case(dir)
3271  case(1)
3272  if ( sw_corner ) then
3273  q1(-1,0) = q1(0,2); q1(0,0) = q1(0,1)
3274  q2(-1,0) = q2(0,2); q2(0,0) = q2(0,1)
3275  endif
3276  if ( se_corner ) then
3277  q1(npx+1,0) = q1(npx,2); q1(npx,0) = q1(npx,1)
3278  q2(npx+1,0) = q2(npx,2); q2(npx,0) = q2(npx,1)
3279  endif
3280  if ( nw_corner ) then
3281  q1(0,npy) = q1(0,npy-1); q1(-1,npy) = q1(0,npy-2)
3282  q2(0,npy) = q2(0,npy-1); q2(-1,npy) = q2(0,npy-2)
3283  endif
3284  if ( ne_corner ) then
3285  q1(npx,npy) = q1(npx,npy-1); q1(npx+1,npy) = q1(npx,npy-2)
3286  q2(npx,npy) = q2(npx,npy-1); q2(npx+1,npy) = q2(npx,npy-2)
3287  endif
3288 
3289  case(2)
3290  if ( sw_corner ) then
3291  q1(0,0) = q1(1,0); q1(0,-1) = q1(2,0)
3292  q2(0,0) = q2(1,0); q2(0,-1) = q2(2,0)
3293  endif
3294  if ( se_corner ) then
3295  q1(npx,0) = q1(npx-1,0); q1(npx,-1) = q1(npx-2,0)
3296  q2(npx,0) = q2(npx-1,0); q2(npx,-1) = q2(npx-2,0)
3297  endif
3298  if ( nw_corner ) then
3299  q1(0,npy) = q1(1,npy); q1(0,npy+1) = q1(2,npy)
3300  q2(0,npy) = q2(1,npy); q2(0,npy+1) = q2(2,npy)
3301  endif
3302  if ( ne_corner ) then
3303  q1(npx,npy) = q1(npx-1,npy); q1(npx,npy+1) = q1(npx-2,npy)
3304  q2(npx,npy) = q2(npx-1,npy); q2(npx,npy+1) = q2(npx-2,npy)
3305  endif
3306 
3307  end select
3308 
3309  end subroutine fill2_4corners
3310 
3312  subroutine fill_4corners(q, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
3313  type(fv_grid_bounds_type), intent(IN) :: bd
3314  integer, intent(in):: dir ! 1: x-dir; 2: y-dir
3315  real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed)
3316  logical, intent(IN) :: sw_corner, se_corner, ne_corner, nw_corner
3317  integer, intent(IN) :: npx, npy
3318 
3319  integer :: is, ie, js, je
3320  integer :: isd, ied, jsd, jed
3321 
3322  is = bd%is
3323  ie = bd%ie
3324  js = bd%js
3325  je = bd%je
3326  isd = bd%isd
3327  ied = bd%ied
3328  jsd = bd%jsd
3329  jed = bd%jed
3330 
3331  select case(dir)
3332  case(1)
3333  if ( sw_corner ) then
3334  q(-1,0) = q(0,2)
3335  q( 0,0) = q(0,1)
3336  endif
3337  if ( se_corner ) then
3338  q(npx+1,0) = q(npx,2)
3339  q(npx, 0) = q(npx,1)
3340  endif
3341  if ( nw_corner ) then
3342  q( 0,npy) = q(0,npy-1)
3343  q(-1,npy) = q(0,npy-2)
3344  endif
3345  if ( ne_corner ) then
3346  q(npx, npy) = q(npx,npy-1)
3347  q(npx+1,npy) = q(npx,npy-2)
3348  endif
3349 
3350  case(2)
3351  if ( sw_corner ) then
3352  q(0, 0) = q(1,0)
3353  q(0,-1) = q(2,0)
3354  endif
3355  if ( se_corner ) then
3356  q(npx, 0) = q(npx-1,0)
3357  q(npx,-1) = q(npx-2,0)
3358  endif
3359  if ( nw_corner ) then
3360  q(0,npy ) = q(1,npy)
3361  q(0,npy+1) = q(2,npy)
3362  endif
3363  if ( ne_corner ) then
3364  q(npx,npy ) = q(npx-1,npy)
3365  q(npx,npy+1) = q(npx-2,npy)
3366  endif
3367 
3368  end select
3369 
3370  end subroutine fill_4corners
3371 
3372  end module sw_core_mod
sw_core_mod::s15
real, parameter s15
Definition: sw_core.F90:61
sw_core_mod::p1
real, parameter p1
0.58333333
Definition: sw_core.F90:71
sw_core_mod::divergence_corner_nest
subroutine, public divergence_corner_nest(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
Definition: sw_core.F90:1803
sw_core_mod::s11
real, parameter s11
Definition: sw_core.F90:61
tp_core_mod::fv_tp_2d
subroutine, public fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, gridstruct, bd, ra_x, ra_y, lim_fac, mfx, mfy, mass, nord, damp_c)
The subroutine &#39;fv_tp_2d&#39; contains the FV advection scheme .
Definition: tp_core.F90:109
sw_core_mod::fill2_4corners
subroutine fill2_4corners(q1, q2, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
The subroutine &#39; fill2_4corners&#39; fills the 4 corners of the scalar fileds only as needed by &#39;c_core&#39;...
Definition: sw_core.F90:3251
fv_mp_mod
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
sw_core_mod::t14
real, parameter t14
Definition: sw_core.F90:60
sw_core_mod::r3
real, parameter r3
Definition: sw_core.F90:59
fv_arrays_mod::fv_grid_type
The type &#39;fv_grid_type&#39; is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
Definition: fv_arrays.F90:123
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1187
test_cases_mod
Definition: test_cases.F90:22
sw_core_mod::c_sw
subroutine, public c_sw(delpc, delp, ptc, pt, u, v, w, uc, vc, ua, va, wc, ut, vt, divg_d, nord, dt2, hydrostatic, dord4, bd, gridstruct, flagstruct)
The subroutine &#39;c_sw&#39; performs a half-timestep advance of the C-grid winds.
Definition: sw_core.F90:105
test_cases_mod::test_case
integer test_case
Definition: test_cases.F90:180
sw_core_mod::b2
real, parameter b2
Definition: sw_core.F90:90
sw_core_mod
The module &#39;sw_core&#39; advances the forward step of the Lagrangian dynamics as described by ...
Definition: sw_core.F90:26
sw_core_mod::b5
real, parameter b5
Definition: sw_core.F90:93
sw_core_mod::b4
real, parameter b4
Definition: sw_core.F90:92
sw_core_mod::t11
real, parameter t11
Definition: sw_core.F90:60
sw_core_mod::divergence_corner
subroutine, public divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
The subroutine &#39;divergence_corner&#39; computes the cell-mean divergence on the "dual grid"...
Definition: sw_core.F90:1696
sw_core_mod::del6_vt_flux
subroutine, public del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd)
The subroutine &#39;del6_vt_flux&#39; applies 2nd, 4th, or 6th-order damping to fluxes ("vorticity damping") ...
Definition: sw_core.F90:1581
sw_core_mod::c3
real, parameter c3
Definition: sw_core.F90:82
sw_core_mod::edge_interpolate4
real function edge_interpolate4(ua, dxa)
Definition: sw_core.F90:3166
sw_core_mod::p2
real, parameter p2
Definition: sw_core.F90:72
sw_core_mod::d2a2c_vect
subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, bd, npx, npy, bounded_domain, grid_type)
Definition: sw_core.F90:2825
tp_core_mod
The module &#39;tp_core&#39; is a collection of routines to support FV transport.
Definition: tp_core.F90:24
a2b_edge_mod
The module &#39;a2b_edge&#39; performs FV-consistent interpolation of pressure to corners.
Definition: a2b_edge.F90:24
sw_core_mod::c1
real, parameter c1
Definition: sw_core.F90:80
sw_core_mod::xtp_u
subroutine xtp_u(is, ie, js, je, isd, ied, jsd, jed, c, u, v, flux, iord, dx, rdx, npx, npy, grid_type, lim_fac, bounded_domain)
Definition: sw_core.F90:2002
sw_core_mod::fill_4corners
subroutine, public fill_4corners(q, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
The subroutine &#39;fill_4corners&#39; fills the 4 corners of the scalar fields only as needed by c_core...
Definition: sw_core.F90:3313
sw_core_mod::a1
real, parameter a1
Definition: sw_core.F90:76
fv_arrays_mod
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
sw_core_mod::big_number
real, parameter big_number
Definition: sw_core.F90:66
tp_core_mod::pert_ppm
subroutine, public pert_ppm(im, a0, al, ar, iv)
Definition: tp_core.F90:1030
sw_core_mod::ytp_v
subroutine ytp_v(is, ie, js, je, isd, ied, jsd, jed, c, u, v, flux, jord, dy, rdy, npx, npy, grid_type, lim_fac, bounded_domain)
Definition: sw_core.F90:2360
sw_core_mod::b1
real, parameter b1
Definition: sw_core.F90:89
sw_core_mod::a2
real, parameter a2
Definition: sw_core.F90:77
sw_core_mod::smag_corner
subroutine smag_corner(dt, u, v, ua, va, smag_c, bd, npx, npy, gridstruct, ng)
The subroutine &#39;smag_corner&#39; computes Smagorinsky damping.
Definition: sw_core.F90:1912
fv_mp_mod::ng
integer, parameter, public ng
Definition: fv_mp_mod.F90:2337
fv_arrays_mod::fv_flags_type
Definition: fv_arrays.F90:248
sw_core_mod::near_zero
real, parameter near_zero
for KE limiter
Definition: sw_core.F90:62
sw_core_mod::t12
real, parameter t12
Definition: sw_core.F90:60
sw_core_mod::d_sw
subroutine, public d_sw(delpc, delp, ptc, pt, u, v, w, uc, vc, ua, va, divg_d, xflux, yflux, cx, cy, crx_adv, cry_adv, xfx_adv, yfx_adv, q_con, z_rat, kgb, heat_source, diss_est, zvir, sphum, nq, q, k, km, inline_q, dt, hord_tr, hord_mt, hord_vt, hord_tm, hord_dp, nord, nord_v, nord_w, nord_t, dddmp, d2_bg, d4_bg, damp_v, damp_w, damp_t, d_con, hydrostatic, gridstruct, flagstruct, bd)
The subroutine &#39;d_sw&#39; peforms a full-timestep advance of the D-grid winds and other prognostic varaia...
Definition: sw_core.F90:526
sw_core_mod::c2
real, parameter c2
Definition: sw_core.F90:81
tp_core_mod::copy_corners
subroutine, public copy_corners(q, npx, npy, dir, bounded_domain, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core.F90:248
sw_core_mod::s13
real, parameter s13
Definition: sw_core.F90:61
a2b_edge_mod::a2b_ord4
subroutine, public a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
Definition: a2b_edge.F90:70
sw_core_mod::s14
real, parameter s14
Definition: sw_core.F90:61
sw_core_mod::t15
real, parameter t15
Definition: sw_core.F90:60
sw_core_mod::fill3_4corners
subroutine fill3_4corners(q1, q2, q3, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
The subroutine &#39;fill3_4corners&#39; fills the 4 corners of the scalar fileds only as needed by &#39;c_core&#39;...
Definition: sw_core.F90:3180
sw_core_mod::t13
real, parameter t13
Definition: sw_core.F90:60
sw_core_mod::b3
real, parameter b3
Definition: sw_core.F90:91