      PROGRAM OCHANL
      PARAMETER(II=101,JJ=101,II1=II-1,JJ1=JJ-1,II2=II-2,JJ2=JJ-2,KK=1)
C
C==========================================================================C
C                                                                          C
C        THIS PROGRAM IS AN ADAPTATION OF THE BLECK AND BOUDRA (1986)      C
C        ISOPYCNIC COORDINATE GENERAL CIRCULATION MODEL. RECTANGULAR       C
C        BASIN VERSION WITH FREE-SLIP OR NO-SLIP BOUNDARY CONDITION        C 
C                                                                          C
C        II  =  NUMBER OF GRID POINTS IN x-DIRECTION                       C
C        JJ  =  NUMBER OF GRID POINTS IN y-DIRECTION                       C
C        KK  =  NUMBER OF LAYERS                                           C
C        II1 = II - 1                                                      C
C        JJ1 = JJ - 1                                                      C
C        DP          --> LAYER THICKNESS IN PRESSURE UNITS                 C
C        PBOT	     --> PRESSURE AT THE BOTTOM  [dyn/(cm**2)]             C
C        UFLUX       --> LAYER VELOCITY TIMES LAYER THICKNESS              C
C                                                                          C
C        ***UNITS ARE CGS                                                  C
C                                                                          C
C==========================================================================C
C
      DIMENSION DPI(II,JJ),PBOT(II)
      REAL DP0(KK),DISVIS(KK),STRESS(KK),XCONT(KK)
      COMMON U(II,JJ,2*KK),V(II,JJ,2*KK),DP(II,JJ,2*KK),PSI(II,JJ,KK)
      COMMON P(II,JJ,KK+1),UFLUX(II,JJ),VFLUX(II,JJ),THTOP(II,JJ)
      COMMON VORT(II,JJ)
      COMMON DDP(II,JJ,KK+2),DDT(II,JJ,KK+2),DDU(II,JJ,KK+2)
      COMMON STRESX(II,JJ),STRESY(II,JJ),CORIO(II,JJ)
     .,POTVOR(II,JJ)
      COMMON DEL2(II,JJ),VISC(II,JJ)
      REAL UOLD(II,JJ,KK),VOLD(II,JJ,KK)
      EQUIVALENCE (UOLD,DDP)
      REAL UTIL1(II,JJ),UTIL2(II,JJ)
      REAL PROXA(II),PROXB(II),PROXC(II),PROXD(II)
C
C
C -------- MISCELLANEOUS CONSTANTS USED IN TIME INTEGRATION ------------
 
      DATA G/980.6/
      DATA BETA/2.E-13/

C ---------- VISCOSITY COEFFICIENT (LAPLACIAN OR BIHARMONIC) -----------
      DATA VISCOS/1.E6/
      DATA EPS/1.E5/
      DATA DELP0/5.E5/,DEFOR/0./
      DATA WGT1,WGT2/.75,.125/

C ---------------- TIME STEP AND GRID SPACING (IN CM) -----------------
      DATA DELT/1200./,SCALE/20.E5/

C --------------------- INITIAL LAYER THICKNESS -----------------------
      DATA DP0/1000.E5/
      DATA DRAG/1.E-7/

C------------------- INITIAL WIND STRESS PARAMETERS -------------------
      DATA STRESSA/1./
      DATA PSTRES/50.E5/
C
C --------------------------- DENSITY JUMP ----------------------------
      DATA DTHETA/0.0007/
C
C ---------------- INITIALIZE PRESSURE AND VELOCITIES -----------------
      DO 10 I=1,II
      DO 10 J=1,JJ
C
C --------------------- BETA-PLANE OR F-PLANE -------------------------
C
        CORIO(I,J)=1.E-4                                 ! F-PLANE
C        CORIO(I,J)=.93E-4+FLOAT(J-JJ1/2)*SCALE*BETA      ! BETA-PLANE
C
C ---------------------------------------------------------------------

      	P(I,J,KK+1)=0.
        UFLUX(I,J)=0.
        VFLUX(I,J)=0.
        PSI(I,J,KK)=0.
