c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine compute_rhs c--------------------------------------------------------------------- c--------------------------------------------------------------------- include 'header.h' integer c, i, j, k, m double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1, > wijk, wp1, wm1 c--------------------------------------------------------------------- c loop over all cells owned by this node c--------------------------------------------------------------------- do c = 1, ncells c--------------------------------------------------------------------- c compute the reciprocal of density, and the kinetic energy, c and the speed of sound. c--------------------------------------------------------------------- do k = -1, cell_size(3,c) do j = -1, cell_size(2,c) do i = -1, cell_size(1,c) rho_inv = 1.0d0/u(i,j,k,1,c) rho_i(i,j,k,c) = rho_inv us(i,j,k,c) = u(i,j,k,2,c) * rho_inv vs(i,j,k,c) = u(i,j,k,3,c) * rho_inv ws(i,j,k,c) = u(i,j,k,4,c) * rho_inv square(i,j,k,c) = 0.5d0* ( > u(i,j,k,2,c)*u(i,j,k,2,c) + > u(i,j,k,3,c)*u(i,j,k,3,c) + > u(i,j,k,4,c)*u(i,j,k,4,c) ) * rho_inv qs(i,j,k,c) = square(i,j,k,c) * rho_inv c--------------------------------------------------------------------- c (don't need speed and ainx until the lhs computation) c--------------------------------------------------------------------- aux = c1c2*rho_inv* (u(i,j,k,5,c) - square(i,j,k,c)) aux = dsqrt(aux) speed(i,j,k,c) = aux ainv(i,j,k,c) = 1.0d0/aux end do end do end do c--------------------------------------------------------------------- c copy the exact forcing term to the right hand side; because c this forcing term is known, we can store it on the whole of every c cell, including the boundary c--------------------------------------------------------------------- do m = 1, 5 do k = 0, cell_size(3,c)-1 do j = 0, cell_size(2,c)-1 do i = 0, cell_size(1,c)-1 rhs(i,j,k,m,c) = forcing(i,j,k,m,c) end do end do end do end do c--------------------------------------------------------------------- c compute xi-direction fluxes 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 uijk = us(i,j,k,c) up1 = us(i+1,j,k,c) um1 = us(i-1,j,k,c) rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dx1tx1 * > (u(i+1,j,k,1,c) - 2.0d0*u(i,j,k,1,c) + > u(i-1,j,k,1,c)) - > tx2 * (u(i+1,j,k,2,c) - u(i-1,j,k,2,c)) rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dx2tx1 * > (u(i+1,j,k,2,c) - 2.0d0*u(i,j,k,2,c) + > u(i-1,j,k,2,c)) + > xxcon2*con43 * (up1 - 2.0d0*uijk + um1) - > tx2 * (u(i+1,j,k,2,c)*up1 - > u(i-1,j,k,2,c)*um1 + > (u(i+1,j,k,5,c)- square(i+1,j,k,c)- > u(i-1,j,k,5,c)+ square(i-1,j,k,c))* > c2) rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dx3tx1 * > (u(i+1,j,k,3,c) - 2.0d0*u(i,j,k,3,c) + > u(i-1,j,k,3,c)) + > xxcon2 * (vs(i+1,j,k,c) - 2.0d0*vs(i,j,k,c) + > vs(i-1,j,k,c)) - > tx2 * (u(i+1,j,k,3,c)*up1 - > u(i-1,j,k,3,c)*um1) rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dx4tx1 * > (u(i+1,j,k,4,c) - 2.0d0*u(i,j,k,4,c) + > u(i-1,j,k,4,c)) + > xxcon2 * (ws(i+1,j,k,c) - 2.0d0*ws(i,j,k,c) + > ws(i-1,j,k,c)) - > tx2 * (u(i+1,j,k,4,c)*up1 - > u(i-1,j,k,4,c)*um1) rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dx5tx1 * > (u(i+1,j,k,5,c) - 2.0d0*u(i,j,k,5,c) + > u(i-1,j,k,5,c)) + > xxcon3 * (qs(i+1,j,k,c) - 2.0d0*qs(i,j,k,c) + > qs(i-1,j,k,c)) + > xxcon4 * (up1*up1 - 2.0d0*uijk*uijk + > um1*um1) + > xxcon5 * (u(i+1,j,k,5,c)*rho_i(i+1,j,k,c) - > 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) + > u(i-1,j,k,5,c)*rho_i(i-1,j,k,c)) - > tx2 * ( (c1*u(i+1,j,k,5,c) - > c2*square(i+1,j,k,c))*up1 - > (c1*u(i-1,j,k,5,c) - > c2*square(i-1,j,k,c))*um1 ) end do end do end do c--------------------------------------------------------------------- c add fourth order xi-direction dissipation c--------------------------------------------------------------------- if (start(1,c) .gt. 0) then i = 1 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp * > ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) + > u(i+2,j,k,m,c)) end do end do end do i = 2 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > (-4.0d0*u(i-1,j,k,m,c) + 6.0d0*u(i,j,k,m,c) - > 4.0d0*u(i+1,j,k,m,c) + u(i+2,j,k,m,c)) end do end do end do endif do m = 1, 5 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 = 3*start(1,c),cell_size(1,c)-3*end(1,c)-1 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) + > 6.0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) + > u(i+2,j,k,m,c) ) end do end do end do end do if (end(1,c) .gt. 0) then i = cell_size(1,c)-3 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) + > 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) ) end do end do end do i = cell_size(1,c)-2 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i-2,j,k,m,c) - 4.d0*u(i-1,j,k,m,c) + > 5.d0*u(i,j,k,m,c) ) end do end do end do endif c--------------------------------------------------------------------- c compute eta-direction fluxes 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 vijk = vs(i,j,k,c) vp1 = vs(i,j+1,k,c) vm1 = vs(i,j-1,k,c) rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dy1ty1 * > (u(i,j+1,k,1,c) - 2.0d0*u(i,j,k,1,c) + > u(i,j-1,k,1,c)) - > ty2 * (u(i,j+1,k,3,c) - u(i,j-1,k,3,c)) rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dy2ty1 * > (u(i,j+1,k,2,c) - 2.0d0*u(i,j,k,2,c) + > u(i,j-1,k,2,c)) + > yycon2 * (us(i,j+1,k,c) - 2.0d0*us(i,j,k,c) + > us(i,j-1,k,c)) - > ty2 * (u(i,j+1,k,2,c)*vp1 - > u(i,j-1,k,2,c)*vm1) rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dy3ty1 * > (u(i,j+1,k,3,c) - 2.0d0*u(i,j,k,3,c) + > u(i,j-1,k,3,c)) + > yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) - > ty2 * (u(i,j+1,k,3,c)*vp1 - > u(i,j-1,k,3,c)*vm1 + > (u(i,j+1,k,5,c) - square(i,j+1,k,c) - > u(i,j-1,k,5,c) + square(i,j-1,k,c)) > *c2) rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dy4ty1 * > (u(i,j+1,k,4,c) - 2.0d0*u(i,j,k,4,c) + > u(i,j-1,k,4,c)) + > yycon2 * (ws(i,j+1,k,c) - 2.0d0*ws(i,j,k,c) + > ws(i,j-1,k,c)) - > ty2 * (u(i,j+1,k,4,c)*vp1 - > u(i,j-1,k,4,c)*vm1) rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dy5ty1 * > (u(i,j+1,k,5,c) - 2.0d0*u(i,j,k,5,c) + > u(i,j-1,k,5,c)) + > yycon3 * (qs(i,j+1,k,c) - 2.0d0*qs(i,j,k,c) + > qs(i,j-1,k,c)) + > yycon4 * (vp1*vp1 - 2.0d0*vijk*vijk + > vm1*vm1) + > yycon5 * (u(i,j+1,k,5,c)*rho_i(i,j+1,k,c) - > 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) + > u(i,j-1,k,5,c)*rho_i(i,j-1,k,c)) - > ty2 * ((c1*u(i,j+1,k,5,c) - > c2*square(i,j+1,k,c)) * vp1 - > (c1*u(i,j-1,k,5,c) - > c2*square(i,j-1,k,c)) * vm1) end do end do end do c--------------------------------------------------------------------- c add fourth order eta-direction dissipation c--------------------------------------------------------------------- if (start(2,c) .gt. 0) then j = 1 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp * > ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) + > u(i,j+2,k,m,c)) end do end do end do j = 2 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > (-4.0d0*u(i,j-1,k,m,c) + 6.0d0*u(i,j,k,m,c) - > 4.0d0*u(i,j+1,k,m,c) + u(i,j+2,k,m,c)) end do end do end do endif do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) + > 6.0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) + > u(i,j+2,k,m,c) ) end do end do end do end do if (end(2,c) .gt. 0) then j = cell_size(2,c)-3 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) + > 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) ) end do end do end do j = cell_size(2,c)-2 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i,j-2,k,m,c) - 4.d0*u(i,j-1,k,m,c) + > 5.d0*u(i,j,k,m,c) ) end do end do end do endif c--------------------------------------------------------------------- c compute zeta-direction fluxes 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 wijk = ws(i,j,k,c) wp1 = ws(i,j,k+1,c) wm1 = ws(i,j,k-1,c) rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dz1tz1 * > (u(i,j,k+1,1,c) - 2.0d0*u(i,j,k,1,c) + > u(i,j,k-1,1,c)) - > tz2 * (u(i,j,k+1,4,c) - u(i,j,k-1,4,c)) rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dz2tz1 * > (u(i,j,k+1,2,c) - 2.0d0*u(i,j,k,2,c) + > u(i,j,k-1,2,c)) + > zzcon2 * (us(i,j,k+1,c) - 2.0d0*us(i,j,k,c) + > us(i,j,k-1,c)) - > tz2 * (u(i,j,k+1,2,c)*wp1 - > u(i,j,k-1,2,c)*wm1) rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dz3tz1 * > (u(i,j,k+1,3,c) - 2.0d0*u(i,j,k,3,c) + > u(i,j,k-1,3,c)) + > zzcon2 * (vs(i,j,k+1,c) - 2.0d0*vs(i,j,k,c) + > vs(i,j,k-1,c)) - > tz2 * (u(i,j,k+1,3,c)*wp1 - > u(i,j,k-1,3,c)*wm1) rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dz4tz1 * > (u(i,j,k+1,4,c) - 2.0d0*u(i,j,k,4,c) + > u(i,j,k-1,4,c)) + > zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) - > tz2 * (u(i,j,k+1,4,c)*wp1 - > u(i,j,k-1,4,c)*wm1 + > (u(i,j,k+1,5,c) - square(i,j,k+1,c) - > u(i,j,k-1,5,c) + square(i,j,k-1,c)) > *c2) rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dz5tz1 * > (u(i,j,k+1,5,c) - 2.0d0*u(i,j,k,5,c) + > u(i,j,k-1,5,c)) + > zzcon3 * (qs(i,j,k+1,c) - 2.0d0*qs(i,j,k,c) + > qs(i,j,k-1,c)) + > zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk + > wm1*wm1) + > zzcon5 * (u(i,j,k+1,5,c)*rho_i(i,j,k+1,c) - > 2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) + > u(i,j,k-1,5,c)*rho_i(i,j,k-1,c)) - > tz2 * ( (c1*u(i,j,k+1,5,c) - > c2*square(i,j,k+1,c))*wp1 - > (c1*u(i,j,k-1,5,c) - > c2*square(i,j,k-1,c))*wm1) end do end do end do c--------------------------------------------------------------------- c add fourth order zeta-direction dissipation c--------------------------------------------------------------------- if (start(3,c) .gt. 0) then k = 1 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp * > ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) + > u(i,j,k+2,m,c)) end do end do end do k = 2 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > (-4.0d0*u(i,j,k-1,m,c) + 6.0d0*u(i,j,k,m,c) - > 4.0d0*u(i,j,k+1,m,c) + u(i,j,k+2,m,c)) end do end do end do endif do m = 1, 5 do k = 3*start(3,c), cell_size(3,c)-3*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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) + > 6.0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) + > u(i,j,k+2,m,c) ) end do end do end do end do if (end(3,c) .gt. 0) then k = cell_size(3,c)-3 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) + > 6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) ) end do end do end do k = cell_size(3,c)-2 do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp * > ( u(i,j,k-2,m,c) - 4.d0*u(i,j,k-1,m,c) + > 5.d0*u(i,j,k,m,c) ) end do end do end do endif do m = 1, 5 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 rhs(i,j,k,m,c) = rhs(i,j,k,m,c) * dt end do end do end do end do end do return end