blob: 9c07a35538449561cefa920939ea89e0d7f4fd04 [file] [log] [blame]
c---------------------------------------------------------------------
c---------------------------------------------------------------------
subroutine lhsy(c)
c---------------------------------------------------------------------
c---------------------------------------------------------------------
c---------------------------------------------------------------------
c This function computes the left hand side for the three y-factors
c---------------------------------------------------------------------
include 'header.h'
double precision ru1
integer i, j, k, c
c---------------------------------------------------------------------
c treat only cell c
c---------------------------------------------------------------------
c---------------------------------------------------------------------
c first fill the lhs for the u-eigenvalue
c---------------------------------------------------------------------
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
do j = start(2,c)-1, cell_size(2,c)-end(2,c)
ru1 = c3c4*rho_i(i,j,k,c)
cv(j) = vs(i,j,k,c)
rhoq(j) = dmax1( dy3 + con43 * ru1,
> dy5 + c1c5*ru1,
> dymax + ru1,
> dy1)
end do
do j = start(2,c), cell_size(2,c)-end(2,c)-1
lhs(i,j,k,1,c) = 0.0d0
lhs(i,j,k,2,c) = -dtty2 * cv(j-1) - dtty1 * rhoq(j-1)
lhs(i,j,k,3,c) = 1.0 + c2dtty1 * rhoq(j)
lhs(i,j,k,4,c) = dtty2 * cv(j+1) - dtty1 * rhoq(j+1)
lhs(i,j,k,5,c) = 0.0d0
end do
end do
end do
c---------------------------------------------------------------------
c add fourth order dissipation
c---------------------------------------------------------------------
if (start(2,c) .gt. 0) then
j = 1
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz5
lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
lhs(i,j+1,k,2,c) = lhs(i,j+1,k,2,c) - comz4
lhs(i,j+1,k,3,c) = lhs(i,j+1,k,3,c) + comz6
lhs(i,j+1,k,4,c) = lhs(i,j+1,k,4,c) - comz4
lhs(i,j+1,k,5,c) = lhs(i,j+1,k,5,c) + comz1
end do
end do
endif
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j=3*start(2,c), cell_size(2,c)-3*end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
end do
end do
end do
if (end(2,c) .gt. 0) then
j = cell_size(2,c)-3
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
lhs(i,j+1,k,1,c) = lhs(i,j+1,k,1,c) + comz1
lhs(i,j+1,k,2,c) = lhs(i,j+1,k,2,c) - comz4
lhs(i,j+1,k,3,c) = lhs(i,j+1,k,3,c) + comz5
end do
end do
endif
c---------------------------------------------------------------------
c subsequently, do the other two factors
c---------------------------------------------------------------------
do k = start(3,c), cell_size(3,c)-end(3,c)-1
do j = start(2,c), cell_size(2,c)-end(2,c)-1
do i = start(1,c), cell_size(1,c)-end(1,c)-1
lhs(i,j,k,1+5,c) = lhs(i,j,k,1,c)
lhs(i,j,k,2+5,c) = lhs(i,j,k,2,c) -
> dtty2 * speed(i,j-1,k,c)
lhs(i,j,k,3+5,c) = lhs(i,j,k,3,c)
lhs(i,j,k,4+5,c) = lhs(i,j,k,4,c) +
> dtty2 * speed(i,j+1,k,c)
lhs(i,j,k,5+5,c) = lhs(i,j,k,5,c)
lhs(i,j,k,1+10,c) = lhs(i,j,k,1,c)
lhs(i,j,k,2+10,c) = lhs(i,j,k,2,c) +
> dtty2 * speed(i,j-1,k,c)
lhs(i,j,k,3+10,c) = lhs(i,j,k,3,c)
lhs(i,j,k,4+10,c) = lhs(i,j,k,4,c) -
> dtty2 * speed(i,j+1,k,c)
lhs(i,j,k,5+10,c) = lhs(i,j,k,5,c)
end do
end do
end do
return
end