C
C --------------------- WIND STRESS PROFILE --------------------------- C
C                                                                       C
          STRESX(I,J)=-STRESSA*COS(FLOAT(J-1)/FLOAT(JJ1)
     .                *6.28318530718)*G/PSTRES
          STRESY(I,J)=0.
C                                                                       C
C	The wind stress is prescribed as a body force and is assumed    C
c	to decrease linearly to the depth corresponding to PSTRES.      C
c	If the upper layer thickness (actually corresponding pressure)  C
c	becomes less than PSTRES, then the forcing is distributed       C
c	to the second layer and so on according to the linear law.      C
C                                                                       C
c --------------------------------------------------------------------- C
C
      DO 10 K=1,KK
      	P(I,J,1)=0.
      	PSI(I,J,KK)=0.
C
C ------ INITIAL LAYER THICKNESS (EXPRESSED AS PRESSURE G*H) ----------
C
      	DP(I,J,K)=DP0(K)
      	DP(I,J,K+KK)=DP0(K)
C
C ------------------------ INITIAL VELOCITIES -------------------------
C
      	U(I,J,K   )=0.
      	U(I,J,K+KK)=0.
      	V(I,J,K   )=0.
 10    	V(I,J,K+KK)=0.   
C
      DELT1=DELT
      X=1./SCALE
      SCALX=SCALE
      SCALY=SCALE
      SCALX2=SCALX**2
      SCALY2=SCALY**2
      WEIGHT=0.
      EPOT=0.
C
C --------------- DETERMINE TYPE OF BOUNDARY CONDITION ----------------
C               (FREE SLIP, IFREE = 1   NO-SLIP, IFREE = 2)
C
      DATA IFREE/1/
C
      if(ifree.eq.1) print *,'Model run with free-slip boundary cond.'
      if(ifree.eq.2) print *,'Model run with no-slip boundary cond.'
C
C ---------------------- 120 DAYS INTEGRATION -------------------------
C
      NSTEP1=0
      NSTEP2=600*86400/INT(DELT)
      PRINT *,NSTEP2
      NSTEP=0
      M=2
      N=1
      IF (NSTEP2.LE.NSTEP1) STOP
      PRINT 113,NSTEP1
113   FORMAT ('MODEL INTEGRATION STARTS AT TIME STEP',I10)
      IF (NSTEP1.EQ.0) GO TO 13
C
C ------------------ RESTART FROM HISTORY FILE ------------------------
C
      NI=1
24    READ (NI) TIME0,TIME,NSTEP
     .    ,(((U(I,J,K),I=1,II),J=1,JJ),K=1   ,KK   )
     .    ,(((U(I,J,K),I=1,II),J=1,JJ),K=1+KK,KK+KK)
     .    ,(((V(I,J,K),I=1,II),J=1,JJ),K=1   ,KK   )
     .    ,(((V(I,J,K),I=1,II),J=1,JJ),K=1+KK,KK+KK)
     .    ,(((DP(I,J,K),I=1,II),J=1,JJ),K=1   ,KK   )
     .    ,(((DP(I,J,K),I=1,II),J=1,JJ),K=1+KK,KK+KK)
      PRINT 111,NSTEP,NSTEP1
111   FORMAT (25X,'TIME STEP IN HISTORY FILE -',I7,9X,'WANTED -',I7)
      IF (NSTEP.NE.NSTEP1) GO TO 24
C
C
C=======================================================================C
C                                                                       C
C              !!!       MAIN LOOP STARTS HERE         !!!              C
C                                                                       C
C=======================================================================C
C
13    MM=(M-1)*KK
      NN=(N-1)*KK
      NSTEP=NSTEP+1
      TIME=NSTEP*DELT/86400.
C
C ------- ENERGETICS CALCULATIONS WILL BE PERFORMED IF INDEX > 0 --------
C
      INDEX=0
      IF (NSTEP.EQ.NSTEP1+1.OR.NSTEP.EQ.NSTEP2) INDEX=1
