364 lines
15 KiB
Fortran
364 lines
15 KiB
Fortran
|
|
! *** generated by SAPFOR with version 2415 and build date: May 4 2025 14:48:40
|
|
! *** Enabled options ***:
|
|
! *** maximum shadow width is 50 percent
|
|
! *** generated by SAPFOR
|
|
|
|
!---------------------------------------------------------------------
|
|
!---------------------------------------------------------------------
|
|
subroutine z_solve_sp ()
|
|
|
|
include 'header_sp.h'
|
|
integer :: i,j,k,k1,k2,m,m1
|
|
double precision :: ru1,fac1,fac2,rhs__(5,0:2)
|
|
double precision :: lhs__(5,0:2),lhsm__(5,0:2),lhsp__(5,0:2)
|
|
double precision :: t1,t2,t3,ac,xvel,yvel,zvel,btuz,ac2u,uzik1
|
|
|
|
!---------------------------------------------------------------------
|
|
!---------------------------------------------------------------------
|
|
!---------------------------------------------------------------------
|
|
! Prepare for z-solve, array redistribution
|
|
!---------------------------------------------------------------------
|
|
|
|
!$SPF PARALLEL_REG r0
|
|
!$SPF ANALYSIS(PRIVATE(LHS__,LHSP__,LHSM__,RHS__))
|
|
! DVM$ PARALLEL (J,I) ON U(*,I,J,*), CUDA_BLOCK (32,4),PRIVATE (M,K1,K2,R
|
|
! DVM$&U1,FAC1,FAC2,K,LHS__,LHSP__,LHSM__,RHS__,T1,T2,T3,AC,XVEL,YVEL,ZVE
|
|
! DVM$&L,BTUZ,AC2U,UZIK1)
|
|
! DVM$ REGION LOCAL (LHS)
|
|
do j = 1,ny2
|
|
do i = 1,nx2
|
|
lhs__(1,0) = 0.0d0
|
|
lhsp__(1,0) = 0.0d0
|
|
lhsm__(1,0) = 0.0d0
|
|
lhs__(2,0) = 0.0d0
|
|
lhsp__(2,0) = 0.0d0
|
|
lhsm__(2,0) = 0.0d0
|
|
lhs__(3,0) = 1.0d0
|
|
lhsp__(3,0) = 1.0d0
|
|
lhsm__(3,0) = 1.0d0
|
|
lhs__(4,0) = 0.0d0
|
|
lhsp__(4,0) = 0.0d0
|
|
lhsm__(4,0) = 0.0d0
|
|
lhs__(5,0) = 0.0d0
|
|
lhsp__(5,0) = 0.0d0
|
|
lhsm__(5,0) = 0.0d0
|
|
lhs__(1,1) = 0.0d0
|
|
ru1 = c3c4 * 1.0d0 / u(1,i,j,0)
|
|
ru1 = dmax1 (dz4 + con43 * ru1,dz5 + c1c5 * ru1,dzmax + ru1,
|
|
&dz1)
|
|
lhs__(2,1) = (-(dttz2)) * ws(i,j,0) - dttz1 * ru1
|
|
ru1 = c3c4 * 1.0d0 / u(1,i,j,1)
|
|
ru1 = dmax1 (dz4 + con43 * ru1,dz5 + c1c5 * ru1,dzmax + ru1,
|
|
&dz1)
|
|
lhs__(3,1) = 1.0d0 + c2dttz1 * ru1
|
|
ru1 = c3c4 * 1.0d0 / u(1,i,j,2)
|
|
ru1 = dmax1 (dz4 + con43 * ru1,dz5 + c1c5 * ru1,dzmax + ru1,
|
|
&dz1)
|
|
lhs__(4,1) = dttz2 * ws(i,j,2) - dttz1 * ru1
|
|
lhs__(5,1) = 0.0d0
|
|
lhs__(3,1) = lhs__(3,1) + comz5
|
|
lhs__(4,1) = lhs__(4,1) - comz4
|
|
lhs__(5,1) = lhs__(5,1) + comz1
|
|
lhsp__(1,1) = lhs__(1,1)
|
|
lhsp__(2,1) = lhs__(2,1) - dttz2 * speed(i,j,1 - 1)
|
|
lhsp__(3,1) = lhs__(3,1)
|
|
lhsp__(4,1) = lhs__(4,1) + dttz2 * speed(i,j,1 + 1)
|
|
lhsp__(5,1) = lhs__(5,1)
|
|
lhsm__(1,1) = lhs__(1,1)
|
|
lhsm__(2,1) = lhs__(2,1) + dttz2 * speed(i,j,1 - 1)
|
|
lhsm__(3,1) = lhs__(3,1)
|
|
lhsm__(4,1) = lhs__(4,1) - dttz2 * speed(i,j,1 + 1)
|
|
lhsm__(5,1) = lhs__(5,1)
|
|
do k = 0,nz2 + 1
|
|
if (k + 2 .lt. nz2 + 1) then
|
|
m = k + 2
|
|
lhs__(1,2) = 0.0d0
|
|
ru1 = c3c4 * 1.0d0 / u(1,i,j,m - 1)
|
|
ru1 = dmax1 (dz4 + con43 * ru1,dz5 + c1c5 * ru1,dzmax
|
|
&+ ru1,dz1)
|
|
lhs__(2,2) = (-(dttz2)) * ws(i,j,m - 1) - dttz1 * ru1
|
|
ru1 = c3c4 * 1.0d0 / u(1,i,j,m)
|
|
ru1 = dmax1 (dz4 + con43 * ru1,dz5 + c1c5 * ru1,dzmax
|
|
&+ ru1,dz1)
|
|
lhs__(3,2) = 1.0d0 + c2dttz1 * ru1
|
|
ru1 = c3c4 * 1.0d0 / u(1,i,j,m + 1)
|
|
ru1 = dmax1 (dz4 + con43 * ru1,dz5 + c1c5 * ru1,dzmax
|
|
&+ ru1,dz1)
|
|
lhs__(4,2) = dttz2 * ws(i,j,m + 1) - dttz1 * ru1
|
|
lhs__(5,2) = 0.0d0
|
|
if (m .eq. 1) then
|
|
lhs__(3,2) = lhs__(3,2) + comz5
|
|
lhs__(4,2) = lhs__(4,2) - comz4
|
|
lhs__(5,2) = lhs__(5,2) + comz1
|
|
else if (m .eq. 2) then
|
|
lhs__(2,2) = lhs__(2,2) - comz4
|
|
lhs__(3,2) = lhs__(3,2) + comz6
|
|
lhs__(4,2) = lhs__(4,2) - comz4
|
|
lhs__(5,2) = lhs__(5,2) + comz1
|
|
else if (m .ge. 3 .and. m .le. nz2 - 2) then
|
|
lhs__(1,2) = lhs__(1,2) + comz1
|
|
lhs__(2,2) = lhs__(2,2) - comz4
|
|
lhs__(3,2) = lhs__(3,2) + comz6
|
|
lhs__(4,2) = lhs__(4,2) - comz4
|
|
lhs__(5,2) = lhs__(5,2) + comz1
|
|
else if (m .eq. nz2 - 1) then
|
|
lhs__(1,2) = lhs__(1,2) + comz1
|
|
lhs__(2,2) = lhs__(2,2) - comz4
|
|
lhs__(3,2) = lhs__(3,2) + comz6
|
|
lhs__(4,2) = lhs__(4,2) - comz4
|
|
else if (m .eq. nz2) then
|
|
lhs__(1,2) = lhs__(1,2) + comz1
|
|
lhs__(2,2) = lhs__(2,2) - comz4
|
|
lhs__(3,2) = lhs__(3,2) + comz5
|
|
endif
|
|
lhsp__(1,2) = lhs__(1,2)
|
|
lhsp__(2,2) = lhs__(2,2) - dttz2 * speed(i,j,m - 1)
|
|
lhsp__(3,2) = lhs__(3,2)
|
|
lhsp__(4,2) = lhs__(4,2) + dttz2 * speed(i,j,m + 1)
|
|
lhsp__(5,2) = lhs__(5,2)
|
|
lhsm__(1,2) = lhs__(1,2)
|
|
lhsm__(2,2) = lhs__(2,2) + dttz2 * speed(i,j,m - 1)
|
|
lhsm__(3,2) = lhs__(3,2)
|
|
lhsm__(4,2) = lhs__(4,2) - dttz2 * speed(i,j,m + 1)
|
|
lhsm__(5,2) = lhs__(5,2)
|
|
else if (k + 2 .eq. nz2 + 1) then
|
|
lhs__(1,2) = 0.0d0
|
|
lhsp__(1,2) = 0.0d0
|
|
lhsm__(1,2) = 0.0d0
|
|
lhs__(2,2) = 0.0d0
|
|
lhsp__(2,2) = 0.0d0
|
|
lhsm__(2,2) = 0.0d0
|
|
lhs__(3,2) = 1.0d0
|
|
lhsp__(3,2) = 1.0d0
|
|
lhsm__(3,2) = 1.0d0
|
|
lhs__(4,2) = 0.0d0
|
|
lhsp__(4,2) = 0.0d0
|
|
lhsm__(4,2) = 0.0d0
|
|
lhs__(5,2) = 0.0d0
|
|
lhsp__(5,2) = 0.0d0
|
|
lhsm__(5,2) = 0.0d0
|
|
endif
|
|
|
|
!********************************** end of init
|
|
k1 = k + 1
|
|
k2 = k + 2
|
|
fac1 = 1.d0 / lhs__(3,0)
|
|
lhs__(4,0) = fac1 * lhs__(4,0)
|
|
lhs__(5,0) = fac1 * lhs__(5,0)
|
|
do m = 1,3
|
|
rhs(m,i,j,k) = fac1 * rhs(m,i,j,k)
|
|
enddo
|
|
if (k .le. nz2 - 1) then
|
|
lhs__(3,1) = lhs__(3,1) - lhs__(2,1) * lhs__(4,0)
|
|
lhs__(4,1) = lhs__(4,1) - lhs__(2,1) * lhs__(5,0)
|
|
lhs__(2,2) = lhs__(2,2) - lhs__(1,2) * lhs__(4,0)
|
|
lhs__(3,2) = lhs__(3,2) - lhs__(1,2) * lhs__(5,0)
|
|
do m = 1,3
|
|
rhs(m,i,j,k1) = rhs(m,i,j,k1) - lhs__(2,1) * rhs(m,
|
|
&i,j,k)
|
|
rhs(m,i,j,k2) = rhs(m,i,j,k2) - lhs__(1,2) * rhs(m,
|
|
&i,j,k)
|
|
enddo
|
|
else
|
|
lhs__(3,1) = lhs__(3,1) - lhs__(2,1) * lhs__(4,0)
|
|
lhs__(4,1) = lhs__(4,1) - lhs__(2,1) * lhs__(5,0)
|
|
fac2 = 1.d0 / lhs__(3,1)
|
|
do m = 1,3
|
|
rhs(m,i,j,k1) = rhs(m,i,j,k1) - lhs__(2,1) * rhs(m,
|
|
&i,j,k)
|
|
rhs(m,i,j,k1) = fac2 * rhs(m,i,j,k1)
|
|
enddo
|
|
endif
|
|
m = 4
|
|
fac1 = 1.d0 / lhsp__(3,0)
|
|
lhsp__(4,0) = fac1 * lhsp__(4,0)
|
|
lhsp__(5,0) = fac1 * lhsp__(5,0)
|
|
rhs(m,i,j,k) = fac1 * rhs(m,i,j,k)
|
|
lhsp__(3,1) = lhsp__(3,1) - lhsp__(2,1) * lhsp__(4,0)
|
|
lhsp__(4,1) = lhsp__(4,1) - lhsp__(2,1) * lhsp__(5,0)
|
|
rhs(m,i,j,k1) = rhs(m,i,j,k1) - lhsp__(2,1) * rhs(m,i,j,k
|
|
&)
|
|
if (k .lt. nz2) then
|
|
lhsp__(2,2) = lhsp__(2,2) - lhsp__(1,2) * lhsp__(4,0)
|
|
lhsp__(3,2) = lhsp__(3,2) - lhsp__(1,2) * lhsp__(5,0)
|
|
rhs(m,i,j,k2) = rhs(m,i,j,k2) - lhsp__(1,2) * rhs(m,i,
|
|
&j,k)
|
|
endif
|
|
m = 5
|
|
fac1 = 1.d0 / lhsm__(3,0)
|
|
lhsm__(4,0) = fac1 * lhsm__(4,0)
|
|
lhsm__(5,0) = fac1 * lhsm__(5,0)
|
|
rhs(m,i,j,k) = fac1 * rhs(m,i,j,k)
|
|
lhsm__(3,1) = lhsm__(3,1) - lhsm__(2,1) * lhsm__(4,0)
|
|
lhsm__(4,1) = lhsm__(4,1) - lhsm__(2,1) * lhsm__(5,0)
|
|
rhs(m,i,j,k1) = rhs(m,i,j,k1) - lhsm__(2,1) * rhs(m,i,j,k
|
|
&)
|
|
if (k .lt. nz2) then
|
|
lhsm__(2,2) = lhsm__(2,2) - lhsm__(1,2) * lhsm__(4,0)
|
|
lhsm__(3,2) = lhsm__(3,2) - lhsm__(1,2) * lhsm__(5,0)
|
|
rhs(m,i,j,k2) = rhs(m,i,j,k2) - lhsm__(1,2) * rhs(m,i,
|
|
&j,k)
|
|
endif
|
|
if (k .eq. nz2) then
|
|
rhs(4,i,j,k1) = rhs(4,i,j,k1) / lhsp__(3,1)
|
|
rhs(5,i,j,k1) = rhs(5,i,j,k1) / lhsm__(3,1)
|
|
do m = 1,3
|
|
rhs(m,i,j,k) = rhs(m,i,j,k) - lhs__(4,0) * rhs(m,i,
|
|
&j,k1)
|
|
enddo
|
|
rhs(4,i,j,k) = rhs(4,i,j,k) - lhsp__(4,0) * rhs(4,i,j,
|
|
&k1)
|
|
rhs(5,i,j,k) = rhs(5,i,j,k) - lhsm__(4,0) * rhs(5,i,j,
|
|
&k1)
|
|
endif
|
|
lhs(0,4,i,j,k) = lhs__(4,0)
|
|
lhs(1,4,i,j,k) = lhsp__(4,0)
|
|
lhs(2,4,i,j,k) = lhsm__(4,0)
|
|
lhs(0,5,i,j,k) = lhs__(5,0)
|
|
lhs(1,5,i,j,k) = lhsp__(5,0)
|
|
lhs(2,5,i,j,k) = lhsm__(5,0)
|
|
lhs__(1,0) = lhs__(1,1)
|
|
lhsp__(1,0) = lhsp__(1,1)
|
|
lhsm__(1,0) = lhsm__(1,1)
|
|
lhs__(1,1) = lhs__(1,2)
|
|
lhsp__(1,1) = lhsp__(1,2)
|
|
lhsm__(1,1) = lhsm__(1,2)
|
|
lhs__(2,0) = lhs__(2,1)
|
|
lhsp__(2,0) = lhsp__(2,1)
|
|
lhsm__(2,0) = lhsm__(2,1)
|
|
lhs__(2,1) = lhs__(2,2)
|
|
lhsp__(2,1) = lhsp__(2,2)
|
|
lhsm__(2,1) = lhsm__(2,2)
|
|
lhs__(3,0) = lhs__(3,1)
|
|
lhsp__(3,0) = lhsp__(3,1)
|
|
lhsm__(3,0) = lhsm__(3,1)
|
|
lhs__(3,1) = lhs__(3,2)
|
|
lhsp__(3,1) = lhsp__(3,2)
|
|
lhsm__(3,1) = lhsm__(3,2)
|
|
lhs__(4,0) = lhs__(4,1)
|
|
lhsp__(4,0) = lhsp__(4,1)
|
|
lhsm__(4,0) = lhsm__(4,1)
|
|
lhs__(4,1) = lhs__(4,2)
|
|
lhsp__(4,1) = lhsp__(4,2)
|
|
lhsm__(4,1) = lhsm__(4,2)
|
|
lhs__(5,0) = lhs__(5,1)
|
|
lhsp__(5,0) = lhsp__(5,1)
|
|
lhsm__(5,0) = lhsm__(5,1)
|
|
lhs__(5,1) = lhs__(5,2)
|
|
lhsp__(5,1) = lhsp__(5,2)
|
|
lhsm__(5,1) = lhsm__(5,2)
|
|
enddo
|
|
k = problem_size - 3
|
|
rhs__(1,2) = rhs(1,i,j,k + 2)
|
|
rhs__(2,2) = rhs(2,i,j,k + 2)
|
|
rhs__(3,2) = rhs(3,i,j,k + 2)
|
|
rhs__(4,2) = rhs(4,i,j,k + 2)
|
|
rhs__(5,2) = rhs(5,i,j,k + 2)
|
|
rhs__(1,1) = rhs(1,i,j,k + 1)
|
|
rhs__(2,1) = rhs(2,i,j,k + 1)
|
|
rhs__(3,1) = rhs(3,i,j,k + 1)
|
|
rhs__(4,1) = rhs(4,i,j,k + 1)
|
|
rhs__(5,1) = rhs(5,i,j,k + 1)
|
|
rhs__(1,0) = rhs(1,i,j,k)
|
|
rhs__(2,0) = rhs(2,i,j,k)
|
|
rhs__(3,0) = rhs(3,i,j,k)
|
|
rhs__(4,0) = rhs(4,i,j,k)
|
|
rhs__(5,0) = rhs(5,i,j,k)
|
|
rhs__(1,0) = rhs__(1,0) - lhs(0,4,i,j,k) * rhs__(1,1) - lhs(
|
|
&0,5,i,j,k) * rhs__(1,2)
|
|
rhs__(2,0) = rhs__(2,0) - lhs(0,4,i,j,k) * rhs__(2,1) - lhs(
|
|
&0,5,i,j,k) * rhs__(2,2)
|
|
rhs__(3,0) = rhs__(3,0) - lhs(0,4,i,j,k) * rhs__(3,1) - lhs(
|
|
&0,5,i,j,k) * rhs__(3,2)
|
|
rhs__(4,0) = rhs__(4,0) - lhs(1,4,i,j,k) * rhs__(4,1) - lhs(
|
|
&1,5,i,j,k) * rhs__(4,2)
|
|
rhs__(5,0) = rhs__(5,0) - lhs(2,4,i,j,k) * rhs__(5,1) - lhs(
|
|
&2,5,i,j,k) * rhs__(5,2)
|
|
rhs__(1,2) = rhs__(1,1)
|
|
rhs__(2,2) = rhs__(2,1)
|
|
rhs__(3,2) = rhs__(3,1)
|
|
rhs__(4,2) = rhs__(4,1)
|
|
rhs__(5,2) = rhs__(5,1)
|
|
rhs__(1,1) = rhs__(1,0)
|
|
rhs__(2,1) = rhs__(2,0)
|
|
rhs__(3,1) = rhs__(3,0)
|
|
rhs__(4,1) = rhs__(4,0)
|
|
rhs__(5,1) = rhs__(5,0)
|
|
do k = problem_size - 4,0,(-(1))
|
|
rhs__(1,0) = rhs(1,i,j,k)
|
|
rhs__(2,0) = rhs(2,i,j,k)
|
|
rhs__(3,0) = rhs(3,i,j,k)
|
|
rhs__(4,0) = rhs(4,i,j,k)
|
|
rhs__(5,0) = rhs(5,i,j,k)
|
|
rhs__(1,0) = rhs__(1,0) - lhs(0,4,i,j,k) * rhs__(1,1) - l
|
|
&hs(0,5,i,j,k) * rhs__(1,2)
|
|
rhs__(2,0) = rhs__(2,0) - lhs(0,4,i,j,k) * rhs__(2,1) - l
|
|
&hs(0,5,i,j,k) * rhs__(2,2)
|
|
rhs__(3,0) = rhs__(3,0) - lhs(0,4,i,j,k) * rhs__(3,1) - l
|
|
&hs(0,5,i,j,k) * rhs__(3,2)
|
|
rhs__(4,0) = rhs__(4,0) - lhs(1,4,i,j,k) * rhs__(4,1) - l
|
|
&hs(1,5,i,j,k) * rhs__(4,2)
|
|
rhs__(5,0) = rhs__(5,0) - lhs(2,4,i,j,k) * rhs__(5,1) - l
|
|
&hs(2,5,i,j,k) * rhs__(5,2)
|
|
xvel = us(i,j,k + 2)
|
|
yvel = vs(i,j,k + 2)
|
|
zvel = ws(i,j,k + 2)
|
|
ac = speed(i,j,k + 2)
|
|
ac2u = ac * ac
|
|
uzik1 = u(1,i,j,k + 2)
|
|
btuz = bt * uzik1
|
|
t1 = btuz / ac * (rhs__(4,2) + rhs__(5,2))
|
|
t2 = rhs__(3,2) + t1
|
|
t3 = btuz * (rhs__(4,2) - rhs__(5,2))
|
|
rhs__(3,2) = uzik1 * rhs__(1,2) + yvel * t2
|
|
rhs__(4,2) = zvel * t2 + t3
|
|
rhs__(5,2) = uzik1 * ((-(xvel)) * rhs__(2,2) + yvel * rhs
|
|
&__(1,2)) + qs(i,j,k + 2) * t2 + c2iv * ac2u * t1 + zvel * t3
|
|
rhs__(1,2) = t2
|
|
rhs__(2,2) = (-(uzik1)) * rhs__(2,2) + xvel * t2
|
|
u(1,i,j,k + 2) = u(1,i,j,k + 2) + rhs__(1,2)
|
|
u(2,i,j,k + 2) = u(2,i,j,k + 2) + rhs__(2,2)
|
|
u(3,i,j,k + 2) = u(3,i,j,k + 2) + rhs__(3,2)
|
|
u(4,i,j,k + 2) = u(4,i,j,k + 2) + rhs__(4,2)
|
|
u(5,i,j,k + 2) = u(5,i,j,k + 2) + rhs__(5,2)
|
|
rhs__(1,2) = rhs__(1,1)
|
|
rhs__(2,2) = rhs__(2,1)
|
|
rhs__(3,2) = rhs__(3,1)
|
|
rhs__(4,2) = rhs__(4,1)
|
|
rhs__(5,2) = rhs__(5,1)
|
|
rhs__(1,1) = rhs__(1,0)
|
|
rhs__(2,1) = rhs__(2,0)
|
|
rhs__(3,1) = rhs__(3,0)
|
|
rhs__(4,1) = rhs__(4,0)
|
|
rhs__(5,1) = rhs__(5,0)
|
|
enddo
|
|
xvel = us(i,j,k + 2)
|
|
yvel = vs(i,j,k + 2)
|
|
zvel = ws(i,j,k + 2)
|
|
ac = speed(i,j,k + 2)
|
|
ac2u = ac * ac
|
|
uzik1 = u(1,i,j,k + 2)
|
|
btuz = bt * uzik1
|
|
t1 = btuz / ac * (rhs__(4,2) + rhs__(5,2))
|
|
t2 = rhs__(3,2) + t1
|
|
t3 = btuz * (rhs__(4,2) - rhs__(5,2))
|
|
rhs__(3,2) = uzik1 * rhs__(1,2) + yvel * t2
|
|
rhs__(4,2) = zvel * t2 + t3
|
|
rhs__(5,2) = uzik1 * ((-(xvel)) * rhs__(2,2) + yvel * rhs__(
|
|
&1,2)) + qs(i,j,k + 2) * t2 + c2iv * ac2u * t1 + zvel * t3
|
|
rhs__(1,2) = t2
|
|
rhs__(2,2) = (-(uzik1)) * rhs__(2,2) + xvel * t2
|
|
u(1,i,j,k + 2) = u(1,i,j,k + 2) + rhs__(1,2)
|
|
u(2,i,j,k + 2) = u(2,i,j,k + 2) + rhs__(2,2)
|
|
u(3,i,j,k + 2) = u(3,i,j,k + 2) + rhs__(3,2)
|
|
u(4,i,j,k + 2) = u(4,i,j,k + 2) + rhs__(4,2)
|
|
u(5,i,j,k + 2) = u(5,i,j,k + 2) + rhs__(5,2)
|
|
enddo
|
|
enddo
|
|
!$SPF END PARALLEL_REG
|
|
! DVM$ END REGION
|
|
return
|
|
end
|
|
|