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