UpdateConsv Subroutine

      subroutine UpdateConsv
      use basicmod
      use fluxmod
      implicit none
      integer::i,j,k

!$acc kernels
!$acc loop collapse(3) independent
      do k=ks,ke
      do j=js,je
      do i=is,ie
         
         d(i,j,k) = d(i,j,k)  &
     & +dt*(                  &
     &  (- nflux1(mden,i+1,j,k) &
     &   + nflux1(mden,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mden,i,j+1,k) &
     &   + nflux2(mden,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(mden,i,j,k+1) &
     &   + nflux3(mden,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

         mv1(i,j,k) = mv1(i,j,k) &
     & +dt*( &
     &      +  grvsrc1(i,j,k)                      &
     & +(- nflux1(mrv1,i+1,j,k) &
     &   + nflux1(mrv1,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mrv1,i,j+1,k) &
     &   + nflux2(mrv1,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(mrv1,i,j,k+1) &
     &   + nflux3(mrv1,i,j,k  ))/(x3a(k+1)-x3a(k)) & 
     &      )
         mv2(i,j,k) = mv2(i,j,k) &
     & +dt*( &
     &      +  grvsrc2(i,j,k)                      &
     & +(- nflux1(mrv2,i+1,j,k) &
     &   + nflux1(mrv2,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mrv2,i,j+1,k) &
     &   + nflux2(mrv2,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(mrv2,i,j,k+1) &
     &   + nflux3(mrv2,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

         mv3(i,j,k) = mv3(i,j,k)  &
     & +dt*( &
     &      +  grvsrc3(i,j,k)                      &
     & +(- nflux1(mrv3,i+1,j,k) &
     &   + nflux1(mrv3,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mrv3,i,j+1,k) &
     &   + nflux2(mrv3,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(mrv3,i,j,k+1) &
     &   + nflux3(mrv3,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

          et(i,j,k) = et(i,j,k) &
     & +dt*( &
     &  (- nflux1(meto,i+1,j,k) &
     &   + nflux1(meto,i  ,j,k))/(x1a(i+1)-x1a(i)) &
     & +(- nflux2(meto,i,j+1,k) &
     &   + nflux2(meto,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(meto,i,j,k+1) &
     &   + nflux3(meto,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

          b1(i,j,k) = b1(i,j,k) &
     & +dt*( &
     &  (- nflux1(mbm1,i+1,j,k) &
     &   + nflux1(mbm1,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mbm1,i,j+1,k) &
     &   + nflux2(mbm1,i,j  ,k))/(x2a(j+1)-x2a(j)) &  
     & +(- nflux3(mbm1,i,j,k+1) &
     &   + nflux3(mbm1,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

          b2(i,j,k) = b2(i,j,k) &
     & +dt*( &
     &  (- nflux1(mbm2,i+1,j,k) &
     &   + nflux1(mbm2,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mbm2,i,j+1,k) &
     &   + nflux2(mbm2,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(mbm2,i,j,k+1) &
     &   + nflux3(mbm2,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

          b3(i,j,k) = b3(i,j,k) &
     & +dt*( &
     &  (- nflux1(mbm3,i+1,j,k) &
     &   + nflux1(mbm3,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mbm3,i,j+1,k) &
     &   + nflux2(mbm3,i,j  ,k))/(x2a(j+1)-x2a(j)) &  
     & +(- nflux3(mbm3,i,j,k+1) &
     &   + nflux3(mbm3,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

          bp(i,j,k) = bp(i,j,k) &
     & +dt*( &
     &  (- nflux1(mbps,i+1,j,k) &
     &   + nflux1(mbps,i  ,j,k))/(x1a(i+1)-x1a(i)) & 
     & +(- nflux2(mbps,i,j+1,k) &
     &   + nflux2(mbps,i,j  ,k))/(x2a(j+1)-x2a(j)) & 
     & +(- nflux3(mbps,i,j,k+1) &
     &   + nflux3(mbps,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )

!          print *, i,j,k,bp(i,j,k)
         DXcomp(1:ncomp,i,j,k) = DXcomp(1:ncomp,i,j,k) &
     & +dt*( &
     & +(- nflux1(mst:med,i+1,j,k) &
     &   + nflux1(mst:med,i  ,j,k))/(x1a(i+1)-x1a(i)) &
     & +(- nflux2(mst:med,i,j+1,k) &
     &   + nflux2(mst:med,i,j  ,k))/(x2a(j+1)-x2a(j)) &
     & +(- nflux3(mst:med,i,j,k+1) &
     &   + nflux3(mst:med,i,j,k  ))/(x3a(k+1)-x3a(k)) &
     &      )
      enddo
      enddo
      enddo
!$acc end kernels

      return
    end subroutine UpdateConsv