c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine rhs(u, rsd, frct, qs, rho_i, nx, nxmax, ny, nz) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c compute the right hand sides c--------------------------------------------------------------------- implicit none include 'header.h' integer nx, nxmax, ny, nz double precision u(5,nxmax,ny,nz), rsd(5,nxmax,ny,nz), $ frct(5,nxmax,ny,nz), qs(nxmax,ny,nz), $ rho_i(nxmax,ny,nz) c--------------------------------------------------------------------- c local variables c--------------------------------------------------------------------- integer i, j, k, m double precision q double precision tmp double precision u21, u31, u41 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(u51im1,u41im1,u31im1,u21im1, !$OMP& u51i,u41i,u31i,u21i,tmp,m,q,u21,i,j,k, !$OMP& u51jm1,u41jm1,u31jm1,u21jm1,u51j,u41j,u31j,u21j,u31, !$OMP& u51km1,u41km1,u31km1,u21km1,u51k,u41k,u31k,u21k,u41) !$OMP& SHARED(dssp,dx5,dx4,dx3,dx2,tx1,dx1,tx3,tx2,nx,ny,nz, !$OMP& dy5,dy4,dy3,dy2,ty1,dy1,ty3,ty2, !$OMP& dz5,dz4,dz3,dz2,tz1,dz1,tz3,tz2) do k = 1, nz !$OMP DO SCHEDULE(STATIC) do j = 1, ny do i = 1, nx do m = 1, 5 rsd(m,i,j,k) = - frct(m,i,j,k) end do tmp = 1.0d+00 / u(1,i,j,k) rho_i(i,j,k) = tmp qs(i,j,k) = 0.50d+00 * ( u(2,i,j,k) * u(2,i,j,k) > + u(3,i,j,k) * u(3,i,j,k) > + u(4,i,j,k) * u(4,i,j,k) ) > * tmp end do end do !$OMP END DO nowait end do !$OMP BARRIER !$OMP MASTER if (timeron) call timer_start(t_rhsx) !$OMP END MASTER 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) = u(2,i,j,k) u21 = u(2,i,j,k) * rho_i(i,j,k) q = qs(i,j,k) flux(2,i) = u(2,i,j,k) * u21 + c2 * > ( u(5,i,j,k) - q ) flux(3,i) = u(3,i,j,k) * u21 flux(4,i) = u(4,i,j,k) * u21 flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21 end do do i = 2, nx-1 do m = 1, 5 rsd(m,i,j,k) = rsd(m,i,j,k) > - tx2 * ( flux(m,i+1) - flux(m,i-1) ) end do end do do i = 2, nx tmp = rho_i(i,j,k) u21i = tmp * u(2,i,j,k) u31i = tmp * u(3,i,j,k) u41i = tmp * u(4,i,j,k) u51i = tmp * u(5,i,j,k) tmp = rho_i(i-1,j,k) u21im1 = tmp * u(2,i-1,j,k) u31im1 = tmp * u(3,i-1,j,k) u41im1 = tmp * u(4,i-1,j,k) u51im1 = tmp * u(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 rsd(1,i,j,k) = rsd(1,i,j,k) > + dx1 * tx1 * ( u(1,i-1,j,k) > - 2.0d+00 * u(1,i,j,k) > + u(1,i+1,j,k) ) rsd(2,i,j,k) = rsd(2,i,j,k) > + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) ) > + dx2 * tx1 * ( u(2,i-1,j,k) > - 2.0d+00 * u(2,i,j,k) > + u(2,i+1,j,k) ) rsd(3,i,j,k) = rsd(3,i,j,k) > + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) ) > + dx3 * tx1 * ( u(3,i-1,j,k) > - 2.0d+00 * u(3,i,j,k) > + u(3,i+1,j,k) ) rsd(4,i,j,k) = rsd(4,i,j,k) > + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) ) > + dx4 * tx1 * ( u(4,i-1,j,k) > - 2.0d+00 * u(4,i,j,k) > + u(4,i+1,j,k) ) rsd(5,i,j,k) = rsd(5,i,j,k) > + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) ) > + dx5 * tx1 * ( u(5,i-1,j,k) > - 2.0d+00 * u(5,i,j,k) > + u(5,i+1,j,k) ) end do end do !$OMP END DO nowait c--------------------------------------------------------------------- c Fourth-order dissipation c--------------------------------------------------------------------- !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do m = 1, 5 rsd(m,2,j,k) = rsd(m,2,j,k) > - dssp * ( + 5.0d+00 * u(m,2,j,k) > - 4.0d+00 * u(m,3,j,k) > + u(m,4,j,k) ) rsd(m,3,j,k) = rsd(m,3,j,k) > - dssp * ( - 4.0d+00 * u(m,2,j,k) > + 6.0d+00 * u(m,3,j,k) > - 4.0d+00 * u(m,4,j,k) > + u(m,5,j,k) ) end do end do !$OMP END DO nowait !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do i = 4, nx - 3 do m = 1, 5 rsd(m,i,j,k) = rsd(m,i,j,k) > - dssp * ( u(m,i-2,j,k) > - 4.0d+00 * u(m,i-1,j,k) > + 6.0d+00 * u(m,i,j,k) > - 4.0d+00 * u(m,i+1,j,k) > + u(m,i+2,j,k) ) end do end do end do !$OMP END DO nowait !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do m = 1, 5 rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k) > - dssp * ( u(m,nx-4,j,k) > - 4.0d+00 * u(m,nx-3,j,k) > + 6.0d+00 * u(m,nx-2,j,k) > - 4.0d+00 * u(m,nx-1,j,k) ) rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k) > - dssp * ( u(m,nx-3,j,k) > - 4.0d+00 * u(m,nx-2,j,k) > + 5.0d+00 * u(m,nx-1,j,k) ) end do end do !$OMP END DO nowait end do !$OMP BARRIER !$OMP MASTER if (timeron) call timer_stop(t_rhsx) if (timeron) call timer_start(t_rhsy) !$OMP END MASTER 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) = u(3,i,j,k) u31 = u(3,i,j,k) * rho_i(i,j,k) q = qs(i,j,k) flux(2,j) = u(2,i,j,k) * u31 flux(3,j) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q) flux(4,j) = u(4,i,j,k) * u31 flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31 end do do j = 2, ny-1 do m = 1, 5 rsd(m,i,j,k) = rsd(m,i,j,k) > - ty2 * ( flux(m,j+1) - flux(m,j-1) ) end do end do do j = 2, ny tmp = rho_i(i,j,k) u21j = tmp * u(2,i,j,k) u31j = tmp * u(3,i,j,k) u41j = tmp * u(4,i,j,k) u51j = tmp * u(5,i,j,k) tmp = rho_i(i,j-1,k) u21jm1 = tmp * u(2,i,j-1,k) u31jm1 = tmp * u(3,i,j-1,k) u41jm1 = tmp * u(4,i,j-1,k) u51jm1 = tmp * u(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 rsd(1,i,j,k) = rsd(1,i,j,k) > + dy1 * ty1 * ( u(1,i,j-1,k) > - 2.0d+00 * u(1,i,j,k) > + u(1,i,j+1,k) ) rsd(2,i,j,k) = rsd(2,i,j,k) > + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) ) > + dy2 * ty1 * ( u(2,i,j-1,k) > - 2.0d+00 * u(2,i,j,k) > + u(2,i,j+1,k) ) rsd(3,i,j,k) = rsd(3,i,j,k) > + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) ) > + dy3 * ty1 * ( u(3,i,j-1,k) > - 2.0d+00 * u(3,i,j,k) > + u(3,i,j+1,k) ) rsd(4,i,j,k) = rsd(4,i,j,k) > + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) ) > + dy4 * ty1 * ( u(4,i,j-1,k) > - 2.0d+00 * u(4,i,j,k) > + u(4,i,j+1,k) ) rsd(5,i,j,k) = rsd(5,i,j,k) > + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) ) > + dy5 * ty1 * ( u(5,i,j-1,k) > - 2.0d+00 * u(5,i,j,k) > + u(5,i,j+1,k) ) end do end do c--------------------------------------------------------------------- c fourth-order dissipation c--------------------------------------------------------------------- do i = 2, nx-1 do m = 1, 5 rsd(m,i,2,k) = rsd(m,i,2,k) > - dssp * ( + 5.0d+00 * u(m,i,2,k) > - 4.0d+00 * u(m,i,3,k) > + u(m,i,4,k) ) rsd(m,i,3,k) = rsd(m,i,3,k) > - dssp * ( - 4.0d+00 * u(m,i,2,k) > + 6.0d+00 * u(m,i,3,k) > - 4.0d+00 * u(m,i,4,k) > + u(m,i,5,k) ) end do end do do j = 4, ny - 3 do i = 2, nx-1 do m = 1, 5 rsd(m,i,j,k) = rsd(m,i,j,k) > - dssp * ( u(m,i,j-2,k) > - 4.0d+00 * u(m,i,j-1,k) > + 6.0d+00 * u(m,i,j,k) > - 4.0d+00 * u(m,i,j+1,k) > + u(m,i,j+2,k) ) end do end do end do do i = 2, nx-1 do m = 1, 5 rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k) > - dssp * ( u(m,i,ny-4,k) > - 4.0d+00 * u(m,i,ny-3,k) > + 6.0d+00 * u(m,i,ny-2,k) > - 4.0d+00 * u(m,i,ny-1,k) ) rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k) > - dssp * ( u(m,i,ny-3,k) > - 4.0d+00 * u(m,i,ny-2,k) > + 5.0d+00 * u(m,i,ny-1,k) ) end do end do end do !$OMP END DO !$OMP MASTER if (timeron) call timer_stop(t_rhsy) if (timeron) call timer_start(t_rhsz) !$OMP END MASTER 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 flux(1,k) = u(4,i,j,k) u41 = u(4,i,j,k) * rho_i(i,j,k) q = qs(i,j,k) flux(2,k) = u(2,i,j,k) * u41 flux(3,k) = u(3,i,j,k) * u41 flux(4,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q) flux(5,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41 end do do k = 2, nz - 1 do m = 1, 5 rsd(m,i,j,k) = rsd(m,i,j,k) > - tz2 * ( flux(m,k+1) - flux(m,k-1) ) end do end do do k = 2, nz tmp = rho_i(i,j,k) u21k = tmp * u(2,i,j,k) u31k = tmp * u(3,i,j,k) u41k = tmp * u(4,i,j,k) u51k = tmp * u(5,i,j,k) tmp = rho_i(i,j,k-1) u21km1 = tmp * u(2,i,j,k-1) u31km1 = tmp * u(3,i,j,k-1) u41km1 = tmp * u(4,i,j,k-1) u51km1 = tmp * u(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 rsd(1,i,j,k) = rsd(1,i,j,k) > + dz1 * tz1 * ( u(1,i,j,k-1) > - 2.0d+00 * u(1,i,j,k) > + u(1,i,j,k+1) ) rsd(2,i,j,k) = rsd(2,i,j,k) > + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) ) > + dz2 * tz1 * ( u(2,i,j,k-1) > - 2.0d+00 * u(2,i,j,k) > + u(2,i,j,k+1) ) rsd(3,i,j,k) = rsd(3,i,j,k) > + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) ) > + dz3 * tz1 * ( u(3,i,j,k-1) > - 2.0d+00 * u(3,i,j,k) > + u(3,i,j,k+1) ) rsd(4,i,j,k) = rsd(4,i,j,k) > + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) ) > + dz4 * tz1 * ( u(4,i,j,k-1) > - 2.0d+00 * u(4,i,j,k) > + u(4,i,j,k+1) ) rsd(5,i,j,k) = rsd(5,i,j,k) > + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) ) > + dz5 * tz1 * ( u(5,i,j,k-1) > - 2.0d+00 * u(5,i,j,k) > + u(5,i,j,k+1) ) end do end do end do !$OMP END DO nowait c--------------------------------------------------------------------- c fourth-order dissipation c--------------------------------------------------------------------- !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do i = 2, nx-1 do m = 1, 5 rsd(m,i,j,2) = rsd(m,i,j,2) > - dssp * ( + 5.0d+00 * u(m,i,j,2) > - 4.0d+00 * u(m,i,j,3) > + u(m,i,j,4) ) rsd(m,i,j,3) = rsd(m,i,j,3) > - dssp * ( - 4.0d+00 * u(m,i,j,2) > + 6.0d+00 * u(m,i,j,3) > - 4.0d+00 * u(m,i,j,4) > + u(m,i,j,5) ) end do end do end do !$OMP END DO nowait do k = 4, nz - 3 !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do i = 2, nx-1 do m = 1, 5 rsd(m,i,j,k) = rsd(m,i,j,k) > - dssp * ( u(m,i,j,k-2) > - 4.0d+00 * u(m,i,j,k-1) > + 6.0d+00 * u(m,i,j,k) > - 4.0d+00 * u(m,i,j,k+1) > + u(m,i,j,k+2) ) end do end do end do !$OMP END DO nowait end do !$OMP DO SCHEDULE(STATIC) do j = 2, ny-1 do i = 2, nx-1 do m = 1, 5 rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2) > - dssp * ( u(m,i,j,nz-4) > - 4.0d+00 * u(m,i,j,nz-3) > + 6.0d+00 * u(m,i,j,nz-2) > - 4.0d+00 * u(m,i,j,nz-1) ) rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1) > - dssp * ( u(m,i,j,nz-3) > - 4.0d+00 * u(m,i,j,nz-2) > + 5.0d+00 * u(m,i,j,nz-1) ) end do end do end do !$OMP END DO nowait !$OMP END PARALLEL if (timeron) call timer_stop(t_rhsz) return end