C
C ------ DIAGNOSTICS AND HISTORY FILE EVERY 72 TIME STEPS (1 DAY) -------
C
      IF (MOD(NSTEP,144).EQ.0) INDEX=1
C
C=======================================================================C
C                           CONTINUITY EQUATION                         C
C=======================================================================C
C
      DO 889 K=1,KK
      KM=K+MM
      KN=K+NN
C 
      DO 32 J=1,JJ1
      DO 32 I=2,II1
32    UFLUX(I,J)=U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))
     .*SCALY
C         
      DO 16 J=2,JJ1
      DO 16 I=1,II1
16    VFLUX(I,J)=V(I,J,KM)*.5*(DP(I,J,KM)+DP(I,J-1,KM))
     .*SCALX
C
      DO 17 J=1,JJ1
      DO 17 I=1,II1
      DDP(I,J,K)=DP(I,J,KN)
17    DP(I,J,KN)=DP(I,J,KN)-(UFLUX(I+1,J)-UFLUX(I,J)
     .                      +VFLUX(I,J+1)-VFLUX(I,J))*DELT1
     ./(SCALX*SCALY)
889   CONTINUE
C
C
C=======================================================================C
C               END OF  CONTINUITY EQUATION CALCULATIONS                C
C=======================================================================C
C
      DO 15 J=1,JJ1
      DO 15 K=1,KK
      DO 15 I=1,II1
 15   P(I,J,K+1)=P(I,J,K)+DP(I,J,K+MM)
C
C  ------------- TIME SMOOTHING PRESSURE THICKNESS FIELD ----------------
C
      DO 666 K=1,KK
      KN=K+NN
      KM=K+MM
      DO 666 J=1,JJ1
      DO 666 I=1,II1
 666  DP(I,J,KM)=DP(I,J,KM)*.98+(DDP(I,J,K)+DP(I,J,KN))*.01
C
C
C=======================================================================C
C                         HYDROSTATIC EQUATION                          C  
C=======================================================================C
C
      DO 6 J=1,JJ1
      DO 6 I=1,II1
      PSI(I,J,KK)=DP(I,J,KK+MM)*DTHETA
 6    CONTINUE
C
      DO 28 K=1,KK
C      
      KM=K+MM
      KN=K+NN
      XCONT(K)=0.
      STRESS(K)=0.
      DISVIS(K)=0.
C
C ------------- SECOND-ORDER FLUXES FOR MOMENTUM EQUATIONS --------------
C
      DO 93 J=1,JJ1
      DO 93 I=2,II1
 93    UFLUX(I,J)=U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))
C
      DO 94 J=2,JJ1         
      DO 94 I=1,II1
 94    VFLUX(I,J)=V(I,J,KM)*.5*(DP(I,J,KM)+DP(I,J-1,KM))
C
C
C=======================================================================C
C             POTENTIAL  VORTICITY COMPUTATION                          C
C=======================================================================C
C
C ----------------------------- Free-slip ------------------------------
C
      DO 180 I=2,II1
      DO 180 J=2,JJ1
180   POTVOR(I,J)=(CORIO(I,J)+(V(I,J,KM)-V(I-1,J,KM))/SCALX-(U(I,J,KM)-
     .U(I,J-1,KM))/SCALY)/(.25*(DP(I,J,KM)+DP(I-1,J,KM)+DP(I,J-1,KM)+
     .DP(I-1,J-1,KM)))
C
C --- Impose the boundary conditions on the calculation of POTVOR ------
C
      DO 810 J=2,JJ1
      POTVOR(1,J)=CORIO(1,J)/(.5*(DP(1,J,KM)+DP(1,J-1,KM)))
 810  POTVOR(II,J)=CORIO(II,J)/(.5*(DP(II1,J,KM)+DP(II1,J-1,KM)))
C
      DO 820 I=2,II1
      POTVOR(I,1)=CORIO(I,1)
     ./(.5*(DP(I,1,KM)+DP(I-1,1,KM)))
 820  POTVOR(I,JJ)=CORIO(I,JJ)/(.5*(DP(I,JJ1,KM)
     .+DP(I-1,JJ1,KM)))
