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