! *** generated by SAPFOR with version 2415 and build date: May  4 2025 14:48:40
! *** Enabled options ***:
! *** maximum shadow width is 50 percent
! *** generated by SAPFOR

!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine compute_rhs_sp (aditional_comp)

       include 'header_sp.h'
      integer ::  i,j,k,m
      double precision ::  aux,rho_inv,uijk,up1,um1,vijk,vp1,vm1,wijk,wp
     &1,wm1,rhs_(5)
      double precision ::  t1,t2,t3,ac,ru1,uu,vv,ww,ac2inv
      integer ::  aditional_comp
      if (timeron)  call timer_start_sp(t_rhs)

!$SPF PARALLEL_REG r0
! DVM$ PARALLEL (K,J,I) ON U(*,I,J,K), PRIVATE (RHO_INV,AUX,M),SHADOW_REN
! DVM$&EW (U(0:0,2:3,2:3,2:3)),SHADOW_COMPUTE 
! DVM$ REGION OUT (US,VS,WS,QS,RHO_I,SPEED,SQUARE)
      do  k = 0,problem_size - 1
         do  j = 0,problem_size - 1
            do  i = 0,problem_size - 1
               rho_inv = 1.0d0 / u(1,i,j,k)
               rho_i(i,j,k) = rho_inv
               us(i,j,k) = u(2,i,j,k) * rho_inv
               vs(i,j,k) = u(3,i,j,k) * rho_inv
               ws(i,j,k) = u(4,i,j,k) * rho_inv
               square(i,j,k) = 0.5d0 * (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)) * rho_inv
               qs(i,j,k) = square(i,j,k) * rho_inv