C
C ------------------------------ No-slip -------------------------------
C
      IF (IFREE.EQ.2) THEN 
C
      DO 81 J=2,JJ1
      POTVOR(1,J)=(CORIO(1,J)+2.*V(1,J,KM)/SCALX)
     ./AMAX1(.5*(DP(1,J,KM)+DP(1,J-1,KM)),EPS)
 81   POTVOR(II,J)=(CORIO(II,J)-2.*V(II1,J,KM)/SCALX)
     ./AMAX1(.5*(DP(II1,J,KM)+DP(II1,J-1,KM)),EPS)
C
      DO 82 I=2,II1
      POTVOR(I,1)=(CORIO(I,1)-2.*U(I,1,KM)/SCALY)
     ./AMAX1(.5*(DP(I,1,KM)+DP(I-1,1,KM)),EPS)
 82   POTVOR(I,JJ)=(CORIO(I,JJ)+2.*U(I,JJ1,KM)/SCALY)
     ./AMAX1(.5*(DP(I,JJ1,KM)+DP(I-1,JJ1,KM)),EPS)
C
      ENDIF
C
C
C======================================================================C
C                              U EQUATION                              C
C======================================================================C
C
C
      DO 37 J=1,JJ1
      DEL2( 1,J)=0.
      DEL2(II,J)=0.
      VISC( 1,J)=DP(  1,J,KM)
      VISC(II,J)=DP(II1,J,KM)
      DO 37 I=2,II1     
      DEL2(I,J)=U(I,J,KN)
      VISC(I,J)=.5*(DP(I,J,KM)+DP(I-1,J,KM))
      UOLD(I,J,K)=U(I,J,KN)
 37   CONTINUE
C
      DO 61 J=1,JJ1
      JA=MAX0(1,J-1)
      JB=MIN0(JJ1,J+1)
C 
C	Set up the free-slip or no-slip boundary conditions
C	(if.ifree.eq.1 ---> free slip boundary condition)
C
         SIG1 = 1.
         SIG2 = 1.
      IF (IFREE.EQ.2) THEN          ! no-slip boundary condition
         IF(J.EQ.  1) SIG1 = -1.
         IF(J.EQ.JJ1) SIG2 = -1.
      ENDIF
C     
      DO 1 I=2,II1
C
C ---------------- HORIZONTAL PRESSURE FORCE ------------------------
C
 1    U(I,J,KN)=U(I,J,KN)-DELT1*((PSI(I,J,K)-PSI(I-1,J,K)
C
C ---------------- KINETIC ENERGY GRADIENT --------------------------
C
     .+.25*(U(I+1,J,KM)**2+V(I,J,KM)**2+V(I,JB,KM)**2
     .     -U(I-1,J,KM)**2-V(I-1,J,KM)**2-V(I-1,JB,KM)**2))/SCALX
C
C ------------------- VORTICITY FLUX ---------------------------------
C
     .-.125*(VFLUX(I,J)+VFLUX(I,JB)+VFLUX(I-1,J)+VFLUX(I-1,JB))*
     .(POTVOR(I,J)+POTVOR(I,JB))
C
     .-.5*(((DEL2(I+1,J)-DEL2( I, J))*(VISC(I,J)+VISC(I+1,J))
     .     -(DEL2( I, J)-DEL2(I-1,J))*(VISC(I,J)+VISC(I-1,J)))/SCALX2
     .    +((SIG2*DEL2( I,JB)-DEL2( I, J))*(VISC(I,J)+VISC( I, JB))
     .     -(DEL2( I, J)-SIG1*DEL2( I,JA))*(VISC(I,J)+VISC( I, JA)))
     ./SCALY2)
     .*VISCOS/AMAX1(.5*(DP(I,J,KM)+DP(I-1,J,KM)),EPS)
C
C -------------------- WIND STRESS FORCING ---------------------------
C
     .-STRESX(I,J)*(AMIN1(P(I,J,K+1)+P(I-1,J,K+1)+1.E3,PSTRES+PSTRES)
     .             -AMIN1(P(I,J,K  )+P(I-1,J,K  )   ,PSTRES+PSTRES))/
     .             (      P(I,J,K+1)+P(I-1,J,K+1)+1.E3
     .                   -P(I,J,K  )-P(I-1,J,K  )))
