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