                             PROGRAM RIGIDSF
C
C       ****************************************************************
C      *  THIS PROGRAM IS AN ADAPTATION OF THE BLECK AND BOUDRA (1986) *
C      *       ISOPYCNIC COORDINATE GENERAL CIRCULATION MODEL.         *
C      *  RECTANGULAR BASIN VERSION  WITH OR WITHOUT NO-SLIP BOUNDARY  *
C      *                      CONDITION                                *
C      *****************************************************************
C
C       II  =  NUMBER OF GRID POINTS IN Y-DIRECTION                 C
C       JJ  =  NUMBER OF GRID POINTS IN X-DIRECTION                 C
C       KK  =  NUMBER OF LAYERS                                     C
C       II1 = II - 1                                                C
C       JJ1 = JJ - 1                                                C
C	DP          --> LAYER THICKNESS
C	PBOT	    --> PRESSURE AT THE BOTTOM  [DYN/(CM**2)]
C	UFLUX       --> LAYER VELOCITY TIMES LAYER THICKNESS
C	UTROP       --> BAROTROPIC VELOCITY
C
C======================================================================C
C======================================================================C
C	RECALL THAT THE # OF GRID POINTS HAS TO BE A PRODUCT OF SMALL
C       PRIMES IN ORDER TO MAKE USE OF THE FFT TO SOLVE THE POISSON
C	EQUATION.
C
      PARAMETER(II=101,JJ=II,KK=2,II2=II-2,JJ2=JJ-2,MAX=II2)
      PARAMETER(NW=350)
C
      REAL U(II,JJ,2*KK) ,V(II,JJ,2*KK) ,DP(II,JJ,2*KK),MONTG(II,JJ,KK)
     .    ,P(II,JJ,KK+1) ,CORIO(II,JJ)  ,TOPO(II,JJ)   ,ABSVOR(II,JJ)
     .    ,UTROP(II,JJ)  ,VTROP(II,JJ)  ,VORT(II,JJ)   ,PS(II,JJ)
     .    ,DEL2U(II,JJ)  ,STRESY(II,JJ) ,STRMF(II,JJ)  ,WSAVE(NW)
     .    ,EIG(II,JJ)    ,STRMAUX(II,JJ),UFLUX(II,JJ)  
     .    ,VFLUX(II,JJ)  ,UOLD(II,JJ,KK) ,DEL2V(II,JJ)
     .    ,VOLD(II,JJ,KK),PBOT(JJ)      ,THETA(KK)     ,DP0(KK)
     .    ,DISVIS(KK)    ,STRESS(KK)    ,XCONT(KK)     ,ZETA(II,JJ)
     .    ,UTIL1(II,JJ)  
     .    ,UTIL2(II,JJ)  
     .    ,PGF(II,JJ)    ,DPAV(II,JJ)
     .    ,GRADX(II,JJ)  ,GRADY(II,JJ)  ,DPU(II,JJ,KK) ,DPV(II,JJ,KK)
     .    ,DEPTHU(II,JJ) ,DEPTHV(II,JJ) ,ULD(II,JJ,KK) ,VLD(II,JJ,KK)
     .    ,PU(II,JJ,KK+1),PV(II,JJ,KK+1),UM(II,JJ,KK), VM(II,JJ,KK)
     .    ,DPM(II,JJ,KK) ,STRESX(II,JJ)
c
c
      EQUIVALENCE   (VORT,STRMF)
C
C ---           FUNCTIONS TO MIMIC CRAY VECTOR MASK FUNCTIONS
        CVMGP(A,B,C)=A*(.5+SIGN(.5,C))+B*(.5-SIGN(.5,C))
C    	 SIGN(A1,A2)= |A1|*SIGN A2
C
C======================================================================C
C                   DEFINE PHYSICAL PARAMETERS                         C
C======================================================================C
C
C --- SUBTRACT FOLLOWING VALUES FROM 1. TO OBTAIN RELATIVE LAYER SPECIFIC VOLUME
C
               DATA THETA/0.,0.002/
C
C
C	RHO=1./(1-THETA) GIVES THE RELATION BETWEEN RHO AND THETA
C
C	EXAMPLE: 3-LAYER CASE WITH 0.02534,0.02629 AND 0.026953
C	CORRESPOND TO RHO=1.026,1.027 AND 1.0277
C
C --- HORIZONTAL GRID SCALE
               DATA SCALE/20.E5/             ! IN CM
C
C --- GRAVITATIONAL ACCELERATION
               DATA G/980.6/
C
C --- VELOCITY TIME SMOOTHING WEIGHTS
               DATA WGT1,WGT2/.75,.125/
C
C --- DEPTH OF WIND STRESS FORCING AND LINEAR BOTTOM DRAG COEFFICIENT
C --- BLB IS BOTTOM BOUNDARY LAYER THICKNESS
               DATA PSTRES/100.E5/,DRAG/1.E-7/,BLB/10.E5/
C
C --- TWO LEAP FROG TIME LEVELS
               DATA M,N/2,1/
C
               DATA VISCOS/1.E6/
C
               DATA EPS/1.E5/,EPS1/5.E5/,EPS2/5.E5/
               DATA EPSIL/1.E-9/
C
C --- 'EPS' IS USED WHEREVER LAYER THICKNESS MUST
C --- BE BOUNDED AWAY FROM ZERO.
C
C --- MODEL IS TO BE INTEGRATED FROM TIME STEP 'NSTEP1' TO 'NSTEP2'
C
      READ(5,*) NDAY1,NDAY2
        TYPE*,'NDAY1= ', NDAY1, ' NDAY2 = ', NDAY2
C
      READ(5,*)DELT               ! TIME STEP OF INTEGRATION
        TYPE*,'OUT1 = ',OUT1,' NSTORE = ',NSTORE
C
      READ(5,*)(DP0(K),K=1,KK)
        TYPE*,'DP0(1) = ',DP0(1),' DP0(2) = ',DP0(2)
C
      READ(5,*)STRESSA,IFREE
        TYPE*,'STRESSA = ', STRESSA
      IF(IFREE.EQ.1) TYPE*,'MODEL RUN WITH FREE-SLIP BOUNDARY COND.'
      IF(IFREE.EQ.2) TYPE*,'MODEL RUN WITH NO-SLIP BOUNDARY COND.'
C 
C	NSTD=IFIX(86400./DELT)
C---
      NSTEP1=NDAY1*86400./DELT + .0001
      NSTEP2=NDAY2*86400./DELT + .0001
C---
	       SCALE2=SCALE**2	
               X=1./SCALE
               II1 = II-1
               JJ1 = JJ-1
               KK1 = KK-1
C
C======================================================================C
C                         I/O ADDRESSES                                C
C======================================================================C
C
C	IF THE RUN STARTS FROM DAY 0, THE PROGRAM WON'T READ ANYTHING
C	FROM THE RESTART_FILE.. HOWEVER, IF IT STARTS FROM ANY OTHER DAY,
C	IT WILL TRY TO OPEN A RESTART FILE.  THIS FILE IN FACT IS THE
C	HISTORY_FILE. RENAMES TO RESTART_20K.101X.  THIS HISTORY FILE
C	WAS CREATED IN THE PREVIOUS RUN.
C	THEN, AT THE END OF THE RUN IT WILL AGAIN
C	CREATE A HISTORY_FILE.  WITH THE FINAL STATE OF THE RUN.
C
c ---   HISTORY FILES
C
        open(20,file='u.1',status='new',form='unformatted')
C
        open(21,file='v.1',status='new',form='unformatted')
C
        open(22,file='dp.1 ',status='new',form='unformatted')
C
C ---   RESTART FILES (to be created)
C
        open(30,file='u.r1',status='new',form='unformatted')
C
        open(31,file='v.r1',status='new',form='unformatted')
C
        open(32,file='dp.r1 ',status='new',form='unformatted')