C
      IF (INDEX.EQ.0) GO TO 61
C
C ----------- EVALUATE TERMS IN KINETIC ENERGY EQUATION -------------
C
      DO 59 I=2,II1
      XCONT(K)=XCONT(K)-U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))
     . *(PSI(I,J,K)-PSI(I-1,J,K))/(G*SCALX)
C
      STRESS(K)=STRESS(K)+U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))
     .*STRESX(I,J)*(AMIN1(P(I,J,K+1)+P(I-1,J,K+1)+1.E3,PSTRES+PSTRES)
     .    -AMIN1(P(I,J,K  )+P(I-1,J,K  )   ,PSTRES+PSTRES))/
     .    (      P(I,J,K+1)+P(I-1,J,K+1)+1.E3
     .          -P(I,J,K  )-P(I-1,J,K  ) )/G
     
      IF (K.EQ.KK) STRESS(K)=STRESS(K)-
     .DRAG*.5*(DP(I,J,KM)+DP(I-1,J,KM))*U(I,J,KM)*U(I,J,KN)/G
C
59    DISVIS(K)=DISVIS(K)-U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))/G
     .*.5*(((DEL2(I+1,J)-DEL2( I, J))*(VISC(I,J)+VISC(I+1,J))
     .     -(DEL2( I, J)-DEL2(I-1,J))*(VISC(I,J)+VISC(I-1,J)))/SCALX2
     .    +((SIG2*DEL2( I,JB)-DEL2( I, J))*(VISC(I,J)+VISC( I,JB ))
     .     -(DEL2( I, J)-SIG1*DEL2( I,JA))*(VISC(I,J)+VISC( I,JA )))
     .      /SCALY2)
     .*VISCOS/AMAX1(.5*(DP(I,J,KM)+DP(I-1,J,KM)),EPS)
61       CONTINUE
C
C
C======================================================================C
C                              V EQUATION                              C
C======================================================================C
C
C              
      DO 34 I=1,II1  
      DEL2(I,1)=0.
      DEL2(I,JJ)=0.
      VISC(I, 1)=DP(I,  1,KM)
      VISC(I,JJ)=DP(I,JJ1,KM)
      DO 34 J=2,JJ1
      DEL2(I,J)=V(I,J,KN)
      VISC(I,J)=.5*(DP(I,J,KM)+DP(I,J-1,KM))
      VOLD(I,J,K)=V(I,J,KN)
 34   CONTINUE
C      
      DO 62 I=1,II1
      IA=MAX0(1,I-1)
      IB=MIN0(II1,I+1)
C
C	Set up the free-slip or no-slip boundary conditions
C	(if.ifree.eq.1 ---> free slip boundary condition)
C
             SIG1 = 1.
             SIG2 = 1.
	IF(IFREE.EQ.2) THEN            ! no-slip boundary condition
             IF(I.EQ.  1) SIG1 = -1.
             IF(I.EQ.II1) SIG2 = -1.
	endif
C    
      DO 2 J=2,JJ1
      JA=MAX0(1,J-1)
      JB=MIN0(JJ1,J+1)
