         PROGRAM SHALLOW



      PARAMETER (N1=513, N2=513)
       REAL,DIMENSION(N1,N2):: U,V,P,UNEW,VNEW,
     1        PNEW,UOLD,VOLD,POLD,
     2        CU,CV,Z,H,PSI


      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL *8 T100,T200,T300
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF
      REAL *8 MFS100,MFS200,MFS300,XFLOPS,TSTART,TEND,TOTTIM,T00
 
          REAL IO_TEMP(N1)
    
C       HPF COMPILER DIRECTIVES
       
CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH U::CU,CV,H,P,PNEW,POLD,PSI,UNEW,UOLD
CHPF$ ALIGN WITH U::V,VNEW,VOLD,Z
CHPF$ DISTRIBUTE U(*,BLOCK) ONTO PROCESS
CHPF$ DISTRIBUTE IO_TEMP(BLOCK) ONTO PROCESS
 
      INTERFACE

      SUBROUTINE INITAL(PSI,P,U,V,UOLD,VOLD,POLD)

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::PSI,P,U,V,UOLD,VOLD,POLD

      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

      
CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::P,POLD,PSI,UOLD
CHPF$ ALIGN WITH *U::V,VOLD
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS

      END SUBROUTINE INITAL

      SUBROUTINE CALC1(CU,CV,Z,H,P,U,V)

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::CU,CV,H,Z,P,U,V

      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF


CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::CU,CV,H,P
CHPF$ ALIGN WITH *U::V,Z
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS

      END SUBROUTINE CALC1
      
      SUBROUTINE CALC2(PNEW,UNEW,VNEW,POLD,UOLD,VOLD,CU,CV,Z,H,U)

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::PNEW,UNEW,VNEW,POLD,VOLD,UOLD,CU,CV,H,Z,U

      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF


CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::CU,CV,H,PNEW,POLD,UNEW,UOLD
CHPF$ ALIGN WITH *U::VNEW,VOLD,Z
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO  *PROCESS

      END SUBROUTINE CALC2 
       
      SUBROUTINE CALC3Z(P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW)

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW
      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF


CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::PNEW,POLD,UNEW,UOLD
CHPF$ ALIGN WITH *U::V,VNEW,VOLD,P
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS

      END SUBROUTINE CALC3Z
      
      SUBROUTINE CALC3(P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW)

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW
      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::P,PNEW,POLD,UNEW,UOLD
CHPF$ ALIGN WITH *U::V,VNEW,VOLD
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS

      END SUBROUTINE CALC3

      END INTERFACE
  
      
      CALL INITAL(PSI,P,U,V,UOLD,VOLD,POLD)

C     PRINT INITIAL VALUES

      WRITE( 6,390 ) N,M,DX,DY,DT,ALPHA,ITMAX
  390 FORMAT('1NUMBER OF POINTS IN THE X DIRECTION', I8/
     1       ' NUMBER OF POINTS IN THE Y DIRECTION', I8/
     2       ' GRID SPACING IN THE X DIRECTION    ', F8.0/
     3       ' GRID SPACING IN THE Y DIRECTION    ', F8.0/
     4       ' TIME STEP                          ', F8.0/
     5       ' TIME FILTER PARAMETER              ', F8.3/
     6       ' NUMBER OF ITERATIONS               ', I8)
      MNMIN = MIN0(M,N)
 
C MOVE SHRUNk ARRAYS TO TEMPORARY FOR IO

      FORALL (I = 1:MNMIN) IO_TEMP(I) = POLD(I,I)
      WRITE( 6,391 ) (IO_TEMP(I),I = 1,MNMIN)
  391 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF P ', //(8E15.7))
      FORALL (I = 1:MNMIN) IO_TEMP(I) = UOLD(I,I)
      WRITE( 6,392 ) (IO_TEMP(I),I = 1,MNMIN)
  392 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF U ', //(8E15.7))
      FORALL (I = 1:MNMIN) IO_TEMP(I) = VOLD(I,I)
      WRITE( 6,393 ) (IO_TEMP(I),I = 1,MNMIN)
  393 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF V ', //(8E15.7))

