40 real,
parameter::
r3 = 1./3.
47 real,
parameter::
b1 = 0.0375
48 real,
parameter::
b2 = -7./30.
49 real,
parameter::
b3 = -23./120.
50 real,
parameter::
b4 = 13./30.
51 real,
parameter::
b5 = -11./240.
54 real,
parameter::
b1 = 1./30.
55 real,
parameter::
b2 = -13./60.
56 real,
parameter::
b3 = -13./60.
57 real,
parameter::
b4 = 0.45
58 real,
parameter::
b5 = -0.05
60 real,
parameter::
t11 = 27./28.,
t12 = -13./28.,
t13=3./7.
61 real,
parameter::
s11 = 11./14.,
s14 = 4./7.,
s15=3./14.
66 real,
parameter::
c1 = -2./14.
67 real,
parameter::
c2 = 11./14.
68 real,
parameter::
c3 = 5./14.
72 real,
parameter::
p1 = 7./12.
73 real,
parameter::
p2 = -1./12.
109 SUBROUTINE fv_tp_2d_adm(q, q_ad, crx, crx_ad, cry, cry_ad, npx, npy, &
110 & hord, fx, fx_ad, fy, fy_ad, xfx, xfx_ad, yfx, yfx_ad, gridstruct, bd&
111 & , ra_x, ra_x_ad, ra_y, ra_y_ad, mfx, mfx_ad, mfy, mfy_ad, mass, &
112 & mass_ad, nord, damp_c)
115 INTEGER,
INTENT(IN) :: npx, npy
116 INTEGER,
INTENT(IN) :: hord
118 REAL,
INTENT(IN) :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed)
119 REAL :: crx_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
121 REAL,
INTENT(IN) :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed)
122 REAL :: xfx_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
124 REAL,
INTENT(IN) :: cry(bd%isd:bd%ied, bd%js:bd%je+1)
125 REAL :: cry_ad(bd%isd:bd%ied, bd%js:bd%je+1)
127 REAL,
INTENT(IN) :: yfx(bd%isd:bd%ied, bd%js:bd%je+1)
128 REAL :: yfx_ad(bd%isd:bd%ied, bd%js:bd%je+1)
129 REAL,
INTENT(IN) :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
130 REAL :: ra_x_ad(bd%is:bd%ie, bd%jsd:bd%jed)
131 REAL,
INTENT(IN) :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
132 REAL :: ra_y_ad(bd%isd:bd%ied, bd%js:bd%je)
134 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
135 REAL,
INTENT(INOUT) :: q_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
137 REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
138 REAL :: fx_ad(bd%is:bd%ie+1, bd%js:bd%je)
140 REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
141 REAL :: fy_ad(bd%is:bd%ie, bd%js:bd%je+1)
145 REAL,
OPTIONAL,
INTENT(IN) :: mfx(bd%is:bd%ie+1, bd%js:bd%je)
146 REAL,
OPTIONAL :: mfx_ad(bd%is:bd%ie+1, bd%js:bd%je)
148 REAL,
OPTIONAL,
INTENT(IN) :: mfy(bd%is:bd%ie, bd%js:bd%je+1)
149 REAL,
OPTIONAL :: mfy_ad(bd%is:bd%ie, bd%js:bd%je+1)
150 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
151 REAL,
OPTIONAL :: mass_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
152 REAL,
OPTIONAL,
INTENT(IN) :: damp_c
153 INTEGER,
OPTIONAL,
INTENT(IN) :: nord
155 INTEGER :: ord_ou, ord_in
156 REAL :: q_i(bd%isd:bd%ied, bd%js:bd%je)
157 REAL :: q_i_ad(bd%isd:bd%ied, bd%js:bd%je)
158 REAL :: q_j(bd%is:bd%ie, bd%jsd:bd%jed)
159 REAL :: q_j_ad(bd%is:bd%ie, bd%jsd:bd%jed)
160 REAL :: fx2(bd%is:bd%ie+1, bd%jsd:bd%jed)
161 REAL :: fx2_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
162 REAL :: fy2(bd%isd:bd%ied, bd%js:bd%je+1)
163 REAL :: fy2_ad(bd%isd:bd%ied, bd%js:bd%je+1)
164 REAL :: fyy(bd%isd:bd%ied, bd%js:bd%je+1)
165 REAL :: fyy_ad(bd%isd:bd%ied, bd%js:bd%je+1)
166 REAL :: fx1(bd%is:bd%ie+1, bd%jsd:bd%jed)
167 REAL :: fx1_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
170 INTEGER :: is, ie, js, je, isd, ied, jsd, jed
187 IF (hord .EQ. 10)
THEN 193 IF (.NOT.gridstruct%nested)
THEN 195 CALL copy_corners(q, npx, npy, 2, gridstruct%nested, bd, &
196 & gridstruct%sw_corner, gridstruct%se_corner, gridstruct&
197 & %nw_corner, gridstruct%ne_corner)
202 CALL yppm(fy2, q, cry, ord_in, isd, ied, isd, ied, js, je, jsd, jed&
203 & , npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%&
207 fyy(i, j) = yfx(i, j)*fy2(i, j)
212 q_i(i, j) = (q(i, j)*gridstruct%area(i, j)+fyy(i, j)-fyy(i, j+1)&
217 CALL xppm(fx, q_i, crx(is:ie+1, js:je), ord_ou, is, ie, isd, ied, js&
218 & , je, jsd, jed, npx, npy, gridstruct%dxa, gridstruct%nested, &
219 & gridstruct%grid_type)
220 IF (.NOT.gridstruct%nested)
THEN 222 CALL copy_corners(q, npx, npy, 1, gridstruct%nested, bd, &
223 & gridstruct%sw_corner, gridstruct%se_corner, gridstruct&
224 & %nw_corner, gridstruct%ne_corner)
229 CALL xppm(fx2, q, crx, ord_in, is, ie, isd, ied, jsd, jed, jsd, jed&
230 & , npx, npy, gridstruct%dxa, gridstruct%nested, gridstruct%&
234 fx1(i, j) = xfx(i, j)*fx2(i, j)
239 q_j(i, j) = (q(i, j)*gridstruct%area(i, j)+fx1(i, j)-fx1(i+1, j)&
244 CALL yppm(fy, q_j, cry, ord_ou, is, ie, isd, ied, js, je, jsd, jed, &
245 & npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%&
250 IF (
PRESENT(mfx) .AND.
PRESENT(mfy))
THEN 251 IF (
PRESENT(nord) .AND.
PRESENT(damp_c) .AND.
PRESENT(mass))
THEN 252 IF (damp_c .GT. 1.e-4)
THEN 253 damp = (damp_c*gridstruct%da_min)**(nord+1)
254 CALL deln_flux_adm(nord, is, ie, js, je, npx, npy, damp, q, &
255 & q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd, mass&
262 temp_ad4 = 0.5*mfy(i, j)*fy_ad(i, j)
263 fy2_ad(i, j) = fy2_ad(i, j) + temp_ad4
264 mfy_ad(i, j) = mfy_ad(i, j) + 0.5*(fy(i, j)+fy2(i, j))*fy_ad(i&
266 fy_ad(i, j) = temp_ad4
272 temp_ad3 = 0.5*mfx(i, j)*fx_ad(i, j)
273 fx2_ad(i, j) = fx2_ad(i, j) + temp_ad3
274 mfx_ad(i, j) = mfx_ad(i, j) + 0.5*(fx(i, j)+fx2(i, j))*fx_ad(i&
276 fx_ad(i, j) = temp_ad3
280 IF (
PRESENT(nord) .AND.
PRESENT(damp_c))
THEN 281 IF (damp_c .GT. 1.e-4)
THEN 282 damp = (damp_c*gridstruct%da_min)**(nord+1)
283 CALL deln_flux_adm(nord, is, ie, js, je, npx, npy, damp, q, &
284 & q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd)
290 temp_ad2 = 0.5*yfx(i, j)*fy_ad(i, j)
291 fy2_ad(i, j) = fy2_ad(i, j) + temp_ad2
292 yfx_ad(i, j) = yfx_ad(i, j) + 0.5*(fy(i, j)+fy2(i, j))*fy_ad(i&
294 fy_ad(i, j) = temp_ad2
300 temp_ad1 = 0.5*xfx(i, j)*fx_ad(i, j)
301 fx2_ad(i, j) = fx2_ad(i, j) + temp_ad1
302 xfx_ad(i, j) = xfx_ad(i, j) + 0.5*(fx(i, j)+fx2(i, j))*fx_ad(i&
304 fx_ad(i, j) = temp_ad1
310 CALL yppm_adm(fy, fy_ad, q_j, q_j_ad, cry, cry_ad, ord_ou, is, ie, &
311 & isd, ied, js, je, jsd, jed, npx, npy, gridstruct%dya, &
312 & gridstruct%nested, gridstruct%grid_type)
316 temp_ad0 = q_j_ad(i, j)/ra_x(i, j)
317 q_ad(i, j) = q_ad(i, j) + gridstruct%area(i, j)*temp_ad0
318 fx1_ad(i, j) = fx1_ad(i, j) + temp_ad0
319 fx1_ad(i+1, j) = fx1_ad(i+1, j) - temp_ad0
320 ra_x_ad(i, j) = ra_x_ad(i, j) - (gridstruct%area(i, j)*q(i, j)+&
321 & fx1(i, j)-fx1(i+1, j))*temp_ad0/ra_x(i, j)
327 xfx_ad(i, j) = xfx_ad(i, j) + fx2(i, j)*fx1_ad(i, j)
328 fx2_ad(i, j) = fx2_ad(i, j) + xfx(i, j)*fx1_ad(i, j)
332 CALL xppm_adm(fx2, fx2_ad, q, q_ad, crx, crx_ad, ord_in, is, ie, isd&
333 & , ied, jsd, jed, jsd, jed, npx, npy, gridstruct%dxa, &
334 & gridstruct%nested, gridstruct%grid_type)
336 IF (branch .EQ. 0)
THEN 339 & , gridstruct%sw_corner, gridstruct%se_corner, &
340 & gridstruct%nw_corner, gridstruct%ne_corner)
344 CALL xppm_adm(fx, fx_ad, q_i, q_i_ad, crx(is:ie+1, js:je), crx_ad(is&
345 & :ie+1, js:je), ord_ou, is, ie, isd, ied, js, je, jsd, jed, &
346 & npx, npy, gridstruct%dxa, gridstruct%nested, gridstruct%&
351 temp_ad = q_i_ad(i, j)/ra_y(i, j)
352 q_ad(i, j) = q_ad(i, j) + gridstruct%area(i, j)*temp_ad
353 fyy_ad(i, j) = fyy_ad(i, j) + temp_ad
354 fyy_ad(i, j+1) = fyy_ad(i, j+1) - temp_ad
355 ra_y_ad(i, j) = ra_y_ad(i, j) - (gridstruct%area(i, j)*q(i, j)+&
356 & fyy(i, j)-fyy(i, j+1))*temp_ad/ra_y(i, j)
362 yfx_ad(i, j) = yfx_ad(i, j) + fy2(i, j)*fyy_ad(i, j)
363 fy2_ad(i, j) = fy2_ad(i, j) + yfx(i, j)*fyy_ad(i, j)
367 CALL yppm_adm(fy2, fy2_ad, q, q_ad, cry, cry_ad, ord_in, isd, ied, &
368 & isd, ied, js, je, jsd, jed, npx, npy, gridstruct%dya, &
369 & gridstruct%nested, gridstruct%grid_type)
371 IF (branch .EQ. 0)
THEN 374 & , gridstruct%sw_corner, gridstruct%se_corner, &
375 & gridstruct%nw_corner, gridstruct%ne_corner)
383 SUBROUTINE fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, &
384 & gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
387 INTEGER,
INTENT(IN) :: npx, npy
388 INTEGER,
INTENT(IN) :: hord
390 REAL,
INTENT(IN) :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed)
392 REAL,
INTENT(IN) :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed)
394 REAL,
INTENT(IN) :: cry(bd%isd:bd%ied, bd%js:bd%je+1)
396 REAL,
INTENT(IN) :: yfx(bd%isd:bd%ied, bd%js:bd%je+1)
397 REAL,
INTENT(IN) :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
398 REAL,
INTENT(IN) :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
400 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
402 REAL,
INTENT(OUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je)
404 REAL,
INTENT(OUT) :: fy(bd%is:bd%ie, bd%js:bd%je+1)
408 REAL,
OPTIONAL,
INTENT(IN) :: mfx(bd%is:bd%ie+1, bd%js:bd%je)
410 REAL,
OPTIONAL,
INTENT(IN) :: mfy(bd%is:bd%ie, bd%js:bd%je+1)
411 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
412 REAL,
OPTIONAL,
INTENT(IN) :: damp_c
413 INTEGER,
OPTIONAL,
INTENT(IN) :: nord
415 INTEGER :: ord_ou, ord_in
416 REAL :: q_i(bd%isd:bd%ied, bd%js:bd%je)
417 REAL :: q_j(bd%is:bd%ie, bd%jsd:bd%jed)
418 REAL :: fx2(bd%is:bd%ie+1, bd%jsd:bd%jed)
419 REAL :: fy2(bd%isd:bd%ied, bd%js:bd%je+1)
420 REAL :: fyy(bd%isd:bd%ied, bd%js:bd%je+1)
421 REAL :: fx1(bd%is:bd%ie+1, bd%jsd:bd%jed)
424 INTEGER :: is, ie, js, je, isd, ied, jsd, jed
434 IF (hord .EQ. 10)
THEN 440 IF (.NOT.gridstruct%nested)
CALL copy_corners(q, npx, npy, 2, &
441 & gridstruct%nested, bd, &
442 & gridstruct%sw_corner, &
443 & gridstruct%se_corner, &
444 & gridstruct%nw_corner, &
445 & gridstruct%ne_corner)
446 CALL yppm(fy2, q, cry, ord_in, isd, ied, isd, ied, js, je, jsd, jed&
447 & , npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%&
451 fyy(i, j) = yfx(i, j)*fy2(i, j)
456 q_i(i, j) = (q(i, j)*gridstruct%area(i, j)+fyy(i, j)-fyy(i, j+1)&
460 CALL xppm(fx, q_i, crx(is:ie+1, js:je), ord_ou, is, ie, isd, ied, js&
461 & , je, jsd, jed, npx, npy, gridstruct%dxa, gridstruct%nested, &
462 & gridstruct%grid_type)
463 IF (.NOT.gridstruct%nested)
CALL copy_corners(q, npx, npy, 1, &
464 & gridstruct%nested, bd, &
465 & gridstruct%sw_corner, &
466 & gridstruct%se_corner, &
467 & gridstruct%nw_corner, &
468 & gridstruct%ne_corner)
469 CALL xppm(fx2, q, crx, ord_in, is, ie, isd, ied, jsd, jed, jsd, jed&
470 & , npx, npy, gridstruct%dxa, gridstruct%nested, gridstruct%&
474 fx1(i, j) = xfx(i, j)*fx2(i, j)
479 q_j(i, j) = (q(i, j)*gridstruct%area(i, j)+fx1(i, j)-fx1(i+1, j)&
483 CALL yppm(fy, q_j, cry, ord_ou, is, ie, isd, ied, js, je, jsd, jed, &
484 & npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%&
489 IF (
PRESENT(mfx) .AND.
PRESENT(mfy))
THEN 495 fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*mfx(i, j)
500 fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*mfy(i, j)
503 IF (
PRESENT(nord) .AND.
PRESENT(damp_c) .AND.
PRESENT(mass))
THEN 504 IF (damp_c .GT. 1.e-4)
THEN 505 damp = (damp_c*gridstruct%da_min)**(nord+1)
506 CALL deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy&
507 & , gridstruct, bd, mass)
516 fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*xfx(i, j)
521 fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*yfx(i, j)
524 IF (
PRESENT(nord) .AND.
PRESENT(damp_c))
THEN 525 IF (damp_c .GT. 1.e-4)
THEN 526 damp = (damp_c*gridstruct%da_min)**(nord+1)
527 CALL deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy&
535 SUBROUTINE copy_corners(q, npx, npy, dir, nested, bd, sw_corner, &
536 & se_corner, nw_corner, ne_corner)
539 INTEGER,
INTENT(IN) :: npx, npy, dir
540 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
541 LOGICAL,
INTENT(IN) :: nested, sw_corner, se_corner, nw_corner, &
546 ELSE IF (dir .EQ. 1)
THEN 558 q(i, j) = q(npy-j, i-npx+1)
565 q(i, j) = q(j, 2*npx-1-i)
572 q(i, j) = q(npy-j, i-1+npx)
576 ELSE IF (dir .EQ. 2)
THEN 588 q(i, j) = q(npy+j-1, npx-i)
595 q(i, j) = q(2*npy-1-j, i)
602 q(i, j) = q(j+1-npx, npy-i)
628 SUBROUTINE xppm_adm(flux, flux_ad, q, q_ad, c, c_ad, iord, is, ie, isd&
629 & , ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
631 INTEGER,
INTENT(IN) :: is, ie, isd, ied, jsd, jed
633 INTEGER,
INTENT(IN) :: jfirst, jlast
634 INTEGER,
INTENT(IN) :: iord
635 INTEGER,
INTENT(IN) :: npx, npy
636 REAL,
INTENT(IN) :: q(isd:ied, jfirst:jlast)
637 REAL :: q_ad(isd:ied, jfirst:jlast)
639 REAL,
INTENT(IN) :: c(is:ie+1, jfirst:jlast)
640 REAL :: c_ad(is:ie+1, jfirst:jlast)
641 REAL,
INTENT(IN) :: dxa(isd:ied, jsd:jed)
642 LOGICAL,
INTENT(IN) :: nested
643 INTEGER,
INTENT(IN) :: grid_type
646 REAL :: flux(is:ie+1, jfirst:jlast)
647 REAL :: flux_ad(is:ie+1, jfirst:jlast)
649 REAL,
DIMENSION(is-1:ie+1) :: bl, br, b0
650 REAL,
DIMENSION(is-1:ie+1) :: bl_ad, br_ad, b0_ad
652 REAL :: q1_ad(isd:ied)
653 REAL,
DIMENSION(is:ie+1) :: fx0, fx1
654 REAL,
DIMENSION(is:ie+1) :: fx0_ad, fx1_ad
655 LOGICAL,
DIMENSION(is-1:ie+1) :: smt5, smt6
656 REAL :: al(is-1:ie+2)
657 REAL :: al_ad(is-1:ie+2)
658 REAL :: dm(is-2:ie+2)
659 REAL :: dm_ad(is-2:ie+2)
660 REAL :: dq(is-3:ie+2)
661 REAL :: dq_ad(is-3:ie+2)
662 INTEGER :: i, j, ie3, is1, ie1
663 REAL :: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
664 REAL :: xt_ad, qtmp_ad, pmp_1_ad, lac_1_ad, pmp_2_ad, lac_2_ad
785 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 786 IF (3 .LT. is - 1)
THEN 791 IF (npx - 2 .GT. ie + 2)
THEN 796 IF (npx - 3 .GT. ie + 1)
THEN 814 IF (iord .LT. 8 .OR. iord .EQ. 333)
THEN 819 al(i) =
p1*(q1(i-1)+q1(i)) +
p2*(q1(i-2)+q1(i+1))
821 IF (iord .EQ. 7)
THEN 823 IF (al(i) .LT. 0.)
THEN 825 al(i) = 0.5*(q1(i-1)+q1(i))
835 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 838 al(0) =
c1*q1(-2) +
c2*q1(-1) +
c3*q1(0)
840 al(1) = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-&
841 & 1))/(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)&
842 & -dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
844 al(2) =
c3*q1(1) +
c2*q1(2) +
c1*q1(3)
845 IF (iord .EQ. 7)
THEN 846 IF (0. .LT. al(0))
THEN 855 IF (0. .LT. al(1))
THEN 864 IF (0. .LT. al(2))
THEN 879 IF (ie + 1 .EQ. npx)
THEN 881 al(npx-1) =
c1*q1(npx-3) +
c2*q1(npx-2) +
c3*q1(npx-1)
883 al(npx) = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-&
884 & dxa(npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((&
885 & 2.*dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1&
886 & ))/(dxa(npx, j)+dxa(npx+1, j)))
888 al(npx+1) =
c3*q1(npx) +
c2*q1(npx+1) +
c1*q1(npx+2)
889 IF (iord .EQ. 7)
THEN 890 IF (0. .LT. al(npx-1))
THEN 892 al(npx-1) = al(npx-1)
899 IF (0. .LT. al(npx))
THEN 908 IF (0. .LT. al(npx+1))
THEN 910 al(npx+1) = al(npx+1)
926 IF (iord .EQ. 1)
THEN 928 IF (c(i, j) .GT. 0.)
THEN 935 ELSE IF (iord .EQ. 2)
THEN 951 ELSE IF (iord .EQ. 333)
THEN 968 ELSE IF (iord .EQ. 3)
THEN 971 bl(i) = al(i) - q1(i)
973 br(i) = al(i+1) - q1(i)
975 b0(i) = bl(i) + br(i)
976 IF (b0(i) .GE. 0.)
THEN 981 IF (bl(i) - br(i) .GE. 0.)
THEN 991 smt6(i) = 3.*x0 .LT. xt
1000 IF (xt .GT. 0.)
THEN 1001 IF (smt6(i-1) .OR. smt5(i))
THEN 1003 fx1(i) = br(i-1) - xt*b0(i-1)
1005 ELSE IF (smt5(i-1))
THEN 1006 IF (bl(i-1) .GE. 0.)
THEN 1013 IF (br(i-1) .GE. 0.)
THEN 1020 IF (x2 .GT. y1)
THEN 1031 fx1(i) = sign(min1, br(i-1))
1036 ELSE IF (smt6(i) .OR. smt5(i-1))
THEN 1038 fx1(i) = bl(i) + xt*b0(i)
1040 ELSE IF (smt5(i))
THEN 1041 IF (bl(i) .GE. 0.)
THEN 1048 IF (br(i) .GE. 0.)
THEN 1055 IF (x3 .GT. y2)
THEN 1065 fx1(i) = sign(min2, bl(i))
1070 IF (xt .GE. 0.)
THEN 1081 ELSE IF (iord .EQ. 4)
THEN 1084 bl(i) = al(i) - q1(i)
1086 br(i) = al(i+1) - q1(i)
1088 b0(i) = bl(i) + br(i)
1089 IF (b0(i) .GE. 0.)
THEN 1094 IF (bl(i) - br(i) .GE. 0.)
THEN 1103 smt5(i) = x0 .LT. xt
1104 smt6(i) = 3.*x0 .LT. xt
1112 IF (c(i, j) .GT. 0.)
THEN 1113 IF (smt6(i-1) .OR. smt5(i))
THEN 1115 fx1(i) = (1.-c(i, j))*(br(i-1)-c(i, j)*b0(i-1))
1120 ELSE IF (smt6(i) .OR. smt5(i-1))
THEN 1122 fx1(i) = (1.+c(i, j))*(bl(i)+c(i, j)*b0(i))
1131 IF (iord .EQ. 5)
THEN 1134 bl(i) = al(i) - q1(i)
1136 br(i) = al(i+1) - q1(i)
1138 b0(i) = bl(i) + br(i)
1139 smt5(i) = bl(i)*br(i) .LT. 0.
1145 bl(i) = al(i) - q1(i)
1147 br(i) = al(i+1) - q1(i)
1149 b0(i) = bl(i) + br(i)
1150 IF (3.*b0(i) .GE. 0.)
THEN 1155 IF (bl(i) - br(i) .GE. 0.)
THEN 1156 abs4 = bl(i) - br(i)
1158 abs4 = -(bl(i)-br(i))
1160 smt5(i) = abs1 .LT. abs4
1166 IF (c(i, j) .GT. 0.)
THEN 1168 fx1(i) = (1.-c(i, j))*(br(i-1)-c(i, j)*b0(i-1))
1172 fx1(i) = (1.+c(i, j))*(bl(i)+c(i, j)*b0(i))
1175 IF (smt5(i-1) .OR. smt5(i))
THEN 1190 xt = 0.25*(q1(i+1)-q1(i-1))
1191 IF (xt .GE. 0.)
THEN 1198 IF (q1(i-1) .LT. q1(i))
THEN 1199 IF (q1(i) .LT. q1(i+1))
THEN 1206 ELSE IF (q1(i-1) .LT. q1(i+1))
THEN 1214 IF (q1(i-1) .GT. q1(i))
THEN 1215 IF (q1(i) .GT. q1(i+1))
THEN 1222 ELSE IF (q1(i-1) .GT. q1(i+1))
THEN 1230 IF (x4 .GT. y3)
THEN 1231 IF (y3 .GT. z1)
THEN 1240 ELSE IF (x4 .GT. z1)
THEN 1249 dm(i) = sign(min3, xt)
1253 al(i) = 0.5*(q1(i-1)+q1(i)) +
r3*(dm(i-1)-dm(i))
1255 IF (iord .EQ. 8)
THEN 1259 IF (xt .GE. 0.)
THEN 1266 IF (al(i) - q1(i) .GE. 0.)
THEN 1273 IF (x5 .GT. y4)
THEN 1283 bl(i) = -sign(min4, xt)
1284 IF (xt .GE. 0.)
THEN 1291 IF (al(i+1) - q1(i) .GE. 0.)
THEN 1292 y5 = al(i+1) - q1(i)
1295 y5 = -(al(i+1)-q1(i))
1298 IF (x6 .GT. y5)
THEN 1308 br(i) = sign(min5, xt)
1311 ELSE IF (iord .EQ. 11)
THEN 1316 IF (xt .GE. 0.)
THEN 1323 IF (al(i) - q1(i) .GE. 0.)
THEN 1330 IF (x7 .GT. y6)
THEN 1340 bl(i) = -sign(min6, xt)
1341 IF (xt .GE. 0.)
THEN 1348 IF (al(i+1) - q1(i) .GE. 0.)
THEN 1349 y7 = al(i+1) - q1(i)
1352 y7 = -(al(i+1)-q1(i))
1355 IF (x8 .GT. y7)
THEN 1365 br(i) = sign(min7, xt)
1370 dq(i) = 2.*(q1(i+1)-q1(i))
1374 bl(i) = al(i) - q1(i)
1376 br(i) = al(i+1) - q1(i)
1377 IF (dm(i-1) .GE. 0.)
THEN 1382 IF (dm(i) .GE. 0.)
THEN 1387 IF (dm(i+1) .GE. 0.)
THEN 1392 IF (abs2 + abs5 + abs7 .LT.
near_zero)
THEN 1397 IF (3.*(bl(i)+br(i)) .GE. 0.)
THEN 1398 abs3 = 3.*(bl(i)+br(i))
1400 abs3 = -(3.*(bl(i)+br(i)))
1402 IF (bl(i) - br(i) .GE. 0.)
THEN 1403 abs6 = bl(i) - br(i)
1405 abs6 = -(bl(i)-br(i))
1407 IF (abs3 .GT. abs6)
THEN 1409 lac_2 = pmp_2 - 0.75*dq(i-2)
1410 IF (0. .LT. pmp_2)
THEN 1411 IF (pmp_2 .LT. lac_2)
THEN 1418 ELSE IF (0. .LT. lac_2)
THEN 1425 IF (0. .GT. pmp_2)
THEN 1426 IF (pmp_2 .GT. lac_2)
THEN 1433 ELSE IF (0. .GT. lac_2)
THEN 1440 IF (br(i) .LT. y14)
THEN 1447 IF (x9 .GT. y8)
THEN 1455 lac_1 = pmp_1 + 0.75*dq(i+1)
1456 IF (0. .LT. pmp_1)
THEN 1457 IF (pmp_1 .LT. lac_1)
THEN 1464 ELSE IF (0. .LT. lac_1)
THEN 1471 IF (0. .GT. pmp_1)
THEN 1472 IF (pmp_1 .GT. lac_1)
THEN 1479 ELSE IF (0. .GT. lac_1)
THEN 1486 IF (bl(i) .LT. y15)
THEN 1493 IF (x10 .GT. y9)
THEN 1508 IF (iord .EQ. 9 .OR. iord .EQ. 13)
THEN 1509 arg1 = ie1 - is1 + 1
1512 CALL pert_ppm(arg1, q1(is1:ie1), bl(is1:ie1), br(is1:ie1), 0)
1517 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 1520 bl(0) =
s14*dm(-1) +
s11*(q1(-1)-q1(0))
1522 xt = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-1))&
1523 & /(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)-&
1524 & dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
1525 IF (q1(1) .GT. q1(2))
THEN 1532 IF (q1(-1) .GT. q1(0))
THEN 1533 IF (q1(0) .GT. z2)
THEN 1540 ELSE IF (q1(-1) .GT. z2)
THEN 1547 IF (xt .LT. y10)
THEN 1554 IF (q1(1) .LT. q1(2))
THEN 1561 IF (q1(-1) .LT. q1(0))
THEN 1562 IF (q1(0) .LT. z3)
THEN 1569 ELSE IF (q1(-1) .LT. z3)
THEN 1576 IF (xt .GT. y11)
THEN 1594 br(2) = al(3) - q1(2)
1597 CALL pert_ppm(3, q1(0:2), bl(0:2), br(0:2), 1)
1602 IF (ie + 1 .EQ. npx)
THEN 1604 bl(npx-2) = al(npx-2) - q1(npx-2)
1606 xt =
s15*q1(npx-1) +
s11*q1(npx-2) +
s14*dm(npx-2)
1608 br(npx-2) = xt - q1(npx-2)
1610 bl(npx-1) = xt - q1(npx-1)
1611 xt = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-dxa(&
1612 & npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((2.*&
1613 & dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1))/&
1614 & (dxa(npx, j)+dxa(npx+1, j)))
1615 IF (q1(npx) .GT. q1(npx+1))
THEN 1622 IF (q1(npx-2) .GT. q1(npx-1))
THEN 1623 IF (q1(npx-1) .GT. z4)
THEN 1630 ELSE IF (q1(npx-2) .GT. z4)
THEN 1637 IF (xt .LT. y12)
THEN 1644 IF (q1(npx) .LT. q1(npx+1))
THEN 1651 IF (q1(npx-2) .LT. q1(npx-1))
THEN 1652 IF (q1(npx-1) .LT. z5)
THEN 1659 ELSE IF (q1(npx-2) .LT. z5)
THEN 1666 IF (xt .GT. y13)
THEN 1675 br(npx-1) = xt - q1(npx-1)
1677 bl(npx) = xt - q1(npx)
1679 br(npx) =
s11*(q1(npx+1)-q1(npx)) -
s14*dm(npx+1)
1682 CALL pert_ppm(3, q1(npx-2:npx), bl(npx-2:npx), br(npx-2:npx)&
1692 IF (c(i, j) .GT. 0.)
THEN 1710 DO j=jlast,jfirst,-1
1712 IF (branch .LT. 3)
THEN 1713 IF (branch .EQ. 0)
THEN 1716 IF (branch .EQ. 0)
THEN 1717 q1_ad(i) = q1_ad(i) + flux_ad(i, j)
1720 q1_ad(i-1) = q1_ad(i-1) + flux_ad(i, j)
1724 ELSE IF (branch .EQ. 1)
THEN 1727 IF (branch .EQ. 0)
THEN 1730 temp0 = al(i) + al(i+1) - 2*qtmp
1731 temp_ad5 = (xt+1.)*flux_ad(i, j)
1732 temp_ad6 = xt*temp_ad5
1733 qtmp_ad = flux_ad(i, j) - temp_ad5 - 2*temp_ad6
1734 xt_ad = temp0*temp_ad5 + (al(i)-qtmp+xt*temp0)*flux_ad(i, &
1736 al_ad(i) = al_ad(i) + temp_ad6 + temp_ad5
1737 al_ad(i+1) = al_ad(i+1) + temp_ad6
1740 q1_ad(i) = q1_ad(i) + qtmp_ad
1744 temp = al(i-1) + al(i) - 2*qtmp
1745 temp_ad3 = (1.-xt)*flux_ad(i, j)
1746 temp_ad4 = -(xt*temp_ad3)
1747 qtmp_ad = flux_ad(i, j) - temp_ad3 - 2*temp_ad4
1748 xt_ad = -(temp*temp_ad3) - (al(i)-qtmp-xt*temp)*flux_ad(i&
1750 al_ad(i) = al_ad(i) + temp_ad4 + temp_ad3
1751 al_ad(i-1) = al_ad(i-1) + temp_ad4
1754 q1_ad(i-1) = q1_ad(i-1) + qtmp_ad
1757 c_ad(i, j) = c_ad(i, j) + xt_ad
1762 IF (branch .EQ. 0)
THEN 1764 temp_ad10 = flux_ad(i, j)/6.0
1765 temp_ad11 = -(0.5*xt*flux_ad(i, j))
1766 temp_ad12 = xt**2*flux_ad(i, j)/6.0
1767 q1_ad(i-1) = q1_ad(i-1) + temp_ad12 - temp_ad11 + 2.0*&
1769 q1_ad(i) = q1_ad(i) + temp_ad11 - 2.0*temp_ad12 + 5.0*&
1771 q1_ad(i+1) = q1_ad(i+1) + temp_ad12 - temp_ad10
1772 xt_ad = ((q1(i+1)-2.0*q1(i)+q1(i-1))*2*xt/6.0-0.5*(q1(i)-&
1773 & q1(i-1)))*flux_ad(i, j)
1777 temp_ad7 = flux_ad(i, j)/6.0
1778 temp_ad8 = -(0.5*xt*flux_ad(i, j))
1779 temp_ad9 = xt**2*flux_ad(i, j)/6.0
1780 q1_ad(i) = q1_ad(i) + temp_ad9 + temp_ad8 + 2.0*temp_ad7
1781 q1_ad(i-1) = q1_ad(i-1) + 5.0*temp_ad7 - temp_ad8 - 2.0*&
1783 q1_ad(i-2) = q1_ad(i-2) + temp_ad9 - temp_ad7
1784 xt_ad = ((q1(i)-2.0*q1(i-1)+q1(i-2))*2*xt/6.0-0.5*(q1(i)-&
1785 & q1(i-1)))*flux_ad(i, j)
1789 c_ad(i, j) = c_ad(i, j) + xt_ad
1792 ELSE IF (branch .LT. 5)
THEN 1793 IF (branch .EQ. 3)
THEN 1795 fx0_ad(i) = fx0_ad(i) + flux_ad(i, j)
1796 abs0_ad = -(fx1(i)*flux_ad(i, j))
1797 fx1_ad(i) = fx1_ad(i) + (1.-abs0)*flux_ad(i, j)
1801 IF (branch .EQ. 0)
THEN 1809 IF (branch .LT. 3)
THEN 1810 IF (branch .NE. 0)
THEN 1811 IF (branch .EQ. 1)
THEN 1813 min2_ad = sign(1.d0, min2*bl(i))*fx1_ad(i)
1816 IF (branch .EQ. 0)
THEN 1826 IF (branch .EQ. 0)
THEN 1827 br_ad(i) = br_ad(i) + y2_ad
1829 br_ad(i) = br_ad(i) - y2_ad
1832 IF (branch .EQ. 0)
THEN 1833 bl_ad(i) = bl_ad(i) + x3_ad
1835 bl_ad(i) = bl_ad(i) - x3_ad
1839 bl_ad(i) = bl_ad(i) + fx1_ad(i)
1840 xt_ad = xt_ad + b0(i)*fx1_ad(i)
1841 b0_ad(i) = b0_ad(i) + xt*fx1_ad(i)
1845 q1_ad(i) = q1_ad(i) + fx0_ad(i)
1848 IF (branch .NE. 3)
THEN 1849 IF (branch .EQ. 4)
THEN 1851 min1_ad = sign(1.d0, min1*br(i-1))*fx1_ad(i)
1854 IF (branch .EQ. 0)
THEN 1864 IF (branch .EQ. 0)
THEN 1865 br_ad(i-1) = br_ad(i-1) + y1_ad
1867 br_ad(i-1) = br_ad(i-1) - y1_ad
1870 IF (branch .EQ. 0)
THEN 1871 bl_ad(i-1) = bl_ad(i-1) + x2_ad
1873 bl_ad(i-1) = bl_ad(i-1) - x2_ad
1877 br_ad(i-1) = br_ad(i-1) + fx1_ad(i)
1878 xt_ad = xt_ad - b0(i-1)*fx1_ad(i)
1879 b0_ad(i-1) = b0_ad(i-1) - xt*fx1_ad(i)
1883 q1_ad(i-1) = q1_ad(i-1) + fx0_ad(i)
1887 c_ad(i, j) = c_ad(i, j) + xt_ad
1895 IF (branch .EQ. 0)
THEN 1901 bl_ad(i) = bl_ad(i) + b0_ad(i)
1902 br_ad(i) = br_ad(i) + b0_ad(i)
1905 al_ad(i+1) = al_ad(i+1) + br_ad(i)
1906 q1_ad(i) = q1_ad(i) - bl_ad(i) - br_ad(i)
1909 al_ad(i) = al_ad(i) + bl_ad(i)
1914 fx0_ad(i) = fx0_ad(i) + flux_ad(i, j)
1915 fx1_ad(i) = fx1_ad(i) + flux_ad(i, j)
1918 IF (branch .LT. 2)
THEN 1919 IF (branch .EQ. 0)
THEN 1921 temp_ad13 = (1.-c(i, j))*fx1_ad(i)
1922 c_ad(i, j) = c_ad(i, j) - b0(i-1)*temp_ad13 - (br(i-1)-c&
1923 & (i, j)*b0(i-1))*fx1_ad(i)
1924 br_ad(i-1) = br_ad(i-1) + temp_ad13
1925 b0_ad(i-1) = b0_ad(i-1) - c(i, j)*temp_ad13
1928 q1_ad(i-1) = q1_ad(i-1) + fx0_ad(i)
1931 IF (branch .EQ. 2)
THEN 1933 temp_ad14 = (c(i, j)+1.)*fx1_ad(i)
1934 c_ad(i, j) = c_ad(i, j) + b0(i)*temp_ad14 + (bl(i)+c(i, &
1935 & j)*b0(i))*fx1_ad(i)
1936 bl_ad(i) = bl_ad(i) + temp_ad14
1937 b0_ad(i) = b0_ad(i) + c(i, j)*temp_ad14
1940 q1_ad(i) = q1_ad(i) + fx0_ad(i)
1950 IF (branch .EQ. 0)
THEN 1956 bl_ad(i) = bl_ad(i) + b0_ad(i)
1957 br_ad(i) = br_ad(i) + b0_ad(i)
1960 al_ad(i+1) = al_ad(i+1) + br_ad(i)
1961 q1_ad(i) = q1_ad(i) - bl_ad(i) - br_ad(i)
1964 al_ad(i) = al_ad(i) + bl_ad(i)
1968 ELSE IF (branch .EQ. 5)
THEN 1971 IF (branch .NE. 0) fx1_ad(i) = fx1_ad(i) + flux_ad(i, j)
1973 IF (branch .EQ. 0)
THEN 1974 q1_ad(i-1) = q1_ad(i-1) + flux_ad(i, j)
1977 temp_ad15 = (1.-c(i, j))*fx1_ad(i)
1978 c_ad(i, j) = c_ad(i, j) - b0(i-1)*temp_ad15 - (br(i-1)-c(i, &
1979 & j)*b0(i-1))*fx1_ad(i)
1980 br_ad(i-1) = br_ad(i-1) + temp_ad15
1981 b0_ad(i-1) = b0_ad(i-1) - c(i, j)*temp_ad15
1984 q1_ad(i) = q1_ad(i) + flux_ad(i, j)
1987 temp_ad16 = (c(i, j)+1.)*fx1_ad(i)
1988 c_ad(i, j) = c_ad(i, j) + b0(i)*temp_ad16 + (bl(i)+c(i, j)*&
1990 bl_ad(i) = bl_ad(i) + temp_ad16
1991 b0_ad(i) = b0_ad(i) + c(i, j)*temp_ad16
1996 IF (branch .EQ. 0)
THEN 1999 bl_ad(i) = bl_ad(i) + b0_ad(i)
2000 br_ad(i) = br_ad(i) + b0_ad(i)
2003 al_ad(i+1) = al_ad(i+1) + br_ad(i)
2004 q1_ad(i) = q1_ad(i) - bl_ad(i) - br_ad(i)
2007 al_ad(i) = al_ad(i) + bl_ad(i)
2013 bl_ad(i) = bl_ad(i) + b0_ad(i)
2014 br_ad(i) = br_ad(i) + b0_ad(i)
2017 al_ad(i+1) = al_ad(i+1) + br_ad(i)
2018 q1_ad(i) = q1_ad(i) - bl_ad(i) - br_ad(i)
2021 al_ad(i) = al_ad(i) + bl_ad(i)
2028 IF (branch .EQ. 0)
THEN 2029 temp_ad23 = (c(i, j)+1.)*flux_ad(i, j)
2030 temp_ad24 = c(i, j)*temp_ad23
2031 q1_ad(i) = q1_ad(i) + flux_ad(i, j)
2032 c_ad(i, j) = c_ad(i, j) + (bl(i)+br(i))*temp_ad23 + (bl(i)+c&
2033 & (i, j)*(bl(i)+br(i)))*flux_ad(i, j)
2034 bl_ad(i) = bl_ad(i) + temp_ad24 + temp_ad23
2035 br_ad(i) = br_ad(i) + temp_ad24
2038 temp1 = bl(i-1) + br(i-1)
2039 temp_ad21 = (1.-c(i, j))*flux_ad(i, j)
2040 temp_ad22 = -(c(i, j)*temp_ad21)
2041 q1_ad(i-1) = q1_ad(i-1) + flux_ad(i, j)
2042 c_ad(i, j) = c_ad(i, j) - temp1*temp_ad21 - (br(i-1)-c(i, j)&
2043 & *temp1)*flux_ad(i, j)
2044 br_ad(i-1) = br_ad(i-1) + temp_ad22 + temp_ad21
2045 bl_ad(i-1) = bl_ad(i-1) + temp_ad22
2050 IF (branch .NE. 0)
THEN 2051 IF (branch .NE. 1)
THEN 2054 CALL pert_ppm_adm(3, q1(npx-2:npx), bl(npx-2:npx), bl_ad(npx&
2055 & -2:npx), br(npx-2:npx), br_ad(npx-2:npx), 1)
2057 q1_ad(npx+1) = q1_ad(npx+1) +
s11*br_ad(npx)
2058 q1_ad(npx) = q1_ad(npx) - bl_ad(npx) -
s11*br_ad(npx)
2059 dm_ad(npx+1) = dm_ad(npx+1) -
s14*br_ad(npx)
2062 xt_ad = br_ad(npx-1) + bl_ad(npx)
2065 q1_ad(npx-1) = q1_ad(npx-1) - br_ad(npx-1)
2068 IF (branch .EQ. 0)
THEN 2075 IF (branch .LT. 2)
THEN 2076 IF (branch .EQ. 0)
THEN 2079 q1_ad(npx-1) = q1_ad(npx-1) + y13_ad
2082 ELSE IF (branch .EQ. 2)
THEN 2085 q1_ad(npx-2) = q1_ad(npx-2) + y13_ad
2089 IF (branch .EQ. 0)
THEN 2090 q1_ad(npx+1) = q1_ad(npx+1) + z5_ad
2092 q1_ad(npx) = q1_ad(npx) + z5_ad
2095 IF (branch .EQ. 0)
THEN 2102 IF (branch .LT. 2)
THEN 2103 IF (branch .EQ. 0)
THEN 2106 q1_ad(npx-1) = q1_ad(npx-1) + y12_ad
2109 ELSE IF (branch .EQ. 2)
THEN 2112 q1_ad(npx-2) = q1_ad(npx-2) + y12_ad
2116 IF (branch .EQ. 0)
THEN 2117 q1_ad(npx+1) = q1_ad(npx+1) + z4_ad
2119 q1_ad(npx) = q1_ad(npx) + z4_ad
2121 temp_ad19 = 0.5*xt_ad/(dxa(npx-2, j)+dxa(npx-1, j))
2122 temp_ad20 = 0.5*xt_ad/(dxa(npx, j)+dxa(npx+1, j))
2123 q1_ad(npx-1) = q1_ad(npx-1) + (dxa(npx-1, j)*2.+dxa(npx-2, j&
2125 q1_ad(npx-2) = q1_ad(npx-2) - dxa(npx-1, j)*temp_ad19
2126 q1_ad(npx) = q1_ad(npx) + (dxa(npx, j)*2.+dxa(npx+1, j))*&
2128 q1_ad(npx+1) = q1_ad(npx+1) - dxa(npx, j)*temp_ad20
2130 xt_ad = br_ad(npx-2) + bl_ad(npx-1)
2131 q1_ad(npx-1) = q1_ad(npx-1) - bl_ad(npx-1)
2134 q1_ad(npx-2) = q1_ad(npx-2) - br_ad(npx-2)
2137 q1_ad(npx-1) = q1_ad(npx-1) +
s15*xt_ad
2138 q1_ad(npx-2) = q1_ad(npx-2) +
s11*xt_ad - bl_ad(npx-2)
2139 dm_ad(npx-2) = dm_ad(npx-2) +
s14*xt_ad
2141 al_ad(npx-2) = al_ad(npx-2) + bl_ad(npx-2)
2145 IF (branch .EQ. 0)
THEN 2148 CALL pert_ppm_adm(3, q1(0:2), bl(0:2), bl_ad(0:2), br(0:2), &
2151 al_ad(3) = al_ad(3) + br_ad(2)
2152 q1_ad(2) = q1_ad(2) - bl_ad(2) - br_ad(2)
2155 xt_ad = br_ad(1) + bl_ad(2)
2158 q1_ad(1) = q1_ad(1) +
s15*xt_ad - br_ad(1)
2160 q1_ad(2) = q1_ad(2) +
s11*xt_ad
2161 dm_ad(2) = dm_ad(2) -
s14*xt_ad
2163 xt_ad = br_ad(0) + bl_ad(1)
2164 q1_ad(1) = q1_ad(1) - bl_ad(1)
2167 q1_ad(0) = q1_ad(0) - br_ad(0)
2170 IF (branch .EQ. 0)
THEN 2177 IF (branch .LT. 2)
THEN 2178 IF (branch .EQ. 0)
THEN 2181 q1_ad(0) = q1_ad(0) + y11_ad
2184 ELSE IF (branch .EQ. 2)
THEN 2187 q1_ad(-1) = q1_ad(-1) + y11_ad
2191 IF (branch .EQ. 0)
THEN 2192 q1_ad(2) = q1_ad(2) + z3_ad
2194 q1_ad(1) = q1_ad(1) + z3_ad
2197 IF (branch .EQ. 0)
THEN 2204 IF (branch .LT. 2)
THEN 2205 IF (branch .EQ. 0)
THEN 2208 q1_ad(0) = q1_ad(0) + y10_ad
2211 ELSE IF (branch .EQ. 2)
THEN 2214 q1_ad(-1) = q1_ad(-1) + y10_ad
2218 IF (branch .EQ. 0)
THEN 2219 q1_ad(2) = q1_ad(2) + z2_ad
2221 q1_ad(1) = q1_ad(1) + z2_ad
2224 temp_ad17 = 0.5*xt_ad/(dxa(-1, j)+dxa(0, j))
2225 temp_ad18 = 0.5*xt_ad/(dxa(1, j)+dxa(2, j))
2226 q1_ad(0) = q1_ad(0) + (dxa(0, j)*2.+dxa(-1, j))*temp_ad17
2227 q1_ad(-1) = q1_ad(-1) - dxa(0, j)*temp_ad17
2228 q1_ad(1) = q1_ad(1) + (dxa(1, j)*2.+dxa(2, j))*temp_ad18
2229 q1_ad(2) = q1_ad(2) - dxa(1, j)*temp_ad18
2231 dm_ad(-1) = dm_ad(-1) +
s14*bl_ad(0)
2232 q1_ad(-1) = q1_ad(-1) +
s11*bl_ad(0)
2233 q1_ad(0) = q1_ad(0) -
s11*bl_ad(0)
2238 IF (branch .EQ. 0)
THEN 2239 arg1 = ie1 - is1 + 1
2242 CALL pert_ppm_adm(arg1, q1(is1:ie1), bl(is1:ie1), bl_ad(is1:&
2243 & ie1), br(is1:ie1), br_ad(is1:ie1), 0)
2246 IF (branch .EQ. 0)
THEN 2249 min5_ad = sign(1.d0, min5*xt)*br_ad(i)
2252 IF (branch .EQ. 0)
THEN 2262 IF (branch .EQ. 0)
THEN 2263 al_ad(i+1) = al_ad(i+1) + y5_ad
2264 q1_ad(i) = q1_ad(i) - y5_ad
2266 q1_ad(i) = q1_ad(i) + y5_ad
2267 al_ad(i+1) = al_ad(i+1) - y5_ad
2270 IF (branch .EQ. 0)
THEN 2276 min4_ad = -(sign(1.d0, min4*xt)*bl_ad(i))
2279 IF (branch .EQ. 0)
THEN 2289 IF (branch .EQ. 0)
THEN 2290 al_ad(i) = al_ad(i) + y4_ad
2291 q1_ad(i) = q1_ad(i) - y4_ad
2293 q1_ad(i) = q1_ad(i) + y4_ad
2294 al_ad(i) = al_ad(i) - y4_ad
2297 IF (branch .EQ. 0)
THEN 2298 xt_ad = xt_ad + x5_ad
2300 xt_ad = xt_ad - x5_ad
2303 dm_ad(i) = dm_ad(i) + 2.*xt_ad
2305 ELSE IF (branch .EQ. 1)
THEN 2308 min7_ad = sign(1.d0, min7*xt)*br_ad(i)
2311 IF (branch .EQ. 0)
THEN 2321 IF (branch .EQ. 0)
THEN 2322 al_ad(i+1) = al_ad(i+1) + y7_ad
2323 q1_ad(i) = q1_ad(i) - y7_ad
2325 q1_ad(i) = q1_ad(i) + y7_ad
2326 al_ad(i+1) = al_ad(i+1) - y7_ad
2329 IF (branch .EQ. 0)
THEN 2335 min6_ad = -(sign(1.d0, min6*xt)*bl_ad(i))
2338 IF (branch .EQ. 0)
THEN 2348 IF (branch .EQ. 0)
THEN 2349 al_ad(i) = al_ad(i) + y6_ad
2350 q1_ad(i) = q1_ad(i) - y6_ad
2352 q1_ad(i) = q1_ad(i) + y6_ad
2353 al_ad(i) = al_ad(i) - y6_ad
2356 IF (branch .EQ. 0)
THEN 2357 xt_ad = xt_ad + x7_ad
2359 xt_ad = xt_ad - x7_ad
2362 dm_ad(i) = dm_ad(i) +
ppm_fac*xt_ad
2367 IF (branch .LT. 2)
THEN 2368 IF (branch .EQ. 0)
THEN 2375 ELSE IF (branch .EQ. 2)
THEN 2385 IF (branch .EQ. 0)
THEN 2388 bl_ad(i) = bl_ad(i) + y9_ad
2392 IF (branch .LT. 2)
THEN 2393 IF (branch .EQ. 0)
THEN 2401 IF (branch .EQ. 2)
THEN 2409 IF (branch .LT. 2)
THEN 2410 IF (branch .EQ. 0)
THEN 2411 lac_1_ad = lac_1_ad + x10_ad
2413 pmp_1_ad = pmp_1_ad + x10_ad
2415 ELSE IF (branch .EQ. 2)
THEN 2416 lac_1_ad = lac_1_ad + x10_ad
2418 pmp_1_ad = pmp_1_ad + lac_1_ad
2419 dq_ad(i+1) = dq_ad(i+1) + 0.75*lac_1_ad
2420 dq_ad(i) = dq_ad(i) - pmp_1_ad
2422 IF (branch .EQ. 0)
THEN 2432 IF (branch .EQ. 0)
THEN 2435 br_ad(i) = br_ad(i) + y8_ad
2439 IF (branch .LT. 2)
THEN 2440 IF (branch .EQ. 0)
THEN 2448 IF (branch .EQ. 2)
THEN 2456 IF (branch .LT. 2)
THEN 2457 IF (branch .EQ. 0)
THEN 2458 lac_2_ad = lac_2_ad + x9_ad
2460 pmp_2_ad = pmp_2_ad + x9_ad
2462 ELSE IF (branch .EQ. 2)
THEN 2463 lac_2_ad = lac_2_ad + x9_ad
2465 pmp_2_ad = pmp_2_ad + lac_2_ad
2466 dq_ad(i-2) = dq_ad(i-2) - 0.75*lac_2_ad
2467 dq_ad(i-1) = dq_ad(i-1) + pmp_2_ad
2469 al_ad(i+1) = al_ad(i+1) + br_ad(i)
2470 q1_ad(i) = q1_ad(i) - bl_ad(i) - br_ad(i)
2473 al_ad(i) = al_ad(i) + bl_ad(i)
2477 q1_ad(i+1) = q1_ad(i+1) + 2.*dq_ad(i)
2478 q1_ad(i) = q1_ad(i) - 2.*dq_ad(i)
2484 q1_ad(i-1) = q1_ad(i-1) + 0.5*al_ad(i)
2485 q1_ad(i) = q1_ad(i) + 0.5*al_ad(i)
2486 dm_ad(i-1) = dm_ad(i-1) +
r3*al_ad(i)
2487 dm_ad(i) = dm_ad(i) -
r3*al_ad(i)
2491 xt = 0.25*(q1(i+1)-q1(i-1))
2492 min3_ad = sign(1.d0, min3*xt)*dm_ad(i)
2495 IF (branch .LT. 2)
THEN 2496 IF (branch .EQ. 0)
THEN 2507 IF (branch .EQ. 2)
THEN 2518 q1_ad(i) = q1_ad(i) + z1_ad
2521 IF (branch .LT. 2)
THEN 2522 IF (branch .EQ. 0)
THEN 2523 q1_ad(i+1) = q1_ad(i+1) + min8_ad
2525 q1_ad(i) = q1_ad(i) + min8_ad
2527 ELSE IF (branch .EQ. 2)
THEN 2528 q1_ad(i+1) = q1_ad(i+1) + min8_ad
2530 q1_ad(i-1) = q1_ad(i-1) + min8_ad
2533 q1_ad(i) = q1_ad(i) - y3_ad
2535 IF (branch .LT. 2)
THEN 2536 IF (branch .EQ. 0)
THEN 2537 q1_ad(i+1) = q1_ad(i+1) + max1_ad
2539 q1_ad(i) = q1_ad(i) + max1_ad
2541 ELSE IF (branch .EQ. 2)
THEN 2542 q1_ad(i+1) = q1_ad(i+1) + max1_ad
2544 q1_ad(i-1) = q1_ad(i-1) + max1_ad
2547 IF (branch .EQ. 0)
THEN 2553 q1_ad(i+1) = q1_ad(i+1) + 0.25*xt_ad
2554 q1_ad(i-1) = q1_ad(i-1) - 0.25*xt_ad
2559 IF (branch .LT. 2)
THEN 2560 IF (branch .NE. 0)
GOTO 110
2561 ELSE IF (branch .EQ. 2)
THEN 2564 IF (branch .EQ. 3)
THEN 2571 IF (branch .EQ. 0)
THEN 2578 IF (branch .EQ. 0)
THEN 2586 q1_ad(npx) = q1_ad(npx) +
c3*al_ad(npx+1)
2587 q1_ad(npx+1) = q1_ad(npx+1) +
c2*al_ad(npx+1)
2588 q1_ad(npx+2) = q1_ad(npx+2) +
c1*al_ad(npx+1)
2591 temp_ad1 = 0.5*al_ad(npx)/(dxa(npx-2, j)+dxa(npx-1, j))
2592 temp_ad2 = 0.5*al_ad(npx)/(dxa(npx, j)+dxa(npx+1, j))
2593 q1_ad(npx-1) = q1_ad(npx-1) + (dxa(npx-1, j)*2.+dxa(npx-2, j))*&
2595 q1_ad(npx-2) = q1_ad(npx-2) - dxa(npx-1, j)*temp_ad1
2596 q1_ad(npx) = q1_ad(npx) + (dxa(npx, j)*2.+dxa(npx+1, j))*temp_ad2
2597 q1_ad(npx+1) = q1_ad(npx+1) - dxa(npx, j)*temp_ad2
2600 q1_ad(npx-3) = q1_ad(npx-3) +
c1*al_ad(npx-1)
2601 q1_ad(npx-2) = q1_ad(npx-2) +
c2*al_ad(npx-1)
2602 q1_ad(npx-1) = q1_ad(npx-1) +
c3*al_ad(npx-1)
2605 IF (branch .LT. 2)
THEN 2606 IF (branch .NE. 0)
GOTO 120
2608 IF (branch .EQ. 2)
THEN 2615 IF (branch .EQ. 0)
THEN 2622 IF (branch .EQ. 0)
THEN 2630 q1_ad(1) = q1_ad(1) +
c3*al_ad(2)
2631 q1_ad(2) = q1_ad(2) +
c2*al_ad(2)
2632 q1_ad(3) = q1_ad(3) +
c1*al_ad(2)
2635 temp_ad = 0.5*al_ad(1)/(dxa(-1, j)+dxa(0, j))
2636 temp_ad0 = 0.5*al_ad(1)/(dxa(1, j)+dxa(2, j))
2637 q1_ad(0) = q1_ad(0) + (dxa(0, j)*2.+dxa(-1, j))*temp_ad
2638 q1_ad(-1) = q1_ad(-1) - dxa(0, j)*temp_ad
2639 q1_ad(1) = q1_ad(1) + (dxa(1, j)*2.+dxa(2, j))*temp_ad0
2640 q1_ad(2) = q1_ad(2) - dxa(1, j)*temp_ad0
2643 q1_ad(-2) = q1_ad(-2) +
c1*al_ad(0)
2644 q1_ad(-1) = q1_ad(-1) +
c2*al_ad(0)
2645 q1_ad(0) = q1_ad(0) +
c3*al_ad(0)
2648 IF (branch .EQ. 0)
THEN 2651 IF (branch .NE. 0)
THEN 2653 q1_ad(i-1) = q1_ad(i-1) + 0.5*al_ad(i)
2654 q1_ad(i) = q1_ad(i) + 0.5*al_ad(i)
2661 q1_ad(i-1) = q1_ad(i-1) +
p1*al_ad(i)
2662 q1_ad(i) = q1_ad(i) +
p1*al_ad(i)
2663 q1_ad(i-2) = q1_ad(i-2) +
p2*al_ad(i)
2664 q1_ad(i+1) = q1_ad(i+1) +
p2*al_ad(i)
2669 q_ad(i, j) = q_ad(i, j) + q1_ad(i)
2675 SUBROUTINE xppm(flux, q, c, iord, is, ie, isd, ied, jfirst, jlast, jsd&
2676 & , jed, npx, npy, dxa, nested, grid_type)
2678 INTEGER,
INTENT(IN) :: is, ie, isd, ied, jsd, jed
2680 INTEGER,
INTENT(IN) :: jfirst, jlast
2681 INTEGER,
INTENT(IN) :: iord
2682 INTEGER,
INTENT(IN) :: npx, npy
2683 REAL,
INTENT(IN) :: q(isd:ied, jfirst:jlast)
2685 REAL,
INTENT(IN) :: c(is:ie+1, jfirst:jlast)
2686 REAL,
INTENT(IN) :: dxa(isd:ied, jsd:jed)
2687 LOGICAL,
INTENT(IN) :: nested
2688 INTEGER,
INTENT(IN) :: grid_type
2691 REAL,
INTENT(OUT) :: flux(is:ie+1, jfirst:jlast)
2693 REAL,
DIMENSION(is-1:ie+1) :: bl, br, b0
2695 REAL,
DIMENSION(is:ie+1) :: fx0, fx1
2696 LOGICAL,
DIMENSION(is-1:ie+1) :: smt5, smt6
2697 REAL :: al(is-1:ie+2)
2698 REAL :: dm(is-2:ie+2)
2699 REAL :: dq(is-3:ie+2)
2700 INTEGER :: i, j, ie3, is1, ie1
2701 REAL :: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
2753 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 2754 IF (3 .LT. is - 1)
THEN 2759 IF (npx - 2 .GT. ie + 2)
THEN 2764 IF (npx - 3 .GT. ie + 1)
THEN 2778 IF (iord .LT. 8 .OR. iord .EQ. 333)
THEN 2782 al(i) =
p1*(q1(i-1)+q1(i)) +
p2*(q1(i-2)+q1(i+1))
2784 IF (iord .EQ. 7)
THEN 2786 IF (al(i) .LT. 0.) al(i) = 0.5*(q1(i-1)+q1(i))
2789 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 2791 al(0) =
c1*q1(-2) +
c2*q1(-1) +
c3*q1(0)
2792 al(1) = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-&
2793 & 1))/(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)&
2794 & -dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
2795 al(2) =
c3*q1(1) +
c2*q1(2) +
c1*q1(3)
2796 IF (iord .EQ. 7)
THEN 2797 IF (0. .LT. al(0))
THEN 2802 IF (0. .LT. al(1))
THEN 2807 IF (0. .LT. al(2))
THEN 2814 IF (ie + 1 .EQ. npx)
THEN 2815 al(npx-1) =
c1*q1(npx-3) +
c2*q1(npx-2) +
c3*q1(npx-1)
2816 al(npx) = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-&
2817 & dxa(npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((&
2818 & 2.*dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1&
2819 & ))/(dxa(npx, j)+dxa(npx+1, j)))
2820 al(npx+1) =
c3*q1(npx) +
c2*q1(npx+1) +
c1*q1(npx+2)
2821 IF (iord .EQ. 7)
THEN 2822 IF (0. .LT. al(npx-1))
THEN 2823 al(npx-1) = al(npx-1)
2827 IF (0. .LT. al(npx))
THEN 2832 IF (0. .LT. al(npx+1))
THEN 2833 al(npx+1) = al(npx+1)
2840 IF (iord .EQ. 1)
THEN 2842 IF (c(i, j) .GT. 0.)
THEN 2843 flux(i, j) = q1(i-1)
2848 ELSE IF (iord .EQ. 2)
THEN 2854 IF (xt .GT. 0.)
THEN 2856 flux(i, j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-&
2860 flux(i, j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-&
2864 ELSE IF (iord .EQ. 333)
THEN 2873 IF (xt .GT. 0.)
THEN 2874 flux(i, j) = (2.0*q1(i)+5.0*q1(i-1)-q1(i-2))/6.0 - 0.5*xt*&
2875 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i)-2.0*q1(i-1)+q1(i-2))
2877 flux(i, j) = (2.0*q1(i-1)+5.0*q1(i)-q1(i+1))/6.0 - 0.5*xt*&
2878 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i+1)-2.0*q1(i)+q1(i-1))
2881 ELSE IF (iord .EQ. 3)
THEN 2883 bl(i) = al(i) - q1(i)
2884 br(i) = al(i+1) - q1(i)
2885 b0(i) = bl(i) + br(i)
2886 IF (b0(i) .GE. 0.)
THEN 2891 IF (bl(i) - br(i) .GE. 0.)
THEN 2896 smt5(i) = x0 .LT. xt
2897 smt6(i) = 3.*x0 .LT. xt
2904 IF (xt .GT. 0.)
THEN 2906 IF (smt6(i-1) .OR. smt5(i))
THEN 2907 fx1(i) = br(i-1) - xt*b0(i-1)
2908 ELSE IF (smt5(i-1))
THEN 2909 IF (bl(i-1) .GE. 0.)
THEN 2914 IF (br(i-1) .GE. 0.)
THEN 2919 IF (x2 .GT. y1)
THEN 2925 fx1(i) = sign(min1, br(i-1))
2929 IF (smt6(i) .OR. smt5(i-1))
THEN 2930 fx1(i) = bl(i) + xt*b0(i)
2931 ELSE IF (smt5(i))
THEN 2932 IF (bl(i) .GE. 0.)
THEN 2937 IF (br(i) .GE. 0.)
THEN 2942 IF (x3 .GT. y2)
THEN 2947 fx1(i) = sign(min2, bl(i))
2950 IF (xt .GE. 0.)
THEN 2955 flux(i, j) = fx0(i) + (1.-abs0)*fx1(i)
2957 ELSE IF (iord .EQ. 4)
THEN 2959 bl(i) = al(i) - q1(i)
2960 br(i) = al(i+1) - q1(i)
2961 b0(i) = bl(i) + br(i)
2962 IF (b0(i) .GE. 0.)
THEN 2967 IF (bl(i) - br(i) .GE. 0.)
THEN 2972 smt5(i) = x0 .LT. xt
2973 smt6(i) = 3.*x0 .LT. xt
2980 IF (c(i, j) .GT. 0.)
THEN 2982 IF (smt6(i-1) .OR. smt5(i)) fx1(i) = (1.-c(i, j))*(br(i-1)&
2986 IF (smt6(i) .OR. smt5(i-1)) fx1(i) = (1.+c(i, j))*(bl(i)+c&
2989 flux(i, j) = fx0(i) + fx1(i)
2993 IF (iord .EQ. 5)
THEN 2995 bl(i) = al(i) - q1(i)
2996 br(i) = al(i+1) - q1(i)
2997 b0(i) = bl(i) + br(i)
2998 smt5(i) = bl(i)*br(i) .LT. 0.
3002 bl(i) = al(i) - q1(i)
3003 br(i) = al(i+1) - q1(i)
3004 b0(i) = bl(i) + br(i)
3005 IF (3.*b0(i) .GE. 0.)
THEN 3010 IF (bl(i) - br(i) .GE. 0.)
THEN 3011 abs4 = bl(i) - br(i)
3013 abs4 = -(bl(i)-br(i))
3015 smt5(i) = abs1 .LT. abs4
3020 IF (c(i, j) .GT. 0.)
THEN 3021 fx1(i) = (1.-c(i, j))*(br(i-1)-c(i, j)*b0(i-1))
3022 flux(i, j) = q1(i-1)
3024 fx1(i) = (1.+c(i, j))*(bl(i)+c(i, j)*b0(i))
3027 IF (smt5(i-1) .OR. smt5(i)) flux(i, j) = flux(i, j) + fx1(i)
3036 xt = 0.25*(q1(i+1)-q1(i-1))
3037 IF (xt .GE. 0.)
THEN 3042 IF (q1(i-1) .LT. q1(i))
THEN 3043 IF (q1(i) .LT. q1(i+1))
THEN 3048 ELSE IF (q1(i-1) .LT. q1(i+1))
THEN 3054 IF (q1(i-1) .GT. q1(i))
THEN 3055 IF (q1(i) .GT. q1(i+1))
THEN 3060 ELSE IF (q1(i-1) .GT. q1(i+1))
THEN 3066 IF (x4 .GT. y3)
THEN 3067 IF (y3 .GT. z1)
THEN 3072 ELSE IF (x4 .GT. z1)
THEN 3077 dm(i) = sign(min3, xt)
3080 al(i) = 0.5*(q1(i-1)+q1(i)) +
r3*(dm(i-1)-dm(i))
3082 IF (iord .EQ. 8)
THEN 3085 IF (xt .GE. 0.)
THEN 3090 IF (al(i) - q1(i) .GE. 0.)
THEN 3095 IF (x5 .GT. y4)
THEN 3100 bl(i) = -sign(min4, xt)
3101 IF (xt .GE. 0.)
THEN 3106 IF (al(i+1) - q1(i) .GE. 0.)
THEN 3107 y5 = al(i+1) - q1(i)
3109 y5 = -(al(i+1)-q1(i))
3111 IF (x6 .GT. y5)
THEN 3116 br(i) = sign(min5, xt)
3118 ELSE IF (iord .EQ. 11)
THEN 3122 IF (xt .GE. 0.)
THEN 3127 IF (al(i) - q1(i) .GE. 0.)
THEN 3132 IF (x7 .GT. y6)
THEN 3137 bl(i) = -sign(min6, xt)
3138 IF (xt .GE. 0.)
THEN 3143 IF (al(i+1) - q1(i) .GE. 0.)
THEN 3144 y7 = al(i+1) - q1(i)
3146 y7 = -(al(i+1)-q1(i))
3148 IF (x8 .GT. y7)
THEN 3153 br(i) = sign(min7, xt)
3157 dq(i) = 2.*(q1(i+1)-q1(i))
3160 bl(i) = al(i) - q1(i)
3161 br(i) = al(i+1) - q1(i)
3162 IF (dm(i-1) .GE. 0.)
THEN 3167 IF (dm(i) .GE. 0.)
THEN 3172 IF (dm(i+1) .GE. 0.)
THEN 3177 IF (abs2 + abs5 + abs7 .LT.
near_zero)
THEN 3181 IF (3.*(bl(i)+br(i)) .GE. 0.)
THEN 3182 abs3 = 3.*(bl(i)+br(i))
3184 abs3 = -(3.*(bl(i)+br(i)))
3186 IF (bl(i) - br(i) .GE. 0.)
THEN 3187 abs6 = bl(i) - br(i)
3189 abs6 = -(bl(i)-br(i))
3191 IF (abs3 .GT. abs6)
THEN 3193 lac_2 = pmp_2 - 0.75*dq(i-2)
3194 IF (0. .LT. pmp_2)
THEN 3195 IF (pmp_2 .LT. lac_2)
THEN 3200 ELSE IF (0. .LT. lac_2)
THEN 3205 IF (0. .GT. pmp_2)
THEN 3206 IF (pmp_2 .GT. lac_2)
THEN 3211 ELSE IF (0. .GT. lac_2)
THEN 3216 IF (br(i) .LT. y14)
THEN 3221 IF (x9 .GT. y8)
THEN 3227 lac_1 = pmp_1 + 0.75*dq(i+1)
3228 IF (0. .LT. pmp_1)
THEN 3229 IF (pmp_1 .LT. lac_1)
THEN 3234 ELSE IF (0. .LT. lac_1)
THEN 3239 IF (0. .GT. pmp_1)
THEN 3240 IF (pmp_1 .GT. lac_1)
THEN 3245 ELSE IF (0. .GT. lac_1)
THEN 3250 IF (bl(i) .LT. y15)
THEN 3255 IF (x10 .GT. y9)
THEN 3265 IF (iord .EQ. 9 .OR. iord .EQ. 13)
THEN 3266 arg1 = ie1 - is1 + 1
3267 CALL pert_ppm(arg1, q1(is1:ie1), bl(is1:ie1), br(is1:ie1), 0)
3269 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 3271 bl(0) =
s14*dm(-1) +
s11*(q1(-1)-q1(0))
3272 xt = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-1))&
3273 & /(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)-&
3274 & dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
3275 IF (q1(1) .GT. q1(2))
THEN 3280 IF (q1(-1) .GT. q1(0))
THEN 3281 IF (q1(0) .GT. z2)
THEN 3286 ELSE IF (q1(-1) .GT. z2)
THEN 3291 IF (xt .LT. y10)
THEN 3296 IF (q1(1) .LT. q1(2))
THEN 3301 IF (q1(-1) .LT. q1(0))
THEN 3302 IF (q1(0) .LT. z3)
THEN 3307 ELSE IF (q1(-1) .LT. z3)
THEN 3312 IF (xt .GT. y11)
THEN 3323 br(2) = al(3) - q1(2)
3324 CALL pert_ppm(3, q1(0:2), bl(0:2), br(0:2), 1)
3326 IF (ie + 1 .EQ. npx)
THEN 3327 bl(npx-2) = al(npx-2) - q1(npx-2)
3328 xt =
s15*q1(npx-1) +
s11*q1(npx-2) +
s14*dm(npx-2)
3329 br(npx-2) = xt - q1(npx-2)
3330 bl(npx-1) = xt - q1(npx-1)
3331 xt = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-dxa(&
3332 & npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((2.*&
3333 & dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1))/&
3334 & (dxa(npx, j)+dxa(npx+1, j)))
3335 IF (q1(npx) .GT. q1(npx+1))
THEN 3340 IF (q1(npx-2) .GT. q1(npx-1))
THEN 3341 IF (q1(npx-1) .GT. z4)
THEN 3346 ELSE IF (q1(npx-2) .GT. z4)
THEN 3351 IF (xt .LT. y12)
THEN 3356 IF (q1(npx) .LT. q1(npx+1))
THEN 3361 IF (q1(npx-2) .LT. q1(npx-1))
THEN 3362 IF (q1(npx-1) .LT. z5)
THEN 3367 ELSE IF (q1(npx-2) .LT. z5)
THEN 3372 IF (xt .GT. y13)
THEN 3378 br(npx-1) = xt - q1(npx-1)
3379 bl(npx) = xt - q1(npx)
3380 br(npx) =
s11*(q1(npx+1)-q1(npx)) -
s14*dm(npx+1)
3381 CALL pert_ppm(3, q1(npx-2:npx), bl(npx-2:npx), br(npx-2:npx)&
3386 IF (c(i, j) .GT. 0.)
THEN 3387 flux(i, j) = q1(i-1) + (1.-c(i, j))*(br(i-1)-c(i, j)*(bl(i-1&
3390 flux(i, j) = q1(i) + (1.+c(i, j))*(bl(i)+c(i, j)*(bl(i)+br(i&
3417 SUBROUTINE yppm_adm(flux, flux_ad, q, q_ad, c, c_ad, jord, ifirst, &
3418 & ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
3421 INTEGER,
INTENT(IN) :: ifirst, ilast
3422 INTEGER,
INTENT(IN) :: isd, ied, js, je, jsd, jed
3423 INTEGER,
INTENT(IN) :: jord
3424 INTEGER,
INTENT(IN) :: npx, npy
3425 REAL,
INTENT(IN) :: q(ifirst:ilast, jsd:jed)
3426 REAL :: q_ad(ifirst:ilast, jsd:jed)
3428 REAL,
INTENT(IN) :: c(isd:ied, js:je+1)
3429 REAL :: c_ad(isd:ied, js:je+1)
3431 REAL :: flux(ifirst:ilast, js:je+1)
3432 REAL :: flux_ad(ifirst:ilast, js:je+1)
3433 REAL,
INTENT(IN) :: dya(isd:ied, jsd:jed)
3434 LOGICAL,
INTENT(IN) :: nested
3435 INTEGER,
INTENT(IN) :: grid_type
3437 REAL :: dm(ifirst:ilast, js-2:je+2)
3438 REAL :: dm_ad(ifirst:ilast, js-2:je+2)
3439 REAL :: al(ifirst:ilast, js-1:je+2)
3440 REAL :: al_ad(ifirst:ilast, js-1:je+2)
3441 REAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: bl, br, b0
3442 REAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: bl_ad, br_ad, b0_ad
3443 REAL :: dq(ifirst:ilast, js-3:je+2)
3444 REAL :: dq_ad(ifirst:ilast, js-3:je+2)
3445 REAL,
DIMENSION(ifirst:ilast) :: fx0, fx1
3446 REAL,
DIMENSION(ifirst:ilast) :: fx0_ad, fx1_ad
3447 LOGICAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: smt5, smt6
3448 REAL :: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
3449 REAL :: xt_ad, qtmp_ad, pmp_1_ad, lac_1_ad, pmp_2_ad, lac_2_ad
3450 INTEGER :: i, j, js1, je3, je1
3572 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 3573 IF (3 .LT. js - 1)
THEN 3578 IF (npy - 2 .GT. je + 2)
THEN 3583 IF (npy - 3 .GT. je + 1)
THEN 3597 IF (jord .LT. 8 .OR. jord .EQ. 333)
THEN 3600 al(i, j) =
p1*(q(i, j-1)+q(i, j)) +
p2*(q(i, j-2)+q(i, j+1))
3603 IF (jord .EQ. 7)
THEN 3606 IF (al(i, j) .LT. 0.)
THEN 3607 al(i, j) = 0.5*(q(i, j)+q(i, j+1))
3618 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 3621 al(i, 0) =
c1*q(i, -2) +
c2*q(i, -1) +
c3*q(i, 0)
3622 al(i, 1) = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)&
3623 & *q(i, -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2)&
3624 & )*q(i, 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
3625 al(i, 2) =
c3*q(i, 1) +
c2*q(i, 2) +
c1*q(i, 3)
3627 IF (jord .EQ. 7)
THEN 3629 IF (0. .LT. al(i, 0))
THEN 3636 IF (0. .LT. al(i, 1))
THEN 3643 IF (0. .LT. al(i, 2))
THEN 3658 IF (je + 1 .EQ. npy)
THEN 3660 al(i, npy-1) =
c1*q(i, npy-3) +
c2*q(i, npy-2) +
c3*q(i, npy&
3662 al(i, npy) = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy&
3663 & -1)-dya(i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1&
3664 & ))+((2.*dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q&
3665 & (i, npy+1))/(dya(i, npy)+dya(i, npy+1)))
3666 al(i, npy+1) =
c3*q(i, npy) +
c2*q(i, npy+1) +
c1*q(i, npy+2&
3669 IF (jord .EQ. 7)
THEN 3671 IF (0. .LT. al(i, npy-1))
THEN 3673 al(i, npy-1) = al(i, npy-1)
3678 IF (0. .LT. al(i, npy))
THEN 3680 al(i, npy) = al(i, npy)
3685 IF (0. .LT. al(i, npy+1))
THEN 3687 al(i, npy+1) = al(i, npy+1)
3703 IF (jord .EQ. 1)
THEN 3706 IF (c(i, j) .GT. 0.)
THEN 3714 DO i=ilast,ifirst,-1
3716 IF (branch .EQ. 0)
THEN 3717 q_ad(i, j) = q_ad(i, j) + flux_ad(i, j)
3720 q_ad(i, j-1) = q_ad(i, j-1) + flux_ad(i, j)
3726 ELSE IF (jord .EQ. 2)
THEN 3733 IF (xt .GT. 0.)
THEN 3742 DO i=ilast,ifirst,-1
3744 IF (branch .EQ. 0)
THEN 3747 temp0 = al(i, j) + al(i, j+1) - 2*qtmp
3748 temp_ad5 = (xt+1.)*flux_ad(i, j)
3749 temp_ad6 = xt*temp_ad5
3750 qtmp_ad = flux_ad(i, j) - temp_ad5 - 2*temp_ad6
3751 xt_ad = temp0*temp_ad5 + (al(i, j)-qtmp+xt*temp0)*flux_ad(&
3753 al_ad(i, j) = al_ad(i, j) + temp_ad6 + temp_ad5
3754 al_ad(i, j+1) = al_ad(i, j+1) + temp_ad6
3756 q_ad(i, j) = q_ad(i, j) + qtmp_ad
3760 temp = al(i, j-1) + al(i, j) - 2*qtmp
3761 temp_ad3 = (1.-xt)*flux_ad(i, j)
3762 temp_ad4 = -(xt*temp_ad3)
3763 qtmp_ad = flux_ad(i, j) - temp_ad3 - 2*temp_ad4
3764 xt_ad = -(temp*temp_ad3) - (al(i, j)-qtmp-xt*temp)*flux_ad&
3766 al_ad(i, j) = al_ad(i, j) + temp_ad4 + temp_ad3
3767 al_ad(i, j-1) = al_ad(i, j-1) + temp_ad4
3769 q_ad(i, j-1) = q_ad(i, j-1) + qtmp_ad
3771 c_ad(i, j) = c_ad(i, j) + xt_ad
3774 ELSE IF (jord .EQ. 333)
THEN 3780 IF (xt .GT. 0.)
THEN 3788 DO i=ilast,ifirst,-1
3790 IF (branch .EQ. 0)
THEN 3792 temp_ad10 = flux_ad(i, j)/6.0
3793 temp_ad11 = -(0.5*xt*flux_ad(i, j))
3794 temp_ad12 = xt**2*flux_ad(i, j)/6.0
3795 q_ad(i, j-1) = q_ad(i, j-1) + temp_ad12 - temp_ad11 + 2.0*&
3797 q_ad(i, j) = q_ad(i, j) + temp_ad11 - 2.0*temp_ad12 + 5.0*&
3799 q_ad(i, j+1) = q_ad(i, j+1) + temp_ad12 - temp_ad10
3800 xt_ad = ((q(i, j+1)-2.0*q(i, j)+q(i, j-1))*2*xt/6.0-0.5*(q&
3801 & (i, j)-q(i, j-1)))*flux_ad(i, j)
3805 temp_ad7 = flux_ad(i, j)/6.0
3806 temp_ad8 = -(0.5*xt*flux_ad(i, j))
3807 temp_ad9 = xt**2*flux_ad(i, j)/6.0
3808 q_ad(i, j) = q_ad(i, j) + temp_ad9 + temp_ad8 + 2.0*&
3810 q_ad(i, j-1) = q_ad(i, j-1) + 5.0*temp_ad7 - temp_ad8 - &
3812 q_ad(i, j-2) = q_ad(i, j-2) + temp_ad9 - temp_ad7
3813 xt_ad = ((q(i, j)-2.0*q(i, j-1)+q(i, j-2))*2*xt/6.0-0.5*(q&
3814 & (i, j)-q(i, j-1)))*flux_ad(i, j)
3817 c_ad(i, j) = c_ad(i, j) + xt_ad
3821 ELSE IF (jord .EQ. 3)
THEN 3824 bl(i, j) = al(i, j) - q(i, j)
3825 br(i, j) = al(i, j+1) - q(i, j)
3826 b0(i, j) = bl(i, j) + br(i, j)
3827 IF (b0(i, j) .GE. 0.)
THEN 3832 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 3833 xt = bl(i, j) - br(i, j)
3835 xt = -(bl(i, j)-br(i, j))
3837 smt5(i, j) = x0 .LT. xt
3838 smt6(i, j) = 3.*x0 .LT. xt
3848 IF (xt .GT. 0.)
THEN 3849 IF (smt6(i, j-1) .OR. smt5(i, j))
THEN 3851 fx1(i) = br(i, j-1) - xt*b0(i, j-1)
3853 ELSE IF (smt5(i, j-1))
THEN 3854 IF (bl(i, j-1) .GE. 0.)
THEN 3861 IF (br(i, j-1) .GE. 0.)
THEN 3868 IF (x1 .GT. y1)
THEN 3879 fx1(i) = sign(min1, br(i, j-1))
3884 ELSE IF (smt6(i, j) .OR. smt5(i, j-1))
THEN 3886 fx1(i) = bl(i, j) + xt*b0(i, j)
3888 ELSE IF (smt5(i, j))
THEN 3889 IF (bl(i, j) .GE. 0.)
THEN 3896 IF (br(i, j) .GE. 0.)
THEN 3903 IF (x2 .GT. y2)
THEN 3913 fx1(i) = sign(min2, bl(i, j))
3918 IF (xt .GE. 0.)
THEN 3935 DO i=ilast,ifirst,-1
3936 fx0_ad(i) = fx0_ad(i) + flux_ad(i, j)
3937 abs0_ad = -(fx1(i)*flux_ad(i, j))
3938 fx1_ad(i) = fx1_ad(i) + (1.-abs0)*flux_ad(i, j)
3942 IF (branch .EQ. 0)
THEN 3950 IF (branch .LT. 3)
THEN 3951 IF (branch .NE. 0)
THEN 3952 IF (branch .EQ. 1)
THEN 3954 min2_ad = sign(1.d0, min2*bl(i, j))*fx1_ad(i)
3957 IF (branch .EQ. 0)
THEN 3967 IF (branch .EQ. 0)
THEN 3968 br_ad(i, j) = br_ad(i, j) + y2_ad
3970 br_ad(i, j) = br_ad(i, j) - y2_ad
3973 IF (branch .EQ. 0)
THEN 3974 bl_ad(i, j) = bl_ad(i, j) + x2_ad
3976 bl_ad(i, j) = bl_ad(i, j) - x2_ad
3980 bl_ad(i, j) = bl_ad(i, j) + fx1_ad(i)
3981 xt_ad = xt_ad + b0(i, j)*fx1_ad(i)
3982 b0_ad(i, j) = b0_ad(i, j) + xt*fx1_ad(i)
3986 q_ad(i, j) = q_ad(i, j) + fx0_ad(i)
3989 IF (branch .NE. 3)
THEN 3990 IF (branch .EQ. 4)
THEN 3992 min1_ad = sign(1.d0, min1*br(i, j-1))*fx1_ad(i)
3995 IF (branch .EQ. 0)
THEN 4005 IF (branch .EQ. 0)
THEN 4006 br_ad(i, j-1) = br_ad(i, j-1) + y1_ad
4008 br_ad(i, j-1) = br_ad(i, j-1) - y1_ad
4011 IF (branch .EQ. 0)
THEN 4012 bl_ad(i, j-1) = bl_ad(i, j-1) + x1_ad
4014 bl_ad(i, j-1) = bl_ad(i, j-1) - x1_ad
4018 br_ad(i, j-1) = br_ad(i, j-1) + fx1_ad(i)
4019 xt_ad = xt_ad - b0(i, j-1)*fx1_ad(i)
4020 b0_ad(i, j-1) = b0_ad(i, j-1) - xt*fx1_ad(i)
4024 q_ad(i, j-1) = q_ad(i, j-1) + fx0_ad(i)
4027 c_ad(i, j) = c_ad(i, j) + xt_ad
4029 DO i=ilast,ifirst,-1
4036 DO i=ilast,ifirst,-1
4037 bl_ad(i, j) = bl_ad(i, j) + b0_ad(i, j)
4038 br_ad(i, j) = br_ad(i, j) + b0_ad(i, j)
4040 al_ad(i, j+1) = al_ad(i, j+1) + br_ad(i, j)
4041 q_ad(i, j) = q_ad(i, j) - bl_ad(i, j) - br_ad(i, j)
4043 al_ad(i, j) = al_ad(i, j) + bl_ad(i, j)
4047 ELSE IF (jord .EQ. 4)
THEN 4050 bl(i, j) = al(i, j) - q(i, j)
4051 br(i, j) = al(i, j+1) - q(i, j)
4052 b0(i, j) = bl(i, j) + br(i, j)
4053 IF (b0(i, j) .GE. 0.)
THEN 4058 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 4059 xt = bl(i, j) - br(i, j)
4061 xt = -(bl(i, j)-br(i, j))
4063 smt5(i, j) = x0 .LT. xt
4064 smt6(i, j) = 3.*x0 .LT. xt
4070 IF (c(i, j) .GT. 0.)
THEN 4071 IF (smt6(i, j-1) .OR. smt5(i, j))
THEN 4076 ELSE IF (smt6(i, j) .OR. smt5(i, j-1))
THEN 4089 DO i=ilast,ifirst,-1
4090 fx0_ad(i) = fx0_ad(i) + flux_ad(i, j)
4091 fx1_ad(i) = fx1_ad(i) + flux_ad(i, j)
4094 IF (branch .LT. 2)
THEN 4095 IF (branch .EQ. 0)
THEN 4096 temp_ad13 = (1.-c(i, j))*fx1_ad(i)
4097 c_ad(i, j) = c_ad(i, j) - b0(i, j-1)*temp_ad13 - (br(i, &
4098 & j-1)-c(i, j)*b0(i, j-1))*fx1_ad(i)
4099 br_ad(i, j-1) = br_ad(i, j-1) + temp_ad13
4100 b0_ad(i, j-1) = b0_ad(i, j-1) - c(i, j)*temp_ad13
4103 q_ad(i, j-1) = q_ad(i, j-1) + fx0_ad(i)
4106 IF (branch .EQ. 2)
THEN 4107 temp_ad14 = (c(i, j)+1.)*fx1_ad(i)
4108 c_ad(i, j) = c_ad(i, j) + b0(i, j)*temp_ad14 + (bl(i, j)&
4109 & +c(i, j)*b0(i, j))*fx1_ad(i)
4110 bl_ad(i, j) = bl_ad(i, j) + temp_ad14
4111 b0_ad(i, j) = b0_ad(i, j) + c(i, j)*temp_ad14
4114 q_ad(i, j) = q_ad(i, j) + fx0_ad(i)
4118 DO i=ilast,ifirst,-1
4124 DO i=ilast,ifirst,-1
4125 bl_ad(i, j) = bl_ad(i, j) + b0_ad(i, j)
4126 br_ad(i, j) = br_ad(i, j) + b0_ad(i, j)
4128 al_ad(i, j+1) = al_ad(i, j+1) + br_ad(i, j)
4129 q_ad(i, j) = q_ad(i, j) - bl_ad(i, j) - br_ad(i, j)
4131 al_ad(i, j) = al_ad(i, j) + bl_ad(i, j)
4137 IF (jord .EQ. 5)
THEN 4140 bl(i, j) = al(i, j) - q(i, j)
4141 br(i, j) = al(i, j+1) - q(i, j)
4142 b0(i, j) = bl(i, j) + br(i, j)
4143 smt5(i, j) = bl(i, j)*br(i, j) .LT. 0.
4150 bl(i, j) = al(i, j) - q(i, j)
4151 br(i, j) = al(i, j+1) - q(i, j)
4152 b0(i, j) = bl(i, j) + br(i, j)
4153 IF (3.*b0(i, j) .GE. 0.)
THEN 4156 abs1 = -(3.*b0(i, j))
4158 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 4159 abs4 = bl(i, j) - br(i, j)
4161 abs4 = -(bl(i, j)-br(i, j))
4163 smt5(i, j) = abs1 .LT. abs4
4171 IF (c(i, j) .GT. 0.)
THEN 4176 IF (smt5(i, j-1) .OR. smt5(i, j))
THEN 4188 DO i=ilast,ifirst,-1
4190 IF (branch .NE. 0) fx1_ad(i) = fx1_ad(i) + flux_ad(i, j)
4192 IF (branch .EQ. 0)
THEN 4193 q_ad(i, j-1) = q_ad(i, j-1) + flux_ad(i, j)
4195 temp_ad15 = (1.-c(i, j))*fx1_ad(i)
4196 c_ad(i, j) = c_ad(i, j) - b0(i, j-1)*temp_ad15 - (br(i, j-&
4197 & 1)-c(i, j)*b0(i, j-1))*fx1_ad(i)
4198 br_ad(i, j-1) = br_ad(i, j-1) + temp_ad15
4199 b0_ad(i, j-1) = b0_ad(i, j-1) - c(i, j)*temp_ad15
4202 q_ad(i, j) = q_ad(i, j) + flux_ad(i, j)
4204 temp_ad16 = (c(i, j)+1.)*fx1_ad(i)
4205 c_ad(i, j) = c_ad(i, j) + b0(i, j)*temp_ad16 + (bl(i, j)+c&
4206 & (i, j)*b0(i, j))*fx1_ad(i)
4207 bl_ad(i, j) = bl_ad(i, j) + temp_ad16
4208 b0_ad(i, j) = b0_ad(i, j) + c(i, j)*temp_ad16
4214 IF (branch .EQ. 0)
THEN 4217 DO i=ilast,ifirst,-1
4218 bl_ad(i, j) = bl_ad(i, j) + b0_ad(i, j)
4219 br_ad(i, j) = br_ad(i, j) + b0_ad(i, j)
4221 al_ad(i, j+1) = al_ad(i, j+1) + br_ad(i, j)
4222 q_ad(i, j) = q_ad(i, j) - bl_ad(i, j) - br_ad(i, j)
4224 al_ad(i, j) = al_ad(i, j) + bl_ad(i, j)
4231 DO i=ilast,ifirst,-1
4232 bl_ad(i, j) = bl_ad(i, j) + b0_ad(i, j)
4233 br_ad(i, j) = br_ad(i, j) + b0_ad(i, j)
4235 al_ad(i, j+1) = al_ad(i, j+1) + br_ad(i, j)
4236 q_ad(i, j) = q_ad(i, j) - bl_ad(i, j) - br_ad(i, j)
4238 al_ad(i, j) = al_ad(i, j) + bl_ad(i, j)
4245 IF (branch .LT. 2)
THEN 4246 IF (branch .EQ. 0)
THEN 4247 DO i=ilast,ifirst,-1
4249 IF (branch .NE. 0) al_ad(i, npy+1) = 0.0
4251 IF (branch .NE. 0) al_ad(i, npy) = 0.0
4253 IF (branch .NE. 0) al_ad(i, npy-1) = 0.0
4256 DO i=ilast,ifirst,-1
4257 q_ad(i, npy) = q_ad(i, npy) +
c3*al_ad(i, npy+1)
4258 q_ad(i, npy+1) = q_ad(i, npy+1) +
c2*al_ad(i, npy+1)
4259 q_ad(i, npy+2) = q_ad(i, npy+2) +
c1*al_ad(i, npy+1)
4260 al_ad(i, npy+1) = 0.0
4261 temp_ad1 = 0.5*al_ad(i, npy)/(dya(i, npy-2)+dya(i, npy-1))
4262 temp_ad2 = 0.5*al_ad(i, npy)/(dya(i, npy)+dya(i, npy+1))
4263 q_ad(i, npy-1) = q_ad(i, npy-1) + (dya(i, npy-1)*2.+dya(i, npy&
4265 q_ad(i, npy-2) = q_ad(i, npy-2) - dya(i, npy-1)*temp_ad1
4266 q_ad(i, npy) = q_ad(i, npy) + (dya(i, npy)*2.+dya(i, npy+1))*&
4268 q_ad(i, npy+1) = q_ad(i, npy+1) - dya(i, npy)*temp_ad2
4270 q_ad(i, npy-3) = q_ad(i, npy-3) +
c1*al_ad(i, npy-1)
4271 q_ad(i, npy-2) = q_ad(i, npy-2) +
c2*al_ad(i, npy-1)
4272 q_ad(i, npy-1) = q_ad(i, npy-1) +
c3*al_ad(i, npy-1)
4273 al_ad(i, npy-1) = 0.0
4275 ELSE IF (branch .NE. 2)
THEN 4279 IF (branch .EQ. 0)
THEN 4280 DO i=ilast,ifirst,-1
4282 IF (branch .NE. 0) al_ad(i, 2) = 0.0
4284 IF (branch .NE. 0) al_ad(i, 1) = 0.0
4286 IF (branch .NE. 0) al_ad(i, 0) = 0.0
4288 ELSE IF (branch .NE. 1)
THEN 4291 DO i=ilast,ifirst,-1
4292 q_ad(i, 1) = q_ad(i, 1) +
c3*al_ad(i, 2)
4293 q_ad(i, 2) = q_ad(i, 2) +
c2*al_ad(i, 2)
4294 q_ad(i, 3) = q_ad(i, 3) +
c1*al_ad(i, 2)
4296 temp_ad = 0.5*al_ad(i, 1)/(dya(i, -1)+dya(i, 0))
4297 temp_ad0 = 0.5*al_ad(i, 1)/(dya(i, 1)+dya(i, 2))
4298 q_ad(i, 0) = q_ad(i, 0) + (dya(i, 0)*2.+dya(i, -1))*temp_ad
4299 q_ad(i, -1) = q_ad(i, -1) - dya(i, 0)*temp_ad
4300 q_ad(i, 1) = q_ad(i, 1) + (dya(i, 1)*2.+dya(i, 2))*temp_ad0
4301 q_ad(i, 2) = q_ad(i, 2) - dya(i, 1)*temp_ad0
4303 q_ad(i, -2) = q_ad(i, -2) +
c1*al_ad(i, 0)
4304 q_ad(i, -1) = q_ad(i, -1) +
c2*al_ad(i, 0)
4305 q_ad(i, 0) = q_ad(i, 0) +
c3*al_ad(i, 0)
4309 IF (branch .EQ. 0)
THEN 4311 DO i=ilast,ifirst,-1
4313 IF (branch .NE. 0)
THEN 4314 q_ad(i, j) = q_ad(i, j) + 0.5*al_ad(i, j)
4315 q_ad(i, j+1) = q_ad(i, j+1) + 0.5*al_ad(i, j)
4322 DO i=ilast,ifirst,-1
4323 q_ad(i, j-1) = q_ad(i, j-1) +
p1*al_ad(i, j)
4324 q_ad(i, j) = q_ad(i, j) +
p1*al_ad(i, j)
4325 q_ad(i, j-2) = q_ad(i, j-2) +
p2*al_ad(i, j)
4326 q_ad(i, j+1) = q_ad(i, j+1) +
p2*al_ad(i, j)
4336 xt = 0.25*(q(i, j+1)-q(i, j-1))
4337 IF (xt .GE. 0.)
THEN 4344 IF (q(i, j-1) .LT. q(i, j))
THEN 4345 IF (q(i, j) .LT. q(i, j+1))
THEN 4352 ELSE IF (q(i, j-1) .LT. q(i, j+1))
THEN 4360 IF (q(i, j-1) .GT. q(i, j))
THEN 4361 IF (q(i, j) .GT. q(i, j+1))
THEN 4368 ELSE IF (q(i, j-1) .GT. q(i, j+1))
THEN 4376 IF (x3 .GT. y3)
THEN 4377 IF (y3 .GT. z1)
THEN 4386 ELSE IF (x3 .GT. z1)
THEN 4395 dm(i, j) = sign(min3, xt)
4400 al(i, j) = 0.5*(q(i, j-1)+q(i, j)) +
r3*(dm(i, j-1)-dm(i, j))
4403 IF (jord .EQ. 8)
THEN 4407 IF (xt .GE. 0.)
THEN 4414 IF (al(i, j) - q(i, j) .GE. 0.)
THEN 4415 y4 = al(i, j) - q(i, j)
4418 y4 = -(al(i, j)-q(i, j))
4421 IF (x4 .GT. y4)
THEN 4430 bl(i, j) = -sign(min4, xt)
4431 IF (xt .GE. 0.)
THEN 4438 IF (al(i, j+1) - q(i, j) .GE. 0.)
THEN 4439 y5 = al(i, j+1) - q(i, j)
4442 y5 = -(al(i, j+1)-q(i, j))
4445 IF (x5 .GT. y5)
THEN 4454 br(i, j) = sign(min5, xt)
4458 ELSE IF (jord .EQ. 11)
THEN 4462 IF (xt .GE. 0.)
THEN 4469 IF (al(i, j) - q(i, j) .GE. 0.)
THEN 4470 y6 = al(i, j) - q(i, j)
4473 y6 = -(al(i, j)-q(i, j))
4476 IF (x6 .GT. y6)
THEN 4485 bl(i, j) = -sign(min6, xt)
4486 IF (xt .GE. 0.)
THEN 4493 IF (al(i, j+1) - q(i, j) .GE. 0.)
THEN 4494 y7 = al(i, j+1) - q(i, j)
4497 y7 = -(al(i, j+1)-q(i, j))
4500 IF (x7 .GT. y7)
THEN 4509 br(i, j) = sign(min7, xt)
4516 dq(i, j) = 2.*(q(i, j+1)-q(i, j))
4521 bl(i, j) = al(i, j) - q(i, j)
4522 br(i, j) = al(i, j+1) - q(i, j)
4523 IF (dm(i, j-1) .GE. 0.)
THEN 4528 IF (dm(i, j) .GE. 0.)
THEN 4533 IF (dm(i, j+1) .GE. 0.)
THEN 4538 IF (abs2 + abs5 + abs7 .LT.
near_zero)
THEN 4543 IF (3.*(bl(i, j)+br(i, j)) .GE. 0.)
THEN 4544 abs3 = 3.*(bl(i, j)+br(i, j))
4546 abs3 = -(3.*(bl(i, j)+br(i, j)))
4548 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 4549 abs6 = bl(i, j) - br(i, j)
4551 abs6 = -(bl(i, j)-br(i, j))
4553 IF (abs3 .GT. abs6)
THEN 4555 lac_2 = pmp_2 - 0.75*dq(i, j-2)
4556 IF (0. .LT. pmp_2)
THEN 4557 IF (pmp_2 .LT. lac_2)
THEN 4564 ELSE IF (0. .LT. lac_2)
THEN 4571 IF (0. .GT. pmp_2)
THEN 4572 IF (pmp_2 .GT. lac_2)
THEN 4579 ELSE IF (0. .GT. lac_2)
THEN 4586 IF (br(i, j) .LT. y14)
THEN 4593 IF (x8 .GT. y8)
THEN 4601 lac_1 = pmp_1 + 0.75*dq(i, j+1)
4602 IF (0. .LT. pmp_1)
THEN 4603 IF (pmp_1 .LT. lac_1)
THEN 4610 ELSE IF (0. .LT. lac_1)
THEN 4617 IF (0. .GT. pmp_1)
THEN 4618 IF (pmp_1 .GT. lac_1)
THEN 4625 ELSE IF (0. .GT. lac_1)
THEN 4632 IF (bl(i, j) .LT. y15)
THEN 4639 IF (x9 .GT. y9)
THEN 4654 IF (jord .EQ. 9 .OR. jord .EQ. 13)
THEN 4657 arg1 = ilast - ifirst + 1
4660 CALL pert_ppm(arg1, q(ifirst:ilast, j), bl(ifirst:ilast, j), &
4661 & br(ifirst:ilast, j), 0)
4667 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 4670 bl(i, 0) =
s14*dm(i, -1) +
s11*(q(i, -1)-q(i, 0))
4671 xt = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)*q(i, &
4672 & -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2))*q(i&
4673 & , 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
4674 IF (q(i, 1) .GT. q(i, 2))
THEN 4681 IF (q(i, -1) .GT. q(i, 0))
THEN 4682 IF (q(i, 0) .GT. z2)
THEN 4689 ELSE IF (q(i, -1) .GT. z2)
THEN 4696 IF (xt .LT. y10)
THEN 4703 IF (q(i, 1) .LT. q(i, 2))
THEN 4710 IF (q(i, -1) .LT. q(i, 0))
THEN 4711 IF (q(i, 0) .LT. z3)
THEN 4718 ELSE IF (q(i, -1) .LT. z3)
THEN 4725 IF (xt .GT. y11)
THEN 4733 br(i, 0) = xt - q(i, 0)
4734 bl(i, 1) = xt - q(i, 1)
4735 xt =
s15*q(i, 1) +
s11*q(i, 2) -
s14*dm(i, 2)
4736 br(i, 1) = xt - q(i, 1)
4737 bl(i, 2) = xt - q(i, 2)
4738 br(i, 2) = al(i, 3) - q(i, 2)
4740 arg1 = 3*(ilast-ifirst+1)
4743 CALL pert_ppm(arg1, q(ifirst:ilast, 0:2), bl(ifirst:ilast, 0:2&
4744 & ), br(ifirst:ilast, 0:2), 1)
4749 IF (je + 1 .EQ. npy)
THEN 4751 bl(i, npy-2) = al(i, npy-2) - q(i, npy-2)
4752 xt =
s15*q(i, npy-1) +
s11*q(i, npy-2) +
s14*dm(i, npy-2)
4753 br(i, npy-2) = xt - q(i, npy-2)
4754 bl(i, npy-1) = xt - q(i, npy-1)
4755 xt = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy-1)-dya(&
4756 & i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1))+((2.*&
4757 & dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q(i, npy+&
4758 & 1))/(dya(i, npy)+dya(i, npy+1)))
4759 IF (q(i, npy) .GT. q(i, npy+1))
THEN 4766 IF (q(i, npy-2) .GT. q(i, npy-1))
THEN 4767 IF (q(i, npy-1) .GT. z4)
THEN 4774 ELSE IF (q(i, npy-2) .GT. z4)
THEN 4781 IF (xt .LT. y12)
THEN 4788 IF (q(i, npy) .LT. q(i, npy+1))
THEN 4795 IF (q(i, npy-2) .LT. q(i, npy-1))
THEN 4796 IF (q(i, npy-1) .LT. z5)
THEN 4803 ELSE IF (q(i, npy-2) .LT. z5)
THEN 4810 IF (xt .GT. y13)
THEN 4818 br(i, npy-1) = xt - q(i, npy-1)
4819 bl(i, npy) = xt - q(i, npy)
4820 br(i, npy) =
s11*(q(i, npy+1)-q(i, npy)) -
s14*dm(i, npy+1)
4822 arg1 = 3*(ilast-ifirst+1)
4825 CALL pert_ppm(arg1, q(ifirst:ilast, npy-2:npy), bl(ifirst:&
4826 & ilast, npy-2:npy), br(ifirst:ilast, npy-2:npy), 1)
4836 IF (c(i, j) .GT. 0.)
THEN 4846 DO i=ilast,ifirst,-1
4848 IF (branch .EQ. 0)
THEN 4849 temp2 = bl(i, j) + br(i, j)
4850 temp_ad23 = (c(i, j)+1.)*flux_ad(i, j)
4851 temp_ad24 = c(i, j)*temp_ad23
4852 q_ad(i, j) = q_ad(i, j) + flux_ad(i, j)
4853 c_ad(i, j) = c_ad(i, j) + temp2*temp_ad23 + (bl(i, j)+c(i, j&
4854 & )*temp2)*flux_ad(i, j)
4855 bl_ad(i, j) = bl_ad(i, j) + temp_ad24 + temp_ad23
4856 br_ad(i, j) = br_ad(i, j) + temp_ad24
4859 temp1 = bl(i, j-1) + br(i, j-1)
4860 temp_ad21 = (1.-c(i, j))*flux_ad(i, j)
4861 temp_ad22 = -(c(i, j)*temp_ad21)
4862 q_ad(i, j-1) = q_ad(i, j-1) + flux_ad(i, j)
4863 c_ad(i, j) = c_ad(i, j) - temp1*temp_ad21 - (br(i, j-1)-c(i&
4864 & , j)*temp1)*flux_ad(i, j)
4865 br_ad(i, j-1) = br_ad(i, j-1) + temp_ad22 + temp_ad21
4866 bl_ad(i, j-1) = bl_ad(i, j-1) + temp_ad22
4872 IF (branch .EQ. 0)
THEN 4876 IF (branch .EQ. 1)
THEN 4882 CALL pert_ppm_adm(arg1, q(ifirst:ilast, npy-2:npy), bl(ifirst:&
4883 & ilast, npy-2:npy), bl_ad(ifirst:ilast, npy-2:npy)&
4884 & , br(ifirst:ilast, npy-2:npy), br_ad(ifirst:ilast&
4888 DO i=ilast,ifirst,-1
4889 q_ad(i, npy+1) = q_ad(i, npy+1) +
s11*br_ad(i, npy)
4890 q_ad(i, npy) = q_ad(i, npy) - bl_ad(i, npy) -
s11*br_ad(i, &
4892 dm_ad(i, npy+1) = dm_ad(i, npy+1) -
s14*br_ad(i, npy)
4894 xt_ad = br_ad(i, npy-1) + bl_ad(i, npy)
4896 q_ad(i, npy-1) = q_ad(i, npy-1) - br_ad(i, npy-1)
4897 br_ad(i, npy-1) = 0.0
4899 IF (branch .EQ. 0)
THEN 4906 IF (branch .LT. 2)
THEN 4907 IF (branch .EQ. 0)
THEN 4910 q_ad(i, npy-1) = q_ad(i, npy-1) + y13_ad
4913 ELSE IF (branch .EQ. 2)
THEN 4916 q_ad(i, npy-2) = q_ad(i, npy-2) + y13_ad
4920 IF (branch .EQ. 0)
THEN 4921 q_ad(i, npy+1) = q_ad(i, npy+1) + z5_ad
4923 q_ad(i, npy) = q_ad(i, npy) + z5_ad
4926 IF (branch .EQ. 0)
THEN 4933 IF (branch .LT. 2)
THEN 4934 IF (branch .EQ. 0)
THEN 4937 q_ad(i, npy-1) = q_ad(i, npy-1) + y12_ad
4940 ELSE IF (branch .EQ. 2)
THEN 4943 q_ad(i, npy-2) = q_ad(i, npy-2) + y12_ad
4947 IF (branch .EQ. 0)
THEN 4948 q_ad(i, npy+1) = q_ad(i, npy+1) + z4_ad
4950 q_ad(i, npy) = q_ad(i, npy) + z4_ad
4952 temp_ad19 = 0.5*xt_ad/(dya(i, npy-2)+dya(i, npy-1))
4953 temp_ad20 = 0.5*xt_ad/(dya(i, npy)+dya(i, npy+1))
4954 q_ad(i, npy-1) = q_ad(i, npy-1) + (dya(i, npy-1)*2.+dya(i, &
4956 q_ad(i, npy-2) = q_ad(i, npy-2) - dya(i, npy-1)*temp_ad19
4957 q_ad(i, npy) = q_ad(i, npy) + (dya(i, npy)*2.+dya(i, npy+1))&
4959 q_ad(i, npy+1) = q_ad(i, npy+1) - dya(i, npy)*temp_ad20
4960 xt_ad = br_ad(i, npy-2) + bl_ad(i, npy-1)
4961 q_ad(i, npy-1) = q_ad(i, npy-1) - bl_ad(i, npy-1)
4962 bl_ad(i, npy-1) = 0.0
4963 q_ad(i, npy-2) = q_ad(i, npy-2) - br_ad(i, npy-2)
4964 br_ad(i, npy-2) = 0.0
4965 q_ad(i, npy-1) = q_ad(i, npy-1) +
s15*xt_ad
4966 q_ad(i, npy-2) = q_ad(i, npy-2) +
s11*xt_ad - bl_ad(i, npy-2&
4968 dm_ad(i, npy-2) = dm_ad(i, npy-2) +
s14*xt_ad
4969 al_ad(i, npy-2) = al_ad(i, npy-2) + bl_ad(i, npy-2)
4970 bl_ad(i, npy-2) = 0.0
4974 IF (branch .EQ. 0)
THEN 4975 arg1 = 3*(ilast-ifirst+1)
4978 CALL pert_ppm_adm(arg1, q(ifirst:ilast, 0:2), bl(ifirst:ilast&
4979 & , 0:2), bl_ad(ifirst:ilast, 0:2), br(ifirst:ilast&
4980 & , 0:2), br_ad(ifirst:ilast, 0:2), 1)
4981 DO i=ilast,ifirst,-1
4982 al_ad(i, 3) = al_ad(i, 3) + br_ad(i, 2)
4983 q_ad(i, 2) = q_ad(i, 2) - bl_ad(i, 2) - br_ad(i, 2)
4985 xt_ad = br_ad(i, 1) + bl_ad(i, 2)
4987 q_ad(i, 1) = q_ad(i, 1) +
s15*xt_ad - br_ad(i, 1)
4989 q_ad(i, 2) = q_ad(i, 2) +
s11*xt_ad
4990 dm_ad(i, 2) = dm_ad(i, 2) -
s14*xt_ad
4991 xt_ad = br_ad(i, 0) + bl_ad(i, 1)
4992 q_ad(i, 1) = q_ad(i, 1) - bl_ad(i, 1)
4994 q_ad(i, 0) = q_ad(i, 0) - br_ad(i, 0)
4997 IF (branch .EQ. 0)
THEN 5004 IF (branch .LT. 2)
THEN 5005 IF (branch .EQ. 0)
THEN 5008 q_ad(i, 0) = q_ad(i, 0) + y11_ad
5011 ELSE IF (branch .EQ. 2)
THEN 5014 q_ad(i, -1) = q_ad(i, -1) + y11_ad
5018 IF (branch .EQ. 0)
THEN 5019 q_ad(i, 2) = q_ad(i, 2) + z3_ad
5021 q_ad(i, 1) = q_ad(i, 1) + z3_ad
5024 IF (branch .EQ. 0)
THEN 5031 IF (branch .LT. 2)
THEN 5032 IF (branch .EQ. 0)
THEN 5035 q_ad(i, 0) = q_ad(i, 0) + y10_ad
5038 ELSE IF (branch .EQ. 2)
THEN 5041 q_ad(i, -1) = q_ad(i, -1) + y10_ad
5045 IF (branch .EQ. 0)
THEN 5046 q_ad(i, 2) = q_ad(i, 2) + z2_ad
5048 q_ad(i, 1) = q_ad(i, 1) + z2_ad
5050 temp_ad17 = 0.5*xt_ad/(dya(i, -1)+dya(i, 0))
5051 temp_ad18 = 0.5*xt_ad/(dya(i, 1)+dya(i, 2))
5052 q_ad(i, 0) = q_ad(i, 0) + (dya(i, 0)*2.+dya(i, -1))*&
5054 q_ad(i, -1) = q_ad(i, -1) - dya(i, 0)*temp_ad17
5055 q_ad(i, 1) = q_ad(i, 1) + (dya(i, 1)*2.+dya(i, 2))*temp_ad18
5056 q_ad(i, 2) = q_ad(i, 2) - dya(i, 1)*temp_ad18
5057 dm_ad(i, -1) = dm_ad(i, -1) +
s14*bl_ad(i, 0)
5058 q_ad(i, -1) = q_ad(i, -1) +
s11*bl_ad(i, 0)
5059 q_ad(i, 0) = q_ad(i, 0) -
s11*bl_ad(i, 0)
5065 IF (branch .EQ. 0)
THEN 5067 arg1 = ilast - ifirst + 1
5070 CALL pert_ppm_adm(arg1, q(ifirst:ilast, j), bl(ifirst:ilast, j&
5071 & ), bl_ad(ifirst:ilast, j), br(ifirst:ilast, j), &
5072 & br_ad(ifirst:ilast, j), 0)
5076 IF (branch .EQ. 0)
THEN 5078 DO i=ilast,ifirst,-1
5080 min5_ad = sign(1.d0, min5*xt)*br_ad(i, j)
5083 IF (branch .EQ. 0)
THEN 5093 IF (branch .EQ. 0)
THEN 5094 al_ad(i, j+1) = al_ad(i, j+1) + y5_ad
5095 q_ad(i, j) = q_ad(i, j) - y5_ad
5097 q_ad(i, j) = q_ad(i, j) + y5_ad
5098 al_ad(i, j+1) = al_ad(i, j+1) - y5_ad
5101 IF (branch .EQ. 0)
THEN 5106 min4_ad = -(sign(1.d0, min4*xt)*bl_ad(i, j))
5109 IF (branch .EQ. 0)
THEN 5119 IF (branch .EQ. 0)
THEN 5120 al_ad(i, j) = al_ad(i, j) + y4_ad
5121 q_ad(i, j) = q_ad(i, j) - y4_ad
5123 q_ad(i, j) = q_ad(i, j) + y4_ad
5124 al_ad(i, j) = al_ad(i, j) - y4_ad
5127 IF (branch .EQ. 0)
THEN 5128 xt_ad = xt_ad + x4_ad
5130 xt_ad = xt_ad - x4_ad
5132 dm_ad(i, j) = dm_ad(i, j) + 2.*xt_ad
5135 ELSE IF (branch .EQ. 1)
THEN 5137 DO i=ilast,ifirst,-1
5139 min7_ad = sign(1.d0, min7*xt)*br_ad(i, j)
5142 IF (branch .EQ. 0)
THEN 5152 IF (branch .EQ. 0)
THEN 5153 al_ad(i, j+1) = al_ad(i, j+1) + y7_ad
5154 q_ad(i, j) = q_ad(i, j) - y7_ad
5156 q_ad(i, j) = q_ad(i, j) + y7_ad
5157 al_ad(i, j+1) = al_ad(i, j+1) - y7_ad
5160 IF (branch .EQ. 0)
THEN 5165 min6_ad = -(sign(1.d0, min6*xt)*bl_ad(i, j))
5168 IF (branch .EQ. 0)
THEN 5178 IF (branch .EQ. 0)
THEN 5179 al_ad(i, j) = al_ad(i, j) + y6_ad
5180 q_ad(i, j) = q_ad(i, j) - y6_ad
5182 q_ad(i, j) = q_ad(i, j) + y6_ad
5183 al_ad(i, j) = al_ad(i, j) - y6_ad
5186 IF (branch .EQ. 0)
THEN 5187 xt_ad = xt_ad + x6_ad
5189 xt_ad = xt_ad - x6_ad
5191 dm_ad(i, j) = dm_ad(i, j) +
ppm_fac*xt_ad
5197 DO i=ilast,ifirst,-1
5199 IF (branch .LT. 2)
THEN 5200 IF (branch .EQ. 0)
THEN 5207 ELSE IF (branch .EQ. 2)
THEN 5217 IF (branch .EQ. 0)
THEN 5220 bl_ad(i, j) = bl_ad(i, j) + y9_ad
5224 IF (branch .LT. 2)
THEN 5225 IF (branch .EQ. 0)
THEN 5233 IF (branch .EQ. 2)
THEN 5241 IF (branch .LT. 2)
THEN 5242 IF (branch .EQ. 0)
THEN 5243 lac_1_ad = lac_1_ad + x9_ad
5245 pmp_1_ad = pmp_1_ad + x9_ad
5247 ELSE IF (branch .EQ. 2)
THEN 5248 lac_1_ad = lac_1_ad + x9_ad
5250 pmp_1_ad = pmp_1_ad + lac_1_ad
5251 dq_ad(i, j+1) = dq_ad(i, j+1) + 0.75*lac_1_ad
5252 dq_ad(i, j) = dq_ad(i, j) - pmp_1_ad
5254 IF (branch .EQ. 0)
THEN 5264 IF (branch .EQ. 0)
THEN 5267 br_ad(i, j) = br_ad(i, j) + y8_ad
5271 IF (branch .LT. 2)
THEN 5272 IF (branch .EQ. 0)
THEN 5280 IF (branch .EQ. 2)
THEN 5288 IF (branch .LT. 2)
THEN 5289 IF (branch .EQ. 0)
THEN 5290 lac_2_ad = lac_2_ad + x8_ad
5292 pmp_2_ad = pmp_2_ad + x8_ad
5294 ELSE IF (branch .EQ. 2)
THEN 5295 lac_2_ad = lac_2_ad + x8_ad
5297 pmp_2_ad = pmp_2_ad + lac_2_ad
5298 dq_ad(i, j-2) = dq_ad(i, j-2) - 0.75*lac_2_ad
5299 dq_ad(i, j-1) = dq_ad(i, j-1) + pmp_2_ad
5300 110 al_ad(i, j+1) = al_ad(i, j+1) + br_ad(i, j)
5301 q_ad(i, j) = q_ad(i, j) - bl_ad(i, j) - br_ad(i, j)
5303 al_ad(i, j) = al_ad(i, j) + bl_ad(i, j)
5308 DO i=ilast,ifirst,-1
5309 q_ad(i, j+1) = q_ad(i, j+1) + 2.*dq_ad(i, j)
5310 q_ad(i, j) = q_ad(i, j) - 2.*dq_ad(i, j)
5316 DO i=ilast,ifirst,-1
5317 q_ad(i, j-1) = q_ad(i, j-1) + 0.5*al_ad(i, j)
5318 q_ad(i, j) = q_ad(i, j) + 0.5*al_ad(i, j)
5319 dm_ad(i, j-1) = dm_ad(i, j-1) +
r3*al_ad(i, j)
5320 dm_ad(i, j) = dm_ad(i, j) -
r3*al_ad(i, j)
5325 DO i=ilast,ifirst,-1
5326 xt = 0.25*(q(i, j+1)-q(i, j-1))
5327 min3_ad = sign(1.d0, min3*xt)*dm_ad(i, j)
5330 IF (branch .LT. 2)
THEN 5331 IF (branch .EQ. 0)
THEN 5342 IF (branch .EQ. 2)
THEN 5353 q_ad(i, j) = q_ad(i, j) + z1_ad
5356 IF (branch .LT. 2)
THEN 5357 IF (branch .EQ. 0)
THEN 5358 q_ad(i, j+1) = q_ad(i, j+1) + min8_ad
5360 q_ad(i, j) = q_ad(i, j) + min8_ad
5362 ELSE IF (branch .EQ. 2)
THEN 5363 q_ad(i, j+1) = q_ad(i, j+1) + min8_ad
5365 q_ad(i, j-1) = q_ad(i, j-1) + min8_ad
5368 q_ad(i, j) = q_ad(i, j) - y3_ad
5370 IF (branch .LT. 2)
THEN 5371 IF (branch .EQ. 0)
THEN 5372 q_ad(i, j+1) = q_ad(i, j+1) + max1_ad
5374 q_ad(i, j) = q_ad(i, j) + max1_ad
5376 ELSE IF (branch .EQ. 2)
THEN 5377 q_ad(i, j+1) = q_ad(i, j+1) + max1_ad
5379 q_ad(i, j-1) = q_ad(i, j-1) + max1_ad
5382 IF (branch .EQ. 0)
THEN 5387 q_ad(i, j+1) = q_ad(i, j+1) + 0.25*xt_ad
5388 q_ad(i, j-1) = q_ad(i, j-1) - 0.25*xt_ad
5394 SUBROUTINE yppm(flux, q, c, jord, ifirst, ilast, isd, ied, js, je, jsd&
5395 & , jed, npx, npy, dya, nested, grid_type)
5398 INTEGER,
INTENT(IN) :: ifirst, ilast
5399 INTEGER,
INTENT(IN) :: isd, ied, js, je, jsd, jed
5400 INTEGER,
INTENT(IN) :: jord
5401 INTEGER,
INTENT(IN) :: npx, npy
5402 REAL,
INTENT(IN) :: q(ifirst:ilast, jsd:jed)
5404 REAL,
INTENT(IN) :: c(isd:ied, js:je+1)
5406 REAL,
INTENT(OUT) :: flux(ifirst:ilast, js:je+1)
5407 REAL,
INTENT(IN) :: dya(isd:ied, jsd:jed)
5408 LOGICAL,
INTENT(IN) :: nested
5409 INTEGER,
INTENT(IN) :: grid_type
5411 REAL :: dm(ifirst:ilast, js-2:je+2)
5412 REAL :: al(ifirst:ilast, js-1:je+2)
5413 REAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: bl, br, b0
5414 REAL :: dq(ifirst:ilast, js-3:je+2)
5415 REAL,
DIMENSION(ifirst:ilast) :: fx0, fx1
5416 LOGICAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: smt5, smt6
5417 REAL :: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
5418 INTEGER :: i, j, js1, je3, je1
5470 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 5471 IF (3 .LT. js - 1)
THEN 5476 IF (npy - 2 .GT. je + 2)
THEN 5481 IF (npy - 3 .GT. je + 1)
THEN 5492 IF (jord .LT. 8 .OR. jord .EQ. 333)
THEN 5495 al(i, j) =
p1*(q(i, j-1)+q(i, j)) +
p2*(q(i, j-2)+q(i, j+1))
5498 IF (jord .EQ. 7)
THEN 5501 IF (al(i, j) .LT. 0.) al(i, j) = 0.5*(q(i, j)+q(i, j+1))
5505 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 5508 al(i, 0) =
c1*q(i, -2) +
c2*q(i, -1) +
c3*q(i, 0)
5509 al(i, 1) = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)&
5510 & *q(i, -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2)&
5511 & )*q(i, 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
5512 al(i, 2) =
c3*q(i, 1) +
c2*q(i, 2) +
c1*q(i, 3)
5514 IF (jord .EQ. 7)
THEN 5516 IF (0. .LT. al(i, 0))
THEN 5521 IF (0. .LT. al(i, 1))
THEN 5526 IF (0. .LT. al(i, 2))
THEN 5534 IF (je + 1 .EQ. npy)
THEN 5536 al(i, npy-1) =
c1*q(i, npy-3) +
c2*q(i, npy-2) +
c3*q(i, npy&
5538 al(i, npy) = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy&
5539 & -1)-dya(i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1&
5540 & ))+((2.*dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q&
5541 & (i, npy+1))/(dya(i, npy)+dya(i, npy+1)))
5542 al(i, npy+1) =
c3*q(i, npy) +
c2*q(i, npy+1) +
c1*q(i, npy+2&
5545 IF (jord .EQ. 7)
THEN 5547 IF (0. .LT. al(i, npy-1))
THEN 5548 al(i, npy-1) = al(i, npy-1)
5552 IF (0. .LT. al(i, npy))
THEN 5553 al(i, npy) = al(i, npy)
5557 IF (0. .LT. al(i, npy+1))
THEN 5558 al(i, npy+1) = al(i, npy+1)
5566 IF (jord .EQ. 1)
THEN 5569 IF (c(i, j) .GT. 0.)
THEN 5570 flux(i, j) = q(i, j-1)
5572 flux(i, j) = q(i, j)
5576 ELSE IF (jord .EQ. 2)
THEN 5583 IF (xt .GT. 0.)
THEN 5585 flux(i, j) = qtmp + (1.-xt)*(al(i, j)-qtmp-xt*(al(i, j-1)+&
5586 & al(i, j)-(qtmp+qtmp)))
5589 flux(i, j) = qtmp + (1.+xt)*(al(i, j)-qtmp+xt*(al(i, j)+al&
5590 & (i, j+1)-(qtmp+qtmp)))
5594 ELSE IF (jord .EQ. 333)
THEN 5600 IF (xt .GT. 0.)
THEN 5601 flux(i, j) = (2.0*q(i, j)+5.0*q(i, j-1)-q(i, j-2))/6.0 - &
5602 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j)-2.0*q(i&
5605 flux(i, j) = (2.0*q(i, j-1)+5.0*q(i, j)-q(i, j+1))/6.0 - &
5606 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j+1)-2.0*q(&
5611 ELSE IF (jord .EQ. 3)
THEN 5614 bl(i, j) = al(i, j) - q(i, j)
5615 br(i, j) = al(i, j+1) - q(i, j)
5616 b0(i, j) = bl(i, j) + br(i, j)
5617 IF (b0(i, j) .GE. 0.)
THEN 5622 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 5623 xt = bl(i, j) - br(i, j)
5625 xt = -(bl(i, j)-br(i, j))
5627 smt5(i, j) = x0 .LT. xt
5628 smt6(i, j) = 3.*x0 .LT. xt
5637 IF (xt .GT. 0.)
THEN 5639 IF (smt6(i, j-1) .OR. smt5(i, j))
THEN 5640 fx1(i) = br(i, j-1) - xt*b0(i, j-1)
5641 ELSE IF (smt5(i, j-1))
THEN 5642 IF (bl(i, j-1) .GE. 0.)
THEN 5647 IF (br(i, j-1) .GE. 0.)
THEN 5652 IF (x1 .GT. y1)
THEN 5658 fx1(i) = sign(min1, br(i, j-1))
5662 IF (smt6(i, j) .OR. smt5(i, j-1))
THEN 5663 fx1(i) = bl(i, j) + xt*b0(i, j)
5664 ELSE IF (smt5(i, j))
THEN 5665 IF (bl(i, j) .GE. 0.)
THEN 5670 IF (br(i, j) .GE. 0.)
THEN 5675 IF (x2 .GT. y2)
THEN 5680 fx1(i) = sign(min2, bl(i, j))
5683 IF (xt .GE. 0.)
THEN 5688 flux(i, j) = fx0(i) + (1.-abs0)*fx1(i)
5691 ELSE IF (jord .EQ. 4)
THEN 5694 bl(i, j) = al(i, j) - q(i, j)
5695 br(i, j) = al(i, j+1) - q(i, j)
5696 b0(i, j) = bl(i, j) + br(i, j)
5697 IF (b0(i, j) .GE. 0.)
THEN 5702 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 5703 xt = bl(i, j) - br(i, j)
5705 xt = -(bl(i, j)-br(i, j))
5707 smt5(i, j) = x0 .LT. xt
5708 smt6(i, j) = 3.*x0 .LT. xt
5717 IF (c(i, j) .GT. 0.)
THEN 5719 IF (smt6(i, j-1) .OR. smt5(i, j)) fx1(i) = (1.-c(i, j))*(&
5720 & br(i, j-1)-c(i, j)*b0(i, j-1))
5723 IF (smt6(i, j) .OR. smt5(i, j-1)) fx1(i) = (1.+c(i, j))*(&
5724 & bl(i, j)+c(i, j)*b0(i, j))
5726 flux(i, j) = fx0(i) + fx1(i)
5731 IF (jord .EQ. 5)
THEN 5734 bl(i, j) = al(i, j) - q(i, j)
5735 br(i, j) = al(i, j+1) - q(i, j)
5736 b0(i, j) = bl(i, j) + br(i, j)
5737 smt5(i, j) = bl(i, j)*br(i, j) .LT. 0.
5743 bl(i, j) = al(i, j) - q(i, j)
5744 br(i, j) = al(i, j+1) - q(i, j)
5745 b0(i, j) = bl(i, j) + br(i, j)
5746 IF (3.*b0(i, j) .GE. 0.)
THEN 5749 abs1 = -(3.*b0(i, j))
5751 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 5752 abs4 = bl(i, j) - br(i, j)
5754 abs4 = -(bl(i, j)-br(i, j))
5756 smt5(i, j) = abs1 .LT. abs4
5763 IF (c(i, j) .GT. 0.)
THEN 5764 fx1(i) = (1.-c(i, j))*(br(i, j-1)-c(i, j)*b0(i, j-1))
5765 flux(i, j) = q(i, j-1)
5767 fx1(i) = (1.+c(i, j))*(bl(i, j)+c(i, j)*b0(i, j))
5768 flux(i, j) = q(i, j)
5770 IF (smt5(i, j-1) .OR. smt5(i, j)) flux(i, j) = flux(i, j) + &
5782 xt = 0.25*(q(i, j+1)-q(i, j-1))
5783 IF (xt .GE. 0.)
THEN 5788 IF (q(i, j-1) .LT. q(i, j))
THEN 5789 IF (q(i, j) .LT. q(i, j+1))
THEN 5794 ELSE IF (q(i, j-1) .LT. q(i, j+1))
THEN 5800 IF (q(i, j-1) .GT. q(i, j))
THEN 5801 IF (q(i, j) .GT. q(i, j+1))
THEN 5806 ELSE IF (q(i, j-1) .GT. q(i, j+1))
THEN 5812 IF (x3 .GT. y3)
THEN 5813 IF (y3 .GT. z1)
THEN 5818 ELSE IF (x3 .GT. z1)
THEN 5823 dm(i, j) = sign(min3, xt)
5828 al(i, j) = 0.5*(q(i, j-1)+q(i, j)) +
r3*(dm(i, j-1)-dm(i, j))
5831 IF (jord .EQ. 8)
THEN 5835 IF (xt .GE. 0.)
THEN 5840 IF (al(i, j) - q(i, j) .GE. 0.)
THEN 5841 y4 = al(i, j) - q(i, j)
5843 y4 = -(al(i, j)-q(i, j))
5845 IF (x4 .GT. y4)
THEN 5850 bl(i, j) = -sign(min4, xt)
5851 IF (xt .GE. 0.)
THEN 5856 IF (al(i, j+1) - q(i, j) .GE. 0.)
THEN 5857 y5 = al(i, j+1) - q(i, j)
5859 y5 = -(al(i, j+1)-q(i, j))
5861 IF (x5 .GT. y5)
THEN 5866 br(i, j) = sign(min5, xt)
5869 ELSE IF (jord .EQ. 11)
THEN 5873 IF (xt .GE. 0.)
THEN 5878 IF (al(i, j) - q(i, j) .GE. 0.)
THEN 5879 y6 = al(i, j) - q(i, j)
5881 y6 = -(al(i, j)-q(i, j))
5883 IF (x6 .GT. y6)
THEN 5888 bl(i, j) = -sign(min6, xt)
5889 IF (xt .GE. 0.)
THEN 5894 IF (al(i, j+1) - q(i, j) .GE. 0.)
THEN 5895 y7 = al(i, j+1) - q(i, j)
5897 y7 = -(al(i, j+1)-q(i, j))
5899 IF (x7 .GT. y7)
THEN 5904 br(i, j) = sign(min7, xt)
5910 dq(i, j) = 2.*(q(i, j+1)-q(i, j))
5915 bl(i, j) = al(i, j) - q(i, j)
5916 br(i, j) = al(i, j+1) - q(i, j)
5917 IF (dm(i, j-1) .GE. 0.)
THEN 5922 IF (dm(i, j) .GE. 0.)
THEN 5927 IF (dm(i, j+1) .GE. 0.)
THEN 5932 IF (abs2 + abs5 + abs7 .LT.
near_zero)
THEN 5936 IF (3.*(bl(i, j)+br(i, j)) .GE. 0.)
THEN 5937 abs3 = 3.*(bl(i, j)+br(i, j))
5939 abs3 = -(3.*(bl(i, j)+br(i, j)))
5941 IF (bl(i, j) - br(i, j) .GE. 0.)
THEN 5942 abs6 = bl(i, j) - br(i, j)
5944 abs6 = -(bl(i, j)-br(i, j))
5946 IF (abs3 .GT. abs6)
THEN 5948 lac_2 = pmp_2 - 0.75*dq(i, j-2)
5949 IF (0. .LT. pmp_2)
THEN 5950 IF (pmp_2 .LT. lac_2)
THEN 5955 ELSE IF (0. .LT. lac_2)
THEN 5960 IF (0. .GT. pmp_2)
THEN 5961 IF (pmp_2 .GT. lac_2)
THEN 5966 ELSE IF (0. .GT. lac_2)
THEN 5971 IF (br(i, j) .LT. y14)
THEN 5976 IF (x8 .GT. y8)
THEN 5982 lac_1 = pmp_1 + 0.75*dq(i, j+1)
5983 IF (0. .LT. pmp_1)
THEN 5984 IF (pmp_1 .LT. lac_1)
THEN 5989 ELSE IF (0. .LT. lac_1)
THEN 5994 IF (0. .GT. pmp_1)
THEN 5995 IF (pmp_1 .GT. lac_1)
THEN 6000 ELSE IF (0. .GT. lac_1)
THEN 6005 IF (bl(i, j) .LT. y15)
THEN 6010 IF (x9 .GT. y9)
THEN 6020 IF (jord .EQ. 9 .OR. jord .EQ. 13)
THEN 6023 arg1 = ilast - ifirst + 1
6024 CALL pert_ppm(arg1, q(ifirst:ilast, j), bl(ifirst:ilast, j), &
6025 & br(ifirst:ilast, j), 0)
6028 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 6031 bl(i, 0) =
s14*dm(i, -1) +
s11*(q(i, -1)-q(i, 0))
6032 xt = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)*q(i, &
6033 & -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2))*q(i&
6034 & , 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
6035 IF (q(i, 1) .GT. q(i, 2))
THEN 6040 IF (q(i, -1) .GT. q(i, 0))
THEN 6041 IF (q(i, 0) .GT. z2)
THEN 6046 ELSE IF (q(i, -1) .GT. z2)
THEN 6051 IF (xt .LT. y10)
THEN 6056 IF (q(i, 1) .LT. q(i, 2))
THEN 6061 IF (q(i, -1) .LT. q(i, 0))
THEN 6062 IF (q(i, 0) .LT. z3)
THEN 6067 ELSE IF (q(i, -1) .LT. z3)
THEN 6072 IF (xt .GT. y11)
THEN 6078 br(i, 0) = xt - q(i, 0)
6079 bl(i, 1) = xt - q(i, 1)
6080 xt =
s15*q(i, 1) +
s11*q(i, 2) -
s14*dm(i, 2)
6081 br(i, 1) = xt - q(i, 1)
6082 bl(i, 2) = xt - q(i, 2)
6083 br(i, 2) = al(i, 3) - q(i, 2)
6085 arg1 = 3*(ilast-ifirst+1)
6086 CALL pert_ppm(arg1, q(ifirst:ilast, 0:2), bl(ifirst:ilast, 0:2&
6087 & ), br(ifirst:ilast, 0:2), 1)
6089 IF (je + 1 .EQ. npy)
THEN 6091 bl(i, npy-2) = al(i, npy-2) - q(i, npy-2)
6092 xt =
s15*q(i, npy-1) +
s11*q(i, npy-2) +
s14*dm(i, npy-2)
6093 br(i, npy-2) = xt - q(i, npy-2)
6094 bl(i, npy-1) = xt - q(i, npy-1)
6095 xt = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy-1)-dya(&
6096 & i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1))+((2.*&
6097 & dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q(i, npy+&
6098 & 1))/(dya(i, npy)+dya(i, npy+1)))
6099 IF (q(i, npy) .GT. q(i, npy+1))
THEN 6104 IF (q(i, npy-2) .GT. q(i, npy-1))
THEN 6105 IF (q(i, npy-1) .GT. z4)
THEN 6110 ELSE IF (q(i, npy-2) .GT. z4)
THEN 6115 IF (xt .LT. y12)
THEN 6120 IF (q(i, npy) .LT. q(i, npy+1))
THEN 6125 IF (q(i, npy-2) .LT. q(i, npy-1))
THEN 6126 IF (q(i, npy-1) .LT. z5)
THEN 6131 ELSE IF (q(i, npy-2) .LT. z5)
THEN 6136 IF (xt .GT. y13)
THEN 6142 br(i, npy-1) = xt - q(i, npy-1)
6143 bl(i, npy) = xt - q(i, npy)
6144 br(i, npy) =
s11*(q(i, npy+1)-q(i, npy)) -
s14*dm(i, npy+1)
6146 arg1 = 3*(ilast-ifirst+1)
6147 CALL pert_ppm(arg1, q(ifirst:ilast, npy-2:npy), bl(ifirst:&
6148 & ilast, npy-2:npy), br(ifirst:ilast, npy-2:npy), 1)
6153 IF (c(i, j) .GT. 0.)
THEN 6154 flux(i, j) = q(i, j-1) + (1.-c(i, j))*(br(i, j-1)-c(i, j)*(&
6155 & bl(i, j-1)+br(i, j-1)))
6157 flux(i, j) = q(i, j) + (1.+c(i, j))*(bl(i, j)+c(i, j)*(bl(i&
6164 SUBROUTINE mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
6165 & kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
6169 INTEGER,
INTENT(IN) :: im, jm, km, nq
6170 INTEGER,
INTENT(IN) :: ifirst, ilast
6171 INTEGER,
INTENT(IN) :: jfirst, jlast
6172 INTEGER,
INTENT(IN) :: kfirst, klast
6174 INTEGER,
INTENT(IN) :: ng_e
6176 INTEGER,
INTENT(IN) :: ng_w
6178 INTEGER,
INTENT(IN) :: ng_s
6180 INTEGER,
INTENT(IN) :: ng_n
6181 REAL,
INTENT(INOUT) :: q_ghst(ifirst-ng_w:ilast+ng_e, jfirst-ng_s:&
6182 & jlast+ng_n, kfirst:klast, nq)
6183 REAL,
OPTIONAL,
INTENT(IN) :: q(ifirst:ilast, jfirst:jlast, kfirst:&
6196 INTEGER :: i, j, k, n
6198 IF (
PRESENT(q)) q_ghst(ifirst:ilast, jfirst:jlast, kfirst:klast, 1:&
6199 & nq) = q(ifirst:ilast, jfirst:jlast, kfirst:klast, 1:nq)
6203 DO j=jfirst-ng_s,jlast+ng_n
6205 q_ghst(ifirst-i, j, k, n) = q_ghst(ilast-i+1, j, k, n)
6208 q_ghst(ilast+i, j, k, n) = q_ghst(ifirst+i-1, j, k, n)
6234 SUBROUTINE pert_ppm_adm(im, a0, al, al_ad, ar, ar_ad, iv)
6236 INTEGER,
INTENT(IN) :: im
6237 INTEGER,
INTENT(IN) :: iv
6238 REAL,
INTENT(IN) :: a0(im)
6239 REAL,
INTENT(INOUT) :: al(im), ar(im)
6240 REAL,
INTENT(INOUT) :: al_ad(im), ar_ad(im)
6242 REAL :: a4, da1, da2, a6da, fmin
6244 REAL,
PARAMETER :: r12=1./12.
6254 IF (a0(i) .LE. 0.)
THEN 6257 a4 = -(3.*(ar(i)+al(i)))
6259 IF (da1 .GE. 0.)
THEN 6264 IF (abs0 .LT. -a4)
THEN 6265 fmin = a0(i) + 0.25/a4*da1**2 + a4*r12
6266 IF (fmin .LT. 0.)
THEN 6267 IF (ar(i) .GT. 0. .AND. al(i) .GT. 0.)
THEN 6269 ELSE IF (da1 .GT. 0.)
THEN 6284 IF (branch .LT. 3)
THEN 6285 IF (branch .NE. 0)
THEN 6286 IF (branch .NE. 1)
THEN 6287 ar_ad(i) = ar_ad(i) - 2.*al_ad(i)
6291 ELSE IF (branch .EQ. 3)
THEN 6292 al_ad(i) = al_ad(i) - 2.*ar_ad(i)
6294 ELSE IF (branch .EQ. 4)
THEN 6305 IF (al(i)*ar(i) .LT. 0.)
THEN 6308 a6da = 3.*(al(i)+ar(i))*da1
6310 IF (a6da .LT. -da2)
THEN 6312 ELSE IF (a6da .GT. da2)
THEN 6323 IF (branch .LT. 2)
THEN 6324 IF (branch .EQ. 0)
THEN 6328 ELSE IF (branch .EQ. 2)
THEN 6329 ar_ad(i) = ar_ad(i) - 2.*al_ad(i)
6332 al_ad(i) = al_ad(i) - 2.*ar_ad(i)
6338 SUBROUTINE pert_ppm(im, a0, al, ar, iv)
6340 INTEGER,
INTENT(IN) :: im
6341 INTEGER,
INTENT(IN) :: iv
6342 REAL,
INTENT(IN) :: a0(im)
6343 REAL,
INTENT(INOUT) :: al(im), ar(im)
6345 REAL :: a4, da1, da2, a6da, fmin
6347 REAL,
PARAMETER :: r12=1./12.
6356 IF (a0(i) .LE. 0.)
THEN 6360 a4 = -(3.*(ar(i)+al(i)))
6362 IF (da1 .GE. 0.)
THEN 6367 IF (abs0 .LT. -a4)
THEN 6368 fmin = a0(i) + 0.25/a4*da1**2 + a4*r12
6369 IF (fmin .LT. 0.)
THEN 6370 IF (ar(i) .GT. 0. .AND. al(i) .GT. 0.)
THEN 6373 ELSE IF (da1 .GT. 0.)
THEN 6385 IF (al(i)*ar(i) .LT. 0.)
THEN 6388 a6da = 3.*(al(i)+ar(i))*da1
6390 IF (a6da .LT. -da2)
THEN 6392 ELSE IF (a6da .GT. da2)
THEN 6423 SUBROUTINE deln_flux_adm(nord, is, ie, js, je, npx, npy, damp, q, q_ad&
6424 & , fx, fx_ad, fy, fy_ad, gridstruct, bd, mass, mass_ad)
6433 TYPE(FV_GRID_BOUNDS_TYPE),
INTENT(IN) :: bd
6435 INTEGER,
INTENT(IN) :: nord
6436 INTEGER,
INTENT(IN) :: is, ie, js, je, npx, npy
6437 REAL,
INTENT(IN) :: damp
6439 REAL,
INTENT(IN) :: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
6440 REAL :: q_ad(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
6441 TYPE(FV_GRID_TYPE),
INTENT(IN),
TARGET :: gridstruct
6443 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
6444 REAL,
OPTIONAL :: mass_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
6446 REAL,
INTENT(INOUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je), fy(bd%is:bd%&
6447 & ie, bd%js:bd%je+1)
6448 REAL,
INTENT(INOUT) :: fx_ad(bd%is:bd%ie+1, bd%js:bd%je), fy_ad(bd%&
6449 & is:bd%ie, bd%js:bd%je+1)
6451 REAL :: fx2(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2(bd%isd:bd%ied, bd%&
6453 REAL :: fx2_ad(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2_ad(bd%isd:bd%ied&
6454 & , bd%jsd:bd%jed+1)
6455 REAL :: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
6456 REAL :: d2_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
6458 INTEGER :: i, j, n, nt, i1, i2, j1, j2
6484 IF (.NOT.
PRESENT(mass))
THEN 6487 d2(i, j) = damp*q(i, j)
6499 IF (nord .GT. 0)
THEN 6500 CALL copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, &
6501 & gridstruct%sw_corner, gridstruct%se_corner, gridstruct&
6502 & %nw_corner, gridstruct%ne_corner)
6507 DO j=js-nord,je+nord
6508 DO i=is-nord,ie+nord+1
6509 fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i-1, j)-d2(i, j))
6512 IF (nord .GT. 0)
THEN 6513 CALL copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, &
6514 & gridstruct%sw_corner, gridstruct%se_corner, gridstruct&
6515 & %nw_corner, gridstruct%ne_corner)
6520 DO j=js-nord,je+nord+1
6521 DO i=is-nord,ie+nord
6522 fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j-1)-d2(i, j))
6525 IF (nord .GT. 0)
THEN 6531 ad_from0 = js - nt - 1
6532 DO j=ad_from0,je+nt+1
6533 ad_from = is - nt - 1
6534 DO i=ad_from,ie+nt+1
6535 d2(i, j) = (fx2(i, j)-fx2(i+1, j)+fy2(i, j)-fy2(i, j+1))*&
6536 & gridstruct%rarea(i, j)
6543 CALL copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, &
6544 & gridstruct%sw_corner, gridstruct%se_corner, &
6545 & gridstruct%nw_corner, gridstruct%ne_corner)
6549 DO i=ad_from1,ie+nt+1
6550 fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i, j)-d2(i-1, j))
6557 CALL copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, &
6558 & gridstruct%sw_corner, gridstruct%se_corner, &
6559 & gridstruct%nw_corner, gridstruct%ne_corner)
6561 DO j=ad_from4,je+nt+1
6564 fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j)-d2(i, j-1))
6579 IF (
PRESENT(mass))
THEN 6585 temp_ad5 = damp2*fy2(i, j)*fy_ad(i, j)
6586 mass_ad(i, j-1) = mass_ad(i, j-1) + temp_ad5
6587 mass_ad(i, j) = mass_ad(i, j) + temp_ad5
6588 fy2_ad(i, j) = fy2_ad(i, j) + damp2*(mass(i, j-1)+mass(i, j))*&
6595 temp_ad4 = damp2*fx2(i, j)*fx_ad(i, j)
6596 mass_ad(i-1, j) = mass_ad(i-1, j) + temp_ad4
6597 mass_ad(i, j) = mass_ad(i, j) + temp_ad4
6598 fx2_ad(i, j) = fx2_ad(i, j) + damp2*(mass(i-1, j)+mass(i, j))*&
6606 fy2_ad(i, j) = fy2_ad(i, j) + fy_ad(i, j)
6612 fx2_ad(i, j) = fx2_ad(i, j) + fx_ad(i, j)
6617 IF (branch .EQ. 0)
THEN 6622 DO j=ad_to4,ad_from4,-1
6625 DO i=ad_to3,ad_from3,-1
6626 temp_ad3 = gridstruct%del6_u(i, j)*fy2_ad(i, j)
6627 d2_ad(i, j) = d2_ad(i, j) + temp_ad3
6628 d2_ad(i, j-1) = d2_ad(i, j-1) - temp_ad3
6633 & , bd, gridstruct%sw_corner, gridstruct%se_corner&
6634 & , gridstruct%nw_corner, gridstruct%ne_corner)
6637 DO j=ad_to2,ad_from2,-1
6640 DO i=ad_to1,ad_from1,-1
6641 temp_ad2 = gridstruct%del6_v(i, j)*fx2_ad(i, j)
6642 d2_ad(i, j) = d2_ad(i, j) + temp_ad2
6643 d2_ad(i-1, j) = d2_ad(i-1, j) - temp_ad2
6648 & , bd, gridstruct%sw_corner, gridstruct%se_corner&
6649 & , gridstruct%nw_corner, gridstruct%ne_corner)
6652 DO j=ad_to0,ad_from0,-1
6655 DO i=ad_to,ad_from,-1
6656 temp_ad1 = gridstruct%rarea(i, j)*d2_ad(i, j)
6657 fx2_ad(i, j) = fx2_ad(i, j) + temp_ad1
6658 fx2_ad(i+1, j) = fx2_ad(i+1, j) - temp_ad1
6659 fy2_ad(i, j) = fy2_ad(i, j) + temp_ad1
6660 fy2_ad(i, j+1) = fy2_ad(i, j+1) - temp_ad1
6668 DO j=je+nord+1,js-nord,-1
6669 DO i=ie+nord,is-nord,-1
6670 temp_ad0 = gridstruct%del6_u(i, j)*fy2_ad(i, j)
6671 d2_ad(i, j-1) = d2_ad(i, j-1) + temp_ad0
6672 d2_ad(i, j) = d2_ad(i, j) - temp_ad0
6678 & gridstruct%nested, bd, gridstruct&
6679 & %sw_corner, gridstruct%se_corner&
6680 & , gridstruct%nw_corner, &
6681 & gridstruct%ne_corner)
6682 DO j=je+nord,js-nord,-1
6683 DO i=ie+nord+1,is-nord,-1
6684 temp_ad = gridstruct%del6_v(i, j)*fx2_ad(i, j)
6685 d2_ad(i-1, j) = d2_ad(i-1, j) + temp_ad
6686 d2_ad(i, j) = d2_ad(i, j) - temp_ad
6692 & gridstruct%nested, bd, gridstruct&
6693 & %sw_corner, gridstruct%se_corner&
6694 & , gridstruct%nw_corner, &
6695 & gridstruct%ne_corner)
6697 IF (branch .EQ. 0)
THEN 6700 q_ad(i, j) = q_ad(i, j) + damp*d2_ad(i, j)
6707 q_ad(i, j) = q_ad(i, j) + d2_ad(i, j)
6736 & sw_corner, se_corner, nw_corner, ne_corner)
6739 INTEGER,
INTENT(IN) :: npx, npy, dir
6740 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
6741 REAL,
INTENT(INOUT) :: q_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
6742 LOGICAL,
INTENT(IN) :: nested, sw_corner, se_corner, nw_corner, &
6762 IF (.NOT.nested)
THEN 6763 IF (dir .EQ. 1)
THEN 6781 DO j=npy+
ng-1,npy,-1
6783 tmp_ad2 = q_ad(i, j)
6785 q_ad(npy-j, i-1+npx) = q_ad(npy-j, i-1+npx) + tmp_ad2
6790 IF (branch .EQ. 0)
THEN 6791 DO j=npy+
ng-1,npy,-1
6792 DO i=npx+
ng-1,npx,-1
6793 tmp_ad1 = q_ad(i, j)
6795 q_ad(j, 2*npx-1-i) = q_ad(j, 2*npx-1-i) + tmp_ad1
6800 IF (branch .EQ. 0)
THEN 6802 DO i=npx+
ng-1,npx,-1
6803 tmp_ad0 = q_ad(i, j)
6805 q_ad(npy-j, i-npx+1) = q_ad(npy-j, i-npx+1) + tmp_ad0
6810 IF (branch .EQ. 0)
THEN 6815 q_ad(j, 1-i) = q_ad(j, 1-i) + tmp_ad
6819 ELSE IF (dir .EQ. 2)
THEN 6837 DO j=npy+
ng-1,npy,-1
6839 tmp_ad6 = q_ad(i, j)
6841 q_ad(j+1-npx, npy-i) = q_ad(j+1-npx, npy-i) + tmp_ad6
6846 IF (branch .EQ. 0)
THEN 6847 DO j=npy+
ng-1,npy,-1
6848 DO i=npx+
ng-1,npx,-1
6849 tmp_ad5 = q_ad(i, j)
6851 q_ad(2*npy-1-j, i) = q_ad(2*npy-1-j, i) + tmp_ad5
6856 IF (branch .EQ. 0)
THEN 6858 DO i=npx+
ng-1,npx,-1
6859 tmp_ad4 = q_ad(i, j)
6861 q_ad(npy+j-1, npx-i) = q_ad(npy+j-1, npx-i) + tmp_ad4
6866 IF (branch .EQ. 0)
THEN 6869 tmp_ad3 = q_ad(i, j)
6871 q_ad(1-j, i) = q_ad(1-j, i) + tmp_ad3
6878 SUBROUTINE deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy, &
6879 & gridstruct, bd, mass)
6888 TYPE(FV_GRID_BOUNDS_TYPE),
INTENT(IN) :: bd
6890 INTEGER,
INTENT(IN) :: nord
6891 INTEGER,
INTENT(IN) :: is, ie, js, je, npx, npy
6892 REAL,
INTENT(IN) :: damp
6894 REAL,
INTENT(IN) :: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
6895 TYPE(FV_GRID_TYPE),
INTENT(IN),
TARGET :: gridstruct
6897 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
6899 REAL,
INTENT(INOUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je), fy(bd%is:bd%&
6900 & ie, bd%js:bd%je+1)
6902 REAL :: fx2(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2(bd%isd:bd%ied, bd%&
6904 REAL :: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
6906 INTEGER :: i, j, n, nt, i1, i2, j1, j2
6912 IF (.NOT.
PRESENT(mass))
THEN 6915 d2(i, j) = damp*q(i, j)
6925 IF (nord .GT. 0)
CALL copy_corners(d2, npx, npy, 1, gridstruct%&
6926 & nested, bd, gridstruct%sw_corner, &
6927 & gridstruct%se_corner, gridstruct%&
6928 & nw_corner, gridstruct%ne_corner)
6929 DO j=js-nord,je+nord
6930 DO i=is-nord,ie+nord+1
6931 fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i-1, j)-d2(i, j))
6934 IF (nord .GT. 0)
CALL copy_corners(d2, npx, npy, 2, gridstruct%&
6935 & nested, bd, gridstruct%sw_corner, &
6936 & gridstruct%se_corner, gridstruct%&
6937 & nw_corner, gridstruct%ne_corner)
6938 DO j=js-nord,je+nord+1
6939 DO i=is-nord,ie+nord
6940 fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j-1)-d2(i, j))
6943 IF (nord .GT. 0)
THEN 6949 DO j=js-nt-1,je+nt+1
6950 DO i=is-nt-1,ie+nt+1
6951 d2(i, j) = (fx2(i, j)-fx2(i+1, j)+fy2(i, j)-fy2(i, j+1))*&
6952 & gridstruct%rarea(i, j)
6955 CALL copy_corners(d2, npx, npy, 1, gridstruct%nested, bd, &
6956 & gridstruct%sw_corner, gridstruct%se_corner, &
6957 & gridstruct%nw_corner, gridstruct%ne_corner)
6960 fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i, j)-d2(i-1, j))
6963 CALL copy_corners(d2, npx, npy, 2, gridstruct%nested, bd, &
6964 & gridstruct%sw_corner, gridstruct%se_corner, &
6965 & gridstruct%nw_corner, gridstruct%ne_corner)
6968 fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j)-d2(i, j-1))
6976 IF (
PRESENT(mass))
THEN 6981 fx(i, j) = fx(i, j) + damp2*(mass(i-1, j)+mass(i, j))*fx2(i, j&
6987 fy(i, j) = fy(i, j) + damp2*(mass(i, j-1)+mass(i, j))*fy2(i, j&
6994 fx(i, j) = fx(i, j) + fx2(i, j)
6999 fy(i, j) = fy(i, j) + fy2(i, j)
7026 SUBROUTINE fv_tp_2d_fwd(q, crx, cry, npx, npy, hord, fx, fy, xfx, &
7027 & yfx, gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
7030 INTEGER,
INTENT(IN) :: npx, npy
7031 INTEGER,
INTENT(IN) :: hord
7033 REAL,
INTENT(IN) :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed)
7035 REAL,
INTENT(IN) :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed)
7037 REAL,
INTENT(IN) :: cry(bd%isd:bd%ied, bd%js:bd%je+1)
7039 REAL,
INTENT(IN) :: yfx(bd%isd:bd%ied, bd%js:bd%je+1)
7040 REAL,
INTENT(IN) :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
7041 REAL,
INTENT(IN) :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
7043 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
7045 REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
7047 REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
7051 REAL,
OPTIONAL,
INTENT(IN) :: mfx(bd%is:bd%ie+1, bd%js:bd%je)
7053 REAL,
OPTIONAL,
INTENT(IN) :: mfy(bd%is:bd%ie, bd%js:bd%je+1)
7054 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
7055 REAL,
OPTIONAL,
INTENT(IN) :: damp_c
7056 INTEGER,
OPTIONAL,
INTENT(IN) :: nord
7058 INTEGER :: ord_ou, ord_in
7059 REAL :: q_i(bd%isd:bd%ied, bd%js:bd%je)
7060 REAL :: q_j(bd%is:bd%ie, bd%jsd:bd%jed)
7061 REAL :: fx2(bd%is:bd%ie+1, bd%jsd:bd%jed)
7062 REAL :: fy2(bd%isd:bd%ied, bd%js:bd%je+1)
7063 REAL :: fyy(bd%isd:bd%ied, bd%js:bd%je+1)
7064 REAL :: fx1(bd%is:bd%ie+1, bd%jsd:bd%jed)
7067 INTEGER :: is, ie, js, je, isd, ied, jsd, jed
7077 IF (hord .EQ. 10)
THEN 7083 IF (.NOT.gridstruct%nested)
THEN 7085 & gridstruct%sw_corner, gridstruct%se_corner, &
7086 & gridstruct%nw_corner, gridstruct%ne_corner)
7091 CALL yppm_fwd(fy2, q, cry, ord_in, isd, ied, isd, ied, js, je, &
7092 & jsd, jed, npx, npy, gridstruct%dya, gridstruct%nested, &
7093 & gridstruct%grid_type)
7096 fyy(i, j) = yfx(i, j)*fy2(i, j)
7101 q_i(i, j) = (q(i, j)*gridstruct%area(i, j)+fyy(i, j)-fyy(i, j+1)&
7105 CALL xppm_fwd(fx, q_i, crx(is:ie+1, js:je), ord_ou, is, ie, isd, &
7106 & ied, js, je, jsd, jed, npx, npy, gridstruct%dxa, &
7107 & gridstruct%nested, gridstruct%grid_type)
7108 IF (.NOT.gridstruct%nested)
THEN 7110 & gridstruct%sw_corner, gridstruct%se_corner, &
7111 & gridstruct%nw_corner, gridstruct%ne_corner)
7116 CALL xppm_fwd(fx2, q, crx, ord_in, is, ie, isd, ied, jsd, jed, &
7117 & jsd, jed, npx, npy, gridstruct%dxa, gridstruct%nested, &
7118 & gridstruct%grid_type)
7122 fx1(i, j) = xfx(i, j)*fx2(i, j)
7127 q_j(i, j) = (q(i, j)*gridstruct%area(i, j)+fx1(i, j)-fx1(i+1, j)&
7131 CALL yppm_fwd(fy, q_j, cry, ord_ou, is, ie, isd, ied, js, je, jsd&
7132 & , jed, npx, npy, gridstruct%dya, gridstruct%nested, &
7133 & gridstruct%grid_type)
7137 IF (
PRESENT(mfx) .AND.
PRESENT(mfy))
THEN 7144 fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*mfx(i, j)
7150 fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*mfy(i, j)
7153 IF (
PRESENT(nord) .AND.
PRESENT(damp_c) .AND.
PRESENT(mass))
THEN 7154 IF (damp_c .GT. 1.e-4)
THEN 7155 damp = (damp_c*gridstruct%da_min)**(nord+1)
7157 & , fx, fy, gridstruct, bd, mass)
7195 fx(i, j) = 0.5*(fx(i, j)+fx2(i, j))*xfx(i, j)
7201 fy(i, j) = 0.5*(fy(i, j)+fy2(i, j))*yfx(i, j)
7204 IF (
PRESENT(nord) .AND.
PRESENT(damp_c))
THEN 7205 IF (damp_c .GT. 1.e-4)
THEN 7206 damp = (damp_c*gridstruct%da_min)**(nord+1)
7208 & , fx, fy, gridstruct, bd)
7263 SUBROUTINE fv_tp_2d_bwd(q, q_ad, crx, crx_ad, cry, cry_ad, npx, npy&
7264 & , hord, fx, fx_ad, fy, fy_ad, xfx, xfx_ad, yfx, yfx_ad, gridstruct, &
7265 & bd, ra_x, ra_x_ad, ra_y, ra_y_ad, mfx, mfx_ad, mfy, mfy_ad, mass, &
7266 & mass_ad, nord, damp_c)
7269 INTEGER,
INTENT(IN) :: npx, npy
7270 INTEGER,
INTENT(IN) :: hord
7271 REAL,
INTENT(IN) :: crx(bd%is:bd%ie+1, bd%jsd:bd%jed)
7272 REAL :: crx_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
7273 REAL,
INTENT(IN) :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed)
7274 REAL :: xfx_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
7275 REAL,
INTENT(IN) :: cry(bd%isd:bd%ied, bd%js:bd%je+1)
7276 REAL :: cry_ad(bd%isd:bd%ied, bd%js:bd%je+1)
7277 REAL,
INTENT(IN) :: yfx(bd%isd:bd%ied, bd%js:bd%je+1)
7278 REAL :: yfx_ad(bd%isd:bd%ied, bd%js:bd%je+1)
7279 REAL,
INTENT(IN) :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
7280 REAL :: ra_x_ad(bd%is:bd%ie, bd%jsd:bd%jed)
7281 REAL,
INTENT(IN) :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
7282 REAL :: ra_y_ad(bd%isd:bd%ied, bd%js:bd%je)
7283 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
7284 REAL,
INTENT(INOUT) :: q_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
7285 REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
7286 REAL :: fx_ad(bd%is:bd%ie+1, bd%js:bd%je)
7287 REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
7288 REAL :: fy_ad(bd%is:bd%ie, bd%js:bd%je+1)
7290 REAL,
OPTIONAL,
INTENT(IN) :: mfx(bd%is:bd%ie+1, bd%js:bd%je)
7291 REAL,
OPTIONAL :: mfx_ad(bd%is:bd%ie+1, bd%js:bd%je)
7292 REAL,
OPTIONAL,
INTENT(IN) :: mfy(bd%is:bd%ie, bd%js:bd%je+1)
7293 REAL,
OPTIONAL :: mfy_ad(bd%is:bd%ie, bd%js:bd%je+1)
7294 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
7295 REAL,
OPTIONAL :: mass_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
7296 REAL,
OPTIONAL,
INTENT(IN) :: damp_c
7297 INTEGER,
OPTIONAL,
INTENT(IN) :: nord
7298 INTEGER :: ord_ou, ord_in
7299 REAL :: q_i(bd%isd:bd%ied, bd%js:bd%je)
7300 REAL :: q_i_ad(bd%isd:bd%ied, bd%js:bd%je)
7301 REAL :: q_j(bd%is:bd%ie, bd%jsd:bd%jed)
7302 REAL :: q_j_ad(bd%is:bd%ie, bd%jsd:bd%jed)
7303 REAL :: fx2(bd%is:bd%ie+1, bd%jsd:bd%jed)
7304 REAL :: fx2_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
7305 REAL :: fy2(bd%isd:bd%ied, bd%js:bd%je+1)
7306 REAL :: fy2_ad(bd%isd:bd%ied, bd%js:bd%je+1)
7307 REAL :: fyy(bd%isd:bd%ied, bd%js:bd%je+1)
7308 REAL :: fyy_ad(bd%isd:bd%ied, bd%js:bd%je+1)
7309 REAL :: fx1(bd%is:bd%ie+1, bd%jsd:bd%jed)
7310 REAL :: fx1_ad(bd%is:bd%ie+1, bd%jsd:bd%jed)
7313 INTEGER :: is, ie, js, je, isd, ied, jsd, jed
7343 IF (branch .LT. 3)
THEN 7344 IF (branch .EQ. 0)
THEN 7346 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7347 CALL poprealarray(q_j, (bd%ie-bd%is+1)*(bd%jed-bd%jsd+1))
7351 CALL poprealarray(fx1, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7352 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7353 CALL poprealarray(fyy, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7354 ELSE IF (branch .EQ. 1)
THEN 7356 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7357 CALL poprealarray(q_j, (bd%ie-bd%is+1)*(bd%jed-bd%jsd+1))
7361 CALL poprealarray(fx1, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7362 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7363 CALL poprealarray(fyy, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7365 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7366 CALL poprealarray(q_j, (bd%ie-bd%is+1)*(bd%jed-bd%jsd+1))
7367 CALL poprealarray(fx1, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7368 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7369 CALL poprealarray(fyy, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7371 damp = (damp_c*gridstruct%da_min)**(nord+1)
7375 CALL deln_flux_bwd(nord, is, ie, js, je, npx, npy, damp, q, &
7376 & q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd, mass&
7383 temp_ad4 = 0.5*mfy(i, j)*fy_ad(i, j)
7384 fy2_ad(i, j) = fy2_ad(i, j) + temp_ad4
7385 mfy_ad(i, j) = mfy_ad(i, j) + 0.5*(fy(i, j)+fy2(i, j))*fy_ad(i&
7387 fy_ad(i, j) = temp_ad4
7394 temp_ad3 = 0.5*mfx(i, j)*fx_ad(i, j)
7395 fx2_ad(i, j) = fx2_ad(i, j) + temp_ad3
7396 mfx_ad(i, j) = mfx_ad(i, j) + 0.5*(fx(i, j)+fx2(i, j))*fx_ad(i&
7398 fx_ad(i, j) = temp_ad3
7402 IF (branch .EQ. 3)
THEN 7404 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7405 CALL poprealarray(q_j, (bd%ie-bd%is+1)*(bd%jed-bd%jsd+1))
7409 CALL poprealarray(fx1, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7410 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7411 CALL poprealarray(fyy, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7412 ELSE IF (branch .EQ. 4)
THEN 7414 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7415 CALL poprealarray(q_j, (bd%ie-bd%is+1)*(bd%jed-bd%jsd+1))
7419 CALL poprealarray(fx1, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7420 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7421 CALL poprealarray(fyy, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7423 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7424 CALL poprealarray(q_j, (bd%ie-bd%is+1)*(bd%jed-bd%jsd+1))
7425 CALL poprealarray(fx1, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7426 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7427 CALL poprealarray(fyy, (bd%ied-bd%isd+1)*(bd%je-bd%js+2))
7429 damp = (damp_c*gridstruct%da_min)**(nord+1)
7433 CALL deln_flux_bwd(nord, is, ie, js, je, npx, npy, damp, q, &
7434 & q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd)
7440 temp_ad2 = 0.5*yfx(i, j)*fy_ad(i, j)
7441 fy2_ad(i, j) = fy2_ad(i, j) + temp_ad2
7442 yfx_ad(i, j) = yfx_ad(i, j) + 0.5*(fy(i, j)+fy2(i, j))*fy_ad(i&
7444 fy_ad(i, j) = temp_ad2
7451 temp_ad1 = 0.5*xfx(i, j)*fx_ad(i, j)
7452 fx2_ad(i, j) = fx2_ad(i, j) + temp_ad1
7453 xfx_ad(i, j) = xfx_ad(i, j) + 0.5*(fx(i, j)+fx2(i, j))*fx_ad(i&
7455 fx_ad(i, j) = temp_ad1
7464 CALL yppm_bwd(fy, fy_ad, q_j, q_j_ad, cry, cry_ad, ord_ou, is, ie&
7465 & , isd, ied, js, je, jsd, jed, npx, npy, gridstruct%dya, &
7466 & gridstruct%nested, gridstruct%grid_type)
7470 temp_ad0 = q_j_ad(i, j)/ra_x(i, j)
7471 q_ad(i, j) = q_ad(i, j) + gridstruct%area(i, j)*temp_ad0
7472 fx1_ad(i, j) = fx1_ad(i, j) + temp_ad0
7473 fx1_ad(i+1, j) = fx1_ad(i+1, j) - temp_ad0
7474 ra_x_ad(i, j) = ra_x_ad(i, j) - (gridstruct%area(i, j)*q(i, j)+&
7475 & fx1(i, j)-fx1(i+1, j))*temp_ad0/ra_x(i, j)
7479 CALL poprealarray(fx2, (bd%ie-bd%is+2)*(bd%jed-bd%jsd+1))
7482 xfx_ad(i, j) = xfx_ad(i, j) + fx2(i, j)*fx1_ad(i, j)
7483 fx2_ad(i, j) = fx2_ad(i, j) + xfx(i, j)*fx1_ad(i, j)
7487 CALL xppm_bwd(fx2, fx2_ad, q, q_ad, crx, crx_ad, ord_in, is, ie, &
7488 & isd, ied, jsd, jed, jsd, jed, npx, npy, gridstruct%dxa, &
7489 & gridstruct%nested, gridstruct%grid_type)
7492 & gridstruct%nested, bd, &
7493 & gridstruct%sw_corner, &
7494 & gridstruct%se_corner, &
7495 & gridstruct%nw_corner, &
7496 & gridstruct%ne_corner)
7498 CALL xppm_bwd(fx, fx_ad, q_i, q_i_ad, crx(is:ie+1, js:je), crx_ad&
7499 & (is:ie+1, js:je), ord_ou, is, ie, isd, ied, js, je, jsd, &
7500 & jed, npx, npy, gridstruct%dxa, gridstruct%nested, &
7501 & gridstruct%grid_type)
7505 temp_ad = q_i_ad(i, j)/ra_y(i, j)
7506 q_ad(i, j) = q_ad(i, j) + gridstruct%area(i, j)*temp_ad
7507 fyy_ad(i, j) = fyy_ad(i, j) + temp_ad
7508 fyy_ad(i, j+1) = fyy_ad(i, j+1) - temp_ad
7509 ra_y_ad(i, j) = ra_y_ad(i, j) - (gridstruct%area(i, j)*q(i, j)+&
7510 & fyy(i, j)-fyy(i, j+1))*temp_ad/ra_y(i, j)
7516 yfx_ad(i, j) = yfx_ad(i, j) + fy2(i, j)*fyy_ad(i, j)
7517 fy2_ad(i, j) = fy2_ad(i, j) + yfx(i, j)*fyy_ad(i, j)
7521 CALL yppm_bwd(fy2, fy2_ad, q, q_ad, cry, cry_ad, ord_in, isd, ied&
7522 & , isd, ied, js, je, jsd, jed, npx, npy, gridstruct%dya, &
7523 & gridstruct%nested, gridstruct%grid_type)
7526 & gridstruct%nested, bd, &
7527 & gridstruct%sw_corner, &
7528 & gridstruct%se_corner, &
7529 & gridstruct%nw_corner, &
7530 & gridstruct%ne_corner)
7552 SUBROUTINE xppm_fwd(flux, q, c, iord, is, ie, isd, ied, jfirst, &
7553 & jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
7555 INTEGER,
INTENT(IN) :: is, ie, isd, ied, jsd, jed
7557 INTEGER,
INTENT(IN) :: jfirst, jlast
7558 INTEGER,
INTENT(IN) :: iord
7559 INTEGER,
INTENT(IN) :: npx, npy
7560 REAL,
INTENT(IN) :: q(isd:ied, jfirst:jlast)
7562 REAL,
INTENT(IN) :: c(is:ie+1, jfirst:jlast)
7563 REAL,
INTENT(IN) :: dxa(isd:ied, jsd:jed)
7564 LOGICAL,
INTENT(IN) :: nested
7565 INTEGER,
INTENT(IN) :: grid_type
7568 REAL :: flux(is:ie+1, jfirst:jlast)
7570 REAL,
DIMENSION(is-1:ie+1) :: bl, br, b0
7572 REAL,
DIMENSION(is:ie+1) :: fx0, fx1
7573 LOGICAL,
DIMENSION(is-1:ie+1) :: smt5, smt6
7574 REAL :: al(is-1:ie+2)
7575 REAL :: dm(is-2:ie+2)
7576 REAL :: dq(is-3:ie+2)
7577 INTEGER :: i, j, ie3, is1, ie1
7578 REAL :: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
7603 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 7604 IF (3 .LT. is - 1)
THEN 7609 IF (npx - 2 .GT. ie + 2)
THEN 7626 IF (iord .LT. 8 .OR. iord .EQ. 333)
THEN 7631 al(i) =
p1*(q1(i-1)+q1(i)) +
p2*(q1(i-2)+q1(i+1))
7633 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 7636 al(0) =
c1*q1(-2) +
c2*q1(-1) +
c3*q1(0)
7638 al(1) = 0.5*(((2.*dxa(0, j)+dxa(-1, j))*q1(0)-dxa(0, j)*q1(-&
7639 & 1))/(dxa(-1, j)+dxa(0, j))+((2.*dxa(1, j)+dxa(2, j))*q1(1)&
7640 & -dxa(1, j)*q1(2))/(dxa(1, j)+dxa(2, j)))
7642 al(2) =
c3*q1(1) +
c2*q1(2) +
c1*q1(3)
7647 IF (ie + 1 .EQ. npx)
THEN 7649 al(npx-1) =
c1*q1(npx-3) +
c2*q1(npx-2) +
c3*q1(npx-1)
7651 al(npx) = 0.5*(((2.*dxa(npx-1, j)+dxa(npx-2, j))*q1(npx-1)-&
7652 & dxa(npx-1, j)*q1(npx-2))/(dxa(npx-2, j)+dxa(npx-1, j))+((&
7653 & 2.*dxa(npx, j)+dxa(npx+1, j))*q1(npx)-dxa(npx, j)*q1(npx+1&
7654 & ))/(dxa(npx, j)+dxa(npx+1, j)))
7656 al(npx+1) =
c3*q1(npx) +
c2*q1(npx+1) +
c1*q1(npx+2)
7664 IF (iord .EQ. 1)
THEN 7666 IF (c(i, j) .GT. 0.)
THEN 7668 flux(i, j) = q1(i-1)
7677 ELSE IF (iord .EQ. 2)
THEN 7682 IF (xt .GT. 0.)
THEN 7686 flux(i, j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-&
7693 flux(i, j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-&
7699 ELSE IF (iord .EQ. 333)
THEN 7704 IF (xt .GT. 0.)
THEN 7706 flux(i, j) = (2.0*q1(i)+5.0*q1(i-1)-q1(i-2))/6.0 - 0.5*xt*&
7707 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i)-2.0*q1(i-1)+q1(i-2))
7711 flux(i, j) = (2.0*q1(i-1)+5.0*q1(i)-q1(i+1))/6.0 - 0.5*xt*&
7712 & (q1(i)-q1(i-1)) + xt*xt/6.0*(q1(i+1)-2.0*q1(i)+q1(i-1))
7750 SUBROUTINE xppm_bwd(flux, flux_ad, q, q_ad, c, c_ad, iord, is, ie, &
7751 & isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
7753 INTEGER,
INTENT(IN) :: is, ie, isd, ied, jsd, jed
7754 INTEGER,
INTENT(IN) :: jfirst, jlast
7755 INTEGER,
INTENT(IN) :: iord
7756 INTEGER,
INTENT(IN) :: npx, npy
7757 REAL,
INTENT(IN) :: q(isd:ied, jfirst:jlast)
7758 REAL :: q_ad(isd:ied, jfirst:jlast)
7759 REAL,
INTENT(IN) :: c(is:ie+1, jfirst:jlast)
7760 REAL :: c_ad(is:ie+1, jfirst:jlast)
7761 REAL,
INTENT(IN) :: dxa(isd:ied, jsd:jed)
7762 LOGICAL,
INTENT(IN) :: nested
7763 INTEGER,
INTENT(IN) :: grid_type
7764 REAL :: flux(is:ie+1, jfirst:jlast)
7765 REAL :: flux_ad(is:ie+1, jfirst:jlast)
7766 REAL,
DIMENSION(is-1:ie+1) :: bl, br, b0
7768 REAL :: q1_ad(isd:ied)
7769 REAL,
DIMENSION(is:ie+1) :: fx0, fx1
7770 LOGICAL,
DIMENSION(is-1:ie+1) :: smt5, smt6
7771 REAL :: al(is-1:ie+2)
7772 REAL :: al_ad(is-1:ie+2)
7773 REAL :: dm(is-2:ie+2)
7774 REAL :: dq(is-3:ie+2)
7775 INTEGER :: i, j, ie3, is1, ie1
7776 REAL :: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2
7777 REAL :: xt_ad, qtmp_ad
7827 DO j=jlast,jfirst,-1
7829 IF (branch .LT. 2)
THEN 7830 IF (branch .EQ. 0)
THEN 7833 IF (branch .EQ. 0)
THEN 7835 q1_ad(i) = q1_ad(i) + flux_ad(i, j)
7839 q1_ad(i-1) = q1_ad(i-1) + flux_ad(i, j)
7846 IF (branch .EQ. 0)
THEN 7850 temp0 = al(i) + al(i+1) - 2*qtmp
7851 temp_ad5 = (xt+1.)*flux_ad(i, j)
7852 temp_ad6 = xt*temp_ad5
7853 qtmp_ad = flux_ad(i, j) - temp_ad5 - 2*temp_ad6
7854 xt_ad = temp0*temp_ad5 + (al(i)-qtmp+xt*temp0)*flux_ad(i, &
7856 al_ad(i) = al_ad(i) + temp_ad6 + temp_ad5
7857 al_ad(i+1) = al_ad(i+1) + temp_ad6
7860 q1_ad(i) = q1_ad(i) + qtmp_ad
7865 temp = al(i-1) + al(i) - 2*qtmp
7866 temp_ad3 = (1.-xt)*flux_ad(i, j)
7867 temp_ad4 = -(xt*temp_ad3)
7868 qtmp_ad = flux_ad(i, j) - temp_ad3 - 2*temp_ad4
7869 xt_ad = -(temp*temp_ad3) - (al(i)-qtmp-xt*temp)*flux_ad(i&
7871 al_ad(i) = al_ad(i) + temp_ad4 + temp_ad3
7872 al_ad(i-1) = al_ad(i-1) + temp_ad4
7875 q1_ad(i-1) = q1_ad(i-1) + qtmp_ad
7877 c_ad(i, j) = c_ad(i, j) + xt_ad
7880 ELSE IF (branch .EQ. 2)
THEN 7883 IF (branch .EQ. 0)
THEN 7886 temp_ad10 = flux_ad(i, j)/6.0
7887 temp_ad11 = -(0.5*xt*flux_ad(i, j))
7888 temp_ad12 = xt**2*flux_ad(i, j)/6.0
7889 q1_ad(i-1) = q1_ad(i-1) + temp_ad12 - temp_ad11 + 2.0*&
7891 q1_ad(i) = q1_ad(i) + temp_ad11 - 2.0*temp_ad12 + 5.0*&
7893 q1_ad(i+1) = q1_ad(i+1) + temp_ad12 - temp_ad10
7894 xt_ad = ((q1(i+1)-2.0*q1(i)+q1(i-1))*2*xt/6.0-0.5*(q1(i)-q1(&
7895 & i-1)))*flux_ad(i, j)
7900 temp_ad7 = flux_ad(i, j)/6.0
7901 temp_ad8 = -(0.5*xt*flux_ad(i, j))
7902 temp_ad9 = xt**2*flux_ad(i, j)/6.0
7903 q1_ad(i) = q1_ad(i) + temp_ad9 + temp_ad8 + 2.0*temp_ad7
7904 q1_ad(i-1) = q1_ad(i-1) + 5.0*temp_ad7 - temp_ad8 - 2.0*&
7906 q1_ad(i-2) = q1_ad(i-2) + temp_ad9 - temp_ad7
7907 xt_ad = ((q1(i)-2.0*q1(i-1)+q1(i-2))*2*xt/6.0-0.5*(q1(i)-q1(&
7908 & i-1)))*flux_ad(i, j)
7911 c_ad(i, j) = c_ad(i, j) + xt_ad
7913 ELSE IF (branch .NE. 3)
THEN 7917 IF (branch .EQ. 0)
THEN 7919 q1_ad(npx) = q1_ad(npx) +
c3*al_ad(npx+1)
7920 q1_ad(npx+1) = q1_ad(npx+1) +
c2*al_ad(npx+1)
7921 q1_ad(npx+2) = q1_ad(npx+2) +
c1*al_ad(npx+1)
7924 temp_ad1 = 0.5*al_ad(npx)/(dxa(npx-2, j)+dxa(npx-1, j))
7925 temp_ad2 = 0.5*al_ad(npx)/(dxa(npx, j)+dxa(npx+1, j))
7926 q1_ad(npx-1) = q1_ad(npx-1) + (dxa(npx-1, j)*2.+dxa(npx-2, j))*&
7928 q1_ad(npx-2) = q1_ad(npx-2) - dxa(npx-1, j)*temp_ad1
7929 q1_ad(npx) = q1_ad(npx) + (dxa(npx, j)*2.+dxa(npx+1, j))*&
7931 q1_ad(npx+1) = q1_ad(npx+1) - dxa(npx, j)*temp_ad2
7934 q1_ad(npx-3) = q1_ad(npx-3) +
c1*al_ad(npx-1)
7935 q1_ad(npx-2) = q1_ad(npx-2) +
c2*al_ad(npx-1)
7936 q1_ad(npx-1) = q1_ad(npx-1) +
c3*al_ad(npx-1)
7938 ELSE IF (branch .NE. 1)
THEN 7942 IF (branch .EQ. 0)
THEN 7944 q1_ad(1) = q1_ad(1) +
c3*al_ad(2)
7945 q1_ad(2) = q1_ad(2) +
c2*al_ad(2)
7946 q1_ad(3) = q1_ad(3) +
c1*al_ad(2)
7949 temp_ad = 0.5*al_ad(1)/(dxa(-1, j)+dxa(0, j))
7950 temp_ad0 = 0.5*al_ad(1)/(dxa(1, j)+dxa(2, j))
7951 q1_ad(0) = q1_ad(0) + (dxa(0, j)*2.+dxa(-1, j))*temp_ad
7952 q1_ad(-1) = q1_ad(-1) - dxa(0, j)*temp_ad
7953 q1_ad(1) = q1_ad(1) + (dxa(1, j)*2.+dxa(2, j))*temp_ad0
7954 q1_ad(2) = q1_ad(2) - dxa(1, j)*temp_ad0
7957 q1_ad(-2) = q1_ad(-2) +
c1*al_ad(0)
7958 q1_ad(-1) = q1_ad(-1) +
c2*al_ad(0)
7959 q1_ad(0) = q1_ad(0) +
c3*al_ad(0)
7964 q1_ad(i-1) = q1_ad(i-1) +
p1*al_ad(i)
7965 q1_ad(i) = q1_ad(i) +
p1*al_ad(i)
7966 q1_ad(i-2) = q1_ad(i-2) +
p2*al_ad(i)
7967 q1_ad(i+1) = q1_ad(i+1) +
p2*al_ad(i)
7972 q_ad(i, j) = q_ad(i, j) + q1_ad(i)
7998 SUBROUTINE yppm_fwd(flux, q, c, jord, ifirst, ilast, isd, ied, js, &
7999 & je, jsd, jed, npx, npy, dya, nested, grid_type)
8002 INTEGER,
INTENT(IN) :: ifirst, ilast
8003 INTEGER,
INTENT(IN) :: isd, ied, js, je, jsd, jed
8004 INTEGER,
INTENT(IN) :: jord
8005 INTEGER,
INTENT(IN) :: npx, npy
8006 REAL,
INTENT(IN) :: q(ifirst:ilast, jsd:jed)
8008 REAL,
INTENT(IN) :: c(isd:ied, js:je+1)
8010 REAL :: flux(ifirst:ilast, js:je+1)
8011 REAL,
INTENT(IN) :: dya(isd:ied, jsd:jed)
8012 LOGICAL,
INTENT(IN) :: nested
8013 INTEGER,
INTENT(IN) :: grid_type
8015 REAL :: dm(ifirst:ilast, js-2:je+2)
8016 REAL :: al(ifirst:ilast, js-1:je+2)
8017 REAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: bl, br, b0
8018 REAL :: dq(ifirst:ilast, js-3:je+2)
8019 REAL,
DIMENSION(ifirst:ilast) :: fx0, fx1
8020 LOGICAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: smt5, smt6
8021 REAL :: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
8022 INTEGER :: i, j, js1, je3, je1
8046 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 8047 IF (3 .LT. js - 1)
THEN 8052 IF (npy - 2 .GT. je + 2)
THEN 8065 IF (jord .LT. 8 .OR. jord .EQ. 333)
THEN 8068 al(i, j) =
p1*(q(i, j-1)+q(i, j)) +
p2*(q(i, j-2)+q(i, j+1))
8071 IF (.NOT.nested .AND.
grid_type .LT. 3)
THEN 8074 al(i, 0) =
c1*q(i, -2) +
c2*q(i, -1) +
c3*q(i, 0)
8075 al(i, 1) = 0.5*(((2.*dya(i, 0)+dya(i, -1))*q(i, 0)-dya(i, 0)&
8076 & *q(i, -1))/(dya(i, -1)+dya(i, 0))+((2.*dya(i, 1)+dya(i, 2)&
8077 & )*q(i, 1)-dya(i, 1)*q(i, 2))/(dya(i, 1)+dya(i, 2)))
8078 al(i, 2) =
c3*q(i, 1) +
c2*q(i, 2) +
c1*q(i, 3)
8084 IF (je + 1 .EQ. npy)
THEN 8086 al(i, npy-1) =
c1*q(i, npy-3) +
c2*q(i, npy-2) +
c3*q(i, npy&
8088 al(i, npy) = 0.5*(((2.*dya(i, npy-1)+dya(i, npy-2))*q(i, npy&
8089 & -1)-dya(i, npy-1)*q(i, npy-2))/(dya(i, npy-2)+dya(i, npy-1&
8090 & ))+((2.*dya(i, npy)+dya(i, npy+1))*q(i, npy)-dya(i, npy)*q&
8091 & (i, npy+1))/(dya(i, npy)+dya(i, npy+1)))
8092 al(i, npy+1) =
c3*q(i, npy) +
c2*q(i, npy+1) +
c1*q(i, npy+2&
8102 IF (jord .EQ. 1)
THEN 8105 IF (c(i, j) .GT. 0.)
THEN 8107 flux(i, j) = q(i, j-1)
8111 flux(i, j) = q(i, j)
8119 ELSE IF (jord .EQ. 2)
THEN 8126 IF (xt .GT. 0.)
THEN 8129 flux(i, j) = qtmp + (1.-xt)*(al(i, j)-qtmp-xt*(al(i, j-1)+&
8130 & al(i, j)-(qtmp+qtmp)))
8135 flux(i, j) = qtmp + (1.+xt)*(al(i, j)-qtmp+xt*(al(i, j)+al&
8136 & (i, j+1)-(qtmp+qtmp)))
8145 ELSE IF (jord .EQ. 333)
THEN 8151 IF (xt .GT. 0.)
THEN 8153 flux(i, j) = (2.0*q(i, j)+5.0*q(i, j-1)-q(i, j-2))/6.0 - &
8154 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j)-2.0*q(i&
8159 flux(i, j) = (2.0*q(i, j-1)+5.0*q(i, j)-q(i, j+1))/6.0 - &
8160 & 0.5*xt*(q(i, j)-q(i, j-1)) + xt*xt/6.0*(q(i, j+1)-2.0*q(&
8198 SUBROUTINE yppm_bwd(flux, flux_ad, q, q_ad, c, c_ad, jord, ifirst, &
8199 & ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
8201 INTEGER,
INTENT(IN) :: ifirst, ilast
8202 INTEGER,
INTENT(IN) :: isd, ied, js, je, jsd, jed
8203 INTEGER,
INTENT(IN) :: jord
8204 INTEGER,
INTENT(IN) :: npx, npy
8205 REAL,
INTENT(IN) :: q(ifirst:ilast, jsd:jed)
8206 REAL :: q_ad(ifirst:ilast, jsd:jed)
8207 REAL,
INTENT(IN) :: c(isd:ied, js:je+1)
8208 REAL :: c_ad(isd:ied, js:je+1)
8209 REAL :: flux(ifirst:ilast, js:je+1)
8210 REAL :: flux_ad(ifirst:ilast, js:je+1)
8211 REAL,
INTENT(IN) :: dya(isd:ied, jsd:jed)
8212 LOGICAL,
INTENT(IN) :: nested
8213 INTEGER,
INTENT(IN) :: grid_type
8214 REAL :: dm(ifirst:ilast, js-2:je+2)
8215 REAL :: al(ifirst:ilast, js-1:je+2)
8216 REAL :: al_ad(ifirst:ilast, js-1:je+2)
8217 REAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: bl, br, b0
8218 REAL :: dq(ifirst:ilast, js-3:je+2)
8219 REAL,
DIMENSION(ifirst:ilast) :: fx0, fx1
8220 LOGICAL,
DIMENSION(ifirst:ilast, js-1:je+1) :: smt5, smt6
8221 REAL :: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
8222 REAL :: xt_ad, qtmp_ad
8223 INTEGER :: i, j, js1, je3, je1
8266 IF (branch .LT. 2)
THEN 8267 IF (branch .EQ. 0)
THEN 8273 DO i=ilast,ifirst,-1
8275 IF (branch .EQ. 0)
THEN 8277 q_ad(i, j) = q_ad(i, j) + flux_ad(i, j)
8281 q_ad(i, j-1) = q_ad(i, j-1) + flux_ad(i, j)
8288 ELSE IF (branch .EQ. 2)
THEN 8294 DO i=ilast,ifirst,-1
8296 IF (branch .EQ. 0)
THEN 8300 temp0 = al(i, j) + al(i, j+1) - 2*qtmp
8301 temp_ad5 = (xt+1.)*flux_ad(i, j)
8302 temp_ad6 = xt*temp_ad5
8303 qtmp_ad = flux_ad(i, j) - temp_ad5 - 2*temp_ad6
8304 xt_ad = temp0*temp_ad5 + (al(i, j)-qtmp+xt*temp0)*flux_ad(i&
8306 al_ad(i, j) = al_ad(i, j) + temp_ad6 + temp_ad5
8307 al_ad(i, j+1) = al_ad(i, j+1) + temp_ad6
8309 q_ad(i, j) = q_ad(i, j) + qtmp_ad
8314 temp = al(i, j-1) + al(i, j) - 2*qtmp
8315 temp_ad3 = (1.-xt)*flux_ad(i, j)
8316 temp_ad4 = -(xt*temp_ad3)
8317 qtmp_ad = flux_ad(i, j) - temp_ad3 - 2*temp_ad4
8318 xt_ad = -(temp*temp_ad3) - (al(i, j)-qtmp-xt*temp)*flux_ad(i&
8320 al_ad(i, j) = al_ad(i, j) + temp_ad4 + temp_ad3
8321 al_ad(i, j-1) = al_ad(i, j-1) + temp_ad4
8323 q_ad(i, j-1) = q_ad(i, j-1) + qtmp_ad
8325 c_ad(i, j) = c_ad(i, j) + xt_ad
8329 IF (branch .EQ. 3)
THEN 8333 DO i=ilast,ifirst,-1
8335 IF (branch .EQ. 0)
THEN 8338 temp_ad10 = flux_ad(i, j)/6.0
8339 temp_ad11 = -(0.5*xt*flux_ad(i, j))
8340 temp_ad12 = xt**2*flux_ad(i, j)/6.0
8341 q_ad(i, j-1) = q_ad(i, j-1) + temp_ad12 - temp_ad11 + 2.0*&
8343 q_ad(i, j) = q_ad(i, j) + temp_ad11 - 2.0*temp_ad12 + 5.0*&
8345 q_ad(i, j+1) = q_ad(i, j+1) + temp_ad12 - temp_ad10
8346 xt_ad = ((q(i, j+1)-2.0*q(i, j)+q(i, j-1))*2*xt/6.0-0.5*(q&
8347 & (i, j)-q(i, j-1)))*flux_ad(i, j)
8352 temp_ad7 = flux_ad(i, j)/6.0
8353 temp_ad8 = -(0.5*xt*flux_ad(i, j))
8354 temp_ad9 = xt**2*flux_ad(i, j)/6.0
8355 q_ad(i, j) = q_ad(i, j) + temp_ad9 + temp_ad8 + 2.0*&
8357 q_ad(i, j-1) = q_ad(i, j-1) + 5.0*temp_ad7 - temp_ad8 - &
8359 q_ad(i, j-2) = q_ad(i, j-2) + temp_ad9 - temp_ad7
8360 xt_ad = ((q(i, j)-2.0*q(i, j-1)+q(i, j-2))*2*xt/6.0-0.5*(q&
8361 & (i, j)-q(i, j-1)))*flux_ad(i, j)
8364 c_ad(i, j) = c_ad(i, j) + xt_ad
8374 IF (branch .EQ. 0)
THEN 8375 DO i=ilast,ifirst,-1
8376 q_ad(i, npy) = q_ad(i, npy) +
c3*al_ad(i, npy+1)
8377 q_ad(i, npy+1) = q_ad(i, npy+1) +
c2*al_ad(i, npy+1)
8378 q_ad(i, npy+2) = q_ad(i, npy+2) +
c1*al_ad(i, npy+1)
8379 al_ad(i, npy+1) = 0.0
8380 temp_ad1 = 0.5*al_ad(i, npy)/(dya(i, npy-2)+dya(i, npy-1))
8381 temp_ad2 = 0.5*al_ad(i, npy)/(dya(i, npy)+dya(i, npy+1))
8382 q_ad(i, npy-1) = q_ad(i, npy-1) + (dya(i, npy-1)*2.+dya(i, npy-2&
8384 q_ad(i, npy-2) = q_ad(i, npy-2) - dya(i, npy-1)*temp_ad1
8385 q_ad(i, npy) = q_ad(i, npy) + (dya(i, npy)*2.+dya(i, npy+1))*&
8387 q_ad(i, npy+1) = q_ad(i, npy+1) - dya(i, npy)*temp_ad2
8389 q_ad(i, npy-3) = q_ad(i, npy-3) +
c1*al_ad(i, npy-1)
8390 q_ad(i, npy-2) = q_ad(i, npy-2) +
c2*al_ad(i, npy-1)
8391 q_ad(i, npy-1) = q_ad(i, npy-1) +
c3*al_ad(i, npy-1)
8392 al_ad(i, npy-1) = 0.0
8394 ELSE IF (branch .NE. 1)
THEN 8398 IF (branch .EQ. 0)
THEN 8399 DO i=ilast,ifirst,-1
8400 q_ad(i, 1) = q_ad(i, 1) +
c3*al_ad(i, 2)
8401 q_ad(i, 2) = q_ad(i, 2) +
c2*al_ad(i, 2)
8402 q_ad(i, 3) = q_ad(i, 3) +
c1*al_ad(i, 2)
8404 temp_ad = 0.5*al_ad(i, 1)/(dya(i, -1)+dya(i, 0))
8405 temp_ad0 = 0.5*al_ad(i, 1)/(dya(i, 1)+dya(i, 2))
8406 q_ad(i, 0) = q_ad(i, 0) + (dya(i, 0)*2.+dya(i, -1))*temp_ad
8407 q_ad(i, -1) = q_ad(i, -1) - dya(i, 0)*temp_ad
8408 q_ad(i, 1) = q_ad(i, 1) + (dya(i, 1)*2.+dya(i, 2))*temp_ad0
8409 q_ad(i, 2) = q_ad(i, 2) - dya(i, 1)*temp_ad0
8411 q_ad(i, -2) = q_ad(i, -2) +
c1*al_ad(i, 0)
8412 q_ad(i, -1) = q_ad(i, -1) +
c2*al_ad(i, 0)
8413 q_ad(i, 0) = q_ad(i, 0) +
c3*al_ad(i, 0)
8418 DO i=ilast,ifirst,-1
8419 q_ad(i, j-1) = q_ad(i, j-1) +
p1*al_ad(i, j)
8420 q_ad(i, j) = q_ad(i, j) +
p1*al_ad(i, j)
8421 q_ad(i, j-2) = q_ad(i, j-2) +
p2*al_ad(i, j)
8422 q_ad(i, j+1) = q_ad(i, j+1) +
p2*al_ad(i, j)
8448 SUBROUTINE deln_flux_fwd(nord, is, ie, js, je, npx, npy, damp, q, &
8449 & fx, fy, gridstruct, bd, mass)
8458 TYPE(FV_GRID_BOUNDS_TYPE),
INTENT(IN) :: bd
8460 INTEGER,
INTENT(IN) :: nord
8461 INTEGER,
INTENT(IN) :: is, ie, js, je, npx, npy
8462 REAL,
INTENT(IN) :: damp
8464 REAL,
INTENT(IN) :: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
8465 TYPE(FV_GRID_TYPE),
INTENT(IN),
TARGET :: gridstruct
8467 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
8469 REAL,
INTENT(INOUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je), fy(bd%is:bd%&
8470 & ie, bd%js:bd%je+1)
8472 REAL :: fx2(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2(bd%isd:bd%ied, bd%&
8474 REAL :: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
8476 INTEGER :: i, j, n, nt, i1, i2, j1, j2
8505 IF (.NOT.
PRESENT(mass))
THEN 8508 d2(i, j) = damp*q(i, j)
8520 IF (nord .GT. 0)
THEN 8522 & gridstruct%sw_corner, gridstruct%se_corner, &
8523 & gridstruct%nw_corner, gridstruct%ne_corner)
8528 DO j=js-nord,je+nord
8529 DO i=is-nord,ie+nord+1
8530 fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i-1, j)-d2(i, j))
8533 IF (nord .GT. 0)
THEN 8535 & gridstruct%sw_corner, gridstruct%se_corner, &
8536 & gridstruct%nw_corner, gridstruct%ne_corner)
8541 DO j=js-nord,je+nord+1
8542 DO i=is-nord,ie+nord
8543 fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j-1)-d2(i, j))
8546 IF (nord .GT. 0)
THEN 8552 ad_from0 = js - nt - 1
8553 DO j=ad_from0,je+nt+1
8554 ad_from = is - nt - 1
8555 DO i=ad_from,ie+nt+1
8556 d2(i, j) = (fx2(i, j)-fx2(i+1, j)+fy2(i, j)-fy2(i, j+1))*&
8557 & gridstruct%rarea(i, j)
8565 & , gridstruct%sw_corner, gridstruct%se_corner&
8566 & , gridstruct%nw_corner, gridstruct%ne_corner)
8570 DO i=ad_from1,ie+nt+1
8571 fx2(i, j) = gridstruct%del6_v(i, j)*(d2(i, j)-d2(i-1, j))
8579 & , gridstruct%sw_corner, gridstruct%se_corner&
8580 & , gridstruct%nw_corner, gridstruct%ne_corner)
8582 DO j=ad_from4,je+nt+1
8585 fy2(i, j) = gridstruct%del6_u(i, j)*(d2(i, j)-d2(i, j-1))
8600 IF (
PRESENT(mass))
THEN 8606 fx(i, j) = fx(i, j) + damp2*(mass(i-1, j)+mass(i, j))*fx2(i, j&
8613 fy(i, j) = fy(i, j) + damp2*(mass(i, j-1)+mass(i, j))*fy2(i, j&
8618 CALL pushrealarray(fx2, (bd%ied-bd%isd+2)*(bd%jed-bd%jsd+1))
8619 CALL pushrealarray(fy2, (bd%ied-bd%isd+1)*(bd%jed-bd%jsd+2))
8626 fx(i, j) = fx(i, j) + fx2(i, j)
8632 fy(i, j) = fy(i, j) + fy2(i, j)
8659 SUBROUTINE deln_flux_bwd(nord, is, ie, js, je, npx, npy, damp, q, &
8660 & q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd, mass, mass_ad)
8662 TYPE(FV_GRID_BOUNDS_TYPE),
INTENT(IN) :: bd
8663 INTEGER,
INTENT(IN) :: nord
8664 INTEGER,
INTENT(IN) :: is, ie, js, je, npx, npy
8665 REAL,
INTENT(IN) :: damp
8666 REAL,
INTENT(IN) :: q(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
8667 REAL :: q_ad(bd%is-ng:bd%ie+ng, bd%js-ng:bd%je+ng)
8668 TYPE(FV_GRID_TYPE),
INTENT(IN),
TARGET :: gridstruct
8669 REAL,
OPTIONAL,
INTENT(IN) :: mass(bd%isd:bd%ied, bd%jsd:bd%jed)
8670 REAL,
OPTIONAL :: mass_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
8671 REAL,
INTENT(INOUT) :: fx(bd%is:bd%ie+1, bd%js:bd%je), fy(bd%is:bd%&
8672 & ie, bd%js:bd%je+1)
8673 REAL,
INTENT(INOUT) :: fx_ad(bd%is:bd%ie+1, bd%js:bd%je), fy_ad(bd%&
8674 & is:bd%ie, bd%js:bd%je+1)
8675 REAL :: fx2(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2(bd%isd:bd%ied, bd%&
8677 REAL :: fx2_ad(bd%isd:bd%ied+1, bd%jsd:bd%jed), fy2_ad(bd%isd:bd%ied&
8678 & , bd%jsd:bd%jed+1)
8679 REAL :: d2(bd%isd:bd%ied, bd%jsd:bd%jed)
8680 REAL :: d2_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
8682 INTEGER :: i, j, n, nt, i1, i2, j1, j2
8729 IF (branch .EQ. 0)
THEN 8731 CALL poprealarray(fy2, (bd%ied-bd%isd+1)*(bd%jed-bd%jsd+2))
8732 CALL poprealarray(fx2, (bd%ied-bd%isd+2)*(bd%jed-bd%jsd+1))
8738 temp_ad5 = damp2*fy2(i, j)*fy_ad(i, j)
8739 mass_ad(i, j-1) = mass_ad(i, j-1) + temp_ad5
8740 mass_ad(i, j) = mass_ad(i, j) + temp_ad5
8741 fy2_ad(i, j) = fy2_ad(i, j) + damp2*(mass(i, j-1)+mass(i, j))*&
8749 temp_ad4 = damp2*fx2(i, j)*fx_ad(i, j)
8750 mass_ad(i-1, j) = mass_ad(i-1, j) + temp_ad4
8751 mass_ad(i, j) = mass_ad(i, j) + temp_ad4
8752 fx2_ad(i, j) = fx2_ad(i, j) + damp2*(mass(i-1, j)+mass(i, j))*&
8762 fy2_ad(i, j) = fy2_ad(i, j) + fy_ad(i, j)
8769 fx2_ad(i, j) = fx2_ad(i, j) + fx_ad(i, j)
8774 IF (branch .EQ. 0)
THEN 8779 DO j=ad_to4,ad_from4,-1
8782 DO i=ad_to3,ad_from3,-1
8783 temp_ad3 = gridstruct%del6_u(i, j)*fy2_ad(i, j)
8784 d2_ad(i, j) = d2_ad(i, j) + temp_ad3
8785 d2_ad(i, j-1) = d2_ad(i, j-1) - temp_ad3
8790 & nested, bd, gridstruct%sw_corner, gridstruct%&
8791 & se_corner, gridstruct%nw_corner, gridstruct%&
8795 DO j=ad_to2,ad_from2,-1
8798 DO i=ad_to1,ad_from1,-1
8799 temp_ad2 = gridstruct%del6_v(i, j)*fx2_ad(i, j)
8800 d2_ad(i, j) = d2_ad(i, j) + temp_ad2
8801 d2_ad(i-1, j) = d2_ad(i-1, j) - temp_ad2
8806 & nested, bd, gridstruct%sw_corner, gridstruct%&
8807 & se_corner, gridstruct%nw_corner, gridstruct%&
8811 DO j=ad_to0,ad_from0,-1
8814 DO i=ad_to,ad_from,-1
8815 temp_ad1 = gridstruct%rarea(i, j)*d2_ad(i, j)
8816 fx2_ad(i, j) = fx2_ad(i, j) + temp_ad1
8817 fx2_ad(i+1, j) = fx2_ad(i+1, j) - temp_ad1
8818 fy2_ad(i, j) = fy2_ad(i, j) + temp_ad1
8819 fy2_ad(i, j+1) = fy2_ad(i, j+1) - temp_ad1
8827 DO j=je+nord+1,js-nord,-1
8828 DO i=ie+nord,is-nord,-1
8829 temp_ad0 = gridstruct%del6_u(i, j)*fy2_ad(i, j)
8830 d2_ad(i, j-1) = d2_ad(i, j-1) + temp_ad0
8831 d2_ad(i, j) = d2_ad(i, j) - temp_ad0
8837 & gridstruct%nested, bd, &
8838 & gridstruct%sw_corner, &
8839 & gridstruct%se_corner, &
8840 & gridstruct%nw_corner, &
8841 & gridstruct%ne_corner)
8842 DO j=je+nord,js-nord,-1
8843 DO i=ie+nord+1,is-nord,-1
8844 temp_ad = gridstruct%del6_v(i, j)*fx2_ad(i, j)
8845 d2_ad(i-1, j) = d2_ad(i-1, j) + temp_ad
8846 d2_ad(i, j) = d2_ad(i, j) - temp_ad
8852 & gridstruct%nested, bd, &
8853 & gridstruct%sw_corner, &
8854 & gridstruct%se_corner, &
8855 & gridstruct%nw_corner, &
8856 & gridstruct%ne_corner)
8861 IF (branch .EQ. 0)
THEN 8864 q_ad(i, j) = q_ad(i, j) + damp*d2_ad(i, j)
8871 q_ad(i, j) = q_ad(i, j) + d2_ad(i, j)
8900 & , se_corner, nw_corner, ne_corner)
8902 TYPE(FV_GRID_BOUNDS_TYPE),
INTENT(IN) :: bd
8903 INTEGER,
INTENT(IN) :: npx, npy, dir
8904 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
8905 LOGICAL,
INTENT(IN) :: nested, sw_corner, se_corner, nw_corner, &
8918 ELSE IF (dir .EQ. 1)
THEN 8935 tmp0 = q(npy-j, i-npx+1)
8947 tmp1 = q(j, 2*npx-1-i)
8959 tmp2 = q(npy-j, i-1+npx)
8968 ELSE IF (dir .EQ. 2)
THEN 8985 tmp4 = q(npy+j-1, npx-i)
8997 tmp5 = q(2*npy-1-j, i)
9009 tmp6 = q(j+1-npx, npy-i)
9045 & sw_corner, se_corner, nw_corner, ne_corner)
9047 TYPE(FV_GRID_BOUNDS_TYPE),
INTENT(IN) :: bd
9048 INTEGER,
INTENT(IN) :: npx, npy, dir
9049 REAL,
INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed)
9050 REAL,
INTENT(INOUT) :: q_ad(bd%isd:bd%ied, bd%jsd:bd%jed)
9051 LOGICAL,
INTENT(IN) :: nested, sw_corner, se_corner, nw_corner, &
9067 IF (branch .LT. 3)
THEN 9068 IF (branch .NE. 0)
THEN 9069 IF (branch .NE. 1)
THEN 9070 DO j=npy+
ng-1,npy,-1
9073 tmp_ad2 = q_ad(i, j)
9075 q_ad(npy-j, i-1+npx) = q_ad(npy-j, i-1+npx) + tmp_ad2
9080 IF (branch .EQ. 0)
THEN 9081 DO j=npy+
ng-1,npy,-1
9082 DO i=npx+
ng-1,npx,-1
9084 tmp_ad1 = q_ad(i, j)
9086 q_ad(j, 2*npx-1-i) = q_ad(j, 2*npx-1-i) + tmp_ad1
9091 IF (branch .EQ. 0)
THEN 9093 DO i=npx+
ng-1,npx,-1
9095 tmp_ad0 = q_ad(i, j)
9097 q_ad(npy-j, i-npx+1) = q_ad(npy-j, i-npx+1) + tmp_ad0
9102 IF (branch .EQ. 0)
THEN 9108 q_ad(j, 1-i) = q_ad(j, 1-i) + tmp_ad
9113 ELSE IF (branch .NE. 3)
THEN 9114 IF (branch .NE. 4)
THEN 9115 DO j=npy+
ng-1,npy,-1
9118 tmp_ad6 = q_ad(i, j)
9120 q_ad(j+1-npx, npy-i) = q_ad(j+1-npx, npy-i) + tmp_ad6
9125 IF (branch .EQ. 0)
THEN 9126 DO j=npy+
ng-1,npy,-1
9127 DO i=npx+
ng-1,npx,-1
9129 tmp_ad5 = q_ad(i, j)
9131 q_ad(2*npy-1-j, i) = q_ad(2*npy-1-j, i) + tmp_ad5
9136 IF (branch .EQ. 0)
THEN 9138 DO i=npx+
ng-1,npx,-1
9140 tmp_ad4 = q_ad(i, j)
9142 q_ad(npy+j-1, npx-i) = q_ad(npy+j-1, npx-i) + tmp_ad4
9147 IF (branch .EQ. 0)
THEN 9151 tmp_ad3 = q_ad(i, j)
9153 q_ad(1-j, i) = q_ad(1-j, i) + tmp_ad3
subroutine xppm(flux, q, c, iord, is, ie, isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
subroutine popinteger4(x)
subroutine popcontrol2b(cc)
subroutine deln_flux_fwd(nord, is, ie, js, je, npx, npy, damp, q, fx, fy, gridstruct, bd, mass)
subroutine, public pushcontrol(ctype, field)
subroutine, public fv_tp_2d_adm(q, q_ad, crx, crx_ad, cry, cry_ad, npx, npy, hord, fx, fx_ad, fy, fy_ad, xfx, xfx_ad, yfx, yfx_ad, gridstruct, bd, ra_x, ra_x_ad, ra_y, ra_y_ad, mfx, mfx_ad, mfy, mfy_ad, mass, mass_ad, nord, damp_c)
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
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)
real, parameter near_zero
subroutine deln_flux_bwd(nord, is, ie, js, je, npx, npy, damp, q, q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd, mass, mass_ad)
subroutine deln_flux_adm(nord, is, ie, js, je, npx, npy, damp, q, q_ad, fx, fx_ad, fy, fy_ad, gridstruct, bd, mass, mass_ad)
subroutine, public fv_tp_2d_fwd(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
subroutine pushcontrol1b(cc)
subroutine xppm_fwd(flux, q, c, iord, is, ie, isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
subroutine deln_flux(nord, is, ie, js, je, npx, npy, damp, q, fx, fy, gridstruct, bd, mass)
subroutine pushcontrol2b(cc)
real, parameter, public big_number
subroutine yppm_bwd(flux, flux_ad, q, q_ad, c, c_ad, jord, ifirst, ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
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)
integer, parameter, public ng
subroutine copy_corners_fwd(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
integer, parameter, public r_grid
subroutine, public fv_tp_2d_bwd(q, q_ad, crx, crx_ad, cry, cry_ad, npx, npy, hord, fx, fx_ad, fy, fy_ad, xfx, xfx_ad, yfx, yfx_ad, gridstruct, bd, ra_x, ra_x_ad, ra_y, ra_y_ad, mfx, mfx_ad, mfy, mfy_ad, mass, mass_ad, nord, damp_c)
subroutine yppm_adm(flux, flux_ad, q, q_ad, c, c_ad, jord, ifirst, ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
subroutine popcontrol3b(cc)
subroutine xppm_bwd(flux, flux_ad, q, q_ad, c, c_ad, iord, is, ie, isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
subroutine yppm(flux, q, c, jord, ifirst, ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
subroutine popcontrol1b(cc)
subroutine, public pert_ppm_adm(im, a0, al, al_ad, ar, ar_ad, iv)
real, parameter ppm_limiter
subroutine, public copy_corners_adm(q, q_ad, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
subroutine xppm_adm(flux, flux_ad, q, q_ad, c, c_ad, iord, is, ie, isd, ied, jfirst, jlast, jsd, jed, npx, npy, dxa, nested, grid_type)
subroutine copy_corners_bwd(q, q_ad, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
subroutine pushcontrol3b(cc)
subroutine yppm_fwd(flux, q, c, jord, ifirst, ilast, isd, ied, js, je, jsd, jed, npx, npy, dya, nested, grid_type)
subroutine, public popcontrol(ctype, field)
Derived type containing the data.
subroutine pushinteger4(x)
subroutine, public pert_ppm(im, a0, al, ar, iv)