C
 2    V(I,J,KN)=V(I,J,KN)-DELT1*((PSI(I,J,K)-PSI(I,JA,K)
C
C ------------------ GRADIENT OF KINETIC ENERGY -----------------------
C
     .+.25*(V(I,JB ,KM)**2+U(I,J  ,KM)**2+U(I+1,J  ,KM)**2
     .     -V(I,JA ,KM)**2-U(I,JA ,KM)**2-U(I+1,JA ,KM)**2))/SCALY
C
C ----------------------- VORTICITY FLUX ------------------------------
C
     .+.125*(UFLUX(I,J)+UFLUX(I+1,J)+UFLUX(I,JA)+UFLUX(I+1,JA))*
     .(POTVOR(I,J)+POTVOR(I+1,J))
C
     .-.5*(((SIG2*DEL2(IB,J)-DEL2(I,J))*(VISC(I,J)+VISC(IB,J))
     .     -(DEL2(I,J)-SIG1*DEL2(IA,J))*(VISC(I,J)+VISC(IA,J)))/SCALX2
     .     +((DEL2(I,JB)-DEL2(I,J))*(VISC(I,J)+VISC(I,JB))
     .     -(DEL2(I,J)-DEL2(I,JA))*(VISC(I,J)+VISC(I,JA)))/SCALY2)
     .*VISCOS/AMAX1(.5*(DP(I,J,KM)+DP(I,JA ,KM)),EPS)
C
C ------------------- WIND STRESS FORCING -----------------------------
C
     .-STRESY(I,J)*(AMIN1(P(I,J,K+1)+P(I,J-1,K+1)+1.E3,PSTRES+PSTRES)
     .             -AMIN1(P(I,J,K  )+P(I,J-1,K  )   ,PSTRES+PSTRES))/
     .             (      P(I,J,K+1)+P(I,J-1,K+1)+1.E3
     .                   -P(I,J,K  )-P(I,J-1,K  )))
C
      IF (INDEX.EQ.0) GO TO 62
C            
C ------------ EVALUATE TERMS IN KINETIC ENERGY EQUATION --------------
C
      DO 60 J=2,JJ1
      JA=MAX0(1,J-1)
      JB=MIN0(JJ1,J+1)
      XCONT(K)=XCONT(K)-V(I,J,KM)*.5*(DP(I,J,KM)+DP(I,JA ,KM))*
     .(PSI(I,J,K)-PSI(I,JA ,K))/(G*SCALY)
C
      STRESS(K)=STRESS(K)+V(I,J,KM)*.5*(DP(I,J,KM)+DP(I,J-1,KM))
     .*STRESY(I,J)*(AMIN1(P(I,J,K+1)+P(I,J-1,K+1)+1.E3,PSTRES+PSTRES)
     .             -AMIN1(P(I,J,K  )+P(I,J-1,K  )   ,PSTRES+PSTRES))/
     .(                   P(I,J,K+1)+P(I,J-1,K+1)+1.E3
     .                   -P(I,J,K  )-P(I,J-1,K  ))/G
C
      IF (K.EQ.KK) STRESS(K)=STRESS(K)-
     .DRAG*.5*(DP(I,J,KM)+DP(I,J-1,KM))*V(I,J,KM)*V(I,J,KN)/G
C
 60   DISVIS(K)=DISVIS(K)-V(I,J,KM)*.5*(DP(I,J,KM)+DP(I,JA ,KM))/G
     .*.5*(((SIG2*DEL2(IB,J)-DEL2( I, J))*(VISC(I,J)+VISC(IB, J))
     .     -(DEL2(I, J)-SIG1*DEL2(IA, J))*(VISC(I,J)+VISC(IA, J)))/SCALX2
     .    +((DEL2(I,JB)-DEL2( I, J))*(VISC(I,J)+VISC( I,JB))
     .     -(DEL2(I, J)-DEL2( I,JA))*(VISC(I,J)+VISC( I,JA)))/SCALY2)
     .*VISCOS/AMAX1(.5*(DP(I,J,KM)+DP(I,JA ,KM)),EPS)
 62   CONTINUE
 28   CONTINUE
C
C======================================================================C
C                        LINEAR  DRAG                                  C
C======================================================================C
C
      KN=KK+NN
      DO 18 J=1,JJ1
      DO 18 I=1,II1
      U(I,J,KN)=U(I,J,KN)*(1.-DRAG*DELT1)
