      SUBROUTINE correc
       
      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'common.txt'


      real flux,flux1,flux2,flux3,sum


c****************************************************************
c
c     Subroutine correc corrects the predicted velocities dUdt, dVdt
c     and dWdt in order to satisfy continuity. All other variables
c     (i.e. all scalars) do not need this correction, which means that
c     the predicted values (subroutine predic) are O.K.
c
c     The predicted velocities, normal to the wall, at the boundary
c     r=R are known (set at end of subroutine predic):
c
c                dUdt(imax,j,k) = 0.0    (r=R).
c
c     Therefore we do not need pressure points outside of the domain
c     for the calculation of the U-velocity on the boundary.
c     For the V- and W-velocity we have to take periodic boundary
c     conditions for the pressure in theta- and z-direction.
c     This means:
c
c                p(i,j1,k) = p(i,1,k)
c                p(i,j,k1) = p(i,j,1)
c
c****************************************************************
c
c
c ***************** CORRECTION  U-VELOCITY **********************


      DO k=1,kmax
        DO j=1,jmax
           DO i=1,imax-1

              umnew(i,j,k) = umnew(i,j,k) -
     +        dt * ( p(i+1,j,k) - p(i,j,k) ) /
     +                ( Rp(i+1) - Rp(i) )

              unew(i,j,k)=umnew(i,j,k)/(0.5*(rhonew(i,j,k)+rhonew(i+1,j
     $             ,k)))


           END DO
        END DO
      END DO


      if(1.eq.1)then

      DO k=1,kmax
         DO j=1,jmax

            umnew(imax,j,k)=umnew(imax,j,k)-dt*(-2*P(imax,j,k))
     $             /(Rp(imax+1)-Rp(imax))

            unew(imax,j,k)=umnew(imax,j,k)/(0.5*(rhonew(imax,j,k)
     $             +rhonew(i1,j,k)))
            
         END DO
      END DO
        
      end if

c ******************** END U-VELOCITY ***************************
c
c
c ******************** CORRECTION  V-VELOCITY *******************


      DO k=1,kmax
         DO j=1,jmax-1
           DO i=1,imax

              vmnew(i,j,k) = vmnew(i,j,k) -
     +        dt * ( p(i,j+1,k) - p(i,j,k) ) /
     +                ( Rp(i) * dphi )

              vnew(i,j,k)=vmnew(i,j,k)/(0.5*(rhonew(i,j,k)+rhonew(i,j+1
     $             ,k)))

           END DO
        END DO
      END DO

c *** Treat j=jmax as a special case: ***


      DO k=1,kmax
        DO i=1,imax

           vmnew(i,jmax,k) = vmnew(i,jmax,k) -
     +       dt * ( p(i,1,k) - p(i,jmax,k) ) /
     +              ( Rp(i) * dphi )

           vnew(i,jmax,k)=vmnew(i,jmax,k)/(0.5*(rhonew(i,jmax,k)
     $          +rhonew(i,j1,k)))

        END DO
      END DO


c ******************** END V-VELOCITY **************************
c
c
c ***************** CORRECTION  W-VELOCITY *********************


      DO k=1,kmax-1
         DO j=1,jmax
            DO i=1,imax
             
               wmnew(i,j,k) = wmnew(i,j,k) -
     +         dt * ( p(i,j,k+1) - p(i,j,k) ) / (zp(k+1)-zp(k))

               wnew(i,j,k)=wmnew(i,j,k)/(0.5*(rhonew(i,j,k)+rhonew(i,j
     $              ,k+1)))

            END DO
         END DO
      END DO


c ******************** END W-VELOCITY **************************


      RETURN
      END