C     DETERMINE OVERHEAD OF TIMING CALLS

      TSTART = CPUTIM()
      T00 = CPUTIM() - TSTART

C     GET TIME STAMP

      CALL TIME0

C     INITIAL CALL TO CPU CLOCK

      TSTART = CPUTIM()
      TIME = 0.
      NCYCLE = 0
   90 NCYCLE = NCYCLE + 1

C     COMPUTE CAPITAL  U, CAPITAL V, Z AND H

      CALL CALC1(CU,CV,Z,H,P,U,V)

C     COMPUTE NEW VALUES U,V AND P

      CALL CALC2(PNEW,UNEW,VNEW,POLD,UOLD,VOLD,CU,CV,Z,H,U)

      TIME = TIME + DT
      IF (MOD(NCYCLE,MPRINT) .NE. 0) GOTO 370
      PTIME = TIME / 3600.
      
      WRITE( 6,350 ) NCYCLE,PTIME
  350 FORMAT(//' CYCLE NUMBER',I5,' MODEL TIME IN  HOURS', F6.2)
      FORALL (I = 1:MNMIN) IO_TEMP(I) = PNEW(I,I)
      WRITE( 6,355 ) (IO_TEMP(I),I = 1,MNMIN)
  355 FORMAT(/' DIAGONAL ELEMENTS OF P ', //(8E15.7))
      FORALL (I = 1:MNMIN) IO_TEMP(I) = UNEW(I,I)
      WRITE( 6,360 ) (IO_TEMP(I),I = 1,MNMIN)
  360 FORMAT(/' DIAGONAL ELEMENTS OF U ', //(8E15.7))
      FORALL (I = 1:MNMIN) IO_TEMP(I) = VNEW(I,I)
      WRITE( 6,360 ) (IO_TEMP(I),I = 1,MNMIN)
  365 FORMAT(/' DIAGONAL ELEMENTS OF V ', //(8E15.7))

C     COMPUTE MFLOPS RATES FOR THE THREE MAJOR LOOPS

      MFS100 = 0
      MFS200 = 0
      MFS300 = 0
      T100 = T100 - T00
      T200 = T200 - T00
      T300 = T300 - T00
      IF (T100 .GT. 0) MFS100 = 24. * M* N * 1.E-6 / T100
      IF (T200 .GT. 0) MFS200 = 26. * M* N * 1.E-6 / T200
      IF (T300 .GT. 0) MFS300 = 15. * M* N * 1.E-6 / T300

      TEND = CPUTIM()
      CTIME = TEND - TSTART - T00
      TCYC = CTIME / FLOAT(NCYCLE)
      WRITE( 6,375 ) NCYCLE,CTIME,TCYC,T100,MFS100,T200,MFS200,T300,
     +    MFS300
  375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', E15.6,
     1       ' TIME PER CYCLE', E15.6 /
     2       ' TIME AND MEGAFLOPS FOR LOOP 100 ', E15.6,F8.3/
     3       ' TIME AND MEGAFLOPS FOR LOOP 200 ', E15.6,F8.3/
     4       ' TIME AND MEGAFLOPS FOR LOOP 300 ', E15.6,F8.3/ )
  370 CONTINUE

C     TEST FOR END OF RUN

       IF (NCYCLE .GE. ITMAX) THEN
       TEND = CPUTIM()
       TOTTIM= (T100+T200+T300)*ITMAX

C       TIME STAMP THE END OF RUN

        CALL TIME0

C        COMPUTE THE MFLOPS FOR THE RUN

        XFLOPS = N * M* ITMAX * 65 * 1.E-6 / TOTTIM
        WRITE( 6,9030 ) TOTTIM,XFLOPS
9030  FORMAT(///'  TOTAL CPU TIME FOR JOB = ', E15.5,/
     1          '  OVERALL MFLOPS RATE    = ', E15.5)
        WRITE( 6,9010 ) 
9010  FORMAT(1H1)

      STOP
      ENDIF

C     TIME SMOOTHING AND UPDATE FOR NEXT CYCLE


      IF (NCYCLE .LE. 1) THEN
      CALL CALC3Z(P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW)

      ELSE

      CALL CALC3(P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW)

      ENDIF

      GOTO 90
      END

      SUBROUTINE TIME0

      RETURN
      END

      SUBROUTINE INITAL(PSI,P,U,V,UOLD,VOLD,POLD)

C     INITIALIZE CONSTANTS AND ARRAYS
           
      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::PSI,P,U,V,UOLD,VOLD,POLD
      
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

C     NOTE BELOW THAT TWO DELTA T (TDT) IS SET TO DT ON THE FIRST
C     CYCLE AFTER WHICH IT IS RESET TO DT+DT.

C     COMPILER DIRECTIVES 

CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::P,POLD,PSI,UOLD
CHPF$ ALIGN WITH *U::V,VOLD
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS
        
      DT = 20.
      TDT = DT
      DX = .25E5
      DY = .25E5
      A = 1.E6
      ALPHA = .001
      ITMAX = 50
      MPRINT = 50

      M= N1 - 1
      N = N2 - 1
      MP1 = M+ 1
      NP1 = N + 1
      EL = FLOAT(N) * DX
      PI = 4. * ATAN(1.)
      TPI = PI + PI
      DI = TPI / FLOAT(M)
      DJ = TPI / FLOAT(N)
      PCF = PI * PI * A * A / (EL * EL)

C     INITIAL VALUES OF THE STREAM FUNCTION AND P

      FORALL (J = 1:NP1, I = 1:MP1) PSI(I,J) = A 
     +    * SIN((FLOAT(I) - .5) * DI) * SIN((FLOAT(J) - .5) * DJ)
      FORALL (J = 1:NP1, I = 1:MP1) P(I,J) = PCF * (
     +    COS(2. * FLOAT(I - 1) * DI) + COS(2. * FLOAT(J - 1) * DJ)) 
     +    + 50000.

C     INITIALIZE VELOCITIES

      U(2:1 + M,:N) = -(PSI(2:1 + M,2:1 + N) - PSI(2:1 + M,:N)) / DY
      V(:M,2:1 + N) = (PSI(2:1 + M,2:1 + N) - PSI(:M,2:1 + N)) / DX

C     PERIODIC CONTINUATION

      U(1,:N) = U(M+ 1,:N)
      V(M+ 1,2:1 + N) = V(1,2:1 + N)
      U(2:1 + M,N + 1) = U(2:1 + M,1)
      V(:M,1) = V(:M,N + 1)
      U(1,N + 1) = U(M+ 1,1)
      V(M+ 1,1) = V(1,N + 1)
      UOLD(:MP1,:NP1) = U(:MP1,:NP1)
      VOLD(:MP1,:NP1) = V(:MP1,:NP1)
      POLD(:MP1,:NP1) = P(:MP1,:NP1)
          
       RETURN
       END

      SUBROUTINE CALC1(CU,CV,Z,H,P,U,V)

C     COMPUTE CAPITAL  U, CAPITAL V, Z AND H

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::CU,CV,H,Z,P,U,V
      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

C      HPF COMPILER DIRECTIVES

CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::CU,CV,H,P
CHPF$ ALIGN WITH *U::V,Z
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS
       
      FSDX = 4. / DX
      FSDY = 4. / DY
      TBEG = CPUTIM()

      CU(2:1 + M,:N) = .5 * (P(2:1 + M,:N) + P(:M,:N)) * U(2:1 + M,:N)
      CV(:M,2:1 + N) = .5 * (P(:M,2:1 + N) + P(:M,:N)) * V(:M,2:1 + N)
      Z(2:1 + M,2:1 + N) = (FSDX * (V(2:1 + M,2:1 + N) - V(:M,2:1 + N))
     +    - FSDY * (U(2:1 + M,2:1 + N) - U(2:1 + M,:N))) / (P(:M,:N)      +    + P(2:1 + M,:N) + P(2:1 + M,2:1 + N) + P(:M,2:1 + N))
      H(:M,:N) = P(:M,:N) + .25 * (U(2:1 + M,:N) * U(2:1 + M,:N) + 
     +    U(:M,:N) * U(:M,:N) + V(:M,2:1 + N) * V(:M,2:1 + N) + V(:M,:N)
     +    *V(:M,:N))
      TEND = CPUTIM()
      T100 = TEND - TBEG

C     PERIODIC CONTINUATION

      CU(1,:N) = CU(M+ 1,:N)
      CV(M+ 1,2:1 + N) = CV(1,2:1 + N)
      Z(1,2:1 + N) = Z(M+ 1,2:1 + N)
      H(M+ 1,:N) = H(1,:N)

      CU(2:1 + M,N + 1) = CU(2:1 + M,1)
      CV(:M,1) = CV(:M,N + 1)
      Z(2:1 + M,1) = Z(2:1 + M,N + 1)
      H(:M,N + 1) = H(:M,1)
      CU(1,N + 1) = CU(M+ 1,1)
      CV(M+ 1,1) = CV(1,N + 1)
      Z(1,1) = Z(M+ 1,N + 1)
      H(M+ 1,N + 1) = H(1,1)

C     END OF LOOP 100 CALCULATIONS

      RETURN
      END


      SUBROUTINE CALC2(PNEW,UNEW,VNEW,POLD,UOLD,VOLD,CU,CV,Z,H,U)

C        COMPUTE NEW VALUES OF U,V,P

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::PNEW,UNEW,VNEW,POLD,VOLD,UOLD,CU,CV,H,Z,U
      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

C     HPF COMPILER DIRECTIVES

CHPF$ PROCESSORS PROCESS(8)

CHPF$ ALIGN WITH *U::CU,CV,H,PNEW,POLD,UNEW,UOLD
CHPF$ ALIGN WITH *U::VNEW,VOLD,Z
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO  *PROCESS
             
      TDTS8 = TDT / 8.
      TDTSDX = TDT / DX
      TDTSDY = TDT / DY
      TBEG = CPUTIM()

      UNEW(2:1 + M,:N) = UOLD(2:1 + M,:N) + TDTS8 * (Z(2:1 + M,2:1 + N)
     +    + Z(2:1 + M,:N)) * (CV(2:1 + M,2:1 + N) + CV(:M,2:1 + N) 
     +    + CV(:M,:N) + CV(2:1 + M,:N)) - TDTSDX * (H(2:1 + M,:N) 
     +    - H(:M,:N))
      VNEW(:M,2:1 + N) = VOLD(:M,2:1 + N) - TDTS8 * (Z(2:1 + M,2:1 + N)
     +    + Z(:M,2:1 + N)) * (CU(2:1 + M,2:1 + N) + CU(:M,2:1 + N) 
     +    + CU(:M,:N) + CU(2:1 + M,:N)) - TDTSDY * (H(:M,2:1 + N) 
     +    - H(:M,:N))
      PNEW(:M,:N) = POLD(:M,:N) - TDTSDX * (CU(2:1 + M,:N) - CU(:M,:N))
     +    - TDTSDY * (CV(:M,2:1 + N) - CV(:M,:N))
      TEND = CPUTIM()
      T200 = TEND - TBEG

C     PERIODIC CONTINUATION

      UNEW(1,:N) = UNEW(M+ 1,:N)
      VNEW(M+ 1,2:1 + N) = VNEW(1,2:1 + N)
      PNEW(M+ 1,:N) = PNEW(1,:N)

      UNEW(2:1 + M,N + 1) = UNEW(2:1 + M,1)
      VNEW(:M,1) = VNEW(:M,N + 1)
      PNEW(:M,N + 1) = PNEW(:M,1)
      UNEW(1,N + 1) = UNEW(M+ 1,1)
      VNEW(M+ 1,1) = VNEW(1,N + 1)
      PNEW(M+ 1,N + 1) = PNEW(1,1)

      RETURN
      END


      SUBROUTINE CALC3(P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW)

C         TIME SMOOTHER

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW
      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

C     HPF COMPILER DIRECTIVES

CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::P,PNEW,POLD,UNEW,UOLD
CHPF$ ALIGN WITH *U::V,VNEW,VOLD
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS
 
C     TIME SMOOTHING AND UPDATE FOR NEXT CYCLE

      TBEG = CPUTIM()

      UOLD(:M,:N) = U(:M,:N) + ALPHA * (UNEW(:M,:N) - 2. * U(:M,:N) 
     +    + UOLD(:M,:N))
      U(:M,:N) = UNEW(:M,:N)
      VOLD(:M,:N) = V(:M,:N) + ALPHA * (VNEW(:M,:N) - 2. * V(:M,:N) 
     +    + VOLD(:M,:N))
      V(:M,:N) = VNEW(:M,:N)
      POLD(:M,:N) = P(:M,:N) + ALPHA * (PNEW(:M,:N) - 2. * P(:M,:N) 
     +    + POLD(:M,:N))
      P(:M,:N) = PNEW(:M,:N)
      TEND = CPUTIM()
      T300 = TEND - TBEG

C     PERIODIC CONTINUATION

      UOLD(M+ 1,:N) = UOLD(1,:N)
      VOLD(M+ 1,:N) = VOLD(1,:N)
      POLD(M+ 1,:N) = POLD(1,:N)
      U(M+ 1,:N) = U(1,:N)
      V(M+ 1,:N) = V(1,:N)
      P(M+ 1,:N) = P(1,:N)

      UOLD(:M,N + 1) = UOLD(:M,1)
      VOLD(:M,N + 1) = VOLD(:M,1)
      POLD(:M,N + 1) = POLD(:M,1)
      U(:M,N + 1) = U(:M,1)
      V(:M,N + 1) = V(:M,1)
      P(:M,N + 1) = P(:M,1)
      UOLD(M+ 1,N + 1) = UOLD(1,1)
      VOLD(M+ 1,N + 1) = VOLD(1,1)
      POLD(M+ 1,N + 1) = POLD(1,1)
      U(M+ 1,N + 1) = U(1,1)
      V(M+ 1,N + 1) = V(1,1)
      P(M+ 1,N + 1) = P(1,1)

      RETURN
      END

      SUBROUTINE CALC3Z(P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW)

C     TIME SMOOTHER FOR FIRST ITERATION

      PARAMETER (N1=513, N2=513)
      REAL,DIMENSION(N1,N2)::P,U,V,POLD,UOLD,VOLD,PNEW,UNEW,VNEW
      COMMON /TIME00/ T100, T200, T300
      REAL * 8 CPUTIM
      REAL * 8 T100,T200,T300,TBEG,TEND
      COMMON /CONS/ DT,TDT,DX,DY,A,ALPHA,ITMAX,MPRINT,M,N,MP1,
     1              NP1,EL,PI,TPI,DI,DJ,PCF

   
C    HPF COMPILER DIRECTIVES

CHPF$ PROCESSORS PROCESS(8)
CHPF$ ALIGN WITH *U::PNEW,POLD,UNEW,UOLD
CHPF$ ALIGN WITH *U::V,VNEW,VOLD,P
CHPF$ DISTRIBUTE U*(*,BLOCK) ONTO *PROCESS
    
  310 TDT = TDT + TDT
      UOLD(:MP1,:NP1) = U(:MP1,:NP1)
      U(:MP1,:NP1) = UNEW(:MP1,:NP1)
      VOLD(:MP1,:NP1) = V(:MP1,:NP1)
      V(:MP1,:NP1) = VNEW(:MP1,:NP1)
      POLD(:MP1,:NP1) = P(:MP1,:NP1)
      P(:MP1,:NP1) = PNEW(:MP1,:NP1)
      RETURN
      END

      FUNCTION CPUTIM()
      REAL * 8 CPUTIM,SECOND

C     FUNCTION TO RETURN THE ELAPSED CPU TIME SINCE THE LAST CALL

      CPUTIM= SECNDS(0.0)
      RETURN
      END      

















