c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_solve(rho_i, us, speed, rhs, nx, nxmax, ny, nz) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c this function performs the solution of the approximate factorization c step in the x-direction for all five matrix components c simultaneously. The Thomas algorithm is employed to solve the c systems for the x-lines. Boundary conditions are non-periodic c--------------------------------------------------------------------- include 'header.h' integer nx, nxmax, ny, nz double precision rho_i( 0:nxmax-1,0:ny-1,0:nz-1), $ us ( 0:nxmax-1,0:ny-1,0:nz-1), $ speed( 0:nxmax-1,0:ny-1,0:nz-1), $ rhs (5,0:nxmax-1,0:ny-1,0:nz-1) integer i, j, k, i1, i2, m double precision ru1, fac1, fac2 c--------------------------------------------------------------------- c--------------------------------------------------------------------- if (timeron) call timer_start(t_xsolve) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(fac2,m,fac1,i2,i1,ru1,i,j,k) !$OMP& SHARED(comz6,comz1,comz4,comz5,c2dttx1,dttx1,dttx2,dx1,dxmax, !$OMP& c1c5,dx5,con43,dx2,c3c4,nx,ny,nz) do k = 1, nz-2 c--------------------------------------------------------------------- c Computes the left hand side for the three x-factors c--------------------------------------------------------------------- c--------------------------------------------------------------------- c zap the whole left hand side for starters c--------------------------------------------------------------------- do j = 1, ny-2 do m = 1, 5 lhs (m,0,j) = 0.0d0 lhsp(m,0,j) = 0.0d0 lhsm(m,0,j) = 0.0d0 lhs (m,nx-1,j) = 0.0d0 lhsp(m,nx-1,j) = 0.0d0 lhsm(m,nx-1,j) = 0.0d0 end do c--------------------------------------------------------------------- c next, set all diagonal values to 1. This is overkill, but c convenient c--------------------------------------------------------------------- lhs (3,0,j) = 1.0d0 lhsp(3,0,j) = 1.0d0 lhsm(3,0,j) = 1.0d0 lhs (3,nx-1,j) = 1.0d0 lhsp(3,nx-1,j) = 1.0d0 lhsm(3,nx-1,j) = 1.0d0 end do c--------------------------------------------------------------------- c first fill the lhs for the u-eigenvalue c--------------------------------------------------------------------- do j = 1, ny-2 do i = 0, nx-1 ru1 = c3c4*rho_i(i,j,k) cv(i) = us(i,j,k) rhon(i) = dmax1(dx2+con43*ru1, > dx5+c1c5*ru1, > dxmax+ru1, > dx1) end do do i = 1, nx-2 lhs(1,i,j) = 0.0d0 lhs(2,i,j) = - dttx2 * cv(i-1) - dttx1 * rhon(i-1) lhs(3,i,j) = 1.0d0 + c2dttx1 * rhon(i) lhs(4,i,j) = dttx2 * cv(i+1) - dttx1 * rhon(i+1) lhs(5,i,j) = 0.0d0 end do end do c--------------------------------------------------------------------- c add fourth order dissipation c--------------------------------------------------------------------- do j = 1, ny-2 i = 1 lhs(3,i,j) = lhs(3,i,j) + comz5 lhs(4,i,j) = lhs(4,i,j) - comz4 lhs(5,i,j) = lhs(5,i,j) + comz1 i = 2 lhs(2,i,j) = lhs(2,i,j) - comz4 lhs(3,i,j) = lhs(3,i,j) + comz6 lhs(4,i,j) = lhs(4,i,j) - comz4 lhs(5,i,j) = lhs(5,i,j) + comz1 end do do j = 1, ny-2 do i=3, nx-4 lhs(1,i,j) = lhs(1,i,j) + comz1 lhs(2,i,j) = lhs(2,i,j) - comz4 lhs(3,i,j) = lhs(3,i,j) + comz6 lhs(4,i,j) = lhs(4,i,j) - comz4 lhs(5,i,j) = lhs(5,i,j) + comz1 end do end do do j = 1, ny-2 i = nx-3 lhs(1,i,j) = lhs(1,i,j) + comz1 lhs(2,i,j) = lhs(2,i,j) - comz4 lhs(3,i,j) = lhs(3,i,j) + comz6 lhs(4,i,j) = lhs(4,i,j) - comz4 i = nx-2 lhs(1,i,j) = lhs(1,i,j) + comz1 lhs(2,i,j) = lhs(2,i,j) - comz4 lhs(3,i,j) = lhs(3,i,j) + comz5 end do c--------------------------------------------------------------------- c subsequently, fill the other factors (u+c), (u-c) by adding to c the first c--------------------------------------------------------------------- do j = 1, ny-2 do i = 1, nx-2 lhsp(1,i,j) = lhs(1,i,j) lhsp(2,i,j) = lhs(2,i,j) - > dttx2 * speed(i-1,j,k) lhsp(3,i,j) = lhs(3,i,j) lhsp(4,i,j) = lhs(4,i,j) + > dttx2 * speed(i+1,j,k) lhsp(5,i,j) = lhs(5,i,j) lhsm(1,i,j) = lhs(1,i,j) lhsm(2,i,j) = lhs(2,i,j) + > dttx2 * speed(i-1,j,k) lhsm(3,i,j) = lhs(3,i,j) lhsm(4,i,j) = lhs(4,i,j) - > dttx2 * speed(i+1,j,k) lhsm(5,i,j) = lhs(5,i,j) end do end do c--------------------------------------------------------------------- c FORWARD ELIMINATION c--------------------------------------------------------------------- c--------------------------------------------------------------------- c perform the Thomas algorithm; first, FORWARD ELIMINATION c--------------------------------------------------------------------- do j = 1, ny-2 do i = 0, nx-3 i1 = i + 1 i2 = i + 2 fac1 = 1.d0/lhs(3,i,j) lhs(4,i,j) = fac1*lhs(4,i,j) lhs(5,i,j) = fac1*lhs(5,i,j) do m = 1, 3 rhs(m,i,j,k) = fac1*rhs(m,i,j,k) end do lhs(3,i1,j) = lhs(3,i1,j) - > lhs(2,i1,j)*lhs(4,i,j) lhs(4,i1,j) = lhs(4,i1,j) - > lhs(2,i1,j)*lhs(5,i,j) do m = 1, 3 rhs(m,i1,j,k) = rhs(m,i1,j,k) - > lhs(2,i1,j)*rhs(m,i,j,k) end do lhs(2,i2,j) = lhs(2,i2,j) - > lhs(1,i2,j)*lhs(4,i,j) lhs(3,i2,j) = lhs(3,i2,j) - > lhs(1,i2,j)*lhs(5,i,j) do m = 1, 3 rhs(m,i2,j,k) = rhs(m,i2,j,k) - > lhs(1,i2,j)*rhs(m,i,j,k) end do end do end do c--------------------------------------------------------------------- c The last two rows in this zone are a bit different, c since they do not have two more rows available for the c elimination of off-diagonal entries c--------------------------------------------------------------------- do j = 1, ny-2 i = nx-2 i1 = nx-1 fac1 = 1.d0/lhs(3,i,j) lhs(4,i,j) = fac1*lhs(4,i,j) lhs(5,i,j) = fac1*lhs(5,i,j) do m = 1, 3 rhs(m,i,j,k) = fac1*rhs(m,i,j,k) end do lhs(3,i1,j) = lhs(3,i1,j) - > lhs(2,i1,j)*lhs(4,i,j) lhs(4,i1,j) = lhs(4,i1,j) - > lhs(2,i1,j)*lhs(5,i,j) do m = 1, 3 rhs(m,i1,j,k) = rhs(m,i1,j,k) - > lhs(2,i1,j)*rhs(m,i,j,k) end do c--------------------------------------------------------------------- c scale the last row immediately c--------------------------------------------------------------------- fac2 = 1.d0/lhs(3,i1,j) do m = 1, 3 rhs(m,i1,j,k) = fac2*rhs(m,i1,j,k) end do end do c--------------------------------------------------------------------- c do the u+c and the u-c factors c--------------------------------------------------------------------- do j = 1, ny-2 do i = 0, nx-3 i1 = i + 1 i2 = i + 2 m = 4 fac1 = 1.d0/lhsp(3,i,j) lhsp(4,i,j) = fac1*lhsp(4,i,j) lhsp(5,i,j) = fac1*lhsp(5,i,j) rhs(m,i,j,k) = fac1*rhs(m,i,j,k) lhsp(3,i1,j) = lhsp(3,i1,j) - > lhsp(2,i1,j)*lhsp(4,i,j) lhsp(4,i1,j) = lhsp(4,i1,j) - > lhsp(2,i1,j)*lhsp(5,i,j) rhs(m,i1,j,k) = rhs(m,i1,j,k) - > lhsp(2,i1,j)*rhs(m,i,j,k) lhsp(2,i2,j) = lhsp(2,i2,j) - > lhsp(1,i2,j)*lhsp(4,i,j) lhsp(3,i2,j) = lhsp(3,i2,j) - > lhsp(1,i2,j)*lhsp(5,i,j) rhs(m,i2,j,k) = rhs(m,i2,j,k) - > lhsp(1,i2,j)*rhs(m,i,j,k) m = 5 fac1 = 1.d0/lhsm(3,i,j) lhsm(4,i,j) = fac1*lhsm(4,i,j) lhsm(5,i,j) = fac1*lhsm(5,i,j) rhs(m,i,j,k) = fac1*rhs(m,i,j,k) lhsm(3,i1,j) = lhsm(3,i1,j) - > lhsm(2,i1,j)*lhsm(4,i,j) lhsm(4,i1,j) = lhsm(4,i1,j) - > lhsm(2,i1,j)*lhsm(5,i,j) rhs(m,i1,j,k) = rhs(m,i1,j,k) - > lhsm(2,i1,j)*rhs(m,i,j,k) lhsm(2,i2,j) = lhsm(2,i2,j) - > lhsm(1,i2,j)*lhsm(4,i,j) lhsm(3,i2,j) = lhsm(3,i2,j) - > lhsm(1,i2,j)*lhsm(5,i,j) rhs(m,i2,j,k) = rhs(m,i2,j,k) - > lhsm(1,i2,j)*rhs(m,i,j,k) end do end do c--------------------------------------------------------------------- c And again the last two rows separately c--------------------------------------------------------------------- do j = 1, ny-2 i = nx-2 i1 = nx-1 m = 4 fac1 = 1.d0/lhsp(3,i,j) lhsp(4,i,j) = fac1*lhsp(4,i,j) lhsp(5,i,j) = fac1*lhsp(5,i,j) rhs(m,i,j,k) = fac1*rhs(m,i,j,k) lhsp(3,i1,j) = lhsp(3,i1,j) - > lhsp(2,i1,j)*lhsp(4,i,j) lhsp(4,i1,j) = lhsp(4,i1,j) - > lhsp(2,i1,j)*lhsp(5,i,j) rhs(m,i1,j,k) = rhs(m,i1,j,k) - > lhsp(2,i1,j)*rhs(m,i,j,k) m = 5 fac1 = 1.d0/lhsm(3,i,j) lhsm(4,i,j) = fac1*lhsm(4,i,j) lhsm(5,i,j) = fac1*lhsm(5,i,j) rhs(m,i,j,k) = fac1*rhs(m,i,j,k) lhsm(3,i1,j) = lhsm(3,i1,j) - > lhsm(2,i1,j)*lhsm(4,i,j) lhsm(4,i1,j) = lhsm(4,i1,j) - > lhsm(2,i1,j)*lhsm(5,i,j) rhs(m,i1,j,k) = rhs(m,i1,j,k) - > lhsm(2,i1,j)*rhs(m,i,j,k) c--------------------------------------------------------------------- c Scale the last row immediately c--------------------------------------------------------------------- rhs(4,i1,j,k) = rhs(4,i1,j,k)/lhsp(3,i1,j) rhs(5,i1,j,k) = rhs(5,i1,j,k)/lhsm(3,i1,j) end do c--------------------------------------------------------------------- c BACKSUBSTITUTION c--------------------------------------------------------------------- do j = 1, ny-2 i = nx-2 i1 = nx-1 do m = 1, 3 rhs(m,i,j,k) = rhs(m,i,j,k) - > lhs(4,i,j)*rhs(m,i1,j,k) end do rhs(4,i,j,k) = rhs(4,i,j,k) - > lhsp(4,i,j)*rhs(4,i1,j,k) rhs(5,i,j,k) = rhs(5,i,j,k) - > lhsm(4,i,j)*rhs(5,i1,j,k) end do c--------------------------------------------------------------------- c The first three factors c--------------------------------------------------------------------- do j = 1, ny-2 do i = nx-3, 0, -1 i1 = i + 1 i2 = i + 2 do m = 1, 3 rhs(m,i,j,k) = rhs(m,i,j,k) - > lhs(4,i,j)*rhs(m,i1,j,k) - > lhs(5,i,j)*rhs(m,i2,j,k) end do c--------------------------------------------------------------------- c And the remaining two c--------------------------------------------------------------------- rhs(4,i,j,k) = rhs(4,i,j,k) - > lhsp(4,i,j)*rhs(4,i1,j,k) - > lhsp(5,i,j)*rhs(4,i2,j,k) rhs(5,i,j,k) = rhs(5,i,j,k) - > lhsm(4,i,j)*rhs(5,i1,j,k) - > lhsm(5,i,j)*rhs(5,i2,j,k) end do end do end do !$OMP END PARALLEL DO if (timeron) call timer_stop(t_xsolve) c--------------------------------------------------------------------- c Do the block-diagonal inversion c--------------------------------------------------------------------- if (timeron) call timer_start(t_ninvr) call ninvr(rhs, nx, nxmax, ny, nz) if (timeron) call timer_stop(t_ninvr) return end