C
C ---   TIME AVERAGE FILES (when desired)
C
C        open(23,file='um.1',status='new',form='unformatted')
C
C        open(24,file='vm.1',status='new',form='unformatted')
C
C        open(25,file='dpm.1',status='new',form='unformatted')
C
C ---   RESTART FILES (to be read)
C
c        open(10,file='u.r2',status='old',form='unformatted')
C
c        open(11,file='v.r2',status='old',form='unformatted')
C
c        open(12,file='dp.r2 ',status='old',form='unformatted')
C
C
C 
C======================================================================C
C                        INITIALIZATION                                C
C======================================================================C
C
C --- SPECIFY INITIAL AND BOUNDARY CONDITIONS
C
      DO 10 I=1,II
          DO 10 J=1,JJ
C
C --- CORIOLIS PARAMETER
	      BETA=2.E-13
          CORIO(I,J)=.83E-4 - FLOAT(I-II1/2)*SCALE*BETA
C
C======================================================================C
C                   ZONAL WIND STRESS PROFILE                          C
C======================================================================C
         STRESX(I,J)=0. 
         STRESY(I,J)=-STRESSA*COS(FLOAT(I)/FLOAT(II1)
     .                *6.28318530718)*G/PSTRES
C======================================================================C
C	THE WIND STRESS IS PRESCRIBED AS A BODY FORCE AND IS ASSUMED
C	TO DECREASE LINEARLY TO THE DEPTH CORRESPONDING TO PSTRES.
C	IF THE UPPER LAYER THICKNESS (ACTUALLY CORRESPONDING PRESSURE)
C	BECOMES LESS THAN PSTRES, THEN THE FORCING IS DISTRIBUTED
C	TO THE SECOND LAYER AND SO ON ACCORDING TO THE LINEAR LAW.
C
          UFLUX(I,J)=0.
          VFLUX(I,J)=0.
          STRMF(I,J)=0.
             DO 10 K=1,KK
             U(I,J,K   )=0.
             U(I,J,K+KK)=0.
             V(J,J,K   )=0.
10           V(I,J,K+KK)=0.
C
C --- HERE INSERT TEST TOPOGRAPHY (A GAUSSIAN-SHAPED RIDGE) WHICH VARIES
C --- ONLY IN THE Y (J) DIRECTION
C --- MEAN DEPTH OF FLUID (H0) AND HEIGHT AND LENGTH SCALE (HM,RL) OF
C --- RIDGE (H0 AND HM IN PRESSURE UNITS, RL IN CM - SO, H0=4000 M
C --- HERE, AND HM=3200 M, RL=250 KM)
      H0=4000.E5
      HM=3200.E5
      RL=250.E5
      Y0=JJ1*SCALE/2.+1.
      DO 63 J=1,JJ1
      Y=J*SCALE-Y0
C --- DEFINE SHAPE OF RIDGE
      DEPTHS=H0-HM*EXP(-Y**2/RL**2)
C --- P ARRAY GIVES DEPTH AT EACH INTERFACE (FROM SURFACE DOWN)
      DO 63 I=1,II1
      P(I,J,1)=0.
      P(I,J,KK+1)=DEPTHS
      DO 52 K=1,KK1
52    P(I,J,K+1)=AMIN1(P(I,J,KK+1),P(I,J,K)+DP0(K))
      DO 63 K=1,KK
      DP(I,J,K   )=P(I,J,K+1)-P(I,J,K)
63    DP(I,J,K+KK)=P(I,J,K+1)-P(I,J,K)
C
      DO 892 J=1,JJ1
      PBOTTM=AMAX1(PBOTTM,P(1,J,KK+1))
  892 PBOT(J)=P(1,J,KK+1)

C
C --- DEFINE DEPTH AND LAYER THICKNESS AT U AND V POINTS
      DO 48 I=2,II1
      DO 48 J=1,JJ1
48    DEPTHU(I,J)=PBOT(J)
      DO 47 J=2,JJ1
      DO 47 I=1,II1
47    DEPTHV(I,J)=UVDEP(PBOT(J),PBOT(J-1))
      DO 347 K=1,KK
      DO 348 I=2,II1
      DO 348 J=1,JJ1
348   DPU(I,J,K)=AMIN1(.5*(DP(I,J,K)+DP(I-1,J,K)),
     .AMAX1(0.,DEPTHU(I,J)-.5*(P(I,J,K)+P(I-1,J,K))))
      DO 349 J=2,JJ1
      DO 349 I=1,II1
349   DPV(I,J,K)=AMIN1(.5*(DP(I,J,K)+DP(I,J-1,K)),
     .AMAX1(0.,DEPTHV(I,J)-.5*(P(I,J,K)+P(I,J-1,K))))
347   CONTINUE
C
C --- CALL INITIALIZATION SUBROUTINE FOR POISSON SOLVER:
C
          CALL POINIT(II-2,JJ-2,II,PBOT,EIG,WSAVE)
C
      DO 893 J=1,JJ
      DO 893 I=1,II
      MONTG(I,J,1)=0.
      DO 894 K=1,KK1
 894  MONTG(I,J,K+1)=MONTG(I,J,K)-P(I,J,K+1)*(THETA(K+1)-THETA(K))
893   TOPO(I,J)=(MONTG(I,J,KK)-P(I,J,KK+1)*(1.-THETA(KK)))/G
C
C
C 
C======================================================================C
                   DELT1=DELT
C
        WRITE(6,551)II
551	FORMAT(2X'NUMBER OF GRID ELEMENTS IN EACH DIRECTION -'I3)
       WRITE(6,552)SCALE
552	FORMAT(2X'GRID SIZE - 'E6.1)
       WRITE(6,553)DELT
553	FORMAT(2X'TIME STEP - 'F8.1)
       WRITE(6,554)VISCOS
554	FORMAT(2X'VISCOSITY - 'E10.4)
C
C
C
C======================================================================C
C --- RESTART FROM RESTART FILES IF COMPUTATION IS NOT TO START         C
C               FROM REST                                              C
C======================================================================C
C
       IF(NSTEP1.EQ.0) GO TO 75
c
C     ---------------------------------------------------------------
C
      DELT1=DELT+DELT
C
C     ---------------------------------------------------------------
C
      NO=10
      READ (NO) NSTEP0,TIME0,U
      NO=11
      READ (NO) NSTEP0,TIME0,V
      NO=12
      READ (NO) NSTEP0,TIME0,DP
C
      IF (MOD(NSTEP0,2).EQ.1) THEN
      L=M
      M=N
      N=L
      ENDIF
C======================================================================C
C
C --- FOLLOWIND DO LOOP NEEDED ONLY IF HISTORY FILE USED TO RESTART
C
C           DO 2 K = 1,KK
C           DO 2 I = 1,II
C           DO 2 J = 1,JJ
C      V(I,J,K+KK) = V(I,J,K)
C      U(I,J,K+KK) = U(I,J,K)
C    2 DP(I,J,K+KK)= DP(I,J,K)
C      DELT1=DELT
C
C======================================================================C
C
       WRITE (6,1112) NSTEP0,NSTEP1
1112    FORMAT (15X'TIME STEP IN HISTORY FILE - 'I7,9X'WANTED - 'I7)
C
C
      IF(NSTEP0.NE.NSTEP1) THEN
      WRITE (6,103) NSTEP0,NSTEP1
103   FORMAT (/' STEP NO.'I5' SHOULD BE'I5)
      ENDIF
C
C
      NN=(N-1)*KK
      DO 791 K=1,KK
      DO 791 I=1,II1
      DO 791 J=1,JJ1
  791 P(I,J,K+1)=P(I,J,K)+DP(I,J,K+NN)
C
      DO 792 J=1,JJ1
  792 PBOT(J)=P(1,J,KK+1)
C
      DO 357 K=1,KK
      DO 358 I=2,II1
      DO 358 J=1,JJ1
358   DPU(I,J,K)=AMIN1(.5*(DP(I,J,K+NN)+DP(I-1,J,K+NN)),
     .AMAX1(0.,DEPTHU(I,J)-.5*(P(I,J,K)+P(I-1,J,K))))
      DO 359 J=2,JJ1
      DO 359 I=1,II1
359   DPV(I,J,K)=AMIN1(.5*(DP(I,J,K+NN)+DP(I,J-1,K+NN)),
     .AMAX1(0.,DEPTHV(I,J)-.5*(P(I,J,K)+P(I,J-1,K))))
