9 integer(kind=intType),
parameter ::
bs = 8
15 integer(kind=intType) ::
nx,
ny,
nz,
il,
jl,
kl,
ie,
je,
ke,
ib,
jb,
kb
24 real(kind=realtype),
dimension(0:bbib, 0:bbjb, 0:bbkb, 1:6) ::
w
25 real(kind=realtype),
dimension(0:bbib, 0:bbjb, 0:bbkb) ::
p,
gamma
26 real(kind=realtype),
dimension(0:bbib, 0:bbjb, 0:bbkb) ::
ss
29 real(kind=realtype),
dimension(0:bbie, 0:bbje, 0:bbke, 3) ::
x
30 real(kind=realtype),
dimension(1:bbie, 1:bbje, 1:bbke) ::
rlv,
rev,
vol,
aa
32 real(kind=realtype),
dimension(1:bbie, 1:bbje, 1:bbke, 3) ::
dss
35 real(kind=realtype),
dimension(2:bbil, 2:bbjl, 2:bbkl) ::
volref,
d2wall
36 integer(kind=intType),
dimension(2:bbil, 2:bbjl, 2:bbkl) ::
iblank
39 integer(kind=porType),
dimension(1:bbil, 2:bbjl, 2:bbkl) ::
pori
40 integer(kind=porType),
dimension(2:bbil, 1:bbjl, 2:bbkl) ::
porj
41 integer(kind=porType),
dimension(2:bbil, 2:bbjl, 1:bbkl) ::
pork
44 real(kind=realtype),
dimension(1:bbie, 1:bbje, 1:bbke, 1:5) ::
fw
45 real(kind=realtype),
dimension(1:bbie, 1:bbje, 1:bbke, 1:6) ::
dw
48 real(kind=realtype),
dimension(0:bbie, 1:bbje, 1:bbke, 3) ::
si
49 real(kind=realtype),
dimension(1:bbie, 0:bbje, 1:bbke, 3) ::
sj
50 real(kind=realtype),
dimension(1:bbie, 1:bbje, 0:bbke, 3) ::
sk
53 real(kind=realtype),
dimension(0:bbie, 1:bbje, 1:bbke) ::
sfacei
54 real(kind=realtype),
dimension(1:bbie, 0:bbje, 1:bbke) ::
sfacej
55 real(kind=realtype),
dimension(1:bbie, 1:bbje, 0:bbke) ::
sfacek
58 real(kind=realtype),
dimension(1:bbil, 1:bbjl, 1:bbkl) ::
ux,
uy,
uz
59 real(kind=realtype),
dimension(1:bbil, 1:bbjl, 1:bbkl) ::
vx,
vy,
vz
60 real(kind=realtype),
dimension(1:bbil, 1:bbjl, 1:bbkl) ::
wx,
wy,
wz
61 real(kind=realtype),
dimension(1:bbil, 1:bbjl, 1:bbkl) ::
qx,
qy,
qz
70 subroutine blocketteres(useDissApprox, useViscApprox, useUpdateIntermed, useFlowRes, useTurbRes, useSpatial, &
71 useStoreWall, famLists, funcValues, forces, bcDataNames, bcDataValues, bcDataFamLists)
103 logical,
intent(in),
optional :: useDissApprox, useViscApprox, useUpdateIntermed, useFlowRes
104 logical,
intent(in),
optional :: useTurbRes, useSpatial, useStoreWall
105 integer(kind=intType),
optional,
dimension(:, :),
intent(in) :: famLists
106 real(kind=realtype),
optional,
dimension(:, :),
intent(out) :: funcvalues
107 character,
optional,
dimension(:, :),
intent(in) :: bcDataNames
108 real(kind=realtype),
optional,
dimension(:),
intent(in) :: bcdatavalues
109 integer(kind=intType),
optional,
dimension(:, :) :: bcDataFamLists
110 real(kind=realtype),
intent(out),
optional,
dimension(:, :, :) :: forces
113 logical :: dissApprox, viscApprox, updateIntermed, flowRes, turbRes, spatial, storeWall
114 integer(kind=intType) :: nn, sps, fSize, lstart, lend, iRegion
115 real(kind=realtype) :: plocal
131 updateintermed = .false.
138 if (
present(usedissapprox))
then
139 dissapprox = usedissapprox
142 if (
present(useviscapprox))
then
143 viscapprox = useviscapprox
146 if (
present(useupdateintermed))
then
147 updateintermed = useupdateintermed
150 if (
present(useflowres))
then
154 if (
present(useturbres))
then
158 if (
present(usespatial))
then
162 if (
present(usestorewall))
then
163 storewall = usestorewall
172 if (
present(bcdatanames))
then
174 call setbcdata(bcdatanames, bcdatavalues, bcdatafamlists, sps, &
175 size(bcdatavalues),
size(bcdatafamlists, 2))
232 if (flowres .and. turbres)
then
236 else if (flowres .and. (.not. turbres))
then
240 else if ((.not. flowres) .and. turbres)
then
246 call whalo2(1_inttype, lstart, lend, .true., .true., .true.)
267 blockloop:
do nn = 1,
ndom
272 call blocketterescore(dissapprox, viscapprox, updateintermed, flowres, turbres, storewall)
274 call blockrescore(dissapprox, viscapprox, updateintermed, flowres, turbres, storewall, nn, sps)
286 if (
present(famlists))
then
290 if (
present(forces))
then
293 fsize =
size(forces, 2)
294 call getforces(forces(:, :, sps), fsize, sps)
299 subroutine blocketterescore(dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall)
306 bnx =>
nx, bny =>
ny, bnz =>
nz, &
307 bil =>
il, bjl =>
jl, bkl =>
kl, &
308 bie =>
ie, bje =>
je, bke =>
ke, &
309 bib =>
ib, bjb =>
jb, bkb =>
kb, &
310 bw =>
w, bp =>
p, bgamma =>
gamma, &
312 bux =>
ux, buy =>
uy, buz =>
uz, &
313 bvx =>
vx, bvy =>
vy, bvz =>
vz, &
314 bwx =>
wx, bwy =>
wy, bwz =>
wz, &
315 bqx =>
qx, bqy =>
qy, bqz =>
qz, &
319 bsi =>
si, bsj =>
sj, bsk =>
sk, &
321 bdtl =>
dtl, baa =>
aa, &
334 logical,
intent(in) :: dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall
337 integer(kind=intType) :: i, j, k, l, lStart, lEnd
340 if (flowres .and. turbres)
then
344 else if (flowres .and. (.not. turbres))
then
348 else if ((.not. flowres) .and. turbres)
then
371 firstblockette:
if (
ii == 2)
then
384 w(i, j, k, 1:
nw) = bw(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1:
nw)
385 p(i, j, k) = bp(i +
ii - 2, j +
jj - 2, k +
kk - 2)
386 gamma(i, j, k) = bgamma(i +
ii - 2, j +
jj - 2, k +
kk - 2)
388 ss(i, j, k) = bshocksensor(i +
ii - 2, j +
jj - 2, k +
kk - 2)
398 rlv(i, j, k) = brlv(i +
ii - 2, j +
jj - 2, k +
kk - 2)
399 rev(i, j, k) = brev(i +
ii - 2, j +
jj - 2, k +
kk - 2)
400 vol(i, j, k) = bvol(i +
ii - 2, j +
jj - 2, k +
kk - 2)
409 x(i, j, k, :) = bx(i +
ii - 2, j +
jj - 2, k +
kk - 2, :)
442 w(i, j, k, 1:
nw) =
w(
bs + i, j, k, 1:
nw)
443 p(i, j, k) =
p(
bs + i, j, k)
445 ss(i, j, k) =
ss(
bs + i, j, k)
464 aa(i, j, k) =
aa(
bs + i, j, k)
465 dss(i, j, k, :) =
dss(
bs + i, j, k, :)
474 x(i, j, k, :) =
x(
bs + i, j, k, :)
482 ux(1, j, k) =
ux(
bs + 1, j, k)
483 uy(1, j, k) =
uy(
bs + 1, j, k)
484 uz(1, j, k) =
uz(
bs + 1, j, k)
486 vx(1, j, k) =
vx(
bs + 1, j, k)
487 vy(1, j, k) =
vy(
bs + 1, j, k)
488 vz(1, j, k) =
vz(
bs + 1, j, k)
490 wx(1, j, k) =
wx(
bs + 1, j, k)
491 wy(1, j, k) =
wy(
bs + 1, j, k)
492 wz(1, j, k) =
wz(
bs + 1, j, k)
494 qx(1, j, k) =
qx(
bs + 1, j, k)
495 qy(1, j, k) =
qy(
bs + 1, j, k)
496 qz(1, j, k) =
qz(
bs + 1, j, k)
499 end if firstblockette
509 w(i, j, k, 1:
nw) = bw(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1:
nw)
510 p(i, j, k) = bp(i +
ii - 2, j +
jj - 2, k +
kk - 2)
511 gamma(i, j, k) = bgamma(i +
ii - 2, j +
jj - 2, k +
kk - 2)
513 ss(i, j, k) = bshocksensor(i +
ii - 2, j +
jj - 2, k +
kk - 2)
523 rlv(i, j, k) = brlv(i +
ii - 2, j +
jj - 2, k +
kk - 2)
524 rev(i, j, k) = brev(i +
ii - 2, j +
jj - 2, k +
kk - 2)
525 vol(i, j, k) = bvol(i +
ii - 2, j +
jj - 2, k +
kk - 2)
534 x(i, j, k, :) = bx(i +
ii - 2, j +
jj - 2, k +
kk - 2, :)
543 iblank(i, j, k) = biblank(i +
ii - 2, j +
jj - 2, k +
kk - 2)
545 d2wall(i, j, k) = bd2wall(i +
ii - 2, j +
jj - 2, k +
kk - 2)
546 volref(i, j, k) = bvolref(i +
ii - 2, j +
jj - 2, k +
kk - 2)
555 pori(i, j, k) = bpori(i +
ii - 2, j +
jj - 2, k +
kk - 2)
563 porj(i, j, k) = bporj(i +
ii - 2, j +
jj - 2, k +
kk - 2)
571 pork(i, j, k) = bpork(i +
ii - 2, j +
jj - 2, k +
kk - 2)
581 sfacei(i, j, k) = bsfacei(i +
ii - 2, j +
jj - 2, k +
kk - 2)
589 sfacej(i, j, k) = bsfacej(i +
ii - 2, j +
jj - 2, k +
kk - 2)
597 sfacek(i, j, k) = bsfacek(i +
ii - 2, j +
jj - 2, k +
kk - 2)
675 bdw(i +
ii - 2, j +
jj - 2, k +
kk - 2, l) =
dw(i, j, k, l)
688 intermed:
if (updateintermed)
then
693 bdtl(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
dtl(i, j, k)
702 bradi(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
radi(i, j, k)
703 bradj(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
radj(i, j, k)
704 bradk(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
radk(i, j, k)
710 visc:
if (
viscous .and. flowres)
then
716 baa(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
aa(i, j, k)
726 bux(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
ux(i, j, k)
727 buy(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
uy(i, j, k)
728 buz(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
uz(i, j, k)
730 bvx(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
vx(i, j, k)
731 bvy(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
vy(i, j, k)
732 bvz(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
vz(i, j, k)
734 bwx(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
wx(i, j, k)
735 bwy(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
wy(i, j, k)
736 bwz(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
wz(i, j, k)
738 bqx(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
qx(i, j, k)
739 bqy(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
qy(i, j, k)
740 bqz(i +
ii - 2, j +
jj - 2, k +
kk - 2) =
qz(i, j, k)
755 subroutine blockrescore(dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall, nn, sps)
778 logical,
intent(in) :: dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall
779 integer(kind=intType),
intent(in) :: nn, sps
782 integer(kind=intType) :: i, j, k, lStart, lEnd
785 if (flowres .and. turbres)
then
789 else if (flowres .and. (.not. turbres))
then
793 else if ((.not. flowres) .and. turbres)
then
819 call inviscidcentralflux_block
823 call invisciddissfluxscalarapprox_block
825 call invisciddissfluxmatrixapprox_block
827 call inviscidupwindflux_block(.true.)
832 call invisciddissfluxscalar_block
834 call invisciddissfluxmatrix_block
836 call inviscidupwindflux_block(.true.)
841 call computespeedofsoundsquared_block
843 call viscousfluxapprox_block
845 call allnodalgradients_block
846 call viscousflux_block
850 call sumdwandfw_block
863 integer(kind=intType) :: i, j, k, l, m, n
864 real(kind=realtype),
dimension(3) :: v1, v2
865 real(kind=realtype) :: fact
881 v1(1) =
x(i, j, n, 1) -
x(i, m, k, 1)
882 v1(2) =
x(i, j, n, 2) -
x(i, m, k, 2)
883 v1(3) =
x(i, j, n, 3) -
x(i, m, k, 3)
885 v2(1) =
x(i, j, k, 1) -
x(i, m, n, 1)
886 v2(2) =
x(i, j, k, 2) -
x(i, m, n, 2)
887 v2(3) =
x(i, j, k, 3) -
x(i, m, n, 3)
893 si(i, j, k, 1) = fact * (v1(2) * v2(3) - v1(3) * v2(2))
894 si(i, j, k, 2) = fact * (v1(3) * v2(1) - v1(1) * v2(3))
895 si(i, j, k, 3) = fact * (v1(1) * v2(2) - v1(2) * v2(1))
911 v1(1) =
x(i, j, n, 1) -
x(l, j, k, 1)
912 v1(2) =
x(i, j, n, 2) -
x(l, j, k, 2)
913 v1(3) =
x(i, j, n, 3) -
x(l, j, k, 3)
915 v2(1) =
x(l, j, n, 1) -
x(i, j, k, 1)
916 v2(2) =
x(l, j, n, 2) -
x(i, j, k, 2)
917 v2(3) =
x(l, j, n, 3) -
x(i, j, k, 3)
923 sj(i, j, k, 1) = fact * (v1(2) * v2(3) - v1(3) * v2(2))
924 sj(i, j, k, 2) = fact * (v1(3) * v2(1) - v1(1) * v2(3))
925 sj(i, j, k, 3) = fact * (v1(1) * v2(2) - v1(2) * v2(1))
941 v1(1) =
x(i, j, k, 1) -
x(l, m, k, 1)
942 v1(2) =
x(i, j, k, 2) -
x(l, m, k, 2)
943 v1(3) =
x(i, j, k, 3) -
x(l, m, k, 3)
945 v2(1) =
x(l, j, k, 1) -
x(i, m, k, 1)
946 v2(2) =
x(l, j, k, 2) -
x(i, m, k, 2)
947 v2(3) =
x(l, j, k, 3) -
x(i, m, k, 3)
953 sk(i, j, k, 1) = fact * (v1(2) * v2(3) - v1(3) * v2(2))
954 sk(i, j, k, 2) = fact * (v1(3) * v2(1) - v1(1) * v2(3))
955 sk(i, j, k, 3) = fact * (v1(1) * v2(2) - v1(2) * v2(1))
970 integer(kind=intType) :: varStart, varEnd
972 dw(:, :, :, varstart:varend) =
zero
993 real(kind=realtype) :: fv1, fv2, ft2
994 real(kind=realtype) ::
sst, nu, dist2inv, chi, chi2, chi3
995 real(kind=realtype) :: rr, gg, gg6, termfw, fwsa, term1, term2
996 real(kind=realtype) :: dfv1, dfv2, dft2, drr, dgg, dfw, sqrtprod
997 real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
998 real(kind=realtype) :: div2, fact, sxx, syy, szz, sxy, sxz, syz
999 real(kind=realtype) :: vortx, vorty, vortz
1000 real(kind=realtype) :: omegax, omegay, omegaz
1001 real(kind=realtype) :: strainmag2, prod
1002 real(kind=realtype),
parameter :: xminn = 1.e-10_realtype
1003 real(kind=realtype),
parameter :: f23 =
two *
third
1004 integer(kind=intType) :: i, j, k
1005 real(kind=realtype) :: term1fact
1032 uux =
w(i + 1, j, k,
ivx) *
si(i, j, k, 1) -
w(i - 1, j, k,
ivx) *
si(i - 1, j, k, 1) &
1033 +
w(i, j + 1, k,
ivx) *
sj(i, j, k, 1) -
w(i, j - 1, k,
ivx) *
sj(i, j - 1, k, 1) &
1034 +
w(i, j, k + 1,
ivx) *
sk(i, j, k, 1) -
w(i, j, k - 1,
ivx) *
sk(i, j, k - 1, 1)
1035 uuy =
w(i + 1, j, k,
ivx) *
si(i, j, k, 2) -
w(i - 1, j, k,
ivx) *
si(i - 1, j, k, 2) &
1036 +
w(i, j + 1, k,
ivx) *
sj(i, j, k, 2) -
w(i, j - 1, k,
ivx) *
sj(i, j - 1, k, 2) &
1037 +
w(i, j, k + 1,
ivx) *
sk(i, j, k, 2) -
w(i, j, k - 1,
ivx) *
sk(i, j, k - 1, 2)
1038 uuz =
w(i + 1, j, k,
ivx) *
si(i, j, k, 3) -
w(i - 1, j, k,
ivx) *
si(i - 1, j, k, 3) &
1039 +
w(i, j + 1, k,
ivx) *
sj(i, j, k, 3) -
w(i, j - 1, k,
ivx) *
sj(i, j - 1, k, 3) &
1040 +
w(i, j, k + 1,
ivx) *
sk(i, j, k, 3) -
w(i, j, k - 1,
ivx) *
sk(i, j, k - 1, 3)
1044 vvx =
w(i + 1, j, k,
ivy) *
si(i, j, k, 1) -
w(i - 1, j, k,
ivy) *
si(i - 1, j, k, 1) &
1045 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 1) -
w(i, j - 1, k,
ivy) *
sj(i, j - 1, k, 1) &
1046 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 1) -
w(i, j, k - 1,
ivy) *
sk(i, j, k - 1, 1)
1047 vvy =
w(i + 1, j, k,
ivy) *
si(i, j, k, 2) -
w(i - 1, j, k,
ivy) *
si(i - 1, j, k, 2) &
1048 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 2) -
w(i, j - 1, k,
ivy) *
sj(i, j - 1, k, 2) &
1049 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 2) -
w(i, j, k - 1,
ivy) *
sk(i, j, k - 1, 2)
1050 vvz =
w(i + 1, j, k,
ivy) *
si(i, j, k, 3) -
w(i - 1, j, k,
ivy) *
si(i - 1, j, k, 3) &
1051 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 3) -
w(i, j - 1, k,
ivy) *
sj(i, j - 1, k, 3) &
1052 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 3) -
w(i, j, k - 1,
ivy) *
sk(i, j, k - 1, 3)
1056 wwx =
w(i + 1, j, k,
ivz) *
si(i, j, k, 1) -
w(i - 1, j, k,
ivz) *
si(i - 1, j, k, 1) &
1057 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 1) -
w(i, j - 1, k,
ivz) *
sj(i, j - 1, k, 1) &
1058 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 1) -
w(i, j, k - 1,
ivz) *
sk(i, j, k - 1, 1)
1059 wwy =
w(i + 1, j, k,
ivz) *
si(i, j, k, 2) -
w(i - 1, j, k,
ivz) *
si(i - 1, j, k, 2) &
1060 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 2) -
w(i, j - 1, k,
ivz) *
sj(i, j - 1, k, 2) &
1061 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 2) -
w(i, j, k - 1,
ivz) *
sk(i, j, k - 1, 2)
1062 wwz =
w(i + 1, j, k,
ivz) *
si(i, j, k, 3) -
w(i - 1, j, k,
ivz) *
si(i - 1, j, k, 3) &
1063 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 3) -
w(i, j - 1, k,
ivz) *
sj(i, j - 1, k, 3) &
1064 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 3) -
w(i, j, k - 1,
ivz) *
sk(i, j, k - 1, 3)
1074 sxx =
two * fact * uux
1075 syy =
two * fact * vvy
1076 szz =
two * fact * wwz
1078 sxy = fact * (uuy + vvx)
1079 sxz = fact * (uuz + wwx)
1080 syz = fact * (vvz + wwy)
1084 div2 = f23 * (sxx + syy + szz)**2
1088 strainmag2 =
two * (sxy**2 + sxz**2 + syz**2) &
1089 + sxx**2 + syy**2 + szz**2
1096 vortx =
two * fact * (wwy - vvz) -
two * omegax
1097 vorty =
two * fact * (uuz - wwx) -
two * omegay
1098 vortz =
two * fact * (vvx - uuy) -
two * omegaz
1101 sqrtprod = sqrt(max(
two * strainmag2 - div2,
eps))
1103 sqrtprod = sqrt(vortx**2 + vorty**2 + vortz**2)
1111 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1113 chi =
w(i, j, k,
itu1) / nu
1116 fv1 = chi3 / (chi3 +
cv13)
1117 fv2 =
one - chi / (
one + chi * fv1)
1125 ft2 = rsact3 * exp(-rsact4 * chi2)
1150 rr = min(rr, 10.0_realtype)
1151 gg = rr + rsacw2 * (rr**6 - rr)
1159 term1 = rsacb1 * (
one - ft2) * sqrtprod * term1fact
1160 term2 = dist2inv * (
kar2inv * rsacb1 * ((
one - ft2) * fv2 + ft2) &
1163 dw(i, j, k,
itu1) =
dw(i, j, k,
itu1) + (term1 + term2 *
w(i, j, k,
itu1)) *
w(i, j, k,
itu1)
1181 real(kind=realtype) :: voli, volmi, volpi, xm, ym, zm, xp, yp, zp
1182 real(kind=realtype) :: xa, ya, za, ttm, ttp, cnud, cam, cap
1183 real(kind=realtype) :: nutm, nutp, num, nup, cdm, cdp
1184 real(kind=realtype) :: c1m, c1p, c10, b1, c1, d1, qs, nu
1185 integer(Kind=intType) :: i, j, k
1204 voli =
one /
vol(i, j, k)
1205 volmi =
two / (
vol(i, j, k) +
vol(i, j, k - 1))
1206 volpi =
two / (
vol(i, j, k) +
vol(i, j, k + 1))
1208 xm =
sk(i, j, k - 1, 1) * volmi
1209 ym =
sk(i, j, k - 1, 2) * volmi
1210 zm =
sk(i, j, k - 1, 3) * volmi
1211 xp =
sk(i, j, k, 1) * volpi
1212 yp =
sk(i, j, k, 2) * volpi
1213 zp =
sk(i, j, k, 3) * volpi
1215 xa =
half * (
sk(i, j, k, 1) +
sk(i, j, k - 1, 1)) * voli
1216 ya =
half * (
sk(i, j, k, 2) +
sk(i, j, k - 1, 2)) * voli
1217 za =
half * (
sk(i, j, k, 3) +
sk(i, j, k - 1, 3)) * voli
1218 ttm = xm * xa + ym * ya + zm * za
1219 ttp = xp * xa + yp * ya + zp * za
1241 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1242 num =
half * (
rlv(i, j, k - 1) /
w(i, j, k - 1,
irho) + nu)
1243 nup =
half * (
rlv(i, j, k + 1) /
w(i, j, k + 1,
irho) + nu)
1244 cdm = (num + (
one + rsacb2) * nutm) * ttm *
cb3inv
1245 cdp = (nup + (
one + rsacb2) * nutp) * ttp *
cb3inv
1247 c1m = max(cdm + cam,
zero)
1248 c1p = max(cdp + cap,
zero)
1255 - c10 *
w(i, j, k,
itu1) + c1p *
w(i, j, k + 1,
itu1)
1269 voli =
one /
vol(i, j, k)
1270 volmi =
two / (
vol(i, j, k) +
vol(i, j - 1, k))
1271 volpi =
two / (
vol(i, j, k) +
vol(i, j + 1, k))
1273 xm =
sj(i, j - 1, k, 1) * volmi
1274 ym =
sj(i, j - 1, k, 2) * volmi
1275 zm =
sj(i, j - 1, k, 3) * volmi
1276 xp =
sj(i, j, k, 1) * volpi
1277 yp =
sj(i, j, k, 2) * volpi
1278 zp =
sj(i, j, k, 3) * volpi
1280 xa =
half * (
sj(i, j, k, 1) +
sj(i, j - 1, k, 1)) * voli
1281 ya =
half * (
sj(i, j, k, 2) +
sj(i, j - 1, k, 2)) * voli
1282 za =
half * (
sj(i, j, k, 3) +
sj(i, j - 1, k, 3)) * voli
1283 ttm = xm * xa + ym * ya + zm * za
1284 ttp = xp * xa + yp * ya + zp * za
1306 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1307 num =
half * (
rlv(i, j - 1, k) /
w(i, j - 1, k,
irho) + nu)
1308 nup =
half * (
rlv(i, j + 1, k) /
w(i, j + 1, k,
irho) + nu)
1309 cdm = (num + (
one + rsacb2) * nutm) * ttm *
cb3inv
1310 cdp = (nup + (
one + rsacb2) * nutp) * ttp *
cb3inv
1312 c1m = max(cdm + cam,
zero)
1313 c1p = max(cdp + cap,
zero)
1320 - c10 *
w(i, j, k,
itu1) + c1p *
w(i, j + 1, k,
itu1)
1335 voli =
one /
vol(i, j, k)
1336 volmi =
two / (
vol(i, j, k) +
vol(i - 1, j, k))
1337 volpi =
two / (
vol(i, j, k) +
vol(i + 1, j, k))
1339 xm =
si(i - 1, j, k, 1) * volmi
1340 ym =
si(i - 1, j, k, 2) * volmi
1341 zm =
si(i - 1, j, k, 3) * volmi
1342 xp =
si(i, j, k, 1) * volpi
1343 yp =
si(i, j, k, 2) * volpi
1344 zp =
si(i, j, k, 3) * volpi
1346 xa =
half * (
si(i, j, k, 1) +
si(i - 1, j, k, 1)) * voli
1347 ya =
half * (
si(i, j, k, 2) +
si(i - 1, j, k, 2)) * voli
1348 za =
half * (
si(i, j, k, 3) +
si(i - 1, j, k, 3)) * voli
1349 ttm = xm * xa + ym * ya + zm * za
1350 ttp = xp * xa + yp * ya + zp * za
1372 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1373 num =
half * (
rlv(i - 1, j, k) /
w(i - 1, j, k,
irho) + nu)
1374 nup =
half * (
rlv(i + 1, j, k) /
w(i + 1, j, k,
irho) + nu)
1375 cdm = (num + (
one + rsacb2) * nutm) * ttm *
cb3inv
1376 cdp = (nup + (
one + rsacb2) * nutp) * ttp *
cb3inv
1378 c1m = max(cdm + cam,
zero)
1379 c1p = max(cdp + cap,
zero)
1386 - c10 *
w(i, j, k,
itu1) + c1p *
w(i + 1, j, k,
itu1)
1403 real(kind=realtype) :: uu, dwt, dwtm1, dwtp1, dwti, dwtj, dwtk, qs
1404 real(kind=realtype) :: voli, xa, ya, za
1405 integer(kind=intType),
parameter :: nAdv = 1
1406 integer(kind=intType) :: offset, i, j, k, ii, jj
1428 xa = (
sk(i, j, k, 1) +
sk(i, j, k - 1, 1)) * voli
1429 ya = (
sk(i, j, k, 2) +
sk(i, j, k - 1, 2)) * voli
1430 za = (
sk(i, j, k, 3) +
sk(i, j, k - 1, 3)) * voli
1432 uu = xa *
w(i, j, k,
ivx) + ya *
w(i, j, k,
ivy) + za *
w(i, j, k,
ivz) - qs
1437 velkdir:
if (uu >
zero)
then
1457 dwtm1 =
w(i, j, k - 1, jj) -
w(i, j, k - 2, jj)
1458 dwt =
w(i, j, k, jj) -
w(i, j, k - 1, jj)
1459 dwtp1 =
w(i, j, k + 1, jj) -
w(i, j, k, jj)
1467 if (dwt * dwtp1 >
zero)
then
1468 if (abs(dwt) < abs(dwtp1))
then
1469 dwtk = dwtk +
half * dwt
1471 dwtk = dwtk +
half * dwtp1
1475 if (dwt * dwtm1 >
zero)
then
1476 if (abs(dwt) < abs(dwtm1))
then
1477 dwtk = dwtk -
half * dwt
1479 dwtk = dwtk -
half * dwtm1
1487 dwtk =
w(i, j, k, jj) -
w(i, j, k - 1, jj)
1495 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtk
1517 dwtm1 =
w(i, j, k, jj) -
w(i, j, k - 1, jj)
1518 dwt =
w(i, j, k + 1, jj) -
w(i, j, k, jj)
1519 dwtp1 =
w(i, j, k + 2, jj) -
w(i, j, k + 1, jj)
1527 if (dwt * dwtp1 >
zero)
then
1528 if (abs(dwt) < abs(dwtp1))
then
1529 dwtk = dwtk -
half * dwt
1531 dwtk = dwtk -
half * dwtp1
1535 if (dwt * dwtm1 >
zero)
then
1536 if (abs(dwt) < abs(dwtm1))
then
1537 dwtk = dwtk +
half * dwt
1539 dwtk = dwtk +
half * dwtm1
1547 dwtk =
w(i, j, k + 1, jj) -
w(i, j, k, jj)
1555 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtk
1582 xa = (
sj(i, j, k, 1) +
sj(i, j - 1, k, 1)) * voli
1583 ya = (
sj(i, j, k, 2) +
sj(i, j - 1, k, 2)) * voli
1584 za = (
sj(i, j, k, 3) +
sj(i, j - 1, k, 3)) * voli
1586 uu = xa *
w(i, j, k,
ivx) + ya *
w(i, j, k,
ivy) + za *
w(i, j, k,
ivz) - qs
1591 veljdir:
if (uu >
zero)
then
1610 dwtm1 =
w(i, j - 1, k, jj) -
w(i, j - 2, k, jj)
1611 dwt =
w(i, j, k, jj) -
w(i, j - 1, k, jj)
1612 dwtp1 =
w(i, j + 1, k, jj) -
w(i, j, k, jj)
1620 if (dwt * dwtp1 >
zero)
then
1621 if (abs(dwt) < abs(dwtp1))
then
1622 dwtj = dwtj +
half * dwt
1624 dwtj = dwtj +
half * dwtp1
1628 if (dwt * dwtm1 >
zero)
then
1629 if (abs(dwt) < abs(dwtm1))
then
1630 dwtj = dwtj -
half * dwt
1632 dwtj = dwtj -
half * dwtm1
1640 dwtj =
w(i, j, k, jj) -
w(i, j - 1, k, jj)
1648 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtj
1670 dwtm1 =
w(i, j, k, jj) -
w(i, j - 1, k, jj)
1671 dwt =
w(i, j + 1, k, jj) -
w(i, j, k, jj)
1672 dwtp1 =
w(i, j + 2, k, jj) -
w(i, j + 1, k, jj)
1680 if (dwt * dwtp1 >
zero)
then
1681 if (abs(dwt) < abs(dwtp1))
then
1682 dwtj = dwtj -
half * dwt
1684 dwtj = dwtj -
half * dwtp1
1688 if (dwt * dwtm1 >
zero)
then
1689 if (abs(dwt) < abs(dwtm1))
then
1690 dwtj = dwtj +
half * dwt
1692 dwtj = dwtj +
half * dwtm1
1700 dwtj =
w(i, j + 1, k, jj) -
w(i, j, k, jj)
1708 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtj
1734 xa = (
si(i, j, k, 1) +
si(i - 1, j, k, 1)) * voli
1735 ya = (
si(i, j, k, 2) +
si(i - 1, j, k, 2)) * voli
1736 za = (
si(i, j, k, 3) +
si(i - 1, j, k, 3)) * voli
1738 uu = xa *
w(i, j, k,
ivx) + ya *
w(i, j, k,
ivy) + za *
w(i, j, k,
ivz) - qs
1743 velidir:
if (uu >
zero)
then
1762 dwtm1 =
w(i - 1, j, k, jj) -
w(i - 2, j, k, jj)
1763 dwt =
w(i, j, k, jj) -
w(i - 1, j, k, jj)
1764 dwtp1 =
w(i + 1, j, k, jj) -
w(i, j, k, jj)
1772 if (dwt * dwtp1 >
zero)
then
1773 if (abs(dwt) < abs(dwtp1))
then
1774 dwti = dwti +
half * dwt
1776 dwti = dwti +
half * dwtp1
1780 if (dwt * dwtm1 >
zero)
then
1781 if (abs(dwt) < abs(dwtm1))
then
1782 dwti = dwti -
half * dwt
1784 dwti = dwti -
half * dwtm1
1792 dwti =
w(i, j, k, jj) -
w(i - 1, j, k, jj)
1800 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwti
1822 dwtm1 =
w(i, j, k, jj) -
w(i - 1, j, k, jj)
1823 dwt =
w(i + 1, j, k, jj) -
w(i, j, k, jj)
1824 dwtp1 =
w(i + 2, j, k, jj) -
w(i + 1, j, k, jj)
1832 if (dwt * dwtp1 >
zero)
then
1833 if (abs(dwt) < abs(dwtp1))
then
1834 dwti = dwti -
half * dwt
1836 dwti = dwti -
half * dwtp1
1840 if (dwt * dwtm1 >
zero)
then
1841 if (abs(dwt) < abs(dwtm1))
then
1842 dwti = dwti +
half * dwt
1844 dwti = dwti +
half * dwtm1
1852 dwti =
w(i + 1, j, k, jj) -
w(i, j, k, jj)
1860 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwti
1885 integer(kind=intType) :: i, j, k, ii
1886 real(kind=realtype) :: rblank
1891 rblank = max(real(
iblank(i, j, k), realtype),
zero)
1915 logical,
intent(in),
optional :: updateDtl
1918 real(kind=realtype),
parameter :: b = 2.0_realtype
1921 real(kind=realtype) :: plim, rlim, clim2
1922 real(kind=realtype) :: cc2, qsi, qsj, qsk, sx, sy, sz, rmu
1923 real(kind=realtype) :: ri, rj, rk, rij, rjk, rki
1924 real(kind=realtype) :: vsi, vsj, vsk, rfl, dpi, dpj, dpk
1925 real(kind=realtype) :: sface, tmp, uux, uuy, uuz
1926 logical :: doScaling, updateDt
1927 integer(kind=intType) :: i, j, k
1930 if (
present(updatedtl))
then
1939 rlim = 0.001_realtype *
rhoinf
1961 uux =
w(i, j, k,
ivx)
1962 uuy =
w(i, j, k,
ivy)
1963 uuz =
w(i, j, k,
ivz)
1964 cc2 =
gamma(i, j, k) *
p(i, j, k) /
w(i, j, k,
irho)
1965 cc2 = max(cc2, clim2)
1976 sx =
si(i - 1, j, k, 1) +
si(i, j, k, 1)
1977 sy =
si(i - 1, j, k, 2) +
si(i, j, k, 2)
1978 sz =
si(i - 1, j, k, 3) +
si(i, j, k, 3)
1980 qsi = uux * sx + uuy * sy + uuz * sz - sface
1982 ri =
half * (abs(qsi) &
1990 sx =
sj(i, j - 1, k, 1) +
sj(i, j, k, 1)
1991 sy =
sj(i, j - 1, k, 2) +
sj(i, j, k, 2)
1992 sz =
sj(i, j - 1, k, 3) +
sj(i, j, k, 3)
1994 qsj = uux * sx + uuy * sy + uuz * sz - sface
1996 rj =
half * (abs(qsj) &
2004 sx =
sk(i, j, k - 1, 1) +
sk(i, j, k, 1)
2005 sy =
sk(i, j, k - 1, 2) +
sk(i, j, k, 2)
2006 sz =
sk(i, j, k - 1, 3) +
sk(i, j, k, 3)
2008 qsk = uux * sx + uuy * sy + uuz * sz - sface
2010 rk =
half * (abs(qsk) &
2015 dtl(i, j, k) = ri + rj + rk
2028 rij = (ri / rj)**
adis
2029 rjk = (rj / rk)**
adis
2030 rki = (rk / ri)**
adis
2047 viscousterm:
if (
viscous)
then
2068 rmu = rmu +
rev(i, j, k)
2074 sx =
si(i, j, k, 1) +
si(i - 1, j, k, 1)
2075 sy =
si(i, j, k, 2) +
si(i - 1, j, k, 2)
2076 sz =
si(i, j, k, 3) +
si(i - 1, j, k, 3)
2078 vsi = rmu * (sx * sx + sy * sy + sz * sz)
2079 dtl(i, j, k) =
dtl(i, j, k) + vsi
2084 sx =
sj(i, j, k, 1) +
sj(i, j - 1, k, 1)
2085 sy =
sj(i, j, k, 2) +
sj(i, j - 1, k, 2)
2086 sz =
sj(i, j, k, 3) +
sj(i, j - 1, k, 3)
2088 vsj = rmu * (sx * sx + sy * sy + sz * sz)
2089 dtl(i, j, k) =
dtl(i, j, k) + vsj
2094 sx =
sk(i, j, k, 1) +
sk(i, j, k - 1, 1)
2095 sy =
sk(i, j, k, 2) +
sk(i, j, k - 1, 2)
2096 sz =
sk(i, j, k, 3) +
sk(i, j, k - 1, 3)
2098 vsk = rmu * (sx * sx + sy * sy + sz * sz)
2099 dtl(i, j, k) =
dtl(i, j, k) + vsk
2120 dtl(i, j, k) =
dtl(i, j, k) + tmp *
vol(i, j, k)
2134 dpi = abs(
p(i + 1, j, k) -
two *
p(i, j, k) +
p(i - 1, j, k)) &
2135 / (
p(i + 1, j, k) +
two *
p(i, j, k) +
p(i - 1, j, k) + plim)
2136 dpj = abs(
p(i, j + 1, k) -
two *
p(i, j, k) +
p(i, j - 1, k)) &
2137 / (
p(i, j + 1, k) +
two *
p(i, j, k) +
p(i, j - 1, k) + plim)
2138 dpk = abs(
p(i, j, k + 1) -
two *
p(i, j, k) +
p(i, j, k - 1)) &
2139 / (
p(i, j, k + 1) +
two *
p(i, j, k) +
p(i, j, k - 1) + plim)
2140 rfl =
one / (
one + b * (dpi + dpj + dpk))
2142 dtl(i, j, k) = rfl /
dtl(i, j, k)
2163 real(kind=realtype) :: qsp, qsm, rqsp, rqsm, porvel, porflux
2164 real(kind=realtype) :: pa, vnp, vnm, fs, sface
2165 integer(kind=intType) :: i, j, k
2166 real(kind=realtype) :: wwx, wwy, wwz, rvol
2179 vnp =
w(i + 1, j, k,
ivx) *
si(i, j, k, 1) &
2180 +
w(i + 1, j, k,
ivy) *
si(i, j, k, 2) &
2181 +
w(i + 1, j, k,
ivz) *
si(i, j, k, 3)
2182 vnm =
w(i, j, k,
ivx) *
si(i, j, k, 1) &
2183 +
w(i, j, k,
ivy) *
si(i, j, k, 2) &
2184 +
w(i, j, k,
ivz) *
si(i, j, k, 3)
2206 porvel = porvel * porflux
2211 qsp = (vnp - sface) * porvel
2212 qsm = (vnm - sface) * porvel
2214 rqsp = qsp *
w(i + 1, j, k,
irho)
2215 rqsm = qsm *
w(i, j, k,
irho)
2221 pa = porflux * (
p(i + 1, j, k) +
p(i, j, k))
2231 fs = rqsp *
w(i + 1, j, k,
ivx) + rqsm *
w(i, j, k,
ivx) &
2232 + pa *
si(i, j, k, 1)
2233 dw(i + 1, j, k,
imx) =
dw(i + 1, j, k,
imx) - fs
2236 fs = rqsp *
w(i + 1, j, k,
ivy) + rqsm *
w(i, j, k,
ivy) &
2237 + pa *
si(i, j, k, 2)
2238 dw(i + 1, j, k,
imy) =
dw(i + 1, j, k,
imy) - fs
2241 fs = rqsp *
w(i + 1, j, k,
ivz) + rqsm *
w(i, j, k,
ivz) &
2242 + pa *
si(i, j, k, 3)
2243 dw(i + 1, j, k,
imz) =
dw(i + 1, j, k,
imz) - fs
2246 fs = qsp *
w(i + 1, j, k,
irhoe) + qsm *
w(i, j, k,
irhoe) &
2247 + porflux * (vnp *
p(i + 1, j, k) + vnm *
p(i, j, k))
2265 vnp =
w(i, j + 1, k,
ivx) *
sj(i, j, k, 1) &
2266 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 2) &
2267 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 3)
2268 vnm =
w(i, j, k,
ivx) *
sj(i, j, k, 1) &
2269 +
w(i, j, k,
ivy) *
sj(i, j, k, 2) &
2270 +
w(i, j, k,
ivz) *
sj(i, j, k, 3)
2293 porvel = porvel * porflux
2298 qsp = (vnp - sface) * porvel
2299 qsm = (vnm - sface) * porvel
2301 rqsp = qsp *
w(i, j + 1, k,
irho)
2302 rqsm = qsm *
w(i, j, k,
irho)
2308 pa = porflux * (
p(i, j + 1, k) +
p(i, j, k))
2318 fs = rqsp *
w(i, j + 1, k,
ivx) + rqsm *
w(i, j, k,
ivx) &
2319 + pa *
sj(i, j, k, 1)
2320 dw(i, j + 1, k,
imx) =
dw(i, j + 1, k,
imx) - fs
2323 fs = rqsp *
w(i, j + 1, k,
ivy) + rqsm *
w(i, j, k,
ivy) &
2324 + pa *
sj(i, j, k, 2)
2325 dw(i, j + 1, k,
imy) =
dw(i, j + 1, k,
imy) - fs
2328 fs = rqsp *
w(i, j + 1, k,
ivz) + rqsm *
w(i, j, k,
ivz) &
2329 + pa *
sj(i, j, k, 3)
2330 dw(i, j + 1, k,
imz) =
dw(i, j + 1, k,
imz) - fs
2333 fs = qsp *
w(i, j + 1, k,
irhoe) + qsm *
w(i, j, k,
irhoe) &
2334 + porflux * (vnp *
p(i, j + 1, k) + vnm *
p(i, j, k))
2352 vnp =
w(i, j, k + 1,
ivx) *
sk(i, j, k, 1) &
2353 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 2) &
2354 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 3)
2355 vnm =
w(i, j, k,
ivx) *
sk(i, j, k, 1) &
2356 +
w(i, j, k,
ivy) *
sk(i, j, k, 2) &
2357 +
w(i, j, k,
ivz) *
sk(i, j, k, 3)
2381 porvel = porvel * porflux
2386 qsp = (vnp - sface) * porvel
2387 qsm = (vnm - sface) * porvel
2389 rqsp = qsp *
w(i, j, k + 1,
irho)
2390 rqsm = qsm *
w(i, j, k,
irho)
2396 pa = porflux * (
p(i, j, k + 1) +
p(i, j, k))
2406 fs = rqsp *
w(i, j, k + 1,
ivx) + rqsm *
w(i, j, k,
ivx) &
2407 + pa *
sk(i, j, k, 1)
2408 dw(i, j, k + 1,
imx) =
dw(i, j, k + 1,
imx) - fs
2411 fs = rqsp *
w(i, j, k + 1,
ivy) + rqsm *
w(i, j, k,
ivy) &
2412 + pa *
sk(i, j, k, 2)
2413 dw(i, j, k + 1,
imy) =
dw(i, j, k + 1,
imy) - fs
2416 fs = rqsp *
w(i, j, k + 1,
ivz) + rqsm *
w(i, j, k,
ivz) &
2417 + pa *
sk(i, j, k, 3)
2418 dw(i, j, k + 1,
imz) =
dw(i, j, k + 1,
imz) - fs
2421 fs = qsp *
w(i, j, k + 1,
irhoe) + qsm *
w(i, j, k,
irhoe) &
2422 + porflux * (vnp *
p(i, j, k + 1) + vnm *
p(i, j, k))
2443 rvol =
w(i, j, k,
irho) *
vol(i, j, k)
2445 + rvol * (wwy *
w(i, j, k,
ivz) - wwz *
w(i, j, k,
ivy))
2447 + rvol * (wwz *
w(i, j, k,
ivx) - wwx *
w(i, j, k,
ivz))
2449 + rvol * (wwx *
w(i, j, k,
ivy) - wwy *
w(i, j, k,
ivx))
2476 real(kind=realtype),
parameter :: dpmax = 0.25_realtype
2477 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
2478 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
2479 real(kind=realtype),
parameter :: omega = 0.5_realtype
2480 real(kind=realtype),
parameter :: oneminomega =
one - omega
2484 integer(kind=intType) :: i, j, k, ind, ii
2486 real(kind=realtype) :: plim, sface
2487 real(kind=realtype) :: sfil, fis2, fis4
2488 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53
2489 real(kind=realtype) :: ppor, rrad, dis2, dis4
2490 real(kind=realtype) :: dp1, dp2, tmp, fs
2491 real(kind=realtype) :: ddw1, ddw2, ddw3, ddw4, ddw5, ddw6
2492 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
2493 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
2494 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
2495 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
2496 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
2497 logical :: correctForK
2530 dss(i, j, k, 1) = abs((
p(i + 1, j, k) -
two *
p(i, j, k) +
p(i - 1, j, k)) &
2531 / (omega * (
p(i + 1, j, k) +
two *
p(i, j, k) +
p(i - 1, j, k)) &
2532 + oneminomega * (abs(
p(i + 1, j, k) -
p(i, j, k)) &
2533 + abs(
p(i, j, k) -
p(i - 1, j, k))) + plim))
2535 dss(i, j, k, 2) = abs((
p(i, j + 1, k) -
two *
p(i, j, k) +
p(i, j - 1, k)) &
2536 / (omega * (
p(i, j + 1, k) +
two *
p(i, j, k) +
p(i, j - 1, k)) &
2537 + oneminomega * (abs(
p(i, j + 1, k) -
p(i, j, k)) &
2538 + abs(
p(i, j, k) -
p(i, j - 1, k))) + plim))
2540 dss(i, j, k, 3) = abs((
p(i, j, k + 1) -
two *
p(i, j, k) +
p(i, j, k - 1)) &
2541 / (omega * (
p(i, j, k + 1) +
two *
p(i, j, k) +
p(i, j, k - 1)) &
2542 + oneminomega * (abs(
p(i, j, k + 1) -
p(i, j, k)) &
2543 + abs(
p(i, j, k) -
p(i, j, k - 1))) + plim))
2558 dis2 = ppor * fis2 * min(dpmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1)))
2559 dis4 = dim(ppor * fis4, dis2)
2564 ddw1 =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
2568 ddw2 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
2571 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivx) &
2574 ddw3 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
2577 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivy) &
2580 ddw4 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
2583 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivz) &
2597 if (correctfork)
then
2598 ddw6 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
2601 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
itu1) &
2611 gm1 = gammaavg -
one
2620 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
2623 area = sqrt(
si(i, j, k, 1)**2 +
si(i, j, k, 2)**2 +
si(i, j, k, 3)**2)
2624 tmp =
one / max(1.e-25_realtype, area)
2625 sx =
si(i, j, k, 1) * tmp
2626 sy =
si(i, j, k, 2) * tmp
2627 sz =
si(i, j, k, 3) * tmp
2629 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2630 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
2632 unavg = uavg * sx + vavg * sy + wavg * sz
2634 ova2avg =
one / a2avg
2639 sface =
sfacei(i, j, k) * tmp
2645 lam1 = abs(unavg - sface + aavg)
2646 lam2 = abs(unavg - sface - aavg)
2647 lam3 = abs(unavg - sface)
2654 lam1 = max(lam1, epsacoustic * rrad) * area
2655 lam2 = max(lam2, epsacoustic * rrad) * area
2656 lam3 = max(lam3, epsshear * rrad) * area
2661 abv1 =
half * (lam1 + lam2)
2662 abv2 =
half * (lam1 - lam2)
2665 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2666 - wavg * drw + dre) - gm53 * drk
2667 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2669 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2670 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2675 fs = lam3 * dr + abv6
2681 fs = lam3 * dru + uavg * abv6 + sx * abv7
2682 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
2687 fs = lam3 * drv + vavg * abv6 + sy * abv7
2688 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
2693 fs = lam3 * drw + wavg * abv6 + sz * abv7
2694 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
2699 fs = lam3 * dre + havg * abv6 + unavg * abv7
2718 dis2 = ppor * fis2 * min(dpmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2)))
2719 dis4 = dim(ppor * fis4, dis2)
2724 ddw1 =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
2728 ddw2 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
2731 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivx) &
2734 ddw3 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
2737 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivy) &
2740 ddw4 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
2743 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivz) &
2757 if (correctfork)
then
2758 ddw6 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
2761 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
itu1) &
2771 gm1 = gammaavg -
one
2780 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
2783 area = sqrt(
sj(i, j, k, 1)**2 +
sj(i, j, k, 2)**2 +
sj(i, j, k, 3)**2)
2784 tmp =
one / max(1.e-25_realtype, area)
2785 sx =
sj(i, j, k, 1) * tmp
2786 sy =
sj(i, j, k, 2) * tmp
2787 sz =
sj(i, j, k, 3) * tmp
2789 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2790 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
2792 unavg = uavg * sx + vavg * sy + wavg * sz
2794 ova2avg =
one / a2avg
2799 sface =
sfacej(i, j, k) * tmp
2805 lam1 = abs(unavg - sface + aavg)
2806 lam2 = abs(unavg - sface - aavg)
2807 lam3 = abs(unavg - sface)
2814 lam1 = max(lam1, epsacoustic * rrad) * area
2815 lam2 = max(lam2, epsacoustic * rrad) * area
2816 lam3 = max(lam3, epsshear * rrad) * area
2821 abv1 =
half * (lam1 + lam2)
2822 abv2 =
half * (lam1 - lam2)
2825 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2826 - wavg * drw + dre) - gm53 * drk
2827 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2829 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2830 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2835 fs = lam3 * dr + abv6
2841 fs = lam3 * dru + uavg * abv6 + sx * abv7
2842 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
2847 fs = lam3 * drv + vavg * abv6 + sy * abv7
2848 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
2853 fs = lam3 * drw + wavg * abv6 + sz * abv7
2854 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
2859 fs = lam3 * dre + havg * abv6 + unavg * abv7
2878 dis2 = ppor * fis2 * min(dpmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3)))
2879 dis4 = dim(ppor * fis4, dis2)
2884 ddw1 =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
2888 ddw2 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
2891 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivx) &
2894 ddw3 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
2897 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivy) &
2900 ddw4 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
2903 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivz) &
2917 if (correctfork)
then
2918 ddw6 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
2921 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
itu1) &
2931 gm1 = gammaavg -
one
2940 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
2943 area = sqrt(
sk(i, j, k, 1)**2 +
sk(i, j, k, 2)**2 +
sk(i, j, k, 3)**2)
2944 tmp =
one / max(1.e-25_realtype, area)
2945 sx =
sk(i, j, k, 1) * tmp
2946 sy =
sk(i, j, k, 2) * tmp
2947 sz =
sk(i, j, k, 3) * tmp
2949 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2950 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
2952 unavg = uavg * sx + vavg * sy + wavg * sz
2954 ova2avg =
one / a2avg
2959 sface =
sfacek(i, j, k) * tmp
2965 lam1 = abs(unavg - sface + aavg)
2966 lam2 = abs(unavg - sface - aavg)
2967 lam3 = abs(unavg - sface)
2974 lam1 = max(lam1, epsacoustic * rrad) * area
2975 lam2 = max(lam2, epsacoustic * rrad) * area
2976 lam3 = max(lam3, epsshear * rrad) * area
2981 abv1 =
half * (lam1 + lam2)
2982 abv2 =
half * (lam1 - lam2)
2985 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2986 - wavg * drw + dre) - gm53 * drk
2987 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2989 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2990 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2995 fs = lam3 * dr + abv6
3001 fs = lam3 * dru + uavg * abv6 + sx * abv7
3002 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
3007 fs = lam3 * drv + vavg * abv6 + sy * abv7
3008 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
3013 fs = lam3 * drw + wavg * abv6 + sz * abv7
3014 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
3019 fs = lam3 * dre + havg * abv6 + unavg * abv7
3044 real(kind=realtype),
parameter :: dssmax = 0.25_realtype
3045 real(kind=realtype) :: sslim, rhoi
3046 real(kind=realtype) :: sfil, fis2, fis4
3047 real(kind=realtype) :: ppor, rrad, dis2, dis4, fs
3048 real(kind=realtype) :: ddw1, ddw2, ddw3, ddw4, ddw5
3049 integer(kind=intType) :: i, j, k
3085 ss(i, j, k) =
p(i, j, k) / (
w(i, j, k,
irho)**
gamma(i, j, k))
3095 dss(i, j, k, 1) = abs((
ss(i + 1, j, k) -
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
3096 / (
ss(i + 1, j, k) +
two *
ss(i, j, k) +
ss(i - 1, j, k) + sslim))
3098 dss(i, j, k, 2) = abs((
ss(i, j + 1, k) -
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
3099 / (
ss(i, j + 1, k) +
two *
ss(i, j, k) +
ss(i, j - 1, k) + sslim))
3101 dss(i, j, k, 3) = abs((
ss(i, j, k + 1) -
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
3102 / (
ss(i, j, k + 1) +
two *
ss(i, j, k) +
ss(i, j, k - 1) + sslim))
3143 rrad = ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
3145 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1)))
3146 dis4 = dim(fis4 * rrad, dis2)
3152 ddw1 =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
3161 ddw2 =
w(i + 1, j, k,
ivx) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
3163 - dis4 * (
w(i + 2, j, k,
ivx) *
w(i + 2, j, k,
irho) - &
3166 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
3171 ddw3 =
w(i + 1, j, k,
ivy) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
3173 - dis4 * (
w(i + 2, j, k,
ivy) *
w(i + 2, j, k,
irho) - &
3176 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
3181 ddw4 =
w(i + 1, j, k,
ivz) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
3183 - dis4 * (
w(i + 2, j, k,
ivz) *
w(i + 2, j, k,
irho) - &
3186 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
3191 ddw5 = (
w(i + 1, j, k,
irhoe) +
p(i + 1, j, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
3193 - dis4 * ((
w(i + 2, j, k,
irhoe) +
p(i + 2, j, k)) - &
3194 (
w(i - 1, j, k,
irhoe) +
p(i - 1, j, k)) -
three * ddw5)
3212 rrad = ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
3214 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2)))
3215 dis4 = dim(fis4 * rrad, dis2)
3221 ddw1 =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
3230 ddw2 =
w(i, j + 1, k,
ivx) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
3232 - dis4 * (
w(i, j + 2, k,
ivx) *
w(i, j + 2, k,
irho) - &
3235 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
3240 ddw3 =
w(i, j + 1, k,
ivy) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
3242 - dis4 * (
w(i, j + 2, k,
ivy) *
w(i, j + 2, k,
irho) - &
3245 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
3250 ddw4 =
w(i, j + 1, k,
ivz) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
3252 - dis4 * (
w(i, j + 2, k,
ivz) *
w(i, j + 2, k,
irho) - &
3255 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
3260 ddw5 = (
w(i, j + 1, k,
irhoe) +
p(i, j + 1, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
3262 - dis4 * ((
w(i, j + 2, k,
irhoe) +
p(i, j + 2, k)) - &
3263 (
w(i, j - 1, k,
irhoe) +
p(i, j - 1, k)) -
three * ddw5)
3281 rrad = ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
3283 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3)))
3284 dis4 = dim(fis4 * rrad, dis2)
3290 ddw1 =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
3299 ddw2 =
w(i, j, k + 1,
ivx) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
3301 - dis4 * (
w(i, j, k + 2,
ivx) *
w(i, j, k + 2,
irho) - &
3304 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
3309 ddw3 =
w(i, j, k + 1,
ivy) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
3311 - dis4 * (
w(i, j, k + 2,
ivy) *
w(i, j, k + 2,
irho) - &
3314 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
3319 ddw4 =
w(i, j, k + 1,
ivz) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
3321 - dis4 * (
w(i, j, k + 2,
ivz) *
w(i, j, k + 2,
irho) - &
3324 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
3329 ddw5 = (
w(i, j, k + 1,
irhoe) +
p(i, j, k + 1)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
3331 - dis4 * ((
w(i, j, k + 2,
irhoe) +
p(i, j, k + 2)) - &
3332 (
w(i, j, k - 1,
irhoe) +
p(i, j, k - 1)) -
three * ddw5)
3366 logical,
intent(in) :: fineGrid
3370 integer(kind=porType) :: por
3372 integer(kind=intType) :: nwInt
3373 integer(kind=intType) :: i, j, k, ind
3374 integer(kind=intType) :: limUsed, riemannUsed
3376 real(kind=realtype) :: sx, sy, sz, omk, opk, sfil, gammaface
3377 real(kind=realtype) :: factminmod, sface
3379 real(kind=realtype),
dimension(nw) :: left, right
3380 real(kind=realtype),
dimension(nw) :: du1, du2, du3
3381 real(kind=realtype),
dimension(nwf) :: flux
3383 logical :: firstOrderK, correctForK, rotationalPeriodic
3389 rotationalperiodic = .true.
3391 rotationalperiodic = .false.
3426 if (finegrid) limused =
limiter
3431 if (finegrid) riemannused =
riemann
3444 if (correctfork)
then
3447 firstorderk = .true.
3450 firstorderk = .false.
3454 firstorderk = .false.
3475 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
3486 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
3489 right(
ivx) =
w(i + 1, j, k,
ivx)
3490 right(
ivy) =
w(i + 1, j, k,
ivy)
3491 right(
ivz) =
w(i + 1, j, k,
ivz)
3492 right(
irhoe) =
p(i + 1, j, k)
3493 if (correctfork) right(
itu1) =
w(i + 1, j, k,
itu1)
3531 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
3542 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
3545 right(
ivx) =
w(i, j + 1, k,
ivx)
3546 right(
ivy) =
w(i, j + 1, k,
ivy)
3547 right(
ivz) =
w(i, j + 1, k,
ivz)
3548 right(
irhoe) =
p(i, j + 1, k)
3549 if (correctfork) right(
itu1) =
w(i, j + 1, k,
itu1)
3586 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
3597 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
3600 right(
ivx) =
w(i, j, k + 1,
ivx)
3601 right(
ivy) =
w(i, j, k + 1,
ivy)
3602 right(
ivz) =
w(i, j, k + 1,
ivz)
3603 right(
irhoe) =
p(i, j, k + 1)
3604 if (correctfork) right(
itu1) =
w(i, j, k + 1,
itu1)
3661 du3(
ivx) =
w(i + 2, j, k,
ivx) -
w(i + 1, j, k,
ivx)
3665 du3(
ivy) =
w(i + 2, j, k,
ivy) -
w(i + 1, j, k,
ivy)
3669 du3(
ivz) =
w(i + 2, j, k,
ivz) -
w(i + 1, j, k,
ivz)
3671 du1(
irhoe) =
p(i, j, k) -
p(i - 1, j, k)
3672 du2(
irhoe) =
p(i + 1, j, k) -
p(i, j, k)
3673 du3(
irhoe) =
p(i + 2, j, k) -
p(i + 1, j, k)
3675 if (correctfork)
then
3697 right(
ivx) = right(
ivx) +
w(i + 1, j, k,
ivx)
3698 right(
ivy) = right(
ivy) +
w(i + 1, j, k,
ivy)
3699 right(
ivz) = right(
ivz) +
w(i + 1, j, k,
ivz)
3702 if (correctfork)
then
3710 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
3756 du3(
ivx) =
w(i, j + 2, k,
ivx) -
w(i, j + 1, k,
ivx)
3760 du3(
ivy) =
w(i, j + 2, k,
ivy) -
w(i, j + 1, k,
ivy)
3764 du3(
ivz) =
w(i, j + 2, k,
ivz) -
w(i, j + 1, k,
ivz)
3766 du1(
irhoe) =
p(i, j, k) -
p(i, j - 1, k)
3767 du2(
irhoe) =
p(i, j + 1, k) -
p(i, j, k)
3768 du3(
irhoe) =
p(i, j + 2, k) -
p(i, j + 1, k)
3770 if (correctfork)
then
3792 right(
ivx) = right(
ivx) +
w(i, j + 1, k,
ivx)
3793 right(
ivy) = right(
ivy) +
w(i, j + 1, k,
ivy)
3794 right(
ivz) = right(
ivz) +
w(i, j + 1, k,
ivz)
3797 if (correctfork)
then
3805 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
3850 du3(
ivx) =
w(i, j, k + 2,
ivx) -
w(i, j, k + 1,
ivx)
3854 du3(
ivy) =
w(i, j, k + 2,
ivy) -
w(i, j, k + 1,
ivy)
3858 du3(
ivz) =
w(i, j, k + 2,
ivz) -
w(i, j, k + 1,
ivz)
3860 du1(
irhoe) =
p(i, j, k) -
p(i, j, k - 1)
3861 du2(
irhoe) =
p(i, j, k + 1) -
p(i, j, k)
3862 du3(
irhoe) =
p(i, j, k + 2) -
p(i, j, k + 1)
3864 if (correctfork)
then
3886 right(
ivx) = right(
ivx) +
w(i, j, k + 1,
ivx)
3887 right(
ivy) = right(
ivy) +
w(i, j, k + 1,
ivy)
3888 right(
ivz) = right(
ivz) +
w(i, j, k + 1,
ivz)
3891 if (correctfork)
then
3899 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
3949 real(kind=realtype),
parameter :: epslim = 1.e-10_realtype
3953 real(kind=realtype),
dimension(:),
intent(inout) :: du1, du2, du3
3954 real(kind=realtype),
dimension(:),
intent(out) :: left, right
3956 real(kind=realtype),
dimension(:, :, :, :, :),
pointer :: rotmatrix
3960 integer(kind=intType) :: l
3962 real(kind=realtype) :: rl1, rl2, rr1, rr2, tmp, dvx, dvy, dvz
3964 real(kind=realtype),
dimension(3, 3) :: rot
3969 if (rotationalperiodic)
then
3974 rot(1, 1) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1, 1)
3975 rot(1, 2) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1, 2)
3976 rot(1, 3) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1, 3)
3978 rot(2, 1) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 2, 1)
3979 rot(2, 2) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 2, 2)
3980 rot(2, 3) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 2, 3)
3982 rot(3, 1) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 3, 1)
3983 rot(3, 2) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 3, 2)
3984 rot(3, 3) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 3, 3)
3989 dvx = du1(
ivx); dvy = du1(
ivy); dvz = du1(
ivz)
3990 du1(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
3991 du1(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
3992 du1(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
3994 dvx = du2(
ivx); dvy = du2(
ivy); dvz = du2(
ivz)
3995 du2(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
3996 du2(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
3997 du2(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
3999 dvx = du3(
ivx); dvy = du3(
ivy); dvz = du3(
ivz)
4000 du3(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
4001 du3(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
4002 du3(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
4008 select case (limused)
4016 left(l) = omk * du1(l) + opk * du2(l)
4017 right(l) = -omk * du3(l) - opk * du2(l)
4032 tmp =
one / sign(max(abs(du2(l)), epslim), du2(l))
4034 du2(l) / sign(max(abs(du1(l)), epslim), du1(l)))
4035 rl2 = max(
zero, du1(l) * tmp)
4037 rr1 = max(
zero, du3(l) * tmp)
4039 du2(l) / sign(max(abs(du3(l)), epslim), du3(l)))
4043 rl1 = rl1 * (rl1 +
one) / (rl1 * rl1 +
one)
4044 rl2 = rl2 * (rl2 +
one) / (rl2 * rl2 +
one)
4045 rr1 = rr1 * (rr1 +
one) / (rr1 * rr1 +
one)
4046 rr2 = rr2 * (rr2 +
one) / (rr2 * rr2 +
one)
4051 left(l) = omk * rl1 * du1(l) + opk * rl2 * du2(l)
4052 right(l) = -opk * rr1 * du2(l) - omk * rr2 * du3(l)
4068 tmp =
one / sign(max(abs(du2(l)), epslim), du2(l))
4070 du2(l) / sign(max(abs(du1(l)), epslim), du1(l)))
4071 rl2 = max(
zero, du1(l) * tmp)
4073 rr1 = max(
zero, du3(l) * tmp)
4075 du2(l) / sign(max(abs(du3(l)), epslim), du3(l)))
4079 rl1 = min(
one, factminmod * rl1)
4080 rl2 = min(
one, factminmod * rl2)
4081 rr1 = min(
one, factminmod * rr1)
4082 rr2 = min(
one, factminmod * rr2)
4087 left(l) = omk * rl1 * du1(l) + opk * rl2 * du2(l)
4088 right(l) = -opk * rr1 * du2(l) - omk * rr2 * du3(l)
4098 if (firstorderk)
then
4107 if (rotationalperiodic)
then
4111 dvx = left(
ivx); dvy = left(
ivy); dvz = left(
ivz)
4112 left(
ivx) = rot(1, 1) * dvx + rot(2, 1) * dvy + rot(3, 1) * dvz
4113 left(
ivy) = rot(1, 2) * dvx + rot(2, 2) * dvy + rot(3, 2) * dvz
4114 left(
ivz) = rot(1, 3) * dvx + rot(2, 3) * dvy + rot(3, 3) * dvz
4118 dvx = right(
ivx); dvy = right(
ivy); dvz = right(
ivz)
4119 right(
ivx) = rot(1, 1) * dvx + rot(2, 1) * dvy + rot(3, 1) * dvz
4120 right(
ivy) = rot(1, 2) * dvx + rot(2, 2) * dvy + rot(3, 2) * dvz
4121 right(
ivz) = rot(1, 3) * dvx + rot(2, 3) * dvy + rot(3, 3) * dvz
4138 real(kind=realtype),
dimension(*),
intent(in) :: left, right
4139 real(kind=realtype),
dimension(*),
intent(out) :: flux
4143 real(kind=realtype) :: porflux, rface
4144 real(kind=realtype) :: etl, etr, z1l, z1r, tmp
4145 real(kind=realtype) :: dr, dru, drv, drw, dre, drk
4146 real(kind=realtype) :: ravg, uavg, vavg, wavg, havg, kavg
4147 real(kind=realtype) :: alphaavg, a2avg, aavg, unavg
4148 real(kind=realtype) :: ovaavg, ova2avg, area, eta
4149 real(kind=realtype) :: gm1, gm53
4150 real(kind=realtype) :: lam1, lam2, lam3
4151 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
4152 real(kind=realtype),
dimension(2) :: ktmp
4162 gm1 = gammaface -
one
4167 select case (riemannused)
4183 z1l = sqrt(left(
irho))
4184 z1r = sqrt(right(
irho))
4185 tmp =
one / (z1l + z1r)
4190 if (correctfork)
then
4195 ktmp(1) = left(
itu1)
4196 ktmp(2) = right(
itu1)
4206 kavg = tmp * (z1l * left(
itu1) + z1r * right(
itu1))
4221 left(
irhoe), ktmp(1), etl, correctfork)
4224 right(
irhoe), ktmp(2), etr, correctfork)
4238 ravg =
fourth * (z1r + z1l)**2
4239 uavg = tmp * (z1l * left(
ivx) + z1r * right(
ivx))
4240 vavg = tmp * (z1l * left(
ivy) + z1r * right(
ivy))
4241 wavg = tmp * (z1l * left(
ivz) + z1r * right(
ivz))
4242 havg = tmp * ((etl + left(
irhoe)) / z1l &
4243 + (etr + right(
irhoe)) / z1r)
4248 area = sqrt(sx**2 + sy**2 + sz**2)
4249 tmp =
one / max(1.e-25_realtype, area)
4258 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4259 a2avg = abs(gm1 * (havg - alphaavg) - gm53 * kavg)
4261 unavg = uavg * sx + vavg * sy + wavg * sz
4264 ova2avg =
one / a2avg
4282 eta =
half * (abs((left(
ivx) - right(
ivx)) * sx &
4283 + (left(
ivy) - right(
ivy)) * sy &
4284 + (left(
ivz) - right(
ivz)) * sz) &
4285 + abs(sqrt(gammaface * left(
irhoe) / left(
irho)) &
4286 - sqrt(gammaface * right(
irhoe) / right(
irho))))
4290 lam1 = abs(unavg - rface + aavg)
4291 lam2 = abs(unavg - rface - aavg)
4292 lam3 = abs(unavg - rface)
4297 if (lam1 < tmp) lam1 = eta +
fourth * lam1 * lam1 / eta
4298 if (lam2 < tmp) lam2 = eta +
fourth * lam2 * lam2 / eta
4299 if (lam3 < tmp) lam3 = eta +
fourth * lam3 * lam3 / eta
4311 abv1 =
half * (lam1 + lam2)
4312 abv2 =
half * (lam1 - lam2)
4315 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4316 - wavg * drw + dre) - gm53 * drk
4317 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4319 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4320 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4326 flux(
irho) = -porflux * (lam3 * dr + abv6)
4327 flux(
imx) = -porflux * (lam3 * dru + uavg * abv6 &
4329 flux(
imy) = -porflux * (lam3 * drv + vavg * abv6 &
4331 flux(
imz) = -porflux * (lam3 * drw + wavg * abv6 &
4333 flux(
irhoe) = -porflux * (lam3 * dre + havg * abv6 &
4347 "Turkel preconditioner not implemented yet")
4351 "choi merkle preconditioner not implemented yet")
4356 call terminate(
"riemannFlux",
"van leer fvs not implemented yet")
4359 call terminate(
"riemannFlux",
"ausmdv fvs not implemented yet")
4382 real(kind=realtype),
parameter :: dssmax = 0.25_realtype
4383 real(kind=realtype) :: sslim, rhoi
4384 real(kind=realtype) :: sfil, fis2, fis4
4385 real(kind=realtype) :: ppor, rrad, dis2, dis4, fs
4386 real(kind=realtype) :: ddw
4387 integer(kind=intType) :: i, j, k
4415 dss(i, j, k, 1) = abs((
ss(i + 1, j, k) -
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
4416 / (
ss(i + 1, j, k) +
two *
ss(i, j, k) +
ss(i - 1, j, k) + sslim))
4418 dss(i, j, k, 2) = abs((
ss(i, j + 1, k) -
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
4419 / (
ss(i, j + 1, k) +
two *
ss(i, j, k) +
ss(i, j - 1, k) + sslim))
4421 dss(i, j, k, 3) = abs((
ss(i, j, k + 1) -
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
4422 / (
ss(i, j, k + 1) +
two *
ss(i, j, k) +
ss(i, j, k - 1) + sslim))
4453 rrad = ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
4455 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1))) +
sigma * fis4 * rrad
4461 ddw =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
4469 ddw =
w(i + 1, j, k,
ivx) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
4472 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
4477 ddw =
w(i + 1, j, k,
ivy) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
4480 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
4484 ddw =
w(i + 1, j, k,
ivz) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
4487 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
4491 ddw = (
w(i + 1, j, k,
irhoe) +
p(i + 1, j, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
4511 rrad = ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
4513 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2))) +
sigma * fis4 * rrad
4519 ddw =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
4527 ddw =
w(i, j + 1, k,
ivx) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
4530 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
4535 ddw =
w(i, j + 1, k,
ivy) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
4538 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
4543 ddw =
w(i, j + 1, k,
ivz) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
4546 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
4551 ddw = (
w(i, j + 1, k,
irhoe) +
p(i, j + 1, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
4570 rrad = ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
4572 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3))) +
sigma * fis4 * rrad
4578 ddw =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
4586 ddw =
w(i, j, k + 1,
ivx) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
4589 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
4594 ddw =
w(i, j, k + 1,
ivy) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
4597 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
4602 ddw =
w(i, j, k + 1,
ivz) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
4605 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
4609 ddw = (
w(i, j, k + 1,
irhoe) +
p(i, j, k + 1)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
4638 real(kind=realtype),
parameter :: dpmax = 0.25_realtype
4639 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
4640 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
4641 real(kind=realtype),
parameter :: omega = 0.5_realtype
4642 real(kind=realtype),
parameter :: oneminomega =
one - omega
4646 integer(kind=intType) :: i, j, k, ind, ii
4648 real(kind=realtype) :: plim, sface
4649 real(kind=realtype) :: sfil, fis2, fis4
4650 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53
4651 real(kind=realtype) :: ppor, rrad, dis2, dis4
4652 real(kind=realtype) :: dp1, dp2, tmp, fs
4653 real(kind=realtype) :: ddw, ddw6
4654 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
4655 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
4656 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
4657 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
4658 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
4659 logical :: correctForK
4692 dss(i, j, k, 1) = abs((
ss(i + 1, j, k) -
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
4693 / (omega * (
ss(i + 1, j, k) +
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
4694 + oneminomega * (abs(
ss(i + 1, j, k) -
ss(i, j, k)) &
4695 + abs(
ss(i, j, k) -
ss(i - 1, j, k))) + plim))
4697 dss(i, j, k, 2) = abs((
ss(i, j + 1, k) -
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
4698 / (omega * (
ss(i, j + 1, k) +
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
4699 + oneminomega * (abs(
ss(i, j + 1, k) -
ss(i, j, k)) &
4700 + abs(
ss(i, j, k) -
ss(i, j - 1, k))) + plim))
4702 dss(i, j, k, 3) = abs((
ss(i, j, k + 1) -
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
4703 / (omega * (
ss(i, j, k + 1) +
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
4704 + oneminomega * (abs(
ss(i, j, k + 1) -
ss(i, j, k)) &
4705 + abs(
ss(i, j, k) -
ss(i, j, k - 1))) + plim))
4721 dis2 = fis2 * ppor * min(dpmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1))) &
4722 +
sigma * fis4 * ppor
4727 ddw =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
4730 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
4734 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
4738 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
4752 if (correctfork)
then
4753 ddw6 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
4756 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
itu1) &
4766 gm1 = gammaavg -
one
4775 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
4778 area = sqrt(
si(i, j, k, 1)**2 +
si(i, j, k, 2)**2 +
si(i, j, k, 3)**2)
4779 tmp =
one / max(1.e-25_realtype, area)
4780 sx =
si(i, j, k, 1) * tmp
4781 sy =
si(i, j, k, 2) * tmp
4782 sz =
si(i, j, k, 3) * tmp
4784 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4785 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4787 unavg = uavg * sx + vavg * sy + wavg * sz
4789 ova2avg =
one / a2avg
4794 sface =
sfacei(i, j, k) * tmp
4800 lam1 = abs(unavg - sface + aavg)
4801 lam2 = abs(unavg - sface - aavg)
4802 lam3 = abs(unavg - sface)
4809 lam1 = max(lam1, epsacoustic * rrad) * area
4810 lam2 = max(lam2, epsacoustic * rrad) * area
4811 lam3 = max(lam3, epsshear * rrad) * area
4816 abv1 =
half * (lam1 + lam2)
4817 abv2 =
half * (lam1 - lam2)
4820 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4821 - wavg * drw + dre) - gm53 * drk
4822 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4824 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4825 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4830 fs = lam3 * dr + abv6
4836 fs = lam3 * dru + uavg * abv6 + sx * abv7
4837 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
4842 fs = lam3 * drv + vavg * abv6 + sy * abv7
4843 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
4848 fs = lam3 * drw + wavg * abv6 + sz * abv7
4849 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
4854 fs = lam3 * dre + havg * abv6 + unavg * abv7
4873 dis2 = fis2 * ppor * min(dpmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2))) &
4874 +
sigma * fis4 * ppor
4879 ddw =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
4882 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
4886 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
4890 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
4904 if (correctfork)
then
4905 ddw6 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
4908 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
itu1) &
4918 gm1 = gammaavg -
one
4927 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
4930 area = sqrt(
sj(i, j, k, 1)**2 +
sj(i, j, k, 2)**2 +
sj(i, j, k, 3)**2)
4931 tmp =
one / max(1.e-25_realtype, area)
4932 sx =
sj(i, j, k, 1) * tmp
4933 sy =
sj(i, j, k, 2) * tmp
4934 sz =
sj(i, j, k, 3) * tmp
4936 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4937 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4939 unavg = uavg * sx + vavg * sy + wavg * sz
4941 ova2avg =
one / a2avg
4946 sface =
sfacej(i, j, k) * tmp
4952 lam1 = abs(unavg - sface + aavg)
4953 lam2 = abs(unavg - sface - aavg)
4954 lam3 = abs(unavg - sface)
4961 lam1 = max(lam1, epsacoustic * rrad) * area
4962 lam2 = max(lam2, epsacoustic * rrad) * area
4963 lam3 = max(lam3, epsshear * rrad) * area
4968 abv1 =
half * (lam1 + lam2)
4969 abv2 =
half * (lam1 - lam2)
4972 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4973 - wavg * drw + dre) - gm53 * drk
4974 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4976 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4977 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4982 fs = lam3 * dr + abv6
4988 fs = lam3 * dru + uavg * abv6 + sx * abv7
4989 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
4994 fs = lam3 * drv + vavg * abv6 + sy * abv7
4995 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
5000 fs = lam3 * drw + wavg * abv6 + sz * abv7
5001 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
5006 fs = lam3 * dre + havg * abv6 + unavg * abv7
5025 dis2 = fis2 * ppor * min(dpmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3))) &
5026 +
sigma * fis4 * ppor
5031 ddw =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
5034 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
5038 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
5042 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
5056 if (correctfork)
then
5057 ddw6 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
5060 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
itu1) &
5070 gm1 = gammaavg -
one
5079 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
5082 area = sqrt(
sk(i, j, k, 1)**2 +
sk(i, j, k, 2)**2 +
sk(i, j, k, 3)**2)
5083 tmp =
one / max(1.e-25_realtype, area)
5084 sx =
sk(i, j, k, 1) * tmp
5085 sy =
sk(i, j, k, 2) * tmp
5086 sz =
sk(i, j, k, 3) * tmp
5088 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
5089 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
5091 unavg = uavg * sx + vavg * sy + wavg * sz
5093 ova2avg =
one / a2avg
5098 sface =
sfacek(i, j, k) * tmp
5104 lam1 = abs(unavg - sface + aavg)
5105 lam2 = abs(unavg - sface - aavg)
5106 lam3 = abs(unavg - sface)
5113 lam1 = max(lam1, epsacoustic * rrad) * area
5114 lam2 = max(lam2, epsacoustic * rrad) * area
5115 lam3 = max(lam3, epsshear * rrad) * area
5120 abv1 =
half * (lam1 + lam2)
5121 abv2 =
half * (lam1 - lam2)
5124 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
5125 - wavg * drw + dre) - gm53 * drk
5126 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
5128 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
5129 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
5134 fs = lam3 * dr + abv6
5140 fs = lam3 * dru + uavg * abv6 + sx * abv7
5141 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
5146 fs = lam3 * drv + vavg * abv6 + sy * abv7
5147 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
5152 fs = lam3 * drw + wavg * abv6 + sz * abv7
5153 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
5158 fs = lam3 * dre + havg * abv6 + unavg * abv7
5177 logical :: correctForK
5178 real(kind=realtype) :: pp
5179 real(kind=realtype),
parameter :: twothird =
two *
third
5180 integer(kind=intType) :: i, j, k
5185 if (correctfork)
then
5189 pp =
p(i, j, k) - twothird *
w(i, j, k,
irho) *
w(i, j, k,
itu1)
5198 aa(i, j, k) =
gamma(i, j, k) *
p(i, j, k) /
w(i, j, k,
irho)
5214 real(kind=realtype) :: a2, ovol, ubar, vbar, wbar, sx, sy, sz
5215 integer(kind=intType) :: i, j, k
5247 sx =
sk(i, j, k - 1, 1) +
sk(i + 1, j, k - 1, 1) &
5248 +
sk(i, j + 1, k - 1, 1) +
sk(i + 1, j + 1, k - 1, 1) &
5249 +
sk(i, j, k, 1) +
sk(i + 1, j, k, 1) &
5250 +
sk(i, j + 1, k, 1) +
sk(i + 1, j + 1, k, 1)
5251 sy =
sk(i, j, k - 1, 2) +
sk(i + 1, j, k - 1, 2) &
5252 +
sk(i, j + 1, k - 1, 2) +
sk(i + 1, j + 1, k - 1, 2) &
5253 +
sk(i, j, k, 2) +
sk(i + 1, j, k, 2) &
5254 +
sk(i, j + 1, k, 2) +
sk(i + 1, j + 1, k, 2)
5255 sz =
sk(i, j, k - 1, 3) +
sk(i + 1, j, k - 1, 3) &
5256 +
sk(i, j + 1, k - 1, 3) +
sk(i + 1, j + 1, k - 1, 3) &
5257 +
sk(i, j, k, 3) +
sk(i + 1, j, k, 3) &
5258 +
sk(i, j + 1, k, 3) +
sk(i + 1, j + 1, k, 3)
5265 +
w(i, j + 1, k,
ivx) +
w(i + 1, j + 1, k,
ivx))
5267 +
w(i, j + 1, k,
ivy) +
w(i + 1, j + 1, k,
ivy))
5269 +
w(i, j + 1, k,
ivz) +
w(i + 1, j + 1, k,
ivz))
5271 a2 =
fourth * (
aa(i, j, k) +
aa(i + 1, j, k) +
aa(i, j + 1, k) +
aa(i + 1, j + 1, k))
5279 ux(i, j, k - 1) =
ux(i, j, k - 1) + ubar * sx
5280 uy(i, j, k - 1) =
uy(i, j, k - 1) + ubar * sy
5281 uz(i, j, k - 1) =
uz(i, j, k - 1) + ubar * sz
5283 vx(i, j, k - 1) =
vx(i, j, k - 1) + vbar * sx
5284 vy(i, j, k - 1) =
vy(i, j, k - 1) + vbar * sy
5285 vz(i, j, k - 1) =
vz(i, j, k - 1) + vbar * sz
5287 wx(i, j, k - 1) =
wx(i, j, k - 1) + wbar * sx
5288 wy(i, j, k - 1) =
wy(i, j, k - 1) + wbar * sy
5289 wz(i, j, k - 1) =
wz(i, j, k - 1) + wbar * sz
5291 qx(i, j, k - 1) =
qx(i, j, k - 1) - a2 * sx
5292 qy(i, j, k - 1) =
qy(i, j, k - 1) - a2 * sy
5293 qz(i, j, k - 1) =
qz(i, j, k - 1) - a2 * sz
5297 ux(i, j, k) =
ux(i, j, k) - ubar * sx
5298 uy(i, j, k) =
uy(i, j, k) - ubar * sy
5299 uz(i, j, k) =
uz(i, j, k) - ubar * sz
5301 vx(i, j, k) =
vx(i, j, k) - vbar * sx
5302 vy(i, j, k) =
vy(i, j, k) - vbar * sy
5303 vz(i, j, k) =
vz(i, j, k) - vbar * sz
5305 wx(i, j, k) =
wx(i, j, k) - wbar * sx
5306 wy(i, j, k) =
wy(i, j, k) - wbar * sy
5307 wz(i, j, k) =
wz(i, j, k) - wbar * sz
5309 qx(i, j, k) =
qx(i, j, k) + a2 * sx
5310 qy(i, j, k) =
qy(i, j, k) + a2 * sy
5311 qz(i, j, k) =
qz(i, j, k) + a2 * sz
5329 sx =
sj(i, j - 1, k, 1) +
sj(i + 1, j - 1, k, 1) &
5330 +
sj(i, j - 1, k + 1, 1) +
sj(i + 1, j - 1, k + 1, 1) &
5331 +
sj(i, j, k, 1) +
sj(i + 1, j, k, 1) &
5332 +
sj(i, j, k + 1, 1) +
sj(i + 1, j, k + 1, 1)
5333 sy =
sj(i, j - 1, k, 2) +
sj(i + 1, j - 1, k, 2) &
5334 +
sj(i, j - 1, k + 1, 2) +
sj(i + 1, j - 1, k + 1, 2) &
5335 +
sj(i, j, k, 2) +
sj(i + 1, j, k, 2) &
5336 +
sj(i, j, k + 1, 2) +
sj(i + 1, j, k + 1, 2)
5337 sz =
sj(i, j - 1, k, 3) +
sj(i + 1, j - 1, k, 3) &
5338 +
sj(i, j - 1, k + 1, 3) +
sj(i + 1, j - 1, k + 1, 3) &
5339 +
sj(i, j, k, 3) +
sj(i + 1, j, k, 3) &
5340 +
sj(i, j, k + 1, 3) +
sj(i + 1, j, k + 1, 3)
5347 +
w(i, j, k + 1,
ivx) +
w(i + 1, j, k + 1,
ivx))
5349 +
w(i, j, k + 1,
ivy) +
w(i + 1, j, k + 1,
ivy))
5351 +
w(i, j, k + 1,
ivz) +
w(i + 1, j, k + 1,
ivz))
5353 a2 =
fourth * (
aa(i, j, k) +
aa(i + 1, j, k) +
aa(i, j, k + 1) +
aa(i + 1, j, k + 1))
5361 ux(i, j - 1, k) =
ux(i, j - 1, k) + ubar * sx
5362 uy(i, j - 1, k) =
uy(i, j - 1, k) + ubar * sy
5363 uz(i, j - 1, k) =
uz(i, j - 1, k) + ubar * sz
5365 vx(i, j - 1, k) =
vx(i, j - 1, k) + vbar * sx
5366 vy(i, j - 1, k) =
vy(i, j - 1, k) + vbar * sy
5367 vz(i, j - 1, k) =
vz(i, j - 1, k) + vbar * sz
5369 wx(i, j - 1, k) =
wx(i, j - 1, k) + wbar * sx
5370 wy(i, j - 1, k) =
wy(i, j - 1, k) + wbar * sy
5371 wz(i, j - 1, k) =
wz(i, j - 1, k) + wbar * sz
5373 qx(i, j - 1, k) =
qx(i, j - 1, k) - a2 * sx
5374 qy(i, j - 1, k) =
qy(i, j - 1, k) - a2 * sy
5375 qz(i, j - 1, k) =
qz(i, j - 1, k) - a2 * sz
5379 ux(i, j, k) =
ux(i, j, k) - ubar * sx
5380 uy(i, j, k) =
uy(i, j, k) - ubar * sy
5381 uz(i, j, k) =
uz(i, j, k) - ubar * sz
5383 vx(i, j, k) =
vx(i, j, k) - vbar * sx
5384 vy(i, j, k) =
vy(i, j, k) - vbar * sy
5385 vz(i, j, k) =
vz(i, j, k) - vbar * sz
5387 wx(i, j, k) =
wx(i, j, k) - wbar * sx
5388 wy(i, j, k) =
wy(i, j, k) - wbar * sy
5389 wz(i, j, k) =
wz(i, j, k) - wbar * sz
5391 qx(i, j, k) =
qx(i, j, k) + a2 * sx
5392 qy(i, j, k) =
qy(i, j, k) + a2 * sy
5393 qz(i, j, k) =
qz(i, j, k) + a2 * sz
5411 sx =
si(i - 1, j, k, 1) +
si(i - 1, j + 1, k, 1) &
5412 +
si(i - 1, j, k + 1, 1) +
si(i - 1, j + 1, k + 1, 1) &
5413 +
si(i, j, k, 1) +
si(i, j + 1, k, 1) &
5414 +
si(i, j, k + 1, 1) +
si(i, j + 1, k + 1, 1)
5415 sy =
si(i - 1, j, k, 2) +
si(i - 1, j + 1, k, 2) &
5416 +
si(i - 1, j, k + 1, 2) +
si(i - 1, j + 1, k + 1, 2) &
5417 +
si(i, j, k, 2) +
si(i, j + 1, k, 2) &
5418 +
si(i, j, k + 1, 2) +
si(i, j + 1, k + 1, 2)
5419 sz =
si(i - 1, j, k, 3) +
si(i - 1, j + 1, k, 3) &
5420 +
si(i - 1, j, k + 1, 3) +
si(i - 1, j + 1, k + 1, 3) &
5421 +
si(i, j, k, 3) +
si(i, j + 1, k, 3) &
5422 +
si(i, j, k + 1, 3) +
si(i, j + 1, k + 1, 3)
5429 +
w(i, j, k + 1,
ivx) +
w(i, j + 1, k + 1,
ivx))
5431 +
w(i, j, k + 1,
ivy) +
w(i, j + 1, k + 1,
ivy))
5433 +
w(i, j, k + 1,
ivz) +
w(i, j + 1, k + 1,
ivz))
5435 a2 =
fourth * (
aa(i, j, k) +
aa(i, j + 1, k) +
aa(i, j, k + 1) +
aa(i, j + 1, k + 1))
5443 ux(i - 1, j, k) =
ux(i - 1, j, k) + ubar * sx
5444 uy(i - 1, j, k) =
uy(i - 1, j, k) + ubar * sy
5445 uz(i - 1, j, k) =
uz(i - 1, j, k) + ubar * sz
5447 vx(i - 1, j, k) =
vx(i - 1, j, k) + vbar * sx
5448 vy(i - 1, j, k) =
vy(i - 1, j, k) + vbar * sy
5449 vz(i - 1, j, k) =
vz(i - 1, j, k) + vbar * sz
5451 wx(i - 1, j, k) =
wx(i - 1, j, k) + wbar * sx
5452 wy(i - 1, j, k) =
wy(i - 1, j, k) + wbar * sy
5453 wz(i - 1, j, k) =
wz(i - 1, j, k) + wbar * sz
5455 qx(i - 1, j, k) =
qx(i - 1, j, k) - a2 * sx
5456 qy(i - 1, j, k) =
qy(i - 1, j, k) - a2 * sy
5457 qz(i - 1, j, k) =
qz(i - 1, j, k) - a2 * sz
5461 ux(i, j, k) =
ux(i, j, k) - ubar * sx
5462 uy(i, j, k) =
uy(i, j, k) - ubar * sy
5463 uz(i, j, k) =
uz(i, j, k) - ubar * sz
5465 vx(i, j, k) =
vx(i, j, k) - vbar * sx
5466 vy(i, j, k) =
vy(i, j, k) - vbar * sy
5467 vz(i, j, k) =
vz(i, j, k) - vbar * sz
5469 wx(i, j, k) =
wx(i, j, k) - wbar * sx
5470 wy(i, j, k) =
wy(i, j, k) - wbar * sy
5471 wz(i, j, k) =
wz(i, j, k) - wbar * sz
5473 qx(i, j, k) =
qx(i, j, k) + a2 * sx
5474 qy(i, j, k) =
qy(i, j, k) + a2 * sy
5475 qz(i, j, k) =
qz(i, j, k) + a2 * sz
5489 ovol =
one / (
vol(i, j, k) +
vol(i, j, k + 1) &
5490 +
vol(i + 1, j, k) +
vol(i + 1, j, k + 1) &
5491 +
vol(i, j + 1, k) +
vol(i, j + 1, k + 1) &
5492 +
vol(i + 1, j + 1, k) +
vol(i + 1, j + 1, k + 1))
5497 ux(i, j, k) =
ux(i, j, k) * ovol
5498 uy(i, j, k) =
uy(i, j, k) * ovol
5499 uz(i, j, k) =
uz(i, j, k) * ovol
5501 vx(i, j, k) =
vx(i, j, k) * ovol
5502 vy(i, j, k) =
vy(i, j, k) * ovol
5503 vz(i, j, k) =
vz(i, j, k) * ovol
5505 wx(i, j, k) =
wx(i, j, k) * ovol
5506 wy(i, j, k) =
wy(i, j, k) * ovol
5507 wz(i, j, k) =
wz(i, j, k) * ovol
5509 qx(i, j, k) =
qx(i, j, k) * ovol
5510 qy(i, j, k) =
qy(i, j, k) * ovol
5511 qz(i, j, k) =
qz(i, j, k) * ovol
5533 logical,
intent(in),
optional :: storeWallTensor
5536 real(kind=realtype) :: rfilv, por, mul, mue, mut, heatcoef
5537 real(kind=realtype) :: gm1, factlamheat, factturbheat
5538 real(kind=realtype) :: u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z
5539 real(kind=realtype) :: q_x, q_y, q_z
5540 real(kind=realtype) :: corr, ssx, ssy, ssz, fracdiv, snrm
5541 real(kind=realtype) :: tauxx, tauyy, tauzz
5542 real(kind=realtype) :: tauxy, tauxz, tauyz
5543 real(kind=realtype) :: tauxxs, tauyys, tauzzs
5544 real(kind=realtype) :: tauxys, tauxzs, tauyzs
5545 real(kind=realtype) :: ubar, vbar, wbar
5546 real(kind=realtype) :: exx, eyy, ezz
5547 real(kind=realtype) :: exy, exz, eyz
5548 real(kind=realtype) :: wxx, wyy, wzz
5549 real(kind=realtype) :: wxy, wxz, wyz, wyx, wzx, wzy
5550 real(kind=realtype) :: den, ccr1
5551 real(kind=realtype) :: fmx, fmy, fmz, frhoe, fact
5552 integer(kind=intType) :: i, j, k, io, jo, ko
5553 real(kind=realtype),
parameter :: xminn = 1.e-10_realtype
5554 real(kind=realtype),
parameter :: twothird =
two *
third
5555 real(kind=realtype),
dimension(9, 2:max(il, jl), 2:max(jl, kl), 2) :: tmpstore
5557 logical :: storeWall
5560 if (
present(storewalltensor))
then
5561 storewall = storewalltensor
5592 mul = por * (
rlv(i, j, k) +
rlv(i, j, k + 1))
5593 mue = por * (
rev(i, j, k) +
rev(i, j, k + 1))
5600 heatcoef = mul * factlamheat + mue * factturbheat
5605 u_x =
fourth * (
ux(i - 1, j - 1, k) +
ux(i, j - 1, k) &
5606 +
ux(i - 1, j, k) +
ux(i, j, k))
5607 u_y =
fourth * (
uy(i - 1, j - 1, k) +
uy(i, j - 1, k) &
5608 +
uy(i - 1, j, k) +
uy(i, j, k))
5609 u_z =
fourth * (
uz(i - 1, j - 1, k) +
uz(i, j - 1, k) &
5610 +
uz(i - 1, j, k) +
uz(i, j, k))
5612 v_x =
fourth * (
vx(i - 1, j - 1, k) +
vx(i, j - 1, k) &
5613 +
vx(i - 1, j, k) +
vx(i, j, k))
5614 v_y =
fourth * (
vy(i - 1, j - 1, k) +
vy(i, j - 1, k) &
5615 +
vy(i - 1, j, k) +
vy(i, j, k))
5616 v_z =
fourth * (
vz(i - 1, j - 1, k) +
vz(i, j - 1, k) &
5617 +
vz(i - 1, j, k) +
vz(i, j, k))
5619 w_x =
fourth * (
wx(i - 1, j - 1, k) +
wx(i, j - 1, k) &
5620 +
wx(i - 1, j, k) +
wx(i, j, k))
5621 w_y =
fourth * (
wy(i - 1, j - 1, k) +
wy(i, j - 1, k) &
5622 +
wy(i - 1, j, k) +
wy(i, j, k))
5623 w_z =
fourth * (
wz(i - 1, j - 1, k) +
wz(i, j - 1, k) &
5624 +
wz(i - 1, j, k) +
wz(i, j, k))
5626 q_x =
fourth * (
qx(i - 1, j - 1, k) +
qx(i, j - 1, k) &
5627 +
qx(i - 1, j, k) +
qx(i, j, k))
5628 q_y =
fourth * (
qy(i - 1, j - 1, k) +
qy(i, j - 1, k) &
5629 +
qy(i - 1, j, k) +
qy(i, j, k))
5630 q_z =
fourth * (
qz(i - 1, j - 1, k) +
qz(i, j - 1, k) &
5631 +
qz(i - 1, j, k) +
qz(i, j, k))
5638 ssx =
eighth * (
x(i - 1, j - 1, k + 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
5639 +
x(i - 1, j, k + 1, 1) -
x(i - 1, j, k - 1, 1) &
5640 +
x(i, j - 1, k + 1, 1) -
x(i, j - 1, k - 1, 1) &
5641 +
x(i, j, k + 1, 1) -
x(i, j, k - 1, 1))
5642 ssy =
eighth * (
x(i - 1, j - 1, k + 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
5643 +
x(i - 1, j, k + 1, 2) -
x(i - 1, j, k - 1, 2) &
5644 +
x(i, j - 1, k + 1, 2) -
x(i, j - 1, k - 1, 2) &
5645 +
x(i, j, k + 1, 2) -
x(i, j, k - 1, 2))
5646 ssz =
eighth * (
x(i - 1, j - 1, k + 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
5647 +
x(i - 1, j, k + 1, 3) -
x(i - 1, j, k - 1, 3) &
5648 +
x(i, j - 1, k + 1, 3) -
x(i, j - 1, k - 1, 3) &
5649 +
x(i, j, k + 1, 3) -
x(i, j, k - 1, 3))
5654 snrm =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
5661 corr = u_x * ssx + u_y * ssy + u_z * ssz &
5662 - (
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)) * snrm
5663 u_x = u_x - corr * ssx
5664 u_y = u_y - corr * ssy
5665 u_z = u_z - corr * ssz
5667 corr = v_x * ssx + v_y * ssy + v_z * ssz &
5668 - (
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)) * snrm
5669 v_x = v_x - corr * ssx
5670 v_y = v_y - corr * ssy
5671 v_z = v_z - corr * ssz
5673 corr = w_x * ssx + w_y * ssy + w_z * ssz &
5674 - (
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)) * snrm
5675 w_x = w_x - corr * ssx
5676 w_y = w_y - corr * ssy
5677 w_z = w_z - corr * ssz
5679 corr = q_x * ssx + q_y * ssy + q_z * ssz &
5680 + (
aa(i, j, k + 1) -
aa(i, j, k)) * snrm
5681 q_x = q_x - corr * ssx
5682 q_y = q_y - corr * ssy
5683 q_z = q_z - corr * ssz
5692 fracdiv = twothird * (u_x + v_y + w_z)
5694 tauxxs =
two * u_x - fracdiv
5695 tauyys =
two * v_y - fracdiv
5696 tauzzs =
two * w_z - fracdiv
5702 q_x = heatcoef * q_x
5703 q_y = heatcoef * q_y
5704 q_z = heatcoef * q_z
5724 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
5725 v_x * v_x + v_y * v_y + v_z * v_z + &
5726 w_x * w_x + w_y * w_y + w_z * w_z)
5730 den = max(den, xminn)
5735 fact = mue * ccr1 / den
5747 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
5748 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
5749 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
5751 exy = fact * (wxy * tauyys + wxz * tauyzs + &
5752 wyx * tauxxs + wyz * tauxzs)
5753 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
5754 wzx * tauxxs + wzy * tauxys)
5755 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
5756 wzx * tauxys + wzy * tauyys)
5759 tauxx = mut * tauxxs - exx
5760 tauyy = mut * tauyys - eyy
5761 tauzz = mut * tauzzs - ezz
5762 tauxy = mut * tauxys - exy
5763 tauxz = mut * tauxzs - exz
5764 tauyz = mut * tauyzs - eyz
5769 tauxx = mut * tauxxs
5770 tauyy = mut * tauyys
5771 tauzz = mut * tauzzs
5772 tauxy = mut * tauxys
5773 tauxz = mut * tauxzs
5774 tauyz = mut * tauyzs
5787 fmx = tauxx *
sk(i, j, k, 1) + tauxy *
sk(i, j, k, 2) &
5788 + tauxz *
sk(i, j, k, 3)
5789 fmy = tauxy *
sk(i, j, k, 1) + tauyy *
sk(i, j, k, 2) &
5790 + tauyz *
sk(i, j, k, 3)
5791 fmz = tauxz *
sk(i, j, k, 1) + tauyz *
sk(i, j, k, 2) &
5792 + tauzz *
sk(i, j, k, 3)
5793 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sk(i, j, k, 1)
5794 frhoe = frhoe + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sk(i, j, k, 2)
5795 frhoe = frhoe + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sk(i, j, k, 3)
5796 frhoe = frhoe - q_x *
sk(i, j, k, 1) - q_y *
sk(i, j, k, 2) - q_z *
sk(i, j, k, 3)
5805 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fmx
5806 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fmy
5807 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fmz
5814 tmpstore(1, i, j, 1) = tauxx
5815 tmpstore(2, i, j, 1) = tauyy
5816 tmpstore(3, i, j, 1) = tauzz
5817 tmpstore(4, i, j, 1) = tauxy
5818 tmpstore(5, i, j, 1) = tauxz
5819 tmpstore(6, i, j, 1) = tauyz
5821 tmpstore(7, i, j, 1) = q_x
5822 tmpstore(8, i, j, 1) = q_y
5823 tmpstore(9, i, j, 1) = q_z
5827 tmpstore(1, i, j, 2) = tauxx
5828 tmpstore(2, i, j, 2) = tauyy
5829 tmpstore(3, i, j, 2) = tauzz
5830 tmpstore(4, i, j, 2) = tauxy
5831 tmpstore(5, i, j, 2) = tauxz
5832 tmpstore(6, i, j, 2) = tauyz
5834 tmpstore(7, i, j, 2) = q_x
5835 tmpstore(8, i, j, 2) = q_y
5836 tmpstore(9, i, j, 2) = q_z
5846 origkmin:
if (
kk - 1 == 1)
then
5860 origkmax:
if (
kk +
nz - 1 == bkl)
then
5892 mul = por * (
rlv(i, j, k) +
rlv(i, j + 1, k))
5893 mue = por * (
rev(i, j, k) +
rev(i, j + 1, k))
5900 heatcoef = mul * factlamheat + mue * factturbheat
5905 u_x =
fourth * (
ux(i - 1, j, k - 1) +
ux(i, j, k - 1) &
5906 +
ux(i - 1, j, k) +
ux(i, j, k))
5907 u_y =
fourth * (
uy(i - 1, j, k - 1) +
uy(i, j, k - 1) &
5908 +
uy(i - 1, j, k) +
uy(i, j, k))
5909 u_z =
fourth * (
uz(i - 1, j, k - 1) +
uz(i, j, k - 1) &
5910 +
uz(i - 1, j, k) +
uz(i, j, k))
5912 v_x =
fourth * (
vx(i - 1, j, k - 1) +
vx(i, j, k - 1) &
5913 +
vx(i - 1, j, k) +
vx(i, j, k))
5914 v_y =
fourth * (
vy(i - 1, j, k - 1) +
vy(i, j, k - 1) &
5915 +
vy(i - 1, j, k) +
vy(i, j, k))
5916 v_z =
fourth * (
vz(i - 1, j, k - 1) +
vz(i, j, k - 1) &
5917 +
vz(i - 1, j, k) +
vz(i, j, k))
5919 w_x =
fourth * (
wx(i - 1, j, k - 1) +
wx(i, j, k - 1) &
5920 +
wx(i - 1, j, k) +
wx(i, j, k))
5921 w_y =
fourth * (
wy(i - 1, j, k - 1) +
wy(i, j, k - 1) &
5922 +
wy(i - 1, j, k) +
wy(i, j, k))
5923 w_z =
fourth * (
wz(i - 1, j, k - 1) +
wz(i, j, k - 1) &
5924 +
wz(i - 1, j, k) +
wz(i, j, k))
5926 q_x =
fourth * (
qx(i - 1, j, k - 1) +
qx(i, j, k - 1) &
5927 +
qx(i - 1, j, k) +
qx(i, j, k))
5928 q_y =
fourth * (
qy(i - 1, j, k - 1) +
qy(i, j, k - 1) &
5929 +
qy(i - 1, j, k) +
qy(i, j, k))
5930 q_z =
fourth * (
qz(i - 1, j, k - 1) +
qz(i, j, k - 1) &
5931 +
qz(i - 1, j, k) +
qz(i, j, k))
5938 ssx =
eighth * (
x(i - 1, j + 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
5939 +
x(i - 1, j + 1, k, 1) -
x(i - 1, j - 1, k, 1) &
5940 +
x(i, j + 1, k - 1, 1) -
x(i, j - 1, k - 1, 1) &
5941 +
x(i, j + 1, k, 1) -
x(i, j - 1, k, 1))
5942 ssy =
eighth * (
x(i - 1, j + 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
5943 +
x(i - 1, j + 1, k, 2) -
x(i - 1, j - 1, k, 2) &
5944 +
x(i, j + 1, k - 1, 2) -
x(i, j - 1, k - 1, 2) &
5945 +
x(i, j + 1, k, 2) -
x(i, j - 1, k, 2))
5946 ssz =
eighth * (
x(i - 1, j + 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
5947 +
x(i - 1, j + 1, k, 3) -
x(i - 1, j - 1, k, 3) &
5948 +
x(i, j + 1, k - 1, 3) -
x(i, j - 1, k - 1, 3) &
5949 +
x(i, j + 1, k, 3) -
x(i, j - 1, k, 3))
5954 snrm =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
5961 corr = u_x * ssx + u_y * ssy + u_z * ssz &
5962 - (
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)) * snrm
5963 u_x = u_x - corr * ssx
5964 u_y = u_y - corr * ssy
5965 u_z = u_z - corr * ssz
5967 corr = v_x * ssx + v_y * ssy + v_z * ssz &
5968 - (
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)) * snrm
5969 v_x = v_x - corr * ssx
5970 v_y = v_y - corr * ssy
5971 v_z = v_z - corr * ssz
5973 corr = w_x * ssx + w_y * ssy + w_z * ssz &
5974 - (
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)) * snrm
5975 w_x = w_x - corr * ssx
5976 w_y = w_y - corr * ssy
5977 w_z = w_z - corr * ssz
5979 corr = q_x * ssx + q_y * ssy + q_z * ssz &
5980 + (
aa(i, j + 1, k) -
aa(i, j, k)) * snrm
5981 q_x = q_x - corr * ssx
5982 q_y = q_y - corr * ssy
5983 q_z = q_z - corr * ssz
5992 fracdiv = twothird * (u_x + v_y + w_z)
5994 tauxxs =
two * u_x - fracdiv
5995 tauyys =
two * v_y - fracdiv
5996 tauzzs =
two * w_z - fracdiv
6002 q_x = heatcoef * q_x
6003 q_y = heatcoef * q_y
6004 q_z = heatcoef * q_z
6024 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
6025 v_x * v_x + v_y * v_y + v_z * v_z + &
6026 w_x * w_x + w_y * w_y + w_z * w_z)
6030 den = max(den, xminn)
6035 fact = mue * ccr1 / den
6047 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
6048 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
6049 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
6051 exy = fact * (wxy * tauyys + wxz * tauyzs + &
6052 wyx * tauxxs + wyz * tauxzs)
6053 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
6054 wzx * tauxxs + wzy * tauxys)
6055 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
6056 wzx * tauxys + wzy * tauyys)
6059 tauxx = mut * tauxxs - exx
6060 tauyy = mut * tauyys - eyy
6061 tauzz = mut * tauzzs - ezz
6062 tauxy = mut * tauxys - exy
6063 tauxz = mut * tauxzs - exz
6064 tauyz = mut * tauyzs - eyz
6069 tauxx = mut * tauxxs
6070 tauyy = mut * tauyys
6071 tauzz = mut * tauzzs
6072 tauxy = mut * tauxys
6073 tauxz = mut * tauxzs
6074 tauyz = mut * tauyzs
6087 fmx = tauxx *
sj(i, j, k, 1) + tauxy *
sj(i, j, k, 2) &
6088 + tauxz *
sj(i, j, k, 3)
6089 fmy = tauxy *
sj(i, j, k, 1) + tauyy *
sj(i, j, k, 2) &
6090 + tauyz *
sj(i, j, k, 3)
6091 fmz = tauxz *
sj(i, j, k, 1) + tauyz *
sj(i, j, k, 2) &
6092 + tauzz *
sj(i, j, k, 3)
6093 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sj(i, j, k, 1) &
6094 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sj(i, j, k, 2) &
6095 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sj(i, j, k, 3) &
6096 - q_x *
sj(i, j, k, 1) - q_y *
sj(i, j, k, 2) - q_z *
sj(i, j, k, 3)
6105 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fmx
6106 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fmy
6107 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fmz
6114 tmpstore(1, i, k, 1) = tauxx
6115 tmpstore(2, i, k, 1) = tauyy
6116 tmpstore(3, i, k, 1) = tauzz
6117 tmpstore(4, i, k, 1) = tauxy
6118 tmpstore(5, i, k, 1) = tauxz
6119 tmpstore(6, i, k, 1) = tauyz
6121 tmpstore(7, i, k, 1) = q_x
6122 tmpstore(8, i, k, 1) = q_y
6123 tmpstore(9, i, k, 1) = q_z
6127 tmpstore(1, i, k, 2) = tauxx
6128 tmpstore(2, i, k, 2) = tauyy
6129 tmpstore(3, i, k, 2) = tauzz
6130 tmpstore(4, i, k, 2) = tauxy
6131 tmpstore(5, i, k, 2) = tauxz
6132 tmpstore(6, i, k, 2) = tauyz
6134 tmpstore(7, i, k, 2) = q_x
6135 tmpstore(8, i, k, 2) = q_y
6136 tmpstore(9, i, k, 2) = q_z
6143 origjmin:
if (
jj - 1 == 1)
then
6157 origjmax:
if (
jj +
ny - 1 == bjl)
then
6188 mul = por * (
rlv(i, j, k) +
rlv(i + 1, j, k))
6189 mue = por * (
rev(i, j, k) +
rev(i + 1, j, k))
6196 heatcoef = mul * factlamheat + mue * factturbheat
6201 u_x =
fourth * (
ux(i, j - 1, k - 1) +
ux(i, j, k - 1) &
6202 +
ux(i, j - 1, k) +
ux(i, j, k))
6203 u_y =
fourth * (
uy(i, j - 1, k - 1) +
uy(i, j, k - 1) &
6204 +
uy(i, j - 1, k) +
uy(i, j, k))
6205 u_z =
fourth * (
uz(i, j - 1, k - 1) +
uz(i, j, k - 1) &
6206 +
uz(i, j - 1, k) +
uz(i, j, k))
6208 v_x =
fourth * (
vx(i, j - 1, k - 1) +
vx(i, j, k - 1) &
6209 +
vx(i, j - 1, k) +
vx(i, j, k))
6210 v_y =
fourth * (
vy(i, j - 1, k - 1) +
vy(i, j, k - 1) &
6211 +
vy(i, j - 1, k) +
vy(i, j, k))
6212 v_z =
fourth * (
vz(i, j - 1, k - 1) +
vz(i, j, k - 1) &
6213 +
vz(i, j - 1, k) +
vz(i, j, k))
6215 w_x =
fourth * (
wx(i, j - 1, k - 1) +
wx(i, j, k - 1) &
6216 +
wx(i, j - 1, k) +
wx(i, j, k))
6217 w_y =
fourth * (
wy(i, j - 1, k - 1) +
wy(i, j, k - 1) &
6218 +
wy(i, j - 1, k) +
wy(i, j, k))
6219 w_z =
fourth * (
wz(i, j - 1, k - 1) +
wz(i, j, k - 1) &
6220 +
wz(i, j - 1, k) +
wz(i, j, k))
6222 q_x =
fourth * (
qx(i, j - 1, k - 1) +
qx(i, j, k - 1) &
6223 +
qx(i, j - 1, k) +
qx(i, j, k))
6224 q_y =
fourth * (
qy(i, j - 1, k - 1) +
qy(i, j, k - 1) &
6225 +
qy(i, j - 1, k) +
qy(i, j, k))
6226 q_z =
fourth * (
qz(i, j - 1, k - 1) +
qz(i, j, k - 1) &
6227 +
qz(i, j - 1, k) +
qz(i, j, k))
6234 ssx =
eighth * (
x(i + 1, j - 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6235 +
x(i + 1, j - 1, k, 1) -
x(i - 1, j - 1, k, 1) &
6236 +
x(i + 1, j, k - 1, 1) -
x(i - 1, j, k - 1, 1) &
6237 +
x(i + 1, j, k, 1) -
x(i - 1, j, k, 1))
6238 ssy =
eighth * (
x(i + 1, j - 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6239 +
x(i + 1, j - 1, k, 2) -
x(i - 1, j - 1, k, 2) &
6240 +
x(i + 1, j, k - 1, 2) -
x(i - 1, j, k - 1, 2) &
6241 +
x(i + 1, j, k, 2) -
x(i - 1, j, k, 2))
6242 ssz =
eighth * (
x(i + 1, j - 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6243 +
x(i + 1, j - 1, k, 3) -
x(i - 1, j - 1, k, 3) &
6244 +
x(i + 1, j, k - 1, 3) -
x(i - 1, j, k - 1, 3) &
6245 +
x(i + 1, j, k, 3) -
x(i - 1, j, k, 3))
6250 snrm =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
6257 corr = u_x * ssx + u_y * ssy + u_z * ssz &
6258 - (
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)) * snrm
6259 u_x = u_x - corr * ssx
6260 u_y = u_y - corr * ssy
6261 u_z = u_z - corr * ssz
6263 corr = v_x * ssx + v_y * ssy + v_z * ssz &
6264 - (
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)) * snrm
6265 v_x = v_x - corr * ssx
6266 v_y = v_y - corr * ssy
6267 v_z = v_z - corr * ssz
6269 corr = w_x * ssx + w_y * ssy + w_z * ssz &
6270 - (
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)) * snrm
6271 w_x = w_x - corr * ssx
6272 w_y = w_y - corr * ssy
6273 w_z = w_z - corr * ssz
6275 corr = q_x * ssx + q_y * ssy + q_z * ssz &
6276 + (
aa(i + 1, j, k) -
aa(i, j, k)) * snrm
6277 q_x = q_x - corr * ssx
6278 q_y = q_y - corr * ssy
6279 q_z = q_z - corr * ssz
6288 fracdiv = twothird * (u_x + v_y + w_z)
6290 tauxxs =
two * u_x - fracdiv
6291 tauyys =
two * v_y - fracdiv
6292 tauzzs =
two * w_z - fracdiv
6298 q_x = heatcoef * q_x
6299 q_y = heatcoef * q_y
6300 q_z = heatcoef * q_z
6320 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
6321 v_x * v_x + v_y * v_y + v_z * v_z + &
6322 w_x * w_x + w_y * w_y + w_z * w_z)
6326 den = max(den, xminn)
6331 fact = mue * ccr1 / den
6343 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
6344 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
6345 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
6347 exy = fact * (wxy * tauyys + wxz * tauyzs + &
6348 wyx * tauxxs + wyz * tauxzs)
6349 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
6350 wzx * tauxxs + wzy * tauxys)
6351 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
6352 wzx * tauxys + wzy * tauyys)
6355 tauxx = mut * tauxxs - exx
6356 tauyy = mut * tauyys - eyy
6357 tauzz = mut * tauzzs - ezz
6358 tauxy = mut * tauxys - exy
6359 tauxz = mut * tauxzs - exz
6360 tauyz = mut * tauyzs - eyz
6365 tauxx = mut * tauxxs
6366 tauyy = mut * tauyys
6367 tauzz = mut * tauzzs
6368 tauxy = mut * tauxys
6369 tauxz = mut * tauxzs
6370 tauyz = mut * tauyzs
6383 fmx = tauxx *
si(i, j, k, 1) + tauxy *
si(i, j, k, 2) &
6384 + tauxz *
si(i, j, k, 3)
6385 fmy = tauxy *
si(i, j, k, 1) + tauyy *
si(i, j, k, 2) &
6386 + tauyz *
si(i, j, k, 3)
6387 fmz = tauxz *
si(i, j, k, 1) + tauyz *
si(i, j, k, 2) &
6388 + tauzz *
si(i, j, k, 3)
6389 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
si(i, j, k, 1) &
6390 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
si(i, j, k, 2) &
6391 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
si(i, j, k, 3) &
6392 - q_x *
si(i, j, k, 1) - q_y *
si(i, j, k, 2) - q_z *
si(i, j, k, 3)
6395 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fmx
6396 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fmy
6397 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fmz
6409 tmpstore(1, j, k, 1) = tauxx
6410 tmpstore(2, j, k, 1) = tauyy
6411 tmpstore(3, j, k, 1) = tauzz
6412 tmpstore(4, j, k, 1) = tauxy
6413 tmpstore(5, j, k, 1) = tauxz
6414 tmpstore(6, j, k, 1) = tauyz
6416 tmpstore(7, j, k, 1) = q_x
6417 tmpstore(8, j, k, 1) = q_y
6418 tmpstore(9, j, k, 1) = q_z
6422 tmpstore(1, j, k, 2) = tauxx
6423 tmpstore(2, j, k, 2) = tauyy
6424 tmpstore(3, j, k, 2) = tauzz
6425 tmpstore(4, j, k, 2) = tauxy
6426 tmpstore(5, j, k, 2) = tauxz
6427 tmpstore(6, j, k, 2) = tauyz
6429 tmpstore(7, j, k, 2) = q_x
6430 tmpstore(8, j, k, 2) = q_y
6431 tmpstore(9, j, k, 2) = q_z
6438 origimin:
if (
ii - 1 == 1)
then
6452 origimax:
if (
ii +
nx - 1 == bil)
then
6477 real(kind=realtype),
parameter :: twothird =
two *
third
6481 integer(kind=intType) :: i, j, k
6482 integer(kind=intType) :: ii, jj, kk
6484 real(kind=realtype) :: rfilv, por, mul, mue, mut, heatcoef
6485 real(kind=realtype) :: gm1, factlamheat, factturbheat
6486 real(kind=realtype) :: u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z
6487 real(kind=realtype) :: q_x, q_y, q_z, ubar, vbar, wbar
6488 real(kind=realtype) :: corr, ssx, ssy, ssz,
ss, fracdiv
6489 real(kind=realtype) :: tauxx, tauyy, tauzz
6490 real(kind=realtype) :: tauxy, tauxz, tauyz
6491 real(kind=realtype) :: fmx, fmy, fmz, frhoe
6492 real(kind=realtype) :: dd
6493 logical :: correctForK
6505 ssx =
eighth * (
x(i + 1, j - 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6506 +
x(i + 1, j - 1, k, 1) -
x(i - 1, j - 1, k, 1) &
6507 +
x(i + 1, j, k - 1, 1) -
x(i - 1, j, k - 1, 1) &
6508 +
x(i + 1, j, k, 1) -
x(i - 1, j, k, 1))
6509 ssy =
eighth * (
x(i + 1, j - 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6510 +
x(i + 1, j - 1, k, 2) -
x(i - 1, j - 1, k, 2) &
6511 +
x(i + 1, j, k - 1, 2) -
x(i - 1, j, k - 1, 2) &
6512 +
x(i + 1, j, k, 2) -
x(i - 1, j, k, 2))
6513 ssz =
eighth * (
x(i + 1, j - 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6514 +
x(i + 1, j - 1, k, 3) -
x(i - 1, j - 1, k, 3) &
6515 +
x(i + 1, j, k - 1, 3) -
x(i - 1, j, k - 1, 3) &
6516 +
x(i + 1, j, k, 3) -
x(i - 1, j, k, 3))
6519 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
6525 dd =
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)
6530 dd =
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)
6535 dd =
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)
6540 dd =
aa(i + 1, j, k) -
aa(i, j, k)
6553 mul = por * (
rlv(i, j, k) +
rlv(i + 1, j, k))
6554 mue = por * (
rev(i, j, k) +
rev(i + 1, j, k))
6561 heatcoef = mul * factlamheat + mue * factturbheat
6565 fracdiv = twothird * (u_x + v_y + w_z)
6567 tauxx = mut * (
two * u_x - fracdiv)
6568 tauyy = mut * (
two * v_y - fracdiv)
6569 tauzz = mut * (
two * w_z - fracdiv)
6571 tauxy = mut * (u_y + v_x)
6572 tauxz = mut * (u_z + w_x)
6573 tauyz = mut * (v_z + w_y)
6575 q_x = heatcoef * q_x
6576 q_y = heatcoef * q_y
6577 q_z = heatcoef * q_z
6588 fmx = tauxx *
si(i, j, k, 1) + tauxy *
si(i, j, k, 2) + tauxz *
si(i, j, k, 3)
6589 fmy = tauxy *
si(i, j, k, 1) + tauyy *
si(i, j, k, 2) + tauyz *
si(i, j, k, 3)
6590 fmz = tauxz *
si(i, j, k, 1) + tauyz *
si(i, j, k, 2) + tauzz *
si(i, j, k, 3)
6591 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
si(i, j, k, 1) &
6592 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
si(i, j, k, 2) &
6593 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
si(i, j, k, 3) &
6594 - q_x *
si(i, j, k, 1) - q_y *
si(i, j, k, 2) - q_z *
si(i, j, k, 3)
6597 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fmx
6598 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fmy
6599 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fmz
6618 ssx =
eighth * (
x(i - 1, j + 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6619 +
x(i - 1, j + 1, k, 1) -
x(i - 1, j - 1, k, 1) &
6620 +
x(i, j + 1, k - 1, 1) -
x(i, j - 1, k - 1, 1) &
6621 +
x(i, j + 1, k, 1) -
x(i, j - 1, k, 1))
6622 ssy =
eighth * (
x(i - 1, j + 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6623 +
x(i - 1, j + 1, k, 2) -
x(i - 1, j - 1, k, 2) &
6624 +
x(i, j + 1, k - 1, 2) -
x(i, j - 1, k - 1, 2) &
6625 +
x(i, j + 1, k, 2) -
x(i, j - 1, k, 2))
6626 ssz =
eighth * (
x(i - 1, j + 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6627 +
x(i - 1, j + 1, k, 3) -
x(i - 1, j - 1, k, 3) &
6628 +
x(i, j + 1, k - 1, 3) -
x(i, j - 1, k - 1, 3) &
6629 +
x(i, j + 1, k, 3) -
x(i, j - 1, k, 3))
6632 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
6638 dd =
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)
6643 dd =
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)
6648 dd =
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)
6653 dd =
aa(i, j + 1, k) -
aa(i, j, k)
6666 mul = por * (
rlv(i, j, k) +
rlv(i, j + 1, k))
6667 mue = por * (
rev(i, j, k) +
rev(i, j + 1, k))
6674 heatcoef = mul * factlamheat + mue * factturbheat
6678 fracdiv = twothird * (u_x + v_y + w_z)
6680 tauxx = mut * (
two * u_x - fracdiv)
6681 tauyy = mut * (
two * v_y - fracdiv)
6682 tauzz = mut * (
two * w_z - fracdiv)
6684 tauxy = mut * (u_y + v_x)
6685 tauxz = mut * (u_z + w_x)
6686 tauyz = mut * (v_z + w_y)
6688 q_x = heatcoef * q_x
6689 q_y = heatcoef * q_y
6690 q_z = heatcoef * q_z
6701 fmx = tauxx *
sj(i, j, k, 1) + tauxy *
sj(i, j, k, 2) + tauxz *
sj(i, j, k, 3)
6702 fmy = tauxy *
sj(i, j, k, 1) + tauyy *
sj(i, j, k, 2) + tauyz *
sj(i, j, k, 3)
6703 fmz = tauxz *
sj(i, j, k, 1) + tauyz *
sj(i, j, k, 2) + tauzz *
sj(i, j, k, 3)
6704 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sj(i, j, k, 1) &
6705 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sj(i, j, k, 2) &
6706 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sj(i, j, k, 3) &
6707 - q_x *
sj(i, j, k, 1) - q_y *
sj(i, j, k, 2) - q_z *
sj(i, j, k, 3)
6716 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fmx
6717 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fmy
6718 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fmz
6732 ssx =
eighth * (
x(i - 1, j - 1, k + 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6733 +
x(i - 1, j, k + 1, 1) -
x(i - 1, j, k - 1, 1) &
6734 +
x(i, j - 1, k + 1, 1) -
x(i, j - 1, k - 1, 1) &
6735 +
x(i, j, k + 1, 1) -
x(i, j, k - 1, 1))
6736 ssy =
eighth * (
x(i - 1, j - 1, k + 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6737 +
x(i - 1, j, k + 1, 2) -
x(i - 1, j, k - 1, 2) &
6738 +
x(i, j - 1, k + 1, 2) -
x(i, j - 1, k - 1, 2) &
6739 +
x(i, j, k + 1, 2) -
x(i, j, k - 1, 2))
6740 ssz =
eighth * (
x(i - 1, j - 1, k + 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6741 +
x(i - 1, j, k + 1, 3) -
x(i - 1, j, k - 1, 3) &
6742 +
x(i, j - 1, k + 1, 3) -
x(i, j - 1, k - 1, 3) &
6743 +
x(i, j, k + 1, 3) -
x(i, j, k - 1, 3))
6745 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
6751 dd =
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)
6756 dd =
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)
6761 dd =
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)
6766 dd =
aa(i, j, k + 1) -
aa(i, j, k)
6779 mul = por * (
rlv(i, j, k) +
rlv(i, j, k + 1))
6780 mue = por * (
rev(i, j, k) +
rev(i, j, k + 1))
6787 heatcoef = mul * factlamheat + mue * factturbheat
6791 fracdiv = twothird * (u_x + v_y + w_z)
6793 tauxx = mut * (
two * u_x - fracdiv)
6794 tauyy = mut * (
two * v_y - fracdiv)
6795 tauzz = mut * (
two * w_z - fracdiv)
6797 tauxy = mut * (u_y + v_x)
6798 tauxz = mut * (u_z + w_x)
6799 tauyz = mut * (v_z + w_y)
6801 q_x = heatcoef * q_x
6802 q_y = heatcoef * q_y
6803 q_z = heatcoef * q_z
6814 fmx = tauxx *
sk(i, j, k, 1) + tauxy *
sk(i, j, k, 2) + tauxz *
sk(i, j, k, 3)
6815 fmy = tauxy *
sk(i, j, k, 1) + tauyy *
sk(i, j, k, 2) + tauyz *
sk(i, j, k, 3)
6816 fmz = tauxz *
sk(i, j, k, 1) + tauyz *
sk(i, j, k, 2) + tauzz *
sk(i, j, k, 3)
6817 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sk(i, j, k, 1) &
6818 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sk(i, j, k, 2) &
6819 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sk(i, j, k, 3) &
6820 - q_x *
sk(i, j, k, 1) - q_y *
sk(i, j, k, 2) - q_z *
sk(i, j, k, 3)
6829 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fmx
6830 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fmy
6831 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fmz
6850 integer(kind=intType) :: nTurb, i, j, k, l
6851 real(kind=realtype) :: ovol, rblank
6858 rblank = max(real(
iblank(i, j, k), realtype),
zero)
6859 dw(i, j, k, l) = (
dw(i, j, k, l) +
fw(i, j, k, l)) * rblank
6874 integer(kind=intType) :: i, j, k, ii, nTurb
6875 real(kind=realtype) :: ovol
6884 dw(i, j, k, 1:
nwf) =
dw(i, j, k, 1:
nwf) * ovol
subroutine riemannflux(left, right, flux)
subroutine leftrightstate(du1, du2, du3, rotmatrix, left, right)
subroutine getforces(forces, npts, sps)
integer(kind=inttype) nactuatorregions
subroutine setbcdatafinegrid(initializationPart)
subroutine setbcdata(bcDataNamesIn, bcDataIn, famLists, sps, nVar, nFamMax)
subroutine applyallbc_block(secondHalo)
integer(kind=inttype) ndom
integer(kind=inttype), parameter bbkb
subroutine viscousflux(storeWallTensor)
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) vol
integer(kind=inttype), dimension(2:bbil, 2:bbjl, 2:bbkl) iblank
integer(kind=inttype), parameter bbkl
integer(kind=inttype), parameter bbjb
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke, 1:6) dw
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) qz
subroutine timestep(updateDtl)
subroutine invisciddissfluxscalar
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) wy
subroutine blocketterescore(dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall)
integer(kind=inttype), parameter bbie
integer(kind=inttype), parameter bs
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) rev
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) vy
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) wz
subroutine initres(varStart, varEnd)
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) vz
subroutine invisciddissfluxmatrixapprox
real(kind=realtype), dimension(0:bbib, 0:bbjb, 0:bbkb) gamma
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) qx
subroutine invisciddissfluxscalarapprox
integer(kind=inttype), parameter bbje
real(kind=realtype), dimension(0:bbie, 0:bbje, 0:bbke, 3) x
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) uz
integer(kind=inttype), parameter bbib
subroutine viscousfluxapprox
integer(kind=inttype) doublehalostart
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke, 1:5) fw
subroutine inviscidcentralflux
real(kind=realtype), dimension(1:bbie, 0:bbje, 1:bbke) sfacej
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) aa
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) radi
subroutine inviscidupwindflux(fineGrid)
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) radj
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) rlv
real(kind=realtype), dimension(0:bbie, 1:bbje, 1:bbke, 3) si
integer(kind=inttype), parameter bbke
real(kind=realtype), dimension(1:bbie, 0:bbje, 1:bbke, 3) sj
subroutine invisciddissfluxmatrix
real(kind=realtype), dimension(2:bbil, 2:bbjl, 2:bbkl) volref
integer(kind=inttype) nodestart
subroutine computespeedofsoundsquared
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) dtl
integer(kind=inttype), parameter bbil
subroutine blockrescore(dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall, nn, sps)
real(kind=realtype), dimension(0:bbib, 0:bbjb, 0:bbkb) p
real(kind=realtype), dimension(0:bbib, 0:bbjb, 0:bbkb, 1:6) w
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke) radk
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) wx
real(kind=realtype), dimension(2:bbil, 2:bbjl, 2:bbkl) d2wall
real(kind=realtype), dimension(1:bbie, 1:bbje, 1:bbke, 3) dss
real(kind=realtype), dimension(1:bbie, 1:bbje, 0:bbke) sfacek
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) uy
subroutine allnodalgradients
integer(kind=inttype) singlehalostart
real(kind=realtype), dimension(0:bbie, 1:bbje, 1:bbke) sfacei
integer(kind=portype), dimension(2:bbil, 2:bbjl, 1:bbkl) pork
real(kind=realtype), dimension(0:bbib, 0:bbjb, 0:bbkb) ss
integer(kind=portype), dimension(1:bbil, 2:bbjl, 2:bbkl) pori
subroutine blocketteres(useDissApprox, useViscApprox, useUpdateIntermed, useFlowRes, useTurbRes, useSpatial, useStoreWall, famLists, funcValues, forces, bcDataNames, bcDataValues, bcDataFamLists)
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) ux
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) vx
integer(kind=portype), dimension(2:bbil, 1:bbjl, 2:bbkl) porj
real(kind=realtype), dimension(1:bbie, 1:bbje, 0:bbke, 3) sk
integer(kind=inttype), parameter bbjl
real(kind=realtype), dimension(1:bbil, 1:bbjl, 1:bbkl) qy
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 viscjmaxpointer
real(kind=realtype), dimension(:, :, :), pointer qy
real(kind=realtype), dimension(:, :, :), pointer aa
real(kind=realtype), dimension(:, :, :), pointer uz
real(kind=realtype), dimension(:, :, :), pointer p
real(kind=realtype), dimension(:, :, :), pointer radj
real(kind=realtype), dimension(:, :, :, :), pointer w
real(kind=realtype), dimension(:, :, :), pointer uy
real(kind=realtype), dimension(:, :, :), pointer sfacei
type(viscsubfacetype), dimension(:), pointer viscsubface
integer(kind=portype), dimension(:, :, :), pointer porj
real(kind=realtype), dimension(:, :, :), pointer d2wall
integer(kind=portype), dimension(:, :, :), pointer pori
integer(kind=inttype), dimension(:, :, :), pointer iblank
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 volref
real(kind=realtype), dimension(:, :, :, :), pointer sj
integer(kind=inttype), dimension(:, :), pointer visckminpointer
integer(kind=inttype) sectionid
real(kind=realtype), dimension(:, :, :), pointer qx
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 dtl
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
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
integer(kind=inttype), parameter strain
integer(kind=inttype), parameter spalartallmaras
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 pi
real(kind=realtype), parameter eps
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=inttype), parameter upwind
integer(kind=portype), parameter boundflux
integer(kind=portype), parameter noflux
integer(kind=inttype), parameter eulerequations
integer(kind=inttype), parameter timespectral
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 dissmatrix
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 updatefast
integer(kind=inttype), parameter nsequations
integer(kind=inttype), parameter secondorder
integer(kind=inttype), parameter dissscalar
integer(kind=inttype), parameter ransequations
real(kind=realtype), parameter sixth
subroutine computepressuresimple(includeHalos)
subroutine computelamviscosity(includeHalos)
subroutine adjustinflowangle()
subroutine etot(rho, u, v, w, p, k, etotal, correctForK)
real(kind=realtype) gammainf
integer(kind=inttype) nt1
real(kind=realtype) pinfcorr
integer(kind=inttype) nwf
real(kind=realtype) rhoinf
real(kind=realtype) timeref
integer(kind=inttype) nt2
subroutine whalo2(level, start, end, commPressure, commGamma, commViscous)
subroutine exchangecoor(level)
subroutine referencestate
real(kind=realtype) totalr0
integer(kind=inttype) currentlevel
real(kind=realtype) totalr
integer(kind=inttype) groundlevel
subroutine updateoversetconnectivity(level, sps)
real(kind=realtype) rsacw1
subroutine sourceterms_block(nn, res, iRegion, pLocal)
subroutine initres_block(varStart, varEnd, nn, sps)
real(kind=realtype) kar2inv
real(kind=realtype) cb3inv
subroutine sa_block(resOnly)
integer(kind=inttype) nsections
type(sectiontype), dimension(:), allocatable sections
subroutine timestep_block(onlyRadii)
subroutine getsolution(famLists, funcValues, globalValues)
subroutine bcturbtreatment
subroutine applyallturbbcthisblock(secondHalo)
subroutine computeeddyviscosity(includeHalos)
subroutine echk(errorcode, file, line)
logical function getcorrectfork()
real(kind=realtype) function mydim(x, y)
subroutine setpointers(nn, mm, ll)
subroutine terminate(routineName, errorMessage)
subroutine updatewalldistancesquickly(nn, level, sps)