12 use blockpointers,
only:
nx,
il,
ie,
ny,
jl,
je,
nz,
kl,
ke,
spectralsol, &
23 integer(kind=intType) :: i, j, k, ind, ii
24 real(kind=realtype) :: qsp, qsm, rqsp, rqsm, porvel, porflux
25 real(kind=realtype) :: pa, fs, sface, vnp, vnm
26 real(kind=realtype) :: wwx, wwy, wwz, rvol
36 #ifdef TAPENADE_REVERSE
38 do ii = 0,
il *
ny *
nz - 1
40 j = mod(ii /
il,
ny) + 2
41 k = ii / (
il *
ny) + 2
54 vnp =
w(i + 1, j, k,
ivx) *
si(i, j, k, 1) &
55 +
w(i + 1, j, k,
ivy) *
si(i, j, k, 2) &
56 +
w(i + 1, j, k,
ivz) *
si(i, j, k, 3)
57 vnm =
w(i, j, k,
ivx) *
si(i, j, k, 1) &
58 +
w(i, j, k,
ivy) *
si(i, j, k, 2) &
59 +
w(i, j, k,
ivz) *
si(i, j, k, 3)
81 porvel = porvel * porflux
86 qsp = (vnp - sface) * porvel
87 qsm = (vnm - sface) * porvel
89 rqsp = qsp *
w(i + 1, j, k,
irho)
90 rqsm = qsm *
w(i, j, k,
irho)
96 pa = porflux * (
p(i + 1, j, k) +
p(i, j, k))
111 fs = rqsp *
w(i + 1, j, k,
ivx) + rqsm *
w(i, j, k,
ivx) &
112 + pa *
si(i, j, k, 1)
113 dw(i + 1, j, k,
imx) =
dw(i + 1, j, k,
imx) - fs
116 fs = rqsp *
w(i + 1, j, k,
ivy) + rqsm *
w(i, j, k,
ivy) &
117 + pa *
si(i, j, k, 2)
118 dw(i + 1, j, k,
imy) =
dw(i + 1, j, k,
imy) - fs
121 fs = rqsp *
w(i + 1, j, k,
ivz) + rqsm *
w(i, j, k,
ivz) &
122 + pa *
si(i, j, k, 3)
123 dw(i + 1, j, k,
imz) =
dw(i + 1, j, k,
imz) - fs
126 fs = qsp *
w(i + 1, j, k,
irhoe) + qsm *
w(i, j, k,
irhoe) &
127 + porflux * (vnp *
p(i + 1, j, k) + vnm *
p(i, j, k))
130 #ifdef TAPENADE_REVERSE
146 #ifdef TAPENADE_REVERSE
148 do ii = 0,
nx *
jl *
nz - 1
150 j = mod(ii /
nx,
jl) + 1
151 k = ii / (
nx *
jl) + 2
164 vnp =
w(i, j + 1, k,
ivx) *
sj(i, j, k, 1) &
165 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 2) &
166 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 3)
167 vnm =
w(i, j, k,
ivx) *
sj(i, j, k, 1) &
168 +
w(i, j, k,
ivy) *
sj(i, j, k, 2) &
169 +
w(i, j, k,
ivz) *
sj(i, j, k, 3)
192 porvel = porvel * porflux
197 qsp = (vnp - sface) * porvel
198 qsm = (vnm - sface) * porvel
200 rqsp = qsp *
w(i, j + 1, k,
irho)
201 rqsm = qsm *
w(i, j, k,
irho)
207 pa = porflux * (
p(i, j + 1, k) +
p(i, j, k))
222 fs = rqsp *
w(i, j + 1, k,
ivx) + rqsm *
w(i, j, k,
ivx) &
223 + pa *
sj(i, j, k, 1)
224 dw(i, j + 1, k,
imx) =
dw(i, j + 1, k,
imx) - fs
227 fs = rqsp *
w(i, j + 1, k,
ivy) + rqsm *
w(i, j, k,
ivy) &
228 + pa *
sj(i, j, k, 2)
229 dw(i, j + 1, k,
imy) =
dw(i, j + 1, k,
imy) - fs
232 fs = rqsp *
w(i, j + 1, k,
ivz) + rqsm *
w(i, j, k,
ivz) &
233 + pa *
sj(i, j, k, 3)
234 dw(i, j + 1, k,
imz) =
dw(i, j + 1, k,
imz) - fs
237 fs = qsp *
w(i, j + 1, k,
irhoe) + qsm *
w(i, j, k,
irhoe) &
238 + porflux * (vnp *
p(i, j + 1, k) + vnm *
p(i, j, k))
242 #ifdef TAPENADE_REVERSE
257 #ifdef TAPENADE_REVERSE
259 do ii = 0,
nx *
ny *
kl - 1
261 j = mod(ii /
nx,
ny) + 2
262 k = ii / (
nx *
ny) + 1
275 vnp =
w(i, j, k + 1,
ivx) *
sk(i, j, k, 1) &
276 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 2) &
277 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 3)
278 vnm =
w(i, j, k,
ivx) *
sk(i, j, k, 1) &
279 +
w(i, j, k,
ivy) *
sk(i, j, k, 2) &
280 +
w(i, j, k,
ivz) *
sk(i, j, k, 3)
304 porvel = porvel * porflux
309 qsp = (vnp - sface) * porvel
310 qsm = (vnm - sface) * porvel
312 rqsp = qsp *
w(i, j, k + 1,
irho)
313 rqsm = qsm *
w(i, j, k,
irho)
319 pa = porflux * (
p(i, j, k + 1) +
p(i, j, k))
334 fs = rqsp *
w(i, j, k + 1,
ivx) + rqsm *
w(i, j, k,
ivx) &
335 + pa *
sk(i, j, k, 1)
336 dw(i, j, k + 1,
imx) =
dw(i, j, k + 1,
imx) - fs
339 fs = rqsp *
w(i, j, k + 1,
ivy) + rqsm *
w(i, j, k,
ivy) &
340 + pa *
sk(i, j, k, 2)
341 dw(i, j, k + 1,
imy) =
dw(i, j, k + 1,
imy) - fs
344 fs = rqsp *
w(i, j, k + 1,
ivz) + rqsm *
w(i, j, k,
ivz) &
345 + pa *
sk(i, j, k, 3)
346 dw(i, j, k + 1,
imz) =
dw(i, j, k + 1,
imz) - fs
349 fs = qsp *
w(i, j, k + 1,
irhoe) + qsm *
w(i, j, k,
irhoe) &
350 + porflux * (vnp *
p(i, j, k + 1) + vnm *
p(i, j, k))
353 #ifdef TAPENADE_REVERSE
383 do ii = 0,
nx *
ny *
nz - 1
385 j = mod(ii /
nx,
ny) + 2
386 k = ii / (
nx *
ny) + 2
387 rvol =
w(i, j, k,
irho) *
vol(i, j, k)
390 + rvol * (wwy *
w(i, j, k,
ivz) - wwz *
w(i, j, k,
ivy))
392 + rvol * (wwz *
w(i, j, k,
ivx) - wwx *
w(i, j, k,
ivz))
394 + rvol * (wwx *
w(i, j, k,
ivy) - wwy *
w(i, j, k,
ivx))
413 use blockpointers,
only:
nx,
ny,
nz,
il,
jl,
kl,
ie,
je,
ke,
ib,
jb,
kb, &
414 w,
p,
pori,
porj,
pork,
fw,
gamma,
si,
sj,
sk, &
427 real(kind=realtype),
parameter :: dpmax = 0.25_realtype
428 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
429 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
430 real(kind=realtype),
parameter :: omega = 0.5_realtype
431 real(kind=realtype),
parameter :: oneminomega =
one - omega
435 integer(kind=intType) :: i, j, k, ind, ii
437 real(kind=realtype) :: plim, sface
438 real(kind=realtype) :: sfil, fis2, fis4
439 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53
440 real(kind=realtype) :: ppor, rrad, dis2, dis4
441 real(kind=realtype) :: dp1, dp2, tmp, fs
442 real(kind=realtype) :: ddw1, ddw2, ddw3, ddw4, ddw5, ddw6
443 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
444 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
445 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
446 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
447 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
448 real(kind=realtype),
dimension(1:ie, 1:je, 1:ke, 3) :: dss
449 logical :: correctForK
484 #ifdef TAPENADE_REVERSE
486 do ii = 0,
ie *
je *
ke - 1
488 j = mod(ii /
ie,
je) + 1
489 k = ii / (
ie *
je) + 1
495 dss(i, j, k, 1) = abs((
p(i + 1, j, k) -
two *
p(i, j, k) +
p(i - 1, j, k)) &
496 / (omega * (
p(i + 1, j, k) +
two *
p(i, j, k) +
p(i - 1, j, k)) &
497 + oneminomega * (abs(
p(i + 1, j, k) -
p(i, j, k)) &
498 + abs(
p(i, j, k) -
p(i - 1, j, k))) + plim))
500 dss(i, j, k, 2) = abs((
p(i, j + 1, k) -
two *
p(i, j, k) +
p(i, j - 1, k)) &
501 / (omega * (
p(i, j + 1, k) +
two *
p(i, j, k) +
p(i, j - 1, k)) &
502 + oneminomega * (abs(
p(i, j + 1, k) -
p(i, j, k)) &
503 + abs(
p(i, j, k) -
p(i, j - 1, k))) + plim))
505 dss(i, j, k, 3) = abs((
p(i, j, k + 1) -
two *
p(i, j, k) +
p(i, j, k - 1)) &
506 / (omega * (
p(i, j, k + 1) +
two *
p(i, j, k) +
p(i, j, k - 1)) &
507 + oneminomega * (abs(
p(i, j, k + 1) -
p(i, j, k)) &
508 + abs(
p(i, j, k) -
p(i, j, k - 1))) + plim))
509 #ifdef TAPENADE_REVERSE
519 #ifdef TAPENADE_REVERSE
521 do ii = 0,
il *
ny *
nz - 1
523 j = mod(ii /
il,
ny) + 2
524 k = ii / (
il *
ny) + 2
534 dis2 = ppor * fis2 * min(dpmax, max(dss(i, j, k, 1), dss(i + 1, j, k, 1)))
535 dis4 =
mydim(ppor * fis4, dis2)
540 ddw1 =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
544 ddw2 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
547 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivx) &
550 ddw3 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
553 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivy) &
556 ddw4 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
559 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivz) &
571 if (correctfork)
then
572 ddw6 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
575 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
itu1) &
597 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
600 area = sqrt(
si(i, j, k, 1)**2 +
si(i, j, k, 2)**2 +
si(i, j, k, 3)**2)
601 tmp =
one / max(1.e-25_realtype, area)
602 sx =
si(i, j, k, 1) * tmp
603 sy =
si(i, j, k, 2) * tmp
604 sz =
si(i, j, k, 3) * tmp
606 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
607 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
609 unavg = uavg * sx + vavg * sy + wavg * sz
611 ova2avg =
one / a2avg
622 lam1 = abs(unavg - sface + aavg)
623 lam2 = abs(unavg - sface - aavg)
624 lam3 = abs(unavg - sface)
631 lam1 = max(lam1, epsacoustic * rrad) * area
632 lam2 = max(lam2, epsacoustic * rrad) * area
633 lam3 = max(lam3, epsshear * rrad) * area
638 abv1 =
half * (lam1 + lam2)
639 abv2 =
half * (lam1 - lam2)
642 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
643 - wavg * drw + dre) - gm53 * drk
644 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
646 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
647 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
652 fs = lam3 * dr + abv6
663 fs = lam3 * dru + uavg * abv6 + sx * abv7
664 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
669 fs = lam3 * drv + vavg * abv6 + sy * abv7
670 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
675 fs = lam3 * drw + wavg * abv6 + sz * abv7
676 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
681 fs = lam3 * dre + havg * abv6 + unavg * abv7
685 #ifdef TAPENADE_REVERSE
695 #ifdef TAPENADE_REVERSE
697 do ii = 0,
nx *
jl *
nz - 1
699 j = mod(ii /
nx,
jl) + 1
700 k = ii / (
nx *
jl) + 2
712 dis2 = ppor * fis2 * min(dpmax, max(dss(i, j, k, 2), dss(i, j + 1, k, 2)))
713 dis4 =
mydim(ppor * fis4, dis2)
718 ddw1 =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
722 ddw2 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
725 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivx) &
728 ddw3 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
731 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivy) &
734 ddw4 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
737 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivz) &
749 if (correctfork)
then
750 ddw6 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
753 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
itu1) &
775 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
778 area = sqrt(
sj(i, j, k, 1)**2 +
sj(i, j, k, 2)**2 +
sj(i, j, k, 3)**2)
779 tmp =
one / max(1.e-25_realtype, area)
780 sx =
sj(i, j, k, 1) * tmp
781 sy =
sj(i, j, k, 2) * tmp
782 sz =
sj(i, j, k, 3) * tmp
784 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
785 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
787 unavg = uavg * sx + vavg * sy + wavg * sz
789 ova2avg =
one / a2avg
800 lam1 = abs(unavg - sface + aavg)
801 lam2 = abs(unavg - sface - aavg)
802 lam3 = abs(unavg - sface)
809 lam1 = max(lam1, epsacoustic * rrad) * area
810 lam2 = max(lam2, epsacoustic * rrad) * area
811 lam3 = max(lam3, epsshear * rrad) * area
816 abv1 =
half * (lam1 + lam2)
817 abv2 =
half * (lam1 - lam2)
820 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
821 - wavg * drw + dre) - gm53 * drk
822 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
824 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
825 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
830 fs = lam3 * dr + abv6
841 fs = lam3 * dru + uavg * abv6 + sx * abv7
842 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
847 fs = lam3 * drv + vavg * abv6 + sy * abv7
848 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
853 fs = lam3 * drw + wavg * abv6 + sz * abv7
854 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
859 fs = lam3 * dre + havg * abv6 + unavg * abv7
863 #ifdef TAPENADE_REVERSE
873 #ifdef TAPENADE_REVERSE
875 do ii = 0,
nx *
ny *
kl - 1
877 j = mod(ii /
nx,
ny) + 2
878 k = ii / (
nx *
ny) + 1
889 dis2 = ppor * fis2 * min(dpmax, max(dss(i, j, k, 3), dss(i, j, k + 1, 3)))
890 dis4 =
mydim(ppor * fis4, dis2)
895 ddw1 =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
899 ddw2 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
902 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivx) &
905 ddw3 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
908 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivy) &
911 ddw4 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
914 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivz) &
926 if (correctfork)
then
927 ddw6 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
930 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
itu1) &
952 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
955 area = sqrt(
sk(i, j, k, 1)**2 +
sk(i, j, k, 2)**2 +
sk(i, j, k, 3)**2)
956 tmp =
one / max(1.e-25_realtype, area)
957 sx =
sk(i, j, k, 1) * tmp
958 sy =
sk(i, j, k, 2) * tmp
959 sz =
sk(i, j, k, 3) * tmp
961 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
962 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
964 unavg = uavg * sx + vavg * sy + wavg * sz
966 ova2avg =
one / a2avg
977 lam1 = abs(unavg - sface + aavg)
978 lam2 = abs(unavg - sface - aavg)
979 lam3 = abs(unavg - sface)
986 lam1 = max(lam1, epsacoustic * rrad) * area
987 lam2 = max(lam2, epsacoustic * rrad) * area
988 lam3 = max(lam3, epsshear * rrad) * area
993 abv1 =
half * (lam1 + lam2)
994 abv2 =
half * (lam1 - lam2)
997 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
998 - wavg * drw + dre) - gm53 * drk
999 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
1001 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
1002 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
1007 fs = lam3 * dr + abv6
1010 #ifndef USE_TAPENADE
1018 fs = lam3 * dru + uavg * abv6 + sx * abv7
1019 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
1024 fs = lam3 * drv + vavg * abv6 + sy * abv7
1025 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
1030 fs = lam3 * drw + wavg * abv6 + sz * abv7
1031 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
1036 fs = lam3 * dre + havg * abv6 + unavg * abv7
1040 #ifdef TAPENADE_REVERSE
1057 use blockpointers,
only:
nx,
ny,
nz,
il,
jl,
kl,
ie,
je,
ke,
ib,
jb,
kb, &
1069 real(kind=realtype),
parameter :: dssmax = 0.25_realtype
1073 integer(kind=intType) :: i, j, k, ind, ii
1075 real(kind=realtype) :: sslim, rhoi
1076 real(kind=realtype) :: sfil, fis2, fis4
1077 real(kind=realtype) :: ppor, rrad, dis2, dis4
1078 real(kind=realtype) :: ddw1, ddw2, ddw3, ddw4, ddw5, fs
1079 real(kind=realtype),
dimension(1:ie, 1:je, 1:ke, 3) :: dss
1080 real(kind=realtype),
dimension(0:ib, 0:jb, 0:kb) :: ss
1119 #ifdef TAPENADE_REVERSE
1121 do ii = 0, (
ib + 1) * (
jb + 1) * (
kb + 1) - 1
1123 j = mod(ii / (
ib + 1),
jb + 1)
1124 k = ii / ((
ib + 1) * (
jb + 1))
1130 ss(i, j, k) =
p(i, j, k) / (
w(i, j, k,
irho)**
gamma(i, j, k))
1131 #ifdef TAPENADE_REVERSE
1141 #ifdef TAPENADE_REVERSE
1143 do ii = 0,
ie *
je *
ke - 1
1145 j = mod(ii /
ie,
je) + 1
1146 k = ii / (
ie *
je) + 1
1152 dss(i, j, k, 1) = abs((ss(i + 1, j, k) -
two * ss(i, j, k) + ss(i - 1, j, k)) &
1153 / (ss(i + 1, j, k) +
two * ss(i, j, k) + ss(i - 1, j, k) + sslim))
1155 dss(i, j, k, 2) = abs((ss(i, j + 1, k) -
two * ss(i, j, k) + ss(i, j - 1, k)) &
1156 / (ss(i, j + 1, k) +
two * ss(i, j, k) + ss(i, j - 1, k) + sslim))
1158 dss(i, j, k, 3) = abs((ss(i, j, k + 1) -
two * ss(i, j, k) + ss(i, j, k - 1)) &
1159 / (ss(i, j, k + 1) +
two * ss(i, j, k) + ss(i, j, k - 1) + sslim))
1160 #ifdef TAPENADE_REVERSE
1197 #ifdef TAPENADE_REVERSE
1199 do ii = 0,
il *
ny *
nz - 1
1201 j = mod(ii /
il,
ny) + 2
1202 k = ii / (
il *
ny) + 2
1212 rrad = ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
1214 dis2 = fis2 * rrad * min(dssmax, max(dss(i, j, k, 1), dss(i + 1, j, k, 1)))
1215 dis4 =
mydim(fis4 * rrad, dis2)
1221 ddw1 =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
1230 ddw2 =
w(i + 1, j, k,
ivx) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
1232 - dis4 * (
w(i + 2, j, k,
ivx) *
w(i + 2, j, k,
irho) - &
1235 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
1240 ddw3 =
w(i + 1, j, k,
ivy) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
1242 - dis4 * (
w(i + 2, j, k,
ivy) *
w(i + 2, j, k,
irho) - &
1245 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
1250 ddw4 =
w(i + 1, j, k,
ivz) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
1252 - dis4 * (
w(i + 2, j, k,
ivz) *
w(i + 2, j, k,
irho) - &
1255 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
1260 ddw5 = (
w(i + 1, j, k,
irhoe) +
p(i + 1, j, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
1262 - dis4 * ((
w(i + 2, j, k,
irhoe) +
p(i + 2, j, k)) - &
1263 (
w(i - 1, j, k,
irhoe) +
p(i - 1, j, k)) -
three * ddw5)
1267 #ifdef TAPENADE_REVERSE
1278 #ifdef TAPENADE_REVERSE
1280 do ii = 0,
nx *
jl *
nz - 1
1282 j = mod(ii /
nx,
jl) + 1
1283 k = ii / (
nx *
jl) + 2
1293 rrad = ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
1295 dis2 = fis2 * rrad * min(dssmax, max(dss(i, j, k, 2), dss(i, j + 1, k, 2)))
1296 dis4 =
mydim(fis4 * rrad, dis2)
1302 ddw1 =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
1311 ddw2 =
w(i, j + 1, k,
ivx) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
1313 - dis4 * (
w(i, j + 2, k,
ivx) *
w(i, j + 2, k,
irho) - &
1316 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
1321 ddw3 =
w(i, j + 1, k,
ivy) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
1323 - dis4 * (
w(i, j + 2, k,
ivy) *
w(i, j + 2, k,
irho) - &
1326 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
1331 ddw4 =
w(i, j + 1, k,
ivz) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
1333 - dis4 * (
w(i, j + 2, k,
ivz) *
w(i, j + 2, k,
irho) - &
1336 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
1341 ddw5 = (
w(i, j + 1, k,
irhoe) +
p(i, j + 1, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
1343 - dis4 * ((
w(i, j + 2, k,
irhoe) +
p(i, j + 2, k)) - &
1344 (
w(i, j - 1, k,
irhoe) +
p(i, j - 1, k)) -
three * ddw5)
1348 #ifdef TAPENADE_REVERSE
1358 #ifdef TAPENADE_REVERSE
1360 do ii = 0,
nx *
ny *
kl - 1
1362 j = mod(ii /
nx,
ny) + 2
1363 k = ii / (
nx *
ny) + 1
1373 rrad = ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
1375 dis2 = fis2 * rrad * min(dssmax, max(dss(i, j, k, 3), dss(i, j, k + 1, 3)))
1376 dis4 =
mydim(fis4 * rrad, dis2)
1382 ddw1 =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
1391 ddw2 =
w(i, j, k + 1,
ivx) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
1393 - dis4 * (
w(i, j, k + 2,
ivx) *
w(i, j, k + 2,
irho) - &
1396 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
1401 ddw3 =
w(i, j, k + 1,
ivy) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
1403 - dis4 * (
w(i, j, k + 2,
ivy) *
w(i, j, k + 2,
irho) - &
1406 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
1411 ddw4 =
w(i, j, k + 1,
ivz) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
1413 - dis4 * (
w(i, j, k + 2,
ivz) *
w(i, j, k + 2,
irho) - &
1416 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
1421 ddw5 = (
w(i, j, k + 1,
irhoe) +
p(i, j, k + 1)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
1423 - dis4 * ((
w(i, j, k + 2,
irhoe) +
p(i, j, k + 2)) - &
1424 (
w(i, j, k - 1,
irhoe) +
p(i, j, k - 1)) -
three * ddw5)
1429 #ifdef TAPENADE_REVERSE
1451 use blockpointers,
only:
il,
jl,
kl,
ie,
je,
ke,
ib,
jb,
kb,
w,
p, &
1468 logical,
intent(in) :: fineGrid
1472 integer(kind=porType) :: por
1474 integer(kind=intType) :: nwInt
1475 integer(kind=intType) :: i, j, k, ind
1476 integer(kind=intType) :: limUsed, riemannUsed
1478 real(kind=realtype) :: sx, sy, sz, omk, opk, sfil, gammaface
1479 real(kind=realtype) :: factminmod, sface
1481 real(kind=realtype),
dimension(nw) :: left, right
1482 real(kind=realtype),
dimension(nw) :: du1, du2, du3
1483 real(kind=realtype),
dimension(nwf) :: flux
1485 logical :: firstOrderK, correctForK, rotationalPeriodic
1496 rotationalperiodic = .true.
1498 rotationalperiodic = .false.
1533 if (finegrid) limused =
limiter
1543 if (finegrid) riemannused =
riemann
1561 if (correctfork)
then
1564 firstorderk = .true.
1567 firstorderk = .false.
1571 firstorderk = .false.
1592 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
1603 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
1606 right(
ivx) =
w(i + 1, j, k,
ivx)
1607 right(
ivy) =
w(i + 1, j, k,
ivy)
1608 right(
ivz) =
w(i + 1, j, k,
ivz)
1609 right(
irhoe) =
p(i + 1, j, k)
1610 if (correctfork) right(
itu1) =
w(i + 1, j, k,
itu1)
1637 #ifndef USE_TAPENADE
1656 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
1667 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
1670 right(
ivx) =
w(i, j + 1, k,
ivx)
1671 right(
ivy) =
w(i, j + 1, k,
ivy)
1672 right(
ivz) =
w(i, j + 1, k,
ivz)
1673 right(
irhoe) =
p(i, j + 1, k)
1674 if (correctfork) right(
itu1) =
w(i, j + 1, k,
itu1)
1701 #ifndef USE_TAPENADE
1720 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
1731 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
1734 right(
ivx) =
w(i, j, k + 1,
ivx)
1735 right(
ivy) =
w(i, j, k + 1,
ivy)
1736 right(
ivz) =
w(i, j, k + 1,
ivz)
1737 right(
irhoe) =
p(i, j, k + 1)
1738 if (correctfork) right(
itu1) =
w(i, j, k + 1,
itu1)
1765 #ifndef USE_TAPENADE
1803 du3(
ivx) =
w(i + 2, j, k,
ivx) -
w(i + 1, j, k,
ivx)
1807 du3(
ivy) =
w(i + 2, j, k,
ivy) -
w(i + 1, j, k,
ivy)
1811 du3(
ivz) =
w(i + 2, j, k,
ivz) -
w(i + 1, j, k,
ivz)
1813 du1(
irhoe) =
p(i, j, k) -
p(i - 1, j, k)
1814 du2(
irhoe) =
p(i + 1, j, k) -
p(i, j, k)
1815 du3(
irhoe) =
p(i + 2, j, k) -
p(i + 1, j, k)
1817 if (correctfork)
then
1839 right(
ivx) = right(
ivx) +
w(i + 1, j, k,
ivx)
1840 right(
ivy) = right(
ivy) +
w(i + 1, j, k,
ivy)
1841 right(
ivz) = right(
ivz) +
w(i + 1, j, k,
ivz)
1844 if (correctfork)
then
1852 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
1881 #ifndef USE_TAPENADE
1906 du3(
ivx) =
w(i, j + 2, k,
ivx) -
w(i, j + 1, k,
ivx)
1910 du3(
ivy) =
w(i, j + 2, k,
ivy) -
w(i, j + 1, k,
ivy)
1914 du3(
ivz) =
w(i, j + 2, k,
ivz) -
w(i, j + 1, k,
ivz)
1916 du1(
irhoe) =
p(i, j, k) -
p(i, j - 1, k)
1917 du2(
irhoe) =
p(i, j + 1, k) -
p(i, j, k)
1918 du3(
irhoe) =
p(i, j + 2, k) -
p(i, j + 1, k)
1920 if (correctfork)
then
1942 right(
ivx) = right(
ivx) +
w(i, j + 1, k,
ivx)
1943 right(
ivy) = right(
ivy) +
w(i, j + 1, k,
ivy)
1944 right(
ivz) = right(
ivz) +
w(i, j + 1, k,
ivz)
1947 if (correctfork)
then
1955 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
1984 #ifndef USE_TAPENADE
2009 du3(
ivx) =
w(i, j, k + 2,
ivx) -
w(i, j, k + 1,
ivx)
2013 du3(
ivy) =
w(i, j, k + 2,
ivy) -
w(i, j, k + 1,
ivy)
2017 du3(
ivz) =
w(i, j, k + 2,
ivz) -
w(i, j, k + 1,
ivz)
2019 du1(
irhoe) =
p(i, j, k) -
p(i, j, k - 1)
2020 du2(
irhoe) =
p(i, j, k + 1) -
p(i, j, k)
2021 du3(
irhoe) =
p(i, j, k + 2) -
p(i, j, k + 1)
2023 if (correctfork)
then
2045 right(
ivx) = right(
ivx) +
w(i, j, k + 1,
ivx)
2046 right(
ivy) = right(
ivy) +
w(i, j, k + 1,
ivy)
2047 right(
ivz) = right(
ivz) +
w(i, j, k + 1,
ivz)
2050 if (correctfork)
then
2058 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
2087 #ifndef USE_TAPENADE
2116 real(kind=realtype),
parameter :: epslim = 1.e-10_realtype
2120 real(kind=realtype),
dimension(:),
intent(inout) :: du1, du2, du3
2121 real(kind=realtype),
dimension(:),
intent(out) :: left, right
2123 real(kind=realtype),
dimension(:, :, :, :, :),
pointer :: rotmatrix
2127 integer(kind=intType) :: l
2129 real(kind=realtype) :: rl1, rl2, rr1, rr2, tmp, dvx, dvy, dvz
2131 real(kind=realtype),
dimension(3, 3) :: rot
2136 if (rotationalperiodic)
then
2141 rot(1, 1) = rotmatrix(i, j, k, 1, 1)
2142 rot(1, 2) = rotmatrix(i, j, k, 1, 2)
2143 rot(1, 3) = rotmatrix(i, j, k, 1, 3)
2145 rot(2, 1) = rotmatrix(i, j, k, 2, 1)
2146 rot(2, 2) = rotmatrix(i, j, k, 2, 2)
2147 rot(2, 3) = rotmatrix(i, j, k, 2, 3)
2149 rot(3, 1) = rotmatrix(i, j, k, 3, 1)
2150 rot(3, 2) = rotmatrix(i, j, k, 3, 2)
2151 rot(3, 3) = rotmatrix(i, j, k, 3, 3)
2156 dvx = du1(
ivx); dvy = du1(
ivy); dvz = du1(
ivz)
2157 du1(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
2158 du1(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
2159 du1(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
2161 dvx = du2(
ivx); dvy = du2(
ivy); dvz = du2(
ivz)
2162 du2(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
2163 du2(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
2164 du2(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
2166 dvx = du3(
ivx); dvy = du3(
ivy); dvz = du3(
ivz)
2167 du3(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
2168 du3(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
2169 du3(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
2175 select case (limused)
2183 left(l) = omk * du1(l) + opk * du2(l)
2184 right(l) = -omk * du3(l) - opk * du2(l)
2199 tmp =
one / sign(max(abs(du2(l)), epslim), du2(l))
2201 du2(l) / sign(max(abs(du1(l)), epslim), du1(l)))
2202 rl2 = max(
zero, du1(l) * tmp)
2204 rr1 = max(
zero, du3(l) * tmp)
2206 du2(l) / sign(max(abs(du3(l)), epslim), du3(l)))
2210 rl1 = rl1 * (rl1 +
one) / (rl1 * rl1 +
one)
2211 rl2 = rl2 * (rl2 +
one) / (rl2 * rl2 +
one)
2212 rr1 = rr1 * (rr1 +
one) / (rr1 * rr1 +
one)
2213 rr2 = rr2 * (rr2 +
one) / (rr2 * rr2 +
one)
2218 left(l) = omk * rl1 * du1(l) + opk * rl2 * du2(l)
2219 right(l) = -opk * rr1 * du2(l) - omk * rr2 * du3(l)
2235 tmp =
one / sign(max(abs(du2(l)), epslim), du2(l))
2237 du2(l) / sign(max(abs(du1(l)), epslim), du1(l)))
2238 rl2 = max(
zero, du1(l) * tmp)
2240 rr1 = max(
zero, du3(l) * tmp)
2242 du2(l) / sign(max(abs(du3(l)), epslim), du3(l)))
2246 rl1 = min(
one, factminmod * rl1)
2247 rl2 = min(
one, factminmod * rl2)
2248 rr1 = min(
one, factminmod * rr1)
2249 rr2 = min(
one, factminmod * rr2)
2254 left(l) = omk * rl1 * du1(l) + opk * rl2 * du2(l)
2255 right(l) = -opk * rr1 * du2(l) - omk * rr2 * du3(l)
2265 if (firstorderk)
then
2274 if (rotationalperiodic)
then
2278 dvx = left(
ivx); dvy = left(
ivy); dvz = left(
ivz)
2279 left(
ivx) = rot(1, 1) * dvx + rot(2, 1) * dvy + rot(3, 1) * dvz
2280 left(
ivy) = rot(1, 2) * dvx + rot(2, 2) * dvy + rot(3, 2) * dvz
2281 left(
ivz) = rot(1, 3) * dvx + rot(2, 3) * dvy + rot(3, 3) * dvz
2285 dvx = right(
ivx); dvy = right(
ivy); dvz = right(
ivz)
2286 right(
ivx) = rot(1, 1) * dvx + rot(2, 1) * dvy + rot(3, 1) * dvz
2287 right(
ivy) = rot(1, 2) * dvx + rot(2, 2) * dvy + rot(3, 2) * dvz
2288 right(
ivz) = rot(1, 3) * dvx + rot(2, 3) * dvy + rot(3, 3) * dvz
2305 real(kind=realtype),
dimension(*),
intent(in) :: left, right
2306 real(kind=realtype),
dimension(*),
intent(out) :: flux
2310 real(kind=realtype) :: porflux, rface
2311 real(kind=realtype) :: etl, etr, z1l, z1r, tmp
2312 real(kind=realtype) :: dr, dru, drv, drw, dre, drk
2313 real(kind=realtype) :: ravg, uavg, vavg, wavg, havg, kavg
2314 real(kind=realtype) :: alphaavg, a2avg, aavg, unavg
2315 real(kind=realtype) :: ovaavg, ova2avg, area, eta
2316 real(kind=realtype) :: gm1, gm53
2317 real(kind=realtype) :: lam1, lam2, lam3
2318 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
2319 real(kind=realtype),
dimension(2) :: ktmp
2329 gm1 = gammaface -
one
2334 select case (riemannused)
2350 z1l = sqrt(left(
irho))
2351 z1r = sqrt(right(
irho))
2352 tmp =
one / (z1l + z1r)
2357 if (correctfork)
then
2362 ktmp(1) = left(
itu1)
2363 ktmp(2) = right(
itu1)
2373 kavg = tmp * (z1l * left(
itu1) + z1r * right(
itu1))
2388 left(
irhoe), ktmp(1), etl, correctfork)
2391 right(
irhoe), ktmp(2), etr, correctfork)
2405 ravg =
fourth * (z1r + z1l)**2
2406 uavg = tmp * (z1l * left(
ivx) + z1r * right(
ivx))
2407 vavg = tmp * (z1l * left(
ivy) + z1r * right(
ivy))
2408 wavg = tmp * (z1l * left(
ivz) + z1r * right(
ivz))
2409 havg = tmp * ((etl + left(
irhoe)) / z1l &
2410 + (etr + right(
irhoe)) / z1r)
2415 area = sqrt(sx**2 + sy**2 + sz**2)
2416 tmp =
one / max(1.e-25_realtype, area)
2425 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2426 a2avg = abs(gm1 * (havg - alphaavg) - gm53 * kavg)
2428 unavg = uavg * sx + vavg * sy + wavg * sz
2431 ova2avg =
one / a2avg
2449 eta =
half * (abs((left(
ivx) - right(
ivx)) * sx &
2450 + (left(
ivy) - right(
ivy)) * sy &
2451 + (left(
ivz) - right(
ivz)) * sz) &
2452 + abs(sqrt(gammaface * left(
irhoe) / left(
irho)) &
2453 - sqrt(gammaface * right(
irhoe) / right(
irho))))
2457 lam1 = abs(unavg - rface + aavg)
2458 lam2 = abs(unavg - rface - aavg)
2459 lam3 = abs(unavg - rface)
2464 if (lam1 < tmp) lam1 = eta +
fourth * lam1 * lam1 / eta
2465 if (lam2 < tmp) lam2 = eta +
fourth * lam2 * lam2 / eta
2466 if (lam3 < tmp) lam3 = eta +
fourth * lam3 * lam3 / eta
2478 abv1 =
half * (lam1 + lam2)
2479 abv2 =
half * (lam1 - lam2)
2482 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2483 - wavg * drw + dre) - gm53 * drk
2484 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2486 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2487 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2493 flux(
irho) = -porflux * (lam3 * dr + abv6)
2494 flux(
imx) = -porflux * (lam3 * dru + uavg * abv6 &
2496 flux(
imy) = -porflux * (lam3 * drv + vavg * abv6 &
2498 flux(
imz) = -porflux * (lam3 * drw + wavg * abv6 &
2500 flux(
irhoe) = -porflux * (lam3 * dre + havg * abv6 &
2514 "Turkel preconditioner not implemented yet")
2518 "choi merkle preconditioner not implemented yet")
2523 call terminate(
"riemannFlux",
"van leer fvs not implemented yet")
2526 call terminate(
"riemannFlux",
"ausmdv fvs not implemented yet")
2546 #ifndef USE_TAPENADE
2553 real(kind=realtype),
parameter :: twothird =
two *
third
2554 real(kind=realtype),
parameter :: xminn = 1.e-14_realtype
2558 integer(kind=intType) :: i, j, k, ii
2560 real(kind=realtype) :: rfilv, por, mul, mue, mut, heatcoef
2561 real(kind=realtype) :: gm1, factlamheat, factturbheat
2562 real(kind=realtype) :: u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z
2563 real(kind=realtype) :: q_x, q_y, q_z, ubar, vbar, wbar
2564 real(kind=realtype) :: corr, ssx, ssy, ssz, ss, fracdiv
2565 real(kind=realtype) :: tauxx, tauyy, tauzz
2566 real(kind=realtype) :: tauxy, tauxz, tauyz
2567 real(kind=realtype) :: tauxxs, tauyys, tauzzs
2568 real(kind=realtype) :: tauxys, tauxzs, tauyzs
2569 real(kind=realtype) :: exx, eyy, ezz
2570 real(kind=realtype) :: exy, exz, eyz
2571 real(kind=realtype) :: wxy, wxz, wyz, wyx, wzx, wzy
2572 real(kind=realtype) :: den, ccr1, fact
2573 real(kind=realtype) :: fmx, fmy, fmz, frhoe
2574 logical :: correctForK, storeWallTensor
2590 storewalltensor = .false.
2592 storewalltensor = .true.
2594 storewalltensor = .true.
2603 #ifdef TAPENADE_REVERSE
2605 do ii = 0,
nx *
ny *
kl - 1
2607 j = mod(ii /
nx,
ny) + 2
2608 k = ii / (
nx *
ny) + 1
2627 mul = por * (
rlv(i, j, k) +
rlv(i, j, k + 1))
2635 heatcoef = mul * factlamheat + mue * factturbheat
2640 u_x =
fourth * (
ux(i - 1, j - 1, k) +
ux(i, j - 1, k) &
2641 +
ux(i - 1, j, k) +
ux(i, j, k))
2642 u_y =
fourth * (
uy(i - 1, j - 1, k) +
uy(i, j - 1, k) &
2643 +
uy(i - 1, j, k) +
uy(i, j, k))
2644 u_z =
fourth * (
uz(i - 1, j - 1, k) +
uz(i, j - 1, k) &
2645 +
uz(i - 1, j, k) +
uz(i, j, k))
2647 v_x =
fourth * (
vx(i - 1, j - 1, k) +
vx(i, j - 1, k) &
2648 +
vx(i - 1, j, k) +
vx(i, j, k))
2649 v_y =
fourth * (
vy(i - 1, j - 1, k) +
vy(i, j - 1, k) &
2650 +
vy(i - 1, j, k) +
vy(i, j, k))
2651 v_z =
fourth * (
vz(i - 1, j - 1, k) +
vz(i, j - 1, k) &
2652 +
vz(i - 1, j, k) +
vz(i, j, k))
2654 w_x =
fourth * (
wx(i - 1, j - 1, k) +
wx(i, j - 1, k) &
2655 +
wx(i - 1, j, k) +
wx(i, j, k))
2656 w_y =
fourth * (
wy(i - 1, j - 1, k) +
wy(i, j - 1, k) &
2657 +
wy(i - 1, j, k) +
wy(i, j, k))
2658 w_z =
fourth * (
wz(i - 1, j - 1, k) +
wz(i, j - 1, k) &
2659 +
wz(i - 1, j, k) +
wz(i, j, k))
2661 q_x =
fourth * (
qx(i - 1, j - 1, k) +
qx(i, j - 1, k) &
2662 +
qx(i - 1, j, k) +
qx(i, j, k))
2663 q_y =
fourth * (
qy(i - 1, j - 1, k) +
qy(i, j - 1, k) &
2664 +
qy(i - 1, j, k) +
qy(i, j, k))
2665 q_z =
fourth * (
qz(i - 1, j - 1, k) +
qz(i, j - 1, k) &
2666 +
qz(i - 1, j, k) +
qz(i, j, k))
2673 ssx =
eighth * (
x(i - 1, j - 1, k + 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
2674 +
x(i - 1, j, k + 1, 1) -
x(i - 1, j, k - 1, 1) &
2675 +
x(i, j - 1, k + 1, 1) -
x(i, j - 1, k - 1, 1) &
2676 +
x(i, j, k + 1, 1) -
x(i, j, k - 1, 1))
2677 ssy =
eighth * (
x(i - 1, j - 1, k + 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
2678 +
x(i - 1, j, k + 1, 2) -
x(i - 1, j, k - 1, 2) &
2679 +
x(i, j - 1, k + 1, 2) -
x(i, j - 1, k - 1, 2) &
2680 +
x(i, j, k + 1, 2) -
x(i, j, k - 1, 2))
2681 ssz =
eighth * (
x(i - 1, j - 1, k + 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
2682 +
x(i - 1, j, k + 1, 3) -
x(i - 1, j, k - 1, 3) &
2683 +
x(i, j - 1, k + 1, 3) -
x(i, j - 1, k - 1, 3) &
2684 +
x(i, j, k + 1, 3) -
x(i, j, k - 1, 3))
2689 ss =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
2696 corr = u_x * ssx + u_y * ssy + u_z * ssz &
2697 - (
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)) * ss
2698 u_x = u_x - corr * ssx
2699 u_y = u_y - corr * ssy
2700 u_z = u_z - corr * ssz
2702 corr = v_x * ssx + v_y * ssy + v_z * ssz &
2703 - (
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)) * ss
2704 v_x = v_x - corr * ssx
2705 v_y = v_y - corr * ssy
2706 v_z = v_z - corr * ssz
2708 corr = w_x * ssx + w_y * ssy + w_z * ssz &
2709 - (
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)) * ss
2710 w_x = w_x - corr * ssx
2711 w_y = w_y - corr * ssy
2712 w_z = w_z - corr * ssz
2714 corr = q_x * ssx + q_y * ssy + q_z * ssz &
2715 + (
aa(i, j, k + 1) -
aa(i, j, k)) * ss
2716 q_x = q_x - corr * ssx
2717 q_y = q_y - corr * ssy
2718 q_z = q_z - corr * ssz
2727 fracdiv = twothird * (u_x + v_y + w_z)
2729 tauxxs =
two * u_x - fracdiv
2730 tauyys =
two * v_y - fracdiv
2731 tauzzs =
two * w_z - fracdiv
2737 q_x = heatcoef * q_x
2738 q_y = heatcoef * q_y
2739 q_z = heatcoef * q_z
2759 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
2760 v_x * v_x + v_y * v_y + v_z * v_z + &
2761 w_x * w_x + w_y * w_y + w_z * w_z)
2765 den = max(den, xminn)
2770 fact = mue * ccr1 / den
2782 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
2783 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
2784 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
2786 exy = fact * (wxy * tauyys + wxz * tauyzs + &
2787 wyx * tauxxs + wyz * tauxzs)
2788 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
2789 wzx * tauxxs + wzy * tauxys)
2790 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
2791 wzx * tauxys + wzy * tauyys)
2794 tauxx = mut * tauxxs - exx
2795 tauyy = mut * tauyys - eyy
2796 tauzz = mut * tauzzs - ezz
2797 tauxy = mut * tauxys - exy
2798 tauxz = mut * tauxzs - exz
2799 tauyz = mut * tauyzs - eyz
2804 tauxx = mut * tauxxs
2805 tauyy = mut * tauyys
2806 tauzz = mut * tauzzs
2807 tauxy = mut * tauxys
2808 tauxz = mut * tauxzs
2809 tauyz = mut * tauyzs
2822 fmx = tauxx *
sk(i, j, k, 1) + tauxy *
sk(i, j, k, 2) &
2823 + tauxz *
sk(i, j, k, 3)
2824 fmy = tauxy *
sk(i, j, k, 1) + tauyy *
sk(i, j, k, 2) &
2825 + tauyz *
sk(i, j, k, 3)
2826 fmz = tauxz *
sk(i, j, k, 1) + tauyz *
sk(i, j, k, 2) &
2827 + tauzz *
sk(i, j, k, 3)
2828 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sk(i, j, k, 1)
2829 frhoe = frhoe + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sk(i, j, k, 2)
2830 frhoe = frhoe + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sk(i, j, k, 3)
2831 frhoe = frhoe - q_x *
sk(i, j, k, 1) - q_y *
sk(i, j, k, 2) - q_z *
sk(i, j, k, 3)
2840 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fmx
2841 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fmy
2842 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fmz
2849 if (k == 1 .and. storewalltensor .and. &
2867 if (k ==
kl .and. storewalltensor .and. &
2880 #ifdef TAPENADE_REVERSE
2896 #ifdef TAPENADE_REVERSE
2898 do ii = 0,
nx *
jl *
nz - 1
2900 j = mod(ii /
nx,
jl) + 1
2901 k = ii / (
nx *
jl) + 2
2920 mul = por * (
rlv(i, j, k) +
rlv(i, j + 1, k))
2928 heatcoef = mul * factlamheat + mue * factturbheat
2933 u_x =
fourth * (
ux(i - 1, j, k - 1) +
ux(i, j, k - 1) &
2934 +
ux(i - 1, j, k) +
ux(i, j, k))
2935 u_y =
fourth * (
uy(i - 1, j, k - 1) +
uy(i, j, k - 1) &
2936 +
uy(i - 1, j, k) +
uy(i, j, k))
2937 u_z =
fourth * (
uz(i - 1, j, k - 1) +
uz(i, j, k - 1) &
2938 +
uz(i - 1, j, k) +
uz(i, j, k))
2940 v_x =
fourth * (
vx(i - 1, j, k - 1) +
vx(i, j, k - 1) &
2941 +
vx(i - 1, j, k) +
vx(i, j, k))
2942 v_y =
fourth * (
vy(i - 1, j, k - 1) +
vy(i, j, k - 1) &
2943 +
vy(i - 1, j, k) +
vy(i, j, k))
2944 v_z =
fourth * (
vz(i - 1, j, k - 1) +
vz(i, j, k - 1) &
2945 +
vz(i - 1, j, k) +
vz(i, j, k))
2947 w_x =
fourth * (
wx(i - 1, j, k - 1) +
wx(i, j, k - 1) &
2948 +
wx(i - 1, j, k) +
wx(i, j, k))
2949 w_y =
fourth * (
wy(i - 1, j, k - 1) +
wy(i, j, k - 1) &
2950 +
wy(i - 1, j, k) +
wy(i, j, k))
2951 w_z =
fourth * (
wz(i - 1, j, k - 1) +
wz(i, j, k - 1) &
2952 +
wz(i - 1, j, k) +
wz(i, j, k))
2954 q_x =
fourth * (
qx(i - 1, j, k - 1) +
qx(i, j, k - 1) &
2955 +
qx(i - 1, j, k) +
qx(i, j, k))
2956 q_y =
fourth * (
qy(i - 1, j, k - 1) +
qy(i, j, k - 1) &
2957 +
qy(i - 1, j, k) +
qy(i, j, k))
2958 q_z =
fourth * (
qz(i - 1, j, k - 1) +
qz(i, j, k - 1) &
2959 +
qz(i - 1, j, k) +
qz(i, j, k))
2966 ssx =
eighth * (
x(i - 1, j + 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
2967 +
x(i - 1, j + 1, k, 1) -
x(i - 1, j - 1, k, 1) &
2968 +
x(i, j + 1, k - 1, 1) -
x(i, j - 1, k - 1, 1) &
2969 +
x(i, j + 1, k, 1) -
x(i, j - 1, k, 1))
2970 ssy =
eighth * (
x(i - 1, j + 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
2971 +
x(i - 1, j + 1, k, 2) -
x(i - 1, j - 1, k, 2) &
2972 +
x(i, j + 1, k - 1, 2) -
x(i, j - 1, k - 1, 2) &
2973 +
x(i, j + 1, k, 2) -
x(i, j - 1, k, 2))
2974 ssz =
eighth * (
x(i - 1, j + 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
2975 +
x(i - 1, j + 1, k, 3) -
x(i - 1, j - 1, k, 3) &
2976 +
x(i, j + 1, k - 1, 3) -
x(i, j - 1, k - 1, 3) &
2977 +
x(i, j + 1, k, 3) -
x(i, j - 1, k, 3))
2982 ss =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
2989 corr = u_x * ssx + u_y * ssy + u_z * ssz &
2990 - (
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)) * ss
2991 u_x = u_x - corr * ssx
2992 u_y = u_y - corr * ssy
2993 u_z = u_z - corr * ssz
2995 corr = v_x * ssx + v_y * ssy + v_z * ssz &
2996 - (
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)) * ss
2997 v_x = v_x - corr * ssx
2998 v_y = v_y - corr * ssy
2999 v_z = v_z - corr * ssz
3001 corr = w_x * ssx + w_y * ssy + w_z * ssz &
3002 - (
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)) * ss
3003 w_x = w_x - corr * ssx
3004 w_y = w_y - corr * ssy
3005 w_z = w_z - corr * ssz
3007 corr = q_x * ssx + q_y * ssy + q_z * ssz &
3008 + (
aa(i, j + 1, k) -
aa(i, j, k)) * ss
3009 q_x = q_x - corr * ssx
3010 q_y = q_y - corr * ssy
3011 q_z = q_z - corr * ssz
3020 fracdiv = twothird * (u_x + v_y + w_z)
3022 tauxxs =
two * u_x - fracdiv
3023 tauyys =
two * v_y - fracdiv
3024 tauzzs =
two * w_z - fracdiv
3030 q_x = heatcoef * q_x
3031 q_y = heatcoef * q_y
3032 q_z = heatcoef * q_z
3052 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
3053 v_x * v_x + v_y * v_y + v_z * v_z + &
3054 w_x * w_x + w_y * w_y + w_z * w_z)
3058 den = max(den, xminn)
3063 fact = mue * ccr1 / den
3075 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
3076 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
3077 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
3079 exy = fact * (wxy * tauyys + wxz * tauyzs + &
3080 wyx * tauxxs + wyz * tauxzs)
3081 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
3082 wzx * tauxxs + wzy * tauxys)
3083 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
3084 wzx * tauxys + wzy * tauyys)
3087 tauxx = mut * tauxxs - exx
3088 tauyy = mut * tauyys - eyy
3089 tauzz = mut * tauzzs - ezz
3090 tauxy = mut * tauxys - exy
3091 tauxz = mut * tauxzs - exz
3092 tauyz = mut * tauyzs - eyz
3097 tauxx = mut * tauxxs
3098 tauyy = mut * tauyys
3099 tauzz = mut * tauzzs
3100 tauxy = mut * tauxys
3101 tauxz = mut * tauxzs
3102 tauyz = mut * tauyzs
3115 fmx = tauxx *
sj(i, j, k, 1) + tauxy *
sj(i, j, k, 2) &
3116 + tauxz *
sj(i, j, k, 3)
3117 fmy = tauxy *
sj(i, j, k, 1) + tauyy *
sj(i, j, k, 2) &
3118 + tauyz *
sj(i, j, k, 3)
3119 fmz = tauxz *
sj(i, j, k, 1) + tauyz *
sj(i, j, k, 2) &
3120 + tauzz *
sj(i, j, k, 3)
3121 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sj(i, j, k, 1) &
3122 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sj(i, j, k, 2) &
3123 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sj(i, j, k, 3) &
3124 - q_x *
sj(i, j, k, 1) - q_y *
sj(i, j, k, 2) - q_z *
sj(i, j, k, 3)
3133 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fmx
3134 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fmy
3135 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fmz
3142 if (j == 1 .and. storewalltensor .and. &
3161 if (j ==
jl .and. storewalltensor .and. &
3174 #ifdef TAPENADE_REVERSE
3190 #ifdef TAPENADE_REVERSE
3192 do ii = 0,
il *
ny *
nz - 1
3194 j = mod(ii /
il,
ny) + 2
3195 k = ii / (
il *
ny) + 2
3214 mul = por * (
rlv(i, j, k) +
rlv(i + 1, j, k))
3222 heatcoef = mul * factlamheat + mue * factturbheat
3227 u_x =
fourth * (
ux(i, j - 1, k - 1) +
ux(i, j, k - 1) &
3228 +
ux(i, j - 1, k) +
ux(i, j, k))
3229 u_y =
fourth * (
uy(i, j - 1, k - 1) +
uy(i, j, k - 1) &
3230 +
uy(i, j - 1, k) +
uy(i, j, k))
3231 u_z =
fourth * (
uz(i, j - 1, k - 1) +
uz(i, j, k - 1) &
3232 +
uz(i, j - 1, k) +
uz(i, j, k))
3234 v_x =
fourth * (
vx(i, j - 1, k - 1) +
vx(i, j, k - 1) &
3235 +
vx(i, j - 1, k) +
vx(i, j, k))
3236 v_y =
fourth * (
vy(i, j - 1, k - 1) +
vy(i, j, k - 1) &
3237 +
vy(i, j - 1, k) +
vy(i, j, k))
3238 v_z =
fourth * (
vz(i, j - 1, k - 1) +
vz(i, j, k - 1) &
3239 +
vz(i, j - 1, k) +
vz(i, j, k))
3241 w_x =
fourth * (
wx(i, j - 1, k - 1) +
wx(i, j, k - 1) &
3242 +
wx(i, j - 1, k) +
wx(i, j, k))
3243 w_y =
fourth * (
wy(i, j - 1, k - 1) +
wy(i, j, k - 1) &
3244 +
wy(i, j - 1, k) +
wy(i, j, k))
3245 w_z =
fourth * (
wz(i, j - 1, k - 1) +
wz(i, j, k - 1) &
3246 +
wz(i, j - 1, k) +
wz(i, j, k))
3248 q_x =
fourth * (
qx(i, j - 1, k - 1) +
qx(i, j, k - 1) &
3249 +
qx(i, j - 1, k) +
qx(i, j, k))
3250 q_y =
fourth * (
qy(i, j - 1, k - 1) +
qy(i, j, k - 1) &
3251 +
qy(i, j - 1, k) +
qy(i, j, k))
3252 q_z =
fourth * (
qz(i, j - 1, k - 1) +
qz(i, j, k - 1) &
3253 +
qz(i, j - 1, k) +
qz(i, j, k))
3260 ssx =
eighth * (
x(i + 1, j - 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
3261 +
x(i + 1, j - 1, k, 1) -
x(i - 1, j - 1, k, 1) &
3262 +
x(i + 1, j, k - 1, 1) -
x(i - 1, j, k - 1, 1) &
3263 +
x(i + 1, j, k, 1) -
x(i - 1, j, k, 1))
3264 ssy =
eighth * (
x(i + 1, j - 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
3265 +
x(i + 1, j - 1, k, 2) -
x(i - 1, j - 1, k, 2) &
3266 +
x(i + 1, j, k - 1, 2) -
x(i - 1, j, k - 1, 2) &
3267 +
x(i + 1, j, k, 2) -
x(i - 1, j, k, 2))
3268 ssz =
eighth * (
x(i + 1, j - 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
3269 +
x(i + 1, j - 1, k, 3) -
x(i - 1, j - 1, k, 3) &
3270 +
x(i + 1, j, k - 1, 3) -
x(i - 1, j, k - 1, 3) &
3271 +
x(i + 1, j, k, 3) -
x(i - 1, j, k, 3))
3276 ss =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
3283 corr = u_x * ssx + u_y * ssy + u_z * ssz &
3284 - (
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)) * ss
3285 u_x = u_x - corr * ssx
3286 u_y = u_y - corr * ssy
3287 u_z = u_z - corr * ssz
3289 corr = v_x * ssx + v_y * ssy + v_z * ssz &
3290 - (
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)) * ss
3291 v_x = v_x - corr * ssx
3292 v_y = v_y - corr * ssy
3293 v_z = v_z - corr * ssz
3295 corr = w_x * ssx + w_y * ssy + w_z * ssz &
3296 - (
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)) * ss
3297 w_x = w_x - corr * ssx
3298 w_y = w_y - corr * ssy
3299 w_z = w_z - corr * ssz
3301 corr = q_x * ssx + q_y * ssy + q_z * ssz &
3302 + (
aa(i + 1, j, k) -
aa(i, j, k)) * ss
3303 q_x = q_x - corr * ssx
3304 q_y = q_y - corr * ssy
3305 q_z = q_z - corr * ssz
3314 fracdiv = twothird * (u_x + v_y + w_z)
3316 tauxxs =
two * u_x - fracdiv
3317 tauyys =
two * v_y - fracdiv
3318 tauzzs =
two * w_z - fracdiv
3324 q_x = heatcoef * q_x
3325 q_y = heatcoef * q_y
3326 q_z = heatcoef * q_z
3346 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
3347 v_x * v_x + v_y * v_y + v_z * v_z + &
3348 w_x * w_x + w_y * w_y + w_z * w_z)
3352 den = max(den, xminn)
3357 fact = mue * ccr1 / den
3369 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
3370 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
3371 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
3373 exy = fact * (wxy * tauyys + wxz * tauyzs + &
3374 wyx * tauxxs + wyz * tauxzs)
3375 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
3376 wzx * tauxxs + wzy * tauxys)
3377 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
3378 wzx * tauxys + wzy * tauyys)
3381 tauxx = mut * tauxxs - exx
3382 tauyy = mut * tauyys - eyy
3383 tauzz = mut * tauzzs - ezz
3384 tauxy = mut * tauxys - exy
3385 tauxz = mut * tauxzs - exz
3386 tauyz = mut * tauyzs - eyz
3391 tauxx = mut * tauxxs
3392 tauyy = mut * tauyys
3393 tauzz = mut * tauzzs
3394 tauxy = mut * tauxys
3395 tauxz = mut * tauxzs
3396 tauyz = mut * tauyzs
3409 fmx = tauxx *
si(i, j, k, 1) + tauxy *
si(i, j, k, 2) &
3410 + tauxz *
si(i, j, k, 3)
3411 fmy = tauxy *
si(i, j, k, 1) + tauyy *
si(i, j, k, 2) &
3412 + tauyz *
si(i, j, k, 3)
3413 fmz = tauxz *
si(i, j, k, 1) + tauyz *
si(i, j, k, 2) &
3414 + tauzz *
si(i, j, k, 3)
3415 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
si(i, j, k, 1) &
3416 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
si(i, j, k, 2) &
3417 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
si(i, j, k, 3) &
3418 - q_x *
si(i, j, k, 1) - q_y *
si(i, j, k, 2) - q_z *
si(i, j, k, 3)
3427 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fmx
3428 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fmy
3429 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fmz
3436 if (i == 1 .and. storewalltensor .and. &
3455 if (i ==
il .and. storewalltensor .and. &
3471 #ifdef TAPENADE_REVERSE
3482 #ifndef USE_TAPENADE
3497 real(kind=realtype),
parameter :: twothird =
two *
third
3501 integer(kind=intType) :: i, j, k
3502 integer(kind=intType) :: ii, jj, kk
3504 real(kind=realtype) :: rfilv, por, mul, mue, mut, heatcoef
3505 real(kind=realtype) :: gm1, factlamheat, factturbheat
3506 real(kind=realtype) :: u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z
3507 real(kind=realtype) :: q_x, q_y, q_z, ubar, vbar, wbar
3508 real(kind=realtype) :: corr, ssx, ssy, ssz, ss, fracdiv
3509 real(kind=realtype) :: tauxx, tauyy, tauzz
3510 real(kind=realtype) :: tauxy, tauxz, tauyz
3511 real(kind=realtype) :: fmx, fmy, fmz, frhoe
3512 real(kind=realtype) :: dd
3513 logical :: correctForK
3525 ssx =
eighth * (
x(i + 1, j - 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
3526 +
x(i + 1, j - 1, k, 1) -
x(i - 1, j - 1, k, 1) &
3527 +
x(i + 1, j, k - 1, 1) -
x(i - 1, j, k - 1, 1) &
3528 +
x(i + 1, j, k, 1) -
x(i - 1, j, k, 1))
3529 ssy =
eighth * (
x(i + 1, j - 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
3530 +
x(i + 1, j - 1, k, 2) -
x(i - 1, j - 1, k, 2) &
3531 +
x(i + 1, j, k - 1, 2) -
x(i - 1, j, k - 1, 2) &
3532 +
x(i + 1, j, k, 2) -
x(i - 1, j, k, 2))
3533 ssz =
eighth * (
x(i + 1, j - 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
3534 +
x(i + 1, j - 1, k, 3) -
x(i - 1, j - 1, k, 3) &
3535 +
x(i + 1, j, k - 1, 3) -
x(i - 1, j, k - 1, 3) &
3536 +
x(i + 1, j, k, 3) -
x(i - 1, j, k, 3))
3539 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
3545 dd =
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)
3550 dd =
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)
3555 dd =
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)
3560 dd =
aa(i + 1, j, k) -
aa(i, j, k)
3573 mul = por * (
rlv(i, j, k) +
rlv(i + 1, j, k))
3581 heatcoef = mul * factlamheat + mue * factturbheat
3585 fracdiv = twothird * (u_x + v_y + w_z)
3587 tauxx = mut * (
two * u_x - fracdiv)
3588 tauyy = mut * (
two * v_y - fracdiv)
3589 tauzz = mut * (
two * w_z - fracdiv)
3591 tauxy = mut * (u_y + v_x)
3592 tauxz = mut * (u_z + w_x)
3593 tauyz = mut * (v_z + w_y)
3595 q_x = heatcoef * q_x
3596 q_y = heatcoef * q_y
3597 q_z = heatcoef * q_z
3608 fmx = tauxx *
si(i, j, k, 1) + tauxy *
si(i, j, k, 2) + tauxz *
si(i, j, k, 3)
3609 fmy = tauxy *
si(i, j, k, 1) + tauyy *
si(i, j, k, 2) + tauyz *
si(i, j, k, 3)
3610 fmz = tauxz *
si(i, j, k, 1) + tauyz *
si(i, j, k, 2) + tauzz *
si(i, j, k, 3)
3611 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
si(i, j, k, 1) &
3612 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
si(i, j, k, 2) &
3613 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
si(i, j, k, 3) &
3614 - q_x *
si(i, j, k, 1) - q_y *
si(i, j, k, 2) - q_z *
si(i, j, k, 3)
3623 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fmx
3624 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fmy
3625 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fmz
3639 ssx =
eighth * (
x(i - 1, j + 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
3640 +
x(i - 1, j + 1, k, 1) -
x(i - 1, j - 1, k, 1) &
3641 +
x(i, j + 1, k - 1, 1) -
x(i, j - 1, k - 1, 1) &
3642 +
x(i, j + 1, k, 1) -
x(i, j - 1, k, 1))
3643 ssy =
eighth * (
x(i - 1, j + 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
3644 +
x(i - 1, j + 1, k, 2) -
x(i - 1, j - 1, k, 2) &
3645 +
x(i, j + 1, k - 1, 2) -
x(i, j - 1, k - 1, 2) &
3646 +
x(i, j + 1, k, 2) -
x(i, j - 1, k, 2))
3647 ssz =
eighth * (
x(i - 1, j + 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
3648 +
x(i - 1, j + 1, k, 3) -
x(i - 1, j - 1, k, 3) &
3649 +
x(i, j + 1, k - 1, 3) -
x(i, j - 1, k - 1, 3) &
3650 +
x(i, j + 1, k, 3) -
x(i, j - 1, k, 3))
3653 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
3659 dd =
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)
3664 dd =
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)
3669 dd =
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)
3674 dd =
aa(i, j + 1, k) -
aa(i, j, k)
3687 mul = por * (
rlv(i, j, k) +
rlv(i, j + 1, k))
3695 heatcoef = mul * factlamheat + mue * factturbheat
3699 fracdiv = twothird * (u_x + v_y + w_z)
3701 tauxx = mut * (
two * u_x - fracdiv)
3702 tauyy = mut * (
two * v_y - fracdiv)
3703 tauzz = mut * (
two * w_z - fracdiv)
3705 tauxy = mut * (u_y + v_x)
3706 tauxz = mut * (u_z + w_x)
3707 tauyz = mut * (v_z + w_y)
3709 q_x = heatcoef * q_x
3710 q_y = heatcoef * q_y
3711 q_z = heatcoef * q_z
3722 fmx = tauxx *
sj(i, j, k, 1) + tauxy *
sj(i, j, k, 2) + tauxz *
sj(i, j, k, 3)
3723 fmy = tauxy *
sj(i, j, k, 1) + tauyy *
sj(i, j, k, 2) + tauyz *
sj(i, j, k, 3)
3724 fmz = tauxz *
sj(i, j, k, 1) + tauyz *
sj(i, j, k, 2) + tauzz *
sj(i, j, k, 3)
3725 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sj(i, j, k, 1) &
3726 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sj(i, j, k, 2) &
3727 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sj(i, j, k, 3) &
3728 - q_x *
sj(i, j, k, 1) - q_y *
sj(i, j, k, 2) - q_z *
sj(i, j, k, 3)
3737 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fmx
3738 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fmy
3739 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fmz
3753 ssx =
eighth * (
x(i - 1, j - 1, k + 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
3754 +
x(i - 1, j, k + 1, 1) -
x(i - 1, j, k - 1, 1) &
3755 +
x(i, j - 1, k + 1, 1) -
x(i, j - 1, k - 1, 1) &
3756 +
x(i, j, k + 1, 1) -
x(i, j, k - 1, 1))
3757 ssy =
eighth * (
x(i - 1, j - 1, k + 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
3758 +
x(i - 1, j, k + 1, 2) -
x(i - 1, j, k - 1, 2) &
3759 +
x(i, j - 1, k + 1, 2) -
x(i, j - 1, k - 1, 2) &
3760 +
x(i, j, k + 1, 2) -
x(i, j, k - 1, 2))
3761 ssz =
eighth * (
x(i - 1, j - 1, k + 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
3762 +
x(i - 1, j, k + 1, 3) -
x(i - 1, j, k - 1, 3) &
3763 +
x(i, j - 1, k + 1, 3) -
x(i, j - 1, k - 1, 3) &
3764 +
x(i, j, k + 1, 3) -
x(i, j, k - 1, 3))
3766 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
3772 dd =
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)
3777 dd =
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)
3782 dd =
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)
3787 dd =
aa(i, j, k + 1) -
aa(i, j, k)
3800 mul = por * (
rlv(i, j, k) +
rlv(i, j, k + 1))
3808 heatcoef = mul * factlamheat + mue * factturbheat
3812 fracdiv = twothird * (u_x + v_y + w_z)
3814 tauxx = mut * (
two * u_x - fracdiv)
3815 tauyy = mut * (
two * v_y - fracdiv)
3816 tauzz = mut * (
two * w_z - fracdiv)
3818 tauxy = mut * (u_y + v_x)
3819 tauxz = mut * (u_z + w_x)
3820 tauyz = mut * (v_z + w_y)
3822 q_x = heatcoef * q_x
3823 q_y = heatcoef * q_y
3824 q_z = heatcoef * q_z
3835 fmx = tauxx *
sk(i, j, k, 1) + tauxy *
sk(i, j, k, 2) + tauxz *
sk(i, j, k, 3)
3836 fmy = tauxy *
sk(i, j, k, 1) + tauyy *
sk(i, j, k, 2) + tauyz *
sk(i, j, k, 3)
3837 fmz = tauxz *
sk(i, j, k, 1) + tauyz *
sk(i, j, k, 2) + tauzz *
sk(i, j, k, 3)
3838 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sk(i, j, k, 1) &
3839 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sk(i, j, k, 2) &
3840 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sk(i, j, k, 3) &
3841 - q_x *
sk(i, j, k, 1) - q_y *
sk(i, j, k, 2) - q_z *
sk(i, j, k, 3)
3850 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fmx
3851 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fmy
3852 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fmz
3880 real(kind=realtype),
parameter :: dssmax = 0.25_realtype
3884 integer(kind=intType) :: i, j, k, ind
3886 real(kind=realtype) :: sslim, rhoi
3887 real(kind=realtype) :: sfil, fis2, fis4
3888 real(kind=realtype) :: ppor, rrad, dis2
3889 real(kind=realtype) :: dss1, dss2, ddw, fs
4047 rrad = ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
4055 dis2 = fis2 * rrad * min(dssmax, max(dss1, dss2)) +
sigma * fis4 * rrad
4061 ddw =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
4069 ddw =
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)
4072 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
4077 ddw =
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)
4080 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
4085 ddw =
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)
4088 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
4133 rrad = ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
4136 dis2 = fis2 * rrad * min(dssmax, max(dss1, dss2)) +
sigma * fis4 * rrad
4142 ddw =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
4150 ddw =
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)
4153 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
4158 ddw =
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)
4161 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
4166 ddw =
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)
4169 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
4214 rrad = ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
4217 dis2 = fis2 * rrad * min(dssmax, max(dss1, dss2)) +
sigma * fis4 * rrad
4223 ddw =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
4231 ddw =
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)
4234 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
4239 ddw =
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)
4242 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
4247 ddw =
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)
4250 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
4278 w(i, j, k,
ivx) =
w(i, j, k,
ivx) * rhoi
4279 w(i, j, k,
ivy) =
w(i, j, k,
ivy) * rhoi
4280 w(i, j, k,
ivz) =
w(i, j, k,
ivz) * rhoi
4289 w(0, j, k,
ivx) =
w(0, j, k,
ivx) * rhoi
4290 w(0, j, k,
ivy) =
w(0, j, k,
ivy) * rhoi
4291 w(0, j, k,
ivz) =
w(0, j, k,
ivz) * rhoi
4295 w(1, j, k,
ivx) =
w(1, j, k,
ivx) * rhoi
4296 w(1, j, k,
ivy) =
w(1, j, k,
ivy) * rhoi
4297 w(1, j, k,
ivz) =
w(1, j, k,
ivz) * rhoi
4317 w(i, 0, k,
ivx) =
w(i, 0, k,
ivx) * rhoi
4318 w(i, 0, k,
ivy) =
w(i, 0, k,
ivy) * rhoi
4319 w(i, 0, k,
ivz) =
w(i, 0, k,
ivz) * rhoi
4323 w(i, 1, k,
ivx) =
w(i, 1, k,
ivx) * rhoi
4324 w(i, 1, k,
ivy) =
w(i, 1, k,
ivy) * rhoi
4325 w(i, 1, k,
ivz) =
w(i, 1, k,
ivz) * rhoi
4365 real(kind=realtype),
parameter :: dpmax = 0.25_realtype
4366 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
4367 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
4368 real(kind=realtype),
parameter :: omega = 0.5_realtype
4369 real(kind=realtype),
parameter :: oneminomega =
one - omega
4373 integer(kind=intType) :: i, j, k, ind
4375 real(kind=realtype) :: plim, sface
4376 real(kind=realtype) :: sfil, fis2, fis4
4377 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53
4378 real(kind=realtype) :: ppor, rrad, dis2
4379 real(kind=realtype) :: dp1, dp2, ddw, tmp, fs
4380 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
4381 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
4382 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
4383 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
4384 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
4385 logical :: correctForK
4462 dis2 = fis2 * ppor * min(dpmax, max(dp1, dp2)) +
sigma * fis4 * ppor
4467 ddw =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
4470 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
4474 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
4478 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
4490 if (correctfork)
then
4491 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
4505 gm1 = gammaavg -
one
4514 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
4517 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
4518 area = sqrt(sx**2 + sy**2 + sz**2)
4519 tmp =
one / max(1.e-25_realtype, area)
4524 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4525 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4527 unavg = uavg * sx + vavg * sy + wavg * sz
4529 ova2avg =
one / a2avg
4540 lam1 = abs(unavg - sface + aavg)
4541 lam2 = abs(unavg - sface - aavg)
4542 lam3 = abs(unavg - sface)
4546 lam1 = max(lam1, epsacoustic * rrad)
4547 lam2 = max(lam2, epsacoustic * rrad)
4548 lam3 = max(lam3, epsshear * rrad)
4560 abv1 =
half * (lam1 + lam2)
4561 abv2 =
half * (lam1 - lam2)
4564 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4565 - wavg * drw + dre) - gm53 * drk
4566 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4568 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4569 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4574 fs = lam3 * dr + abv6
4580 fs = lam3 * dru + uavg * abv6 + sx * abv7
4581 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
4586 fs = lam3 * drv + vavg * abv6 + sy * abv7
4587 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
4592 fs = lam3 * drw + wavg * abv6 + sz * abv7
4593 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
4598 fs = lam3 * dre + havg * abv6 + unavg * abv7
4641 dis2 = fis2 * ppor * min(dpmax, max(dp1, dp2)) +
sigma * fis4 * ppor
4646 ddw =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
4649 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
4653 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
4657 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
4669 if (correctfork)
then
4670 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
4684 gm1 = gammaavg -
one
4693 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
4696 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
4697 area = sqrt(sx**2 + sy**2 + sz**2)
4698 tmp =
one / max(1.e-25_realtype, area)
4703 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4704 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4706 unavg = uavg * sx + vavg * sy + wavg * sz
4708 ova2avg =
one / a2avg
4719 lam1 = abs(unavg - sface + aavg)
4720 lam2 = abs(unavg - sface - aavg)
4721 lam3 = abs(unavg - sface)
4725 lam1 = max(lam1, epsacoustic * rrad)
4726 lam2 = max(lam2, epsacoustic * rrad)
4727 lam3 = max(lam3, epsshear * rrad)
4739 abv1 =
half * (lam1 + lam2)
4740 abv2 =
half * (lam1 - lam2)
4743 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4744 - wavg * drw + dre) - gm53 * drk
4745 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4747 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4748 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4753 fs = lam3 * dr + abv6
4759 fs = lam3 * dru + uavg * abv6 + sx * abv7
4760 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
4765 fs = lam3 * drv + vavg * abv6 + sy * abv7
4766 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
4771 fs = lam3 * drw + wavg * abv6 + sz * abv7
4772 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
4777 fs = lam3 * dre + havg * abv6 + unavg * abv7
4820 dis2 = fis2 * ppor * min(dpmax, max(dp1, dp2)) +
sigma * fis4 * ppor
4825 ddw =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
4828 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
4832 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
4836 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
4848 if (correctfork)
then
4849 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
4862 gm1 = gammaavg -
one
4871 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
4874 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
4875 area = sqrt(sx**2 + sy**2 + sz**2)
4876 tmp =
one / max(1.e-25_realtype, area)
4881 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4882 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4884 unavg = uavg * sx + vavg * sy + wavg * sz
4886 ova2avg =
one / a2avg
4897 lam1 = abs(unavg - sface + aavg)
4898 lam2 = abs(unavg - sface - aavg)
4899 lam3 = abs(unavg - sface)
4903 lam1 = max(lam1, epsacoustic * rrad)
4904 lam2 = max(lam2, epsacoustic * rrad)
4905 lam3 = max(lam3, epsshear * rrad)
4917 abv1 =
half * (lam1 + lam2)
4918 abv2 =
half * (lam1 - lam2)
4921 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4922 - wavg * drw + dre) - gm53 * drk
4923 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4925 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4926 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4931 fs = lam3 * dr + abv6
4937 fs = lam3 * dru + uavg * abv6 + sx * abv7
4938 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
4943 fs = lam3 * drv + vavg * abv6 + sy * abv7
4944 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
4949 fs = lam3 * drw + wavg * abv6 + sz * abv7
4950 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
4955 fs = lam3 * dre + havg * abv6 + unavg * abv7
4975 #ifndef USE_TAPENADE
4995 integer(kind=intType) :: i, j, k
4997 real(kind=realtype) :: sfil, fis0, dis0, ppor, fs, rhoi
5053 dis0 = fis0 * ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
5058 fs = dis0 * (
w(i + 1, j, k,
irho) -
w(i, j, k,
irho))
5064 fs = dis0 * (
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx))
5065 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
5070 fs = dis0 * (
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy))
5071 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
5076 fs = dis0 * (
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz))
5077 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
5082 fs = dis0 * (
w(i + 1, j, k,
irhoe) -
w(i, j, k,
irhoe))
5101 dis0 = fis0 * ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
5106 fs = dis0 * (
w(i, j + 1, k,
irho) -
w(i, j, k,
irho))
5112 fs = dis0 * (
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx))
5113 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
5118 fs = dis0 * (
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy))
5119 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
5124 fs = dis0 * (
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz))
5125 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
5130 fs = dis0 * (
w(i, j + 1, k,
irhoe) -
w(i, j, k,
irhoe))
5149 dis0 = fis0 * ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
5154 fs = dis0 * (
w(i, j, k + 1,
irho) -
w(i, j, k,
irho))
5160 fs = dis0 * (
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx))
5161 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
5166 fs = dis0 * (
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy))
5167 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
5172 fs = dis0 * (
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz))
5173 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
5178 fs = dis0 * (
w(i, j, k + 1,
irhoe) -
w(i, j, k,
irhoe))
5194 w(i, j, k,
ivx) =
w(i, j, k,
ivx) * rhoi
5195 w(i, j, k,
ivy) =
w(i, j, k,
ivy) * rhoi
5196 w(i, j, k,
ivz) =
w(i, j, k,
ivz) * rhoi
5216 use blockpointers,
only:
il,
jl,
kl,
ie,
je,
ke,
ib,
jb,
kb,
w,
p, &
5227 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
5228 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
5232 integer(kind=intType) :: i, j, k
5234 real(kind=realtype) :: sfil, fis0, dis0, ppor, rrad, sface
5235 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53, tmp, fs
5236 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
5237 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
5238 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
5239 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
5240 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
5241 logical :: correctForK
5295 dr = dis0 * (
w(i + 1, j, k,
irho) -
w(i, j, k,
irho))
5296 dru = dis0 * (
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
5298 drv = dis0 * (
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
5300 drw = dis0 * (
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
5302 dre = dis0 * (
w(i + 1, j, k,
irhoe) -
w(i, j, k,
irhoe))
5309 if (correctfork)
then
5310 drk = dis0 * (
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
5322 gm1 = gammaavg -
one
5331 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
5334 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
5335 area = sqrt(sx**2 + sy**2 + sz**2)
5336 tmp =
one / max(1.e-25_realtype, area)
5341 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
5342 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
5344 unavg = uavg * sx + vavg * sy + wavg * sz
5346 ova2avg =
one / a2avg
5357 lam1 = abs(unavg - sface + aavg)
5358 lam2 = abs(unavg - sface - aavg)
5359 lam3 = abs(unavg - sface)
5363 lam1 = max(lam1, epsacoustic * rrad)
5364 lam2 = max(lam2, epsacoustic * rrad)
5365 lam3 = max(lam3, epsshear * rrad)
5377 abv1 =
half * (lam1 + lam2)
5378 abv2 =
half * (lam1 - lam2)
5381 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
5382 - wavg * drw + dre) - gm53 * drk
5383 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
5385 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
5386 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
5391 fs = lam3 * dr + abv6
5397 fs = lam3 * dru + uavg * abv6 + sx * abv7
5398 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
5403 fs = lam3 * drv + vavg * abv6 + sy * abv7
5404 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
5409 fs = lam3 * drw + wavg * abv6 + sz * abv7
5410 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
5415 fs = lam3 * dre + havg * abv6 + unavg * abv7
5439 dr = dis0 * (
w(i, j + 1, k,
irho) -
w(i, j, k,
irho))
5440 dru = dis0 * (
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
5442 drv = dis0 * (
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
5444 drw = dis0 * (
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
5446 dre = dis0 * (
w(i, j + 1, k,
irhoe) -
w(i, j, k,
irhoe))
5453 if (correctfork)
then
5454 drk = dis0 * (
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
5466 gm1 = gammaavg -
one
5475 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
5478 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
5479 area = sqrt(sx**2 + sy**2 + sz**2)
5480 tmp =
one / max(1.e-25_realtype, area)
5485 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
5486 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
5488 unavg = uavg * sx + vavg * sy + wavg * sz
5490 ova2avg =
one / a2avg
5501 lam1 = abs(unavg - sface + aavg)
5502 lam2 = abs(unavg - sface - aavg)
5503 lam3 = abs(unavg - sface)
5507 lam1 = max(lam1, epsacoustic * rrad)
5508 lam2 = max(lam2, epsacoustic * rrad)
5509 lam3 = max(lam3, epsshear * rrad)
5521 abv1 =
half * (lam1 + lam2)
5522 abv2 =
half * (lam1 - lam2)
5525 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
5526 - wavg * drw + dre) - gm53 * drk
5527 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
5529 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
5530 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
5535 fs = lam3 * dr + abv6
5541 fs = lam3 * dru + uavg * abv6 + sx * abv7
5542 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
5547 fs = lam3 * drv + vavg * abv6 + sy * abv7
5548 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
5553 fs = lam3 * drw + wavg * abv6 + sz * abv7
5554 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
5559 fs = lam3 * dre + havg * abv6 + unavg * abv7
5583 dr = dis0 * (
w(i, j, k + 1,
irho) -
w(i, j, k,
irho))
5584 dru = dis0 * (
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
5586 drv = dis0 * (
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
5588 drw = dis0 * (
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
5590 dre = dis0 * (
w(i, j, k + 1,
irhoe) -
w(i, j, k,
irhoe))
5597 if (correctfork)
then
5598 drk = dis0 * (
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
5610 gm1 = gammaavg -
one
5619 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
5622 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
5623 area = sqrt(sx**2 + sy**2 + sz**2)
5624 tmp =
one / max(1.e-25_realtype, area)
5629 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
5630 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
5632 unavg = uavg * sx + vavg * sy + wavg * sz
5634 ova2avg =
one / a2avg
5645 lam1 = abs(unavg - sface + aavg)
5646 lam2 = abs(unavg - sface - aavg)
5647 lam3 = abs(unavg - sface)
5651 lam1 = max(lam1, epsacoustic * rrad)
5652 lam2 = max(lam2, epsacoustic * rrad)
5653 lam3 = max(lam3, epsshear * rrad)
5665 abv1 =
half * (lam1 + lam2)
5666 abv2 =
half * (lam1 - lam2)
5669 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
5670 - wavg * drw + dre) - gm53 * drk
5671 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
5673 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
5674 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
5679 fs = lam3 * dr + abv6
5685 fs = lam3 * dru + uavg * abv6 + sx * abv7
5686 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
5691 fs = lam3 * drv + vavg * abv6 + sy * abv7
5692 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
5697 fs = lam3 * drw + wavg * abv6 + sz * abv7
5698 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
5703 fs = lam3 * dre + havg * abv6 + unavg * abv7
subroutine riemannflux(left, right, flux)
subroutine leftrightstate(du1, du2, du3, rotmatrix, left, right)
real(kind=realtype), dimension(:, :, :), pointer sfacek
integer(kind=inttype), dimension(:, :), pointer viscjminpointer
real(kind=realtype), dimension(:, :, :), pointer gamma
real(kind=realtype), dimension(:, :, :), pointer qz
logical addgridvelocities
real(kind=realtype), dimension(:, :, :), pointer radk
integer(kind=portype), dimension(:, :, :), pointer pork
integer(kind=inttype), dimension(:, :, :), pointer indfamilyj
integer(kind=inttype), dimension(:, :), pointer viscjmaxpointer
real(kind=realtype), dimension(:, :, :), pointer qy
real(kind=realtype), dimension(:, :, :), pointer aa
real(kind=realtype), dimension(:, :, :), pointer uz
integer(kind=inttype), dimension(:, :, :), pointer factfamilyj
integer(kind=inttype) spectralsol
real(kind=realtype), dimension(:, :, :), pointer p
real(kind=realtype), dimension(:, :, :), pointer radj
real(kind=realtype), dimension(:, :, :, :), pointer w
real(kind=realtype), dimension(:, :, :), pointer uy
integer(kind=inttype), dimension(:, :, :), pointer indfamilyk
real(kind=realtype), dimension(:, :, :), pointer sfacei
type(viscsubfacetype), dimension(:), pointer viscsubface
integer(kind=portype), dimension(:, :, :), pointer porj
integer(kind=portype), dimension(:, :, :), pointer pori
real(kind=realtype), dimension(:, :, :), pointer wx
integer(kind=inttype) nbkglobal
integer(kind=inttype), dimension(:, :), pointer viscimaxpointer
real(kind=realtype), dimension(:, :, :), pointer rlv
real(kind=realtype), dimension(:, :, :, :), pointer si
real(kind=realtype), dimension(:, :, :, :), pointer sj
integer(kind=inttype), dimension(:, :), pointer visckminpointer
integer(kind=inttype), dimension(:, :, :), pointer factfamilyi
real(kind=realtype), dimension(:, :, :), pointer qx
integer(kind=inttype), dimension(:, :, :), pointer factfamilyk
real(kind=realtype), dimension(:, :, :), pointer vz
real(kind=realtype), dimension(:, :, :, :, :), pointer rotmatrixj
real(kind=realtype), dimension(:, :, :), pointer rev
real(kind=realtype), dimension(:, :, :, :), pointer dw
real(kind=realtype), dimension(:, :, :, :), pointer sk
real(kind=realtype), dimension(:, :, :), pointer ux
real(kind=realtype), dimension(:, :, :), pointer shocksensor
real(kind=realtype), dimension(:, :, :), pointer wy
real(kind=realtype), dimension(:, :, :), pointer wz
real(kind=realtype), dimension(:, :, :), pointer vol
real(kind=realtype), dimension(:, :, :), pointer sfacej
real(kind=realtype), dimension(:, :, :), pointer vy
real(kind=realtype), dimension(:, :, :, :), pointer fw
real(kind=realtype), dimension(:, :, :), pointer vx
real(kind=realtype), dimension(:, :, :), pointer radi
integer(kind=inttype), dimension(:, :, :), pointer indfamilyi
real(kind=realtype), dimension(:, :, :, :), pointer x
integer(kind=inttype), dimension(:, :), pointer visckmaxpointer
integer(kind=inttype), dimension(:, :), pointer visciminpointer
real(kind=realtype), dimension(:, :, :, :, :), pointer rotmatrixi
real(kind=realtype), dimension(:, :, :, :, :), pointer rotmatrixk
type(cgnsblockinfotype), dimension(:), allocatable cgnsdoms
real(kind=realtype), dimension(:, :), allocatable massflowfamilydiss
real(kind=realtype), dimension(:, :), allocatable massflowfamilyinv
integer(kind=inttype), parameter roe
integer(kind=inttype), parameter vanleer
integer(kind=inttype), parameter firstorder
real(kind=realtype), parameter zero
real(kind=realtype), parameter three
real(kind=realtype), parameter third
real(kind=realtype), parameter thresholdreal
integer(kind=inttype), parameter ausmdv
real(kind=realtype), parameter eighth
integer(kind=inttype), parameter turkel
integer(kind=inttype), parameter choimerkle
integer(kind=inttype), parameter vanalbeda
integer(kind=portype), parameter boundflux
integer(kind=portype), parameter noflux
integer(kind=inttype), parameter eulerequations
integer(kind=portype), parameter normalflux
real(kind=realtype), parameter five
integer(kind=inttype), parameter noprecond
real(kind=realtype), parameter one
real(kind=realtype), parameter half
integer(kind=inttype), parameter steady
real(kind=realtype), parameter two
integer(kind=inttype), parameter minmod
real(kind=realtype), parameter fourth
integer(kind=inttype), parameter nolimiter
integer(kind=inttype), parameter nsequations
integer(kind=inttype), parameter ransequations
subroutine etot(rho, u, v, w, p, k, etotal, correctForK)
real(kind=realtype) gammainf
real(kind=realtype) pinfcorr
integer(kind=inttype) nwf
real(kind=realtype) rhoinf
real(kind=realtype) timeref
subroutine invisciddissfluxmatrixcoarse
subroutine invisciddissfluxmatrixapprox
subroutine invisciddissfluxscalarcoarse
subroutine invisciddissfluxscalar
subroutine invisciddissfluxmatrix
subroutine inviscidcentralflux
subroutine viscousfluxapprox
subroutine inviscidupwindflux(fineGrid)
subroutine invisciddissfluxscalarapprox
real(kind=realtype) totalr0
integer(kind=inttype) currentlevel
real(kind=realtype) totalr
integer(kind=inttype) groundlevel
integer(kind=inttype) rkstage
logical function getcorrectfork()
real(kind=realtype) function mydim(x, y)
subroutine terminate(routineName, errorMessage)