357   CONTINUE
C
C
75       	CONTINUE
      IREC=0
      NSTEP=NSTEP0
      WRITE (6,102) NSTEP
102   FORMAT (/' MODEL INTEGRATION STARTS FROM TIME STEP'I6)
C     
C
C
C=======================================================================C
C                                                                       C
C                     MAIN LOOP STARTS HERE                             C
C                     =====================                             C
C                           =========                                   C
C=======================================================================C
C
 
15    MM=(M-1)*KK
      NN=(N-1)*KK
c
c
      NSTEP=NSTEP+1
      TIME=NSTEP*DELT/86400.
      INDEX=0
      IF (AMOD(TIME+.0001,10.).LT..0002) INDEX=1
C
C=======================================================================C
C           CONTINUITY EQUATION (LAYER THICKNESS PREDICTION)            C
C=======================================================================C
C
C --- ADVECTIVE EQUATION (USING 'FLUX CORRECTED TRANSPORT' ALGORITHM)
C
      DO 39 J=1,JJ1
      DO 39 I=1,II1
      UTROP(I,J)=0.
39    VTROP(I,J)=0.
      DO 1 K=1,KK
      KM=K+MM
      KN=K+NN
C
C --- COMPUTE LOW-ORDER (DIFFUSIVE) FLUXES
C
      DO 11 J=1,JJ1
      DO 11 I=2,II1
11    UFLUX(I,J)=U(I,J,KM)*CVMGP(
     .AMIN1(DP(I-1,J,KN),AMAX1(0.,DEPTHU(I,J)-P(I-1,J,K))),
     .AMIN1(DP(I  ,J,KN),AMAX1(0.,DEPTHU(I,J)-P(I  ,J,K))),
     .U(I,J,KM))
C
      DO 12 J=2,JJ1
      DO 12 I=1,II1
12    VFLUX(I,J)=V(I,J,KM)*CVMGP(
     .AMIN1(DP(I,J-1,KN),AMAX1(0.,DEPTHV(I,J)-P(I,J-1,K))),
     .AMIN1(DP(I,J  ,KN),AMAX1(0.,DEPTHV(I,J)-P(I,J  ,K))),
     .V(I,J,KM))
C
C --- ADVANCE  DP  FIELD USING LOW-ORDER (DIFFUSIVE) FLUX VALUES
C
      DO 19 J=1,JJ1
      DO 19 I=1,II1
      UOLD(I,J,K)=DP(I,J,KN)
19    DP(I,J,KN)=DP(I,J,KN)-(UFLUX(I+1,J)-UFLUX(I,J)
     .                      +VFLUX(I,J+1)-VFLUX(I,J))*X*DELT1
C
C --- COMPUTE 'ANTIDIFFUSIVE' (I.E., HIGH-ORDER MINUS LOW-ORDER) FLUXES
C
C
C --- SECOND-ORDER FLUXES FOR CONTINUITY EQUATION
C
      DO 16 J=1,JJ1 
      UFLUX(  2,J)=U(  2,J,KM)*DPU(  2,J,K)-UFLUX(  2,J)
16    UFLUX(II1,J)=U(II1,J,KM)*DPU(II1,J,K)-UFLUX(II1,J)
      DO 17 I=1,II1
      VFLUX(I,  2)=V(I,  2,KM)*DPV(I,  2,K)-VFLUX(I,  2)
17    VFLUX(I,JJ1)=V(I,JJ1,KM)*DPV(I,JJ1,K)-VFLUX(I,JJ1)
      DO 20 J=1,JJ1
      DO 20 I=3,II2
   20 UFLUX(I,J)=U(I,J,KM)*DPU(I,J,K)-UFLUX(I,J)
      DO 21 J=3,JJ2
      DO 21 I=1,II1
   21 VFLUX(I,J)=V(I,J,KM)*DPV(I,J,K)-VFLUX(I,J)
C
C --- AT EACH GRID POINT, DETERMINE THE RATIO OF THE LARGEST PERMISSIBLE
C --- POS. (NEG.) CHANGE IN  DP  TO THE SUM OF ALL INCOMING (OUTGOING) FLUXES
C
C --- FIRST, THE FOUR CORNER POINTS...
      UTIL1(  1,  1)=-(DP(  1,  1,KN)-AMAX1(
     .DP(  1,  1,KN),DP(  2,  1,KN),DP(  1,  2,KN)))
     ./((-AMIN1(0.,UFLUX(  2,  1))
     .   -AMIN1(0.,VFLUX(  1,  2))+EPSIL)*X*DELT1)
      UTIL2(  1,  1)=-(DP(  1,  1,KN)-AMAX1(0.,AMIN1(
     .DP(  1,  1,KN),DP(  2,  1,KN),DP(  1,  2,KN))))
     ./((-AMAX1(0.,UFLUX(  2,  1))
     .   -AMAX1(0.,VFLUX(  1,  2))-EPSIL)*X*DELT1)
C
      UTIL1(  1,JJ1)=-(DP(  1,JJ1,KN)-AMAX1(
     .DP(  1,JJ1,KN),DP(  2,JJ1,KN),DP(  1,JJ2,KN)))
     ./((-AMIN1(0.,UFLUX(  2,JJ1))
     .   +AMAX1(0.,VFLUX(  1,JJ1))+EPSIL)*X*DELT1)
      UTIL2(  1,JJ1)=-(DP(  1,JJ1,KN)-AMAX1(0.,AMIN1(
     .DP(  1,JJ1,KN),DP(  2,JJ1,KN),DP(  1,JJ2,KN))))
     ./((-AMAX1(0.,UFLUX(  2,JJ1))
     .   +AMIN1(0.,VFLUX(  1,JJ1))-EPSIL)*X*DELT1)
C
      UTIL1(II1,  1)=-(DP(II1,  1,KN)-AMAX1(
     .DP(II1,  1,KN),DP(II2,  1,KN),DP(II1,  2,KN)))
     ./(( AMAX1(0.,UFLUX(II1,  1))
     .   -AMIN1(0.,VFLUX(II1,  2))+EPSIL)*X*DELT1)
      UTIL2(II1,  1)=-(DP(II1,  1,KN)-AMAX1(0.,AMIN1(
     .DP(II1,  1,KN),DP(II2,  1,KN),DP(II1,  2,KN))))
     ./(( AMIN1(0.,UFLUX(II1,  1))
     .   -AMAX1(0.,VFLUX(II1,  2))-EPSIL)*X*DELT1)
C
      UTIL1(II1,JJ1)=-(DP(II1,JJ1,KN)-AMAX1(
     .DP(II1,JJ1,KN),DP(II2,JJ1,KN),DP(II1,JJ2,KN)))
     ./(( AMAX1(0.,UFLUX(II1,JJ1))
     .   +AMAX1(0.,VFLUX(II1,JJ1))+EPSIL)*X*DELT1)
      UTIL2(II1,JJ1)=-(DP(II1,JJ1,KN)-AMAX1(0.,AMIN1(
     .DP(II1,JJ1,KN),DP(II2,JJ1,KN),DP(II1,JJ2,KN))))
     ./(( AMIN1(0.,UFLUX(II1,JJ1))
     .   +AMIN1(0.,VFLUX(II1,JJ1))-EPSIL)*X*DELT1)
