c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine erhs(frct, rsd, nx, nxmax, ny, nz) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c c compute the right hand side based on exact solution c c--------------------------------------------------------------------- implicit none include 'header.h' integer nx, nxmax, ny, nz double precision frct(5,nxmax,ny,nz), rsd(5,nxmax,ny,nz) 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 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(u51km1,u41km1,u31km1,u21km1,u51k, !$OMP& u41k,u31k,u21k,u41,u51jm1,u41jm1,u31jm1,u21jm1,u51j,u41j,u31j, !$OMP& u21j,u31,u51im1,u41im1,u31im1,u21im1,u51i,u41i,u31i,u21i,tmp,q, !$OMP& u21,xi,eta,zeta,m,i,j,k) !$OMP& SHARED(dz5,dz4,dz3,dz2,tz1,dz1,tz3,tz2,dy5,dy4,dy3,dy2,ty1,dy1, !$OMP& ty3,ty2,dssp,dx5,dx4,dx3,dx2,tx1,dx1,tx3,tx2,nx,ny,nz) do k = 1, nz !$OMP DO SCHEDULE(STATIC) 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 !$OMP END DO nowait end do do k = 1, nz zeta = ( dble(k-1) ) / ( nz - 1 ) !$OMP DO SCHEDULE(STATIC) do j = 1, ny eta = ( dble(j-1) ) / ( ny - 1 ) do i = 1, nx xi = ( dble(i-1) ) / ( nx - 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 !$OMP END DO nowait end do !$OMP BARRIER c--------------------------------------------------------------------- c xi-direction flux differences c--------------------------------------------------------------------- do k = 2, nz - 1 !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 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 = 2, nx-1 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 = 2, 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 = 2, nx-1 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 !$OMP END DO nowait end do !$OMP BARRIER c--------------------------------------------------------------------- c eta-direction flux differences c--------------------------------------------------------------------- !$OMP DO SCHEDULE(STATIC) do k = 2, nz - 1 do i = 2, nx-1 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 = 2, ny-1 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 = 2, 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 = 2, ny-1 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 !$OMP END DO c--------------------------------------------------------------------- c zeta-direction flux differences c--------------------------------------------------------------------- !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do i = 2, nx-1 do k = 1, nz rtmp(1,k) = rsd(1,i,j,k) rtmp(2,k) = rsd(2,i,j,k) rtmp(3,k) = rsd(3,i,j,k) rtmp(4,k) = rsd(4,i,j,k) rtmp(5,k) = rsd(5,i,j,k) end do do k = 1, nz flux(1,k) = rtmp(4,k) u41 = rtmp(4,k) / rtmp(1,k) q = 0.50d+00 * ( rtmp(2,k) * rtmp(2,k) > + rtmp(3,k) * rtmp(3,k) > + rtmp(4,k) * rtmp(4,k) ) > / rtmp(1,k) flux(2,k) = rtmp(2,k) * u41 flux(3,k) = rtmp(3,k) * u41 flux(4,k) = rtmp(4,k) * u41 + c2 * > ( rtmp(5,k) - q ) flux(5,k) = ( c1 * rtmp(5,k) - c2 * q ) * u41 end do do k = 2, nz - 1 do m = 1, 5 utmp(m,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 / rtmp(1,k) u21k = tmp * rtmp(2,k) u31k = tmp * rtmp(3,k) u41k = tmp * rtmp(4,k) u51k = tmp * rtmp(5,k) tmp = 1.0d+00 / rtmp(1,k-1) u21km1 = tmp * rtmp(2,k-1) u31km1 = tmp * rtmp(3,k-1) u41km1 = tmp * rtmp(4,k-1) u51km1 = tmp * rtmp(5,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 utmp(1,k) = utmp(1,k) > + dz1 * tz1 * ( rtmp(1,k+1) > - 2.0d+00 * rtmp(1,k) > + rtmp(1,k-1) ) utmp(2,k) = utmp(2,k) > + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) ) > + dz2 * tz1 * ( rtmp(2,k+1) > - 2.0d+00 * rtmp(2,k) > + rtmp(2,k-1) ) utmp(3,k) = utmp(3,k) > + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) ) > + dz3 * tz1 * ( rtmp(3,k+1) > - 2.0d+00 * rtmp(3,k) > + rtmp(3,k-1) ) utmp(4,k) = utmp(4,k) > + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) ) > + dz4 * tz1 * ( rtmp(4,k+1) > - 2.0d+00 * rtmp(4,k) > + rtmp(4,k-1) ) utmp(5,k) = utmp(5,k) > + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) ) > + dz5 * tz1 * ( rtmp(5,k+1) > - 2.0d+00 * rtmp(5,k) > + rtmp(5,k-1) ) end do c--------------------------------------------------------------------- c fourth-order dissipation c--------------------------------------------------------------------- do m = 1, 5 frct(m,i,j,2) = utmp(m,2) > - dssp * ( + 5.0d+00 * rtmp(m,2) > - 4.0d+00 * rtmp(m,3) > + rtmp(m,4) ) frct(m,i,j,3) = utmp(m,3) > - dssp * (- 4.0d+00 * rtmp(m,2) > + 6.0d+00 * rtmp(m,3) > - 4.0d+00 * rtmp(m,4) > + rtmp(m,5) ) end do do k = 4, nz - 3 do m = 1, 5 frct(m,i,j,k) = utmp(m,k) > - dssp * ( rtmp(m,k-2) > - 4.0d+00 * rtmp(m,k-1) > + 6.0d+00 * rtmp(m,k) > - 4.0d+00 * rtmp(m,k+1) > + rtmp(m,k+2) ) end do end do do m = 1, 5 frct(m,i,j,nz-2) = utmp(m,nz-2) > - dssp * ( rtmp(m,nz-4) > - 4.0d+00 * rtmp(m,nz-3) > + 6.0d+00 * rtmp(m,nz-2) > - 4.0d+00 * rtmp(m,nz-1) ) frct(m,i,j,nz-1) = utmp(m,nz-1) > - dssp * ( rtmp(m,nz-3) > - 4.0d+00 * rtmp(m,nz-2) > + 5.0d+00 * rtmp(m,nz-1) ) end do end do end do !$OMP END DO nowait !$OMP END PARALLEL return end