18    V(I,J,KN)=V(I,J,KN)*(1.-DRAG*DELT1)
C
C ------------- SMOOTH PROGNOSTIC VARIABLES IN TIME -------------------C
C
      DO 22 K=1,KK
      KM=K+MM
      KN=K+NN
      DO 324 J=1,JJ1
      DO 324 I=2,II1
      U(I,J,KM)=U(I,J,KM)*WGT1+(UOLD(I,J,K)+U(I,J,KN))*WGT2
324   CONTINUE
C
      DO 222 J=2,JJ1
      DO 222 I=1,II1
      V(I,J,KM)=V(I,J,KM)*WGT1+(VOLD(I,J,K)+V(I,J,KN))*WGT2
222	CONTINUE
C
22	CONTINUE
C
C======================================================================C
C                OUTPUT AND DIAGNOSTIC CALCULATIONS                    C
C======================================================================C
C
      IF (INDEX.EQ.0) GO TO 25
      PRINT 100,NSTEP,TIME,SCALX,SCALY
100   FORMAT (' TIME STEP',I9,25X,'DAY',F8.3,25X,
     .'GRID SIZE (KM)',-5P2F9.3)
C
C
C
C ------------ ENERGY DIAGNOSTICS ----------------
C
      DPBAR=DP0(KK)
      SUMWGT=0.
      SUMPOT=0.
      SUMKIN=0.
      SUMXGR=0.
      SUMSTR=0.
      SUMDIS=0.
      SUMEKT=0.
      DO 53 K=1,KK
         KM=K+MM
         EPOT=0.
         EKIN=0.
         WEIGHT=0.
      DO 51 J=1,JJ1
      DO 51 I=1,II1
         EPOT=EPOT+.5*(DP(I,J,KM)-DPBAR)**2*DTHETA/G
         EKIN=EKIN+.25*DP(I,J,KM)*(U(I,J,KM)**2+U(I+1,J,KM)**2
     .                            +V(I,J,KM)**2+V(I,J+1 ,KM)**2)/G
51    WEIGHT=WEIGHT+DP(I,J,KM)
      
      EPOT=EPOT/(II1*JJ1)
      EKIN=EKIN/(II1*JJ1)
      DISVIS(K)=-DISVIS(K)/(II1*JJ1)
      XCONT(K)=XCONT(K)/(II1*JJ1)
      STRESS(K)=STRESS(K)/(II1*JJ1)
      WEIGHT=WEIGHT/(II1*JJ1)
      EKINT=STRESS(K)+XCONT(K)+DISVIS(K)
      PRINT 101,WEIGHT,EPOT,EKIN,EKINT,XCONT(K),STRESS(K),DISVIS(K),K
101   FORMAT (' MASS',1PE10.2,'  EPOT',E10.2,'  EKIN',E10.2,'  EKINT'
     .,E10.2,'  XCONT',E10.2,'  STRESS',E10.2,'  DISSIP',E10.2,I4)
      SUMWGT=SUMWGT+WEIGHT
      SUMKIN=SUMKIN+EKIN
      SUMPOT=SUMPOT+EPOT
      SUMXGR=SUMXGR+XCONT(K)
      SUMSTR=SUMSTR+STRESS(K)
      SUMDIS=SUMDIS+DISVIS(K)
      SUMEKT=SUMEKT+EKINT
53    DPBAR=DPBAR+DP0(K)
      PRINT 101,SUMWGT,SUMPOT,SUMKIN,SUMEKT,SUMXGR,SUMSTR,SUMDIS
      
      TIMOUT=TIME
C
25    CONTINUE
C
      IF(INDEX.EQ.1)
     .WRITE(1)M,N,(((U(I,J,K),I=1,II),J=1,JJ),K=1,KK),
     . (((V(I,J,K),I=1,II),J=1,JJ),K=1,KK),
     . (((DP(I,J,K),I=1,II),J=1,JJ),K=1,KK),NSTEP,DELT1      

      IF(NSTEP.GE.NSTEP2) STOP
      L=M
      M=N
      N=L
      DELT1=2.*DELT

      GO TO 13
      PRINT*, "NSTEP = ",NSTEP
      END