C
C --- NOW THE REMAINING LATERAL BOUNDARY POINTS...
      DO 25 I=2,II2
      UTIL1(I,  1)=-(DP(I,  1,KN)-AMAX1(
     .DP(I,  1,KN),DP(I-1,  1,KN),DP(I+1,  1,KN),DP(I,  2,KN)))
     ./((AMAX1(0.,UFLUX(I,  1))-AMIN1(0.,UFLUX(I+1,  1))
     .                         -AMIN1(0.,VFLUX(I  ,  2))+EPSIL)*X*DELT1)
      UTIL2(I,  1)=-(DP(I,  1,KN)-AMAX1(0.,AMIN1(
     .DP(I,  1,KN),DP(I-1,  1,KN),DP(I+1,  1,KN),DP(I,  2,KN))))
     ./((AMIN1(0.,UFLUX(I,  1))-AMAX1(0.,UFLUX(I+1,  1))
     .                         -AMAX1(0.,VFLUX(I  ,  2))-EPSIL)*X*DELT1)
      UTIL1(I,JJ1)=-(DP(I,JJ1,KN)-AMAX1(
     .DP(I,JJ1,KN),DP(I-1,JJ1,KN),DP(I+1,JJ1,KN),DP(I,JJ2,KN)))
     ./((AMAX1(0.,UFLUX(I,JJ1))-AMIN1(0.,UFLUX(I+1,JJ1))
     .  +AMAX1(0.,VFLUX(I,JJ1))                         +EPSIL)*X*DELT1)
 25   UTIL2(I,JJ1)=-(DP(I,JJ1,KN)-AMAX1(0.,AMIN1(
     .DP(I,JJ1,KN),DP(I-1,JJ1,KN),DP(I+1,JJ1,KN),DP(I,JJ2,KN))))
     ./((AMIN1(0.,UFLUX(I,JJ1))-AMAX1(0.,UFLUX(I+1,JJ1))
     .  +AMIN1(0.,VFLUX(I,JJ1))                         -EPSIL)*X*DELT1)
C
      DO 26 J=2,JJ2
      UTIL1(  1,J)=-(DP(  1,J,KN)-AMAX1(
     .DP(  1,J,KN),DP(  2,J,KN),DP(  1,J-1,KN),DP(  1,J+1,KN)))
     ./((                      -AMIN1(0.,UFLUX(  2,J  ))
     .  +AMAX1(0.,VFLUX(  1,J))-AMIN1(0.,VFLUX(  1,J+1))+EPSIL)*X*DELT1)
      UTIL2(  1,J)=-(DP(  1,J,KN)-AMAX1(0.,AMIN1(
     .DP(  1,J,KN),DP(  2,J,KN),DP(  1,J-1,KN),DP(  1,J+1,KN))))
     ./((                      -AMAX1(0.,UFLUX(  2,J  ))
     .  +AMIN1(0.,VFLUX(  1,J))-AMAX1(0.,VFLUX(  1,J+1))-EPSIL)*X*DELT1)
      UTIL1(II1,J)=-(DP(II1,J,KN)-AMAX1(
     .DP(II1,J,KN),DP(II2,J,KN),DP(II1,J-1,KN),DP(II1,J+1,KN)))
     ./((AMAX1(0.,UFLUX(II1,J))
     .  +AMAX1(0.,VFLUX(II1,J))-AMIN1(0.,VFLUX(II1,J+1))+EPSIL)*X*DELT1)
 26   UTIL2(II1,J)=-(DP(II1,J,KN)-AMAX1(0.,AMIN1(
     .DP(II1,J,KN),DP(II2,J,KN),DP(II1,J-1,KN),DP(II1,J+1,KN))))
     ./((AMIN1(0.,UFLUX(II1,J))
     .  +AMIN1(0.,VFLUX(II1,J))-AMAX1(0.,VFLUX(II1,J+1))-EPSIL)*X*DELT1)
C
C --- FINALLY, THE INTERIOR GRID POINTS...
      DO 27 J=2,JJ2
      DO 27 I=2,II2
      UTIL1(I,J)=-(DP(I,J,KN)-AMAX1(
     .DP(I,J,KN),DP(I-1,J,KN),DP(I+1,J,KN),DP(I,J-1,KN),DP(I,J+1,KN)))
     ./((AMAX1(0.,UFLUX(I,J))-AMIN1(0.,UFLUX(I+1,J))
     .  +AMAX1(0.,VFLUX(I,J))-AMIN1(0.,VFLUX(I,J+1))+EPSIL)*X*DELT1)
 27   UTIL2(I,J)=-(DP(I,J,KN)-AMAX1(0.,AMIN1(
     .DP(I,J,KN),DP(I-1,J,KN),DP(I+1,J,KN),DP(I,J-1,KN),DP(I,J+1,KN))))
     ./((AMIN1(0.,UFLUX(I,J))-AMAX1(0.,UFLUX(I+1,J))
     .  +AMIN1(0.,VFLUX(I,J))-AMAX1(0.,VFLUX(I,J+1))-EPSIL)*X*DELT1)
C
C --- LIMIT ANTIDIFFUSIVE FLUXES
C --- (RETAIN INFORMATION ABOUT FLUX CLIPPING IN -UTROP,VTROP-. THIS
C --- WILL ALLOW US LATER TO RESTORE NONDIVERGENCE OF BAROTROPIC FLOW.)
C
      DO 28 J=1,JJ1
      DO 28 I=2,II1
      UTROP(I,J)=UTROP(I,J)+UFLUX(I,J)*
     .(1.-CVMGP(AMIN1(1.,UTIL1(I,J),UTIL2(I-1,J)),
     .          AMIN1(1.,UTIL2(I,J),UTIL1(I-1,J)),UFLUX(I,J)))
28    UFLUX(I,J)=UFLUX(I,J)*CVMGP(AMIN1(1.,UTIL1(I,J),UTIL2(I-1,J)),
     .                            AMIN1(1.,UTIL2(I,J),UTIL1(I-1,J)),
     .UFLUX(I,J))
      DO 29 J=2,JJ1
      DO 29 I=1,II1
      VTROP(I,J)=VTROP(I,J)+VFLUX(I,J)*
     .(1.-CVMGP(AMIN1(1.,UTIL1(I,J),UTIL2(I,J-1)),
     .          AMIN1(1.,UTIL2(I,J),UTIL1(I,J-1)),VFLUX(I,J)))
29    VFLUX(I,J)=VFLUX(I,J)*CVMGP(AMIN1(1.,UTIL1(I,J),UTIL2(I,J-1)),
     .                            AMIN1(1.,UTIL2(I,J),UTIL1(I,J-1)),
     .VFLUX(I,J))
C
C --- ADD ANTIDIFFUSIVE FLUXES TO  DP  FIELD
C
      DO 1 J=1,JJ1
      DO 1 I=1,II1
1     DP(I,J,KN)=DP(I,J,KN)-(UFLUX(I+1,J)-UFLUX(I,J)
     .                      +VFLUX(I,J+1)-VFLUX(I,J))*X*DELT1
C
C --- RESTORE NONDIVERGENCE OF VERTICALLY INTEGRATED FLOW
C
      DO 5 J=1,JJ1
      DO 5 K=1,KK
      DO 5 I=1,II1
5     P(I,J,K+1)=P(I,J,K)+DP(I,J,K+NN)
      DO 14 K=1,KK
      KN=K+NN
      DO 44 J=1,JJ1
      DO 44 I=2,II1
44    UFLUX(I,J)=UTROP(I,J)*
     .   CVMGP(DP(I-1,J,KN)/PBOT(J),DP(I,J,KN)/PBOT(J),UTROP(I,J))

      DO 45 J=2,JJ1
      DO 45 I=1,II1
45    VFLUX(I,J)=VTROP(I,J)*
     .   CVMGP(DP(I,J-1,KN)/PBOT(J-1),DP(I,J,KN)/PBOT(J),VTROP(I,J))

      DO 14 J=1,JJ1
      DO 14 I=1,II1
14    DP(I,J,KN)=DP(I,J,KN)-(UFLUX(I+1,J)-UFLUX(I,J)
     .                      +VFLUX(I,J+1)-VFLUX(I,J))*X*DELT1
c
C
C ---                   END OF  F C T  CALCULATIONS
C
      DO 120 J=1,JJ1
      DO 120 I=1,II1
120   MONTG(I,J,KK)=G*TOPO(I,J)+P(I,J,KK+1)*(1.-THETA(KK))
C
C=======================================================================C
C                       HYDROSTATIC EQUATION                            C
C=======================================================================C
C
      DO 83 J=1,JJ1
      DO 8 K=1,KK
      DO 8 I=1,II1
    8      P(I,J,K+1)=P(I,J,K)+DP(I,J,K+MM)
      DO 13 KI=2,KK
      K=KK+1-KI
      DO 13 I=1,II1
   13    MONTG(I,J,K)=MONTG(I,J,K+1)+P(I,J,K+1)*(THETA(K+1)-THETA(K))