!---------------------------------------------------------------------
!               (don't need speed and ainx until the lhs computation)
!---------------------------------------------------------------------
               aux = c1c2 * rho_inv * (u(5,i,j,k) - square(i,j,k))
               speed(i,j,k) = dsqrt (aux)
               do  m = 1,5
                  rhs(m,i,j,k) = forcing(m,i,j,k)
               enddo  
            enddo  
         enddo  
      enddo  

!$SPF ANALYSIS(PRIVATE(rhs_))
! DVM$ PARALLEL (K,J,I) ON RHS(*,I,J,K), PRIVATE (UIJK,UP1,UM1,M,VIJK,VP1
! DVM$&,VM1,WIJK,WP1,WM1,RHS_,T1,T2,T3,AC,RU1,UU,VV,WW,AC2INV),CUDA_BLOCK
! DVM$& (32,4)
      do  k = 1,nz2
         do  j = 1,ny2
            do  i = 1,nx2
               uijk = us(i,j,k)
               up1 = us(i + 1,j,k)
               um1 = us(i - 1,j,k)
               rhs_(1) = rhs(1,i,j,k)
               rhs_(2) = rhs(2,i,j,k)
               rhs_(3) = rhs(3,i,j,k)
               rhs_(4) = rhs(4,i,j,k)
               rhs_(5) = rhs(5,i,j,k)
               rhs_(1) = rhs_(1) + dx1tx1 * (u(1,i + 1,j,k) - 2.0d0 * u(
     &1,i,j,k) + u(1,i - 1,j,k)) - tx2 * (u(2,i + 1,j,k) - u(2,i - 1,j,k
     &))
               rhs_(2) = rhs_(2) + dx2tx1 * (u(2,i + 1,j,k) - 2.0d0 * u(
     &2,i,j,k) + u(2,i - 1,j,k)) + xxcon2 * con43 * (up1 - 2.0d0 * uijk 
     &+ um1) - tx2 * (u(2,i + 1,j,k) * up1 - u(2,i - 1,j,k) * um1 + (u(5
     &,i + 1,j,k) - square(i + 1,j,k) - u(5,i - 1,j,k) + square(i - 1,j,
     &k)) * c2)
               rhs_(3) = rhs_(3) + dx3tx1 * (u(3,i + 1,j,k) - 2.0d0 * u(
     &3,i,j,k) + u(3,i - 1,j,k)) + xxcon2 * (vs(i + 1,j,k) - 2.0d0 * vs(
     &i,j,k) + vs(i - 1,j,k)) - tx2 * (u(3,i + 1,j,k) * up1 - u(3,i - 1,
     &j,k) * um1)
               rhs_(4) = rhs_(4) + dx4tx1 * (u(4,i + 1,j,k) - 2.0d0 * u(
     &4,i,j,k) + u(4,i - 1,j,k)) + xxcon2 * (ws(i + 1,j,k) - 2.0d0 * ws(
     &i,j,k) + ws(i - 1,j,k)) - tx2 * (u(4,i + 1,j,k) * up1 - u(4,i - 1,
     &j,k) * um1)
               rhs_(5) = rhs_(5) + dx5tx1 * (u(5,i + 1,j,k) - 2.0d0 * u(
     &5,i,j,k) + u(5,i - 1,j,k)) + xxcon3 * (qs(i + 1,j,k) - 2.0d0 * qs(
     &i,j,k) + qs(i - 1,j,k)) + xxcon4 * (up1 * up1 - 2.0d0 * uijk * uij
     &k + um1 * um1) + xxcon5 * (u(5,i + 1,j,k) * rho_i(i + 1,j,k) - 2.0
     &d0 * u(5,i,j,k) * rho_i(i,j,k) + u(5,i - 1,j,k) * rho_i(i - 1,j,k)
     &) - tx2 * ((c1 * u(5,i + 1,j,k) - c2 * square(i + 1,j,k)) * up1 - 
     &(c1 * u(5,i - 1,j,k) - c2 * square(i - 1,j,k)) * um1)
               if (i .eq. 1) then
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (5.0d0 * u(m,i,j,k) - 4.
     &0d0 * u(m,i + 1,j,k) + u(m,i + 2,j,k))
                  enddo  
               else if (i .eq. 2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * ((-(4.0d0)) * u(m,i - 1,
     &j,k) + 6.0d0 * u(m,i,j,k) - 4.0d0 * u(m,i + 1,j,k) + u(m,i + 2,j,k
     &))
                  enddo  
               else if (i .ge. 3 .and. i .le. nx2 - 2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i - 2,j,k) - 4.0d0 
     &* u(m,i - 1,j,k) + 6.0 * u(m,i,j,k) - 4.0d0 * u(m,i + 1,j,k) + u(m
     &,i + 2,j,k))
                  enddo  
               else if (i .eq. nx2 - 1) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i - 2,j,k) - 4.0d0 
     &* u(m,i - 1,j,k) + 6.0d0 * u(m,i,j,k) - 4.0d0 * u(m,i + 1,j,k))
                  enddo  
               else if (i .eq. nx2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i - 2,j,k) - 4.d0 *
     & u(m,i - 1,j,k) + 5.d0 * u(m,i,j,k))
                  enddo  
               endif  
               vijk = vs(i,j,k)
               vp1 = vs(i,j + 1,k)
               vm1 = vs(i,j - 1,k)
               rhs_(1) = rhs_(1) + dy1ty1 * (u(1,i,j + 1,k) - 2.0d0 * u(
     &1,i,j,k) + u(1,i,j - 1,k)) - ty2 * (u(3,i,j + 1,k) - u(3,i,j - 1,k
     &))
               rhs_(2) = rhs_(2) + dy2ty1 * (u(2,i,j + 1,k) - 2.0d0 * u(
     &2,i,j,k) + u(2,i,j - 1,k)) + yycon2 * (us(i,j + 1,k) - 2.0d0 * us(
     &i,j,k) + us(i,j - 1,k)) - ty2 * (u(2,i,j + 1,k) * vp1 - u(2,i,j - 
     &1,k) * vm1)
               rhs_(3) = rhs_(3) + dy3ty1 * (u(3,i,j + 1,k) - 2.0d0 * u(
     &3,i,j,k) + u(3,i,j - 1,k)) + yycon2 * con43 * (vp1 - 2.0d0 * vijk 
     &+ vm1) - ty2 * (u(3,i,j + 1,k) * vp1 - u(3,i,j - 1,k) * vm1 + (u(5
     &,i,j + 1,k) - square(i,j + 1,k) - u(5,i,j - 1,k) + square(i,j - 1,
     &k)) * c2)
               rhs_(4) = rhs_(4) + dy4ty1 * (u(4,i,j + 1,k) - 2.0d0 * u(
     &4,i,j,k) + u(4,i,j - 1,k)) + yycon2 * (ws(i,j + 1,k) - 2.0d0 * ws(
     &i,j,k) + ws(i,j - 1,k)) - ty2 * (u(4,i,j + 1,k) * vp1 - u(4,i,j - 
     &1,k) * vm1)
               rhs_(5) = rhs_(5) + dy5ty1 * (u(5,i,j + 1,k) - 2.0d0 * u(
     &5,i,j,k) + u(5,i,j - 1,k)) + yycon3 * (qs(i,j + 1,k) - 2.0d0 * qs(
     &i,j,k) + qs(i,j - 1,k)) + yycon4 * (vp1 * vp1 - 2.0d0 * vijk * vij
     &k + vm1 * vm1) + yycon5 * (u(5,i,j + 1,k) * rho_i(i,j + 1,k) - 2.0
     &d0 * u(5,i,j,k) * rho_i(i,j,k) + u(5,i,j - 1,k) * rho_i(i,j - 1,k)
     &) - ty2 * ((c1 * u(5,i,j + 1,k) - c2 * square(i,j + 1,k)) * vp1 - 
     &(c1 * u(5,i,j - 1,k) - c2 * square(i,j - 1,k)) * vm1)
               if (j .eq. 1) then
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (5.0d0 * u(m,i,j,k) - 4.
     &0d0 * u(m,i,j + 1,k) + u(m,i,j + 2,k))
                  enddo  
               else if (j .eq. 2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * ((-(4.0d0)) * u(m,i,j - 
     &1,k) + 6.0d0 * u(m,i,j,k) - 4.0d0 * u(m,i,j + 1,k) + u(m,i,j + 2,k
     &))
                  enddo  
               else if (j .ge. 3 .and. j .le. ny2 - 2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i,j - 2,k) - 4.0d0 
     &* u(m,i,j - 1,k) + 6.0 * u(m,i,j,k) - 4.0d0 * u(m,i,j + 1,k) + u(m
     &,i,j + 2,k))
                  enddo  
               else if (j .eq. ny2 - 1) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i,j - 2,k) - 4.0d0 
     &* u(m,i,j - 1,k) + 6.0d0 * u(m,i,j,k) - 4.0d0 * u(m,i,j + 1,k))
                  enddo  
               else if (j .eq. ny2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i,j - 2,k) - 4.d0 *
     & u(m,i,j - 1,k) + 5.d0 * u(m,i,j,k))
                  enddo  
               endif  
               wijk = ws(i,j,k)
               wp1 = ws(i,j,k + 1)
               wm1 = ws(i,j,k - 1)
               rhs_(1) = rhs_(1) + dz1tz1 * (u(1,i,j,k + 1) - 2.0d0 * u(
     &1,i,j,k) + u(1,i,j,k - 1)) - tz2 * (u(4,i,j,k + 1) - u(4,i,j,k - 1
     &))
               rhs_(2) = rhs_(2) + dz2tz1 * (u(2,i,j,k + 1) - 2.0d0 * u(
     &2,i,j,k) + u(2,i,j,k - 1)) + zzcon2 * (us(i,j,k + 1) - 2.0d0 * us(
     &i,j,k) + us(i,j,k - 1)) - tz2 * (u(2,i,j,k + 1) * wp1 - u(2,i,j,k 
     &- 1) * wm1)
               rhs_(3) = rhs_(3) + dz3tz1 * (u(3,i,j,k + 1) - 2.0d0 * u(
     &3,i,j,k) + u(3,i,j,k - 1)) + zzcon2 * (vs(i,j,k + 1) - 2.0d0 * vs(
     &i,j,k) + vs(i,j,k - 1)) - tz2 * (u(3,i,j,k + 1) * wp1 - u(3,i,j,k 
     &- 1) * wm1)
               rhs_(4) = rhs_(4) + dz4tz1 * (u(4,i,j,k + 1) - 2.0d0 * u(
     &4,i,j,k) + u(4,i,j,k - 1)) + zzcon2 * con43 * (wp1 - 2.0d0 * wijk 
     &+ wm1) - tz2 * (u(4,i,j,k + 1) * wp1 - u(4,i,j,k - 1) * wm1 + (u(5
     &,i,j,k + 1) - square(i,j,k + 1) - u(5,i,j,k - 1) + square(i,j,k - 
     &1)) * c2)
               rhs_(5) = rhs_(5) + dz5tz1 * (u(5,i,j,k + 1) - 2.0d0 * u(
     &5,i,j,k) + u(5,i,j,k - 1)) + zzcon3 * (qs(i,j,k + 1) - 2.0d0 * qs(
     &i,j,k) + qs(i,j,k - 1)) + zzcon4 * (wp1 * wp1 - 2.0d0 * wijk * wij
     &k + wm1 * wm1) + zzcon5 * (u(5,i,j,k + 1) * rho_i(i,j,k + 1) - 2.0
     &d0 * u(5,i,j,k) * rho_i(i,j,k) + u(5,i,j,k - 1) * rho_i(i,j,k - 1)
     &) - tz2 * ((c1 * u(5,i,j,k + 1) - c2 * square(i,j,k + 1)) * wp1 - 
     &(c1 * u(5,i,j,k - 1) - c2 * square(i,j,k - 1)) * wm1)
               if (k .eq. 1) then
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (5.0d0 * u(m,i,j,k) - 4.
     &0d0 * u(m,i,j,k + 1) + u(m,i,j,k + 2))
                  enddo  
               else if (k .eq. 2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * ((-(4.0d0)) * u(m,i,j,k 
     &- 1) + 6.0d0 * u(m,i,j,k) - 4.0d0 * u(m,i,j,k + 1) + u(m,i,j,k + 2
     &))
                  enddo  
               else if (k .ge. 3 .and. k .le. nz2 - 2) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i,j,k - 2) - 4.0d0 
     &* u(m,i,j,k - 1) + 6.0 * u(m,i,j,k) - 4.0d0 * u(m,i,j,k + 1) + u(m
     &,i,j,k + 2))
                  enddo  
               else if (k .eq. nz2 - 1) then  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i,j,k - 2) - 4.0d0 
     &* u(m,i,j,k - 1) + 6.0d0 * u(m,i,j,k) - 4.0d0 * u(m,i,j,k + 1))
                  enddo  
               else  
                  do  m = 1,5
                     rhs_(m) = rhs_(m) - dssp * (u(m,i,j,k - 2) - 4.d0 *
     & u(m,i,j,k - 1) + 5.d0 * u(m,i,j,k))
                  enddo  
               endif  
               rhs_(1) = rhs_(1) * dt
               rhs_(2) = rhs_(2) * dt
               rhs_(3) = rhs_(3) * dt
               rhs_(4) = rhs_(4) * dt
               rhs_(5) = rhs_(5) * dt
               rhs(1,i,j,k) = rhs_(1)
               rhs(2,i,j,k) = rhs_(2)
               rhs(3,i,j,k) = rhs_(3)
               rhs(4,i,j,k) = rhs_(4)
               rhs(5,i,j,k) = rhs_(5)
               if (aditional_comp .eq. 1) then
                  ru1 = rho_i(i,j,k)
                  uu = us(i,j,k)
                  vv = vs(i,j,k)
                  ww = ws(i,j,k)
                  ac = speed(i,j,k)
                  ac2inv = ac * ac
                  t1 = c2 / ac2inv * (qs(i,j,k) * rhs_(1) - uu * rhs_(2)
     & - vv * rhs_(3) - ww * rhs_(4) + rhs_(5))
                  t2 = bt * ru1 * (uu * rhs_(1) - rhs_(2))
                  t3 = bt * ru1 * ac * t1
                  rhs(1,i,j,k) = rhs_(1) - t1
                  rhs(2,i,j,k) = (-(ru1)) * (ww * rhs_(1) - rhs_(4))
                  rhs(3,i,j,k) = ru1 * (vv * rhs_(1) - rhs_(3))
                  rhs(4,i,j,k) = (-(t2)) + t3
                  rhs(5,i,j,k) = t2 + t3
               endif  
            enddo  
         enddo  
      enddo  
!$SPF END PARALLEL_REG
! DVM$ END REGION
      if (timeron)  call timer_stop_sp(t_rhs)
      return 
      end