c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine erhs c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c c compute the right hand side based on exact solution c c--------------------------------------------------------------------- implicit none include 'applu.incl' c--------------------------------------------------------------------- c local variables c--------------------------------------------------------------------- integer i, j, k, m double precision xi, eta, zeta double precision q double precision u21, u31, u41 double precision tmp double precision u21i, u31i, u41i, u51i double precision u21j, u31j, u41j, u51j double precision u21k, u31k, u41k, u51k double precision u21im1, u31im1, u41im1, u51im1 double precision u21jm1, u31jm1, u41jm1, u51jm1 double precision u21km1, u31km1, u41km1, u51km1 do k = 1, nz do j = 1, ny do i = 1, nx do m = 1, 5 frct( m, i, j, k ) = 0.0d+00 end do end do end do end do do k = 1, nz zeta = ( dble(k-1) ) / ( nz - 1 ) do j = 1, ny eta = ( dble(j-1) ) / ( ny0 - 1 ) do i = 1, nx xi = ( dble(i-1) ) / ( nx0 - 1 ) do m = 1, 5 rsd(m,i,j,k) = ce(m,1) > + (ce(m,2) > + (ce(m,5) > + (ce(m,8) > + ce(m,11) * xi) * xi) * xi) * xi > + (ce(m,3) > + (ce(m,6) > + (ce(m,9) > + ce(m,12) * eta) * eta) * eta) * eta > + (ce(m,4) > + (ce(m,7) > + (ce(m,10) > + ce(m,13) * zeta) * zeta) * zeta) * zeta end do end do end do end do c--------------------------------------------------------------------- c xi-direction flux differences c--------------------------------------------------------------------- do k = 2, nz - 1 do j = jst, jend do i = 1, nx flux(1,i) = rsd(2,i,j,k) u21 = rsd(2,i,j,k) / rsd(1,i,j,k) q = 0.50d+00 * ( rsd(2,i,j,k) * rsd(2,i,j,k) > + rsd(3,i,j,k) * rsd(3,i,j,k) > + rsd(4,i,j,k) * rsd(4,i,j,k) ) > / rsd(1,i,j,k) flux(2,i) = rsd(2,i,j,k) * u21 + c2 * > ( rsd(5,i,j,k) - q ) flux(3,i) = rsd(3,i,j,k) * u21 flux(4,i) = rsd(4,i,j,k) * u21 flux(5,i) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u21 end do do i = ist, iend do m = 1, 5 frct(m,i,j,k) = frct(m,i,j,k) > - tx2 * ( flux(m,i+1) - flux(m,i-1) ) end do end do do i = ist, nx tmp = 1.0d+00 / rsd(1,i,j,k) u21i = tmp * rsd(2,i,j,k) u31i = tmp * rsd(3,i,j,k) u41i = tmp * rsd(4,i,j,k) u51i = tmp * rsd(5,i,j,k) tmp = 1.0d+00 / rsd(1,i-1,j,k) u21im1 = tmp * rsd(2,i-1,j,k) u31im1 = tmp * rsd(3,i-1,j,k) u41im1 = tmp * rsd(4,i-1,j,k) u51im1 = tmp * rsd(5,i-1,j,k) flux(2,i) = (4.0d+00/3.0d+00) * tx3 * > ( u21i - u21im1 ) flux(3,i) = tx3 * ( u31i - u31im1 ) flux(4,i) = tx3 * ( u41i - u41im1 ) flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 ) > * tx3 * ( ( u21i **2 + u31i **2 + u41i **2 ) > - ( u21im1**2 + u31im1**2 + u41im1**2 ) ) > + (1.0d+00/6.0d+00) > * tx3 * ( u21i**2 - u21im1**2 ) > + c1 * c5 * tx3 * ( u51i - u51im1 ) end do do i = ist, iend frct(1,i,j,k) = frct(1,i,j,k) > + dx1 * tx1 * ( rsd(1,i-1,j,k) > - 2.0d+00 * rsd(1,i,j,k) > + rsd(1,i+1,j,k) ) frct(2,i,j,k) = frct(2,i,j,k) > + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) ) > + dx2 * tx1 * ( rsd(2,i-1,j,k) > - 2.0d+00 * rsd(2,i,j,k) > + rsd(2,i+1,j,k) ) frct(3,i,j,k) = frct(3,i,j,k) > + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) ) > + dx3 * tx1 * ( rsd(3,i-1,j,k) > - 2.0d+00 * rsd(3,i,j,k) > + rsd(3,i+1,j,k) ) frct(4,i,j,k) = frct(4,i,j,k) > + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) ) > + dx4 * tx1 * ( rsd(4,i-1,j,k) > - 2.0d+00 * rsd(4,i,j,k) > + rsd(4,i+1,j,k) ) frct(5,i,j,k) = frct(5,i,j,k) > + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) ) > + dx5 * tx1 * ( rsd(5,i-1,j,k) > - 2.0d+00 * rsd(5,i,j,k) > + rsd(5,i+1,j,k) ) end do c--------------------------------------------------------------------- c Fourth-order dissipation c--------------------------------------------------------------------- do m = 1, 5 frct(m,2,j,k) = frct(m,2,j,k) > - dssp * ( + 5.0d+00 * rsd(m,2,j,k) > - 4.0d+00 * rsd(m,3,j,k) > + rsd(m,4,j,k) ) frct(m,3,j,k) = frct(m,3,j,k) > - dssp * ( - 4.0d+00 * rsd(m,2,j,k) > + 6.0d+00 * rsd(m,3,j,k) > - 4.0d+00 * rsd(m,4,j,k) > + rsd(m,5,j,k) ) end do do i = 4, nx - 3 do m = 1, 5 frct(m,i,j,k) = frct(m,i,j,k) > - dssp * ( rsd(m,i-2,j,k) > - 4.0d+00 * rsd(m,i-1,j,k) > + 6.0d+00 * rsd(m,i,j,k) > - 4.0d+00 * rsd(m,i+1,j,k) > + rsd(m,i+2,j,k) ) end do end do do m = 1, 5 frct(m,nx-2,j,k) = frct(m,nx-2,j,k) > - dssp * ( rsd(m,nx-4,j,k) > - 4.0d+00 * rsd(m,nx-3,j,k) > + 6.0d+00 * rsd(m,nx-2,j,k) > - 4.0d+00 * rsd(m,nx-1,j,k) ) frct(m,nx-1,j,k) = frct(m,nx-1,j,k) > - dssp * ( rsd(m,nx-3,j,k) > - 4.0d+00 * rsd(m,nx-2,j,k) > + 5.0d+00 * rsd(m,nx-1,j,k) ) end do end do end do c--------------------------------------------------------------------- c eta-direction flux differences c--------------------------------------------------------------------- do k = 2, nz - 1 do i = ist, iend do j = 1, ny flux(1,j) = rsd(3,i,j,k) u31 = rsd(3,i,j,k) / rsd(1,i,j,k) q = 0.50d+00 * ( rsd(2,i,j,k) * rsd(2,i,j,k) > + rsd(3,i,j,k) * rsd(3,i,j,k) > + rsd(4,i,j,k) * rsd(4,i,j,k) ) > / rsd(1,i,j,k) flux(2,j) = rsd(2,i,j,k) * u31 flux(3,j) = rsd(3,i,j,k) * u31 + c2 * > ( rsd(5,i,j,k) - q ) flux(4,j) = rsd(4,i,j,k) * u31 flux(5,j) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u31 end do do j = jst, jend do m = 1, 5 frct(m,i,j,k) = frct(m,i,j,k) > - ty2 * ( flux(m,j+1) - flux(m,j-1) ) end do end do do j = jst, ny tmp = 1.0d+00 / rsd(1,i,j,k) u21j = tmp * rsd(2,i,j,k) u31j = tmp * rsd(3,i,j,k) u41j = tmp * rsd(4,i,j,k) u51j = tmp * rsd(5,i,j,k) tmp = 1.0d+00 / rsd(1,i,j-1,k) u21jm1 = tmp * rsd(2,i,j-1,k) u31jm1 = tmp * rsd(3,i,j-1,k) u41jm1 = tmp * rsd(4,i,j-1,k) u51jm1 = tmp * rsd(5,i,j-1,k) flux(2,j) = ty3 * ( u21j - u21jm1 ) flux(3,j) = (4.0d+00/3.0d+00) * ty3 * > ( u31j - u31jm1 ) flux(4,j) = ty3 * ( u41j - u41jm1 ) flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 ) > * ty3 * ( ( u21j **2 + u31j **2 + u41j **2 ) > - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) ) > + (1.0d+00/6.0d+00) > * ty3 * ( u31j**2 - u31jm1**2 ) > + c1 * c5 * ty3 * ( u51j - u51jm1 ) end do do j = jst, jend frct(1,i,j,k) = frct(1,i,j,k) > + dy1 * ty1 * ( rsd(1,i,j-1,k) > - 2.0d+00 * rsd(1,i,j,k) > + rsd(1,i,j+1,k) ) frct(2,i,j,k) = frct(2,i,j,k) > + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) ) > + dy2 * ty1 * ( rsd(2,i,j-1,k) > - 2.0d+00 * rsd(2,i,j,k) > + rsd(2,i,j+1,k) ) frct(3,i,j,k) = frct(3,i,j,k) > + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) ) > + dy3 * ty1 * ( rsd(3,i,j-1,k) > - 2.0d+00 * rsd(3,i,j,k) > + rsd(3,i,j+1,k) ) frct(4,i,j,k) = frct(4,i,j,k) > + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) ) > + dy4 * ty1 * ( rsd(4,i,j-1,k) > - 2.0d+00 * rsd(4,i,j,k) > + rsd(4,i,j+1,k) ) frct(5,i,j,k) = frct(5,i,j,k) > + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) ) > + dy5 * ty1 * ( rsd(5,i,j-1,k) > - 2.0d+00 * rsd(5,i,j,k) > + rsd(5,i,j+1,k) ) end do c--------------------------------------------------------------------- c fourth-order dissipation c--------------------------------------------------------------------- do m = 1, 5 frct(m,i,2,k) = frct(m,i,2,k) > - dssp * ( + 5.0d+00 * rsd(m,i,2,k) > - 4.0d+00 * rsd(m,i,3,k) > + rsd(m,i,4,k) ) frct(m,i,3,k) = frct(m,i,3,k) > - dssp * ( - 4.0d+00 * rsd(m,i,2,k) > + 6.0d+00 * rsd(m,i,3,k) > - 4.0d+00 * rsd(m,i,4,k) > + rsd(m,i,5,k) ) end do do j = 4, ny - 3 do m = 1, 5 frct(m,i,j,k) = frct(m,i,j,k) > - dssp * ( rsd(m,i,j-2,k) > - 4.0d+00 * rsd(m,i,j-1,k) > + 6.0d+00 * rsd(m,i,j,k) > - 4.0d+00 * rsd(m,i,j+1,k) > + rsd(m,i,j+2,k) ) end do end do do m = 1, 5 frct(m,i,ny-2,k) = frct(m,i,ny-2,k) > - dssp * ( rsd(m,i,ny-4,k) > - 4.0d+00 * rsd(m,i,ny-3,k) > + 6.0d+00 * rsd(m,i,ny-2,k) > - 4.0d+00 * rsd(m,i,ny-1,k) ) frct(m,i,ny-1,k) = frct(m,i,ny-1,k) > - dssp * ( rsd(m,i,ny-3,k) > - 4.0d+00 * rsd(m,i,ny-2,k) > + 5.0d+00 * rsd(m,i,ny-1,k) ) end do end do end do c--------------------------------------------------------------------- c zeta-direction flux differences c--------------------------------------------------------------------- do j = jst, jend do i = ist, iend do k = 1, nz flux(1,k) = rsd(4,i,j,k) u41 = rsd(4,i,j,k) / rsd(1,i,j,k) q = 0.50d+00 * ( rsd(2,i,j,k) * rsd(2,i,j,k) > + rsd(3,i,j,k) * rsd(3,i,j,k) > + rsd(4,i,j,k) * rsd(4,i,j,k) ) > / rsd(1,i,j,k) flux(2,k) = rsd(2,i,j,k) * u41 flux(3,k) = rsd(3,i,j,k) * u41 flux(4,k) = rsd(4,i,j,k) * u41 + c2 * > ( rsd(5,i,j,k) - q ) flux(5,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u41 end do do k = 2, nz - 1 do m = 1, 5 frct(m,i,j,k) = frct(m,i,j,k) > - tz2 * ( flux(m,k+1) - flux(m,k-1) ) end do end do do k = 2, nz tmp = 1.0d+00 / rsd(1,i,j,k) u21k = tmp * rsd(2,i,j,k) u31k = tmp * rsd(3,i,j,k) u41k = tmp * rsd(4,i,j,k) u51k = tmp * rsd(5,i,j,k) tmp = 1.0d+00 / rsd(1,i,j,k-1) u21km1 = tmp * rsd(2,i,j,k-1) u31km1 = tmp * rsd(3,i,j,k-1) u41km1 = tmp * rsd(4,i,j,k-1) u51km1 = tmp * rsd(5,i,j,k-1) flux(2,k) = tz3 * ( u21k - u21km1 ) flux(3,k) = tz3 * ( u31k - u31km1 ) flux(4,k) = (4.0d+00/3.0d+00) * tz3 * ( u41k > - u41km1 ) flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 ) > * tz3 * ( ( u21k **2 + u31k **2 + u41k **2 ) > - ( u21km1**2 + u31km1**2 + u41km1**2 ) ) > + (1.0d+00/6.0d+00) > * tz3 * ( u41k**2 - u41km1**2 ) > + c1 * c5 * tz3 * ( u51k - u51km1 ) end do do k = 2, nz - 1 frct(1,i,j,k) = frct(1,i,j,k) > + dz1 * tz1 * ( rsd(1,i,j,k+1) > - 2.0d+00 * rsd(1,i,j,k) > + rsd(1,i,j,k-1) ) frct(2,i,j,k) = frct(2,i,j,k) > + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) ) > + dz2 * tz1 * ( rsd(2,i,j,k+1) > - 2.0d+00 * rsd(2,i,j,k) > + rsd(2,i,j,k-1) ) frct(3,i,j,k) = frct(3,i,j,k) > + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) ) > + dz3 * tz1 * ( rsd(3,i,j,k+1) > - 2.0d+00 * rsd(3,i,j,k) > + rsd(3,i,j,k-1) ) frct(4,i,j,k) = frct(4,i,j,k) > + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) ) > + dz4 * tz1 * ( rsd(4,i,j,k+1) > - 2.0d+00 * rsd(4,i,j,k) > + rsd(4,i,j,k-1) ) frct(5,i,j,k) = frct(5,i,j,k) > + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) ) > + dz5 * tz1 * ( rsd(5,i,j,k+1) > - 2.0d+00 * rsd(5,i,j,k) > + rsd(5,i,j,k-1) ) end do c--------------------------------------------------------------------- c fourth-order dissipation c--------------------------------------------------------------------- do m = 1, 5 frct(m,i,j,2) = frct(m,i,j,2) > - dssp * ( + 5.0d+00 * rsd(m,i,j,2) > - 4.0d+00 * rsd(m,i,j,3) > + rsd(m,i,j,4) ) frct(m,i,j,3) = frct(m,i,j,3) > - dssp * (- 4.0d+00 * rsd(m,i,j,2) > + 6.0d+00 * rsd(m,i,j,3) > - 4.0d+00 * rsd(m,i,j,4) > + rsd(m,i,j,5) ) end do do k = 4, nz - 3 do m = 1, 5 frct(m,i,j,k) = frct(m,i,j,k) > - dssp * ( rsd(m,i,j,k-2) > - 4.0d+00 * rsd(m,i,j,k-1) > + 6.0d+00 * rsd(m,i,j,k) > - 4.0d+00 * rsd(m,i,j,k+1) > + rsd(m,i,j,k+2) ) end do end do do m = 1, 5 frct(m,i,j,nz-2) = frct(m,i,j,nz-2) > - dssp * ( rsd(m,i,j,nz-4) > - 4.0d+00 * rsd(m,i,j,nz-3) > + 6.0d+00 * rsd(m,i,j,nz-2) > - 4.0d+00 * rsd(m,i,j,nz-1) ) frct(m,i,j,nz-1) = frct(m,i,j,nz-1) > - dssp * ( rsd(m,i,j,nz-3) > - 4.0d+00 * rsd(m,i,j,nz-2) > + 5.0d+00 * rsd(m,i,j,nz-1) ) end do end do end do return end