C
C --- TIME SMOOTHING OF THICKNESS FIELD
C
      DO 83 K=1,KK
      DO 83 I=1,II1
   83      DP(I,J,K+MM)=DP(I,J,K+MM)*.98+(UOLD(I,J,K)+DP(I,J,K+NN))*.01
C
C=======================================================================C
C                       MOMENTUM EQUATIONS                              C
C=======================================================================C
C
      DO 9 K=1,KK
      KM=K+MM
      KN=K+NN
      XCONT(K)=0.
      DISVIS(K)=0.
      STRESS(K)=0.
C
C=======================================================================C
C                        VORTICITY COMPUTATION                          C
C                       FOR A RECTANGULAR BASIN:                        C
C=======================================================================C
C
      DO 80 I=2,II1
      DO 80 J=2,JJ1
   80 ABSVOR(I,J)=(CORIO(I,J)+(V(I,J,KM)-V(I-1,J,KM)-U(I,J,KM)
     .+U(I,J-1,KM))*X)
C
C     IMPOSE THE BOUNDARY CONDITIONS ON THE CALCULATION OF ABSVOR
C	NEAR THE BOUNDARIES
C
	      IF(IFREE.EQ.2) THEN
C
      DO 81 J=2,JJ1
      ABSVOR(1,J)=(CORIO(1,J)+2.*V(1,J,KM)*X)
   81 ABSVOR(II,J)=(CORIO(II,J)-2.*V(II1,J,KM)*X)
      DO 82 I=2,II1
      ABSVOR(I,1)=(CORIO(I,1)-2.*U(I,1,KM)*X)
 82   ABSVOR(I,JJ)=(CORIO(I,JJ)+2.*U(I,JJ1,KM)*X)
      ABSVOR(1, 1) = CORIO(1, 1)
      ABSVOR(1,JJ) = CORIO(1,JJ)
      ABSVOR(II,1) = CORIO(II,1)
      ABSVOR(II,JJ)= CORIO(II,JJ)
C
	      ELSE                                  ! FOR FREE SLIP
C
      DO 881 J=2,JJ1
      ABSVOR(1,J)=CORIO(1,J)
  881 ABSVOR(II,J)=CORIO(II,J)
      DO 882 I=2,II1
      ABSVOR(I,1)=CORIO(I,1)
882   ABSVOR(I,JJ)=CORIO(I,JJ)
      ABSVOR(1, 1) = CORIO(1, 1)
      ABSVOR(1,JJ) = CORIO(1,JJ)
      ABSVOR(II,1) = CORIO(II,1)
      ABSVOR(II,JJ)= CORIO(II,JJ)
	      ENDIF
C
C======================================================================C
C                            U EQUATION                                 C
C=======================================================================C
C
C=======================================================================C
C                THICKNESS FIELD FOR LATERAL DIFFUSION                  C
C=======================================================================C
C 
      DO 412 J=1,JJ1
      PS(1,J)=DPU(1,J,K)
      PS(II,J)=DPU(II1,J,K)
      UOLD(1,J,K)=0.
      UOLD(II,J,K)=0.
C
      DO 412 I=2,II1
      PS(I,J)=DPU(I,J,K)
      UOLD(I,J,K)=U(I,J,KN)
  412 CONTINUE
C
      DO 37 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
             IF(J.EQ.  1) SIG1 = -1.
             IF(J.EQ.JJ1) SIG2 = -1.
	      ELSE
             IF(J.EQ.  1) SIG1 = 1.
             IF(J.EQ.JJ1) SIG2 = 1.
	      ENDIF
C
C --- RECTANGULAR BASIN
C
C	...THIS DEL2U TERM IS NECCESARY FOR THE LAPLACIAN VISCOSITY TERM...
C
      DO 37 I=2,II1
   37 DEL2U(I,J)=
     .((PS(I+1,J)+PS(I,J))*(UOLD(I+1,J,K)-UOLD(I,J,K))-
     . (PS(I,J)+PS(I-1,J))*(UOLD(I,J,K)-UOLD(I-1,J,K))+
     . (PS(I,JB)+PS(I,J))*(SIG2*UOLD(I,JB,K)-UOLD(I,J,K))-
     . (PS(I,J )+PS(I,JA))*(UOLD(I,J,K)-SIG1*UOLD(I,JA,K)))*.5
C
      DO 130 J=1,JJ
      DO 130 I=1,II
      PGF(I,J)=0.
  130 GRADX(I,J)=0.
C
C --- FIND PRESSURE GRADIENT FOR U EQUATION
C
      DO 820 I=2,II1
      DO 820 J=1,JJ1
      PGF(I,J)=(MONTG(I,J,K)-MONTG(I-1,J,K))/(1.-THETA(K))
 820  GRADX(I,J)=PGF(I,J)
C --- 'SCHEME 2' FROM APPENDIX A IN BLECK-SMITH PAPER
C --- NOTE THIS IS PRIMARILY NEEDED IN CASES WITH STEEP TOPOGRAPHY
      IF (K.EQ.1) GO TO 822
      DO 861 I=2,II1
      DO 861 J=1,JJ1
861   DPAV(I,J)=DPU(I,J,K)
      DO 851 J=1,JJ1
      JA=MAX0(1  ,J-1)
      JB=MIN0(JJ1,J+1)
      DO 851 I=2,II1
      GRADX(I,J)=CVMGP(PGF(I,J),
     .(AMIN1(DPAV(I,J),EPS1)*PGF(I,J)+(EPS1-AMIN1(DPAV(I,J),EPS1))*
     . (PGF(I-1,J )*AMIN1(DPAV(I-1,J ),EPS1)
     . +PGF(I+1,J )*AMIN1(DPAV(I+1,J ),EPS1)
     . +PGF(I  ,JA)*AMIN1(DPAV(I  ,JA),EPS1)
     . +PGF(I  ,JB)*AMIN1(DPAV(I  ,JB),EPS1))/
     .(AMIN1(DPAV(I-1,J ),EPS1)+AMIN1(DPAV(I+1,J ),EPS1)
     .+AMIN1(DPAV(I  ,JA),EPS1)+AMIN1(DPAV(I  ,JB),EPS1)+.001))/EPS1,
     .EPS1-AMIN1(P(I,J,K),P(I-1,J,K)))
 851  CONTINUE
822   CONTINUE
C
C=======================================================================C
C            BEGIN INTEGRATION OF U MOMENTUM EQUATION                   C
C=======================================================================C
      DO 61 J=1,JJ1
      JA=MAX0(  1,J-1)
      JB=MIN0(JJ1,J+1)
      DO 6 I=2,II1
6     U(I,J,KN)=
C
C --- HORIZONTAL PRESSURE FORCE (X-DIRECTION)
     .U(I,J,KN)-DELT1*((GRADX(I,J)
C
C --- GRADIENT OF KINETIC ENERGY
     .+.25*(U(I+1,J,KM)**2+V(I  ,J,KM)**2+V(I  ,J+1,KM)**2
     .     -U(I-1,J,KM)**2-V(I-1,J,KM)**2-V(I-1,J+1,KM)**2))*X
C
C --- WIND FORCING STRESX (AS FOR V EQUATION) IF DESIRED
C
C --- ABSOLUTE VORTICITY FLUX
     .-.125*(V(I,J,KM)+V(I,J+1,KM)+V(I-1,J,KM)+V(I-1,J+1,KM))*
     .(ABSVOR(I,J)+ABSVOR(I,J+1))
C
C --- VISCOUS DIFFUSION
     .-DEL2U(I,J)*VISCOS*X*X/AMAX1(DPU(I,J,K),EPS1))
C
      IF (INDEX.EQ.0) GO TO 61
C --- EVALUATE TERMS IN KINETIC ENERGY EQUATION
      DO 59 I=2,II1
C
      XCONT(K)=XCONT(K)-U(I,J,KM)*DPU(I,J,K)*GRADX(I,J)/(G*SCALE)
C
      IF (K.EQ.KK) STRESS(K)=STRESS(K)-
     .DRAG*DPU(I,J,K)*U(I,J,KM)*U(I,J,KN)/G
