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
1031 uux =
w(i + 1, j, k,
ivx) *
si(i, j, k, 1) -
w(i - 1, j, k,
ivx) *
si(i - 1, j, k, 1) &
1032 +
w(i, j + 1, k,
ivx) *
sj(i, j, k, 1) -
w(i, j - 1, k,
ivx) *
sj(i, j - 1, k, 1) &
1033 +
w(i, j, k + 1,
ivx) *
sk(i, j, k, 1) -
w(i, j, k - 1,
ivx) *
sk(i, j, k - 1, 1)
1034 uuy =
w(i + 1, j, k,
ivx) *
si(i, j, k, 2) -
w(i - 1, j, k,
ivx) *
si(i - 1, j, k, 2) &
1035 +
w(i, j + 1, k,
ivx) *
sj(i, j, k, 2) -
w(i, j - 1, k,
ivx) *
sj(i, j - 1, k, 2) &
1036 +
w(i, j, k + 1,
ivx) *
sk(i, j, k, 2) -
w(i, j, k - 1,
ivx) *
sk(i, j, k - 1, 2)
1037 uuz =
w(i + 1, j, k,
ivx) *
si(i, j, k, 3) -
w(i - 1, j, k,
ivx) *
si(i - 1, j, k, 3) &
1038 +
w(i, j + 1, k,
ivx) *
sj(i, j, k, 3) -
w(i, j - 1, k,
ivx) *
sj(i, j - 1, k, 3) &
1039 +
w(i, j, k + 1,
ivx) *
sk(i, j, k, 3) -
w(i, j, k - 1,
ivx) *
sk(i, j, k - 1, 3)
1043 vvx =
w(i + 1, j, k,
ivy) *
si(i, j, k, 1) -
w(i - 1, j, k,
ivy) *
si(i - 1, j, k, 1) &
1044 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 1) -
w(i, j - 1, k,
ivy) *
sj(i, j - 1, k, 1) &
1045 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 1) -
w(i, j, k - 1,
ivy) *
sk(i, j, k - 1, 1)
1046 vvy =
w(i + 1, j, k,
ivy) *
si(i, j, k, 2) -
w(i - 1, j, k,
ivy) *
si(i - 1, j, k, 2) &
1047 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 2) -
w(i, j - 1, k,
ivy) *
sj(i, j - 1, k, 2) &
1048 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 2) -
w(i, j, k - 1,
ivy) *
sk(i, j, k - 1, 2)
1049 vvz =
w(i + 1, j, k,
ivy) *
si(i, j, k, 3) -
w(i - 1, j, k,
ivy) *
si(i - 1, j, k, 3) &
1050 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 3) -
w(i, j - 1, k,
ivy) *
sj(i, j - 1, k, 3) &
1051 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 3) -
w(i, j, k - 1,
ivy) *
sk(i, j, k - 1, 3)
1055 wwx =
w(i + 1, j, k,
ivz) *
si(i, j, k, 1) -
w(i - 1, j, k,
ivz) *
si(i - 1, j, k, 1) &
1056 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 1) -
w(i, j - 1, k,
ivz) *
sj(i, j - 1, k, 1) &
1057 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 1) -
w(i, j, k - 1,
ivz) *
sk(i, j, k - 1, 1)
1058 wwy =
w(i + 1, j, k,
ivz) *
si(i, j, k, 2) -
w(i - 1, j, k,
ivz) *
si(i - 1, j, k, 2) &
1059 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 2) -
w(i, j - 1, k,
ivz) *
sj(i, j - 1, k, 2) &
1060 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 2) -
w(i, j, k - 1,
ivz) *
sk(i, j, k - 1, 2)
1061 wwz =
w(i + 1, j, k,
ivz) *
si(i, j, k, 3) -
w(i - 1, j, k,
ivz) *
si(i - 1, j, k, 3) &
1062 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 3) -
w(i, j - 1, k,
ivz) *
sj(i, j - 1, k, 3) &
1063 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 3) -
w(i, j, k - 1,
ivz) *
sk(i, j, k - 1, 3)
1073 sxx =
two * fact * uux
1074 syy =
two * fact * vvy
1075 szz =
two * fact * wwz
1077 sxy = fact * (uuy + vvx)
1078 sxz = fact * (uuz + wwx)
1079 syz = fact * (vvz + wwy)
1083 div2 = f23 * (sxx + syy + szz)**2
1087 strainmag2 =
two * (sxy**2 + sxz**2 + syz**2) &
1088 + sxx**2 + syy**2 + szz**2
1095 vortx =
two * fact * (wwy - vvz) -
two * omegax
1096 vorty =
two * fact * (uuz - wwx) -
two * omegay
1097 vortz =
two * fact * (vvx - uuy) -
two * omegaz
1100 sqrtprod = sqrt(max(
two * strainmag2 - div2,
eps))
1102 sqrtprod = sqrt(vortx**2 + vorty**2 + vortz**2)
1110 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1112 chi =
w(i, j, k,
itu1) / nu
1115 fv1 = chi3 / (chi3 +
cv13)
1116 fv2 =
one - chi / (
one + chi * fv1)
1149 rr = min(rr, 10.0_realtype)
1150 gg = rr +
rsacw2 * (rr**6 - rr)
1158 term1 =
rsacb1 * (
one - ft2) * sqrtprod * term1fact
1162 dw(i, j, k,
itu1) =
dw(i, j, k,
itu1) + (term1 + term2 *
w(i, j, k,
itu1)) *
w(i, j, k,
itu1)
1180 real(kind=realtype) :: voli, volmi, volpi, xm, ym, zm, xp, yp, zp
1181 real(kind=realtype) :: xa, ya, za, ttm, ttp, cnud, cam, cap
1182 real(kind=realtype) :: nutm, nutp, num, nup, cdm, cdp
1183 real(kind=realtype) :: c1m, c1p, c10, b1, c1, d1, qs, nu
1184 integer(Kind=intType) :: i, j, k
1202 voli =
one /
vol(i, j, k)
1203 volmi =
two / (
vol(i, j, k) +
vol(i, j, k - 1))
1204 volpi =
two / (
vol(i, j, k) +
vol(i, j, k + 1))
1206 xm =
sk(i, j, k - 1, 1) * volmi
1207 ym =
sk(i, j, k - 1, 2) * volmi
1208 zm =
sk(i, j, k - 1, 3) * volmi
1209 xp =
sk(i, j, k, 1) * volpi
1210 yp =
sk(i, j, k, 2) * volpi
1211 zp =
sk(i, j, k, 3) * volpi
1213 xa =
half * (
sk(i, j, k, 1) +
sk(i, j, k - 1, 1)) * voli
1214 ya =
half * (
sk(i, j, k, 2) +
sk(i, j, k - 1, 2)) * voli
1215 za =
half * (
sk(i, j, k, 3) +
sk(i, j, k - 1, 3)) * voli
1216 ttm = xm * xa + ym * ya + zm * za
1217 ttp = xp * xa + yp * ya + zp * za
1239 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1240 num =
half * (
rlv(i, j, k - 1) /
w(i, j, k - 1,
irho) + nu)
1241 nup =
half * (
rlv(i, j, k + 1) /
w(i, j, k + 1,
irho) + nu)
1245 c1m = max(cdm + cam,
zero)
1246 c1p = max(cdp + cap,
zero)
1253 - c10 *
w(i, j, k,
itu1) + c1p *
w(i, j, k + 1,
itu1)
1267 voli =
one /
vol(i, j, k)
1268 volmi =
two / (
vol(i, j, k) +
vol(i, j - 1, k))
1269 volpi =
two / (
vol(i, j, k) +
vol(i, j + 1, k))
1271 xm =
sj(i, j - 1, k, 1) * volmi
1272 ym =
sj(i, j - 1, k, 2) * volmi
1273 zm =
sj(i, j - 1, k, 3) * volmi
1274 xp =
sj(i, j, k, 1) * volpi
1275 yp =
sj(i, j, k, 2) * volpi
1276 zp =
sj(i, j, k, 3) * volpi
1278 xa =
half * (
sj(i, j, k, 1) +
sj(i, j - 1, k, 1)) * voli
1279 ya =
half * (
sj(i, j, k, 2) +
sj(i, j - 1, k, 2)) * voli
1280 za =
half * (
sj(i, j, k, 3) +
sj(i, j - 1, k, 3)) * voli
1281 ttm = xm * xa + ym * ya + zm * za
1282 ttp = xp * xa + yp * ya + zp * za
1304 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1305 num =
half * (
rlv(i, j - 1, k) /
w(i, j - 1, k,
irho) + nu)
1306 nup =
half * (
rlv(i, j + 1, k) /
w(i, j + 1, k,
irho) + nu)
1310 c1m = max(cdm + cam,
zero)
1311 c1p = max(cdp + cap,
zero)
1318 - c10 *
w(i, j, k,
itu1) + c1p *
w(i, j + 1, k,
itu1)
1333 voli =
one /
vol(i, j, k)
1334 volmi =
two / (
vol(i, j, k) +
vol(i - 1, j, k))
1335 volpi =
two / (
vol(i, j, k) +
vol(i + 1, j, k))
1337 xm =
si(i - 1, j, k, 1) * volmi
1338 ym =
si(i - 1, j, k, 2) * volmi
1339 zm =
si(i - 1, j, k, 3) * volmi
1340 xp =
si(i, j, k, 1) * volpi
1341 yp =
si(i, j, k, 2) * volpi
1342 zp =
si(i, j, k, 3) * volpi
1344 xa =
half * (
si(i, j, k, 1) +
si(i - 1, j, k, 1)) * voli
1345 ya =
half * (
si(i, j, k, 2) +
si(i - 1, j, k, 2)) * voli
1346 za =
half * (
si(i, j, k, 3) +
si(i - 1, j, k, 3)) * voli
1347 ttm = xm * xa + ym * ya + zm * za
1348 ttp = xp * xa + yp * ya + zp * za
1370 nu =
rlv(i, j, k) /
w(i, j, k,
irho)
1371 num =
half * (
rlv(i - 1, j, k) /
w(i - 1, j, k,
irho) + nu)
1372 nup =
half * (
rlv(i + 1, j, k) /
w(i + 1, j, k,
irho) + nu)
1376 c1m = max(cdm + cam,
zero)
1377 c1p = max(cdp + cap,
zero)
1384 - c10 *
w(i, j, k,
itu1) + c1p *
w(i + 1, j, k,
itu1)
1401 real(kind=realtype) :: uu, dwt, dwtm1, dwtp1, dwti, dwtj, dwtk, qs
1402 real(kind=realtype) :: voli, xa, ya, za
1403 integer(kind=intType),
parameter :: nAdv = 1
1404 integer(kind=intType) :: offset, i, j, k, ii, jj
1426 xa = (
sk(i, j, k, 1) +
sk(i, j, k - 1, 1)) * voli
1427 ya = (
sk(i, j, k, 2) +
sk(i, j, k - 1, 2)) * voli
1428 za = (
sk(i, j, k, 3) +
sk(i, j, k - 1, 3)) * voli
1430 uu = xa *
w(i, j, k,
ivx) + ya *
w(i, j, k,
ivy) + za *
w(i, j, k,
ivz) - qs
1435 velkdir:
if (uu >
zero)
then
1455 dwtm1 =
w(i, j, k - 1, jj) -
w(i, j, k - 2, jj)
1456 dwt =
w(i, j, k, jj) -
w(i, j, k - 1, jj)
1457 dwtp1 =
w(i, j, k + 1, jj) -
w(i, j, k, jj)
1465 if (dwt * dwtp1 >
zero)
then
1466 if (abs(dwt) < abs(dwtp1))
then
1467 dwtk = dwtk +
half * dwt
1469 dwtk = dwtk +
half * dwtp1
1473 if (dwt * dwtm1 >
zero)
then
1474 if (abs(dwt) < abs(dwtm1))
then
1475 dwtk = dwtk -
half * dwt
1477 dwtk = dwtk -
half * dwtm1
1485 dwtk =
w(i, j, k, jj) -
w(i, j, k - 1, jj)
1493 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtk
1515 dwtm1 =
w(i, j, k, jj) -
w(i, j, k - 1, jj)
1516 dwt =
w(i, j, k + 1, jj) -
w(i, j, k, jj)
1517 dwtp1 =
w(i, j, k + 2, jj) -
w(i, j, k + 1, jj)
1525 if (dwt * dwtp1 >
zero)
then
1526 if (abs(dwt) < abs(dwtp1))
then
1527 dwtk = dwtk -
half * dwt
1529 dwtk = dwtk -
half * dwtp1
1533 if (dwt * dwtm1 >
zero)
then
1534 if (abs(dwt) < abs(dwtm1))
then
1535 dwtk = dwtk +
half * dwt
1537 dwtk = dwtk +
half * dwtm1
1545 dwtk =
w(i, j, k + 1, jj) -
w(i, j, k, jj)
1553 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtk
1580 xa = (
sj(i, j, k, 1) +
sj(i, j - 1, k, 1)) * voli
1581 ya = (
sj(i, j, k, 2) +
sj(i, j - 1, k, 2)) * voli
1582 za = (
sj(i, j, k, 3) +
sj(i, j - 1, k, 3)) * voli
1584 uu = xa *
w(i, j, k,
ivx) + ya *
w(i, j, k,
ivy) + za *
w(i, j, k,
ivz) - qs
1589 veljdir:
if (uu >
zero)
then
1608 dwtm1 =
w(i, j - 1, k, jj) -
w(i, j - 2, k, jj)
1609 dwt =
w(i, j, k, jj) -
w(i, j - 1, k, jj)
1610 dwtp1 =
w(i, j + 1, k, jj) -
w(i, j, k, jj)
1618 if (dwt * dwtp1 >
zero)
then
1619 if (abs(dwt) < abs(dwtp1))
then
1620 dwtj = dwtj +
half * dwt
1622 dwtj = dwtj +
half * dwtp1
1626 if (dwt * dwtm1 >
zero)
then
1627 if (abs(dwt) < abs(dwtm1))
then
1628 dwtj = dwtj -
half * dwt
1630 dwtj = dwtj -
half * dwtm1
1638 dwtj =
w(i, j, k, jj) -
w(i, j - 1, k, jj)
1646 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtj
1668 dwtm1 =
w(i, j, k, jj) -
w(i, j - 1, k, jj)
1669 dwt =
w(i, j + 1, k, jj) -
w(i, j, k, jj)
1670 dwtp1 =
w(i, j + 2, k, jj) -
w(i, j + 1, k, jj)
1678 if (dwt * dwtp1 >
zero)
then
1679 if (abs(dwt) < abs(dwtp1))
then
1680 dwtj = dwtj -
half * dwt
1682 dwtj = dwtj -
half * dwtp1
1686 if (dwt * dwtm1 >
zero)
then
1687 if (abs(dwt) < abs(dwtm1))
then
1688 dwtj = dwtj +
half * dwt
1690 dwtj = dwtj +
half * dwtm1
1698 dwtj =
w(i, j + 1, k, jj) -
w(i, j, k, jj)
1706 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwtj
1732 xa = (
si(i, j, k, 1) +
si(i - 1, j, k, 1)) * voli
1733 ya = (
si(i, j, k, 2) +
si(i - 1, j, k, 2)) * voli
1734 za = (
si(i, j, k, 3) +
si(i - 1, j, k, 3)) * voli
1736 uu = xa *
w(i, j, k,
ivx) + ya *
w(i, j, k,
ivy) + za *
w(i, j, k,
ivz) - qs
1741 velidir:
if (uu >
zero)
then
1760 dwtm1 =
w(i - 1, j, k, jj) -
w(i - 2, j, k, jj)
1761 dwt =
w(i, j, k, jj) -
w(i - 1, j, k, jj)
1762 dwtp1 =
w(i + 1, j, k, jj) -
w(i, j, k, jj)
1770 if (dwt * dwtp1 >
zero)
then
1771 if (abs(dwt) < abs(dwtp1))
then
1772 dwti = dwti +
half * dwt
1774 dwti = dwti +
half * dwtp1
1778 if (dwt * dwtm1 >
zero)
then
1779 if (abs(dwt) < abs(dwtm1))
then
1780 dwti = dwti -
half * dwt
1782 dwti = dwti -
half * dwtm1
1790 dwti =
w(i, j, k, jj) -
w(i - 1, j, k, jj)
1798 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwti
1820 dwtm1 =
w(i, j, k, jj) -
w(i - 1, j, k, jj)
1821 dwt =
w(i + 1, j, k, jj) -
w(i, j, k, jj)
1822 dwtp1 =
w(i + 2, j, k, jj) -
w(i + 1, j, k, jj)
1830 if (dwt * dwtp1 >
zero)
then
1831 if (abs(dwt) < abs(dwtp1))
then
1832 dwti = dwti -
half * dwt
1834 dwti = dwti -
half * dwtp1
1838 if (dwt * dwtm1 >
zero)
then
1839 if (abs(dwt) < abs(dwtm1))
then
1840 dwti = dwti +
half * dwt
1842 dwti = dwti +
half * dwtm1
1850 dwti =
w(i + 1, j, k, jj) -
w(i, j, k, jj)
1858 dw(i, j, k,
itu1 + ii - 1) =
dw(i, j, k,
itu1 + ii - 1) - uu * dwti
1883 integer(kind=intType) :: i, j, k, ii
1884 real(kind=realtype) :: rblank
1889 rblank = max(real(
iblank(i, j, k), realtype),
zero)
1913 logical,
intent(in),
optional :: updateDtl
1916 real(kind=realtype),
parameter :: b = 2.0_realtype
1919 real(kind=realtype) :: plim, rlim, clim2
1920 real(kind=realtype) :: cc2, qsi, qsj, qsk, sx, sy, sz, rmu
1921 real(kind=realtype) :: ri, rj, rk, rij, rjk, rki
1922 real(kind=realtype) :: vsi, vsj, vsk, rfl, dpi, dpj, dpk
1923 real(kind=realtype) :: sface, tmp, uux, uuy, uuz
1924 logical :: doScaling, updateDt
1925 integer(kind=intType) :: i, j, k
1928 if (
present(updatedtl))
then
1937 rlim = 0.001_realtype *
rhoinf
1959 uux =
w(i, j, k,
ivx)
1960 uuy =
w(i, j, k,
ivy)
1961 uuz =
w(i, j, k,
ivz)
1962 cc2 =
gamma(i, j, k) *
p(i, j, k) /
w(i, j, k,
irho)
1963 cc2 = max(cc2, clim2)
1974 sx =
si(i - 1, j, k, 1) +
si(i, j, k, 1)
1975 sy =
si(i - 1, j, k, 2) +
si(i, j, k, 2)
1976 sz =
si(i - 1, j, k, 3) +
si(i, j, k, 3)
1978 qsi = uux * sx + uuy * sy + uuz * sz - sface
1980 ri =
half * (abs(qsi) &
1988 sx =
sj(i, j - 1, k, 1) +
sj(i, j, k, 1)
1989 sy =
sj(i, j - 1, k, 2) +
sj(i, j, k, 2)
1990 sz =
sj(i, j - 1, k, 3) +
sj(i, j, k, 3)
1992 qsj = uux * sx + uuy * sy + uuz * sz - sface
1994 rj =
half * (abs(qsj) &
2002 sx =
sk(i, j, k - 1, 1) +
sk(i, j, k, 1)
2003 sy =
sk(i, j, k - 1, 2) +
sk(i, j, k, 2)
2004 sz =
sk(i, j, k - 1, 3) +
sk(i, j, k, 3)
2006 qsk = uux * sx + uuy * sy + uuz * sz - sface
2008 rk =
half * (abs(qsk) &
2013 dtl(i, j, k) = ri + rj + rk
2026 rij = (ri / rj)**
adis
2027 rjk = (rj / rk)**
adis
2028 rki = (rk / ri)**
adis
2045 viscousterm:
if (
viscous)
then
2066 rmu = rmu +
rev(i, j, k)
2072 sx =
si(i, j, k, 1) +
si(i - 1, j, k, 1)
2073 sy =
si(i, j, k, 2) +
si(i - 1, j, k, 2)
2074 sz =
si(i, j, k, 3) +
si(i - 1, j, k, 3)
2076 vsi = rmu * (sx * sx + sy * sy + sz * sz)
2077 dtl(i, j, k) =
dtl(i, j, k) + vsi
2082 sx =
sj(i, j, k, 1) +
sj(i, j - 1, k, 1)
2083 sy =
sj(i, j, k, 2) +
sj(i, j - 1, k, 2)
2084 sz =
sj(i, j, k, 3) +
sj(i, j - 1, k, 3)
2086 vsj = rmu * (sx * sx + sy * sy + sz * sz)
2087 dtl(i, j, k) =
dtl(i, j, k) + vsj
2092 sx =
sk(i, j, k, 1) +
sk(i, j, k - 1, 1)
2093 sy =
sk(i, j, k, 2) +
sk(i, j, k - 1, 2)
2094 sz =
sk(i, j, k, 3) +
sk(i, j, k - 1, 3)
2096 vsk = rmu * (sx * sx + sy * sy + sz * sz)
2097 dtl(i, j, k) =
dtl(i, j, k) + vsk
2118 dtl(i, j, k) =
dtl(i, j, k) + tmp *
vol(i, j, k)
2132 dpi = abs(
p(i + 1, j, k) -
two *
p(i, j, k) +
p(i - 1, j, k)) &
2133 / (
p(i + 1, j, k) +
two *
p(i, j, k) +
p(i - 1, j, k) + plim)
2134 dpj = abs(
p(i, j + 1, k) -
two *
p(i, j, k) +
p(i, j - 1, k)) &
2135 / (
p(i, j + 1, k) +
two *
p(i, j, k) +
p(i, j - 1, k) + plim)
2136 dpk = abs(
p(i, j, k + 1) -
two *
p(i, j, k) +
p(i, j, k - 1)) &
2137 / (
p(i, j, k + 1) +
two *
p(i, j, k) +
p(i, j, k - 1) + plim)
2138 rfl =
one / (
one + b * (dpi + dpj + dpk))
2140 dtl(i, j, k) = rfl /
dtl(i, j, k)
2161 real(kind=realtype) :: qsp, qsm, rqsp, rqsm, porvel, porflux
2162 real(kind=realtype) :: pa, vnp, vnm, fs, sface
2163 integer(kind=intType) :: i, j, k
2164 real(kind=realtype) :: wwx, wwy, wwz, rvol
2177 vnp =
w(i + 1, j, k,
ivx) *
si(i, j, k, 1) &
2178 +
w(i + 1, j, k,
ivy) *
si(i, j, k, 2) &
2179 +
w(i + 1, j, k,
ivz) *
si(i, j, k, 3)
2180 vnm =
w(i, j, k,
ivx) *
si(i, j, k, 1) &
2181 +
w(i, j, k,
ivy) *
si(i, j, k, 2) &
2182 +
w(i, j, k,
ivz) *
si(i, j, k, 3)
2204 porvel = porvel * porflux
2209 qsp = (vnp - sface) * porvel
2210 qsm = (vnm - sface) * porvel
2212 rqsp = qsp *
w(i + 1, j, k,
irho)
2213 rqsm = qsm *
w(i, j, k,
irho)
2219 pa = porflux * (
p(i + 1, j, k) +
p(i, j, k))
2229 fs = rqsp *
w(i + 1, j, k,
ivx) + rqsm *
w(i, j, k,
ivx) &
2230 + pa *
si(i, j, k, 1)
2231 dw(i + 1, j, k,
imx) =
dw(i + 1, j, k,
imx) - fs
2234 fs = rqsp *
w(i + 1, j, k,
ivy) + rqsm *
w(i, j, k,
ivy) &
2235 + pa *
si(i, j, k, 2)
2236 dw(i + 1, j, k,
imy) =
dw(i + 1, j, k,
imy) - fs
2239 fs = rqsp *
w(i + 1, j, k,
ivz) + rqsm *
w(i, j, k,
ivz) &
2240 + pa *
si(i, j, k, 3)
2241 dw(i + 1, j, k,
imz) =
dw(i + 1, j, k,
imz) - fs
2244 fs = qsp *
w(i + 1, j, k,
irhoe) + qsm *
w(i, j, k,
irhoe) &
2245 + porflux * (vnp *
p(i + 1, j, k) + vnm *
p(i, j, k))
2263 vnp =
w(i, j + 1, k,
ivx) *
sj(i, j, k, 1) &
2264 +
w(i, j + 1, k,
ivy) *
sj(i, j, k, 2) &
2265 +
w(i, j + 1, k,
ivz) *
sj(i, j, k, 3)
2266 vnm =
w(i, j, k,
ivx) *
sj(i, j, k, 1) &
2267 +
w(i, j, k,
ivy) *
sj(i, j, k, 2) &
2268 +
w(i, j, k,
ivz) *
sj(i, j, k, 3)
2291 porvel = porvel * porflux
2296 qsp = (vnp - sface) * porvel
2297 qsm = (vnm - sface) * porvel
2299 rqsp = qsp *
w(i, j + 1, k,
irho)
2300 rqsm = qsm *
w(i, j, k,
irho)
2306 pa = porflux * (
p(i, j + 1, k) +
p(i, j, k))
2316 fs = rqsp *
w(i, j + 1, k,
ivx) + rqsm *
w(i, j, k,
ivx) &
2317 + pa *
sj(i, j, k, 1)
2318 dw(i, j + 1, k,
imx) =
dw(i, j + 1, k,
imx) - fs
2321 fs = rqsp *
w(i, j + 1, k,
ivy) + rqsm *
w(i, j, k,
ivy) &
2322 + pa *
sj(i, j, k, 2)
2323 dw(i, j + 1, k,
imy) =
dw(i, j + 1, k,
imy) - fs
2326 fs = rqsp *
w(i, j + 1, k,
ivz) + rqsm *
w(i, j, k,
ivz) &
2327 + pa *
sj(i, j, k, 3)
2328 dw(i, j + 1, k,
imz) =
dw(i, j + 1, k,
imz) - fs
2331 fs = qsp *
w(i, j + 1, k,
irhoe) + qsm *
w(i, j, k,
irhoe) &
2332 + porflux * (vnp *
p(i, j + 1, k) + vnm *
p(i, j, k))
2350 vnp =
w(i, j, k + 1,
ivx) *
sk(i, j, k, 1) &
2351 +
w(i, j, k + 1,
ivy) *
sk(i, j, k, 2) &
2352 +
w(i, j, k + 1,
ivz) *
sk(i, j, k, 3)
2353 vnm =
w(i, j, k,
ivx) *
sk(i, j, k, 1) &
2354 +
w(i, j, k,
ivy) *
sk(i, j, k, 2) &
2355 +
w(i, j, k,
ivz) *
sk(i, j, k, 3)
2379 porvel = porvel * porflux
2384 qsp = (vnp - sface) * porvel
2385 qsm = (vnm - sface) * porvel
2387 rqsp = qsp *
w(i, j, k + 1,
irho)
2388 rqsm = qsm *
w(i, j, k,
irho)
2394 pa = porflux * (
p(i, j, k + 1) +
p(i, j, k))
2404 fs = rqsp *
w(i, j, k + 1,
ivx) + rqsm *
w(i, j, k,
ivx) &
2405 + pa *
sk(i, j, k, 1)
2406 dw(i, j, k + 1,
imx) =
dw(i, j, k + 1,
imx) - fs
2409 fs = rqsp *
w(i, j, k + 1,
ivy) + rqsm *
w(i, j, k,
ivy) &
2410 + pa *
sk(i, j, k, 2)
2411 dw(i, j, k + 1,
imy) =
dw(i, j, k + 1,
imy) - fs
2414 fs = rqsp *
w(i, j, k + 1,
ivz) + rqsm *
w(i, j, k,
ivz) &
2415 + pa *
sk(i, j, k, 3)
2416 dw(i, j, k + 1,
imz) =
dw(i, j, k + 1,
imz) - fs
2419 fs = qsp *
w(i, j, k + 1,
irhoe) + qsm *
w(i, j, k,
irhoe) &
2420 + porflux * (vnp *
p(i, j, k + 1) + vnm *
p(i, j, k))
2441 rvol =
w(i, j, k,
irho) *
vol(i, j, k)
2443 + rvol * (wwy *
w(i, j, k,
ivz) - wwz *
w(i, j, k,
ivy))
2445 + rvol * (wwz *
w(i, j, k,
ivx) - wwx *
w(i, j, k,
ivz))
2447 + rvol * (wwx *
w(i, j, k,
ivy) - wwy *
w(i, j, k,
ivx))
2474 real(kind=realtype),
parameter :: dpmax = 0.25_realtype
2475 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
2476 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
2477 real(kind=realtype),
parameter :: omega = 0.5_realtype
2478 real(kind=realtype),
parameter :: oneminomega =
one - omega
2482 integer(kind=intType) :: i, j, k, ind, ii
2484 real(kind=realtype) :: plim, sface
2485 real(kind=realtype) :: sfil, fis2, fis4
2486 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53
2487 real(kind=realtype) :: ppor, rrad, dis2, dis4
2488 real(kind=realtype) :: dp1, dp2, tmp, fs
2489 real(kind=realtype) :: ddw1, ddw2, ddw3, ddw4, ddw5, ddw6
2490 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
2491 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
2492 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
2493 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
2494 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
2495 logical :: correctForK
2528 dss(i, j, k, 1) = abs((
p(i + 1, j, k) -
two *
p(i, j, k) +
p(i - 1, j, k)) &
2529 / (omega * (
p(i + 1, j, k) +
two *
p(i, j, k) +
p(i - 1, j, k)) &
2530 + oneminomega * (abs(
p(i + 1, j, k) -
p(i, j, k)) &
2531 + abs(
p(i, j, k) -
p(i - 1, j, k))) + plim))
2533 dss(i, j, k, 2) = abs((
p(i, j + 1, k) -
two *
p(i, j, k) +
p(i, j - 1, k)) &
2534 / (omega * (
p(i, j + 1, k) +
two *
p(i, j, k) +
p(i, j - 1, k)) &
2535 + oneminomega * (abs(
p(i, j + 1, k) -
p(i, j, k)) &
2536 + abs(
p(i, j, k) -
p(i, j - 1, k))) + plim))
2538 dss(i, j, k, 3) = abs((
p(i, j, k + 1) -
two *
p(i, j, k) +
p(i, j, k - 1)) &
2539 / (omega * (
p(i, j, k + 1) +
two *
p(i, j, k) +
p(i, j, k - 1)) &
2540 + oneminomega * (abs(
p(i, j, k + 1) -
p(i, j, k)) &
2541 + abs(
p(i, j, k) -
p(i, j, k - 1))) + plim))
2556 dis2 = ppor * fis2 * min(dpmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1)))
2557 dis4 = dim(ppor * fis4, dis2)
2562 ddw1 =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
2566 ddw2 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
2569 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivx) &
2572 ddw3 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
2575 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivy) &
2578 ddw4 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
2581 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
ivz) &
2595 if (correctfork)
then
2596 ddw6 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
2599 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
itu1) &
2609 gm1 = gammaavg -
one
2618 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
2621 area = sqrt(
si(i, j, k, 1)**2 +
si(i, j, k, 2)**2 +
si(i, j, k, 3)**2)
2622 tmp =
one / max(1.e-25_realtype, area)
2623 sx =
si(i, j, k, 1) * tmp
2624 sy =
si(i, j, k, 2) * tmp
2625 sz =
si(i, j, k, 3) * tmp
2627 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2628 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
2630 unavg = uavg * sx + vavg * sy + wavg * sz
2632 ova2avg =
one / a2avg
2637 sface =
sfacei(i, j, k) * tmp
2643 lam1 = abs(unavg - sface + aavg)
2644 lam2 = abs(unavg - sface - aavg)
2645 lam3 = abs(unavg - sface)
2652 lam1 = max(lam1, epsacoustic * rrad) * area
2653 lam2 = max(lam2, epsacoustic * rrad) * area
2654 lam3 = max(lam3, epsshear * rrad) * area
2659 abv1 =
half * (lam1 + lam2)
2660 abv2 =
half * (lam1 - lam2)
2663 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2664 - wavg * drw + dre) - gm53 * drk
2665 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2667 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2668 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2673 fs = lam3 * dr + abv6
2679 fs = lam3 * dru + uavg * abv6 + sx * abv7
2680 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
2685 fs = lam3 * drv + vavg * abv6 + sy * abv7
2686 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
2691 fs = lam3 * drw + wavg * abv6 + sz * abv7
2692 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
2697 fs = lam3 * dre + havg * abv6 + unavg * abv7
2716 dis2 = ppor * fis2 * min(dpmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2)))
2717 dis4 = dim(ppor * fis4, dis2)
2722 ddw1 =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
2726 ddw2 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
2729 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivx) &
2732 ddw3 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
2735 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivy) &
2738 ddw4 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
2741 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
ivz) &
2755 if (correctfork)
then
2756 ddw6 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
2759 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
itu1) &
2769 gm1 = gammaavg -
one
2778 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
2781 area = sqrt(
sj(i, j, k, 1)**2 +
sj(i, j, k, 2)**2 +
sj(i, j, k, 3)**2)
2782 tmp =
one / max(1.e-25_realtype, area)
2783 sx =
sj(i, j, k, 1) * tmp
2784 sy =
sj(i, j, k, 2) * tmp
2785 sz =
sj(i, j, k, 3) * tmp
2787 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2788 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
2790 unavg = uavg * sx + vavg * sy + wavg * sz
2792 ova2avg =
one / a2avg
2797 sface =
sfacej(i, j, k) * tmp
2803 lam1 = abs(unavg - sface + aavg)
2804 lam2 = abs(unavg - sface - aavg)
2805 lam3 = abs(unavg - sface)
2812 lam1 = max(lam1, epsacoustic * rrad) * area
2813 lam2 = max(lam2, epsacoustic * rrad) * area
2814 lam3 = max(lam3, epsshear * rrad) * area
2819 abv1 =
half * (lam1 + lam2)
2820 abv2 =
half * (lam1 - lam2)
2823 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2824 - wavg * drw + dre) - gm53 * drk
2825 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2827 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2828 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2833 fs = lam3 * dr + abv6
2839 fs = lam3 * dru + uavg * abv6 + sx * abv7
2840 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
2845 fs = lam3 * drv + vavg * abv6 + sy * abv7
2846 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
2851 fs = lam3 * drw + wavg * abv6 + sz * abv7
2852 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
2857 fs = lam3 * dre + havg * abv6 + unavg * abv7
2876 dis2 = ppor * fis2 * min(dpmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3)))
2877 dis4 = dim(ppor * fis4, dis2)
2882 ddw1 =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
2886 ddw2 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
2889 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivx) &
2892 ddw3 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
2895 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivy) &
2898 ddw4 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
2901 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
ivz) &
2915 if (correctfork)
then
2916 ddw6 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
2919 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
itu1) &
2929 gm1 = gammaavg -
one
2938 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
2941 area = sqrt(
sk(i, j, k, 1)**2 +
sk(i, j, k, 2)**2 +
sk(i, j, k, 3)**2)
2942 tmp =
one / max(1.e-25_realtype, area)
2943 sx =
sk(i, j, k, 1) * tmp
2944 sy =
sk(i, j, k, 2) * tmp
2945 sz =
sk(i, j, k, 3) * tmp
2947 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
2948 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
2950 unavg = uavg * sx + vavg * sy + wavg * sz
2952 ova2avg =
one / a2avg
2957 sface =
sfacek(i, j, k) * tmp
2963 lam1 = abs(unavg - sface + aavg)
2964 lam2 = abs(unavg - sface - aavg)
2965 lam3 = abs(unavg - sface)
2972 lam1 = max(lam1, epsacoustic * rrad) * area
2973 lam2 = max(lam2, epsacoustic * rrad) * area
2974 lam3 = max(lam3, epsshear * rrad) * area
2979 abv1 =
half * (lam1 + lam2)
2980 abv2 =
half * (lam1 - lam2)
2983 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
2984 - wavg * drw + dre) - gm53 * drk
2985 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
2987 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
2988 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
2993 fs = lam3 * dr + abv6
2999 fs = lam3 * dru + uavg * abv6 + sx * abv7
3000 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
3005 fs = lam3 * drv + vavg * abv6 + sy * abv7
3006 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
3011 fs = lam3 * drw + wavg * abv6 + sz * abv7
3012 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
3017 fs = lam3 * dre + havg * abv6 + unavg * abv7
3042 real(kind=realtype),
parameter :: dssmax = 0.25_realtype
3043 real(kind=realtype) :: sslim, rhoi
3044 real(kind=realtype) :: sfil, fis2, fis4
3045 real(kind=realtype) :: ppor, rrad, dis2, dis4, fs
3046 real(kind=realtype) :: ddw1, ddw2, ddw3, ddw4, ddw5
3047 integer(kind=intType) :: i, j, k
3083 ss(i, j, k) =
p(i, j, k) / (
w(i, j, k,
irho)**
gamma(i, j, k))
3093 dss(i, j, k, 1) = abs((
ss(i + 1, j, k) -
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
3094 / (
ss(i + 1, j, k) +
two *
ss(i, j, k) +
ss(i - 1, j, k) + sslim))
3096 dss(i, j, k, 2) = abs((
ss(i, j + 1, k) -
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
3097 / (
ss(i, j + 1, k) +
two *
ss(i, j, k) +
ss(i, j - 1, k) + sslim))
3099 dss(i, j, k, 3) = abs((
ss(i, j, k + 1) -
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
3100 / (
ss(i, j, k + 1) +
two *
ss(i, j, k) +
ss(i, j, k - 1) + sslim))
3141 rrad = ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
3143 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1)))
3144 dis4 = dim(fis4 * rrad, dis2)
3150 ddw1 =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
3159 ddw2 =
w(i + 1, j, k,
ivx) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
3161 - dis4 * (
w(i + 2, j, k,
ivx) *
w(i + 2, j, k,
irho) - &
3164 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
3169 ddw3 =
w(i + 1, j, k,
ivy) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
3171 - dis4 * (
w(i + 2, j, k,
ivy) *
w(i + 2, j, k,
irho) - &
3174 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
3179 ddw4 =
w(i + 1, j, k,
ivz) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
3181 - dis4 * (
w(i + 2, j, k,
ivz) *
w(i + 2, j, k,
irho) - &
3184 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
3189 ddw5 = (
w(i + 1, j, k,
irhoe) +
p(i + 1, j, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
3191 - dis4 * ((
w(i + 2, j, k,
irhoe) +
p(i + 2, j, k)) - &
3192 (
w(i - 1, j, k,
irhoe) +
p(i - 1, j, k)) -
three * ddw5)
3210 rrad = ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
3212 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2)))
3213 dis4 = dim(fis4 * rrad, dis2)
3219 ddw1 =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
3228 ddw2 =
w(i, j + 1, k,
ivx) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
3230 - dis4 * (
w(i, j + 2, k,
ivx) *
w(i, j + 2, k,
irho) - &
3233 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
3238 ddw3 =
w(i, j + 1, k,
ivy) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
3240 - dis4 * (
w(i, j + 2, k,
ivy) *
w(i, j + 2, k,
irho) - &
3243 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
3248 ddw4 =
w(i, j + 1, k,
ivz) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
3250 - dis4 * (
w(i, j + 2, k,
ivz) *
w(i, j + 2, k,
irho) - &
3253 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
3258 ddw5 = (
w(i, j + 1, k,
irhoe) +
p(i, j + 1, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
3260 - dis4 * ((
w(i, j + 2, k,
irhoe) +
p(i, j + 2, k)) - &
3261 (
w(i, j - 1, k,
irhoe) +
p(i, j - 1, k)) -
three * ddw5)
3279 rrad = ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
3281 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3)))
3282 dis4 = dim(fis4 * rrad, dis2)
3288 ddw1 =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
3297 ddw2 =
w(i, j, k + 1,
ivx) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
3299 - dis4 * (
w(i, j, k + 2,
ivx) *
w(i, j, k + 2,
irho) - &
3302 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
3307 ddw3 =
w(i, j, k + 1,
ivy) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
3309 - dis4 * (
w(i, j, k + 2,
ivy) *
w(i, j, k + 2,
irho) - &
3312 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
3317 ddw4 =
w(i, j, k + 1,
ivz) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
3319 - dis4 * (
w(i, j, k + 2,
ivz) *
w(i, j, k + 2,
irho) - &
3322 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
3327 ddw5 = (
w(i, j, k + 1,
irhoe) +
p(i, j, k + 1)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
3329 - dis4 * ((
w(i, j, k + 2,
irhoe) +
p(i, j, k + 2)) - &
3330 (
w(i, j, k - 1,
irhoe) +
p(i, j, k - 1)) -
three * ddw5)
3364 logical,
intent(in) :: fineGrid
3368 integer(kind=porType) :: por
3370 integer(kind=intType) :: nwInt
3371 integer(kind=intType) :: i, j, k, ind
3372 integer(kind=intType) :: limUsed, riemannUsed
3374 real(kind=realtype) :: sx, sy, sz, omk, opk, sfil, gammaface
3375 real(kind=realtype) :: factminmod, sface
3377 real(kind=realtype),
dimension(nw) :: left, right
3378 real(kind=realtype),
dimension(nw) :: du1, du2, du3
3379 real(kind=realtype),
dimension(nwf) :: flux
3381 logical :: firstOrderK, correctForK, rotationalPeriodic
3387 rotationalperiodic = .true.
3389 rotationalperiodic = .false.
3424 if (finegrid) limused =
limiter
3429 if (finegrid) riemannused =
riemann
3442 if (correctfork)
then
3445 firstorderk = .true.
3448 firstorderk = .false.
3452 firstorderk = .false.
3473 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
3484 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
3487 right(
ivx) =
w(i + 1, j, k,
ivx)
3488 right(
ivy) =
w(i + 1, j, k,
ivy)
3489 right(
ivz) =
w(i + 1, j, k,
ivz)
3490 right(
irhoe) =
p(i + 1, j, k)
3491 if (correctfork) right(
itu1) =
w(i + 1, j, k,
itu1)
3529 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
3540 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
3543 right(
ivx) =
w(i, j + 1, k,
ivx)
3544 right(
ivy) =
w(i, j + 1, k,
ivy)
3545 right(
ivz) =
w(i, j + 1, k,
ivz)
3546 right(
irhoe) =
p(i, j + 1, k)
3547 if (correctfork) right(
itu1) =
w(i, j + 1, k,
itu1)
3584 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
3595 if (correctfork) left(
itu1) =
w(i, j, k,
itu1)
3598 right(
ivx) =
w(i, j, k + 1,
ivx)
3599 right(
ivy) =
w(i, j, k + 1,
ivy)
3600 right(
ivz) =
w(i, j, k + 1,
ivz)
3601 right(
irhoe) =
p(i, j, k + 1)
3602 if (correctfork) right(
itu1) =
w(i, j, k + 1,
itu1)
3659 du3(
ivx) =
w(i + 2, j, k,
ivx) -
w(i + 1, j, k,
ivx)
3663 du3(
ivy) =
w(i + 2, j, k,
ivy) -
w(i + 1, j, k,
ivy)
3667 du3(
ivz) =
w(i + 2, j, k,
ivz) -
w(i + 1, j, k,
ivz)
3669 du1(
irhoe) =
p(i, j, k) -
p(i - 1, j, k)
3670 du2(
irhoe) =
p(i + 1, j, k) -
p(i, j, k)
3671 du3(
irhoe) =
p(i + 2, j, k) -
p(i + 1, j, k)
3673 if (correctfork)
then
3695 right(
ivx) = right(
ivx) +
w(i + 1, j, k,
ivx)
3696 right(
ivy) = right(
ivy) +
w(i + 1, j, k,
ivy)
3697 right(
ivz) = right(
ivz) +
w(i + 1, j, k,
ivz)
3700 if (correctfork)
then
3708 sx =
si(i, j, k, 1); sy =
si(i, j, k, 2); sz =
si(i, j, k, 3)
3754 du3(
ivx) =
w(i, j + 2, k,
ivx) -
w(i, j + 1, k,
ivx)
3758 du3(
ivy) =
w(i, j + 2, k,
ivy) -
w(i, j + 1, k,
ivy)
3762 du3(
ivz) =
w(i, j + 2, k,
ivz) -
w(i, j + 1, k,
ivz)
3764 du1(
irhoe) =
p(i, j, k) -
p(i, j - 1, k)
3765 du2(
irhoe) =
p(i, j + 1, k) -
p(i, j, k)
3766 du3(
irhoe) =
p(i, j + 2, k) -
p(i, j + 1, k)
3768 if (correctfork)
then
3790 right(
ivx) = right(
ivx) +
w(i, j + 1, k,
ivx)
3791 right(
ivy) = right(
ivy) +
w(i, j + 1, k,
ivy)
3792 right(
ivz) = right(
ivz) +
w(i, j + 1, k,
ivz)
3795 if (correctfork)
then
3803 sx =
sj(i, j, k, 1); sy =
sj(i, j, k, 2); sz =
sj(i, j, k, 3)
3848 du3(
ivx) =
w(i, j, k + 2,
ivx) -
w(i, j, k + 1,
ivx)
3852 du3(
ivy) =
w(i, j, k + 2,
ivy) -
w(i, j, k + 1,
ivy)
3856 du3(
ivz) =
w(i, j, k + 2,
ivz) -
w(i, j, k + 1,
ivz)
3858 du1(
irhoe) =
p(i, j, k) -
p(i, j, k - 1)
3859 du2(
irhoe) =
p(i, j, k + 1) -
p(i, j, k)
3860 du3(
irhoe) =
p(i, j, k + 2) -
p(i, j, k + 1)
3862 if (correctfork)
then
3884 right(
ivx) = right(
ivx) +
w(i, j, k + 1,
ivx)
3885 right(
ivy) = right(
ivy) +
w(i, j, k + 1,
ivy)
3886 right(
ivz) = right(
ivz) +
w(i, j, k + 1,
ivz)
3889 if (correctfork)
then
3897 sx =
sk(i, j, k, 1); sy =
sk(i, j, k, 2); sz =
sk(i, j, k, 3)
3947 real(kind=realtype),
parameter :: epslim = 1.e-10_realtype
3951 real(kind=realtype),
dimension(:),
intent(inout) :: du1, du2, du3
3952 real(kind=realtype),
dimension(:),
intent(out) :: left, right
3954 real(kind=realtype),
dimension(:, :, :, :, :),
pointer :: rotmatrix
3958 integer(kind=intType) :: l
3960 real(kind=realtype) :: rl1, rl2, rr1, rr2, tmp, dvx, dvy, dvz
3962 real(kind=realtype),
dimension(3, 3) :: rot
3967 if (rotationalperiodic)
then
3972 rot(1, 1) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1, 1)
3973 rot(1, 2) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1, 2)
3974 rot(1, 3) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 1, 3)
3976 rot(2, 1) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 2, 1)
3977 rot(2, 2) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 2, 2)
3978 rot(2, 3) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 2, 3)
3980 rot(3, 1) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 3, 1)
3981 rot(3, 2) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 3, 2)
3982 rot(3, 3) = rotmatrix(i +
ii - 2, j +
jj - 2, k +
kk - 2, 3, 3)
3987 dvx = du1(
ivx); dvy = du1(
ivy); dvz = du1(
ivz)
3988 du1(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
3989 du1(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
3990 du1(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
3992 dvx = du2(
ivx); dvy = du2(
ivy); dvz = du2(
ivz)
3993 du2(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
3994 du2(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
3995 du2(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
3997 dvx = du3(
ivx); dvy = du3(
ivy); dvz = du3(
ivz)
3998 du3(
ivx) = rot(1, 1) * dvx + rot(1, 2) * dvy + rot(1, 3) * dvz
3999 du3(
ivy) = rot(2, 1) * dvx + rot(2, 2) * dvy + rot(2, 3) * dvz
4000 du3(
ivz) = rot(3, 1) * dvx + rot(3, 2) * dvy + rot(3, 3) * dvz
4006 select case (limused)
4014 left(l) = omk * du1(l) + opk * du2(l)
4015 right(l) = -omk * du3(l) - opk * du2(l)
4030 tmp =
one / sign(max(abs(du2(l)), epslim), du2(l))
4032 du2(l) / sign(max(abs(du1(l)), epslim), du1(l)))
4033 rl2 = max(
zero, du1(l) * tmp)
4035 rr1 = max(
zero, du3(l) * tmp)
4037 du2(l) / sign(max(abs(du3(l)), epslim), du3(l)))
4041 rl1 = rl1 * (rl1 +
one) / (rl1 * rl1 +
one)
4042 rl2 = rl2 * (rl2 +
one) / (rl2 * rl2 +
one)
4043 rr1 = rr1 * (rr1 +
one) / (rr1 * rr1 +
one)
4044 rr2 = rr2 * (rr2 +
one) / (rr2 * rr2 +
one)
4049 left(l) = omk * rl1 * du1(l) + opk * rl2 * du2(l)
4050 right(l) = -opk * rr1 * du2(l) - omk * rr2 * du3(l)
4066 tmp =
one / sign(max(abs(du2(l)), epslim), du2(l))
4068 du2(l) / sign(max(abs(du1(l)), epslim), du1(l)))
4069 rl2 = max(
zero, du1(l) * tmp)
4071 rr1 = max(
zero, du3(l) * tmp)
4073 du2(l) / sign(max(abs(du3(l)), epslim), du3(l)))
4077 rl1 = min(
one, factminmod * rl1)
4078 rl2 = min(
one, factminmod * rl2)
4079 rr1 = min(
one, factminmod * rr1)
4080 rr2 = min(
one, factminmod * rr2)
4085 left(l) = omk * rl1 * du1(l) + opk * rl2 * du2(l)
4086 right(l) = -opk * rr1 * du2(l) - omk * rr2 * du3(l)
4096 if (firstorderk)
then
4105 if (rotationalperiodic)
then
4109 dvx = left(
ivx); dvy = left(
ivy); dvz = left(
ivz)
4110 left(
ivx) = rot(1, 1) * dvx + rot(2, 1) * dvy + rot(3, 1) * dvz
4111 left(
ivy) = rot(1, 2) * dvx + rot(2, 2) * dvy + rot(3, 2) * dvz
4112 left(
ivz) = rot(1, 3) * dvx + rot(2, 3) * dvy + rot(3, 3) * dvz
4116 dvx = right(
ivx); dvy = right(
ivy); dvz = right(
ivz)
4117 right(
ivx) = rot(1, 1) * dvx + rot(2, 1) * dvy + rot(3, 1) * dvz
4118 right(
ivy) = rot(1, 2) * dvx + rot(2, 2) * dvy + rot(3, 2) * dvz
4119 right(
ivz) = rot(1, 3) * dvx + rot(2, 3) * dvy + rot(3, 3) * dvz
4136 real(kind=realtype),
dimension(*),
intent(in) :: left, right
4137 real(kind=realtype),
dimension(*),
intent(out) :: flux
4141 real(kind=realtype) :: porflux, rface
4142 real(kind=realtype) :: etl, etr, z1l, z1r, tmp
4143 real(kind=realtype) :: dr, dru, drv, drw, dre, drk
4144 real(kind=realtype) :: ravg, uavg, vavg, wavg, havg, kavg
4145 real(kind=realtype) :: alphaavg, a2avg, aavg, unavg
4146 real(kind=realtype) :: ovaavg, ova2avg, area, eta
4147 real(kind=realtype) :: gm1, gm53
4148 real(kind=realtype) :: lam1, lam2, lam3
4149 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
4150 real(kind=realtype),
dimension(2) :: ktmp
4160 gm1 = gammaface -
one
4165 select case (riemannused)
4181 z1l = sqrt(left(
irho))
4182 z1r = sqrt(right(
irho))
4183 tmp =
one / (z1l + z1r)
4188 if (correctfork)
then
4193 ktmp(1) = left(
itu1)
4194 ktmp(2) = right(
itu1)
4204 kavg = tmp * (z1l * left(
itu1) + z1r * right(
itu1))
4219 left(
irhoe), ktmp(1), etl, correctfork)
4222 right(
irhoe), ktmp(2), etr, correctfork)
4236 ravg =
fourth * (z1r + z1l)**2
4237 uavg = tmp * (z1l * left(
ivx) + z1r * right(
ivx))
4238 vavg = tmp * (z1l * left(
ivy) + z1r * right(
ivy))
4239 wavg = tmp * (z1l * left(
ivz) + z1r * right(
ivz))
4240 havg = tmp * ((etl + left(
irhoe)) / z1l &
4241 + (etr + right(
irhoe)) / z1r)
4246 area = sqrt(sx**2 + sy**2 + sz**2)
4247 tmp =
one / max(1.e-25_realtype, area)
4256 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4257 a2avg = abs(gm1 * (havg - alphaavg) - gm53 * kavg)
4259 unavg = uavg * sx + vavg * sy + wavg * sz
4262 ova2avg =
one / a2avg
4280 eta =
half * (abs((left(
ivx) - right(
ivx)) * sx &
4281 + (left(
ivy) - right(
ivy)) * sy &
4282 + (left(
ivz) - right(
ivz)) * sz) &
4283 + abs(sqrt(gammaface * left(
irhoe) / left(
irho)) &
4284 - sqrt(gammaface * right(
irhoe) / right(
irho))))
4288 lam1 = abs(unavg - rface + aavg)
4289 lam2 = abs(unavg - rface - aavg)
4290 lam3 = abs(unavg - rface)
4295 if (lam1 < tmp) lam1 = eta +
fourth * lam1 * lam1 / eta
4296 if (lam2 < tmp) lam2 = eta +
fourth * lam2 * lam2 / eta
4297 if (lam3 < tmp) lam3 = eta +
fourth * lam3 * lam3 / eta
4309 abv1 =
half * (lam1 + lam2)
4310 abv2 =
half * (lam1 - lam2)
4313 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4314 - wavg * drw + dre) - gm53 * drk
4315 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4317 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4318 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4324 flux(
irho) = -porflux * (lam3 * dr + abv6)
4325 flux(
imx) = -porflux * (lam3 * dru + uavg * abv6 &
4327 flux(
imy) = -porflux * (lam3 * drv + vavg * abv6 &
4329 flux(
imz) = -porflux * (lam3 * drw + wavg * abv6 &
4331 flux(
irhoe) = -porflux * (lam3 * dre + havg * abv6 &
4345 "Turkel preconditioner not implemented yet")
4349 "choi merkle preconditioner not implemented yet")
4354 call terminate(
"riemannFlux",
"van leer fvs not implemented yet")
4357 call terminate(
"riemannFlux",
"ausmdv fvs not implemented yet")
4380 real(kind=realtype),
parameter :: dssmax = 0.25_realtype
4381 real(kind=realtype) :: sslim, rhoi
4382 real(kind=realtype) :: sfil, fis2, fis4
4383 real(kind=realtype) :: ppor, rrad, dis2, dis4, fs
4384 real(kind=realtype) :: ddw
4385 integer(kind=intType) :: i, j, k
4413 dss(i, j, k, 1) = abs((
ss(i + 1, j, k) -
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
4414 / (
ss(i + 1, j, k) +
two *
ss(i, j, k) +
ss(i - 1, j, k) + sslim))
4416 dss(i, j, k, 2) = abs((
ss(i, j + 1, k) -
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
4417 / (
ss(i, j + 1, k) +
two *
ss(i, j, k) +
ss(i, j - 1, k) + sslim))
4419 dss(i, j, k, 3) = abs((
ss(i, j, k + 1) -
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
4420 / (
ss(i, j, k + 1) +
two *
ss(i, j, k) +
ss(i, j, k - 1) + sslim))
4451 rrad = ppor * (
radi(i, j, k) +
radi(i + 1, j, k))
4453 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1))) +
sigma * fis4 * rrad
4459 ddw =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
4467 ddw =
w(i + 1, j, k,
ivx) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
4470 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
4475 ddw =
w(i + 1, j, k,
ivy) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
4478 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
4482 ddw =
w(i + 1, j, k,
ivz) *
w(i + 1, j, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
4485 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
4489 ddw = (
w(i + 1, j, k,
irhoe) +
p(i + 1, j, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
4509 rrad = ppor * (
radj(i, j, k) +
radj(i, j + 1, k))
4511 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2))) +
sigma * fis4 * rrad
4517 ddw =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
4525 ddw =
w(i, j + 1, k,
ivx) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
4528 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
4533 ddw =
w(i, j + 1, k,
ivy) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
4536 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
4541 ddw =
w(i, j + 1, k,
ivz) *
w(i, j + 1, k,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
4544 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
4549 ddw = (
w(i, j + 1, k,
irhoe) +
p(i, j + 1, k)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
4568 rrad = ppor * (
radk(i, j, k) +
radk(i, j, k + 1))
4570 dis2 = fis2 * rrad * min(dssmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3))) +
sigma * fis4 * rrad
4576 ddw =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
4584 ddw =
w(i, j, k + 1,
ivx) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivx) *
w(i, j, k,
irho)
4587 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
4592 ddw =
w(i, j, k + 1,
ivy) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivy) *
w(i, j, k,
irho)
4595 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
4600 ddw =
w(i, j, k + 1,
ivz) *
w(i, j, k + 1,
irho) -
w(i, j, k,
ivz) *
w(i, j, k,
irho)
4603 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
4607 ddw = (
w(i, j, k + 1,
irhoe) +
p(i, j, k + 1)) - (
w(i, j, k,
irhoe) +
p(i, j, k))
4636 real(kind=realtype),
parameter :: dpmax = 0.25_realtype
4637 real(kind=realtype),
parameter :: epsacoustic = 0.25_realtype
4638 real(kind=realtype),
parameter :: epsshear = 0.025_realtype
4639 real(kind=realtype),
parameter :: omega = 0.5_realtype
4640 real(kind=realtype),
parameter :: oneminomega =
one - omega
4644 integer(kind=intType) :: i, j, k, ind, ii
4646 real(kind=realtype) :: plim, sface
4647 real(kind=realtype) :: sfil, fis2, fis4
4648 real(kind=realtype) :: gammaavg, gm1, ovgm1, gm53
4649 real(kind=realtype) :: ppor, rrad, dis2, dis4
4650 real(kind=realtype) :: dp1, dp2, tmp, fs
4651 real(kind=realtype) :: ddw, ddw6
4652 real(kind=realtype) :: dr, dru, drv, drw, dre, drk, sx, sy, sz
4653 real(kind=realtype) :: uavg, vavg, wavg, a2avg, aavg, havg
4654 real(kind=realtype) :: alphaavg, unavg, ovaavg, ova2avg
4655 real(kind=realtype) :: kavg, lam1, lam2, lam3, area
4656 real(kind=realtype) :: abv1, abv2, abv3, abv4, abv5, abv6, abv7
4657 logical :: correctForK
4690 dss(i, j, k, 1) = abs((
ss(i + 1, j, k) -
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
4691 / (omega * (
ss(i + 1, j, k) +
two *
ss(i, j, k) +
ss(i - 1, j, k)) &
4692 + oneminomega * (abs(
ss(i + 1, j, k) -
ss(i, j, k)) &
4693 + abs(
ss(i, j, k) -
ss(i - 1, j, k))) + plim))
4695 dss(i, j, k, 2) = abs((
ss(i, j + 1, k) -
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
4696 / (omega * (
ss(i, j + 1, k) +
two *
ss(i, j, k) +
ss(i, j - 1, k)) &
4697 + oneminomega * (abs(
ss(i, j + 1, k) -
ss(i, j, k)) &
4698 + abs(
ss(i, j, k) -
ss(i, j - 1, k))) + plim))
4700 dss(i, j, k, 3) = abs((
ss(i, j, k + 1) -
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
4701 / (omega * (
ss(i, j, k + 1) +
two *
ss(i, j, k) +
ss(i, j, k - 1)) &
4702 + oneminomega * (abs(
ss(i, j, k + 1) -
ss(i, j, k)) &
4703 + abs(
ss(i, j, k) -
ss(i, j, k - 1))) + plim))
4719 dis2 = fis2 * ppor * min(dpmax, max(
dss(i, j, k, 1),
dss(i + 1, j, k, 1))) &
4720 +
sigma * fis4 * ppor
4725 ddw =
w(i + 1, j, k,
irho) -
w(i, j, k,
irho)
4728 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivx) &
4732 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivy) &
4736 ddw =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
ivz) &
4750 if (correctfork)
then
4751 ddw6 =
w(i + 1, j, k,
irho) *
w(i + 1, j, k,
itu1) &
4754 - dis4 * (
w(i + 2, j, k,
irho) *
w(i + 2, j, k,
itu1) &
4764 gm1 = gammaavg -
one
4773 a2avg =
half * (
gamma(i + 1, j, k) *
p(i + 1, j, k) /
w(i + 1, j, k,
irho) &
4776 area = sqrt(
si(i, j, k, 1)**2 +
si(i, j, k, 2)**2 +
si(i, j, k, 3)**2)
4777 tmp =
one / max(1.e-25_realtype, area)
4778 sx =
si(i, j, k, 1) * tmp
4779 sy =
si(i, j, k, 2) * tmp
4780 sz =
si(i, j, k, 3) * tmp
4782 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4783 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4785 unavg = uavg * sx + vavg * sy + wavg * sz
4787 ova2avg =
one / a2avg
4792 sface =
sfacei(i, j, k) * tmp
4798 lam1 = abs(unavg - sface + aavg)
4799 lam2 = abs(unavg - sface - aavg)
4800 lam3 = abs(unavg - sface)
4807 lam1 = max(lam1, epsacoustic * rrad) * area
4808 lam2 = max(lam2, epsacoustic * rrad) * area
4809 lam3 = max(lam3, epsshear * rrad) * area
4814 abv1 =
half * (lam1 + lam2)
4815 abv2 =
half * (lam1 - lam2)
4818 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4819 - wavg * drw + dre) - gm53 * drk
4820 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4822 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4823 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4828 fs = lam3 * dr + abv6
4834 fs = lam3 * dru + uavg * abv6 + sx * abv7
4835 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fs
4840 fs = lam3 * drv + vavg * abv6 + sy * abv7
4841 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fs
4846 fs = lam3 * drw + wavg * abv6 + sz * abv7
4847 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fs
4852 fs = lam3 * dre + havg * abv6 + unavg * abv7
4871 dis2 = fis2 * ppor * min(dpmax, max(
dss(i, j, k, 2),
dss(i, j + 1, k, 2))) &
4872 +
sigma * fis4 * ppor
4877 ddw =
w(i, j + 1, k,
irho) -
w(i, j, k,
irho)
4880 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivx) &
4884 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivy) &
4888 ddw =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
ivz) &
4902 if (correctfork)
then
4903 ddw6 =
w(i, j + 1, k,
irho) *
w(i, j + 1, k,
itu1) &
4906 - dis4 * (
w(i, j + 2, k,
irho) *
w(i, j + 2, k,
itu1) &
4916 gm1 = gammaavg -
one
4925 a2avg =
half * (
gamma(i, j + 1, k) *
p(i, j + 1, k) /
w(i, j + 1, k,
irho) &
4928 area = sqrt(
sj(i, j, k, 1)**2 +
sj(i, j, k, 2)**2 +
sj(i, j, k, 3)**2)
4929 tmp =
one / max(1.e-25_realtype, area)
4930 sx =
sj(i, j, k, 1) * tmp
4931 sy =
sj(i, j, k, 2) * tmp
4932 sz =
sj(i, j, k, 3) * tmp
4934 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
4935 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
4937 unavg = uavg * sx + vavg * sy + wavg * sz
4939 ova2avg =
one / a2avg
4944 sface =
sfacej(i, j, k) * tmp
4950 lam1 = abs(unavg - sface + aavg)
4951 lam2 = abs(unavg - sface - aavg)
4952 lam3 = abs(unavg - sface)
4959 lam1 = max(lam1, epsacoustic * rrad) * area
4960 lam2 = max(lam2, epsacoustic * rrad) * area
4961 lam3 = max(lam3, epsshear * rrad) * area
4966 abv1 =
half * (lam1 + lam2)
4967 abv2 =
half * (lam1 - lam2)
4970 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
4971 - wavg * drw + dre) - gm53 * drk
4972 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
4974 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
4975 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
4980 fs = lam3 * dr + abv6
4986 fs = lam3 * dru + uavg * abv6 + sx * abv7
4987 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fs
4992 fs = lam3 * drv + vavg * abv6 + sy * abv7
4993 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fs
4998 fs = lam3 * drw + wavg * abv6 + sz * abv7
4999 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fs
5004 fs = lam3 * dre + havg * abv6 + unavg * abv7
5023 dis2 = fis2 * ppor * min(dpmax, max(
dss(i, j, k, 3),
dss(i, j, k + 1, 3))) &
5024 +
sigma * fis4 * ppor
5029 ddw =
w(i, j, k + 1,
irho) -
w(i, j, k,
irho)
5032 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivx) &
5036 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivy) &
5040 ddw =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
ivz) &
5054 if (correctfork)
then
5055 ddw6 =
w(i, j, k + 1,
irho) *
w(i, j, k + 1,
itu1) &
5058 - dis4 * (
w(i, j, k + 2,
irho) *
w(i, j, k + 2,
itu1) &
5068 gm1 = gammaavg -
one
5077 a2avg =
half * (
gamma(i, j, k + 1) *
p(i, j, k + 1) /
w(i, j, k + 1,
irho) &
5080 area = sqrt(
sk(i, j, k, 1)**2 +
sk(i, j, k, 2)**2 +
sk(i, j, k, 3)**2)
5081 tmp =
one / max(1.e-25_realtype, area)
5082 sx =
sk(i, j, k, 1) * tmp
5083 sy =
sk(i, j, k, 2) * tmp
5084 sz =
sk(i, j, k, 3) * tmp
5086 alphaavg =
half * (uavg**2 + vavg**2 + wavg**2)
5087 havg = alphaavg + ovgm1 * (a2avg - gm53 * kavg)
5089 unavg = uavg * sx + vavg * sy + wavg * sz
5091 ova2avg =
one / a2avg
5096 sface =
sfacek(i, j, k) * tmp
5102 lam1 = abs(unavg - sface + aavg)
5103 lam2 = abs(unavg - sface - aavg)
5104 lam3 = abs(unavg - sface)
5111 lam1 = max(lam1, epsacoustic * rrad) * area
5112 lam2 = max(lam2, epsacoustic * rrad) * area
5113 lam3 = max(lam3, epsshear * rrad) * area
5118 abv1 =
half * (lam1 + lam2)
5119 abv2 =
half * (lam1 - lam2)
5122 abv4 = gm1 * (alphaavg * dr - uavg * dru - vavg * drv &
5123 - wavg * drw + dre) - gm53 * drk
5124 abv5 = sx * dru + sy * drv + sz * drw - unavg * dr
5126 abv6 = abv3 * abv4 * ova2avg + abv2 * abv5 * ovaavg
5127 abv7 = abv2 * abv4 * ovaavg + abv3 * abv5
5132 fs = lam3 * dr + abv6
5138 fs = lam3 * dru + uavg * abv6 + sx * abv7
5139 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fs
5144 fs = lam3 * drv + vavg * abv6 + sy * abv7
5145 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fs
5150 fs = lam3 * drw + wavg * abv6 + sz * abv7
5151 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fs
5156 fs = lam3 * dre + havg * abv6 + unavg * abv7
5175 logical :: correctForK
5176 real(kind=realtype) :: pp
5177 real(kind=realtype),
parameter :: twothird =
two *
third
5178 integer(kind=intType) :: i, j, k
5183 if (correctfork)
then
5187 pp =
p(i, j, k) - twothird *
w(i, j, k,
irho) *
w(i, j, k,
itu1)
5196 aa(i, j, k) =
gamma(i, j, k) *
p(i, j, k) /
w(i, j, k,
irho)
5212 real(kind=realtype) :: a2, ovol, ubar, vbar, wbar, sx, sy, sz
5213 integer(kind=intType) :: i, j, k
5245 sx =
sk(i, j, k - 1, 1) +
sk(i + 1, j, k - 1, 1) &
5246 +
sk(i, j + 1, k - 1, 1) +
sk(i + 1, j + 1, k - 1, 1) &
5247 +
sk(i, j, k, 1) +
sk(i + 1, j, k, 1) &
5248 +
sk(i, j + 1, k, 1) +
sk(i + 1, j + 1, k, 1)
5249 sy =
sk(i, j, k - 1, 2) +
sk(i + 1, j, k - 1, 2) &
5250 +
sk(i, j + 1, k - 1, 2) +
sk(i + 1, j + 1, k - 1, 2) &
5251 +
sk(i, j, k, 2) +
sk(i + 1, j, k, 2) &
5252 +
sk(i, j + 1, k, 2) +
sk(i + 1, j + 1, k, 2)
5253 sz =
sk(i, j, k - 1, 3) +
sk(i + 1, j, k - 1, 3) &
5254 +
sk(i, j + 1, k - 1, 3) +
sk(i + 1, j + 1, k - 1, 3) &
5255 +
sk(i, j, k, 3) +
sk(i + 1, j, k, 3) &
5256 +
sk(i, j + 1, k, 3) +
sk(i + 1, j + 1, k, 3)
5263 +
w(i, j + 1, k,
ivx) +
w(i + 1, j + 1, k,
ivx))
5265 +
w(i, j + 1, k,
ivy) +
w(i + 1, j + 1, k,
ivy))
5267 +
w(i, j + 1, k,
ivz) +
w(i + 1, j + 1, k,
ivz))
5269 a2 =
fourth * (
aa(i, j, k) +
aa(i + 1, j, k) +
aa(i, j + 1, k) +
aa(i + 1, j + 1, k))
5277 ux(i, j, k - 1) =
ux(i, j, k - 1) + ubar * sx
5278 uy(i, j, k - 1) =
uy(i, j, k - 1) + ubar * sy
5279 uz(i, j, k - 1) =
uz(i, j, k - 1) + ubar * sz
5281 vx(i, j, k - 1) =
vx(i, j, k - 1) + vbar * sx
5282 vy(i, j, k - 1) =
vy(i, j, k - 1) + vbar * sy
5283 vz(i, j, k - 1) =
vz(i, j, k - 1) + vbar * sz
5285 wx(i, j, k - 1) =
wx(i, j, k - 1) + wbar * sx
5286 wy(i, j, k - 1) =
wy(i, j, k - 1) + wbar * sy
5287 wz(i, j, k - 1) =
wz(i, j, k - 1) + wbar * sz
5289 qx(i, j, k - 1) =
qx(i, j, k - 1) - a2 * sx
5290 qy(i, j, k - 1) =
qy(i, j, k - 1) - a2 * sy
5291 qz(i, j, k - 1) =
qz(i, j, k - 1) - a2 * sz
5295 ux(i, j, k) =
ux(i, j, k) - ubar * sx
5296 uy(i, j, k) =
uy(i, j, k) - ubar * sy
5297 uz(i, j, k) =
uz(i, j, k) - ubar * sz
5299 vx(i, j, k) =
vx(i, j, k) - vbar * sx
5300 vy(i, j, k) =
vy(i, j, k) - vbar * sy
5301 vz(i, j, k) =
vz(i, j, k) - vbar * sz
5303 wx(i, j, k) =
wx(i, j, k) - wbar * sx
5304 wy(i, j, k) =
wy(i, j, k) - wbar * sy
5305 wz(i, j, k) =
wz(i, j, k) - wbar * sz
5307 qx(i, j, k) =
qx(i, j, k) + a2 * sx
5308 qy(i, j, k) =
qy(i, j, k) + a2 * sy
5309 qz(i, j, k) =
qz(i, j, k) + a2 * sz
5327 sx =
sj(i, j - 1, k, 1) +
sj(i + 1, j - 1, k, 1) &
5328 +
sj(i, j - 1, k + 1, 1) +
sj(i + 1, j - 1, k + 1, 1) &
5329 +
sj(i, j, k, 1) +
sj(i + 1, j, k, 1) &
5330 +
sj(i, j, k + 1, 1) +
sj(i + 1, j, k + 1, 1)
5331 sy =
sj(i, j - 1, k, 2) +
sj(i + 1, j - 1, k, 2) &
5332 +
sj(i, j - 1, k + 1, 2) +
sj(i + 1, j - 1, k + 1, 2) &
5333 +
sj(i, j, k, 2) +
sj(i + 1, j, k, 2) &
5334 +
sj(i, j, k + 1, 2) +
sj(i + 1, j, k + 1, 2)
5335 sz =
sj(i, j - 1, k, 3) +
sj(i + 1, j - 1, k, 3) &
5336 +
sj(i, j - 1, k + 1, 3) +
sj(i + 1, j - 1, k + 1, 3) &
5337 +
sj(i, j, k, 3) +
sj(i + 1, j, k, 3) &
5338 +
sj(i, j, k + 1, 3) +
sj(i + 1, j, k + 1, 3)
5345 +
w(i, j, k + 1,
ivx) +
w(i + 1, j, k + 1,
ivx))
5347 +
w(i, j, k + 1,
ivy) +
w(i + 1, j, k + 1,
ivy))
5349 +
w(i, j, k + 1,
ivz) +
w(i + 1, j, k + 1,
ivz))
5351 a2 =
fourth * (
aa(i, j, k) +
aa(i + 1, j, k) +
aa(i, j, k + 1) +
aa(i + 1, j, k + 1))
5359 ux(i, j - 1, k) =
ux(i, j - 1, k) + ubar * sx
5360 uy(i, j - 1, k) =
uy(i, j - 1, k) + ubar * sy
5361 uz(i, j - 1, k) =
uz(i, j - 1, k) + ubar * sz
5363 vx(i, j - 1, k) =
vx(i, j - 1, k) + vbar * sx
5364 vy(i, j - 1, k) =
vy(i, j - 1, k) + vbar * sy
5365 vz(i, j - 1, k) =
vz(i, j - 1, k) + vbar * sz
5367 wx(i, j - 1, k) =
wx(i, j - 1, k) + wbar * sx
5368 wy(i, j - 1, k) =
wy(i, j - 1, k) + wbar * sy
5369 wz(i, j - 1, k) =
wz(i, j - 1, k) + wbar * sz
5371 qx(i, j - 1, k) =
qx(i, j - 1, k) - a2 * sx
5372 qy(i, j - 1, k) =
qy(i, j - 1, k) - a2 * sy
5373 qz(i, j - 1, k) =
qz(i, j - 1, k) - a2 * sz
5377 ux(i, j, k) =
ux(i, j, k) - ubar * sx
5378 uy(i, j, k) =
uy(i, j, k) - ubar * sy
5379 uz(i, j, k) =
uz(i, j, k) - ubar * sz
5381 vx(i, j, k) =
vx(i, j, k) - vbar * sx
5382 vy(i, j, k) =
vy(i, j, k) - vbar * sy
5383 vz(i, j, k) =
vz(i, j, k) - vbar * sz
5385 wx(i, j, k) =
wx(i, j, k) - wbar * sx
5386 wy(i, j, k) =
wy(i, j, k) - wbar * sy
5387 wz(i, j, k) =
wz(i, j, k) - wbar * sz
5389 qx(i, j, k) =
qx(i, j, k) + a2 * sx
5390 qy(i, j, k) =
qy(i, j, k) + a2 * sy
5391 qz(i, j, k) =
qz(i, j, k) + a2 * sz
5409 sx =
si(i - 1, j, k, 1) +
si(i - 1, j + 1, k, 1) &
5410 +
si(i - 1, j, k + 1, 1) +
si(i - 1, j + 1, k + 1, 1) &
5411 +
si(i, j, k, 1) +
si(i, j + 1, k, 1) &
5412 +
si(i, j, k + 1, 1) +
si(i, j + 1, k + 1, 1)
5413 sy =
si(i - 1, j, k, 2) +
si(i - 1, j + 1, k, 2) &
5414 +
si(i - 1, j, k + 1, 2) +
si(i - 1, j + 1, k + 1, 2) &
5415 +
si(i, j, k, 2) +
si(i, j + 1, k, 2) &
5416 +
si(i, j, k + 1, 2) +
si(i, j + 1, k + 1, 2)
5417 sz =
si(i - 1, j, k, 3) +
si(i - 1, j + 1, k, 3) &
5418 +
si(i - 1, j, k + 1, 3) +
si(i - 1, j + 1, k + 1, 3) &
5419 +
si(i, j, k, 3) +
si(i, j + 1, k, 3) &
5420 +
si(i, j, k + 1, 3) +
si(i, j + 1, k + 1, 3)
5427 +
w(i, j, k + 1,
ivx) +
w(i, j + 1, k + 1,
ivx))
5429 +
w(i, j, k + 1,
ivy) +
w(i, j + 1, k + 1,
ivy))
5431 +
w(i, j, k + 1,
ivz) +
w(i, j + 1, k + 1,
ivz))
5433 a2 =
fourth * (
aa(i, j, k) +
aa(i, j + 1, k) +
aa(i, j, k + 1) +
aa(i, j + 1, k + 1))
5441 ux(i - 1, j, k) =
ux(i - 1, j, k) + ubar * sx
5442 uy(i - 1, j, k) =
uy(i - 1, j, k) + ubar * sy
5443 uz(i - 1, j, k) =
uz(i - 1, j, k) + ubar * sz
5445 vx(i - 1, j, k) =
vx(i - 1, j, k) + vbar * sx
5446 vy(i - 1, j, k) =
vy(i - 1, j, k) + vbar * sy
5447 vz(i - 1, j, k) =
vz(i - 1, j, k) + vbar * sz
5449 wx(i - 1, j, k) =
wx(i - 1, j, k) + wbar * sx
5450 wy(i - 1, j, k) =
wy(i - 1, j, k) + wbar * sy
5451 wz(i - 1, j, k) =
wz(i - 1, j, k) + wbar * sz
5453 qx(i - 1, j, k) =
qx(i - 1, j, k) - a2 * sx
5454 qy(i - 1, j, k) =
qy(i - 1, j, k) - a2 * sy
5455 qz(i - 1, j, k) =
qz(i - 1, j, k) - a2 * sz
5459 ux(i, j, k) =
ux(i, j, k) - ubar * sx
5460 uy(i, j, k) =
uy(i, j, k) - ubar * sy
5461 uz(i, j, k) =
uz(i, j, k) - ubar * sz
5463 vx(i, j, k) =
vx(i, j, k) - vbar * sx
5464 vy(i, j, k) =
vy(i, j, k) - vbar * sy
5465 vz(i, j, k) =
vz(i, j, k) - vbar * sz
5467 wx(i, j, k) =
wx(i, j, k) - wbar * sx
5468 wy(i, j, k) =
wy(i, j, k) - wbar * sy
5469 wz(i, j, k) =
wz(i, j, k) - wbar * sz
5471 qx(i, j, k) =
qx(i, j, k) + a2 * sx
5472 qy(i, j, k) =
qy(i, j, k) + a2 * sy
5473 qz(i, j, k) =
qz(i, j, k) + a2 * sz
5487 ovol =
one / (
vol(i, j, k) +
vol(i, j, k + 1) &
5488 +
vol(i + 1, j, k) +
vol(i + 1, j, k + 1) &
5489 +
vol(i, j + 1, k) +
vol(i, j + 1, k + 1) &
5490 +
vol(i + 1, j + 1, k) +
vol(i + 1, j + 1, k + 1))
5495 ux(i, j, k) =
ux(i, j, k) * ovol
5496 uy(i, j, k) =
uy(i, j, k) * ovol
5497 uz(i, j, k) =
uz(i, j, k) * ovol
5499 vx(i, j, k) =
vx(i, j, k) * ovol
5500 vy(i, j, k) =
vy(i, j, k) * ovol
5501 vz(i, j, k) =
vz(i, j, k) * ovol
5503 wx(i, j, k) =
wx(i, j, k) * ovol
5504 wy(i, j, k) =
wy(i, j, k) * ovol
5505 wz(i, j, k) =
wz(i, j, k) * ovol
5507 qx(i, j, k) =
qx(i, j, k) * ovol
5508 qy(i, j, k) =
qy(i, j, k) * ovol
5509 qz(i, j, k) =
qz(i, j, k) * ovol
5531 logical,
intent(in),
optional :: storeWallTensor
5534 real(kind=realtype) :: rfilv, por, mul, mue, mut, heatcoef
5535 real(kind=realtype) :: gm1, factlamheat, factturbheat
5536 real(kind=realtype) :: u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z
5537 real(kind=realtype) :: q_x, q_y, q_z
5538 real(kind=realtype) :: corr, ssx, ssy, ssz, fracdiv, snrm
5539 real(kind=realtype) :: tauxx, tauyy, tauzz
5540 real(kind=realtype) :: tauxy, tauxz, tauyz
5541 real(kind=realtype) :: tauxxs, tauyys, tauzzs
5542 real(kind=realtype) :: tauxys, tauxzs, tauyzs
5543 real(kind=realtype) :: ubar, vbar, wbar
5544 real(kind=realtype) :: exx, eyy, ezz
5545 real(kind=realtype) :: exy, exz, eyz
5546 real(kind=realtype) :: wxx, wyy, wzz
5547 real(kind=realtype) :: wxy, wxz, wyz, wyx, wzx, wzy
5548 real(kind=realtype) :: den, ccr1
5549 real(kind=realtype) :: fmx, fmy, fmz, frhoe, fact
5550 integer(kind=intType) :: i, j, k, io, jo, ko
5551 real(kind=realtype),
parameter :: xminn = 1.e-10_realtype
5552 real(kind=realtype),
parameter :: twothird =
two *
third
5553 real(kind=realtype),
dimension(9, 2:max(il, jl), 2:max(jl, kl), 2) :: tmpstore
5555 logical :: storeWall
5558 if (
present(storewalltensor))
then
5559 storewall = storewalltensor
5590 mul = por * (
rlv(i, j, k) +
rlv(i, j, k + 1))
5591 mue = por * (
rev(i, j, k) +
rev(i, j, k + 1))
5598 heatcoef = mul * factlamheat + mue * factturbheat
5603 u_x =
fourth * (
ux(i - 1, j - 1, k) +
ux(i, j - 1, k) &
5604 +
ux(i - 1, j, k) +
ux(i, j, k))
5605 u_y =
fourth * (
uy(i - 1, j - 1, k) +
uy(i, j - 1, k) &
5606 +
uy(i - 1, j, k) +
uy(i, j, k))
5607 u_z =
fourth * (
uz(i - 1, j - 1, k) +
uz(i, j - 1, k) &
5608 +
uz(i - 1, j, k) +
uz(i, j, k))
5610 v_x =
fourth * (
vx(i - 1, j - 1, k) +
vx(i, j - 1, k) &
5611 +
vx(i - 1, j, k) +
vx(i, j, k))
5612 v_y =
fourth * (
vy(i - 1, j - 1, k) +
vy(i, j - 1, k) &
5613 +
vy(i - 1, j, k) +
vy(i, j, k))
5614 v_z =
fourth * (
vz(i - 1, j - 1, k) +
vz(i, j - 1, k) &
5615 +
vz(i - 1, j, k) +
vz(i, j, k))
5617 w_x =
fourth * (
wx(i - 1, j - 1, k) +
wx(i, j - 1, k) &
5618 +
wx(i - 1, j, k) +
wx(i, j, k))
5619 w_y =
fourth * (
wy(i - 1, j - 1, k) +
wy(i, j - 1, k) &
5620 +
wy(i - 1, j, k) +
wy(i, j, k))
5621 w_z =
fourth * (
wz(i - 1, j - 1, k) +
wz(i, j - 1, k) &
5622 +
wz(i - 1, j, k) +
wz(i, j, k))
5624 q_x =
fourth * (
qx(i - 1, j - 1, k) +
qx(i, j - 1, k) &
5625 +
qx(i - 1, j, k) +
qx(i, j, k))
5626 q_y =
fourth * (
qy(i - 1, j - 1, k) +
qy(i, j - 1, k) &
5627 +
qy(i - 1, j, k) +
qy(i, j, k))
5628 q_z =
fourth * (
qz(i - 1, j - 1, k) +
qz(i, j - 1, k) &
5629 +
qz(i - 1, j, k) +
qz(i, j, k))
5636 ssx =
eighth * (
x(i - 1, j - 1, k + 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
5637 +
x(i - 1, j, k + 1, 1) -
x(i - 1, j, k - 1, 1) &
5638 +
x(i, j - 1, k + 1, 1) -
x(i, j - 1, k - 1, 1) &
5639 +
x(i, j, k + 1, 1) -
x(i, j, k - 1, 1))
5640 ssy =
eighth * (
x(i - 1, j - 1, k + 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
5641 +
x(i - 1, j, k + 1, 2) -
x(i - 1, j, k - 1, 2) &
5642 +
x(i, j - 1, k + 1, 2) -
x(i, j - 1, k - 1, 2) &
5643 +
x(i, j, k + 1, 2) -
x(i, j, k - 1, 2))
5644 ssz =
eighth * (
x(i - 1, j - 1, k + 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
5645 +
x(i - 1, j, k + 1, 3) -
x(i - 1, j, k - 1, 3) &
5646 +
x(i, j - 1, k + 1, 3) -
x(i, j - 1, k - 1, 3) &
5647 +
x(i, j, k + 1, 3) -
x(i, j, k - 1, 3))
5652 snrm =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
5659 corr = u_x * ssx + u_y * ssy + u_z * ssz &
5660 - (
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)) * snrm
5661 u_x = u_x - corr * ssx
5662 u_y = u_y - corr * ssy
5663 u_z = u_z - corr * ssz
5665 corr = v_x * ssx + v_y * ssy + v_z * ssz &
5666 - (
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)) * snrm
5667 v_x = v_x - corr * ssx
5668 v_y = v_y - corr * ssy
5669 v_z = v_z - corr * ssz
5671 corr = w_x * ssx + w_y * ssy + w_z * ssz &
5672 - (
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)) * snrm
5673 w_x = w_x - corr * ssx
5674 w_y = w_y - corr * ssy
5675 w_z = w_z - corr * ssz
5677 corr = q_x * ssx + q_y * ssy + q_z * ssz &
5678 + (
aa(i, j, k + 1) -
aa(i, j, k)) * snrm
5679 q_x = q_x - corr * ssx
5680 q_y = q_y - corr * ssy
5681 q_z = q_z - corr * ssz
5690 fracdiv = twothird * (u_x + v_y + w_z)
5692 tauxxs =
two * u_x - fracdiv
5693 tauyys =
two * v_y - fracdiv
5694 tauzzs =
two * w_z - fracdiv
5700 q_x = heatcoef * q_x
5701 q_y = heatcoef * q_y
5702 q_z = heatcoef * q_z
5722 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
5723 v_x * v_x + v_y * v_y + v_z * v_z + &
5724 w_x * w_x + w_y * w_y + w_z * w_z)
5728 den = max(den, xminn)
5733 fact = mue * ccr1 / den
5745 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
5746 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
5747 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
5749 exy = fact * (wxy * tauyys + wxz * tauyzs + &
5750 wyx * tauxxs + wyz * tauxzs)
5751 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
5752 wzx * tauxxs + wzy * tauxys)
5753 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
5754 wzx * tauxys + wzy * tauyys)
5757 tauxx = mut * tauxxs - exx
5758 tauyy = mut * tauyys - eyy
5759 tauzz = mut * tauzzs - ezz
5760 tauxy = mut * tauxys - exy
5761 tauxz = mut * tauxzs - exz
5762 tauyz = mut * tauyzs - eyz
5767 tauxx = mut * tauxxs
5768 tauyy = mut * tauyys
5769 tauzz = mut * tauzzs
5770 tauxy = mut * tauxys
5771 tauxz = mut * tauxzs
5772 tauyz = mut * tauyzs
5785 fmx = tauxx *
sk(i, j, k, 1) + tauxy *
sk(i, j, k, 2) &
5786 + tauxz *
sk(i, j, k, 3)
5787 fmy = tauxy *
sk(i, j, k, 1) + tauyy *
sk(i, j, k, 2) &
5788 + tauyz *
sk(i, j, k, 3)
5789 fmz = tauxz *
sk(i, j, k, 1) + tauyz *
sk(i, j, k, 2) &
5790 + tauzz *
sk(i, j, k, 3)
5791 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sk(i, j, k, 1)
5792 frhoe = frhoe + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sk(i, j, k, 2)
5793 frhoe = frhoe + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sk(i, j, k, 3)
5794 frhoe = frhoe - q_x *
sk(i, j, k, 1) - q_y *
sk(i, j, k, 2) - q_z *
sk(i, j, k, 3)
5803 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fmx
5804 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fmy
5805 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fmz
5812 tmpstore(1, i, j, 1) = tauxx
5813 tmpstore(2, i, j, 1) = tauyy
5814 tmpstore(3, i, j, 1) = tauzz
5815 tmpstore(4, i, j, 1) = tauxy
5816 tmpstore(5, i, j, 1) = tauxz
5817 tmpstore(6, i, j, 1) = tauyz
5819 tmpstore(7, i, j, 1) = q_x
5820 tmpstore(8, i, j, 1) = q_y
5821 tmpstore(9, i, j, 1) = q_z
5825 tmpstore(1, i, j, 2) = tauxx
5826 tmpstore(2, i, j, 2) = tauyy
5827 tmpstore(3, i, j, 2) = tauzz
5828 tmpstore(4, i, j, 2) = tauxy
5829 tmpstore(5, i, j, 2) = tauxz
5830 tmpstore(6, i, j, 2) = tauyz
5832 tmpstore(7, i, j, 2) = q_x
5833 tmpstore(8, i, j, 2) = q_y
5834 tmpstore(9, i, j, 2) = q_z
5844 origkmin:
if (
kk - 1 == 1)
then
5858 origkmax:
if (
kk +
nz - 1 == bkl)
then
5890 mul = por * (
rlv(i, j, k) +
rlv(i, j + 1, k))
5891 mue = por * (
rev(i, j, k) +
rev(i, j + 1, k))
5898 heatcoef = mul * factlamheat + mue * factturbheat
5903 u_x =
fourth * (
ux(i - 1, j, k - 1) +
ux(i, j, k - 1) &
5904 +
ux(i - 1, j, k) +
ux(i, j, k))
5905 u_y =
fourth * (
uy(i - 1, j, k - 1) +
uy(i, j, k - 1) &
5906 +
uy(i - 1, j, k) +
uy(i, j, k))
5907 u_z =
fourth * (
uz(i - 1, j, k - 1) +
uz(i, j, k - 1) &
5908 +
uz(i - 1, j, k) +
uz(i, j, k))
5910 v_x =
fourth * (
vx(i - 1, j, k - 1) +
vx(i, j, k - 1) &
5911 +
vx(i - 1, j, k) +
vx(i, j, k))
5912 v_y =
fourth * (
vy(i - 1, j, k - 1) +
vy(i, j, k - 1) &
5913 +
vy(i - 1, j, k) +
vy(i, j, k))
5914 v_z =
fourth * (
vz(i - 1, j, k - 1) +
vz(i, j, k - 1) &
5915 +
vz(i - 1, j, k) +
vz(i, j, k))
5917 w_x =
fourth * (
wx(i - 1, j, k - 1) +
wx(i, j, k - 1) &
5918 +
wx(i - 1, j, k) +
wx(i, j, k))
5919 w_y =
fourth * (
wy(i - 1, j, k - 1) +
wy(i, j, k - 1) &
5920 +
wy(i - 1, j, k) +
wy(i, j, k))
5921 w_z =
fourth * (
wz(i - 1, j, k - 1) +
wz(i, j, k - 1) &
5922 +
wz(i - 1, j, k) +
wz(i, j, k))
5924 q_x =
fourth * (
qx(i - 1, j, k - 1) +
qx(i, j, k - 1) &
5925 +
qx(i - 1, j, k) +
qx(i, j, k))
5926 q_y =
fourth * (
qy(i - 1, j, k - 1) +
qy(i, j, k - 1) &
5927 +
qy(i - 1, j, k) +
qy(i, j, k))
5928 q_z =
fourth * (
qz(i - 1, j, k - 1) +
qz(i, j, k - 1) &
5929 +
qz(i - 1, j, k) +
qz(i, j, k))
5936 ssx =
eighth * (
x(i - 1, j + 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
5937 +
x(i - 1, j + 1, k, 1) -
x(i - 1, j - 1, k, 1) &
5938 +
x(i, j + 1, k - 1, 1) -
x(i, j - 1, k - 1, 1) &
5939 +
x(i, j + 1, k, 1) -
x(i, j - 1, k, 1))
5940 ssy =
eighth * (
x(i - 1, j + 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
5941 +
x(i - 1, j + 1, k, 2) -
x(i - 1, j - 1, k, 2) &
5942 +
x(i, j + 1, k - 1, 2) -
x(i, j - 1, k - 1, 2) &
5943 +
x(i, j + 1, k, 2) -
x(i, j - 1, k, 2))
5944 ssz =
eighth * (
x(i - 1, j + 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
5945 +
x(i - 1, j + 1, k, 3) -
x(i - 1, j - 1, k, 3) &
5946 +
x(i, j + 1, k - 1, 3) -
x(i, j - 1, k - 1, 3) &
5947 +
x(i, j + 1, k, 3) -
x(i, j - 1, k, 3))
5952 snrm =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
5959 corr = u_x * ssx + u_y * ssy + u_z * ssz &
5960 - (
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)) * snrm
5961 u_x = u_x - corr * ssx
5962 u_y = u_y - corr * ssy
5963 u_z = u_z - corr * ssz
5965 corr = v_x * ssx + v_y * ssy + v_z * ssz &
5966 - (
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)) * snrm
5967 v_x = v_x - corr * ssx
5968 v_y = v_y - corr * ssy
5969 v_z = v_z - corr * ssz
5971 corr = w_x * ssx + w_y * ssy + w_z * ssz &
5972 - (
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)) * snrm
5973 w_x = w_x - corr * ssx
5974 w_y = w_y - corr * ssy
5975 w_z = w_z - corr * ssz
5977 corr = q_x * ssx + q_y * ssy + q_z * ssz &
5978 + (
aa(i, j + 1, k) -
aa(i, j, k)) * snrm
5979 q_x = q_x - corr * ssx
5980 q_y = q_y - corr * ssy
5981 q_z = q_z - corr * ssz
5990 fracdiv = twothird * (u_x + v_y + w_z)
5992 tauxxs =
two * u_x - fracdiv
5993 tauyys =
two * v_y - fracdiv
5994 tauzzs =
two * w_z - fracdiv
6000 q_x = heatcoef * q_x
6001 q_y = heatcoef * q_y
6002 q_z = heatcoef * q_z
6022 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
6023 v_x * v_x + v_y * v_y + v_z * v_z + &
6024 w_x * w_x + w_y * w_y + w_z * w_z)
6028 den = max(den, xminn)
6033 fact = mue * ccr1 / den
6045 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
6046 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
6047 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
6049 exy = fact * (wxy * tauyys + wxz * tauyzs + &
6050 wyx * tauxxs + wyz * tauxzs)
6051 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
6052 wzx * tauxxs + wzy * tauxys)
6053 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
6054 wzx * tauxys + wzy * tauyys)
6057 tauxx = mut * tauxxs - exx
6058 tauyy = mut * tauyys - eyy
6059 tauzz = mut * tauzzs - ezz
6060 tauxy = mut * tauxys - exy
6061 tauxz = mut * tauxzs - exz
6062 tauyz = mut * tauyzs - eyz
6067 tauxx = mut * tauxxs
6068 tauyy = mut * tauyys
6069 tauzz = mut * tauzzs
6070 tauxy = mut * tauxys
6071 tauxz = mut * tauxzs
6072 tauyz = mut * tauyzs
6085 fmx = tauxx *
sj(i, j, k, 1) + tauxy *
sj(i, j, k, 2) &
6086 + tauxz *
sj(i, j, k, 3)
6087 fmy = tauxy *
sj(i, j, k, 1) + tauyy *
sj(i, j, k, 2) &
6088 + tauyz *
sj(i, j, k, 3)
6089 fmz = tauxz *
sj(i, j, k, 1) + tauyz *
sj(i, j, k, 2) &
6090 + tauzz *
sj(i, j, k, 3)
6091 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sj(i, j, k, 1) &
6092 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sj(i, j, k, 2) &
6093 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sj(i, j, k, 3) &
6094 - q_x *
sj(i, j, k, 1) - q_y *
sj(i, j, k, 2) - q_z *
sj(i, j, k, 3)
6103 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fmx
6104 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fmy
6105 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fmz
6112 tmpstore(1, i, k, 1) = tauxx
6113 tmpstore(2, i, k, 1) = tauyy
6114 tmpstore(3, i, k, 1) = tauzz
6115 tmpstore(4, i, k, 1) = tauxy
6116 tmpstore(5, i, k, 1) = tauxz
6117 tmpstore(6, i, k, 1) = tauyz
6119 tmpstore(7, i, k, 1) = q_x
6120 tmpstore(8, i, k, 1) = q_y
6121 tmpstore(9, i, k, 1) = q_z
6125 tmpstore(1, i, k, 2) = tauxx
6126 tmpstore(2, i, k, 2) = tauyy
6127 tmpstore(3, i, k, 2) = tauzz
6128 tmpstore(4, i, k, 2) = tauxy
6129 tmpstore(5, i, k, 2) = tauxz
6130 tmpstore(6, i, k, 2) = tauyz
6132 tmpstore(7, i, k, 2) = q_x
6133 tmpstore(8, i, k, 2) = q_y
6134 tmpstore(9, i, k, 2) = q_z
6141 origjmin:
if (
jj - 1 == 1)
then
6155 origjmax:
if (
jj +
ny - 1 == bjl)
then
6186 mul = por * (
rlv(i, j, k) +
rlv(i + 1, j, k))
6187 mue = por * (
rev(i, j, k) +
rev(i + 1, j, k))
6194 heatcoef = mul * factlamheat + mue * factturbheat
6199 u_x =
fourth * (
ux(i, j - 1, k - 1) +
ux(i, j, k - 1) &
6200 +
ux(i, j - 1, k) +
ux(i, j, k))
6201 u_y =
fourth * (
uy(i, j - 1, k - 1) +
uy(i, j, k - 1) &
6202 +
uy(i, j - 1, k) +
uy(i, j, k))
6203 u_z =
fourth * (
uz(i, j - 1, k - 1) +
uz(i, j, k - 1) &
6204 +
uz(i, j - 1, k) +
uz(i, j, k))
6206 v_x =
fourth * (
vx(i, j - 1, k - 1) +
vx(i, j, k - 1) &
6207 +
vx(i, j - 1, k) +
vx(i, j, k))
6208 v_y =
fourth * (
vy(i, j - 1, k - 1) +
vy(i, j, k - 1) &
6209 +
vy(i, j - 1, k) +
vy(i, j, k))
6210 v_z =
fourth * (
vz(i, j - 1, k - 1) +
vz(i, j, k - 1) &
6211 +
vz(i, j - 1, k) +
vz(i, j, k))
6213 w_x =
fourth * (
wx(i, j - 1, k - 1) +
wx(i, j, k - 1) &
6214 +
wx(i, j - 1, k) +
wx(i, j, k))
6215 w_y =
fourth * (
wy(i, j - 1, k - 1) +
wy(i, j, k - 1) &
6216 +
wy(i, j - 1, k) +
wy(i, j, k))
6217 w_z =
fourth * (
wz(i, j - 1, k - 1) +
wz(i, j, k - 1) &
6218 +
wz(i, j - 1, k) +
wz(i, j, k))
6220 q_x =
fourth * (
qx(i, j - 1, k - 1) +
qx(i, j, k - 1) &
6221 +
qx(i, j - 1, k) +
qx(i, j, k))
6222 q_y =
fourth * (
qy(i, j - 1, k - 1) +
qy(i, j, k - 1) &
6223 +
qy(i, j - 1, k) +
qy(i, j, k))
6224 q_z =
fourth * (
qz(i, j - 1, k - 1) +
qz(i, j, k - 1) &
6225 +
qz(i, j - 1, k) +
qz(i, j, k))
6232 ssx =
eighth * (
x(i + 1, j - 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6233 +
x(i + 1, j - 1, k, 1) -
x(i - 1, j - 1, k, 1) &
6234 +
x(i + 1, j, k - 1, 1) -
x(i - 1, j, k - 1, 1) &
6235 +
x(i + 1, j, k, 1) -
x(i - 1, j, k, 1))
6236 ssy =
eighth * (
x(i + 1, j - 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6237 +
x(i + 1, j - 1, k, 2) -
x(i - 1, j - 1, k, 2) &
6238 +
x(i + 1, j, k - 1, 2) -
x(i - 1, j, k - 1, 2) &
6239 +
x(i + 1, j, k, 2) -
x(i - 1, j, k, 2))
6240 ssz =
eighth * (
x(i + 1, j - 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6241 +
x(i + 1, j - 1, k, 3) -
x(i - 1, j - 1, k, 3) &
6242 +
x(i + 1, j, k - 1, 3) -
x(i - 1, j, k - 1, 3) &
6243 +
x(i + 1, j, k, 3) -
x(i - 1, j, k, 3))
6248 snrm =
one / sqrt(ssx * ssx + ssy * ssy + ssz * ssz)
6255 corr = u_x * ssx + u_y * ssy + u_z * ssz &
6256 - (
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)) * snrm
6257 u_x = u_x - corr * ssx
6258 u_y = u_y - corr * ssy
6259 u_z = u_z - corr * ssz
6261 corr = v_x * ssx + v_y * ssy + v_z * ssz &
6262 - (
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)) * snrm
6263 v_x = v_x - corr * ssx
6264 v_y = v_y - corr * ssy
6265 v_z = v_z - corr * ssz
6267 corr = w_x * ssx + w_y * ssy + w_z * ssz &
6268 - (
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)) * snrm
6269 w_x = w_x - corr * ssx
6270 w_y = w_y - corr * ssy
6271 w_z = w_z - corr * ssz
6273 corr = q_x * ssx + q_y * ssy + q_z * ssz &
6274 + (
aa(i + 1, j, k) -
aa(i, j, k)) * snrm
6275 q_x = q_x - corr * ssx
6276 q_y = q_y - corr * ssy
6277 q_z = q_z - corr * ssz
6286 fracdiv = twothird * (u_x + v_y + w_z)
6288 tauxxs =
two * u_x - fracdiv
6289 tauyys =
two * v_y - fracdiv
6290 tauzzs =
two * w_z - fracdiv
6296 q_x = heatcoef * q_x
6297 q_y = heatcoef * q_y
6298 q_z = heatcoef * q_z
6318 den = sqrt(u_x * u_x + u_y * u_y + u_z * u_z + &
6319 v_x * v_x + v_y * v_y + v_z * v_z + &
6320 w_x * w_x + w_y * w_y + w_z * w_z)
6324 den = max(den, xminn)
6329 fact = mue * ccr1 / den
6341 exx = fact * (wxy * tauxys + wxz * tauxzs) *
two
6342 eyy = fact * (wyx * tauxys + wyz * tauyzs) *
two
6343 ezz = fact * (wzx * tauxzs + wzy * tauyzs) *
two
6345 exy = fact * (wxy * tauyys + wxz * tauyzs + &
6346 wyx * tauxxs + wyz * tauxzs)
6347 exz = fact * (wxy * tauyzs + wxz * tauzzs + &
6348 wzx * tauxxs + wzy * tauxys)
6349 eyz = fact * (wyx * tauxzs + wyz * tauzzs + &
6350 wzx * tauxys + wzy * tauyys)
6353 tauxx = mut * tauxxs - exx
6354 tauyy = mut * tauyys - eyy
6355 tauzz = mut * tauzzs - ezz
6356 tauxy = mut * tauxys - exy
6357 tauxz = mut * tauxzs - exz
6358 tauyz = mut * tauyzs - eyz
6363 tauxx = mut * tauxxs
6364 tauyy = mut * tauyys
6365 tauzz = mut * tauzzs
6366 tauxy = mut * tauxys
6367 tauxz = mut * tauxzs
6368 tauyz = mut * tauyzs
6381 fmx = tauxx *
si(i, j, k, 1) + tauxy *
si(i, j, k, 2) &
6382 + tauxz *
si(i, j, k, 3)
6383 fmy = tauxy *
si(i, j, k, 1) + tauyy *
si(i, j, k, 2) &
6384 + tauyz *
si(i, j, k, 3)
6385 fmz = tauxz *
si(i, j, k, 1) + tauyz *
si(i, j, k, 2) &
6386 + tauzz *
si(i, j, k, 3)
6387 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
si(i, j, k, 1) &
6388 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
si(i, j, k, 2) &
6389 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
si(i, j, k, 3) &
6390 - q_x *
si(i, j, k, 1) - q_y *
si(i, j, k, 2) - q_z *
si(i, j, k, 3)
6393 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fmx
6394 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fmy
6395 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fmz
6407 tmpstore(1, j, k, 1) = tauxx
6408 tmpstore(2, j, k, 1) = tauyy
6409 tmpstore(3, j, k, 1) = tauzz
6410 tmpstore(4, j, k, 1) = tauxy
6411 tmpstore(5, j, k, 1) = tauxz
6412 tmpstore(6, j, k, 1) = tauyz
6414 tmpstore(7, j, k, 1) = q_x
6415 tmpstore(8, j, k, 1) = q_y
6416 tmpstore(9, j, k, 1) = q_z
6420 tmpstore(1, j, k, 2) = tauxx
6421 tmpstore(2, j, k, 2) = tauyy
6422 tmpstore(3, j, k, 2) = tauzz
6423 tmpstore(4, j, k, 2) = tauxy
6424 tmpstore(5, j, k, 2) = tauxz
6425 tmpstore(6, j, k, 2) = tauyz
6427 tmpstore(7, j, k, 2) = q_x
6428 tmpstore(8, j, k, 2) = q_y
6429 tmpstore(9, j, k, 2) = q_z
6436 origimin:
if (
ii - 1 == 1)
then
6450 origimax:
if (
ii +
nx - 1 == bil)
then
6475 real(kind=realtype),
parameter :: twothird =
two *
third
6479 integer(kind=intType) :: i, j, k
6480 integer(kind=intType) :: ii, jj, kk
6482 real(kind=realtype) :: rfilv, por, mul, mue, mut, heatcoef
6483 real(kind=realtype) :: gm1, factlamheat, factturbheat
6484 real(kind=realtype) :: u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z
6485 real(kind=realtype) :: q_x, q_y, q_z, ubar, vbar, wbar
6486 real(kind=realtype) :: corr, ssx, ssy, ssz,
ss, fracdiv
6487 real(kind=realtype) :: tauxx, tauyy, tauzz
6488 real(kind=realtype) :: tauxy, tauxz, tauyz
6489 real(kind=realtype) :: fmx, fmy, fmz, frhoe
6490 real(kind=realtype) :: dd
6491 logical :: correctForK
6503 ssx =
eighth * (
x(i + 1, j - 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6504 +
x(i + 1, j - 1, k, 1) -
x(i - 1, j - 1, k, 1) &
6505 +
x(i + 1, j, k - 1, 1) -
x(i - 1, j, k - 1, 1) &
6506 +
x(i + 1, j, k, 1) -
x(i - 1, j, k, 1))
6507 ssy =
eighth * (
x(i + 1, j - 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6508 +
x(i + 1, j - 1, k, 2) -
x(i - 1, j - 1, k, 2) &
6509 +
x(i + 1, j, k - 1, 2) -
x(i - 1, j, k - 1, 2) &
6510 +
x(i + 1, j, k, 2) -
x(i - 1, j, k, 2))
6511 ssz =
eighth * (
x(i + 1, j - 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6512 +
x(i + 1, j - 1, k, 3) -
x(i - 1, j - 1, k, 3) &
6513 +
x(i + 1, j, k - 1, 3) -
x(i - 1, j, k - 1, 3) &
6514 +
x(i + 1, j, k, 3) -
x(i - 1, j, k, 3))
6517 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
6523 dd =
w(i + 1, j, k,
ivx) -
w(i, j, k,
ivx)
6528 dd =
w(i + 1, j, k,
ivy) -
w(i, j, k,
ivy)
6533 dd =
w(i + 1, j, k,
ivz) -
w(i, j, k,
ivz)
6538 dd =
aa(i + 1, j, k) -
aa(i, j, k)
6551 mul = por * (
rlv(i, j, k) +
rlv(i + 1, j, k))
6552 mue = por * (
rev(i, j, k) +
rev(i + 1, j, k))
6559 heatcoef = mul * factlamheat + mue * factturbheat
6563 fracdiv = twothird * (u_x + v_y + w_z)
6565 tauxx = mut * (
two * u_x - fracdiv)
6566 tauyy = mut * (
two * v_y - fracdiv)
6567 tauzz = mut * (
two * w_z - fracdiv)
6569 tauxy = mut * (u_y + v_x)
6570 tauxz = mut * (u_z + w_x)
6571 tauyz = mut * (v_z + w_y)
6573 q_x = heatcoef * q_x
6574 q_y = heatcoef * q_y
6575 q_z = heatcoef * q_z
6586 fmx = tauxx *
si(i, j, k, 1) + tauxy *
si(i, j, k, 2) + tauxz *
si(i, j, k, 3)
6587 fmy = tauxy *
si(i, j, k, 1) + tauyy *
si(i, j, k, 2) + tauyz *
si(i, j, k, 3)
6588 fmz = tauxz *
si(i, j, k, 1) + tauyz *
si(i, j, k, 2) + tauzz *
si(i, j, k, 3)
6589 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
si(i, j, k, 1) &
6590 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
si(i, j, k, 2) &
6591 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
si(i, j, k, 3) &
6592 - q_x *
si(i, j, k, 1) - q_y *
si(i, j, k, 2) - q_z *
si(i, j, k, 3)
6595 fw(i + 1, j, k,
imx) =
fw(i + 1, j, k,
imx) + fmx
6596 fw(i + 1, j, k,
imy) =
fw(i + 1, j, k,
imy) + fmy
6597 fw(i + 1, j, k,
imz) =
fw(i + 1, j, k,
imz) + fmz
6616 ssx =
eighth * (
x(i - 1, j + 1, k - 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6617 +
x(i - 1, j + 1, k, 1) -
x(i - 1, j - 1, k, 1) &
6618 +
x(i, j + 1, k - 1, 1) -
x(i, j - 1, k - 1, 1) &
6619 +
x(i, j + 1, k, 1) -
x(i, j - 1, k, 1))
6620 ssy =
eighth * (
x(i - 1, j + 1, k - 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6621 +
x(i - 1, j + 1, k, 2) -
x(i - 1, j - 1, k, 2) &
6622 +
x(i, j + 1, k - 1, 2) -
x(i, j - 1, k - 1, 2) &
6623 +
x(i, j + 1, k, 2) -
x(i, j - 1, k, 2))
6624 ssz =
eighth * (
x(i - 1, j + 1, k - 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6625 +
x(i - 1, j + 1, k, 3) -
x(i - 1, j - 1, k, 3) &
6626 +
x(i, j + 1, k - 1, 3) -
x(i, j - 1, k - 1, 3) &
6627 +
x(i, j + 1, k, 3) -
x(i, j - 1, k, 3))
6630 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
6636 dd =
w(i, j + 1, k,
ivx) -
w(i, j, k,
ivx)
6641 dd =
w(i, j + 1, k,
ivy) -
w(i, j, k,
ivy)
6646 dd =
w(i, j + 1, k,
ivz) -
w(i, j, k,
ivz)
6651 dd =
aa(i, j + 1, k) -
aa(i, j, k)
6664 mul = por * (
rlv(i, j, k) +
rlv(i, j + 1, k))
6665 mue = por * (
rev(i, j, k) +
rev(i, j + 1, k))
6672 heatcoef = mul * factlamheat + mue * factturbheat
6676 fracdiv = twothird * (u_x + v_y + w_z)
6678 tauxx = mut * (
two * u_x - fracdiv)
6679 tauyy = mut * (
two * v_y - fracdiv)
6680 tauzz = mut * (
two * w_z - fracdiv)
6682 tauxy = mut * (u_y + v_x)
6683 tauxz = mut * (u_z + w_x)
6684 tauyz = mut * (v_z + w_y)
6686 q_x = heatcoef * q_x
6687 q_y = heatcoef * q_y
6688 q_z = heatcoef * q_z
6699 fmx = tauxx *
sj(i, j, k, 1) + tauxy *
sj(i, j, k, 2) + tauxz *
sj(i, j, k, 3)
6700 fmy = tauxy *
sj(i, j, k, 1) + tauyy *
sj(i, j, k, 2) + tauyz *
sj(i, j, k, 3)
6701 fmz = tauxz *
sj(i, j, k, 1) + tauyz *
sj(i, j, k, 2) + tauzz *
sj(i, j, k, 3)
6702 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sj(i, j, k, 1) &
6703 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sj(i, j, k, 2) &
6704 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sj(i, j, k, 3) &
6705 - q_x *
sj(i, j, k, 1) - q_y *
sj(i, j, k, 2) - q_z *
sj(i, j, k, 3)
6714 fw(i, j + 1, k,
imx) =
fw(i, j + 1, k,
imx) + fmx
6715 fw(i, j + 1, k,
imy) =
fw(i, j + 1, k,
imy) + fmy
6716 fw(i, j + 1, k,
imz) =
fw(i, j + 1, k,
imz) + fmz
6730 ssx =
eighth * (
x(i - 1, j - 1, k + 1, 1) -
x(i - 1, j - 1, k - 1, 1) &
6731 +
x(i - 1, j, k + 1, 1) -
x(i - 1, j, k - 1, 1) &
6732 +
x(i, j - 1, k + 1, 1) -
x(i, j - 1, k - 1, 1) &
6733 +
x(i, j, k + 1, 1) -
x(i, j, k - 1, 1))
6734 ssy =
eighth * (
x(i - 1, j - 1, k + 1, 2) -
x(i - 1, j - 1, k - 1, 2) &
6735 +
x(i - 1, j, k + 1, 2) -
x(i - 1, j, k - 1, 2) &
6736 +
x(i, j - 1, k + 1, 2) -
x(i, j - 1, k - 1, 2) &
6737 +
x(i, j, k + 1, 2) -
x(i, j, k - 1, 2))
6738 ssz =
eighth * (
x(i - 1, j - 1, k + 1, 3) -
x(i - 1, j - 1, k - 1, 3) &
6739 +
x(i - 1, j, k + 1, 3) -
x(i - 1, j, k - 1, 3) &
6740 +
x(i, j - 1, k + 1, 3) -
x(i, j - 1, k - 1, 3) &
6741 +
x(i, j, k + 1, 3) -
x(i, j, k - 1, 3))
6743 ss =
one / (ssx * ssx + ssy * ssy + ssz * ssz)
6749 dd =
w(i, j, k + 1,
ivx) -
w(i, j, k,
ivx)
6754 dd =
w(i, j, k + 1,
ivy) -
w(i, j, k,
ivy)
6759 dd =
w(i, j, k + 1,
ivz) -
w(i, j, k,
ivz)
6764 dd =
aa(i, j, k + 1) -
aa(i, j, k)
6777 mul = por * (
rlv(i, j, k) +
rlv(i, j, k + 1))
6778 mue = por * (
rev(i, j, k) +
rev(i, j, k + 1))
6785 heatcoef = mul * factlamheat + mue * factturbheat
6789 fracdiv = twothird * (u_x + v_y + w_z)
6791 tauxx = mut * (
two * u_x - fracdiv)
6792 tauyy = mut * (
two * v_y - fracdiv)
6793 tauzz = mut * (
two * w_z - fracdiv)
6795 tauxy = mut * (u_y + v_x)
6796 tauxz = mut * (u_z + w_x)
6797 tauyz = mut * (v_z + w_y)
6799 q_x = heatcoef * q_x
6800 q_y = heatcoef * q_y
6801 q_z = heatcoef * q_z
6812 fmx = tauxx *
sk(i, j, k, 1) + tauxy *
sk(i, j, k, 2) + tauxz *
sk(i, j, k, 3)
6813 fmy = tauxy *
sk(i, j, k, 1) + tauyy *
sk(i, j, k, 2) + tauyz *
sk(i, j, k, 3)
6814 fmz = tauxz *
sk(i, j, k, 1) + tauyz *
sk(i, j, k, 2) + tauzz *
sk(i, j, k, 3)
6815 frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) *
sk(i, j, k, 1) &
6816 + (ubar * tauxy + vbar * tauyy + wbar * tauyz) *
sk(i, j, k, 2) &
6817 + (ubar * tauxz + vbar * tauyz + wbar * tauzz) *
sk(i, j, k, 3) &
6818 - q_x *
sk(i, j, k, 1) - q_y *
sk(i, j, k, 2) - q_z *
sk(i, j, k, 3)
6827 fw(i, j, k + 1,
imx) =
fw(i, j, k + 1,
imx) + fmx
6828 fw(i, j, k + 1,
imy) =
fw(i, j, k + 1,
imy) + fmy
6829 fw(i, j, k + 1,
imz) =
fw(i, j, k + 1,
imz) + fmz
6848 integer(kind=intType) :: nTurb, i, j, k, l
6849 real(kind=realtype) :: ovol, rblank
6856 rblank = max(real(
iblank(i, j, k), realtype),
zero)
6857 dw(i, j, k, l) = (
dw(i, j, k, l) +
fw(i, j, k, l)) * rblank
6872 integer(kind=intType) :: i, j, k, ii, nTurb
6873 real(kind=realtype) :: ovol
6882 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), parameter rsacw1
real(kind=realtype), parameter rsak
real(kind=realtype), parameter rsact4
real(kind=realtype), parameter rsacrot
real(kind=realtype), parameter rsacb2
real(kind=realtype), parameter rsact3
real(kind=realtype), parameter rsacw3
real(kind=realtype), parameter rsacw2
real(kind=realtype), parameter rsacv1
real(kind=realtype), parameter rsacb3
real(kind=realtype), parameter rsacb1
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)