FV3 Bundle
tp_core_tlm.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_tlm_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, r_grid
28 
29  implicit none
30 
31  private
32  public fv_tp_2d, pert_ppm, copy_corners
33  public fv_tp_2d_tlm, copy_corners_tlm
34 
35  real, parameter:: ppm_fac = 1.5 ! nonlinear scheme limiter: between 1 and 2
36  real, parameter:: r3 = 1./3.
37  real, parameter:: near_zero = 1.e-25
38  real, parameter:: ppm_limiter = 2.0
39 
40 #ifdef WAVE_FORM
41 ! Suresh & Huynh scheme 2.2 (purtabation form)
42 ! The wave-form is more diffusive than scheme 2.1
43  real, parameter:: b1 = 0.0375
44  real, parameter:: b2 = -7./30.
45  real, parameter:: b3 = -23./120.
46  real, parameter:: b4 = 13./30.
47  real, parameter:: b5 = -11./240.
48 #else
49 ! scheme 2.1: perturbation form
50  real, parameter:: b1 = 1./30.
51  real, parameter:: b2 = -13./60.
52  real, parameter:: b3 = -13./60.
53  real, parameter:: b4 = 0.45
54  real, parameter:: b5 = -0.05
55 #endif
56  real, parameter:: t11 = 27./28., t12 = -13./28., t13=3./7.
57  real, parameter:: s11 = 11./14., s14 = 4./7., s15=3./14.
58 !----------------------------------------------------
59 ! volume-conserving cubic with 2nd drv=0 at end point:
60 !----------------------------------------------------
61 ! Non-monotonic
62  real, parameter:: c1 = -2./14.
63  real, parameter:: c2 = 11./14.
64  real, parameter:: c3 = 5./14.
65 !----------------------
66 ! PPM volume mean form:
67 !----------------------
68  real, parameter:: p1 = 7./12. ! 0.58333333
69  real, parameter:: p2 = -1./12.
70 ! q(i+0.5) = p1*(q(i-1)+q(i)) + p2*(q(i-2)+q(i+1))
71 ! integer:: is, ie, js, je, isd, ied, jsd, jed
72 
73 !
74 !EOP
75 !-----------------------------------------------------------------------
76 
77 CONTAINS
78 ! q(i+0.5) = p1*(q(i-1)+q(i)) + p2*(q(i-2)+q(i+1))
79 ! integer:: is, ie, js, je, isd, ied, jsd, jed
80 !
81 !EOP
82 !-----------------------------------------------------------------------
83  SUBROUTINE fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, &
84 & gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
85  IMPLICIT NONE
86  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
87  INTEGER, INTENT(IN) :: npx, npy
88  INTEGER, INTENT(IN) :: hord
89 !
90  REAL, INTENT(IN) :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed)
91 !
92  REAL, INTENT(IN) :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed)
93 !
94  REAL, INTENT(IN) :: cry(bd%isd:bd%ied, bd%js:bd%je+1)
95 !
96  REAL, INTENT(IN) :: yfx(bd%isd:bd%ied, bd%js:bd%je+1)
97  REAL, INTENT(IN) :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
98  REAL, INTENT(IN) :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
99 ! transported scalar
100  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
101 ! Flux in x ( E )
102  REAL, INTENT(OUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je)
103 ! Flux in y ( N )
104  REAL, INTENT(OUT) :: fy(bd%is:bd%ie, bd%js:bd%je+1)
105  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
106 ! optional Arguments:
107 ! Mass Flux X-Dir
108  REAL, OPTIONAL, INTENT(IN) :: mfx(bd%is:bd%ie+1, bd%js:bd%je)
109 ! Mass Flux Y-Dir
110  REAL, OPTIONAL, INTENT(IN) :: mfy(bd%is:bd%ie, bd%js:bd%je+1)
111  REAL, OPTIONAL, INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
112  REAL, OPTIONAL, INTENT(IN) :: damp_c
113  INTEGER, OPTIONAL, INTENT(IN) :: nord
114 ! Local:
115  INTEGER :: ord_ou, ord_in
116  REAL :: q_i(bd%isd:bd%ied, bd%js:bd%je)
117  REAL :: q_j(bd%is:bd%ie, bd%jsd:bd%jed)
118  REAL :: fx2(bd%is:bd%ie+1, bd%jsd:bd%jed)
119  REAL :: fy2(bd%isd:bd%ied, bd%js:bd%je+1)
120  REAL :: fyy(bd%isd:bd%ied, bd%js:bd%je+1)
121  REAL :: fx1(bd%is:bd%ie+1)
122  REAL :: damp
123  INTEGER :: i, j
124  INTEGER :: is, ie, js, je, isd, ied, jsd, jed
125  INTRINSIC PRESENT
126  REAL*8 :: pwx1
127  INTEGER :: pwy1
128  is = bd%is
129  ie = bd%ie
130  js = bd%js
131  je = bd%je
132  isd = bd%isd
133  ied = bd%ied
134  jsd = bd%jsd
135  jed = bd%jed
136  IF (hord .EQ. 10) THEN
137  ord_in = 8
138  ELSE
139  ord_in = hord
140  END IF
141  ord_ou = hord
142  IF (.NOT.gridstruct%nested) CALL copy_corners(q, npx, npy, 2, &
143 & gridstruct%nested, bd, &
144 & gridstruct%sw_corner, &
145 & gridstruct%se_corner, &
146 & gridstruct%nw_corner, &
147 & gridstruct%ne_corner)
148  CALL yppm(fy2, q, cry, ord_in, isd, ied, isd, ied, js, je, jsd, jed&
149 & , npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%&
150 & grid_type)
151  DO j=js,je+1
152  DO i=isd,ied
153  fyy(i, j) = yfx(i, j)*fy2(i, j)
154  END DO
155  END DO
156  DO j=js,je
157  DO i=isd,ied
158  q_i(i, j) = (q(i, j)*gridstruct%area(i, j)+fyy(i, j)-fyy(i, j+1)&
159 & )/ra_y(i, j)
160  END DO
161  END DO
162  CALL xppm(fx, q_i, crx(is:ie+1, js:je), ord_ou, is, ie, isd, ied, js&
163 & , je, jsd, jed, npx, npy, gridstruct%dxa, gridstruct%nested, &
164 & gridstruct%grid_type)
165  IF (.NOT.gridstruct%nested) CALL copy_corners(q, npx, npy, 1, &
166 & gridstruct%nested, bd, &
167 & gridstruct%sw_corner, &
168 & gridstruct%se_corner, &
169 & gridstruct%nw_corner, &
170 & gridstruct%ne_corner)
171  CALL xppm(fx2, q, crx, ord_in, is, ie, isd, ied, jsd, jed, jsd, jed&
172 & , npx, npy, gridstruct%dxa, gridstruct%nested, gridstruct%&
173 & grid_type)
174  DO j=jsd,jed
175  DO i=is,ie+1
176  fx1(i) = xfx(i, j)*fx2(i, j)
177  END DO
178  DO i=is,ie
179  q_j(i, j) = (q(i, j)*gridstruct%area(i, j)+fx1(i)-fx1(i+1))/ra_x&
180 & (i, j)
181  END DO
182  END DO
183  CALL yppm(fy, q_j, cry, ord_ou, is, ie, isd, ied, js, je, jsd, jed, &
184 & npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%&
185 & grid_type)
186 !----------------
187 ! Flux averaging:
188 !----------------
189  IF (PRESENT(mfx) .AND. PRESENT(mfy)) THEN
190 !---------------------------------
191 ! For transport of pt and tracers
192 !---------------------------------
193  DO j=js,je
194  DO i=is,ie+1
195  fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*mfx(i, j)
196  END DO
197  END DO
198  DO j=js,je+1
199  DO i=is,ie
200  fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*mfy(i, j)
201  END DO
202  END DO
203  IF (PRESENT(nord) .AND. PRESENT(damp_c) .AND. PRESENT(mass)) THEN
204  IF (damp_c .GT. 1.e-4) THEN
205  pwx1 = damp_c*gridstruct%da_min
206  pwy1 = nord + 1
207  damp = pwx1**pwy1
208  CALL deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy&
209 & , gridstruct, bd, mass)
210  END IF
211  END IF
212  ELSE
213 !---------------------------------
214 ! For transport of delp, vorticity
215 !---------------------------------
216  DO j=js,je
217  DO i=is,ie+1
218  fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*xfx(i, j)
219  END DO
220  END DO
221  DO j=js,je+1
222  DO i=is,ie
223  fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*yfx(i, j)
224  END DO
225  END DO
226  IF (PRESENT(nord) .AND. PRESENT(damp_c)) THEN
227  IF (damp_c .GT. 1.e-4) THEN
228  pwx1 = damp_c*gridstruct%da_min
229  pwy1 = nord + 1
230  damp = pwx1**pwy1
231  CALL deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy&
232 & , gridstruct, bd)
233  END IF
234  END IF
235  END IF
236  END SUBROUTINE fv_tp_2d
237  SUBROUTINE xppm(flux, q, c, iord, is, ie, isd, ied, jfirst, jlast, jsd&
238 & , jed, npx, npy, dxa, nested, grid_type)
239  IMPLICIT NONE
240  INTEGER, INTENT(IN) :: is, ie, isd, ied, jsd, jed
241 ! compute domain
242  INTEGER, INTENT(IN) :: jfirst, jlast
243  INTEGER, INTENT(IN) :: iord
244  INTEGER, INTENT(IN) :: npx, npy
245  REAL, INTENT(IN) :: q(isd:ied, jfirst:jlast)
246 ! Courant N (like FLUX)
247  REAL, INTENT(IN) :: c(is:ie+1, jfirst:jlast)
248  REAL, INTENT(IN) :: dxa(isd:ied, jsd:jed)
249  LOGICAL, INTENT(IN) :: nested
250  INTEGER, INTENT(IN) :: grid_type
251 ! !OUTPUT PARAMETERS:
252 ! Flux
253  REAL, INTENT(OUT) :: flux(is:ie+1, jfirst:jlast)
254 ! Local
255  REAL, DIMENSION(is-1:ie+1) :: bl, br, b0
256  REAL :: q1(isd:ied)
257  REAL, DIMENSION(is:ie+1) :: fx0, fx1
258  LOGICAL, DIMENSION(is-1:ie+1) :: smt5, smt6
259  REAL :: al(is-1:ie+2)
260  REAL :: dm(is-2:ie+2)
261  REAL :: dq(is-3:ie+2)
262  INTEGER :: i, j, ie3, is1, ie1
263  REAL :: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
264  INTRINSIC max
265  INTRINSIC min
266  INTRINSIC abs
267  INTRINSIC sign
268  REAL :: min1
269  REAL :: min2
270  REAL :: abs0
271  REAL :: abs1
272  REAL :: min3
273  REAL :: min4
274  REAL :: min5
275  REAL :: min6
276  REAL :: min7
277  REAL :: abs2
278  REAL :: abs3
279  REAL :: abs4
280  REAL :: max1
281  REAL :: min8
282  REAL :: abs5
283  REAL :: abs6
284  REAL :: abs7
285  REAL :: x10
286  REAL :: x9
287  REAL :: x8
288  REAL :: x7
289  REAL :: x6
290  REAL :: x5
291  REAL :: x4
292  REAL :: x3
293  REAL :: x2
294  REAL :: y15
295  REAL :: y14
296  REAL :: y13
297  REAL :: y12
298  REAL :: y11
299  REAL :: y10
300  REAL :: z5
301  REAL :: z4
302  REAL :: z3
303  REAL :: z2
304  REAL :: z1
305  REAL :: y9
306  REAL :: y8
307  REAL :: y7
308  REAL :: y6
309  REAL :: y5
310  REAL :: y4
311  REAL :: y3
312  REAL :: y2
313  REAL :: y1
314  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
315  IF (3 .LT. is - 1) THEN
316  is1 = is - 1
317  ELSE
318  is1 = 3
319  END IF
320  IF (npx - 2 .GT. ie + 2) THEN
321  ie3 = ie + 2
322  ELSE
323  ie3 = npx - 2
324  END IF
325  IF (npx - 3 .GT. ie + 1) THEN
326  ie1 = ie + 1
327  ELSE
328  ie1 = npx - 3
329  END IF
330  ELSE
331  is1 = is - 1
332  ie3 = ie + 2
333  ie1 = ie + 1
334  END IF
335  DO j=jfirst,jlast
336  DO i=isd,ied
337  q1(i) = q(i, j)
338  END DO
339  IF (iord .LT. 8 .OR. iord .EQ. 333) THEN
340 ! ord = 2: perfectly linear ppm scheme
341 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
342  DO i=is1,ie3
343  al(i) = p1*(q1(i-1)+q1(i)) + p2*(q1(i-2)+q1(i+1))
344  END DO
345  IF (iord .EQ. 7) THEN
346  DO i=is1,ie3
347  IF (al(i) .LT. 0.) al(i) = 0.5*(q1(i-1)+q1(i))
348  END DO
349  END IF
350  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
351  IF (is .EQ. 1) THEN
352  al(0) = c1*q1(-2) + c2*q1(-1) + c3*q1(0)
353  al(1) = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-&
354 & 1))/(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)&
355 & -dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
356  al(2) = c3*q1(1) + c2*q1(2) + c1*q1(3)
357  IF (iord .EQ. 7) THEN
358  IF (0. .LT. al(0)) THEN
359  al(0) = al(0)
360  ELSE
361  al(0) = 0.
362  END IF
363  IF (0. .LT. al(1)) THEN
364  al(1) = al(1)
365  ELSE
366  al(1) = 0.
367  END IF
368  IF (0. .LT. al(2)) THEN
369  al(2) = al(2)
370  ELSE
371  al(2) = 0.
372  END IF
373  END IF
374  END IF
375  IF (ie + 1 .EQ. npx) THEN
376  al(npx-1) = c1*q1(npx-3) + c2*q1(npx-2) + c3*q1(npx-1)
377  al(npx) = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-&
378 & dxa(npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((&
379 & 2.*dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1&
380 & ))/(dxa(npx, j)+dxa(npx+1, j)))
381  al(npx+1) = c3*q1(npx) + c2*q1(npx+1) + c1*q1(npx+2)
382  IF (iord .EQ. 7) THEN
383  IF (0. .LT. al(npx-1)) THEN
384  al(npx-1) = al(npx-1)
385  ELSE
386  al(npx-1) = 0.
387  END IF
388  IF (0. .LT. al(npx)) THEN
389  al(npx) = al(npx)
390  ELSE
391  al(npx) = 0.
392  END IF
393  IF (0. .LT. al(npx+1)) THEN
394  al(npx+1) = al(npx+1)
395  ELSE
396  al(npx+1) = 0.
397  END IF
398  END IF
399  END IF
400  END IF
401  IF (iord .EQ. 1) THEN
402  DO i=is,ie+1
403  IF (c(i, j) .GT. 0.) THEN
404  flux(i, j) = q1(i-1)
405  ELSE
406  flux(i, j) = q1(i)
407  END IF
408  END DO
409  ELSE IF (iord .EQ. 2) THEN
410 ! perfectly linear scheme
411 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7
412 !DEC$ VECTOR ALWAYS
413  DO i=is,ie+1
414  xt = c(i, j)
415  IF (xt .GT. 0.) THEN
416  qtmp = q1(i-1)
417  flux(i, j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-&
418 & (qtmp+qtmp)))
419  ELSE
420  qtmp = q1(i)
421  flux(i, j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-&
422 & (qtmp+qtmp)))
423  END IF
424  END DO
425  ELSE IF (iord .EQ. 333) THEN
426 ! x0 = sign(dim(xt, 0.), 1.)
427 ! x1 = sign(dim(0., xt), 1.)
428 ! flux(i,j) = x0*(q1(i-1)+(1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-(qtmp+qtmp)))) &
429 ! + x1*(q1(i) +(1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-(qtmp+qtmp))))
430 ! Perfectly linear scheme, more diffusive than ord=2 (HoldawayKent-2015-TellusA)
431 !DEC$ VECTOR ALWAYS
432  DO i=is,ie+1
433  xt = c(i, j)
434  IF (xt .GT. 0.) THEN
435  flux(i, j) = (2.0*q1(i)+5.0*q1(i-1)-q1(i-2))/6.0 - 0.5*xt*&
436 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i)-2.0*q1(i-1)+q1(i-2))
437  ELSE
438  flux(i, j) = (2.0*q1(i-1)+5.0*q1(i)-q1(i+1))/6.0 - 0.5*xt*&
439 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i+1)-2.0*q1(i)+q1(i-1))
440  END IF
441  END DO
442  ELSE IF (iord .EQ. 3) THEN
443  DO i=is-1,ie+1
444  bl(i) = al(i) - q1(i)
445  br(i) = al(i+1) - q1(i)
446  b0(i) = bl(i) + br(i)
447  IF (b0(i) .GE. 0.) THEN
448  x0 = b0(i)
449  ELSE
450  x0 = -b0(i)
451  END IF
452  IF (bl(i) - br(i) .GE. 0.) THEN
453  xt = bl(i) - br(i)
454  ELSE
455  xt = -(bl(i)-br(i))
456  END IF
457  smt5(i) = x0 .LT. xt
458  smt6(i) = 3.*x0 .LT. xt
459  END DO
460  DO i=is,ie+1
461  fx1(i) = 0.
462  END DO
463  DO i=is,ie+1
464  xt = c(i, j)
465  IF (xt .GT. 0.) THEN
466  fx0(i) = q1(i-1)
467  IF (smt6(i-1) .OR. smt5(i)) THEN
468  fx1(i) = br(i-1) - xt*b0(i-1)
469  ELSE IF (smt5(i-1)) THEN
470  IF (bl(i-1) .GE. 0.) THEN
471  x2 = bl(i-1)
472  ELSE
473  x2 = -bl(i-1)
474  END IF
475  IF (br(i-1) .GE. 0.) THEN
476  y1 = br(i-1)
477  ELSE
478  y1 = -br(i-1)
479  END IF
480  IF (x2 .GT. y1) THEN
481  min1 = y1
482  ELSE
483  min1 = x2
484  END IF
485 ! 2nd order, piece-wise linear
486  fx1(i) = sign(min1, br(i-1))
487  END IF
488  ELSE
489  fx0(i) = q1(i)
490  IF (smt6(i) .OR. smt5(i-1)) THEN
491  fx1(i) = bl(i) + xt*b0(i)
492  ELSE IF (smt5(i)) THEN
493  IF (bl(i) .GE. 0.) THEN
494  x3 = bl(i)
495  ELSE
496  x3 = -bl(i)
497  END IF
498  IF (br(i) .GE. 0.) THEN
499  y2 = br(i)
500  ELSE
501  y2 = -br(i)
502  END IF
503  IF (x3 .GT. y2) THEN
504  min2 = y2
505  ELSE
506  min2 = x3
507  END IF
508  fx1(i) = sign(min2, bl(i))
509  END IF
510  END IF
511  IF (xt .GE. 0.) THEN
512  abs0 = xt
513  ELSE
514  abs0 = -xt
515  END IF
516  flux(i, j) = fx0(i) + (1.-abs0)*fx1(i)
517  END DO
518  ELSE IF (iord .EQ. 4) THEN
519  DO i=is-1,ie+1
520  bl(i) = al(i) - q1(i)
521  br(i) = al(i+1) - q1(i)
522  b0(i) = bl(i) + br(i)
523  IF (b0(i) .GE. 0.) THEN
524  x0 = b0(i)
525  ELSE
526  x0 = -b0(i)
527  END IF
528  IF (bl(i) - br(i) .GE. 0.) THEN
529  xt = bl(i) - br(i)
530  ELSE
531  xt = -(bl(i)-br(i))
532  END IF
533  smt5(i) = x0 .LT. xt
534  smt6(i) = 3.*x0 .LT. xt
535  END DO
536  DO i=is,ie+1
537  fx1(i) = 0.
538  END DO
539 !DEC$ VECTOR ALWAYS
540  DO i=is,ie+1
541  IF (c(i, j) .GT. 0.) THEN
542  fx0(i) = q1(i-1)
543  IF (smt6(i-1) .OR. smt5(i)) fx1(i) = (1.-c(i, j))*(br(i-1)&
544 & -c(i, j)*b0(i-1))
545  ELSE
546  fx0(i) = q1(i)
547  IF (smt6(i) .OR. smt5(i-1)) fx1(i) = (1.+c(i, j))*(bl(i)+c&
548 & (i, j)*b0(i))
549  END IF
550  flux(i, j) = fx0(i) + fx1(i)
551  END DO
552  ELSE
553 ! iord = 5 & 6
554  IF (iord .EQ. 5) THEN
555  DO i=is-1,ie+1
556  bl(i) = al(i) - q1(i)
557  br(i) = al(i+1) - q1(i)
558  b0(i) = bl(i) + br(i)
559  smt5(i) = bl(i)*br(i) .LT. 0.
560  END DO
561  ELSE
562  DO i=is-1,ie+1
563  bl(i) = al(i) - q1(i)
564  br(i) = al(i+1) - q1(i)
565  b0(i) = bl(i) + br(i)
566  IF (3.*b0(i) .GE. 0.) THEN
567  abs1 = 3.*b0(i)
568  ELSE
569  abs1 = -(3.*b0(i))
570  END IF
571  IF (bl(i) - br(i) .GE. 0.) THEN
572  abs4 = bl(i) - br(i)
573  ELSE
574  abs4 = -(bl(i)-br(i))
575  END IF
576  smt5(i) = abs1 .LT. abs4
577  END DO
578  END IF
579 !DEC$ VECTOR ALWAYS
580  DO i=is,ie+1
581  IF (c(i, j) .GT. 0.) THEN
582  fx1(i) = (1.-c(i, j))*(br(i-1)-c(i, j)*b0(i-1))
583  flux(i, j) = q1(i-1)
584  ELSE
585  fx1(i) = (1.+c(i, j))*(bl(i)+c(i, j)*b0(i))
586  flux(i, j) = q1(i)
587  END IF
588  IF (smt5(i-1) .OR. smt5(i)) flux(i, j) = flux(i, j) + fx1(i)
589  END DO
590  END IF
591  ELSE
592 ! Monotonic constraints:
593 ! ord = 8: PPM with Lin's PPM fast monotone constraint
594 ! ord = 10: PPM with Lin's modification of Huynh 2nd constraint
595 ! ord = 13: 10 plus positive definite constraint
596  DO i=is-2,ie+2
597  xt = 0.25*(q1(i+1)-q1(i-1))
598  IF (xt .GE. 0.) THEN
599  x4 = xt
600  ELSE
601  x4 = -xt
602  END IF
603  IF (q1(i-1) .LT. q1(i)) THEN
604  IF (q1(i) .LT. q1(i+1)) THEN
605  max1 = q1(i+1)
606  ELSE
607  max1 = q1(i)
608  END IF
609  ELSE IF (q1(i-1) .LT. q1(i+1)) THEN
610  max1 = q1(i+1)
611  ELSE
612  max1 = q1(i-1)
613  END IF
614  y3 = max1 - q1(i)
615  IF (q1(i-1) .GT. q1(i)) THEN
616  IF (q1(i) .GT. q1(i+1)) THEN
617  min8 = q1(i+1)
618  ELSE
619  min8 = q1(i)
620  END IF
621  ELSE IF (q1(i-1) .GT. q1(i+1)) THEN
622  min8 = q1(i+1)
623  ELSE
624  min8 = q1(i-1)
625  END IF
626  z1 = q1(i) - min8
627  IF (x4 .GT. y3) THEN
628  IF (y3 .GT. z1) THEN
629  min3 = z1
630  ELSE
631  min3 = y3
632  END IF
633  ELSE IF (x4 .GT. z1) THEN
634  min3 = z1
635  ELSE
636  min3 = x4
637  END IF
638  dm(i) = sign(min3, xt)
639  END DO
640  DO i=is1,ie1+1
641  al(i) = 0.5*(q1(i-1)+q1(i)) + r3*(dm(i-1)-dm(i))
642  END DO
643  IF (iord .EQ. 8) THEN
644  DO i=is1,ie1
645  xt = 2.*dm(i)
646  IF (xt .GE. 0.) THEN
647  x5 = xt
648  ELSE
649  x5 = -xt
650  END IF
651  IF (al(i) - q1(i) .GE. 0.) THEN
652  y4 = al(i) - q1(i)
653  ELSE
654  y4 = -(al(i)-q1(i))
655  END IF
656  IF (x5 .GT. y4) THEN
657  min4 = y4
658  ELSE
659  min4 = x5
660  END IF
661  bl(i) = -sign(min4, xt)
662  IF (xt .GE. 0.) THEN
663  x6 = xt
664  ELSE
665  x6 = -xt
666  END IF
667  IF (al(i+1) - q1(i) .GE. 0.) THEN
668  y5 = al(i+1) - q1(i)
669  ELSE
670  y5 = -(al(i+1)-q1(i))
671  END IF
672  IF (x6 .GT. y5) THEN
673  min5 = y5
674  ELSE
675  min5 = x6
676  END IF
677  br(i) = sign(min5, xt)
678  END DO
679  ELSE IF (iord .EQ. 11) THEN
680 ! This is emulation of 2nd van Leer scheme using PPM codes
681  DO i=is1,ie1
682  xt = ppm_fac*dm(i)
683  IF (xt .GE. 0.) THEN
684  x7 = xt
685  ELSE
686  x7 = -xt
687  END IF
688  IF (al(i) - q1(i) .GE. 0.) THEN
689  y6 = al(i) - q1(i)
690  ELSE
691  y6 = -(al(i)-q1(i))
692  END IF
693  IF (x7 .GT. y6) THEN
694  min6 = y6
695  ELSE
696  min6 = x7
697  END IF
698  bl(i) = -sign(min6, xt)
699  IF (xt .GE. 0.) THEN
700  x8 = xt
701  ELSE
702  x8 = -xt
703  END IF
704  IF (al(i+1) - q1(i) .GE. 0.) THEN
705  y7 = al(i+1) - q1(i)
706  ELSE
707  y7 = -(al(i+1)-q1(i))
708  END IF
709  IF (x8 .GT. y7) THEN
710  min7 = y7
711  ELSE
712  min7 = x8
713  END IF
714  br(i) = sign(min7, xt)
715  END DO
716  ELSE
717  DO i=is1-2,ie1+1
718  dq(i) = 2.*(q1(i+1)-q1(i))
719  END DO
720  DO i=is1,ie1
721  bl(i) = al(i) - q1(i)
722  br(i) = al(i+1) - q1(i)
723  IF (dm(i-1) .GE. 0.) THEN
724  abs2 = dm(i-1)
725  ELSE
726  abs2 = -dm(i-1)
727  END IF
728  IF (dm(i) .GE. 0.) THEN
729  abs5 = dm(i)
730  ELSE
731  abs5 = -dm(i)
732  END IF
733  IF (dm(i+1) .GE. 0.) THEN
734  abs7 = dm(i+1)
735  ELSE
736  abs7 = -dm(i+1)
737  END IF
738  IF (abs2 + abs5 + abs7 .LT. near_zero) THEN
739  bl(i) = 0.
740  br(i) = 0.
741  ELSE
742  IF (3.*(bl(i)+br(i)) .GE. 0.) THEN
743  abs3 = 3.*(bl(i)+br(i))
744  ELSE
745  abs3 = -(3.*(bl(i)+br(i)))
746  END IF
747  IF (bl(i) - br(i) .GE. 0.) THEN
748  abs6 = bl(i) - br(i)
749  ELSE
750  abs6 = -(bl(i)-br(i))
751  END IF
752  IF (abs3 .GT. abs6) THEN
753  pmp_2 = dq(i-1)
754  lac_2 = pmp_2 - 0.75*dq(i-2)
755  IF (0. .LT. pmp_2) THEN
756  IF (pmp_2 .LT. lac_2) THEN
757  x9 = lac_2
758  ELSE
759  x9 = pmp_2
760  END IF
761  ELSE IF (0. .LT. lac_2) THEN
762  x9 = lac_2
763  ELSE
764  x9 = 0.
765  END IF
766  IF (0. .GT. pmp_2) THEN
767  IF (pmp_2 .GT. lac_2) THEN
768  y14 = lac_2
769  ELSE
770  y14 = pmp_2
771  END IF
772  ELSE IF (0. .GT. lac_2) THEN
773  y14 = lac_2
774  ELSE
775  y14 = 0.
776  END IF
777  IF (br(i) .LT. y14) THEN
778  y8 = y14
779  ELSE
780  y8 = br(i)
781  END IF
782  IF (x9 .GT. y8) THEN
783  br(i) = y8
784  ELSE
785  br(i) = x9
786  END IF
787  pmp_1 = -dq(i)
788  lac_1 = pmp_1 + 0.75*dq(i+1)
789  IF (0. .LT. pmp_1) THEN
790  IF (pmp_1 .LT. lac_1) THEN
791  x10 = lac_1
792  ELSE
793  x10 = pmp_1
794  END IF
795  ELSE IF (0. .LT. lac_1) THEN
796  x10 = lac_1
797  ELSE
798  x10 = 0.
799  END IF
800  IF (0. .GT. pmp_1) THEN
801  IF (pmp_1 .GT. lac_1) THEN
802  y15 = lac_1
803  ELSE
804  y15 = pmp_1
805  END IF
806  ELSE IF (0. .GT. lac_1) THEN
807  y15 = lac_1
808  ELSE
809  y15 = 0.
810  END IF
811  IF (bl(i) .LT. y15) THEN
812  y9 = y15
813  ELSE
814  y9 = bl(i)
815  END IF
816  IF (x10 .GT. y9) THEN
817  bl(i) = y9
818  ELSE
819  bl(i) = x10
820  END IF
821  END IF
822  END IF
823  END DO
824  END IF
825 ! Positive definite constraint:
826  IF (iord .EQ. 9 .OR. iord .EQ. 13) CALL pert_ppm(ie1 - is1 + 1, &
827 & q1(is1:ie1), bl(is1:&
828 & ie1), br(is1:ie1), 0)
829  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
830  IF (is .EQ. 1) THEN
831  bl(0) = s14*dm(-1) + s11*(q1(-1)-q1(0))
832  xt = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-1))&
833 & /(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)-&
834 & dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
835  IF (q1(1) .GT. q1(2)) THEN
836  z2 = q1(2)
837  ELSE
838  z2 = q1(1)
839  END IF
840  IF (q1(-1) .GT. q1(0)) THEN
841  IF (q1(0) .GT. z2) THEN
842  y10 = z2
843  ELSE
844  y10 = q1(0)
845  END IF
846  ELSE IF (q1(-1) .GT. z2) THEN
847  y10 = z2
848  ELSE
849  y10 = q1(-1)
850  END IF
851  IF (xt .LT. y10) THEN
852  xt = y10
853  ELSE
854  xt = xt
855  END IF
856  IF (q1(1) .LT. q1(2)) THEN
857  z3 = q1(2)
858  ELSE
859  z3 = q1(1)
860  END IF
861  IF (q1(-1) .LT. q1(0)) THEN
862  IF (q1(0) .LT. z3) THEN
863  y11 = z3
864  ELSE
865  y11 = q1(0)
866  END IF
867  ELSE IF (q1(-1) .LT. z3) THEN
868  y11 = z3
869  ELSE
870  y11 = q1(-1)
871  END IF
872  IF (xt .GT. y11) THEN
873  xt = y11
874  ELSE
875  xt = xt
876  END IF
877 ! endif
878  br(0) = xt - q1(0)
879  bl(1) = xt - q1(1)
880  xt = s15*q1(1) + s11*q1(2) - s14*dm(2)
881  br(1) = xt - q1(1)
882  bl(2) = xt - q1(2)
883  br(2) = al(3) - q1(2)
884  CALL pert_ppm(3, q1(0:2), bl(0:2), br(0:2), 1)
885  END IF
886  IF (ie + 1 .EQ. npx) THEN
887  bl(npx-2) = al(npx-2) - q1(npx-2)
888  xt = s15*q1(npx-1) + s11*q1(npx-2) + s14*dm(npx-2)
889  br(npx-2) = xt - q1(npx-2)
890  bl(npx-1) = xt - q1(npx-1)
891  xt = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-dxa(&
892 & npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((2.*&
893 & dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1))/&
894 & (dxa(npx, j)+dxa(npx+1, j)))
895  IF (q1(npx) .GT. q1(npx+1)) THEN
896  z4 = q1(npx+1)
897  ELSE
898  z4 = q1(npx)
899  END IF
900  IF (q1(npx-2) .GT. q1(npx-1)) THEN
901  IF (q1(npx-1) .GT. z4) THEN
902  y12 = z4
903  ELSE
904  y12 = q1(npx-1)
905  END IF
906  ELSE IF (q1(npx-2) .GT. z4) THEN
907  y12 = z4
908  ELSE
909  y12 = q1(npx-2)
910  END IF
911  IF (xt .LT. y12) THEN
912  xt = y12
913  ELSE
914  xt = xt
915  END IF
916  IF (q1(npx) .LT. q1(npx+1)) THEN
917  z5 = q1(npx+1)
918  ELSE
919  z5 = q1(npx)
920  END IF
921  IF (q1(npx-2) .LT. q1(npx-1)) THEN
922  IF (q1(npx-1) .LT. z5) THEN
923  y13 = z5
924  ELSE
925  y13 = q1(npx-1)
926  END IF
927  ELSE IF (q1(npx-2) .LT. z5) THEN
928  y13 = z5
929  ELSE
930  y13 = q1(npx-2)
931  END IF
932  IF (xt .GT. y13) THEN
933  xt = y13
934  ELSE
935  xt = xt
936  END IF
937 ! endif
938  br(npx-1) = xt - q1(npx-1)
939  bl(npx) = xt - q1(npx)
940  br(npx) = s11*(q1(npx+1)-q1(npx)) - s14*dm(npx+1)
941  CALL pert_ppm(3, q1(npx-2:npx), bl(npx-2:npx), br(npx-2:npx)&
942 & , 1)
943  END IF
944  END IF
945  DO i=is,ie+1
946  IF (c(i, j) .GT. 0.) THEN
947  flux(i, j) = q1(i-1) + (1.-c(i, j))*(br(i-1)-c(i, j)*(bl(i-1&
948 & )+br(i-1)))
949  ELSE
950  flux(i, j) = q1(i) + (1.+c(i, j))*(bl(i)+c(i, j)*(bl(i)+br(i&
951 & )))
952  END IF
953  END DO
954  END IF
955  END DO
956  END SUBROUTINE xppm
957  SUBROUTINE yppm(flux, q, c, jord, ifirst, ilast, isd, ied, js, je, jsd&
958 & , jed, npx, npy, dya, nested, grid_type)
959  IMPLICIT NONE
960 ! Compute domain
961  INTEGER, INTENT(IN) :: ifirst, ilast
962  INTEGER, INTENT(IN) :: isd, ied, js, je, jsd, jed
963  INTEGER, INTENT(IN) :: jord
964  INTEGER, INTENT(IN) :: npx, npy
965  REAL, INTENT(IN) :: q(ifirst:ilast, jsd:jed)
966 ! Courant number
967  REAL, INTENT(IN) :: c(isd:ied, js:je+1)
968 ! Flux
969  REAL, INTENT(OUT) :: flux(ifirst:ilast, js:je+1)
970  REAL, INTENT(IN) :: dya(isd:ied, jsd:jed)
971  LOGICAL, INTENT(IN) :: nested
972  INTEGER, INTENT(IN) :: grid_type
973 ! Local:
974  REAL :: dm(ifirst:ilast, js-2:je+2)
975  REAL :: al(ifirst:ilast, js-1:je+2)
976  REAL, DIMENSION(ifirst:ilast, js-1:je+1) :: bl, br, b0
977  REAL :: dq(ifirst:ilast, js-3:je+2)
978  REAL, DIMENSION(ifirst:ilast) :: fx0, fx1
979  LOGICAL, DIMENSION(ifirst:ilast, js-1:je+1) :: smt5, smt6
980  REAL :: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
981  INTEGER :: i, j, js1, je3, je1
982  INTRINSIC max
983  INTRINSIC min
984  INTRINSIC abs
985  INTRINSIC sign
986  REAL :: min1
987  REAL :: min2
988  REAL :: abs0
989  REAL :: abs1
990  REAL :: min3
991  REAL :: min4
992  REAL :: min5
993  REAL :: min6
994  REAL :: min7
995  REAL :: abs2
996  REAL :: abs3
997  REAL :: abs4
998  REAL :: max1
999  REAL :: min8
1000  REAL :: abs5
1001  REAL :: abs6
1002  REAL :: abs7
1003  REAL :: x9
1004  REAL :: x8
1005  REAL :: x7
1006  REAL :: x6
1007  REAL :: x5
1008  REAL :: x4
1009  REAL :: x3
1010  REAL :: x2
1011  REAL :: x1
1012  REAL :: y15
1013  REAL :: y14
1014  REAL :: y13
1015  REAL :: y12
1016  REAL :: y11
1017  REAL :: y10
1018  REAL :: z5
1019  REAL :: z4
1020  REAL :: z3
1021  REAL :: z2
1022  REAL :: z1
1023  REAL :: y9
1024  REAL :: y8
1025  REAL :: y7
1026  REAL :: y6
1027  REAL :: y5
1028  REAL :: y4
1029  REAL :: y3
1030  REAL :: y2
1031  REAL :: y1
1032  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
1033  IF (3 .LT. js - 1) THEN
1034  js1 = js - 1
1035  ELSE
1036  js1 = 3
1037  END IF
1038  IF (npy - 2 .GT. je + 2) THEN
1039  je3 = je + 2
1040  ELSE
1041  je3 = npy - 2
1042  END IF
1043  IF (npy - 3 .GT. je + 1) THEN
1044  je1 = je + 1
1045  ELSE
1046  je1 = npy - 3
1047  END IF
1048  ELSE
1049 ! Nested grid OR Doubly periodic domain:
1050  js1 = js - 1
1051  je3 = je + 2
1052  je1 = je + 1
1053  END IF
1054  IF (jord .LT. 8 .OR. jord .EQ. 333) THEN
1055  DO j=js1,je3
1056  DO i=ifirst,ilast
1057  al(i, j) = p1*(q(i, j-1)+q(i, j)) + p2*(q(i, j-2)+q(i, j+1))
1058  END DO
1059  END DO
1060  IF (jord .EQ. 7) THEN
1061  DO j=js1,je3
1062  DO i=ifirst,ilast
1063  IF (al(i, j) .LT. 0.) al(i, j) = 0.5*(q(i, j)+q(i, j+1))
1064  END DO
1065  END DO
1066  END IF
1067  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
1068  IF (js .EQ. 1) THEN
1069  DO i=ifirst,ilast
1070  al(i, 0) = c1*q(i, -2) + c2*q(i, -1) + c3*q(i, 0)
1071  al(i, 1) = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)&
1072 & *q(i, -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2)&
1073 & )*q(i, 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
1074  al(i, 2) = c3*q(i, 1) + c2*q(i, 2) + c1*q(i, 3)
1075  END DO
1076  IF (jord .EQ. 7) THEN
1077  DO i=ifirst,ilast
1078  IF (0. .LT. al(i, 0)) THEN
1079  al(i, 0) = al(i, 0)
1080  ELSE
1081  al(i, 0) = 0.
1082  END IF
1083  IF (0. .LT. al(i, 1)) THEN
1084  al(i, 1) = al(i, 1)
1085  ELSE
1086  al(i, 1) = 0.
1087  END IF
1088  IF (0. .LT. al(i, 2)) THEN
1089  al(i, 2) = al(i, 2)
1090  ELSE
1091  al(i, 2) = 0.
1092  END IF
1093  END DO
1094  END IF
1095  END IF
1096  IF (je + 1 .EQ. npy) THEN
1097  DO i=ifirst,ilast
1098  al(i, npy-1) = c1*q(i, npy-3) + c2*q(i, npy-2) + c3*q(i, npy&
1099 & -1)
1100  al(i, npy) = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy&
1101 & -1)-dya(i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1&
1102 & ))+((2.*dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q&
1103 & (i, npy+1))/(dya(i, npy)+dya(i, npy+1)))
1104  al(i, npy+1) = c3*q(i, npy) + c2*q(i, npy+1) + c1*q(i, npy+2&
1105 & )
1106  END DO
1107  IF (jord .EQ. 7) THEN
1108  DO i=ifirst,ilast
1109  IF (0. .LT. al(i, npy-1)) THEN
1110  al(i, npy-1) = al(i, npy-1)
1111  ELSE
1112  al(i, npy-1) = 0.
1113  END IF
1114  IF (0. .LT. al(i, npy)) THEN
1115  al(i, npy) = al(i, npy)
1116  ELSE
1117  al(i, npy) = 0.
1118  END IF
1119  IF (0. .LT. al(i, npy+1)) THEN
1120  al(i, npy+1) = al(i, npy+1)
1121  ELSE
1122  al(i, npy+1) = 0.
1123  END IF
1124  END DO
1125  END IF
1126  END IF
1127  END IF
1128  IF (jord .EQ. 1) THEN
1129  DO j=js,je+1
1130  DO i=ifirst,ilast
1131  IF (c(i, j) .GT. 0.) THEN
1132  flux(i, j) = q(i, j-1)
1133  ELSE
1134  flux(i, j) = q(i, j)
1135  END IF
1136  END DO
1137  END DO
1138  ELSE IF (jord .EQ. 2) THEN
1139 ! Perfectly linear scheme
1140 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7
1141  DO j=js,je+1
1142 !DEC$ VECTOR ALWAYS
1143  DO i=ifirst,ilast
1144  xt = c(i, j)
1145  IF (xt .GT. 0.) THEN
1146  qtmp = q(i, j-1)
1147  flux(i, j) = qtmp + (1.-xt)*(al(i, j)-qtmp-xt*(al(i, j-1)+&
1148 & al(i, j)-(qtmp+qtmp)))
1149  ELSE
1150  qtmp = q(i, j)
1151  flux(i, j) = qtmp + (1.+xt)*(al(i, j)-qtmp+xt*(al(i, j)+al&
1152 & (i, j+1)-(qtmp+qtmp)))
1153  END IF
1154  END DO
1155  END DO
1156  ELSE IF (jord .EQ. 333) THEN
1157 ! Perfectly linear scheme, more diffusive than ord=2 (HoldawayKent-2015-TellusA)
1158  DO j=js,je+1
1159 !DEC$ VECTOR ALWAYS
1160  DO i=ifirst,ilast
1161  xt = c(i, j)
1162  IF (xt .GT. 0.) THEN
1163  flux(i, j) = (2.0*q(i, j)+5.0*q(i, j-1)-q(i, j-2))/6.0 - &
1164 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j)-2.0*q(i&
1165 & , j-1)+q(i, j-2))
1166  ELSE
1167  flux(i, j) = (2.0*q(i, j-1)+5.0*q(i, j)-q(i, j+1))/6.0 - &
1168 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j+1)-2.0*q(&
1169 & i, j)+q(i, j-1))
1170  END IF
1171  END DO
1172  END DO
1173  ELSE IF (jord .EQ. 3) THEN
1174  DO j=js-1,je+1
1175  DO i=ifirst,ilast
1176  bl(i, j) = al(i, j) - q(i, j)
1177  br(i, j) = al(i, j+1) - q(i, j)
1178  b0(i, j) = bl(i, j) + br(i, j)
1179  IF (b0(i, j) .GE. 0.) THEN
1180  x0 = b0(i, j)
1181  ELSE
1182  x0 = -b0(i, j)
1183  END IF
1184  IF (bl(i, j) - br(i, j) .GE. 0.) THEN
1185  xt = bl(i, j) - br(i, j)
1186  ELSE
1187  xt = -(bl(i, j)-br(i, j))
1188  END IF
1189  smt5(i, j) = x0 .LT. xt
1190  smt6(i, j) = 3.*x0 .LT. xt
1191  END DO
1192  END DO
1193  DO j=js,je+1
1194  DO i=ifirst,ilast
1195  fx1(i) = 0.
1196  END DO
1197  DO i=ifirst,ilast
1198  xt = c(i, j)
1199  IF (xt .GT. 0.) THEN
1200  fx0(i) = q(i, j-1)
1201  IF (smt6(i, j-1) .OR. smt5(i, j)) THEN
1202  fx1(i) = br(i, j-1) - xt*b0(i, j-1)
1203  ELSE IF (smt5(i, j-1)) THEN
1204  IF (bl(i, j-1) .GE. 0.) THEN
1205  x1 = bl(i, j-1)
1206  ELSE
1207  x1 = -bl(i, j-1)
1208  END IF
1209  IF (br(i, j-1) .GE. 0.) THEN
1210  y1 = br(i, j-1)
1211  ELSE
1212  y1 = -br(i, j-1)
1213  END IF
1214  IF (x1 .GT. y1) THEN
1215  min1 = y1
1216  ELSE
1217  min1 = x1
1218  END IF
1219 ! both up-downwind sides are noisy; 2nd order, piece-wise linear
1220  fx1(i) = sign(min1, br(i, j-1))
1221  END IF
1222  ELSE
1223  fx0(i) = q(i, j)
1224  IF (smt6(i, j) .OR. smt5(i, j-1)) THEN
1225  fx1(i) = bl(i, j) + xt*b0(i, j)
1226  ELSE IF (smt5(i, j)) THEN
1227  IF (bl(i, j) .GE. 0.) THEN
1228  x2 = bl(i, j)
1229  ELSE
1230  x2 = -bl(i, j)
1231  END IF
1232  IF (br(i, j) .GE. 0.) THEN
1233  y2 = br(i, j)
1234  ELSE
1235  y2 = -br(i, j)
1236  END IF
1237  IF (x2 .GT. y2) THEN
1238  min2 = y2
1239  ELSE
1240  min2 = x2
1241  END IF
1242  fx1(i) = sign(min2, bl(i, j))
1243  END IF
1244  END IF
1245  IF (xt .GE. 0.) THEN
1246  abs0 = xt
1247  ELSE
1248  abs0 = -xt
1249  END IF
1250  flux(i, j) = fx0(i) + (1.-abs0)*fx1(i)
1251  END DO
1252  END DO
1253  ELSE IF (jord .EQ. 4) THEN
1254  DO j=js-1,je+1
1255  DO i=ifirst,ilast
1256  bl(i, j) = al(i, j) - q(i, j)
1257  br(i, j) = al(i, j+1) - q(i, j)
1258  b0(i, j) = bl(i, j) + br(i, j)
1259  IF (b0(i, j) .GE. 0.) THEN
1260  x0 = b0(i, j)
1261  ELSE
1262  x0 = -b0(i, j)
1263  END IF
1264  IF (bl(i, j) - br(i, j) .GE. 0.) THEN
1265  xt = bl(i, j) - br(i, j)
1266  ELSE
1267  xt = -(bl(i, j)-br(i, j))
1268  END IF
1269  smt5(i, j) = x0 .LT. xt
1270  smt6(i, j) = 3.*x0 .LT. xt
1271  END DO
1272  END DO
1273  DO j=js,je+1
1274  DO i=ifirst,ilast
1275  fx1(i) = 0.
1276  END DO
1277 !DEC$ VECTOR ALWAYS
1278  DO i=ifirst,ilast
1279  IF (c(i, j) .GT. 0.) THEN
1280  fx0(i) = q(i, j-1)
1281  IF (smt6(i, j-1) .OR. smt5(i, j)) fx1(i) = (1.-c(i, j))*(&
1282 & br(i, j-1)-c(i, j)*b0(i, j-1))
1283  ELSE
1284  fx0(i) = q(i, j)
1285  IF (smt6(i, j) .OR. smt5(i, j-1)) fx1(i) = (1.+c(i, j))*(&
1286 & bl(i, j)+c(i, j)*b0(i, j))
1287  END IF
1288  flux(i, j) = fx0(i) + fx1(i)
1289  END DO
1290  END DO
1291  ELSE
1292 ! jord=5,6,7
1293  IF (jord .EQ. 5) THEN
1294  DO j=js-1,je+1
1295  DO i=ifirst,ilast
1296  bl(i, j) = al(i, j) - q(i, j)
1297  br(i, j) = al(i, j+1) - q(i, j)
1298  b0(i, j) = bl(i, j) + br(i, j)
1299  smt5(i, j) = bl(i, j)*br(i, j) .LT. 0.
1300  END DO
1301  END DO
1302  ELSE
1303  DO j=js-1,je+1
1304  DO i=ifirst,ilast
1305  bl(i, j) = al(i, j) - q(i, j)
1306  br(i, j) = al(i, j+1) - q(i, j)
1307  b0(i, j) = bl(i, j) + br(i, j)
1308  IF (3.*b0(i, j) .GE. 0.) THEN
1309  abs1 = 3.*b0(i, j)
1310  ELSE
1311  abs1 = -(3.*b0(i, j))
1312  END IF
1313  IF (bl(i, j) - br(i, j) .GE. 0.) THEN
1314  abs4 = bl(i, j) - br(i, j)
1315  ELSE
1316  abs4 = -(bl(i, j)-br(i, j))
1317  END IF
1318  smt5(i, j) = abs1 .LT. abs4
1319  END DO
1320  END DO
1321  END IF
1322  DO j=js,je+1
1323 !DEC$ VECTOR ALWAYS
1324  DO i=ifirst,ilast
1325  IF (c(i, j) .GT. 0.) THEN
1326  fx1(i) = (1.-c(i, j))*(br(i, j-1)-c(i, j)*b0(i, j-1))
1327  flux(i, j) = q(i, j-1)
1328  ELSE
1329  fx1(i) = (1.+c(i, j))*(bl(i, j)+c(i, j)*b0(i, j))
1330  flux(i, j) = q(i, j)
1331  END IF
1332  IF (smt5(i, j-1) .OR. smt5(i, j)) flux(i, j) = flux(i, j) + &
1333 & fx1(i)
1334  END DO
1335  END DO
1336  END IF
1337  RETURN
1338  ELSE
1339 ! Monotonic constraints:
1340 ! ord = 8: PPM with Lin's PPM fast monotone constraint
1341 ! ord > 8: PPM with Lin's modification of Huynh 2nd constraint
1342  DO j=js-2,je+2
1343  DO i=ifirst,ilast
1344  xt = 0.25*(q(i, j+1)-q(i, j-1))
1345  IF (xt .GE. 0.) THEN
1346  x3 = xt
1347  ELSE
1348  x3 = -xt
1349  END IF
1350  IF (q(i, j-1) .LT. q(i, j)) THEN
1351  IF (q(i, j) .LT. q(i, j+1)) THEN
1352  max1 = q(i, j+1)
1353  ELSE
1354  max1 = q(i, j)
1355  END IF
1356  ELSE IF (q(i, j-1) .LT. q(i, j+1)) THEN
1357  max1 = q(i, j+1)
1358  ELSE
1359  max1 = q(i, j-1)
1360  END IF
1361  y3 = max1 - q(i, j)
1362  IF (q(i, j-1) .GT. q(i, j)) THEN
1363  IF (q(i, j) .GT. q(i, j+1)) THEN
1364  min8 = q(i, j+1)
1365  ELSE
1366  min8 = q(i, j)
1367  END IF
1368  ELSE IF (q(i, j-1) .GT. q(i, j+1)) THEN
1369  min8 = q(i, j+1)
1370  ELSE
1371  min8 = q(i, j-1)
1372  END IF
1373  z1 = q(i, j) - min8
1374  IF (x3 .GT. y3) THEN
1375  IF (y3 .GT. z1) THEN
1376  min3 = z1
1377  ELSE
1378  min3 = y3
1379  END IF
1380  ELSE IF (x3 .GT. z1) THEN
1381  min3 = z1
1382  ELSE
1383  min3 = x3
1384  END IF
1385  dm(i, j) = sign(min3, xt)
1386  END DO
1387  END DO
1388  DO j=js1,je1+1
1389  DO i=ifirst,ilast
1390  al(i, j) = 0.5*(q(i, j-1)+q(i, j)) + r3*(dm(i, j-1)-dm(i, j))
1391  END DO
1392  END DO
1393  IF (jord .EQ. 8) THEN
1394  DO j=js1,je1
1395  DO i=ifirst,ilast
1396  xt = 2.*dm(i, j)
1397  IF (xt .GE. 0.) THEN
1398  x4 = xt
1399  ELSE
1400  x4 = -xt
1401  END IF
1402  IF (al(i, j) - q(i, j) .GE. 0.) THEN
1403  y4 = al(i, j) - q(i, j)
1404  ELSE
1405  y4 = -(al(i, j)-q(i, j))
1406  END IF
1407  IF (x4 .GT. y4) THEN
1408  min4 = y4
1409  ELSE
1410  min4 = x4
1411  END IF
1412  bl(i, j) = -sign(min4, xt)
1413  IF (xt .GE. 0.) THEN
1414  x5 = xt
1415  ELSE
1416  x5 = -xt
1417  END IF
1418  IF (al(i, j+1) - q(i, j) .GE. 0.) THEN
1419  y5 = al(i, j+1) - q(i, j)
1420  ELSE
1421  y5 = -(al(i, j+1)-q(i, j))
1422  END IF
1423  IF (x5 .GT. y5) THEN
1424  min5 = y5
1425  ELSE
1426  min5 = x5
1427  END IF
1428  br(i, j) = sign(min5, xt)
1429  END DO
1430  END DO
1431  ELSE IF (jord .EQ. 11) THEN
1432  DO j=js1,je1
1433  DO i=ifirst,ilast
1434  xt = ppm_fac*dm(i, j)
1435  IF (xt .GE. 0.) THEN
1436  x6 = xt
1437  ELSE
1438  x6 = -xt
1439  END IF
1440  IF (al(i, j) - q(i, j) .GE. 0.) THEN
1441  y6 = al(i, j) - q(i, j)
1442  ELSE
1443  y6 = -(al(i, j)-q(i, j))
1444  END IF
1445  IF (x6 .GT. y6) THEN
1446  min6 = y6
1447  ELSE
1448  min6 = x6
1449  END IF
1450  bl(i, j) = -sign(min6, xt)
1451  IF (xt .GE. 0.) THEN
1452  x7 = xt
1453  ELSE
1454  x7 = -xt
1455  END IF
1456  IF (al(i, j+1) - q(i, j) .GE. 0.) THEN
1457  y7 = al(i, j+1) - q(i, j)
1458  ELSE
1459  y7 = -(al(i, j+1)-q(i, j))
1460  END IF
1461  IF (x7 .GT. y7) THEN
1462  min7 = y7
1463  ELSE
1464  min7 = x7
1465  END IF
1466  br(i, j) = sign(min7, xt)
1467  END DO
1468  END DO
1469  ELSE
1470  DO j=js1-2,je1+1
1471  DO i=ifirst,ilast
1472  dq(i, j) = 2.*(q(i, j+1)-q(i, j))
1473  END DO
1474  END DO
1475  DO j=js1,je1
1476  DO i=ifirst,ilast
1477  bl(i, j) = al(i, j) - q(i, j)
1478  br(i, j) = al(i, j+1) - q(i, j)
1479  IF (dm(i, j-1) .GE. 0.) THEN
1480  abs2 = dm(i, j-1)
1481  ELSE
1482  abs2 = -dm(i, j-1)
1483  END IF
1484  IF (dm(i, j) .GE. 0.) THEN
1485  abs5 = dm(i, j)
1486  ELSE
1487  abs5 = -dm(i, j)
1488  END IF
1489  IF (dm(i, j+1) .GE. 0.) THEN
1490  abs7 = dm(i, j+1)
1491  ELSE
1492  abs7 = -dm(i, j+1)
1493  END IF
1494  IF (abs2 + abs5 + abs7 .LT. near_zero) THEN
1495  bl(i, j) = 0.
1496  br(i, j) = 0.
1497  ELSE
1498  IF (3.*(bl(i, j)+br(i, j)) .GE. 0.) THEN
1499  abs3 = 3.*(bl(i, j)+br(i, j))
1500  ELSE
1501  abs3 = -(3.*(bl(i, j)+br(i, j)))
1502  END IF
1503  IF (bl(i, j) - br(i, j) .GE. 0.) THEN
1504  abs6 = bl(i, j) - br(i, j)
1505  ELSE
1506  abs6 = -(bl(i, j)-br(i, j))
1507  END IF
1508  IF (abs3 .GT. abs6) THEN
1509  pmp_2 = dq(i, j-1)
1510  lac_2 = pmp_2 - 0.75*dq(i, j-2)
1511  IF (0. .LT. pmp_2) THEN
1512  IF (pmp_2 .LT. lac_2) THEN
1513  x8 = lac_2
1514  ELSE
1515  x8 = pmp_2
1516  END IF
1517  ELSE IF (0. .LT. lac_2) THEN
1518  x8 = lac_2
1519  ELSE
1520  x8 = 0.
1521  END IF
1522  IF (0. .GT. pmp_2) THEN
1523  IF (pmp_2 .GT. lac_2) THEN
1524  y14 = lac_2
1525  ELSE
1526  y14 = pmp_2
1527  END IF
1528  ELSE IF (0. .GT. lac_2) THEN
1529  y14 = lac_2
1530  ELSE
1531  y14 = 0.
1532  END IF
1533  IF (br(i, j) .LT. y14) THEN
1534  y8 = y14
1535  ELSE
1536  y8 = br(i, j)
1537  END IF
1538  IF (x8 .GT. y8) THEN
1539  br(i, j) = y8
1540  ELSE
1541  br(i, j) = x8
1542  END IF
1543  pmp_1 = -dq(i, j)
1544  lac_1 = pmp_1 + 0.75*dq(i, j+1)
1545  IF (0. .LT. pmp_1) THEN
1546  IF (pmp_1 .LT. lac_1) THEN
1547  x9 = lac_1
1548  ELSE
1549  x9 = pmp_1
1550  END IF
1551  ELSE IF (0. .LT. lac_1) THEN
1552  x9 = lac_1
1553  ELSE
1554  x9 = 0.
1555  END IF
1556  IF (0. .GT. pmp_1) THEN
1557  IF (pmp_1 .GT. lac_1) THEN
1558  y15 = lac_1
1559  ELSE
1560  y15 = pmp_1
1561  END IF
1562  ELSE IF (0. .GT. lac_1) THEN
1563  y15 = lac_1
1564  ELSE
1565  y15 = 0.
1566  END IF
1567  IF (bl(i, j) .LT. y15) THEN
1568  y9 = y15
1569  ELSE
1570  y9 = bl(i, j)
1571  END IF
1572  IF (x9 .GT. y9) THEN
1573  bl(i, j) = y9
1574  ELSE
1575  bl(i, j) = x9
1576  END IF
1577  END IF
1578  END IF
1579  END DO
1580  END DO
1581  END IF
1582  IF (jord .EQ. 9 .OR. jord .EQ. 13) THEN
1583 ! Positive definite constraint:
1584  DO j=js1,je1
1585  CALL pert_ppm(ilast - ifirst + 1, q(ifirst:ilast, j), bl(&
1586 & ifirst:ilast, j), br(ifirst:ilast, j), 0)
1587  END DO
1588  END IF
1589  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
1590  IF (js .EQ. 1) THEN
1591  DO i=ifirst,ilast
1592  bl(i, 0) = s14*dm(i, -1) + s11*(q(i, -1)-q(i, 0))
1593  xt = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)*q(i, &
1594 & -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2))*q(i&
1595 & , 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
1596  IF (q(i, 1) .GT. q(i, 2)) THEN
1597  z2 = q(i, 2)
1598  ELSE
1599  z2 = q(i, 1)
1600  END IF
1601  IF (q(i, -1) .GT. q(i, 0)) THEN
1602  IF (q(i, 0) .GT. z2) THEN
1603  y10 = z2
1604  ELSE
1605  y10 = q(i, 0)
1606  END IF
1607  ELSE IF (q(i, -1) .GT. z2) THEN
1608  y10 = z2
1609  ELSE
1610  y10 = q(i, -1)
1611  END IF
1612  IF (xt .LT. y10) THEN
1613  xt = y10
1614  ELSE
1615  xt = xt
1616  END IF
1617  IF (q(i, 1) .LT. q(i, 2)) THEN
1618  z3 = q(i, 2)
1619  ELSE
1620  z3 = q(i, 1)
1621  END IF
1622  IF (q(i, -1) .LT. q(i, 0)) THEN
1623  IF (q(i, 0) .LT. z3) THEN
1624  y11 = z3
1625  ELSE
1626  y11 = q(i, 0)
1627  END IF
1628  ELSE IF (q(i, -1) .LT. z3) THEN
1629  y11 = z3
1630  ELSE
1631  y11 = q(i, -1)
1632  END IF
1633  IF (xt .GT. y11) THEN
1634  xt = y11
1635  ELSE
1636  xt = xt
1637  END IF
1638 ! endif
1639  br(i, 0) = xt - q(i, 0)
1640  bl(i, 1) = xt - q(i, 1)
1641  xt = s15*q(i, 1) + s11*q(i, 2) - s14*dm(i, 2)
1642  br(i, 1) = xt - q(i, 1)
1643  bl(i, 2) = xt - q(i, 2)
1644  br(i, 2) = al(i, 3) - q(i, 2)
1645  END DO
1646  CALL pert_ppm(3*(ilast-ifirst+1), q(ifirst:ilast, 0:2), bl(&
1647 & ifirst:ilast, 0:2), br(ifirst:ilast, 0:2), 1)
1648  END IF
1649  IF (je + 1 .EQ. npy) THEN
1650  DO i=ifirst,ilast
1651  bl(i, npy-2) = al(i, npy-2) - q(i, npy-2)
1652  xt = s15*q(i, npy-1) + s11*q(i, npy-2) + s14*dm(i, npy-2)
1653  br(i, npy-2) = xt - q(i, npy-2)
1654  bl(i, npy-1) = xt - q(i, npy-1)
1655  xt = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy-1)-dya(&
1656 & i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1))+((2.*&
1657 & dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q(i, npy+&
1658 & 1))/(dya(i, npy)+dya(i, npy+1)))
1659  IF (q(i, npy) .GT. q(i, npy+1)) THEN
1660  z4 = q(i, npy+1)
1661  ELSE
1662  z4 = q(i, npy)
1663  END IF
1664  IF (q(i, npy-2) .GT. q(i, npy-1)) THEN
1665  IF (q(i, npy-1) .GT. z4) THEN
1666  y12 = z4
1667  ELSE
1668  y12 = q(i, npy-1)
1669  END IF
1670  ELSE IF (q(i, npy-2) .GT. z4) THEN
1671  y12 = z4
1672  ELSE
1673  y12 = q(i, npy-2)
1674  END IF
1675  IF (xt .LT. y12) THEN
1676  xt = y12
1677  ELSE
1678  xt = xt
1679  END IF
1680  IF (q(i, npy) .LT. q(i, npy+1)) THEN
1681  z5 = q(i, npy+1)
1682  ELSE
1683  z5 = q(i, npy)
1684  END IF
1685  IF (q(i, npy-2) .LT. q(i, npy-1)) THEN
1686  IF (q(i, npy-1) .LT. z5) THEN
1687  y13 = z5
1688  ELSE
1689  y13 = q(i, npy-1)
1690  END IF
1691  ELSE IF (q(i, npy-2) .LT. z5) THEN
1692  y13 = z5
1693  ELSE
1694  y13 = q(i, npy-2)
1695  END IF
1696  IF (xt .GT. y13) THEN
1697  xt = y13
1698  ELSE
1699  xt = xt
1700  END IF
1701 ! endif
1702  br(i, npy-1) = xt - q(i, npy-1)
1703  bl(i, npy) = xt - q(i, npy)
1704  br(i, npy) = s11*(q(i, npy+1)-q(i, npy)) - s14*dm(i, npy+1)
1705  END DO
1706  CALL pert_ppm(3*(ilast-ifirst+1), q(ifirst:ilast, npy-2:npy), &
1707 & bl(ifirst:ilast, npy-2:npy), br(ifirst:ilast, npy-2:&
1708 & npy), 1)
1709  END IF
1710  END IF
1711  DO j=js,je+1
1712  DO i=ifirst,ilast
1713  IF (c(i, j) .GT. 0.) THEN
1714  flux(i, j) = q(i, j-1) + (1.-c(i, j))*(br(i, j-1)-c(i, j)*(&
1715 & bl(i, j-1)+br(i, j-1)))
1716  ELSE
1717  flux(i, j) = q(i, j) + (1.+c(i, j))*(bl(i, j)+c(i, j)*(bl(i&
1718 & , j)+br(i, j)))
1719  END IF
1720  END DO
1721  END DO
1722  END IF
1723  END SUBROUTINE yppm
1724  SUBROUTINE mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
1725 & kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
1726  IMPLICIT NONE
1727 !
1728 ! !INPUT PARAMETERS:
1729  INTEGER, INTENT(IN) :: im, jm, km, nq
1730  INTEGER, INTENT(IN) :: ifirst, ilast
1731  INTEGER, INTENT(IN) :: jfirst, jlast
1732  INTEGER, INTENT(IN) :: kfirst, klast
1733 ! eastern zones to ghost
1734  INTEGER, INTENT(IN) :: ng_e
1735 ! western zones to ghost
1736  INTEGER, INTENT(IN) :: ng_w
1737 ! southern zones to ghost
1738  INTEGER, INTENT(IN) :: ng_s
1739 ! northern zones to ghost
1740  INTEGER, INTENT(IN) :: ng_n
1741  REAL, INTENT(INOUT) :: q_ghst(ifirst-ng_w:ilast+ng_e, jfirst-ng_s:&
1742 & jlast+ng_n, kfirst:klast, nq)
1743  REAL, OPTIONAL, INTENT(IN) :: q(ifirst:ilast, jfirst:jlast, kfirst:&
1744 & klast, nq)
1745 !
1746 ! !DESCRIPTION:
1747 !
1748 ! Ghost 4d east/west
1749 !
1750 ! !REVISION HISTORY:
1751 ! 2005.08.22 Putman
1752 !
1753 !EOP
1754 !------------------------------------------------------------------------------
1755 !BOC
1756  INTEGER :: i, j, k, n
1757  INTRINSIC PRESENT
1758  IF (PRESENT(q)) q_ghst(ifirst:ilast, jfirst:jlast, kfirst:klast, 1:&
1759 & nq) = q(ifirst:ilast, jfirst:jlast, kfirst:klast, 1:nq)
1760 ! Assume Periodicity in X-dir and not overlapping
1761  DO n=1,nq
1762  DO k=kfirst,klast
1763  DO j=jfirst-ng_s,jlast+ng_n
1764  DO i=1,ng_w
1765  q_ghst(ifirst-i, j, k, n) = q_ghst(ilast-i+1, j, k, n)
1766  END DO
1767  DO i=1,ng_e
1768  q_ghst(ilast+i, j, k, n) = q_ghst(ifirst+i-1, j, k, n)
1769  END DO
1770  END DO
1771  END DO
1772  END DO
1773  END SUBROUTINE mp_ghost_ew
1774 ! Differentiation of pert_ppm in forward (tangent) mode:
1775 ! variations of useful results: al ar
1776 ! with respect to varying inputs: al ar
1777  SUBROUTINE pert_ppm_tlm(im, a0, al, al_tl, ar, ar_tl, iv)
1778  IMPLICIT NONE
1779  INTEGER, INTENT(IN) :: im
1780  INTEGER, INTENT(IN) :: iv
1781  REAL, INTENT(IN) :: a0(im)
1782  REAL, INTENT(INOUT) :: al(im), ar(im)
1783  REAL, INTENT(INOUT) :: al_tl(im), ar_tl(im)
1784 ! Local:
1785  REAL :: a4, da1, da2, a6da, fmin
1786  INTEGER :: i
1787  REAL, PARAMETER :: r12=1./12.
1788  INTRINSIC abs
1789  REAL :: abs0
1790 !-----------------------------------
1791 ! Optimized PPM in perturbation form:
1792 !-----------------------------------
1793  IF (iv .EQ. 0) THEN
1794 ! Positive definite constraint
1795  DO i=1,im
1796  IF (a0(i) .LE. 0.) THEN
1797  al_tl(i) = 0.0
1798  al(i) = 0.
1799  ar_tl(i) = 0.0
1800  ar(i) = 0.
1801  ELSE
1802  a4 = -(3.*(ar(i)+al(i)))
1803  da1 = ar(i) - al(i)
1804  IF (da1 .GE. 0.) THEN
1805  abs0 = da1
1806  ELSE
1807  abs0 = -da1
1808  END IF
1809  IF (abs0 .LT. -a4) THEN
1810  fmin = a0(i) + 0.25/a4*da1**2 + a4*r12
1811  IF (fmin .LT. 0.) THEN
1812  IF (ar(i) .GT. 0. .AND. al(i) .GT. 0.) THEN
1813  ar_tl(i) = 0.0
1814  ar(i) = 0.
1815  al_tl(i) = 0.0
1816  al(i) = 0.
1817  ELSE IF (da1 .GT. 0.) THEN
1818  ar_tl(i) = -(2.*al_tl(i))
1819  ar(i) = -(2.*al(i))
1820  ELSE
1821  al_tl(i) = -(2.*ar_tl(i))
1822  al(i) = -(2.*ar(i))
1823  END IF
1824  END IF
1825  END IF
1826  END IF
1827  END DO
1828  ELSE
1829 ! Standard PPM constraint
1830  DO i=1,im
1831  IF (al(i)*ar(i) .LT. 0.) THEN
1832  da1 = al(i) - ar(i)
1833  da2 = da1**2
1834  a6da = 3.*(al(i)+ar(i))*da1
1835 ! abs(a6da) > da2 --> 3.*abs(al+ar) > abs(al-ar)
1836  IF (a6da .LT. -da2) THEN
1837  ar_tl(i) = -(2.*al_tl(i))
1838  ar(i) = -(2.*al(i))
1839  ELSE IF (a6da .GT. da2) THEN
1840  al_tl(i) = -(2.*ar_tl(i))
1841  al(i) = -(2.*ar(i))
1842  END IF
1843  ELSE
1844 ! effect of dm=0 included here
1845  al_tl(i) = 0.0
1846  al(i) = 0.
1847  ar_tl(i) = 0.0
1848  ar(i) = 0.
1849  END IF
1850  END DO
1851  END IF
1852  END SUBROUTINE pert_ppm_tlm
1853  SUBROUTINE pert_ppm(im, a0, al, ar, iv)
1854  IMPLICIT NONE
1855  INTEGER, INTENT(IN) :: im
1856  INTEGER, INTENT(IN) :: iv
1857  REAL, INTENT(IN) :: a0(im)
1858  REAL, INTENT(INOUT) :: al(im), ar(im)
1859 ! Local:
1860  REAL :: a4, da1, da2, a6da, fmin
1861  INTEGER :: i
1862  REAL, PARAMETER :: r12=1./12.
1863  INTRINSIC abs
1864  REAL :: abs0
1865 !-----------------------------------
1866 ! Optimized PPM in perturbation form:
1867 !-----------------------------------
1868  IF (iv .EQ. 0) THEN
1869 ! Positive definite constraint
1870  DO i=1,im
1871  IF (a0(i) .LE. 0.) THEN
1872  al(i) = 0.
1873  ar(i) = 0.
1874  ELSE
1875  a4 = -(3.*(ar(i)+al(i)))
1876  da1 = ar(i) - al(i)
1877  IF (da1 .GE. 0.) THEN
1878  abs0 = da1
1879  ELSE
1880  abs0 = -da1
1881  END IF
1882  IF (abs0 .LT. -a4) THEN
1883  fmin = a0(i) + 0.25/a4*da1**2 + a4*r12
1884  IF (fmin .LT. 0.) THEN
1885  IF (ar(i) .GT. 0. .AND. al(i) .GT. 0.) THEN
1886  ar(i) = 0.
1887  al(i) = 0.
1888  ELSE IF (da1 .GT. 0.) THEN
1889  ar(i) = -(2.*al(i))
1890  ELSE
1891  al(i) = -(2.*ar(i))
1892  END IF
1893  END IF
1894  END IF
1895  END IF
1896  END DO
1897  ELSE
1898 ! Standard PPM constraint
1899  DO i=1,im
1900  IF (al(i)*ar(i) .LT. 0.) THEN
1901  da1 = al(i) - ar(i)
1902  da2 = da1**2
1903  a6da = 3.*(al(i)+ar(i))*da1
1904 ! abs(a6da) > da2 --> 3.*abs(al+ar) > abs(al-ar)
1905  IF (a6da .LT. -da2) THEN
1906  ar(i) = -(2.*al(i))
1907  ELSE IF (a6da .GT. da2) THEN
1908  al(i) = -(2.*ar(i))
1909  END IF
1910  ELSE
1911 ! effect of dm=0 included here
1912  al(i) = 0.
1913  ar(i) = 0.
1914  END IF
1915  END DO
1916  END IF
1917  END SUBROUTINE pert_ppm
1918  SUBROUTINE deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy, &
1919 & gridstruct, bd, mass)
1920  IMPLICIT NONE
1921 ! Del-n damping for the cell-mean values (A grid)
1922 !------------------
1923 ! nord = 0: del-2
1924 ! nord = 1: del-4
1925 ! nord = 2: del-6
1926 ! nord = 3: del-8 --> requires more ghosting than current
1927 !------------------
1928  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
1929 ! del-n
1930  INTEGER, INTENT(IN) :: nord
1931  INTEGER, INTENT(IN) :: is, ie, js, je, npx, npy
1932  REAL, INTENT(IN) :: damp
1933 ! q ghosted on input
1934  REAL, INTENT(IN) :: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
1935  TYPE(FV_GRID_TYPE), INTENT(IN), TARGET :: gridstruct
1936 ! q ghosted on input
1937  REAL, OPTIONAL, INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
1938 ! diffusive fluxes:
1939  REAL, INTENT(INOUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je), fy(bd%is:bd%&
1940 & ie, bd%js:bd%je+1)
1941 ! local:
1942  REAL :: fx2(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2(bd%isd:bd%ied, bd%&
1943 & jsd:bd%jed+1)
1944  REAL :: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
1945  REAL :: damp2
1946  INTEGER :: i, j, n, nt, i1, i2, j1, j2
1947  INTRINSIC PRESENT
1948  i1 = is - 1 - nord
1949  i2 = ie + 1 + nord
1950  j1 = js - 1 - nord
1951  j2 = je + 1 + nord
1952  IF (.NOT.PRESENT(mass)) THEN
1953  DO j=j1,j2
1954  DO i=i1,i2
1955  d2(i, j) = damp*q(i, j)
1956  END DO
1957  END DO
1958  ELSE
1959  DO j=j1,j2
1960  DO i=i1,i2
1961  d2(i, j) = q(i, j)
1962  END DO
1963  END DO
1964  END IF
1965  IF (nord .GT. 0) CALL copy_corners(d2, npx, npy, 1, gridstruct%&
1966 & nested, bd, gridstruct%sw_corner, &
1967 & gridstruct%se_corner, gridstruct%&
1968 & nw_corner, gridstruct%ne_corner)
1969  DO j=js-nord,je+nord
1970  DO i=is-nord,ie+nord+1
1971  fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i-1, j)-d2(i, j))
1972  END DO
1973  END DO
1974  IF (nord .GT. 0) CALL copy_corners(d2, npx, npy, 2, gridstruct%&
1975 & nested, bd, gridstruct%sw_corner, &
1976 & gridstruct%se_corner, gridstruct%&
1977 & nw_corner, gridstruct%ne_corner)
1978  DO j=js-nord,je+nord+1
1979  DO i=is-nord,ie+nord
1980  fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j-1)-d2(i, j))
1981  END DO
1982  END DO
1983  IF (nord .GT. 0) THEN
1984 !----------
1985 ! high-order
1986 !----------
1987  DO n=1,nord
1988  nt = nord - n
1989  DO j=js-nt-1,je+nt+1
1990  DO i=is-nt-1,ie+nt+1
1991  d2(i, j) = (fx2(i, j)-fx2(i+1, j)+fy2(i, j)-fy2(i, j+1))*&
1992 & gridstruct%rarea(i, j)
1993  END DO
1994  END DO
1995  CALL copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, &
1996 & gridstruct%sw_corner, gridstruct%se_corner, &
1997 & gridstruct%nw_corner, gridstruct%ne_corner)
1998  DO j=js-nt,je+nt
1999  DO i=is-nt,ie+nt+1
2000  fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i, j)-d2(i-1, j))
2001  END DO
2002  END DO
2003  CALL copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, &
2004 & gridstruct%sw_corner, gridstruct%se_corner, &
2005 & gridstruct%nw_corner, gridstruct%ne_corner)
2006  DO j=js-nt,je+nt+1
2007  DO i=is-nt,ie+nt
2008  fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j)-d2(i, j-1))
2009  END DO
2010  END DO
2011  END DO
2012  END IF
2013 !---------------------------------------------
2014 ! Add the diffusive fluxes to the flux arrays:
2015 !---------------------------------------------
2016  IF (PRESENT(mass)) THEN
2017 ! Apply mass weighting to diffusive fluxes:
2018  damp2 = 0.5*damp
2019  DO j=js,je
2020  DO i=is,ie+1
2021  fx(i, j) = fx(i, j) + damp2*(mass(i-1, j)+mass(i, j))*fx2(i, j&
2022 & )
2023  END DO
2024  END DO
2025  DO j=js,je+1
2026  DO i=is,ie
2027  fy(i, j) = fy(i, j) + damp2*(mass(i, j-1)+mass(i, j))*fy2(i, j&
2028 & )
2029  END DO
2030  END DO
2031  ELSE
2032  DO j=js,je
2033  DO i=is,ie+1
2034  fx(i, j) = fx(i, j) + fx2(i, j)
2035  END DO
2036  END DO
2037  DO j=js,je+1
2038  DO i=is,ie
2039  fy(i, j) = fy(i, j) + fy2(i, j)
2040  END DO
2041  END DO
2042  END IF
2043  END SUBROUTINE deln_flux
2044 !Weird arguments are because this routine is called in a lot of
2045 !places outside of tp_core, sometimes very deeply nested in the call tree.
2046  SUBROUTINE copy_corners(q, npx, npy, dir, nested, bd, sw_corner, &
2047 & se_corner, nw_corner, ne_corner)
2048  IMPLICIT NONE
2049  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
2050  INTEGER, INTENT(IN) :: npx, npy, dir
2051  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
2052  LOGICAL, INTENT(IN) :: nested, sw_corner, se_corner, nw_corner, &
2053 & ne_corner
2054  INTEGER :: i, j
2055  IF (nested) THEN
2056  RETURN
2057  ELSE IF (dir .EQ. 1) THEN
2058 ! XDir:
2059  IF (sw_corner) THEN
2060  DO j=1-ng,0
2061  DO i=1-ng,0
2062  q(i, j) = q(j, 1-i)
2063  END DO
2064  END DO
2065  END IF
2066  IF (se_corner) THEN
2067  DO j=1-ng,0
2068  DO i=npx,npx+ng-1
2069  q(i, j) = q(npy-j, i-npx+1)
2070  END DO
2071  END DO
2072  END IF
2073  IF (ne_corner) THEN
2074  DO j=npy,npy+ng-1
2075  DO i=npx,npx+ng-1
2076  q(i, j) = q(j, 2*npx-1-i)
2077  END DO
2078  END DO
2079  END IF
2080  IF (nw_corner) THEN
2081  DO j=npy,npy+ng-1
2082  DO i=1-ng,0
2083  q(i, j) = q(npy-j, i-1+npx)
2084  END DO
2085  END DO
2086  END IF
2087  ELSE IF (dir .EQ. 2) THEN
2088 ! YDir:
2089  IF (sw_corner) THEN
2090  DO j=1-ng,0
2091  DO i=1-ng,0
2092  q(i, j) = q(1-j, i)
2093  END DO
2094  END DO
2095  END IF
2096  IF (se_corner) THEN
2097  DO j=1-ng,0
2098  DO i=npx,npx+ng-1
2099  q(i, j) = q(npy+j-1, npx-i)
2100  END DO
2101  END DO
2102  END IF
2103  IF (ne_corner) THEN
2104  DO j=npy,npy+ng-1
2105  DO i=npx,npx+ng-1
2106  q(i, j) = q(2*npy-1-j, i)
2107  END DO
2108  END DO
2109  END IF
2110  IF (nw_corner) THEN
2111  DO j=npy,npy+ng-1
2112  DO i=1-ng,0
2113  q(i, j) = q(j+1-npx, npy-i)
2114  END DO
2115  END DO
2116  END IF
2117  END IF
2118  END SUBROUTINE copy_corners
2119 ! Differentiation of fv_tp_2d in forward (tangent) mode:
2120 ! variations of useful results: q fx fy
2121 ! with respect to varying inputs: xfx q mass mfx mfy ra_x ra_y
2122 ! yfx fx fy crx cry
2123  SUBROUTINE fv_tp_2d_tlm(q, q_tl, crx, crx_tl, cry, cry_tl, npx, npy&
2124 & , hord, fx, fx_tl, fy, fy_tl, xfx, xfx_tl, yfx, yfx_tl, gridstruct, &
2125 & bd, ra_x, ra_x_tl, ra_y, ra_y_tl, mfx, mfx_tl, mfy, mfy_tl, mass, &
2126 & mass_tl, nord, damp_c)
2127  IMPLICIT NONE
2128  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
2129  INTEGER, INTENT(IN) :: npx, npy
2130  INTEGER, INTENT(IN) :: hord
2131 !
2132  REAL, INTENT(IN) :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed)
2133  REAL, INTENT(IN) :: crx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed)
2134 !
2135  REAL, INTENT(IN) :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed)
2136  REAL, INTENT(IN) :: xfx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed)
2137 !
2138  REAL, INTENT(IN) :: cry(bd%isd:bd%ied, bd%js:bd%je+1)
2139  REAL, INTENT(IN) :: cry_tl(bd%isd:bd%ied, bd%js:bd%je+1)
2140 !
2141  REAL, INTENT(IN) :: yfx(bd%isd:bd%ied, bd%js:bd%je+1)
2142  REAL, INTENT(IN) :: yfx_tl(bd%isd:bd%ied, bd%js:bd%je+1)
2143  REAL, INTENT(IN) :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
2144  REAL, INTENT(IN) :: ra_x_tl(bd%is:bd%ie, bd%jsd:bd%jed)
2145  REAL, INTENT(IN) :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
2146  REAL, INTENT(IN) :: ra_y_tl(bd%isd:bd%ied, bd%js:bd%je)
2147 ! transported scalar
2148  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
2149  REAL, INTENT(INOUT) :: q_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
2150 ! Flux in x ( E )
2151  REAL, INTENT(OUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je)
2152  REAL, INTENT(OUT) :: fx_tl(bd%is:bd%ie+1, bd%js:bd%je)
2153 ! Flux in y ( N )
2154  REAL, INTENT(OUT) :: fy(bd%is:bd%ie, bd%js:bd%je+1)
2155  REAL, INTENT(OUT) :: fy_tl(bd%is:bd%ie, bd%js:bd%je+1)
2156  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
2157 ! optional Arguments:
2158 ! Mass Flux X-Dir
2159  REAL, OPTIONAL, INTENT(IN) :: mfx(bd%is:bd%ie+1, bd%js:bd%je)
2160  REAL, OPTIONAL, INTENT(IN) :: mfx_tl(bd%is:bd%ie+1, bd%js:bd%je)
2161 ! Mass Flux Y-Dir
2162  REAL, OPTIONAL, INTENT(IN) :: mfy(bd%is:bd%ie, bd%js:bd%je+1)
2163  REAL, OPTIONAL, INTENT(IN) :: mfy_tl(bd%is:bd%ie, bd%js:bd%je+1)
2164  REAL, OPTIONAL, INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
2165  REAL, OPTIONAL, INTENT(IN) :: mass_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
2166  REAL, OPTIONAL, INTENT(IN) :: damp_c
2167  INTEGER, OPTIONAL, INTENT(IN) :: nord
2168 ! Local:
2169  INTEGER :: ord_ou, ord_in
2170  REAL :: q_i(bd%isd:bd%ied, bd%js:bd%je)
2171  REAL :: q_i_tl(bd%isd:bd%ied, bd%js:bd%je)
2172  REAL :: q_j(bd%is:bd%ie, bd%jsd:bd%jed)
2173  REAL :: q_j_tl(bd%is:bd%ie, bd%jsd:bd%jed)
2174  REAL :: fx2(bd%is:bd%ie+1, bd%jsd:bd%jed)
2175  REAL :: fx2_tl(bd%is:bd%ie+1, bd%jsd:bd%jed)
2176  REAL :: fy2(bd%isd:bd%ied, bd%js:bd%je+1)
2177  REAL :: fy2_tl(bd%isd:bd%ied, bd%js:bd%je+1)
2178  REAL :: fyy(bd%isd:bd%ied, bd%js:bd%je+1)
2179  REAL :: fyy_tl(bd%isd:bd%ied, bd%js:bd%je+1)
2180  REAL :: fx1(bd%is:bd%ie+1)
2181  REAL :: fx1_tl(bd%is:bd%ie+1)
2182  REAL :: damp
2183  INTEGER :: i, j
2184  INTEGER :: is, ie, js, je, isd, ied, jsd, jed
2185  INTRINSIC PRESENT
2186  REAL*8 :: pwx1
2187  INTEGER :: pwy1
2188  is = bd%is
2189  ie = bd%ie
2190  js = bd%js
2191  je = bd%je
2192  isd = bd%isd
2193  ied = bd%ied
2194  jsd = bd%jsd
2195  jed = bd%jed
2196  IF (hord .EQ. 10) THEN
2197  ord_in = 8
2198  ELSE
2199  ord_in = hord
2200  END IF
2201  ord_ou = hord
2202  IF (.NOT.gridstruct%nested) CALL copy_corners_tlm(q, q_tl, npx, &
2203 & npy, 2, gridstruct%&
2204 & nested, bd, &
2205 & gridstruct%sw_corner&
2206 & , gridstruct%&
2207 & se_corner, gridstruct&
2208 & %nw_corner, &
2209 & gridstruct%ne_corner)
2210  fy2_tl = 0.0
2211  CALL yppm_tlm(fy2, fy2_tl, q, q_tl, cry, cry_tl, ord_in, isd, ied&
2212 & , isd, ied, js, je, jsd, jed, npx, npy, gridstruct%dya, &
2213 & gridstruct%nested, gridstruct%grid_type)
2214  fyy_tl = 0.0
2215  DO j=js,je+1
2216  DO i=isd,ied
2217  fyy_tl(i, j) = yfx_tl(i, j)*fy2(i, j) + yfx(i, j)*fy2_tl(i, j)
2218  fyy(i, j) = yfx(i, j)*fy2(i, j)
2219  END DO
2220  END DO
2221  q_i_tl = 0.0
2222  DO j=js,je
2223  DO i=isd,ied
2224  q_i_tl(i, j) = ((gridstruct%area(i, j)*q_tl(i, j)+fyy_tl(i, j)-&
2225 & fyy_tl(i, j+1))*ra_y(i, j)-(q(i, j)*gridstruct%area(i, j)+fyy(&
2226 & i, j)-fyy(i, j+1))*ra_y_tl(i, j))/ra_y(i, j)**2
2227  q_i(i, j) = (q(i, j)*gridstruct%area(i, j)+fyy(i, j)-fyy(i, j+1)&
2228 & )/ra_y(i, j)
2229  END DO
2230  END DO
2231  CALL xppm_tlm(fx, fx_tl, q_i, q_i_tl, crx(is:ie+1, js:je), crx_tl&
2232 & (is:ie+1, js:je), ord_ou, is, ie, isd, ied, js, je, jsd, &
2233 & jed, npx, npy, gridstruct%dxa, gridstruct%nested, &
2234 & gridstruct%grid_type)
2235  IF (.NOT.gridstruct%nested) CALL copy_corners_tlm(q, q_tl, npx, &
2236 & npy, 1, gridstruct%&
2237 & nested, bd, &
2238 & gridstruct%sw_corner&
2239 & , gridstruct%&
2240 & se_corner, gridstruct&
2241 & %nw_corner, &
2242 & gridstruct%ne_corner)
2243  fx2_tl = 0.0
2244  CALL xppm_tlm(fx2, fx2_tl, q, q_tl, crx, crx_tl, ord_in, is, ie, &
2245 & isd, ied, jsd, jed, jsd, jed, npx, npy, gridstruct%dxa, &
2246 & gridstruct%nested, gridstruct%grid_type)
2247  q_j_tl = 0.0
2248  fx1_tl = 0.0
2249  DO j=jsd,jed
2250  DO i=is,ie+1
2251  fx1_tl(i) = xfx_tl(i, j)*fx2(i, j) + xfx(i, j)*fx2_tl(i, j)
2252  fx1(i) = xfx(i, j)*fx2(i, j)
2253  END DO
2254  DO i=is,ie
2255  q_j_tl(i, j) = ((gridstruct%area(i, j)*q_tl(i, j)+fx1_tl(i)-&
2256 & fx1_tl(i+1))*ra_x(i, j)-(q(i, j)*gridstruct%area(i, j)+fx1(i)-&
2257 & fx1(i+1))*ra_x_tl(i, j))/ra_x(i, j)**2
2258  q_j(i, j) = (q(i, j)*gridstruct%area(i, j)+fx1(i)-fx1(i+1))/ra_x&
2259 & (i, j)
2260  END DO
2261  END DO
2262  CALL yppm_tlm(fy, fy_tl, q_j, q_j_tl, cry, cry_tl, ord_ou, is, ie&
2263 & , isd, ied, js, je, jsd, jed, npx, npy, gridstruct%dya, &
2264 & gridstruct%nested, gridstruct%grid_type)
2265 !----------------
2266 ! Flux averaging:
2267 !----------------
2268  IF (PRESENT(mfx) .AND. PRESENT(mfy)) THEN
2269 !---------------------------------
2270 ! For transport of pt and tracers
2271 !---------------------------------
2272  DO j=js,je
2273  DO i=is,ie+1
2274  fx_tl(i, j) = 0.5*((fx_tl(i, j)+fx2_tl(i, j))*mfx(i, j)+(fx(i&
2275 & , j)+fx2(i, j))*mfx_tl(i, j))
2276  fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*mfx(i, j)
2277  END DO
2278  END DO
2279  DO j=js,je+1
2280  DO i=is,ie
2281  fy_tl(i, j) = 0.5*((fy_tl(i, j)+fy2_tl(i, j))*mfy(i, j)+(fy(i&
2282 & , j)+fy2(i, j))*mfy_tl(i, j))
2283  fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*mfy(i, j)
2284  END DO
2285  END DO
2286  IF (PRESENT(nord) .AND. PRESENT(damp_c) .AND. PRESENT(mass)) THEN
2287  IF (damp_c .GT. 1.e-4) THEN
2288  pwx1 = damp_c*gridstruct%da_min
2289  pwy1 = nord + 1
2290  damp = pwx1**pwy1
2291  CALL deln_flux_tlm(nord, is, ie, js, je, npx, npy, damp, q&
2292 & , q_tl, fx, fx_tl, fy, fy_tl, gridstruct, bd, &
2293 & mass, mass_tl)
2294  END IF
2295  END IF
2296  ELSE
2297 !---------------------------------
2298 ! For transport of delp, vorticity
2299 !---------------------------------
2300  DO j=js,je
2301  DO i=is,ie+1
2302  fx_tl(i, j) = 0.5*((fx_tl(i, j)+fx2_tl(i, j))*xfx(i, j)+(fx(i&
2303 & , j)+fx2(i, j))*xfx_tl(i, j))
2304  fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*xfx(i, j)
2305  END DO
2306  END DO
2307  DO j=js,je+1
2308  DO i=is,ie
2309  fy_tl(i, j) = 0.5*((fy_tl(i, j)+fy2_tl(i, j))*yfx(i, j)+(fy(i&
2310 & , j)+fy2(i, j))*yfx_tl(i, j))
2311  fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*yfx(i, j)
2312  END DO
2313  END DO
2314  IF (PRESENT(nord) .AND. PRESENT(damp_c)) THEN
2315  IF (damp_c .GT. 1.e-4) THEN
2316  pwx1 = damp_c*gridstruct%da_min
2317  pwy1 = nord + 1
2318  damp = pwx1**pwy1
2319  CALL deln_flux_tlm(nord, is, ie, js, je, npx, npy, damp, q&
2320 & , q_tl, fx, fx_tl, fy, fy_tl, gridstruct, bd)
2321  END IF
2322  END IF
2323  END IF
2324  END SUBROUTINE fv_tp_2d_tlm
2325 ! Differentiation of xppm in forward (tangent) mode:
2326 ! variations of useful results: flux
2327 ! with respect to varying inputs: q flux c
2328  SUBROUTINE xppm_tlm(flux, flux_tl, q, q_tl, c, c_tl, iord, is, ie, &
2329 & isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
2330  IMPLICIT NONE
2331  INTEGER, INTENT(IN) :: is, ie, isd, ied, jsd, jed
2332 ! compute domain
2333  INTEGER, INTENT(IN) :: jfirst, jlast
2334  INTEGER, INTENT(IN) :: iord
2335  INTEGER, INTENT(IN) :: npx, npy
2336  REAL, INTENT(IN) :: q(isd:ied, jfirst:jlast)
2337  REAL, INTENT(IN) :: q_tl(isd:ied, jfirst:jlast)
2338 ! Courant N (like FLUX)
2339  REAL, INTENT(IN) :: c(is:ie+1, jfirst:jlast)
2340  REAL, INTENT(IN) :: c_tl(is:ie+1, jfirst:jlast)
2341  REAL, INTENT(IN) :: dxa(isd:ied, jsd:jed)
2342  LOGICAL, INTENT(IN) :: nested
2343  INTEGER, INTENT(IN) :: grid_type
2344 ! !OUTPUT PARAMETERS:
2345 ! Flux
2346  REAL, INTENT(OUT) :: flux(is:ie+1, jfirst:jlast)
2347  REAL, INTENT(OUT) :: flux_tl(is:ie+1, jfirst:jlast)
2348 ! Local
2349  REAL, DIMENSION(is-1:ie+1) :: bl, br, b0
2350  REAL :: q1(isd:ied)
2351  REAL :: q1_tl(isd:ied)
2352  REAL, DIMENSION(is:ie+1) :: fx0, fx1
2353  LOGICAL, DIMENSION(is-1:ie+1) :: smt5, smt6
2354  REAL :: al(is-1:ie+2)
2355  REAL :: al_tl(is-1:ie+2)
2356  REAL :: dm(is-2:ie+2)
2357  REAL :: dq(is-3:ie+2)
2358  INTEGER :: i, j, ie3, is1, ie1
2359  REAL :: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
2360  REAL :: xt_tl, qtmp_tl
2361  INTRINSIC max
2362  INTRINSIC min
2363  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
2364  IF (3 .LT. is - 1) THEN
2365  is1 = is - 1
2366  ELSE
2367  is1 = 3
2368  END IF
2369  IF (npx - 2 .GT. ie + 2) THEN
2370  ie3 = ie + 2
2371  ELSE
2372  ie3 = npx - 2
2373  END IF
2374  IF (npx - 3 .GT. ie + 1) THEN
2375  ie1 = ie + 1
2376  ELSE
2377  ie1 = npx - 3
2378  END IF
2379  al_tl = 0.0
2380  q1_tl = 0.0
2381  ELSE
2382  is1 = is - 1
2383  ie3 = ie + 2
2384  ie1 = ie + 1
2385  al_tl = 0.0
2386  q1_tl = 0.0
2387  END IF
2388  DO j=jfirst,jlast
2389  DO i=isd,ied
2390  q1_tl(i) = q_tl(i, j)
2391  q1(i) = q(i, j)
2392  END DO
2393  IF (iord .LT. 8 .OR. iord .EQ. 333) THEN
2394 ! ord = 2: perfectly linear ppm scheme
2395 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
2396  DO i=is1,ie3
2397  al_tl(i) = p1*(q1_tl(i-1)+q1_tl(i)) + p2*(q1_tl(i-2)+q1_tl(i+1&
2398 & ))
2399  al(i) = p1*(q1(i-1)+q1(i)) + p2*(q1(i-2)+q1(i+1))
2400  END DO
2401  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
2402  IF (is .EQ. 1) THEN
2403  al_tl(0) = c1*q1_tl(-2) + c2*q1_tl(-1) + c3*q1_tl(0)
2404  al(0) = c1*q1(-2) + c2*q1(-1) + c3*q1(0)
2405  al_tl(1) = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1_tl(0)-dxa(0, j&
2406 & )*q1_tl(-1))/(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, &
2407 & j))*q1_tl(1)-dxa(1, j)*q1_tl(2))/(dxa(1, j)+dxa(2, j)))
2408  al(1) = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-&
2409 & 1))/(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)&
2410 & -dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
2411  al_tl(2) = c3*q1_tl(1) + c2*q1_tl(2) + c1*q1_tl(3)
2412  al(2) = c3*q1(1) + c2*q1(2) + c1*q1(3)
2413  END IF
2414  IF (ie + 1 .EQ. npx) THEN
2415  al_tl(npx-1) = c1*q1_tl(npx-3) + c2*q1_tl(npx-2) + c3*q1_tl(&
2416 & npx-1)
2417  al(npx-1) = c1*q1(npx-3) + c2*q1(npx-2) + c3*q1(npx-1)
2418  al_tl(npx) = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1_tl(&
2419 & npx-1)-dxa(npx-1, j)*q1_tl(npx-2))/(dxa(npx-2, j)+dxa(npx-&
2420 & 1, j))+((2.*dxa(npx, j)+dxa(npx+1, j))*q1_tl(npx)-dxa(npx&
2421 & , j)*q1_tl(npx+1))/(dxa(npx, j)+dxa(npx+1, j)))
2422  al(npx) = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-&
2423 & dxa(npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((&
2424 & 2.*dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1&
2425 & ))/(dxa(npx, j)+dxa(npx+1, j)))
2426  al_tl(npx+1) = c3*q1_tl(npx) + c2*q1_tl(npx+1) + c1*q1_tl(&
2427 & npx+2)
2428  al(npx+1) = c3*q1(npx) + c2*q1(npx+1) + c1*q1(npx+2)
2429  END IF
2430  END IF
2431  IF (iord .EQ. 1) THEN
2432  DO i=is,ie+1
2433  IF (c(i, j) .GT. 0.) THEN
2434  flux_tl(i, j) = q1_tl(i-1)
2435  flux(i, j) = q1(i-1)
2436  ELSE
2437  flux_tl(i, j) = q1_tl(i)
2438  flux(i, j) = q1(i)
2439  END IF
2440  END DO
2441  ELSE IF (iord .EQ. 2) THEN
2442 ! perfectly linear scheme
2443 !DEC$ VECTOR ALWAYS
2444  DO i=is,ie+1
2445  xt_tl = c_tl(i, j)
2446  xt = c(i, j)
2447  IF (xt .GT. 0.) THEN
2448  qtmp_tl = q1_tl(i-1)
2449  qtmp = q1(i-1)
2450  flux_tl(i, j) = qtmp_tl + (1.-xt)*(al_tl(i)-qtmp_tl-xt_tl*&
2451 & (al(i-1)+al(i)-(qtmp+qtmp))-xt*(al_tl(i-1)+al_tl(i)-2*&
2452 & qtmp_tl)) - xt_tl*(al(i)-qtmp-xt*(al(i-1)+al(i)-(qtmp+&
2453 & qtmp)))
2454  flux(i, j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-&
2455 & (qtmp+qtmp)))
2456  ELSE
2457  qtmp_tl = q1_tl(i)
2458  qtmp = q1(i)
2459  flux_tl(i, j) = qtmp_tl + xt_tl*(al(i)-qtmp+xt*(al(i)+al(i&
2460 & +1)-(qtmp+qtmp))) + (1.+xt)*(al_tl(i)-qtmp_tl+xt_tl*(al(&
2461 & i)+al(i+1)-(qtmp+qtmp))+xt*(al_tl(i)+al_tl(i+1)-2*&
2462 & qtmp_tl))
2463  flux(i, j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-&
2464 & (qtmp+qtmp)))
2465  END IF
2466  END DO
2467  ELSE IF (iord .EQ. 333) THEN
2468 ! Perfectly linear scheme, more diffusive than ord=2 (HoldawayKent-2015-TellusA)
2469 !DEC$ VECTOR ALWAYS
2470  DO i=is,ie+1
2471  xt_tl = c_tl(i, j)
2472  xt = c(i, j)
2473  IF (xt .GT. 0.) THEN
2474  flux_tl(i, j) = (2.0*q1_tl(i)+5.0*q1_tl(i-1)-q1_tl(i-2))/&
2475 & 6.0 - 0.5*(xt_tl*(q1(i)-q1(i-1))+xt*(q1_tl(i)-q1_tl(i-1)&
2476 & )) + (xt_tl*xt+xt*xt_tl)*(q1(i)-2.0*q1(i-1)+q1(i-2))/6.0&
2477 & + xt**2*(q1_tl(i)-2.0*q1_tl(i-1)+q1_tl(i-2))/6.0
2478  flux(i, j) = (2.0*q1(i)+5.0*q1(i-1)-q1(i-2))/6.0 - 0.5*xt*&
2479 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i)-2.0*q1(i-1)+q1(i-2))
2480  ELSE
2481  flux_tl(i, j) = (2.0*q1_tl(i-1)+5.0*q1_tl(i)-q1_tl(i+1))/&
2482 & 6.0 - 0.5*(xt_tl*(q1(i)-q1(i-1))+xt*(q1_tl(i)-q1_tl(i-1)&
2483 & )) + (xt_tl*xt+xt*xt_tl)*(q1(i+1)-2.0*q1(i)+q1(i-1))/6.0&
2484 & + xt**2*(q1_tl(i+1)-2.0*q1_tl(i)+q1_tl(i-1))/6.0
2485  flux(i, j) = (2.0*q1(i-1)+5.0*q1(i)-q1(i+1))/6.0 - 0.5*xt*&
2486 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i+1)-2.0*q1(i)+q1(i-1))
2487  END IF
2488  END DO
2489  END IF
2490  END IF
2491  END DO
2492  END SUBROUTINE xppm_tlm
2493 ! Differentiation of yppm in forward (tangent) mode:
2494 ! variations of useful results: flux
2495 ! with respect to varying inputs: q flux c
2496  SUBROUTINE yppm_tlm(flux, flux_tl, q, q_tl, c, c_tl, jord, ifirst, &
2497 & ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
2498  IMPLICIT NONE
2499 ! Compute domain
2500  INTEGER, INTENT(IN) :: ifirst, ilast
2501  INTEGER, INTENT(IN) :: isd, ied, js, je, jsd, jed
2502  INTEGER, INTENT(IN) :: jord
2503  INTEGER, INTENT(IN) :: npx, npy
2504  REAL, INTENT(IN) :: q(ifirst:ilast, jsd:jed)
2505  REAL, INTENT(IN) :: q_tl(ifirst:ilast, jsd:jed)
2506 ! Courant number
2507  REAL, INTENT(IN) :: c(isd:ied, js:je+1)
2508  REAL, INTENT(IN) :: c_tl(isd:ied, js:je+1)
2509 ! Flux
2510  REAL, INTENT(OUT) :: flux(ifirst:ilast, js:je+1)
2511  REAL, INTENT(OUT) :: flux_tl(ifirst:ilast, js:je+1)
2512  REAL, INTENT(IN) :: dya(isd:ied, jsd:jed)
2513  LOGICAL, INTENT(IN) :: nested
2514  INTEGER, INTENT(IN) :: grid_type
2515 ! Local:
2516  REAL :: dm(ifirst:ilast, js-2:je+2)
2517  REAL :: al(ifirst:ilast, js-1:je+2)
2518  REAL :: al_tl(ifirst:ilast, js-1:je+2)
2519  REAL, DIMENSION(ifirst:ilast, js-1:je+1) :: bl, br, b0
2520  REAL :: dq(ifirst:ilast, js-3:je+2)
2521  REAL, DIMENSION(ifirst:ilast) :: fx0, fx1
2522  LOGICAL, DIMENSION(ifirst:ilast, js-1:je+1) :: smt5, smt6
2523  REAL :: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
2524  REAL :: xt_tl, qtmp_tl
2525  INTEGER :: i, j, js1, je3, je1
2526  INTRINSIC max
2527  INTRINSIC min
2528  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
2529  IF (3 .LT. js - 1) THEN
2530  js1 = js - 1
2531  ELSE
2532  js1 = 3
2533  END IF
2534  IF (npy - 2 .GT. je + 2) THEN
2535  je3 = je + 2
2536  ELSE
2537  je3 = npy - 2
2538  END IF
2539  IF (npy - 3 .GT. je + 1) THEN
2540  je1 = je + 1
2541  ELSE
2542  je1 = npy - 3
2543  END IF
2544  ELSE
2545 ! Nested grid OR Doubly periodic domain:
2546  js1 = js - 1
2547  je3 = je + 2
2548  je1 = je + 1
2549  END IF
2550  IF (jord .LT. 8 .OR. jord .EQ. 333) THEN
2551  al_tl = 0.0
2552  DO j=js1,je3
2553  DO i=ifirst,ilast
2554  al_tl(i, j) = p1*(q_tl(i, j-1)+q_tl(i, j)) + p2*(q_tl(i, j-2)+&
2555 & q_tl(i, j+1))
2556  al(i, j) = p1*(q(i, j-1)+q(i, j)) + p2*(q(i, j-2)+q(i, j+1))
2557  END DO
2558  END DO
2559  IF (.NOT.nested .AND. grid_type .LT. 3) THEN
2560  IF (js .EQ. 1) THEN
2561  DO i=ifirst,ilast
2562  al_tl(i, 0) = c1*q_tl(i, -2) + c2*q_tl(i, -1) + c3*q_tl(i, 0&
2563 & )
2564  al(i, 0) = c1*q(i, -2) + c2*q(i, -1) + c3*q(i, 0)
2565  al_tl(i, 1) = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q_tl(i, 0)-dya&
2566 & (i, 0)*q_tl(i, -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+&
2567 & dya(i, 2))*q_tl(i, 1)-dya(i, 1)*q_tl(i, 2))/(dya(i, 1)+dya&
2568 & (i, 2)))
2569  al(i, 1) = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)&
2570 & *q(i, -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2)&
2571 & )*q(i, 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
2572  al_tl(i, 2) = c3*q_tl(i, 1) + c2*q_tl(i, 2) + c1*q_tl(i, 3)
2573  al(i, 2) = c3*q(i, 1) + c2*q(i, 2) + c1*q(i, 3)
2574  END DO
2575  END IF
2576  IF (je + 1 .EQ. npy) THEN
2577  DO i=ifirst,ilast
2578  al_tl(i, npy-1) = c1*q_tl(i, npy-3) + c2*q_tl(i, npy-2) + c3&
2579 & *q_tl(i, npy-1)
2580  al(i, npy-1) = c1*q(i, npy-3) + c2*q(i, npy-2) + c3*q(i, npy&
2581 & -1)
2582  al_tl(i, npy) = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q_tl(&
2583 & i, npy-1)-dya(i, npy-1)*q_tl(i, npy-2))/(dya(i, npy-2)+dya&
2584 & (i, npy-1))+((2.*dya(i, npy)+dya(i, npy+1))*q_tl(i, npy)-&
2585 & dya(i, npy)*q_tl(i, npy+1))/(dya(i, npy)+dya(i, npy+1)))
2586  al(i, npy) = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy&
2587 & -1)-dya(i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1&
2588 & ))+((2.*dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q&
2589 & (i, npy+1))/(dya(i, npy)+dya(i, npy+1)))
2590  al_tl(i, npy+1) = c3*q_tl(i, npy) + c2*q_tl(i, npy+1) + c1*&
2591 & q_tl(i, npy+2)
2592  al(i, npy+1) = c3*q(i, npy) + c2*q(i, npy+1) + c1*q(i, npy+2&
2593 & )
2594  END DO
2595  END IF
2596  END IF
2597  IF (jord .EQ. 1) THEN
2598  DO j=js,je+1
2599  DO i=ifirst,ilast
2600  IF (c(i, j) .GT. 0.) THEN
2601  flux_tl(i, j) = q_tl(i, j-1)
2602  flux(i, j) = q(i, j-1)
2603  ELSE
2604  flux_tl(i, j) = q_tl(i, j)
2605  flux(i, j) = q(i, j)
2606  END IF
2607  END DO
2608  END DO
2609  ELSE IF (jord .EQ. 2) THEN
2610 ! Perfectly linear scheme
2611 ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7
2612  DO j=js,je+1
2613 !DEC$ VECTOR ALWAYS
2614  DO i=ifirst,ilast
2615  xt_tl = c_tl(i, j)
2616  xt = c(i, j)
2617  IF (xt .GT. 0.) THEN
2618  qtmp_tl = q_tl(i, j-1)
2619  qtmp = q(i, j-1)
2620  flux_tl(i, j) = qtmp_tl + (1.-xt)*(al_tl(i, j)-qtmp_tl-&
2621 & xt_tl*(al(i, j-1)+al(i, j)-(qtmp+qtmp))-xt*(al_tl(i, j-1&
2622 & )+al_tl(i, j)-2*qtmp_tl)) - xt_tl*(al(i, j)-qtmp-xt*(al(&
2623 & i, j-1)+al(i, j)-(qtmp+qtmp)))
2624  flux(i, j) = qtmp + (1.-xt)*(al(i, j)-qtmp-xt*(al(i, j-1)+&
2625 & al(i, j)-(qtmp+qtmp)))
2626  ELSE
2627  qtmp_tl = q_tl(i, j)
2628  qtmp = q(i, j)
2629  flux_tl(i, j) = qtmp_tl + xt_tl*(al(i, j)-qtmp+xt*(al(i, j&
2630 & )+al(i, j+1)-(qtmp+qtmp))) + (1.+xt)*(al_tl(i, j)-&
2631 & qtmp_tl+xt_tl*(al(i, j)+al(i, j+1)-(qtmp+qtmp))+xt*(&
2632 & al_tl(i, j)+al_tl(i, j+1)-2*qtmp_tl))
2633  flux(i, j) = qtmp + (1.+xt)*(al(i, j)-qtmp+xt*(al(i, j)+al&
2634 & (i, j+1)-(qtmp+qtmp)))
2635  END IF
2636  END DO
2637  END DO
2638  ELSE IF (jord .EQ. 333) THEN
2639 ! Perfectly linear scheme, more diffusive than ord=2 (HoldawayKent-2015-TellusA)
2640  DO j=js,je+1
2641 !DEC$ VECTOR ALWAYS
2642  DO i=ifirst,ilast
2643  xt_tl = c_tl(i, j)
2644  xt = c(i, j)
2645  IF (xt .GT. 0.) THEN
2646  flux_tl(i, j) = (2.0*q_tl(i, j)+5.0*q_tl(i, j-1)-q_tl(i, j&
2647 & -2))/6.0 - 0.5*(xt_tl*(q(i, j)-q(i, j-1))+xt*(q_tl(i, j)&
2648 & -q_tl(i, j-1))) + (xt_tl*xt+xt*xt_tl)*(q(i, j)-2.0*q(i, &
2649 & j-1)+q(i, j-2))/6.0 + xt**2*(q_tl(i, j)-2.0*q_tl(i, j-1)&
2650 & +q_tl(i, j-2))/6.0
2651  flux(i, j) = (2.0*q(i, j)+5.0*q(i, j-1)-q(i, j-2))/6.0 - &
2652 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j)-2.0*q(i&
2653 & , j-1)+q(i, j-2))
2654  ELSE
2655  flux_tl(i, j) = (2.0*q_tl(i, j-1)+5.0*q_tl(i, j)-q_tl(i, j&
2656 & +1))/6.0 - 0.5*(xt_tl*(q(i, j)-q(i, j-1))+xt*(q_tl(i, j)&
2657 & -q_tl(i, j-1))) + (xt_tl*xt+xt*xt_tl)*(q(i, j+1)-2.0*q(i&
2658 & , j)+q(i, j-1))/6.0 + xt**2*(q_tl(i, j+1)-2.0*q_tl(i, j)&
2659 & +q_tl(i, j-1))/6.0
2660  flux(i, j) = (2.0*q(i, j-1)+5.0*q(i, j)-q(i, j+1))/6.0 - &
2661 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j+1)-2.0*q(&
2662 & i, j)+q(i, j-1))
2663  END IF
2664  END DO
2665  END DO
2666  END IF
2667  RETURN
2668  END IF
2669  END SUBROUTINE yppm_tlm
2670 ! Differentiation of deln_flux in forward (tangent) mode:
2671 ! variations of useful results: fx fy
2672 ! with respect to varying inputs: q mass fx fy
2673  SUBROUTINE deln_flux_tlm(nord, is, ie, js, je, npx, npy, damp, q, &
2674 & q_tl, fx, fx_tl, fy, fy_tl, gridstruct, bd, mass, mass_tl)
2675  IMPLICIT NONE
2676 ! Del-n damping for the cell-mean values (A grid)
2677 !------------------
2678 ! nord = 0: del-2
2679 ! nord = 1: del-4
2680 ! nord = 2: del-6
2681 ! nord = 3: del-8 --> requires more ghosting than current
2682 !------------------
2683  TYPE(FV_GRID_BOUNDS_TYPE), INTENT(IN) :: bd
2684 ! del-n
2685  INTEGER, INTENT(IN) :: nord
2686  INTEGER, INTENT(IN) :: is, ie, js, je, npx, npy
2687  REAL, INTENT(IN) :: damp
2688 ! q ghosted on input
2689  REAL, INTENT(IN) :: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
2690  REAL, INTENT(IN) :: q_tl(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
2691  TYPE(FV_GRID_TYPE), INTENT(IN), TARGET :: gridstruct
2692 ! q ghosted on input
2693  REAL, OPTIONAL, INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
2694  REAL, OPTIONAL, INTENT(IN) :: mass_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
2695 ! diffusive fluxes:
2696  REAL, INTENT(INOUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je), fy(bd%is:bd%&
2697 & ie, bd%js:bd%je+1)
2698  REAL, INTENT(INOUT) :: fx_tl(bd%is:bd%ie+1, bd%js:bd%je), fy_tl(bd%&
2699 & is:bd%ie, bd%js:bd%je+1)
2700 ! local:
2701  REAL :: fx2(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2(bd%isd:bd%ied, bd%&
2702 & jsd:bd%jed+1)
2703  REAL :: fx2_tl(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2_tl(bd%isd:bd%ied&
2704 & , bd%jsd:bd%jed+1)
2705  REAL :: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
2706  REAL :: d2_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
2707  REAL :: damp2
2708  INTEGER :: i, j, n, nt, i1, i2, j1, j2
2709  INTRINSIC PRESENT
2710  i1 = is - 1 - nord
2711  i2 = ie + 1 + nord
2712  j1 = js - 1 - nord
2713  j2 = je + 1 + nord
2714  IF (.NOT.PRESENT(mass)) THEN
2715  d2_tl = 0.0
2716  DO j=j1,j2
2717  DO i=i1,i2
2718  d2_tl(i, j) = damp*q_tl(i, j)
2719  d2(i, j) = damp*q(i, j)
2720  END DO
2721  END DO
2722  ELSE
2723  d2_tl = 0.0
2724  DO j=j1,j2
2725  DO i=i1,i2
2726  d2_tl(i, j) = q_tl(i, j)
2727  d2(i, j) = q(i, j)
2728  END DO
2729  END DO
2730  END IF
2731  IF (nord .GT. 0) THEN
2732  CALL copy_corners_tlm(d2, d2_tl, npx, npy, 1, gridstruct%nested&
2733 & , bd, gridstruct%sw_corner, gridstruct%&
2734 & se_corner, gridstruct%nw_corner, gridstruct%&
2735 & ne_corner)
2736  fx2_tl = 0.0
2737  ELSE
2738  fx2_tl = 0.0
2739  END IF
2740  DO j=js-nord,je+nord
2741  DO i=is-nord,ie+nord+1
2742  fx2_tl(i, j) = gridstruct%del6_v(i, j)*(d2_tl(i-1, j)-d2_tl(i, j&
2743 & ))
2744  fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i-1, j)-d2(i, j))
2745  END DO
2746  END DO
2747  IF (nord .GT. 0) THEN
2748  CALL copy_corners_tlm(d2, d2_tl, npx, npy, 2, gridstruct%nested&
2749 & , bd, gridstruct%sw_corner, gridstruct%&
2750 & se_corner, gridstruct%nw_corner, gridstruct%&
2751 & ne_corner)
2752  fy2_tl = 0.0
2753  ELSE
2754  fy2_tl = 0.0
2755  END IF
2756  DO j=js-nord,je+nord+1
2757  DO i=is-nord,ie+nord
2758  fy2_tl(i, j) = gridstruct%del6_u(i, j)*(d2_tl(i, j-1)-d2_tl(i, j&
2759 & ))
2760  fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j-1)-d2(i, j))
2761  END DO
2762  END DO
2763  IF (nord .GT. 0) THEN
2764 !----------
2765 ! high-order
2766 !----------
2767  DO n=1,nord
2768  nt = nord - n
2769  DO j=js-nt-1,je+nt+1
2770  DO i=is-nt-1,ie+nt+1
2771  d2_tl(i, j) = gridstruct%rarea(i, j)*(fx2_tl(i, j)-fx2_tl(i+&
2772 & 1, j)+fy2_tl(i, j)-fy2_tl(i, j+1))
2773  d2(i, j) = (fx2(i, j)-fx2(i+1, j)+fy2(i, j)-fy2(i, j+1))*&
2774 & gridstruct%rarea(i, j)
2775  END DO
2776  END DO
2777  CALL copy_corners_tlm(d2, d2_tl, npx, npy, 1, gridstruct%&
2778 & nested, bd, gridstruct%sw_corner, gridstruct%&
2779 & se_corner, gridstruct%nw_corner, gridstruct%&
2780 & ne_corner)
2781  DO j=js-nt,je+nt
2782  DO i=is-nt,ie+nt+1
2783  fx2_tl(i, j) = gridstruct%del6_v(i, j)*(d2_tl(i, j)-d2_tl(i-&
2784 & 1, j))
2785  fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i, j)-d2(i-1, j))
2786  END DO
2787  END DO
2788  CALL copy_corners_tlm(d2, d2_tl, npx, npy, 2, gridstruct%&
2789 & nested, bd, gridstruct%sw_corner, gridstruct%&
2790 & se_corner, gridstruct%nw_corner, gridstruct%&
2791 & ne_corner)
2792  DO j=js-nt,je+nt+1
2793  DO i=is-nt,ie+nt
2794  fy2_tl(i, j) = gridstruct%del6_u(i, j)*(d2_tl(i, j)-d2_tl(i&
2795 & , j-1))
2796  fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j)-d2(i, j-1))
2797  END DO
2798  END DO
2799  END DO
2800  END IF
2801 !---------------------------------------------
2802 ! Add the diffusive fluxes to the flux arrays:
2803 !---------------------------------------------
2804  IF (PRESENT(mass)) THEN
2805 ! Apply mass weighting to diffusive fluxes:
2806  damp2 = 0.5*damp
2807  DO j=js,je
2808  DO i=is,ie+1
2809  fx_tl(i, j) = fx_tl(i, j) + damp2*((mass_tl(i-1, j)+mass_tl(i&
2810 & , j))*fx2(i, j)+(mass(i-1, j)+mass(i, j))*fx2_tl(i, j))
2811  fx(i, j) = fx(i, j) + damp2*(mass(i-1, j)+mass(i, j))*fx2(i, j&
2812 & )
2813  END DO
2814  END DO
2815  DO j=js,je+1
2816  DO i=is,ie
2817  fy_tl(i, j) = fy_tl(i, j) + damp2*((mass_tl(i, j-1)+mass_tl(i&
2818 & , j))*fy2(i, j)+(mass(i, j-1)+mass(i, j))*fy2_tl(i, j))
2819  fy(i, j) = fy(i, j) + damp2*(mass(i, j-1)+mass(i, j))*fy2(i, j&
2820 & )
2821  END DO
2822  END DO
2823  ELSE
2824  DO j=js,je
2825  DO i=is,ie+1
2826  fx_tl(i, j) = fx_tl(i, j) + fx2_tl(i, j)
2827  fx(i, j) = fx(i, j) + fx2(i, j)
2828  END DO
2829  END DO
2830  DO j=js,je+1
2831  DO i=is,ie
2832  fy_tl(i, j) = fy_tl(i, j) + fy2_tl(i, j)
2833  fy(i, j) = fy(i, j) + fy2(i, j)
2834  END DO
2835  END DO
2836  END IF
2837  END SUBROUTINE deln_flux_tlm
2838 ! Differentiation of copy_corners in forward (tangent) mode:
2839 ! variations of useful results: q
2840 ! with respect to varying inputs: q
2841 !Weird arguments are because this routine is called in a lot of
2842 !places outside of tp_core, sometimes very deeply nested in the call tree.
2843  SUBROUTINE copy_corners_tlm(q, q_tl, npx, npy, dir, nested, bd, &
2844 & sw_corner, se_corner, nw_corner, ne_corner)
2845  IMPLICIT NONE
2846  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
2847  INTEGER, INTENT(IN) :: npx, npy, dir
2848  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
2849  REAL, INTENT(INOUT) :: q_tl(bd%isd:bd%ied, bd%jsd:bd%jed)
2850  LOGICAL, INTENT(IN) :: nested, sw_corner, se_corner, nw_corner, &
2851 & ne_corner
2852  INTEGER :: i, j
2853  IF (nested) THEN
2854  RETURN
2855  ELSE IF (dir .EQ. 1) THEN
2856 ! XDir:
2857  IF (sw_corner) THEN
2858  DO j=1-ng,0
2859  DO i=1-ng,0
2860  q_tl(i, j) = q_tl(j, 1-i)
2861  q(i, j) = q(j, 1-i)
2862  END DO
2863  END DO
2864  END IF
2865  IF (se_corner) THEN
2866  DO j=1-ng,0
2867  DO i=npx,npx+ng-1
2868  q_tl(i, j) = q_tl(npy-j, i-npx+1)
2869  q(i, j) = q(npy-j, i-npx+1)
2870  END DO
2871  END DO
2872  END IF
2873  IF (ne_corner) THEN
2874  DO j=npy,npy+ng-1
2875  DO i=npx,npx+ng-1
2876  q_tl(i, j) = q_tl(j, 2*npx-1-i)
2877  q(i, j) = q(j, 2*npx-1-i)
2878  END DO
2879  END DO
2880  END IF
2881  IF (nw_corner) THEN
2882  DO j=npy,npy+ng-1
2883  DO i=1-ng,0
2884  q_tl(i, j) = q_tl(npy-j, i-1+npx)
2885  q(i, j) = q(npy-j, i-1+npx)
2886  END DO
2887  END DO
2888  END IF
2889  ELSE IF (dir .EQ. 2) THEN
2890 ! YDir:
2891  IF (sw_corner) THEN
2892  DO j=1-ng,0
2893  DO i=1-ng,0
2894  q_tl(i, j) = q_tl(1-j, i)
2895  q(i, j) = q(1-j, i)
2896  END DO
2897  END DO
2898  END IF
2899  IF (se_corner) THEN
2900  DO j=1-ng,0
2901  DO i=npx,npx+ng-1
2902  q_tl(i, j) = q_tl(npy+j-1, npx-i)
2903  q(i, j) = q(npy+j-1, npx-i)
2904  END DO
2905  END DO
2906  END IF
2907  IF (ne_corner) THEN
2908  DO j=npy,npy+ng-1
2909  DO i=npx,npx+ng-1
2910  q_tl(i, j) = q_tl(2*npy-1-j, i)
2911  q(i, j) = q(2*npy-1-j, i)
2912  END DO
2913  END DO
2914  END IF
2915  IF (nw_corner) THEN
2916  DO j=npy,npy+ng-1
2917  DO i=1-ng,0
2918  q_tl(i, j) = q_tl(j+1-npx, npy-i)
2919  q(i, j) = q(j+1-npx, npy-i)
2920  END DO
2921  END DO
2922  END IF
2923  END IF
2924  END SUBROUTINE copy_corners_tlm
2925 !Weird arguments are because this routine is called in a lot of
2926 !places outside of tp_core, sometimes very deeply nested in the call tree.
2927 end module tp_core_tlm_mod
tp_core_tlm_mod::t12
real, parameter t12
Definition: tp_core_tlm.F90:56
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
tp_core_tlm_mod::ppm_fac
real, parameter ppm_fac
Definition: tp_core_tlm.F90:35
tp_core_tlm_mod::xppm_tlm
subroutine xppm_tlm(flux, flux_tl, q, q_tl, c, c_tl, iord, is, ie, isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
Definition: tp_core_tlm.F90:2330
tp_core_tlm_mod::ppm_limiter
real, parameter ppm_limiter
Definition: tp_core_tlm.F90:38
tp_core_tlm_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_tlm.F90:239
tp_core_tlm_mod::p2
real, parameter p2
Definition: tp_core_tlm.F90:69
tp_core_tlm_mod::fv_tp_2d_tlm
subroutine, public fv_tp_2d_tlm(q, q_tl, crx, crx_tl, cry, cry_tl, npx, npy, hord, fx, fx_tl, fy, fy_tl, xfx, xfx_tl, yfx, yfx_tl, gridstruct, bd, ra_x, ra_x_tl, ra_y, ra_y_tl, mfx, mfx_tl, mfy, mfy_tl, mass, mass_tl, nord, damp_c)
Definition: tp_core_tlm.F90:2127
tp_core_tlm_mod::b4
real, parameter b4
Definition: tp_core_tlm.F90:53
fv_arrays_nlm_mod::fv_grid_type
Definition: fv_arrays_nlm.F90:115
tp_core_tlm_mod::yppm_tlm
subroutine yppm_tlm(flux, flux_tl, q, q_tl, c, c_tl, jord, ifirst, ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
Definition: tp_core_tlm.F90:2498
tp_core_tlm_mod::b2
real, parameter b2
Definition: tp_core_tlm.F90:51
tp_core_tlm_mod
Definition: tp_core_tlm.F90:20
tp_core_tlm_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_tlm.F90:959
fv_grid_utils_nlm_mod
Definition: fv_grid_utils_nlm.F90:20
tp_core_tlm_mod::pert_ppm
subroutine, public pert_ppm(im, a0, al, ar, iv)
Definition: tp_core_tlm.F90:1854
tp_core_tlm_mod::r3
real, parameter r3
Definition: tp_core_tlm.F90:36
fv_grid_utils_nlm_mod::big_number
real, parameter, public big_number
Definition: fv_grid_utils_nlm.F90:49
tp_core_tlm_mod::c1
real, parameter c1
Definition: tp_core_tlm.F90:62
fv_mp_nlm_mod::ng
integer, parameter, public ng
Definition: fv_mp_nlm_mod.F90:2465
tp_core_tlm_mod::c2
real, parameter c2
Definition: tp_core_tlm.F90:63
tp_core_tlm_mod::b3
real, parameter b3
Definition: tp_core_tlm.F90:52
fv_arrays_nlm_mod
Definition: fv_arrays_nlm.F90:20
tp_core_tlm_mod::near_zero
real, parameter near_zero
Definition: tp_core_tlm.F90:37
tp_core_tlm_mod::copy_corners
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core_tlm.F90:2048
fv_arrays_nlm_mod::r_grid
integer, parameter, public r_grid
Definition: fv_arrays_nlm.F90:32
tp_core_tlm_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_tlm.F90:1726
tp_core_tlm_mod::b5
real, parameter b5
Definition: tp_core_tlm.F90:54
tp_core_tlm_mod::t11
real, parameter t11
Definition: tp_core_tlm.F90:56
tp_core_tlm_mod::s11
real, parameter s11
Definition: tp_core_tlm.F90:57
tp_core_tlm_mod::c3
real, parameter c3
Definition: tp_core_tlm.F90:64
tp_core_tlm_mod::p1
real, parameter p1
Definition: tp_core_tlm.F90:68
tp_core_tlm_mod::s14
real, parameter s14
Definition: tp_core_tlm.F90:57
max
#define max(a, b)
Definition: mosaic_util.h:33
tp_core_tlm_mod::s15
real, parameter s15
Definition: tp_core_tlm.F90:57
tp_core_tlm_mod::deln_flux_tlm
subroutine deln_flux_tlm(nord, is, ie, js, je, npx, npy, damp, q, q_tl, fx, fx_tl, fy, fy_tl, gridstruct, bd, mass, mass_tl)
Definition: tp_core_tlm.F90:2675
tp_core_tlm_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_tlm.F90:85
tp_core_tlm_mod::b1
real, parameter b1
Definition: tp_core_tlm.F90:50
min
#define min(a, b)
Definition: mosaic_util.h:32
tp_core_tlm_mod::copy_corners_tlm
subroutine, public copy_corners_tlm(q, q_tl, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core_tlm.F90:2845
tp_core_tlm_mod::t13
real, parameter t13
Definition: tp_core_tlm.F90:56
tp_core_tlm_mod::deln_flux
subroutine deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy, gridstruct, bd, mass)
Definition: tp_core_tlm.F90:1920
fv_arrays_nlm_mod::fv_grid_bounds_type
Definition: fv_arrays_nlm.F90:601
unstructured_grid_mod::grid_type
Derived type containing the data.
Definition: unstructured_grid_mod.F90:23
tp_core_tlm_mod::pert_ppm_tlm
subroutine pert_ppm_tlm(im, a0, al, al_tl, ar, ar_tl, iv)
Definition: tp_core_tlm.F90:1778