C
59    DISVIS(K)=DISVIS(K)-U(I,J,KM)*DPU(I,J,K)/G
     .*DEL2U(I,J)*VISCOS*X*X
     ./AMAX1(DPU(I,J,K),EPS1)
61    CONTINUE
C
C=======================================================================C
C                             V EQUATION                                C
C=======================================================================C
C
      DO 76 I=1,II1
      PS(I,1)=DPV(I,1,K)
      PS(I,JJ)=DPV(I,JJ1,K)
      VOLD(I,1,K)=0.
      VOLD(I,JJ,K)=0.
      DO 76 J=2,JJ1
      PS(I,J)=DPV(I,J,K)
      VOLD(I,J,K)=V(I,J,KN)
   76 CONTINUE
C
C --- RECTANGULAR BASIN
C
      DO 420 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
             IF(I.EQ.  1) SIG1 = -1.
             IF(I.EQ.II1) SIG2 = -1.
	      ELSE
             IF(I.EQ.  1) SIG1 = 1.
             IF(I.EQ.II1) SIG2 = 1.
	      ENDIF
C
      DO 420 J=2,JJ1
  420 DEL2V(I,J)=
     .((PS(IB,J)+PS(I,J))*(SIG2*VOLD(IB,J,K)-VOLD(I,J,K))-
     . (PS(I,J)+PS(IA,J))*(VOLD(I,J,K)-SIG1*VOLD(IA,J,K))+
     . (PS(I,J+1)+PS(I,J))*(VOLD(I,J+1,K)-VOLD(I,J,K))-
     . (PS(I,J)+PS(I,J-1))*(VOLD(I,J,K)-VOLD(I,J-1,K)))*.5
C
      DO 140 J=1,JJ
      DO 140 I=1,II
      PGF(I,J)=0.
  140 GRADY(I,J)=0.
C
C --- FIND PRESSURE GRADIENT FOR V EQUATION
C
      DO 821 J=2,JJ1
      DO 821 I=1,II1
      PGF(I,J)=(MONTG(I,J,K)-MONTG(I,J-1,K))/(1.-THETA(K))
 821  GRADY(I,J)=PGF(I,J)
C --- 'SCHEME 2' FROM APPENDIX A IN BLECK-SMITH PAPER
      IF (K.EQ.1) GO TO 824
      DO 862 J=2,JJ1
      DO 862 I=1,II1
862   DPAV(I,J)=DPV(I,J,K)
      DO 860 I=1,II1
      IA=MAX0(  1,I-1)
      IB=MIN0(II1,I+1)
      DO 860 J=2,JJ1
      GRADY(I,J)=CVMGP(PGF(I,J),
     .(AMIN1(DPAV(I,J),EPS1)*PGF(I,J)+(EPS1-AMIN1(DPAV(I,J),EPS1))*
     .(PGF(I,J-1)*AMIN1(DPAV(I,J-1),EPS1)
     .+PGF(I,J+1)*AMIN1(DPAV(I,J+1),EPS1)
     .+PGF(IA,J)*AMIN1(DPAV(IA,J),EPS1)
     .+PGF(IB,J)*AMIN1(DPAV(IB,J),EPS1))/
     .(AMIN1(DPAV(I,J-1),EPS1)+AMIN1(DPAV(I,J+1),EPS1)
     .+AMIN1(DPAV(IA,J),EPS1)+AMIN1(DPAV(IB,J),EPS1)+.001))/EPS1,
     .EPS1-AMIN1(P(I,J,K),P(I,J-1,K)))
 860  CONTINUE
 824  CONTINUE
C
C======================================================================C
C                 BEGIN INTEGRATION OF V-MOMENTUM-EQUATION             C
C======================================================================C
C
      DO 62 I=1,II1
      IA=MAX0(  1,I-1)
      IB=MIN0(II1,I+1)
      DO 7 J=2,JJ1
