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