+++ /dev/null
-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
- integer iglob, jglob
- integer iex
- integer L1, L2
- integer ist1, iend1
- integer jst1, jend1
- double precision dsspm
- 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
-
- dsspm = dssp
-
-
- 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
- jglob = jpt + j
- eta = ( dble(jglob-1) ) / ( ny0 - 1 )
- do i = 1, nx
- iglob = ipt + i
- xi = ( dble(iglob-1) ) / ( nx0 - 1 )
- do m = 1, 5
- rsd(m,i,j,k) = ce(m,1)
- > + ce(m,2) * xi
- > + ce(m,3) * eta
- > + ce(m,4) * zeta
- > + ce(m,5) * xi * xi
- > + ce(m,6) * eta * eta
- > + ce(m,7) * zeta * zeta
- > + ce(m,8) * xi * xi * xi
- > + ce(m,9) * eta * eta * eta
- > + ce(m,10) * zeta * zeta * zeta
- > + ce(m,11) * xi * xi * xi * xi
- > + ce(m,12) * eta * eta * eta * eta
- > + ce(m,13) * zeta * zeta * zeta * zeta
- end do
- end do
- end do
- end do
-
-c---------------------------------------------------------------------
-c xi-direction flux differences
-c---------------------------------------------------------------------
-c
-c iex = flag : iex = 0 north/south communication
-c : iex = 1 east/west communication
-c
-c---------------------------------------------------------------------
- iex = 0
-
-c---------------------------------------------------------------------
-c communicate and receive/send two rows of data
-c---------------------------------------------------------------------
- call exchange_3 (rsd,iex)
-
- L1 = 0
- if (north.eq.-1) L1 = 1
- L2 = nx + 1
- if (south.eq.-1) L2 = nx
-
- ist1 = 1
- iend1 = nx
- if (north.eq.-1) ist1 = 4
- if (south.eq.-1) iend1 = nx - 3
-
- do k = 2, nz - 1
- do j = jst, jend
- do i = L1, L2
- flux(1,i,j,k) = 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,j,k) = rsd(2,i,j,k) * u21 + c2 *
- > ( rsd(5,i,j,k) - q )
- flux(3,i,j,k) = rsd(3,i,j,k) * u21
- flux(4,i,j,k) = rsd(4,i,j,k) * u21
- flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u21
- end do
- end do
- end do
-
- do k = 2, nz - 1
- do j = jst, jend
- do i = ist, iend
- do m = 1, 5
- frct(m,i,j,k) = frct(m,i,j,k)
- > - tx2 * ( flux(m,i+1,j,k) - flux(m,i-1,j,k) )
- end do
- end do
- do i = ist, L2
- 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,j,k) = (4.0d+00/3.0d+00) * tx3 *
- > ( u21i - u21im1 )
- flux(3,i,j,k) = tx3 * ( u31i - u31im1 )
- flux(4,i,j,k) = tx3 * ( u41i - u41im1 )
- flux(5,i,j,k) = 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,j,k) - flux(2,i,j,k) )
- > + 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,j,k) - flux(3,i,j,k) )
- > + 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,j,k) - flux(4,i,j,k) )
- > + 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,j,k) - flux(5,i,j,k) )
- > + 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---------------------------------------------------------------------
- IF (north.eq.-1) then
- do m = 1, 5
- frct(m,2,j,k) = frct(m,2,j,k)
- > - dsspm * ( + 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)
- > - dsspm * ( - 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
- END IF
-
- do i = ist1,iend1
- do m = 1, 5
- frct(m,i,j,k) = frct(m,i,j,k)
- > - dsspm * ( 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
-
- IF (south.eq.-1) then
- do m = 1, 5
- frct(m,nx-2,j,k) = frct(m,nx-2,j,k)
- > - dsspm * ( 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)
- > - dsspm * ( 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 IF
-
- end do
- end do
-
-c---------------------------------------------------------------------
-c eta-direction flux differences
-c---------------------------------------------------------------------
-c
-c iex = flag : iex = 0 north/south communication
-c : iex = 1 east/west communication
-c
-c---------------------------------------------------------------------
- iex = 1
-
-c---------------------------------------------------------------------
-c communicate and receive/send two rows of data
-c---------------------------------------------------------------------
- call exchange_3 (rsd,iex)
-
- L1 = 0
- if (west.eq.-1) L1 = 1
- L2 = ny + 1
- if (east.eq.-1) L2 = ny
-
- jst1 = 1
- jend1 = ny
- if (west.eq.-1) jst1 = 4
- if (east.eq.-1) jend1 = ny - 3
-
- do k = 2, nz - 1
- do j = L1, L2
- do i = ist, iend
- flux(1,i,j,k) = 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,i,j,k) = rsd(2,i,j,k) * u31
- flux(3,i,j,k) = rsd(3,i,j,k) * u31 + c2 *
- > ( rsd(5,i,j,k) - q )
- flux(4,i,j,k) = rsd(4,i,j,k) * u31
- flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u31
- end do
- end do
- end do
-
- do k = 2, nz - 1
- do i = ist, iend
- do j = jst, jend
- do m = 1, 5
- frct(m,i,j,k) = frct(m,i,j,k)
- > - ty2 * ( flux(m,i,j+1,k) - flux(m,i,j-1,k) )
- end do
- end do
- end do
-
- do j = jst, L2
- do i = ist, iend
- 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,i,j,k) = ty3 * ( u21j - u21jm1 )
- flux(3,i,j,k) = (4.0d+00/3.0d+00) * ty3 *
- > ( u31j - u31jm1 )
- flux(4,i,j,k) = ty3 * ( u41j - u41jm1 )
- flux(5,i,j,k) = 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
- end do
-
- do j = jst, jend
- do i = ist, iend
- 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,i,j+1,k) - flux(2,i,j,k) )
- > + 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,i,j+1,k) - flux(3,i,j,k) )
- > + 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,i,j+1,k) - flux(4,i,j,k) )
- > + 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,i,j+1,k) - flux(5,i,j,k) )
- > + dy5 * ty1 * ( rsd(5,i,j-1,k)
- > - 2.0d+00 * rsd(5,i,j,k)
- > + rsd(5,i,j+1,k) )
- end do
- end do
-
-c---------------------------------------------------------------------
-c fourth-order dissipation
-c---------------------------------------------------------------------
- IF (west.eq.-1) then
- do i = ist, iend
- do m = 1, 5
- frct(m,i,2,k) = frct(m,i,2,k)
- > - dsspm * ( + 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)
- > - dsspm * ( - 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
- end do
- END IF
-
- do j = jst1, jend1
- do i = ist, iend
- do m = 1, 5
- frct(m,i,j,k) = frct(m,i,j,k)
- > - dsspm * ( 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
- end do
-
- IF (east.eq.-1) then
- do i = ist, iend
- do m = 1, 5
- frct(m,i,ny-2,k) = frct(m,i,ny-2,k)
- > - dsspm * ( 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)
- > - dsspm * ( 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 IF
-
- end do
-
-c---------------------------------------------------------------------
-c zeta-direction flux differences
-c---------------------------------------------------------------------
- do k = 1, nz
- do j = jst, jend
- do i = ist, iend
- flux(1,i,j,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,i,j,k) = rsd(2,i,j,k) * u41
- flux(3,i,j,k) = rsd(3,i,j,k) * u41
- flux(4,i,j,k) = rsd(4,i,j,k) * u41 + c2 *
- > ( rsd(5,i,j,k) - q )
- flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u41
- end do
- end do
- end do
-
- do k = 2, nz - 1
- do j = jst, jend
- do i = ist, iend
- do m = 1, 5
- frct(m,i,j,k) = frct(m,i,j,k)
- > - tz2 * ( flux(m,i,j,k+1) - flux(m,i,j,k-1) )
- end do
- end do
- end do
- end do
-
- do k = 2, nz
- do j = jst, jend
- do i = ist, iend
- 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,i,j,k) = tz3 * ( u21k - u21km1 )
- flux(3,i,j,k) = tz3 * ( u31k - u31km1 )
- flux(4,i,j,k) = (4.0d+00/3.0d+00) * tz3 * ( u41k
- > - u41km1 )
- flux(5,i,j,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
- end do
- end do
-
- do k = 2, nz - 1
- do j = jst, jend
- do i = ist, iend
- 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,i,j,k+1) - flux(2,i,j,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,i,j,k+1) - flux(3,i,j,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,i,j,k+1) - flux(4,i,j,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,i,j,k+1) - flux(5,i,j,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
- end do
- end do
-
-c---------------------------------------------------------------------
-c fourth-order dissipation
-c---------------------------------------------------------------------
- do j = jst, jend
- do i = ist, iend
- do m = 1, 5
- frct(m,i,j,2) = frct(m,i,j,2)
- > - dsspm * ( + 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)
- > - dsspm * (- 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
- end do
- end do
-
- do k = 4, nz - 3
- do j = jst, jend
- do i = ist, iend
- do m = 1, 5
- frct(m,i,j,k) = frct(m,i,j,k)
- > - dsspm * ( 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
- end do
- end do
-
- do j = jst, jend
- do i = ist, iend
- do m = 1, 5
- frct(m,i,j,nz-2) = frct(m,i,j,nz-2)
- > - dsspm * ( 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)
- > - dsspm * ( 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