7     V(I,J,KN)=
C
C --- HORIZONTAL PRESSURE FORCE (Y-DIRECTION)
     .V(I,J,KN)-DELT1*((GRADY(I,J)
C
C --- GRADIENT OF KINETIC ENERGY
     .+.25*(V(I,J+1,KM)**2+U(I,J  ,KM)**2+U(I+1,J  ,KM)**2
     .     -V(I,J-1,KM)**2-U(I,J-1,KM)**2-U(I+1,J-1,KM)**2))*X
C
C --- ABSOLUTE VORTICITY FLUX
     .+.125*(U(I,J,KM)+U(I+1,J,KM)+U(I,J-1,KM)+U(I+1,J-1,KM))*
     .(ABSVOR(I,J)+ABSVOR(I+1,J))
C
C --- WIND STRESS FORCING
     .-STRESY(I,J)*(AMIN1(P(I,J,K+1)+P(I,J-1,K+1)+1.,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.
     .                   -P(I,J,K  )-P(I,J-1,K  ))
C
C --- VISCOUS DIFFUSION
     .-DEL2V(I,J)*X*X*VISCOS/AMAX1(DPV(I,J,K),EPS1))
      IF (INDEX.EQ.0) GO TO 62
C
C --- EVALUATE TERMS IN KINETIC ENERGY EQUATION
      DO 60 J=2,JJ1
      XCONT(K)=XCONT(K)-V(I,J,KM)*DPV(I,J,K)*GRADY(I,J)/(G*SCALE)
C
      STRESS(K)=STRESS(K)+V(I,J,KM)*DPV(I,J,K)*
     .STRESY(I,J)*(AMIN1(P(I,J,K+1)+P(I,J-1,K+1)+1.,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.
     .    -P(I,J,K  )-P(I,J-1,K  ))/G
C
      IF (K.EQ.KK) STRESS(K)=STRESS(K)-
     .DRAG*DPV(I,J,K)*V(I,J,KM)*V(I,J,KN)/G
C
60    DISVIS(K)=DISVIS(K)-V(I,J,KM)*DPV(I,J,K)/G
     .*DEL2V(I,J)*X*X*VISCOS
     ./AMAX1(DPV(I,J,K),EPS1)
62    CONTINUE
9     CONTINUE
C
      DO 330 K=1,KK
      KN=K+NN
      DO 330 I=1,II
      DO 330 J=1,JJ
  330 P(I,J,K+1)=P(I,J,K)+DP(I,J,KN)
C
C --- SUBSTITUTE DEPTH-WEIGHTED AVERAGES FOR (U,V) AT GRID POINTS WITH
C --- THICKNESS LESS THAN 5 DB  (I.E., 5 M)
C
      DO 771 K=1,KK+1
      DO 771 J=1,JJ1
      DO 771 I=2,II1
  771 PU(I,J,K)=AMIN1(.5*(P(I,J,K)+P(I-1,J,K)),DEPTHU(I,J))
C
      DO 772 K=1,KK+1
      DO 772 J=2,JJ1
      DO 772 I=1,II1
  772 PV(I,J,K)=AMIN1(.5*(P(I,J,K)+P(I,J-1,K)),DEPTHV(I,J))
C
      DO 774 K=1,KK
      DO 668 J=2,JJ1
      DO 668 I=1,II1
  668 VLD(I,J,K)=V(I,J,K+NN)*(PV(I,J,K+1)-PV(I,J,K))
C
      DO 773 K1=1,K-1
      DO 773 J=2,JJ1
      DO 773 I=1,II1
  773 VLD(I,J,K)=VLD(I,J,K)+V(I,J,K1+NN)*
     .(AMIN1(PV(I,J,K  ),AMAX1(PV(I,J,K+1)-EPS2,PV(I,J,K1+1)))
     .-AMIN1(PV(I,J,K  ),AMAX1(PV(I,J,K+1)-EPS2,PV(I,J,K1  ))))
C
      DO 774 K1=K+1,KK
      DO 774 J=2,JJ1
      DO 774 I=1,II1
  774 VLD(I,J,K)=VLD(I,J,K)+V(I,J,K1+NN)*
     .(AMAX1(PV(I,J,K+1),AMIN1(PV(I,J,K  )+EPS2,PV(I,J,K1+1)))
     .-AMAX1(PV(I,J,K+1),AMIN1(PV(I,J,K  )+EPS2,PV(I,J,K1  ))))
C
      DO 663 K=1,KK
      DO 663 J=2,JJ1
      DO 663 I=1,II1
  663 V(I,J,K+NN)=VLD(I,J,K)/
     .(AMAX1(PV(I,J,K+1),AMIN1(PV(I,J,K  )+EPS2,PV(I,J,KK+1)))
     .-AMIN1(PV(I,J,K  ),AMAX1(PV(I,J,K+1)-EPS2,PV(I,J,   1))))
C
      DO 863 J=1,JJ1
      DO 874 K=1,KK
      DO 888 I=2,II1
  888 ULD(I,J,K)=U(I,J,K+NN)*(PU(I,J,K+1)-PU(I,J,K))
      DO 873 K1=1,K-1
      DO 873 I=2,II1
  873 ULD(I,J,K)=ULD(I,J,K)+U(I,J,K1+NN)*
     .(AMIN1(PU(I,J,K  ),AMAX1(PU(I,J,K+1)-EPS2,PU(I,J,K1+1)))
     .-AMIN1(PU(I,J,K  ),AMAX1(PU(I,J,K+1)-EPS2,PU(I,J,K1  ))))
C
      DO 874 K1=K+1,KK
      DO 874 I=2,II1
  874 ULD(I,J,K)=ULD(I,J,K)+U(I,J,K1+NN)*
     .(AMAX1(PU(I,J,K+1),AMIN1(PU(I,J,K  )+EPS2,PU(I,J,K1+1)))
     .-AMAX1(PU(I,J,K+1),AMIN1(PU(I,J,K  )+EPS2,PU(I,J,K1  ))))
C
      DO 863 K=1,KK
      DO 863 I=2,II1
  863 U(I,J,K+NN)=ULD(I,J,K)/
     .(AMAX1(PU(I,J,K+1),AMIN1(PU(I,J,K  )+EPS2,PU(I,J,KK+1)))
     .-AMIN1(PU(I,J,K  ),AMAX1(PU(I,J,K+1)-EPS2,PU(I,J,1   ))))
C
      DO 367 K=1,KK
      KN=K+NN
      DO 368 J=1,JJ1
      DO 368 I=2,II1
368   DPU(I,J,K)=AMIN1(.5*(DP(I,J,KN)+DP(I-1,J,KN)),
     .AMAX1(0.,DEPTHU(I,J)-.5*(P(I,J,K)+P(I-1,J,K))))
      DO 369 J=2,JJ1
      DO 369 I=1,II1
369   DPV(I,J,K)=AMIN1(.5*(DP(I,J,KN)+DP(I,J-1,KN)),
     .AMAX1(0.,DEPTHV(I,J)-.5*(P(I,J,K)+P(I,J-1,K))))
367   CONTINUE
C
C======================================================================C
C                         BOTTOM DRAG                                  C
C======================================================================C
C
      DO 18 K=1,KK
        KN=K+NN
      DO 188 J=1,JJ1
      DO 188 I=2,II1
 188  U(I,J,KN)=CVMGP(U(I,J,KN)*(1.-DRAG*DELT1),U(I,J,KN),
     . BLB+0.5*(P(I,J,K+1)+P(I-1,J,K+1)-P(I,J,KK+1)-P(I-1,J,KK+1)))
      DO 18 J=2,JJ1
      DO 18 I=1,II1
 18   V(I,J,KN)=CVMGP(V(I,J,KN)*(1.-DRAG*DELT1),V(I,J,KN),
     . BLB+0.5*(P(I,J,K+1)+P(I,J-1,K+1)-P(I,J,KK+1)-P(I,J-1,KK+1)))
C
C======================================================================C
C REMOVE DIVERGENT COMPONENT FROM MEAN MOTION (RIGID LID APPROXIMAT.)  C
C======================================================================C
C
      DO 36 J=1,JJ1
      DO 36 I=1,II1
        P(I,J,KK+1)=0.
         UTROP(I,J)=0.
36        VTROP(I,J)=0.
      DO 32 K=1,KK
       KN=K+NN
        DO 30 J=1,JJ1
         DO 30 I=1,II1
30        P(I,J,KK+1)=P(I,J,KK+1)+DP(I,J,KN)
      DO 31 J=1,JJ1
      DO 31 I=2,II1
   31 UTROP(I,J)=UTROP(I,J)+U(I,J,KN)*DPU(I,J,K)
      DO 32 J=2,JJ1
      DO 32 I=1,II1
   32 VTROP(I,J)=VTROP(I,J)+V(I,J,KN)*DPV(I,J,K)
C
      DO 34 J=1,JJ1
      DO 34 I=2,II1
34    UTROP(I,J)=UTROP(I,J)/DEPTHU(I,J)
      DO 35 J=2,JJ1
      DO 35 I=1,II1
35    VTROP(I,J)=VTROP(I,J)/DEPTHV(I,J)
      DO 33 J=2,JJ1
      DO 33 I=2,II1
33    VORT(I,J)=(VTROP(I,J)-VTROP(I-1,J)-UTROP(I,J)+UTROP(I,J-1))
C
C --- SOLVE POISSON EQUATION TO DETERMINE ROTATIONAL PART OF MEAN MOTION
C --- FOR TESTING THE FFT ROUTINE, ACTIVATE THE TWO STATEMENTS ..PRINT 110...
C --- AND THE LOOP DO 562.... THE RESULTS SHOULD AGREE
C       IF (INDEX.GT.0) WRITE(6,110) (VORT(I,I),I=2,13)
C 110   FORMAT (1X,1P12E10.3)
C
       CALL POISSON(VORT(2,2),II,II2,JJ2,PBOT,EIG,WSAVE)
C
C      DO 562 I=2,II1
C      DO 562 J=2,JJ1
C  562 UTIL1(I,J)=
C     .2.*((STRMF(I+1,J)-STRMF(I,J))/(P(I  ,J,KK+1)+P(I  ,J-1,KK+1))
C     .   +(STRMF(I-1,J)-STRMF(I,J))/(P(I-1,J,KK+1)+P(I-1,J-1,KK+1))
C     .   +(STRMF(I,J+1)-STRMF(I,J))/(P(I,J  ,KK+1)+P(I-1,J  ,KK+1))
C     .   +(STRMF(I,J-1)-STRMF(I,J))/(P(I,J-1,KK+1)+P(I-1,J-1,KK+1)))
C
C=======================================================================C
C     COMPUTE THE NONDIVERGENT BAROTROPIC VELOCITY FIELD (UTROP) FROM   C
C     THE  STREMFUNCTION AND SUBTRACT THIS FROM THE TOTAL TO FIND THE   C
C     DIVERGENT PART OF THE BAROTROPIC FLOW.                            C
C=======================================================================C
C
      DO 43 J=1,JJ1
       DO 43 I=2,II1
  43  UTROP(I,J)=UTROP(I,J)-(STRMF(I,J)-STRMF(I,J+1))/DEPTHU(I,J)
C
      DO 40 J=2,JJ1
       DO 40 I=1,II1
  40  VTROP(I,J)=VTROP(I,J)-(STRMF(I+1,J)-STRMF(I,J))/DEPTHV(I,J)
C
      DO 41 K=1,KK
        KN=K+NN
         DO 42 J=1,JJ1
          DO 42 I=2,II1
42         U(I,J,KN)=U(I,J,KN)-UTROP(I,J)
            DO 41 J=2,JJ1
             DO 41 I=1,II1
41    V(I,J,KN)=V(I,J,KN)-VTROP(I,J)
C
C=======================================================================C
C                 SMOOTH  U  AND  V  FIELDS IN TIME                     C
C=======================================================================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
      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
22	  CONTINUE
C
C=======================================================================C
C                                                                       C
C                  OUTPUT AND DIAGNOSTIC CALCULATIONS                   C
C                  ==================================                   C
C                                                                       C
C=======================================================================C
C
      IF (INDEX.EQ.0) GO TO 23
      PRINT 100,NSTEP,TIME
100   FORMAT (' T I M E   S T E P'I9,25X'D A Y'F8.1)
C
C --- ENERGY DIAGNOSTICS
C
      PBAR=0.
      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=PBOTTM
      DO 51 I=1,II1
      DO 51 J=1,JJ1
      IF (K.EQ.1)
     .EPOT=EPOT+.5*(STRMF(I,J)*CORIO(I,J)*SCALE/P(I,J,KK+1))**2/G
      IF (K.GT.1)
     .EPOT=EPOT+.5*(P(I,J,K)-PBAR)**2*(THETA(K)-THETA(K-1))/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
C --- CORRECT KINETIC ENERGY BUDGET FOR REMOVAL OF DIVERGENT PART OF MEAN MOTION
      XCONT(K)=XCONT(K)-.5*DP(I,J,KM)*
     .(U(I  ,J,KM)*UTROP(I  ,J)+V(I,J  ,KM)*VTROP(I,J  )
     .+U(I+1,J,KM)*UTROP(I+1,J)+V(I,J+1,KM)*VTROP(I,J+1))/(G*DELT1)
51    WEIGHT=AMIN1(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)
      EKINT=STRESS(K)+XCONT(K)+DISVIS(K)
      PRINT 101,WEIGHT,EPOT,EKIN,EKINT,XCONT(K),STRESS(K),DISVIS(K),K
101   FORMAT (1P' THKN'E10.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    PBAR=PBAR+DP0(K)
      PRINT 101,SUMWGT,SUMPOT,SUMKIN,SUMEKT,SUMXGR,SUMSTR,SUMDIS
C
C --- OUTPUT TO RESTART FILE
C
      NO=20
      REWIND (NO)
      WRITE (NO) NSTEP,DAY,U
      NO=21
      REWIND (NO)
      WRITE (NO) NSTEP,DAY,V
      NO=22
      REWIND (NO)
      WRITE (NO) NSTEP,DAY,DP
C
C --- OUTPUT TO HISTORY FILE
C      
      NO=30
      WRITE (NO) (((U(I,J,K),I=1,II),J=1,JJ),K=1,KK),NSTEP
      NO=31
      WRITE (NO) (((V(I,J,K),I=1,II),J=1,JJ),K=1,KK),NSTEP
      NO=32
      WRITE (NO) (((DP(I,J,K),I=1,II),J=1,JJ),K=1,KK),NSTEP
C
C --- OUTPUT TO AVERAGED FILES
C
C
C
      IREC=IREC+1
      DO 180 K=1,KK
      DO 180 J=1,JJ
      DO 180 I=1,II
      DPM(I,J,K)=DP(I,J,K+NN)+DPM(I,J,K)
      UM(I,J,K)=U(I,J,K+NN)+UM(I,J,K)
 180  VM(I,J,K)=V(I,J,K+NN)+VM(I,J,K)
C
      NO=23
C      REWIND (NO)
C      WRITE (NO) UM,IREC
      NO=24
C      REWIND (NO)
C      WRITE (NO) VM,IREC
      NO=25
C      REWIND (NO)
C      WRITE (NO) DPM,IREC
C
C
23    L=M
      M=N
      N=L
      IF (NSTEP.GE.NSTEP2) STOP '(NORMAL)'
      DELT1=DELT+DELT
      GO TO 15
      END
C
      FUNCTION UVDEP(A,B)
C
C --- FUNCTION SUBROUTINE TO CHOOSE FLUID DEPTH AT U AND V POINTS AS
C --- A FUNCTION OF DEPTH AT NEIGHBORING GRID POINTS
C --- SEE BLECK AND SMITH (JGR, 1990) FOR DISCUSSION OF CHOOSING
C --- MINIMUM RATHER THAN AVERAGE
C
      UVDEP=AMIN1(A,B)
      RETURN
      END
C
          SUBROUTINE POINIT(M,N,L,PBOT,EIG,WSAVE)
C         =======================================                       C
C    INITIALIZATION ROUTINE FOR FFTPSSN:                                C
C    THIS VERSION OF POINIT PERMITS PBOTTOM TO VARY IN Y DIRECTION      C
C
          REAL EIG(L,N), WORK(1), PBOT(1)
           PI = 4.*ATAN(1.)
C--- INITIALIZE FFT ROUTINE
           CALL SINTI(M,WSAVE)
C--- EIGENVALUES OF TRIDIAGONAL MATRIX
           DO 1 J=1,M
           DO 1 K=1,N
    1      EIG(J,K) = -1./PBOT(K) - 1./PBOT(K+1)
     .         -8.*SIN(J*PI/2./(M+1))**2/(2.*UVDEP(PBOT(K),PBOT(K+1)))
C--- GAUSSIAN ELIMINATION OF TRI-DIAGONAL SYSTEM - LEFT HAND SIDE
           DO 2 K = 2,N
           DO 2 J = 1,M
    2      EIG(J,K) = EIG(J,K) - 1./(PBOT(K)**2*EIG(J,K-1))
C
           RETURN
           END
C
       SUBROUTINE POISSON(ZETA,II,M,N,PBOT,EIG,WSAVE)
C                                                                 C
C       THIS SUBROUTINE USES A FAST SINE TRANSFORM ALGORITHM      C
C       TO SOLVE  THE  POISSON  EQUATION:                         C
C                                                                 C
C                      PSI  + PSI  = ZETA                         C
C                         XX     YY                               C
C                                                                 C
C       USING  THE DIRICHLET'S BOUNDARY  CONDITIONS:              C
C                                                                 C
C         PSI(0,Y) = PSI(LX,Y) = PSI(X,0) = PSI(X,LY) = 0         C
C                                                                 C
C       IN A RECTANGULAR DOMAIN  WITH DIMENSIONS LX, LY.          C
C                      ---------------------------                C
C       NOTE THAT THE THE ARRAY ZETA INPUTS THE RIGHT HAND        C
C       SIDE OF THE EQUATION AND  RETURNS THE SOLUTION, PSI.      C
C                                                                 C
C       WSAVE IS AN WORKING ARRAY OF DIMENSION INT(2.5*MAX+15)    C
C           IF II2=JJ2 THEN WSAVE DOES NOT NEED TO BE MODIFIED    C
C             BY CALLING SINTI AGAIN.                             C
C                    ********************                         C
C       THIS SUBROUTINE CALLS SUBROUTINES SINT AND SINTI,  FROM   C
C       NCAR'S PACKAGE: FFTPACK (LINK WITH LIBRARY MYLIBRY.OLB)   C
C 
        DIMENSION ZETA(II,1), EIG(II,1), PBOT(1), WSAVE(1)
C 
          PI = 4.*ATAN(1.)
          N1 = N+1
          M1 = M+1
C=================================================================C
C         PART I - COMPUTES TRANSFORM OF ZETA(I,J)                C
C=================================================================C
C 
          DO 60 K = 1,N
          DO 50 J = 1,M
   50     ZETA(J,K) = ZETA(J,K)/(2.*(M1))
C--- FORWARD TRANSFORM:
   60     CALL SINT(M,ZETA(1,K),WSAVE)
C 
C--- GAUSSIAN ELIMINATION - RIGHT HAND SIDE
          DO 2 K=2,N
          DO 2 J=1,M
    2     ZETA(J,K) = ZETA(J,K)-ZETA(J,K-1)/(PBOT(K)*EIG(J,K-1))
          DO 3 J=1,M
    3     ZETA(J,N) = ZETA(J,N)/EIG(J,N)
          DO 4 K=N-1,1,-1
          DO 4 J=1,M
    4     ZETA(J,K)=(ZETA(J,K)-ZETA(J,K+1)/PBOT(K+1))/EIG(J,K)
C--- BACK TRANSFORM:
          DO 6 K=1,N
          CALL SINT(M,ZETA(1,K),WSAVE)
    6     ZETA(M1,K) = 0.
          RETURN
          END

c=========================================================================
c=========================================================================
c		INPUT FILE NEEDED WITH PARAMETERS
c=========================================================================
c
c 0,5                         ! nstep1, nstep2
c 1152                        ! delt
c 500.e5,4500.e5              ! dp(1),dp(2)
c 1.0,2                       !  wind stress, 1 --> free slip, 2 --> no-slip
c 5.e6                        !  viscosity in cgs units
c
c=========================================================================

