FV3 Bundle
tp_core_nlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
20 module tp_core_nlm_mod
21 !BOP
22 !
23 ! !MODULE: tp_core --- A collection of routines to support FV transport
24 !
25  use fv_mp_nlm_mod, only: ng
26  use fv_grid_utils_nlm_mod, only: big_number
27  use fv_arrays_nlm_mod, only: fv_grid_type, fv_grid_bounds_type
28 
29  implicit none
30 
31  private
32  public fv_tp_2d, pert_ppm, copy_corners
33 
34  real, parameter:: ppm_fac = 1.5 ! nonlinear scheme limiter: between 1 and 2
35  real, parameter:: r3 = 1./3.
36  real, parameter:: near_zero = 1.e-25
37  real, parameter:: ppm_limiter = 2.0
38 
39 #ifdef WAVE_FORM
40 ! Suresh & Huynh scheme 2.2 (purtabation form)
41 ! The wave-form is more diffusive than scheme 2.1
42  real, parameter:: b1 = 0.0375
43  real, parameter:: b2 = -7./30.
44  real, parameter:: b3 = -23./120.
45  real, parameter:: b4 = 13./30.
46  real, parameter:: b5 = -11./240.
47 #else
48 ! scheme 2.1: perturbation form
49  real, parameter:: b1 = 1./30.
50  real, parameter:: b2 = -13./60.
51  real, parameter:: b3 = -13./60.
52  real, parameter:: b4 = 0.45
53  real, parameter:: b5 = -0.05
54 #endif
55  real, parameter:: t11 = 27./28., t12 = -13./28., t13=3./7.
56  real, parameter:: s11 = 11./14., s14 = 4./7., s15=3./14.
57 !----------------------------------------------------
58 ! volume-conserving cubic with 2nd drv=0 at end point:
59 !----------------------------------------------------
60 ! Non-monotonic
61  real, parameter:: c1 = -2./14.
62  real, parameter:: c2 = 11./14.
63  real, parameter:: c3 = 5./14.
64 !----------------------
65 ! PPM volume mean form:
66 !----------------------
67  real, parameter:: p1 = 7./12. ! 0.58333333
68  real, parameter:: p2 = -1./12.
69 ! q(i+0.5) = p1*(q(i-1)+q(i)) + p2*(q(i-2)+q(i+1))
70 ! integer:: is, ie, js, je, isd, ied, jsd, jed
71 
72 !
73 !EOP
74 !-----------------------------------------------------------------------
75 
76 contains
77 
78  subroutine fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, &
79  gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
80  type(fv_grid_bounds_type), intent(IN) :: bd
81  integer, intent(in):: npx, npy
82  integer, intent(in)::hord
83 
84  real, intent(in):: crx(bd%is:bd%ie+1,bd%jsd:bd%jed) !
85  real, intent(in):: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed) !
86  real, intent(in):: cry(bd%isd:bd%ied,bd%js:bd%je+1 ) !
87  real, intent(in):: yfx(bd%isd:bd%ied,bd%js:bd%je+1 ) !
88  real, intent(in):: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
89  real, intent(in):: ra_y(bd%isd:bd%ied,bd%js:bd%je)
90  real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed) ! transported scalar
91  real, intent(out)::fx(bd%is:bd%ie+1 ,bd%js:bd%je) ! Flux in x ( E )
92  real, intent(out)::fy(bd%is:bd%ie, bd%js:bd%je+1 ) ! Flux in y ( N )
93 
94  type(fv_grid_type), intent(IN), target :: gridstruct
95 ! optional Arguments:
96  real, OPTIONAL, intent(in):: mfx(bd%is:bd%ie+1,bd%js:bd%je ) ! Mass Flux X-Dir
97  real, OPTIONAL, intent(in):: mfy(bd%is:bd%ie ,bd%js:bd%je+1) ! Mass Flux Y-Dir
98  real, OPTIONAL, intent(in):: mass(bd%isd:bd%ied,bd%jsd:bd%jed)
99  real, OPTIONAL, intent(in):: damp_c
100  integer, OPTIONAL, intent(in):: nord
101 ! Local:
102  integer ord_ou, ord_in
103  real q_i(bd%isd:bd%ied,bd%js:bd%je)
104  real q_j(bd%is:bd%ie,bd%jsd:bd%jed)
105  real fx2(bd%is:bd%ie+1,bd%jsd:bd%jed)
106  real fy2(bd%isd:bd%ied,bd%js:bd%je+1)
107  real fyy(bd%isd:bd%ied,bd%js:bd%je+1)
108  real fx1(bd%is:bd%ie+1)
109  real damp
110  integer i, j
111 
112  integer:: is, ie, js, je, isd, ied, jsd, jed
113 
114  is = bd%is
115  ie = bd%ie
116  js = bd%js
117  je = bd%je
118  isd = bd%isd
119  ied = bd%ied
120  jsd = bd%jsd
121  jed = bd%jed
122 
123  if ( hord == 10 ) then
124  ord_in = 8
125  else
126  ord_in = hord
127  endif
128  ord_ou = hord
129 
130  if (.not. gridstruct%nested) call copy_corners(q, npx, npy, 2, gridstruct%nested, bd, &
131  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
132 
133  call yppm(fy2, q, cry, ord_in, isd,ied,isd,ied, js,je,jsd,jed, npx,npy, gridstruct%dya, gridstruct%nested, gridstruct%grid_type)
134 
135  do j=js,je+1
136  do i=isd,ied
137  fyy(i,j) = yfx(i,j) * fy2(i,j)
138  enddo
139  enddo
140  do j=js,je
141  do i=isd,ied
142  q_i(i,j) = (q(i,j)*gridstruct%area(i,j) + fyy(i,j)-fyy(i,j+1))/ra_y(i,j)
143  enddo
144  enddo
145 
146  call xppm(fx, q_i, crx(is,js), ord_ou, is,ie,isd,ied, js,je,jsd,jed, npx,npy, gridstruct%dxa, gridstruct%nested, gridstruct%grid_type)
147 
148  if (.not. gridstruct%nested) call copy_corners(q, npx, npy, 1, gridstruct%nested, bd, &
149  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
150 
151  call xppm(fx2, q, crx, ord_in, is,ie,isd,ied, jsd,jed,jsd,jed, npx,npy, gridstruct%dxa, gridstruct%nested, gridstruct%grid_type)
152 
153  do j=jsd,jed
154  do i=is,ie+1
155  fx1(i) = xfx(i,j) * fx2(i,j)
156  enddo
157  do i=is,ie
158  q_j(i,j) = (q(i,j)*gridstruct%area(i,j) + fx1(i)-fx1(i+1))/ra_x(i,j)
159  enddo
160  enddo
161 
162  call yppm(fy, q_j, cry, ord_ou, is,ie,isd,ied, js,je,jsd,jed, npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%grid_type)
163 
164 !----------------
165 ! Flux averaging:
166 !----------------
167 
168  if ( present(mfx) .and. present(mfy) ) then
169 !---------------------------------
170 ! For transport of pt and tracers
171 !---------------------------------
172  do j=js,je
173  do i=is,ie+1
174  fx(i,j) = 0.5*(fx(i,j) + fx2(i,j)) * mfx(i,j)
175  enddo
176  enddo
177  do j=js,je+1
178  do i=is,ie
179  fy(i,j) = 0.5*(fy(i,j) + fy2(i,j)) * mfy(i,j)
180  enddo
181  enddo
182  if ( present(nord) .and. present(damp_c) .and. present(mass) ) then
183  if ( damp_c > 1.e-4 ) then
184  damp = (damp_c * gridstruct%da_min)**(nord+1)
185  call deln_flux(nord, is,ie,js,je, npx, npy, damp, q, fx, fy, gridstruct, bd, mass )
186  endif
187  endif
188  else
189 !---------------------------------
190 ! For transport of delp, vorticity
191 !---------------------------------
192  do j=js,je
193  do i=is,ie+1
194  fx(i,j) = 0.5*(fx(i,j) + fx2(i,j)) * xfx(i,j)
195  enddo
196  enddo
197  do j=js,je+1
198  do i=is,ie
199  fy(i,j) = 0.5*(fy(i,j) + fy2(i,j)) * yfx(i,j)
200  enddo
201  enddo
202  if ( present(nord) .and. present(damp_c) ) then
203  if ( damp_c > 1.e-4 ) then
204  damp = (damp_c * gridstruct%da_min)**(nord+1)
205  call deln_flux(nord, is,ie,js,je, npx, npy, damp, q, fx, fy, gridstruct, bd)
206  endif
207  endif
208  endif
209 
210  end subroutine fv_tp_2d
211 
212  !Weird arguments are because this routine is called in a lot of
213  !places outside of tp_core, sometimes very deeply nested in the call tree.
214  subroutine copy_corners(q, npx, npy, dir, nested, bd, &
215  sw_corner, se_corner, nw_corner, ne_corner)
216  type(fv_grid_bounds_type), intent(IN) :: bd
217  integer, intent(in):: npx, npy, dir
218  real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed)
219  logical, intent(IN) :: nested, sw_corner, se_corner, nw_corner, ne_corner
220  integer i,j
221 
222  if (nested) return
223 
224  if ( dir == 1 ) then
225 ! XDir:
226  if ( sw_corner ) then
227  do j=1-ng,0
228  do i=1-ng,0
229  q(i,j) = q(j,1-i)
230  enddo
231  enddo
232  endif
233  if ( se_corner ) then
234  do j=1-ng,0
235  do i=npx,npx+ng-1
236  q(i,j) = q(npy-j,i-npx+1)
237  enddo
238  enddo
239  endif
240  if ( ne_corner ) then
241  do j=npy,npy+ng-1
242  do i=npx,npx+ng-1
243  q(i,j) = q(j,2*npx-1-i)
244  enddo
245  enddo
246  endif
247  if ( nw_corner ) then
248  do j=npy,npy+ng-1
249  do i=1-ng,0
250  q(i,j) = q(npy-j,i-1+npx)
251  enddo
252  enddo
253  endif
254 
255  elseif ( dir == 2 ) then
256 ! YDir:
257 
258  if ( sw_corner ) then
259  do j=1-ng,0
260  do i=1-ng,0
261  q(i,j) = q(1-j,i)
262  enddo
263  enddo
264  endif
265  if ( se_corner ) then
266  do j=1-ng,0
267  do i=npx,npx+ng-1
268  q(i,j) = q(npy+j-1,npx-i)
269  enddo
270  enddo
271  endif
272  if ( ne_corner ) then
273  do j=npy,npy+ng-1
274  do i=npx,npx+ng-1
275  q(i,j) = q(2*npy-1-j,i)
276  enddo
277  enddo
278  endif
279  if ( nw_corner ) then
280  do j=npy,npy+ng-1
281  do i=1-ng,0
282  q(i,j) = q(j+1-npx,npy-i)
283  enddo
284  enddo
285  endif
286 
287  endif
288 
289  end subroutine copy_corners
290 
291  subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, dxa, nested, grid_type)
292  integer, INTENT(IN) :: is, ie, isd, ied, jsd, jed
293  integer, INTENT(IN) :: jfirst, jlast ! compute domain
294  integer, INTENT(IN) :: iord
295  integer, INTENT(IN) :: npx, npy
296  real , INTENT(IN) :: q(isd:ied,jfirst:jlast)
297  real , INTENT(IN) :: c(is:ie+1,jfirst:jlast) ! Courant N (like FLUX)
298  real , intent(IN) :: dxa(isd:ied,jsd:jed)
299  logical, intent(IN) :: nested
300  integer, intent(IN) :: grid_type
301 ! !OUTPUT PARAMETERS:
302  real , INTENT(OUT) :: flux(is:ie+1,jfirst:jlast) ! Flux
303 ! Local
304  real, dimension(is-1:ie+1):: bl, br, b0
305  real:: q1(isd:ied)
306  real, dimension(is:ie+1):: fx0, fx1
307  logical, dimension(is-1:ie+1):: smt5, smt6
308  real al(is-1:ie+2)
309  real dm(is-2:ie+2)
310  real dq(is-3:ie+2)
311  integer:: i, j, ie3, is1, ie1
312  real:: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
313 
314  if ( .not. nested .and. grid_type<3 ) then
315  is1 = max(3,is-1); ie3 = min(npx-2,ie+2)
316  ie1 = min(npx-3,ie+1)
317  else
318  is1 = is-1; ie3 = ie+2
319  ie1 = ie+1
320  end if
321 
322  do 666 j=jfirst,jlast
323 
324  do i=isd, ied
325  q1(i) = q(i,j)
326  enddo
327 
328  if ( iord < 8 ) then
329 ! ord = 2: perfectly linear ppm scheme
330 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
331 
332  do i=is1, ie3
333  al(i) = p1*(q1(i-1)+q1(i)) + p2*(q1(i-2)+q1(i+1))
334  enddo
335  if ( iord==7 ) then
336  do i=is1, ie3
337  if ( al(i)<0. ) al(i) = 0.5*(q1(i-1)+q1(i))
338  enddo
339  endif
340 
341  if ( .not.nested .and. grid_type<3 ) then
342  if ( is==1 ) then
343  al(0) = c1*q1(-2) + c2*q1(-1) + c3*q1(0)
344  al(1) = 0.5*(((2.*dxa(0,j)+dxa(-1,j))*q1(0)-dxa(0,j)*q1(-1))/(dxa(-1,j)+dxa(0,j)) &
345  + ((2.*dxa(1,j)+dxa( 2,j))*q1(1)-dxa(1,j)*q1( 2))/(dxa(1, j)+dxa(2,j)))
346  al(2) = c3*q1(1) + c2*q1(2) +c1*q1(3)
347  if(iord==7) then
348  al(0) = max(0., al(0))
349  al(1) = max(0., al(1))
350  al(2) = max(0., al(2))
351  endif
352  endif
353  if ( (ie+1)==npx ) then
354  al(npx-1) = c1*q1(npx-3) + c2*q1(npx-2) + c3*q1(npx-1)
355  al(npx) = 0.5*(((2.*dxa(npx-1,j)+dxa(npx-2,j))*q1(npx-1)-dxa(npx-1,j)*q1(npx-2))/(dxa(npx-2,j)+dxa(npx-1,j)) &
356  + ((2.*dxa(npx, j)+dxa(npx+1,j))*q1(npx )-dxa(npx, j)*q1(npx+1))/(dxa(npx, j)+dxa(npx+1,j)))
357  al(npx+1) = c3*q1(npx) + c2*q1(npx+1) + c1*q1(npx+2)
358  if(iord==7) then
359  al(npx-1) = max(0., al(npx-1))
360  al(npx ) = max(0., al(npx ))
361  al(npx+1) = max(0., al(npx+1))
362  endif
363  endif
364  endif
365 
366  if ( iord==2 ) then ! perfectly linear scheme
367 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7
368 
369 !DEC$ VECTOR ALWAYS
370  do i=is,ie+1
371  xt = c(i,j)
372  if ( xt > 0. ) then
373  qtmp = q1(i-1)
374  flux(i,j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-(qtmp+qtmp)))
375  else
376  qtmp = q1(i)
377  flux(i,j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-(qtmp+qtmp)))
378  endif
379 ! x0 = sign(dim(xt, 0.), 1.)
380 ! x1 = sign(dim(0., xt), 1.)
381 ! flux(i,j) = x0*(q1(i-1)+(1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-(qtmp+qtmp)))) &
382 ! + x1*(q1(i) +(1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-(qtmp+qtmp))))
383  enddo
384 
385  elseif ( iord==3 ) then
386  do i=is-1,ie+1
387  bl(i) = al(i) - q1(i)
388  br(i) = al(i+1) - q1(i)
389  b0(i) = bl(i) + br(i)
390  x0 = abs(b0(i))
391  xt = abs(bl(i)-br(i))
392  smt5(i) = x0 < xt
393  smt6(i) = 3.*x0 < xt
394  enddo
395  do i=is,ie+1
396  fx1(i) = 0.
397  enddo
398  do i=is,ie+1
399  xt = c(i,j)
400  if ( xt > 0. ) then
401  fx0(i) = q1(i-1)
402  if ( smt6(i-1).or.smt5(i) ) then
403  fx1(i) = br(i-1) - xt*b0(i-1)
404  elseif ( smt5(i-1) ) then ! 2nd order, piece-wise linear
405  fx1(i) = sign(min(abs(bl(i-1)),abs(br(i-1))), br(i-1))
406  endif
407  else
408  fx0(i) = q1(i)
409  if ( smt6(i).or.smt5(i-1) ) then
410  fx1(i) = bl(i) + xt*b0(i)
411  elseif ( smt5(i) ) then
412  fx1(i) = sign(min(abs(bl(i)), abs(br(i))), bl(i))
413  endif
414  endif
415  flux(i,j) = fx0(i) + (1.-abs(xt))*fx1(i)
416  enddo
417  elseif ( iord==4 ) then
418  do i=is-1,ie+1
419  bl(i) = al(i) - q1(i)
420  br(i) = al(i+1) - q1(i)
421  b0(i) = bl(i) + br(i)
422  x0 = abs(b0(i))
423  xt = abs(bl(i)-br(i))
424  smt5(i) = x0 < xt
425  smt6(i) = 3.*x0 < xt
426  enddo
427  do i=is,ie+1
428  fx1(i) = 0.
429  enddo
430 !DEC$ VECTOR ALWAYS
431  do i=is,ie+1
432  if ( c(i,j) > 0. ) then
433  fx0(i) = q1(i-1)
434  if ( smt6(i-1).or.smt5(i) ) fx1(i) = (1.-c(i,j))*(br(i-1) - c(i,j)*b0(i-1))
435  else
436  fx0(i) = q1(i)
437  if ( smt6(i).or.smt5(i-1) ) fx1(i) = (1.+c(i,j))*(bl(i) + c(i,j)*b0(i))
438  endif
439  flux(i,j) = fx0(i) + fx1(i)
440  enddo
441  else
442 ! iord = 5 & 6
443  if ( iord==5 ) then
444  do i=is-1,ie+1
445  bl(i) = al(i) - q1(i)
446  br(i) = al(i+1) - q1(i)
447  b0(i) = bl(i) + br(i)
448  smt5(i) = bl(i)*br(i) < 0.
449  enddo
450  else
451  do i=is-1,ie+1
452  bl(i) = al(i) - q1(i)
453  br(i) = al(i+1) - q1(i)
454  b0(i) = bl(i) + br(i)
455  smt5(i) = abs(3.*b0(i)) < abs(bl(i)-br(i))
456  enddo
457  endif
458 !DEC$ VECTOR ALWAYS
459  do i=is,ie+1
460  if ( c(i,j) > 0. ) then
461  fx1(i) = (1.-c(i,j))*(br(i-1) - c(i,j)*b0(i-1))
462  flux(i,j) = q1(i-1)
463  else
464  fx1(i) = (1.+c(i,j))*(bl(i) + c(i,j)*b0(i))
465  flux(i,j) = q1(i)
466  endif
467  if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx1(i)
468  enddo
469  endif
470  goto 666
471 
472  else
473 
474 ! Monotonic constraints:
475 ! ord = 8: PPM with Lin's PPM fast monotone constraint
476 ! ord = 10: PPM with Lin's modification of Huynh 2nd constraint
477 ! ord = 13: 10 plus positive definite constraint
478 
479  do i=is-2,ie+2
480  xt = 0.25*(q1(i+1) - q1(i-1))
481  dm(i) = sign(min(abs(xt), max(q1(i-1), q1(i), q1(i+1)) - q1(i), &
482  q1(i) - min(q1(i-1), q1(i), q1(i+1))), xt)
483  enddo
484  do i=is1,ie1+1
485  al(i) = 0.5*(q1(i-1)+q1(i)) + r3*(dm(i-1)-dm(i))
486  enddo
487 
488  if ( iord==8 ) then
489  do i=is1, ie1
490  xt = 2.*dm(i)
491  bl(i) = -sign(min(abs(xt), abs(al(i )-q1(i))), xt)
492  br(i) = sign(min(abs(xt), abs(al(i+1)-q1(i))), xt)
493  enddo
494  elseif ( iord==11 ) then
495 ! This is emulation of 2nd van Leer scheme using PPM codes
496  do i=is1, ie1
497  xt = ppm_fac*dm(i)
498  bl(i) = -sign(min(abs(xt), abs(al(i )-q1(i))), xt)
499  br(i) = sign(min(abs(xt), abs(al(i+1)-q1(i))), xt)
500  enddo
501  else
502  do i=is1-2, ie1+1
503  dq(i) = 2.*(q1(i+1) - q1(i))
504  enddo
505  do i=is1, ie1
506  bl(i) = al(i ) - q1(i)
507  br(i) = al(i+1) - q1(i)
508  if ( abs(dm(i-1))+abs(dm(i))+abs(dm(i+1)) < near_zero ) then
509  bl(i) = 0.
510  br(i) = 0.
511  elseif( abs(3.*(bl(i)+br(i))) > abs(bl(i)-br(i)) ) then
512  pmp_2 = dq(i-1)
513  lac_2 = pmp_2 - 0.75*dq(i-2)
514  br(i) = min( max(0., pmp_2, lac_2), max(br(i), min(0., pmp_2, lac_2)) )
515  pmp_1 = -dq(i)
516  lac_1 = pmp_1 + 0.75*dq(i+1)
517  bl(i) = min( max(0., pmp_1, lac_1), max(bl(i), min(0., pmp_1, lac_1)) )
518  endif
519  enddo
520  endif
521 ! Positive definite constraint:
522  if(iord==9 .or. iord==13) call pert_ppm(ie1-is1+1, q1(is1), bl(is1), br(is1), 0)
523 
524  if (.not. nested .and. grid_type<3) then
525  if ( is==1 ) then
526  bl(0) = s14*dm(-1) + s11*(q1(-1)-q1(0))
527 
528  xt = 0.5*(((2.*dxa(0,j)+dxa(-1,j))*q1(0)-dxa(0,j)*q1(-1))/(dxa(-1,j)+dxa(0,j)) &
529  + ((2.*dxa(1,j)+dxa( 2,j))*q1(1)-dxa(1,j)*q1( 2))/(dxa(1, j)+dxa(2,j)))
530 ! if ( iord==8 .or. iord==10 ) then
531  xt = max(xt, min(q1(-1),q1(0),q1(1),q1(2)))
532  xt = min(xt, max(q1(-1),q1(0),q1(1),q1(2)))
533 ! endif
534  br(0) = xt - q1(0)
535  bl(1) = xt - q1(1)
536  xt = s15*q1(1) + s11*q1(2) - s14*dm(2)
537  br(1) = xt - q1(1)
538  bl(2) = xt - q1(2)
539 
540  br(2) = al(3) - q1(2)
541  call pert_ppm(3, q1(0), bl(0), br(0), 1)
542  endif
543  if ( (ie+1)==npx ) then
544  bl(npx-2) = al(npx-2) - q1(npx-2)
545 
546  xt = s15*q1(npx-1) + s11*q1(npx-2) + s14*dm(npx-2)
547  br(npx-2) = xt - q1(npx-2)
548  bl(npx-1) = xt - q1(npx-1)
549 
550  xt = 0.5*(((2.*dxa(npx-1,j)+dxa(npx-2,j))*q1(npx-1)-dxa(npx-1,j)*q1(npx-2))/(dxa(npx-2,j)+dxa(npx-1,j)) &
551  + ((2.*dxa(npx, j)+dxa(npx+1,j))*q1(npx )-dxa(npx, j)*q1(npx+1))/(dxa(npx, j)+dxa(npx+1,j)))
552 ! if ( iord==8 .or. iord==10 ) then
553  xt = max(xt, min(q1(npx-2),q1(npx-1),q1(npx),q1(npx+1)))
554  xt = min(xt, max(q1(npx-2),q1(npx-1),q1(npx),q1(npx+1)))
555 ! endif
556  br(npx-1) = xt - q1(npx-1)
557  bl(npx ) = xt - q1(npx )
558 
559  br(npx) = s11*(q1(npx+1)-q1(npx)) - s14*dm(npx+1)
560  call pert_ppm(3, q1(npx-2), bl(npx-2), br(npx-2), 1)
561  endif
562  endif
563 
564  endif
565 
566  do i=is,ie+1
567  if( c(i,j)>0. ) then
568  flux(i,j) = q1(i-1) + (1.-c(i,j))*(br(i-1)-c(i,j)*(bl(i-1)+br(i-1)))
569  else
570  flux(i,j) = q1(i ) + (1.+c(i,j))*(bl(i )+c(i,j)*(bl(i)+br(i)))
571  endif
572  enddo
573 
574 666 continue
575 
576  end subroutine xppm
577 
578 
579  subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy, dya, nested, grid_type)
580  integer, INTENT(IN) :: ifirst,ilast ! Compute domain
581  integer, INTENT(IN) :: isd,ied, js,je,jsd,jed
582  integer, INTENT(IN) :: jord
583  integer, INTENT(IN) :: npx, npy
584  real , INTENT(IN) :: q(ifirst:ilast,jsd:jed)
585  real , intent(in) :: c(isd:ied,js:je+1 ) ! Courant number
586  real , INTENT(OUT):: flux(ifirst:ilast,js:je+1) ! Flux
587  real , intent(IN) :: dya(isd:ied,jsd:jed)
588  logical, intent(IN) :: nested
589  integer, intent(IN) :: grid_type
590 ! Local:
591  real:: dm(ifirst:ilast,js-2:je+2)
592  real:: al(ifirst:ilast,js-1:je+2)
593  real, dimension(ifirst:ilast,js-1:je+1):: bl, br, b0
594  real:: dq(ifirst:ilast,js-3:je+2)
595  real, dimension(ifirst:ilast):: fx0, fx1
596  logical, dimension(ifirst:ilast,js-1:je+1):: smt5, smt6
597  real:: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
598  integer:: i, j, js1, je3, je1
599 
600  if ( .not.nested .and. grid_type < 3 ) then
601 ! Cubed-sphere:
602  js1 = max(3,js-1); je3 = min(npy-2,je+2)
603  je1 = min(npy-3,je+1)
604  else
605 ! Nested grid OR Doubly periodic domain:
606  js1 = js-1; je3 = je+2
607  je1 = je+1
608  endif
609 
610 if ( jord < 8 ) then
611 
612  do j=js1, je3
613  do i=ifirst,ilast
614  al(i,j) = p1*(q(i,j-1)+q(i,j)) + p2*(q(i,j-2)+q(i,j+1))
615  enddo
616  enddo
617  if ( jord==7 ) then
618  do j=js1, je3
619  do i=ifirst,ilast
620  if ( al(i,j)<0. ) al(i,j) = 0.5*(q(i,j)+q(i,j+1))
621  enddo
622  enddo
623  endif
624 
625  if ( .not. nested .and. grid_type<3 ) then
626  if( js==1 ) then
627  do i=ifirst,ilast
628  al(i,0) = c1*q(i,-2) + c2*q(i,-1) + c3*q(i,0)
629  al(i,1) = 0.5*(((2.*dya(i,0)+dya(i,-1))*q(i,0)-dya(i,0)*q(i,-1))/(dya(i,-1)+dya(i,0)) &
630  + ((2.*dya(i,1)+dya(i,2))*q(i,1)-dya(i,1)*q(i,2))/(dya(i,1)+dya(i,2)))
631  al(i,2) = c3*q(i,1) + c2*q(i,2) + c1*q(i,3)
632  enddo
633  if ( jord==7 ) then
634  do i=ifirst,ilast
635  al(i,0) = max(0., al(i,0))
636  al(i,1) = max(0., al(i,1))
637  al(i,2) = max(0., al(i,2))
638  enddo
639  endif
640  endif
641  if( (je+1)==npy ) then
642  do i=ifirst,ilast
643  al(i,npy-1) = c1*q(i,npy-3) + c2*q(i,npy-2) + c3*q(i,npy-1)
644  al(i,npy) = 0.5*(((2.*dya(i,npy-1)+dya(i,npy-2))*q(i,npy-1)-dya(i,npy-1)*q(i,npy-2))/(dya(i,npy-2)+dya(i,npy-1)) &
645  + ((2.*dya(i,npy)+dya(i,npy+1))*q(i,npy)-dya(i,npy)*q(i,npy+1))/(dya(i,npy)+dya(i,npy+1)))
646  al(i,npy+1) = c3*q(i,npy) + c2*q(i,npy+1) + c1*q(i,npy+2)
647  enddo
648  if (jord==7 ) then
649  do i=ifirst,ilast
650  al(i,npy-1) = max(0., al(i,npy-1))
651  al(i,npy ) = max(0., al(i,npy ))
652  al(i,npy+1) = max(0., al(i,npy+1))
653  enddo
654  endif
655  endif
656  endif
657 
658  if ( jord==2 ) then ! Perfectly linear scheme
659 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7
660 
661  do j=js,je+1
662 !DEC$ VECTOR ALWAYS
663  do i=ifirst,ilast
664  xt = c(i,j)
665  if ( xt > 0. ) then
666  qtmp = q(i,j-1)
667  flux(i,j) = qtmp + (1.-xt)*(al(i,j)-qtmp-xt*(al(i,j-1)+al(i,j)-(qtmp+qtmp)))
668  else
669  qtmp = q(i,j)
670  flux(i,j) = qtmp + (1.+xt)*(al(i,j)-qtmp+xt*(al(i,j)+al(i,j+1)-(qtmp+qtmp)))
671  endif
672  enddo
673  enddo
674 
675  elseif ( jord==3 ) then
676  do j=js-1,je+1
677  do i=ifirst,ilast
678  bl(i,j) = al(i,j ) - q(i,j)
679  br(i,j) = al(i,j+1) - q(i,j)
680  b0(i,j) = bl(i,j) + br(i,j)
681  x0 = abs(b0(i,j))
682  xt = abs(bl(i,j)-br(i,j))
683  smt5(i,j) = x0 < xt
684  smt6(i,j) = 3.*x0 < xt
685  enddo
686  enddo
687  do j=js,je+1
688  do i=ifirst,ilast
689  fx1(i) = 0.
690  enddo
691  do i=ifirst,ilast
692  xt = c(i,j)
693  if ( xt > 0. ) then
694  fx0(i) = q(i,j-1)
695  if( smt6(i,j-1).or.smt5(i,j) ) then
696  fx1(i) = br(i,j-1) - xt*b0(i,j-1)
697  elseif ( smt5(i,j-1) ) then ! both up-downwind sides are noisy; 2nd order, piece-wise linear
698  fx1(i) = sign(min(abs(bl(i,j-1)),abs(br(i,j-1))),br(i,j-1))
699  endif
700  else
701  fx0(i) = q(i,j)
702  if( smt6(i,j).or.smt5(i,j-1) ) then
703  fx1(i) = bl(i,j) + xt*b0(i,j)
704  elseif ( smt5(i,j) ) then
705  fx1(i) = sign(min(abs(bl(i,j)),abs(br(i,j))), bl(i,j))
706  endif
707  endif
708  flux(i,j) = fx0(i) + (1.-abs(xt))*fx1(i)
709  enddo
710  enddo
711 
712  elseif ( jord==4 ) then
713  do j=js-1,je+1
714  do i=ifirst,ilast
715  bl(i,j) = al(i,j ) - q(i,j)
716  br(i,j) = al(i,j+1) - q(i,j)
717  b0(i,j) = bl(i,j) + br(i,j)
718  x0 = abs(b0(i,j))
719  xt = abs(bl(i,j)-br(i,j))
720  smt5(i,j) = x0 < xt
721  smt6(i,j) = 3.*x0 < xt
722  enddo
723  enddo
724  do j=js,je+1
725  do i=ifirst,ilast
726  fx1(i) = 0.
727  enddo
728 !DEC$ VECTOR ALWAYS
729  do i=ifirst,ilast
730  if ( c(i,j) > 0. ) then
731  fx0(i) = q(i,j-1)
732  if( smt6(i,j-1).or.smt5(i,j) ) fx1(i) = (1.-c(i,j))*(br(i,j-1) - c(i,j)*b0(i,j-1))
733  else
734  fx0(i) = q(i,j)
735  if( smt6(i,j).or.smt5(i,j-1) ) fx1(i) = (1.+c(i,j))*(bl(i,j) + c(i,j)*b0(i,j))
736  endif
737  flux(i,j) = fx0(i) + fx1(i)
738  enddo
739  enddo
740 
741  else ! jord=5,6,7
742  if ( jord==5 ) then
743  do j=js-1,je+1
744  do i=ifirst,ilast
745  bl(i,j) = al(i,j ) - q(i,j)
746  br(i,j) = al(i,j+1) - q(i,j)
747  b0(i,j) = bl(i,j) + br(i,j)
748  smt5(i,j) = bl(i,j)*br(i,j) < 0.
749  enddo
750  enddo
751  else
752  do j=js-1,je+1
753  do i=ifirst,ilast
754  bl(i,j) = al(i,j ) - q(i,j)
755  br(i,j) = al(i,j+1) - q(i,j)
756  b0(i,j) = bl(i,j) + br(i,j)
757  smt5(i,j) = abs(3.*b0(i,j)) < abs(bl(i,j)-br(i,j))
758  enddo
759  enddo
760  endif
761  do j=js,je+1
762 !DEC$ VECTOR ALWAYS
763  do i=ifirst,ilast
764  if ( c(i,j) > 0. ) then
765  fx1(i) = (1.-c(i,j))*(br(i,j-1) - c(i,j)*b0(i,j-1))
766  flux(i,j) = q(i,j-1)
767  else
768  fx1(i) = (1.+c(i,j))*(bl(i,j) + c(i,j)*b0(i,j))
769  flux(i,j) = q(i,j)
770  endif
771  if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx1(i)
772  enddo
773  enddo
774  endif
775  return
776 
777 else
778 ! Monotonic constraints:
779 ! ord = 8: PPM with Lin's PPM fast monotone constraint
780 ! ord > 8: PPM with Lin's modification of Huynh 2nd constraint
781 
782  do j=js-2,je+2
783  do i=ifirst,ilast
784  xt = 0.25*(q(i,j+1) - q(i,j-1))
785  dm(i,j) = sign(min(abs(xt), max(q(i,j-1), q(i,j), q(i,j+1)) - q(i,j), &
786  q(i,j) - min(q(i,j-1), q(i,j), q(i,j+1))), xt)
787  enddo
788  enddo
789  do j=js1,je1+1
790  do i=ifirst,ilast
791  al(i,j) = 0.5*(q(i,j-1)+q(i,j)) + r3*(dm(i,j-1) - dm(i,j))
792  enddo
793  enddo
794 
795  if ( jord==8 ) then
796  do j=js1,je1
797  do i=ifirst,ilast
798  xt = 2.*dm(i,j)
799  bl(i,j) = -sign(min(abs(xt), abs(al(i,j)-q(i,j))), xt)
800  br(i,j) = sign(min(abs(xt), abs(al(i,j+1)-q(i,j))), xt)
801  enddo
802  enddo
803  elseif ( jord==11 ) then
804  do j=js1,je1
805  do i=ifirst,ilast
806  xt = ppm_fac*dm(i,j)
807  bl(i,j) = -sign(min(abs(xt), abs(al(i,j)-q(i,j))), xt)
808  br(i,j) = sign(min(abs(xt), abs(al(i,j+1)-q(i,j))), xt)
809  enddo
810  enddo
811  else
812  do j=js1-2,je1+1
813  do i=ifirst,ilast
814  dq(i,j) = 2.*(q(i,j+1) - q(i,j))
815  enddo
816  enddo
817  do j=js1,je1
818  do i=ifirst,ilast
819  bl(i,j) = al(i,j ) - q(i,j)
820  br(i,j) = al(i,j+1) - q(i,j)
821  if ( abs(dm(i,j-1))+abs(dm(i,j))+abs(dm(i,j+1)) < near_zero ) then
822  bl(i,j) = 0.
823  br(i,j) = 0.
824  elseif( abs(3.*(bl(i,j)+br(i,j))) > abs(bl(i,j)-br(i,j)) ) then
825  pmp_2 = dq(i,j-1)
826  lac_2 = pmp_2 - 0.75*dq(i,j-2)
827  br(i,j) = min(max(0.,pmp_2,lac_2), max(br(i,j), min(0.,pmp_2,lac_2)))
828  pmp_1 = -dq(i,j)
829  lac_1 = pmp_1 + 0.75*dq(i,j+1)
830  bl(i,j) = min(max(0.,pmp_1,lac_1), max(bl(i,j), min(0.,pmp_1,lac_1)))
831  endif
832  enddo
833  enddo
834  endif
835  if ( jord==9 .or. jord==13 ) then
836 ! Positive definite constraint:
837  do j=js1,je1
838  call pert_ppm(ilast-ifirst+1, q(ifirst,j), bl(ifirst,j), br(ifirst,j), 0)
839  enddo
840  endif
841 
842  if (.not. nested .and. grid_type<3) then
843  if( js==1 ) then
844  do i=ifirst,ilast
845  bl(i,0) = s14*dm(i,-1) + s11*(q(i,-1)-q(i,0))
846 
847  xt = 0.5*(((2.*dya(i,0)+dya(i,-1))*q(i,0)-dya(i,0)*q(i,-1))/(dya(i,-1)+dya(i,0)) &
848  + ((2.*dya(i,1)+dya(i,2))*q(i,1)-dya(i,1)*q(i,2))/(dya(i,1)+dya(i,2)))
849 ! if ( jord==8 .or. jord==10 ) then
850  xt = max(xt, min(q(i,-1),q(i,0),q(i,1),q(i,2)))
851  xt = min(xt, max(q(i,-1),q(i,0),q(i,1),q(i,2)))
852 ! endif
853  br(i,0) = xt - q(i,0)
854  bl(i,1) = xt - q(i,1)
855 
856  xt = s15*q(i,1) + s11*q(i,2) - s14*dm(i,2)
857  br(i,1) = xt - q(i,1)
858  bl(i,2) = xt - q(i,2)
859 
860  br(i,2) = al(i,3) - q(i,2)
861  enddo
862  call pert_ppm(3*(ilast-ifirst+1), q(ifirst,0), bl(ifirst,0), br(ifirst,0), 1)
863  endif
864  if( (je+1)==npy ) then
865  do i=ifirst,ilast
866  bl(i,npy-2) = al(i,npy-2) - q(i,npy-2)
867 
868  xt = s15*q(i,npy-1) + s11*q(i,npy-2) + s14*dm(i,npy-2)
869  br(i,npy-2) = xt - q(i,npy-2)
870  bl(i,npy-1) = xt - q(i,npy-1)
871 
872  xt = 0.5*(((2.*dya(i,npy-1)+dya(i,npy-2))*q(i,npy-1)-dya(i,npy-1)*q(i,npy-2))/(dya(i,npy-2)+dya(i,npy-1)) &
873  + ((2.*dya(i,npy)+dya(i,npy+1))*q(i,npy)-dya(i,npy)*q(i,npy+1))/(dya(i,npy)+dya(i,npy+1)))
874 ! if ( jord==8 .or. jord==10 ) then
875  xt = max(xt, min(q(i,npy-2),q(i,npy-1),q(i,npy),q(i,npy+1)))
876  xt = min(xt, max(q(i,npy-2),q(i,npy-1),q(i,npy),q(i,npy+1)))
877 ! endif
878  br(i,npy-1) = xt - q(i,npy-1)
879  bl(i,npy ) = xt - q(i,npy)
880 
881  br(i,npy) = s11*(q(i,npy+1)-q(i,npy)) - s14*dm(i,npy+1)
882  enddo
883  call pert_ppm(3*(ilast-ifirst+1), q(ifirst,npy-2), bl(ifirst,npy-2), br(ifirst,npy-2), 1)
884  endif
885  end if
886 
887 endif
888 
889  do j=js,je+1
890  do i=ifirst,ilast
891  if( c(i,j)>0. ) then
892  flux(i,j) = q(i,j-1) + (1.-c(i,j))*(br(i,j-1)-c(i,j)*(bl(i,j-1)+br(i,j-1)))
893  else
894  flux(i,j) = q(i,j ) + (1.+c(i,j))*(bl(i,j )+c(i,j)*(bl(i,j)+br(i,j)))
895  endif
896  enddo
897  enddo
898 
899  end subroutine yppm
900 
901 
902 
903  subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
904  kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
905 !
906 ! !INPUT PARAMETERS:
907  integer, intent(in):: im, jm, km, nq
908  integer, intent(in):: ifirst, ilast
909  integer, intent(in):: jfirst, jlast
910  integer, intent(in):: kfirst, klast
911  integer, intent(in):: ng_e ! eastern zones to ghost
912  integer, intent(in):: ng_w ! western zones to ghost
913  integer, intent(in):: ng_s ! southern zones to ghost
914  integer, intent(in):: ng_n ! northern zones to ghost
915  real, intent(inout):: q_ghst(ifirst-ng_w:ilast+ng_e,jfirst-ng_s:jlast+ng_n,kfirst:klast,nq)
916  real, optional, intent(in):: q(ifirst:ilast,jfirst:jlast,kfirst:klast,nq)
917 !
918 ! !DESCRIPTION:
919 !
920 ! Ghost 4d east/west
921 !
922 ! !REVISION HISTORY:
923 ! 2005.08.22 Putman
924 !
925 !EOP
926 !------------------------------------------------------------------------------
927 !BOC
928  integer :: i,j,k,n
929 
930  if (present(q)) then
931  q_ghst(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq) = &
932  q(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq)
933  endif
934 
935 ! Assume Periodicity in X-dir and not overlapping
936  do n=1,nq
937  do k=kfirst,klast
938  do j=jfirst-ng_s,jlast+ng_n
939  do i=1, ng_w
940  q_ghst(ifirst-i,j,k,n) = q_ghst(ilast-i+1,j,k,n)
941  enddo
942  do i=1, ng_e
943  q_ghst(ilast+i,j,k,n) = q_ghst(ifirst+i-1,j,k,n)
944  enddo
945  enddo
946  enddo
947  enddo
948 
949  end subroutine mp_ghost_ew
950 
951 
952 
953  subroutine pert_ppm(im, a0, al, ar, iv)
954  integer, intent(in):: im
955  integer, intent(in):: iv
956  real, intent(in) :: a0(im)
957  real, intent(inout):: al(im), ar(im)
958 ! Local:
959  real a4, da1, da2, a6da, fmin
960  integer i
961  real, parameter:: r12 = 1./12.
962 
963 !-----------------------------------
964 ! Optimized PPM in perturbation form:
965 !-----------------------------------
966 
967  if ( iv==0 ) then
968 ! Positive definite constraint
969  do i=1,im
970  if ( a0(i) <= 0. ) then
971  al(i) = 0.
972  ar(i) = 0.
973  else
974  a4 = -3.*(ar(i) + al(i))
975  da1 = ar(i) - al(i)
976  if( abs(da1) < -a4 ) then
977  fmin = a0(i) + 0.25/a4*da1**2 + a4*r12
978  if( fmin < 0. ) then
979  if( ar(i)>0. .and. al(i)>0. ) then
980  ar(i) = 0.
981  al(i) = 0.
982  elseif( da1 > 0. ) then
983  ar(i) = -2.*al(i)
984  else
985  al(i) = -2.*ar(i)
986  endif
987  endif
988  endif
989  endif
990  enddo
991  else
992 ! Standard PPM constraint
993  do i=1,im
994  if ( al(i)*ar(i) < 0. ) then
995  da1 = al(i) - ar(i)
996  da2 = da1**2
997  a6da = 3.*(al(i)+ar(i))*da1
998 ! abs(a6da) > da2 --> 3.*abs(al+ar) > abs(al-ar)
999  if( a6da < -da2 ) then
1000  ar(i) = -2.*al(i)
1001  elseif( a6da > da2 ) then
1002  al(i) = -2.*ar(i)
1003  endif
1004  else
1005 ! effect of dm=0 included here
1006  al(i) = 0.
1007  ar(i) = 0.
1008  endif
1009  enddo
1010  endif
1011 
1012  end subroutine pert_ppm
1013 
1014 
1015  subroutine deln_flux(nord,is,ie,js,je, npx, npy, damp, q, fx, fy, gridstruct, bd, mass )
1016 ! Del-n damping for the cell-mean values (A grid)
1017 !------------------
1018 ! nord = 0: del-2
1019 ! nord = 1: del-4
1020 ! nord = 2: del-6
1021 ! nord = 3: del-8 --> requires more ghosting than current
1022 !------------------
1023  type(fv_grid_bounds_type), intent(IN) :: bd
1024  integer, intent(in):: nord ! del-n
1025  integer, intent(in):: is,ie,js,je, npx, npy
1026  real, intent(in):: damp
1027  real, intent(in):: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng) ! q ghosted on input
1028  type(fv_grid_type), intent(IN), target :: gridstruct
1029  real, optional, intent(in):: mass(bd%isd:bd%ied, bd%jsd:bd%jed) ! q ghosted on input
1030 ! diffusive fluxes:
1031  real, intent(inout):: fx(bd%is:bd%ie+1,bd%js:bd%je), fy(bd%is:bd%ie,bd%js:bd%je+1)
1032 ! local:
1033  real fx2(bd%isd:bd%ied+1,bd%jsd:bd%jed), fy2(bd%isd:bd%ied,bd%jsd:bd%jed+1)
1034  real d2(bd%isd:bd%ied,bd%jsd:bd%jed)
1035  real damp2
1036  integer i,j, n, nt, i1, i2, j1, j2
1037 
1038 #ifdef USE_SG
1039  real, pointer, dimension(:,:) :: dx, dy, rdxc, rdyc
1040  real, pointer, dimension(:,:,:) :: sin_sg
1041  dx => gridstruct%dx
1042  dy => gridstruct%dy
1043  rdxc => gridstruct%rdxc
1044  rdyc => gridstruct%rdyc
1045  sin_sg => gridstruct%sin_sg
1046 #endif
1047 
1048  i1 = is-1-nord; i2 = ie+1+nord
1049  j1 = js-1-nord; j2 = je+1+nord
1050 
1051  if ( .not. present(mass) ) then
1052  do j=j1, j2
1053  do i=i1,i2
1054  d2(i,j) = damp*q(i,j)
1055  enddo
1056  enddo
1057  else
1058  do j=j1, j2
1059  do i=i1,i2
1060  d2(i,j) = q(i,j)
1061  enddo
1062  enddo
1063  endif
1064 
1065  if( nord>0 ) call copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, &
1066  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1067 
1068  do j=js-nord,je+nord
1069  do i=is-nord,ie+nord+1
1070 #ifdef USE_SG
1071  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)
1072 #else
1073  fx2(i,j) = gridstruct%del6_v(i,j)*(d2(i-1,j)-d2(i,j))
1074 #endif
1075  enddo
1076  enddo
1077 
1078  if( nord>0 ) call copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, &
1079  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1080  do j=js-nord,je+nord+1
1081  do i=is-nord,ie+nord
1082 #ifdef USE_SG
1083  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)
1084 #else
1085  fy2(i,j) = gridstruct%del6_u(i,j)*(d2(i,j-1)-d2(i,j))
1086 #endif
1087  enddo
1088  enddo
1089 
1090  if ( nord>0 ) then
1091 
1092 !----------
1093 ! high-order
1094 !----------
1095 
1096  do n=1, nord
1097 
1098  nt = nord-n
1099 
1100  do j=js-nt-1,je+nt+1
1101  do i=is-nt-1,ie+nt+1
1102  d2(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j)
1103  enddo
1104  enddo
1105 
1106  call copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, &
1107  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1108  do j=js-nt,je+nt
1109  do i=is-nt,ie+nt+1
1110 #ifdef USE_SG
1111  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)
1112 #else
1113  fx2(i,j) = gridstruct%del6_v(i,j)*(d2(i,j)-d2(i-1,j))
1114 #endif
1115  enddo
1116  enddo
1117 
1118  call copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, &
1119  gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
1120  do j=js-nt,je+nt+1
1121  do i=is-nt,ie+nt
1122 #ifdef USE_SG
1123  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)
1124 #else
1125  fy2(i,j) = gridstruct%del6_u(i,j)*(d2(i,j)-d2(i,j-1))
1126 #endif
1127  enddo
1128  enddo
1129  enddo
1130 
1131  endif
1132 
1133 !---------------------------------------------
1134 ! Add the diffusive fluxes to the flux arrays:
1135 !---------------------------------------------
1136 
1137  if ( present(mass) ) then
1138 ! Apply mass weighting to diffusive fluxes:
1139  damp2 = 0.5*damp
1140  do j=js,je
1141  do i=is,ie+1
1142  fx(i,j) = fx(i,j) + damp2*(mass(i-1,j)+mass(i,j))*fx2(i,j)
1143  enddo
1144  enddo
1145  do j=js,je+1
1146  do i=is,ie
1147  fy(i,j) = fy(i,j) + damp2*(mass(i,j-1)+mass(i,j))*fy2(i,j)
1148  enddo
1149  enddo
1150  else
1151  do j=js,je
1152  do i=is,ie+1
1153  fx(i,j) = fx(i,j) + fx2(i,j)
1154  enddo
1155  enddo
1156  do j=js,je+1
1157  do i=is,ie
1158  fy(i,j) = fy(i,j) + fy2(i,j)
1159  enddo
1160  enddo
1161  endif
1162 
1163  end subroutine deln_flux
1164 
1165 
1166 end module tp_core_nlm_mod
tp_core_nlm_mod::b2
real, parameter b2
Definition: tp_core_nlm.F90:50
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
tp_core_nlm_mod::b4
real, parameter b4
Definition: tp_core_nlm.F90:52
fv_arrays_nlm_mod::fv_grid_type
Definition: fv_arrays_nlm.F90:115
tp_core_nlm_mod::r3
real, parameter r3
Definition: tp_core_nlm.F90:35
tp_core_nlm_mod::mp_ghost_ew
subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
Definition: tp_core_nlm.F90:905
tp_core_nlm_mod::pert_ppm
subroutine, public pert_ppm(im, a0, al, ar, iv)
Definition: tp_core_nlm.F90:954
tp_core_nlm_mod::ppm_limiter
real, parameter ppm_limiter
Definition: tp_core_nlm.F90:37
tp_core_nlm_mod::near_zero
real, parameter near_zero
Definition: tp_core_nlm.F90:36
tp_core_nlm_mod::p2
real, parameter p2
Definition: tp_core_nlm.F90:68
tp_core_nlm_mod::b3
real, parameter b3
Definition: tp_core_nlm.F90:51
tp_core_nlm_mod::c3
real, parameter c3
Definition: tp_core_nlm.F90:63
tp_core_nlm_mod::t11
real, parameter t11
Definition: tp_core_nlm.F90:55
tp_core_nlm_mod::t13
real, parameter t13
Definition: tp_core_nlm.F90:55
tp_core_nlm_mod::c1
real, parameter c1
Definition: tp_core_nlm.F90:61
tp_core_nlm_mod::b1
real, parameter b1
Definition: tp_core_nlm.F90:49
fv_grid_utils_nlm_mod
Definition: fv_grid_utils_nlm.F90:20
tp_core_nlm_mod::deln_flux
subroutine deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy, gridstruct, bd, mass)
Definition: tp_core_nlm.F90:1016
fv_grid_utils_nlm_mod::big_number
real, parameter, public big_number
Definition: fv_grid_utils_nlm.F90:49
tp_core_nlm_mod::ppm_fac
real, parameter ppm_fac
Definition: tp_core_nlm.F90:34
fv_mp_nlm_mod::ng
integer, parameter, public ng
Definition: fv_mp_nlm_mod.F90:2465
tp_core_nlm_mod::s11
real, parameter s11
Definition: tp_core_nlm.F90:56
fv_arrays_nlm_mod
Definition: fv_arrays_nlm.F90:20
tp_core_nlm_mod::b5
real, parameter b5
Definition: tp_core_nlm.F90:53
tp_core_nlm_mod::c2
real, parameter c2
Definition: tp_core_nlm.F90:62
tp_core_nlm_mod::s14
real, parameter s14
Definition: tp_core_nlm.F90:56
tp_core_nlm_mod::s15
real, parameter s15
Definition: tp_core_nlm.F90:56
tp_core_nlm_mod::p1
real, parameter p1
Definition: tp_core_nlm.F90:67
max
#define max(a, b)
Definition: mosaic_util.h:33
tp_core_nlm_mod::xppm
subroutine xppm(flux, q, c, iord, is, ie, isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
Definition: tp_core_nlm.F90:292
tp_core_nlm_mod::t12
real, parameter t12
Definition: tp_core_nlm.F90:55
tp_core_nlm_mod
Definition: tp_core_nlm.F90:20
tp_core_nlm_mod::copy_corners
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core_nlm.F90:216
min
#define min(a, b)
Definition: mosaic_util.h:32
tp_core_nlm_mod::fv_tp_2d
subroutine, public fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
Definition: tp_core_nlm.F90:80
fv_arrays_nlm_mod::fv_grid_bounds_type
Definition: fv_arrays_nlm.F90:601
tp_core_nlm_mod::yppm
subroutine yppm(flux, q, c, jord, ifirst, ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
Definition: tp_core_nlm.F90:580
unstructured_grid_mod::grid_type
Derived type containing the data.
Definition: unstructured_grid_mod.F90:23