                             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	UTIL3       --> SEA LEVEL ELEVATION, OR VISCOSITY
C       UTIL5       --> KINETIC ENERGY (BAROTROPIC) [g/(s**2)]
C						OR ERG/(CM**2)	
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=201,JJ=II,KK=2,II2=II-2,JJ2=JJ-2,MAX=II2)
      PARAMETER(NW=INT(2.5*MAX+II))
C
C======================================================================C
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)    ,UFLUX(II,JJ)  ,MSTRMF(II,JJ) 
     .    ,VFLUX(II,JJ)  ,UOLD(II,JJ,KK),NPTS(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)  ,UTIL3(II,JJ)  ,UTIL5(II,JJ)
     .    ,U0(II,JJ) , U1(II,JJ) ,U2(II,JJ) ,ULD(II,JJ,KK)
     .    ,V0(II,JJ) , V1(II,JJ) , V2(II,JJ) ,VLD(II,JJ,KK)
     .    ,PSI(II,JJ,KK),VG(II,JJ),VGR(II,JJ),PU(II,JJ,KK+1)
     .    ,PV(II,JJ,KK+1)
c
C
      EQUIVALENCE   (VORT,STRMF)
C     ----------------------------------------------------------------------
C
c	-------------------------------------------------------------
c
C ---           FUNCTIONS TO MIMIC CRAY VECTOR MASK FUNCTIONS
c
              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
C======================================================================C
C                   DEFINE PHYSICAL PARAMETERS                         C
C======================================================================C
C
C --- SUBTRACT FOLLOWING VALUES FROM 1. TO OBTAIN LAYER SPECIFIC VOLUME
c
               DATA THETA/0.,0.003/   
c
C
c	RHO=1./(1-THETA) GIVES THE RELATION BETWEEN RHO AND THETA
c
C
C --- HORIZONTAL GRID SCALE

               DATA SCALE/10.E5/             ! in cm
C
c
C --- PRESSURE AND DP's  ARE GIVEN IN UNITS OF dyn/(cm**2)
C --- EAST-WEST DIRECTION ARRAY OF BOTTOM PRESSURES

               DATA PBOT/JJ*5000.E5/
	       DATA PBOTTM/5000.E5/         

C
C --- GRAVITATIONAL ACCELERATION
               DATA G/980.6/
C
C
C --- VELOCITY TIME SMOOTHING WEIGHTS
               DATA WGT1,WGT2/.75,.125/
C
C
C ---  LINEAR BOTTOM DRAG COEFFICIENT
               DATA DRAG/0.E-7/
C
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/,EPSIL/1.E-9/
C
C
C
C --- 'EPS' IS USED WHEREVER LAYER THICKNESS MUST 
C --- BE BOUNDED AWAY FROM ZERO.
C

c
C --- MODEL IS TO BE INTEGRATED FROM TIME STEP 'NSTEP1' TO 'NSTEP2'
C
        READ(5,*) NDAY1,NDAY2
          type*,'nday1= ', nday1, ' nday2 = ', nday2
	READ(5,*)DELT               ! TIME STEP OF INTEGRATION
        READ(5,*)(DP0(K),K=1,KK)
          type*,'dp0(1) = ',dp0(1),' dp0(2) = ',dp0(2)
        READ(5,*)IFREE
	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
      NSTEP1=NDAY1*86400./DELT+.001
      NSTEP2=NDAY2*86400./DELT+.001
c
C---
	       SCALE2=SCALE**2	
               X=1./SCALE
               II1 = II-1
               JJ1 = JJ-1
C
C
C
c
c
C
c
c
c
C
C======================================================================C
C                         I/O ADDRESSES                                C
C======================================================================C
C
c        open(20,file='u.7',status='new',form='unformatted')
C
c        open(21,file='v.7',status='new',form='unformatted')
C
c        open(22,file='dp.7',status='new',form='unformatted')
C
        open(30,file='u.7',status='new',form='unformatted')
C
        open(31,file='v.7',status='new',form='unformatted')
C
        open(32,file='dp.7',status='new',form='unformatted')
C
c        open(10,file='ul4.1',status='old',form='unformatted')
C
c        open(11,file='vl4.1',status='old',form='unformatted')
C
c        open(12,file='dpl4.1 ',status='old',form='unformatted')
C
C ---------------------------------------------------------------------C
      
C======================================================================C
C                        INITIALIZATION                                C
C======================================================================C
C
C --- CALL INITIALIZATION SUBROUTINE FOR POISSON SOLVER:
C
          CALL POINIT(II-2,JJ-2,II,PBOT,EIG,WSAVE)                   
C
C --- SPECIFY INITIAL AND BOUNDARY CONDITIONS
C
      DO 73 K=1,KK
          DO 73 I=1,II
             DO 73 J=1,JJ
   73        TOPO(I,J)=TOPO(I,J)-(1.-THETA(K))*DP0(K)/G
C
      DO 10 I=1,II
          DO 10 J=1,JJ
C
C --- CORIOLIS PARAMETER
	BETA=2.E-13
          CORIO(I,J)=.93E-4 - FLOAT(I-II1/2)*SCALE*BETA
C
C
C
          UFLUX(I,J)=0.
          VFLUX(I,J)=0.
          STRMF(I,J)=0.
          MONTG(I,J,KK)=0.
          P(I,J,1)=0.
C
             DO 10 K=1,KK
             U(I,J,K   )=0.
             U(I,J,K+KK)=0.
             V(J,J,K   )=0.
             V(I,J,K+KK)=0.
             DP(I,J,K   )=DP0(K)
10           DP(I,J,K+KK)=DP0(K)
C
C======================================================================C
C======================================================================C
C
C
C
C
C
C    INITIAL PRESSURE FIELD for the eddy simulation
C
C    initial COORDINATES OF THE CENTER
      I0=100
      J0=100
C    MAXIMUM INTERFACE DISPLACEMENT AT THE CENTER (in meters*e5)
C    POS. ANTICYCLONE NEG. CYCLONE
      CD=300.E5
C    RADIUS OF MAXIMUM VELOCITY (in km*e5)
      RMAX=65.E5
      DO 3004 J=1,JJ
      DO 3004 I=1,II
3004  P(I,J,2)=DP0(1)
      DO 3071 J=60,140
      DO 3071 I=60,140
C    GAUSSIAN VORTEX
      P(I,J,2)=DP0(1)+CD*EXP(-(SCALE**2*((I-I0)**2+
     .(J-J0)**2))/(2*RMAX**2))
3071    CONTINUE
      DO 3019 J=1,JJ
      DO 3019 I=1,II
      P(I,J,3)=PBOTTM
      DP(I,J,1)=P(I,J,2)
3019    DP(I,J,2)=PBOTTM-DP(I,J,1)
C
C
C   INITIAL VELOCITY FIELD (upper layer only)
C
C
      DO 3064 J=1,JJ1
      DO 3064 I=1,II1
      DO 3064 K=KK-1,1,-1
3064    PSI(I,J,K)=P(I,J,K+1)*(THETA(K+1)-THETA(K))
      K=1
      DO 3047 I=2,II2
      DO 3047 J=2,JJ2
      R=SQRT(FLOAT((I-I0)**2+(J-J0)**2))
      VG(I,J)=G*(THETA(K+1)-THETA(K))*R*SCALE*(DP(I,J,1)-DP0(1))
     ./(RMAX**2*.5*(CORIO(I,J)+CORIO(I+1,J))*1000.)
      DEN=(.5*(CORIO(I,J)+CORIO(I+1,J)))**2-4.*G*(THETA(K+1)-
     .THETA(K))*(DP(I,J,1)-DP0(1))/(1000.*RMAX**2)
      if (den.lt.0.) stop 'negative den'
      VGR(I,J)=-.5*R*SCALE*(-.5*(CORIO(I,J)+CORIO(I+1,J))+SQRT(DEN))
3047   CONTINUE
C
      DO 930 J=2,JJ2
      DO 930 I=2,II2
      IF(J.EQ.J0) GOTO 931
      UI=I
      UJ=J
      RU=SQRT((UI-I0)**2+(J-J0)**2)
      U1(I,J)=.5*(VG(I,J)+VG(I,J))*(UJ-J0)/RU
      U2(I,J)=.5*(VGR(I,J)+VGR(I,J))*(UJ-J0)/RU
C
931    IF(I.EQ.I0) GOTO 930
      VJ=J
      VI=I
      RV=SQRT((I-I0)**2+(VJ-J0)**2)
      V1(I,J)=-.5*(VG(I,J)+VG(I,J))*(VI-I0)/RV
      V2(I,J)=-.5*(VGR(I,J)+VGR(I,J))*(VI-I0)/RV
930    CONTINUE
C
      DO 3048 J=2,JJ2
      DO 3048 I=2,II2
      U0(I,J)=(PSI(I,J-1,K)+PSI(I-1,J-1,K)-
     .PSI(I,J+1,K)-PSI(I-1,J+1,K))/(SCALE*4.*CORIO(I,J))
      V0(I,J)=(PSI(I+1,J,K)+PSI(I+1,J-1,K)-PSI(I-1,J,K)-PSI(I-1,J-1,K))
     ./((CORIO(I,J)+CORIO(I+1,J))*2.*SCALE)
      U1(I,J)=.5*(U1(I,J)+U1(I-1,J))
      U2(I,J)=.5*(U2(I,J)+U2(I-1,J))
      V1(I,J)=.5*(V1(I,J)+V1(I,J-1))
      V2(I,J)=.5*(V2(I,J)+V2(I,J-1))
      U(I,J,1)=U2(I,J)
      V(I,J,1)=V2(I,J)
c      U(I,J,2)=U0(I,J)
c      V(I,J,2)=V0(I,J)
c      DP(I,J,1)=DP0(1)
c      DP(I,J,2)=DP0(2)
3048  CONTINUE

C                      CHOICE OF 
C
C     U0 AND V0 ARE GEOSTROPHIC COMPUTED FROM THE ACTUAL DISPLACEMENT
C     U1 AND V1 ARE GEOSTROPHIC FROM V=G'/F dH/dr
C     U2 AND V2 ARE GRADIENT BALANCE  
C
C                      FOR U,V (DEFINED BEFORE STATEMENT 3048)
C
244   DO 2 K = 1,KK
      DO 2 I = 1,II
      DO 2 J = 1,JJ
      V(I,J,K+KK) = V(I,J,K)
      U(I,J,K+KK) = U(I,J,K)
    2 DP(I,J,K+KK)= DP(I,J,K)
C
C======================================================================C
C --- RESTART FROM HISTORY FILE IF COMPUTATION IS NOT TO START         C
C               FROM REST                                              C
C======================================================================C
C
       DELT1=DELT
       IF(NSTEP1.EQ.0) GO TO 75
c
C
C --- RESTART FROM HISTORY FILE
C
      NO=10
      READ (NO) NSTEP0,DAY,U
      NO=11
      READ (NO) NSTEP0,DAY,V
      NO=12
      READ (NO) NSTEP0,DAY,DP
C
      DELT1=DELT+DELT
C
      WRITE (6,111) NSTEP0,NSTEP1
 111    FORMAT (25X'TIME STEP IN HISTORY FILE -'I7,9X'WANTED -'I7)
C
C
      NSTEP1=NSTEP0
      DAY0=DELT*NSTEP1/86400.
C




C
      IF (MOD(NSTEP0,2).EQ.1) THEN
      L=M
      M=N
      N=L
      ENDIF
C
75       	CONTINUE
      NSTEP=NSTEP1
      WRITE (6,102) NSTEP
 102    FORMAT (/' MODEL INTEGRATION STARTS FROM TIME STEP'I7)

C     
C=======================================================================
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
      NSTEP=NSTEP+1
c
c
      TIME=NSTEP*DELT/86400.
      INDEX=0
      IF (AMOD(TIME+.0001,10.).LT..0002) INDEX=1
C
C
C=======================================================================C
C           CONTINUITY EQUATION (LAYER THICKNESS PREDICTION)            C
C=======================================================================C
C
C     ****************************************************************
C --- ----------------------------------------------------------------
C
C ---    (USING 'FLUX CORRECTED TRANSPORT' ALGORITHM)
C
C --- ----------------------------------------------------------------
C     ****************************************************************
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(DP(I-1,J,KN),DP(I,J,KN),U(I,J,KM))
      DO 12 J=2,JJ1
      DO 12 I=1,II1
12    VFLUX(I,J)=V(I,J,KM)*CVMGP(DP(I,J-1,KN),DP(I,J,KN),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)*.5*(DP(  2,J,KM)+DP(  1,J,KM))-
     .UFLUX(  2,J)
16    UFLUX(II1,J)=U(II1,J,KM)*.5*(DP(II2,J,KM)+DP(II1,J,KM))-
     .UFLUX(II1,J)
      DO 17 I=1,II1
      VFLUX(I,  2)=V(I,  2,KM)*.5*(DP(I,  2,KM)+DP(I,  1,KM))-
     .VFLUX(I,  2)
17    VFLUX(I,JJ1)=V(I,JJ1,KM)*.5*(DP(I,JJ2,KM)+DP(I,JJ1,KM))-
     .VFLUX(I,JJ1)
      DO 20 J=1,JJ1
      DO 20 I=3,II2
   20 UFLUX(I,J)=U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))-
     .UFLUX(I,J)
      DO 21 J=3,JJ2
      DO 21 I=1,II1
   21 VFLUX(I,J)=V(I,J,KM)*.5*(DP(I,J,KM)+DP(I,J-1,KM))-
     .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)/P(I-1,J,KK+1),
     .                            DP(I  ,J,KN)/P(I  ,J,KK+1),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)/P(I,J-1,KK+1),
     .                            DP(I,J  ,KN)/P(I,J  ,KK+1),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    
C ---                   END OF  F C T  CALCULATIONS
C
C     ****************************************************************
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
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 --- SECOND-ORDER FLUXES FOR MOMENTUM EQUATIONS
C
      DO 190 J=1,JJ1
      UFLUX(  2,J)=U(  2,J,KM)*.5*(DP(  2,J,KM)+DP(  1,J,KM))
 190  UFLUX(II1,J)=U(II1,J,KM)*.5*(DP(II2,J,KM)+DP(II1,J,KM))
      DO 92 I=1,II1
      VFLUX(I,  2)=V(I,  2,KM)*.5*(DP(I,  2,KM)+DP(I,  1,KM))
92    VFLUX(I,JJ1)=V(I,JJ1,KM)*.5*(DP(I,JJ2,KM)+DP(I,JJ1,KM))
      DO 93 J=1,JJ1
      DO 93 I=3,II2
93    UFLUX(I,J)=U(I,J,KM)*.5*(DP(I,J,KM)+DP(I-1,J,KM))
      DO 94 J=3,JJ2
      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
C                        VORTICITY COMPUTATION                          C
C=======================================================================C
C                                                                       C
C                                                                       C
C=======================================================================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
C
	if(ifree.eq.2) then 

      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)
C
      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)

	else                                  ! for free slip

      DO 881 J=2,JJ1
      ABSVOR(1,J)=CORIO(1,J)
  881 ABSVOR(II,J)=CORIO(II,J)
C
      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
C                       OOOOOOOOOOOOOOOOOOOOO                           C
C=======================================================================C
C                            U EQUATION                                 C
C=======================================================================C
C
C=======================================================================C
C                THICKNESS FIELD FOR LATERAL DIFFUSION                  C
C=======================================================================C

      DO 412 J=1,JJ1
      PS(1,J)=DP(1,J,KM)
      PS(II,J)=DP(II1,J,KM)
      UOLD(1,J,K)=0.
      UOLD(II,J,K)=0.
C
      DO 412 I=2,II1
      PS(I,J)=.5*(DP(I,J,KM)+DP(I-1,J,KM))
      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=======================================================================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*((MONTG(I,J,K)-MONTG(I-1,J,K)
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 --- 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(.5*(DP(I,J,KM)+DP(I-1,J,KM)),EPS))
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)*.5*(DP(I,J,KM)+DP(I-1,J,KM))*
     .(MONTG(I,J,K)-MONTG(I-1,J,K))/(G*SCALE)
C
      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
     .*DEL2U(I,J)*VISCOS*X*X
     ./AMAX1(.5*(DP(I,J,KM)+DP(I-1,J,KM)),EPS)
61    CONTINUE
C
C=======================================================================C
C                       OOOOOOOOOOOOOOOOOOOOO                           C
C=======================================================================C
C                             V EQUATION                                C
C=======================================================================C
C
      DO 76 I=1,II1
      PS(I,1)=DP(I,1,KM)
      PS(I,JJ)=DP(I,JJ1,KM)
      VOLD(I,1,K)=0.
      VOLD(I,JJ,K)=0.
      DO 76 J=2,JJ1
      PS(I,J)=.5*(DP(I,J,KM)+DP(I,J-1,KM))
      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======================================================================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*((MONTG(I,J,K)-MONTG(I,J-1,K)
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
C --- VISCOUS DIFFUSION
     .-DEL2V(I,J)*X*X*VISCOS
     ./AMAX1(.5*(DP(I,J,KM)+DP(I,J-1,KM)),EPS))
      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)*.5*(DP(I,J,KM)+DP(I,J-1,KM))*
     .(MONTG(I,J,K)-MONTG(I,J-1,K))/(G*SCALE)
C
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,J-1,KM))/G
     .*DEL2V(I,J)*X*X*VISCOS
     ./AMAX1(.5*(DP(I,J,KM)+DP(I,J-1,KM)),EPS)
62    CONTINUE
9     CONTINUE
C
C======================================================================C
C                         BOTTOM DRAG                                  C
C======================================================================C
C
c      KN=KK+NN
c      DO 18 J=1,JJ1
c      DO 18 I=1,II1
c      U(I,J,KN)=U(I,J,KN)*(1.-DRAG*DELT1)
c18    V(I,J,KN)=V(I,J,KN)*(1.-DRAG*DELT1)
c

C
C --- SUBSTITUTE DEPTH-WEIGHTED AVERAGES FOR (U,V) AT MASSLESS GRID POINTS
C
      DO 771 K=1,KK+1
      DO 771 J=1,JJ1
      DO 771 I=2,II1
  771 PU(I,J,K)=.5*(P(I,J,K)+P(I-1,J,K))
C
      DO 772 K=1,KK+1
      DO 772 J=2,JJ1
      DO 772 I=1,II1
  772 PV(I,J,K)=.5*(P(I,J,K)+P(I,J-1,K))
C
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)-EPS1,PV(I,J,K1+1)))
     .-AMIN1(PV(I,J,K  ),AMAX1(PV(I,J,K+1)-EPS1,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  )+EPS1,PV(I,J,K1+1)))
     .-AMAX1(PV(I,J,K+1),AMIN1(PV(I,J,K  )+EPS1,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  )+EPS1,PV(I,J,KK+1)))
     .-AMIN1(PV(I,J,K  ),AMAX1(PV(I,J,K+1)-EPS1,PV(I,J,   1))))
C
      DO 863 J=1,JJ1
C
      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)-EPS1,PU(I,J,K1+1)))
     .-AMIN1(PU(I,J,K  ),AMAX1(PU(I,J,K+1)-EPS1,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  )+EPS1,PU(I,J,K1+1)))
     .-AMAX1(PU(I,J,K+1),AMIN1(PU(I,J,K  )+EPS1,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  )+EPS1,PU(I,J,KK+1)))
     .-AMIN1(PU(I,J,K  ),AMAX1(PU(I,J,K+1)-EPS1,PU(I,J,1   ))))

c
c
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 46 J=1,JJ1
      UTROP(  2,J)=UTROP(  2,J)+U(  2,J,KN)*(DP(  2,J,KN)+DP(  1,J,KN))
46    UTROP(II1,J)=UTROP(II1,J)+U(II1,J,KN)*(DP(II2,J,KN)+DP(II1,J,KN))
      DO 31 J=1,JJ1
      DO 31 I=3,II2
   31 UTROP(I,J)=UTROP(I,J)+U(I,J,KN)*(DP(I,J,KN)+DP(I-1,J,KN))
      DO 47 I=1,II1
      VTROP(I,  2)=VTROP(I,  2)+V(I,  2,KN)*(DP(I,  2,KN)+DP(I,  1,KN))
47    VTROP(I,JJ1)=VTROP(I,JJ1)+V(I,JJ1,KN)*(DP(I,JJ2,KN)+DP(I,JJ1,KN))
      DO 32 J=3,JJ2
      DO 32 I=1,II1
   32 VTROP(I,J)=VTROP(I,J)+V(I,J,KN)*(DP(I,J,KN)+DP(I,J-1,KN))
C
C
      DO 34 J=1,JJ1
      DO 34 I=2,II1
34    UTROP(I,J)=UTROP(I,J)/(P(I,J,KK+1)+P(I-1,J,KK+1))
      DO 35 J=2,JJ1
      DO 35 I=1,II1
35    VTROP(I,J)=VTROP(I,J)/(P(I,J,KK+1)+P(I,J-1,KK+1))
      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=======================================================================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 52.... FOR PBOTTM=CONST., THE RESULTS SHOULD AGREE
C      IF (INDEX.GT.0) WRITE(12,110) (VORT(I,I),I=2,13)
C110   FORMAT (1X,1P12E10.3)

       CALL POISSON(VORT(2,2),II,II2,JJ2,PBOT,EIG,WSAVE)
C
c         DO 111 I = 1,II
c           DO 111 J = 1,JJ
c        STRMAUX(J,I) = STRMF(II+1-I,J)*SCALE/(G*1.E12)
c        IF(NSTEP.LE.NSTEP2/2) GO TO 111
c        MSTRMF( J,I) = MSTRMF(J,I) + STRMAUX(J,I)/FLOAT(NSTEP2/2)
c  111   CONTINUE 
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))
     ./(.5*(P(I,J,KK+1)+P(I-1,J,KK+1)))
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))
     ./(.5*(P(I,J,KK+1)+P(I,J-1,KK+1)))
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
C
22	CONTINUE
c
C -----------------------------------------------------------------------------
C
C --- OUTPUT AND DIAGNOSTIC CALCULATIONS
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
      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 
C --- 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 HISTORY FILE
C
      NO=20
c      REWIND (NO)
c      WRITE (NO) NSTEP,DAY,U
      NO=21
c      REWIND (NO)
c      WRITE (NO) NSTEP,DAY,V
      NO=22
c      REWIND (NO)
c      WRITE (NO) NSTEP,DAY,DP
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
C
C
C
C
 23       L=M
      M=N
      N=L
      IF (NSTEP.GE.NSTEP2) STOP '(NORMAL)'
      DELT1=DELT+DELT
      GO TO 15
      END
C
C=======================================================================C
C                                                                       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=======================================================================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/(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
C=================================================================C
C                                                                 C
       SUBROUTINE POISSON(ZETA,II,M,N,PBOT,EIG,WSAVE)
C      ==============================================             C  
C                                                                 C
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=================================================================C
        
        DIMENSION ZETA(II,1), EIG(II,1), PBOT(1), WSAVE(1)

          PI = 4.*ATAN(1.)
          N1 = N+1
          M1 = M+1
C=================================================================C
C         PART I - COMPUTES TRANSFORM OF ZETA(I,J)                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--- 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 ********************************************************************** 
C                     FFTPACK.FOR                                      C
C         Contains the Fortran codes for FFT subroutines in            C
C                      Library FFTPACK.LIB                             C
c                                                                      c
c	This is a reduced version of what is contained in sint.stuff,  c
c        so that only those subroutines that are being used in the     c
c	POINIT were kept.                                              c
C======================================================================C

      SUBROUTINE SINTI (N,WSAVE)                                         
      DIMENSION       WSAVE(1)                                           
C      LOGICAL Q8Q4                                                       
C      SAVE Q8Q4                                                          
      DATA PI /3.14159265358979/                                         
C      DATA Q8Q4 /.TRUE./                                                 
C      IF (Q8Q4) THEN                                                     
C          CALL Q8QST4 ('LOCLIB', 'FFTPACK', 'SINTI',  'VERSION  4')      
C          Q8Q4 = .FALSE.                                                 
C      ENDIF                                                              
      IF (N .LE. 1) RETURN                                               
      NS2 = N/2                                                          
      NP1 = N+1                                                          
      DT = PI/FLOAT(NP1)                                                 
      DO 101 K=1,NS2                                                     
         WSAVE(K) = 2.*SIN(K*DT)                                         
  101 CONTINUE                                                           
      CALL RFFTI (NP1,WSAVE(NS2+1))                                      
      RETURN                                                             
      END                                                                
c
c
      SUBROUTINE SINT (N,X,WSAVE)                                        
      DIMENSION       X(1)       ,WSAVE(1)                               
C      LOGICAL Q8Q4                                                       
C      SAVE Q8Q4                                                          
C      DATA Q8Q4 /.TRUE./                                                 
C      IF (Q8Q4) THEN                                                     
C          CALL Q8QST4 ('LOCLIB', 'FFTPACK', 'SINT',   'VERSION  4')      
C          Q8Q4 = .FALSE.                                                 
C      ENDIF                                                              
      NP1 = N+1                                                          
      IW1 = N/2+1                                                        
      IW2 = IW1+NP1                                                      
      IW3 = IW2+NP1                                                      
      CALL SINT1(N,X,WSAVE,WSAVE(IW1),WSAVE(IW2),WSAVE(IW3))             
      RETURN                                                             
      END                                                                
c
c
      SUBROUTINE SINT1(N,WAR,WAS,XH,X,IFAC)                              
      DIMENSION WAR(1),WAS(1),X(1),XH(1),IFAC(1)                         
      DATA SQRT3 /1.73205080756888/                                      
      DO 100 I=1,N                                                       
      XH(I) = WAR(I)                                                     
      WAR(I) = X(I)                                                      
  100 CONTINUE                                                           
      IF (N-2) 101,102,103                                               
  101 XH(1) = XH(1)+XH(1)                                                
      GO TO 106                                                          
  102 XHOLD = SQRT3*(XH(1)+XH(2))                                        
      XH(2) = SQRT3*(XH(1)-XH(2))                                        
      XH(1) = XHOLD                                                      
      GO TO 106                                                          
  103 NP1 = N+1                                                          
      NS2 = N/2                                                          
      X(1) = 0.                                                          
      DO 104 K=1,NS2                                                     
         KC = NP1-K                                                      
         T1 = XH(K)-XH(KC)                                               
         T2 = WAS(K)*(XH(K)+XH(KC))                                      
         X(K+1) = T1+T2                                                  
         X(KC+1) = T2-T1                                                 
  104 CONTINUE                                                           
      MODN = MOD(N,2)                                                    
      IF (MODN .NE. 0) X(NS2+2) = 4.*XH(NS2+1)                           
      CALL RFFTF1 (NP1,X,XH,WAR,IFAC)                                    
      XH(1) = .5*X(1)                                                    
      DO 105 I=3,N,2                                                     
         XH(I-1) = -X(I)                                                 
         XH(I) = XH(I-2)+X(I-1)                                          
  105 CONTINUE                                                           
      IF (MODN .NE. 0) GO TO 106                                         
      XH(N) = -X(N+1)                                                    
  106 DO 107 I=1,N                                                       
      X(I) = WAR(I)                                                      
      WAR(I) = XH(I)                                                     
  107 CONTINUE                                                           
      RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RFFTI (N,WSAVE)                                         
      DIMENSION       WSAVE(1)                                           
C      LOGICAL Q8Q4                                                       
C      SAVE Q8Q4                                                          
C      DATA Q8Q4 /.TRUE./                                                 
C      IF (Q8Q4) THEN                                                     
C          CALL Q8QST4 ('LOCLIB', 'FFTPACK', 'RFFTI',  'VERSION  4')      
C          Q8Q4 = .FALSE.                                                 
C      ENDIF                                                              
      IF (N .EQ. 1) RETURN                                               
      CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))                            
      RETURN                                                             
      END                                                                
      SUBROUTINE RFFTI1 (N,WA,IFAC)                                      
      DIMENSION       WA(1)      ,IFAC(1)    ,NTRYH(4)                   
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/                  
      NL = N                                                             
      NF = 0                                                             
      J = 0                                                              
  101 J = J+1                                                            
      IF (J-4) 102,102,103                                               
  102 NTRY = NTRYH(J)                                                    
      GO TO 104                                                          
  103 NTRY = NTRY+2                                                      
  104 NQ = NL/NTRY                                                       
      NR = NL-NTRY*NQ                                                    
      IF (NR) 101,105,101                                                
  105 NF = NF+1                                                          
      IFAC(NF+2) = NTRY                                                  
      NL = NQ                                                            
      IF (NTRY .NE. 2) GO TO 107                                         
      IF (NF .EQ. 1) GO TO 107                                           
      DO 106 I=2,NF                                                      
         IB = NF-I+2                                                     
         IFAC(IB+2) = IFAC(IB+1)                                         
  106 CONTINUE                                                           
      IFAC(3) = 2                                                        
  107 IF (NL .NE. 1) GO TO 104                                           
      IFAC(1) = N                                                        
      IFAC(2) = NF                                                       
      TPI = 6.28318530717959                                             
      ARGH = TPI/FLOAT(N)                                                
      IS = 0                                                             
      NFM1 = NF-1                                                        
      L1 = 1                                                             
      IF (NFM1 .EQ. 0) RETURN                                            
      DO 110 K1=1,NFM1                                                   
         IP = IFAC(K1+2)                                                 
         LD = 0                                                          
         L2 = L1*IP                                                      
         IDO = N/L2                                                      
         IPM = IP-1                                                      
         DO 109 J=1,IPM                                                  
            LD = LD+L1                                                   
            I = IS                                                       
            ARGLD = FLOAT(LD)*ARGH                                       
            FI = 0.                                                      
            DO 108 II=3,IDO,2                                            
               I = I+2                                                   
               FI = FI+1.                                                
               ARG = FI*ARGLD                                            
               WA(I-1) = COS(ARG)                                        
               WA(I) = SIN(ARG)                                          
  108       CONTINUE                                                     
            IS = IS+IDO                                                  
  109    CONTINUE                                                        
         L1 = L2                                                         
  110 CONTINUE                                                           
      RETURN                                                             
      END                                                                
c
c
c
	SUBROUTINE RFFTF (N,R,WSAVE)                                       
      DIMENSION       R(1)       ,WSAVE(1)                               
C      LOGICAL Q8Q4                                                       
C      SAVE Q8Q4                                                          
C      DATA Q8Q4 /.TRUE./                                                 
C      IF (Q8Q4) THEN                                                     
C          CALL Q8QST4 ('LOCLIB', 'FFTPACK', 'RFFTF',  'VERSION  4')      
C          Q8Q4 = .FALSE.                                                 
C      ENDIF                                                              
      IF (N .EQ. 1) RETURN                                               
      CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))                    
      RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC)                                 
      DIMENSION       CH(1)      ,C(1)       ,WA(1)      ,IFAC(1)        
      NF = IFAC(2)                                                       
      NA = 1                                                             
      L2 = N                                                             
      IW = N                                                             
      DO 111 K1=1,NF                                                     
         KH = NF-K1                                                      
         IP = IFAC(KH+3)                                                 
         L1 = L2/IP                                                      
         IDO = N/L2                                                      
         IDL1 = IDO*L1                                                   
         IW = IW-(IP-1)*IDO                                              
         NA = 1-NA                                                       
         IF (IP .NE. 4) GO TO 102                                        
         IX2 = IW+IDO                                                    
         IX3 = IX2+IDO                                                   
         IF (NA .NE. 0) GO TO 101                                        
         CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))                 
         GO TO 110                                                       
  101    CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))                 
         GO TO 110                                                       
  102    IF (IP .NE. 2) GO TO 104                                        
         IF (NA .NE. 0) GO TO 103                                        
         CALL RADF2 (IDO,L1,C,CH,WA(IW))                                 
         GO TO 110                                                       
  103    CALL RADF2 (IDO,L1,CH,C,WA(IW))                                 
         GO TO 110                                                       
  104    IF (IP .NE. 3) GO TO 106                                        
         IX2 = IW+IDO                                                    
         IF (NA .NE. 0) GO TO 105                                        
         CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))                         
         GO TO 110                                                       
  105    CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))                         
         GO TO 110                                                       
  106    IF (IP .NE. 5) GO TO 108                                        
         IX2 = IW+IDO                                                    
         IX3 = IX2+IDO                                                   
         IX4 = IX3+IDO                                                   
         IF (NA .NE. 0) GO TO 107                                        
         CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))         
         GO TO 110                                                       
  107    CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))         
         GO TO 110                                                       
  108    IF (IDO .EQ. 1) NA = 1-NA                                       
         IF (NA .NE. 0) GO TO 109                                        
         CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))                  
         NA = 1                                                          
         GO TO 110                                                       
  109    CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))                 
         NA = 0                                                          
  110    L2 = L1                                                         
  111 CONTINUE                                                           
      IF (NA .EQ. 1) RETURN                                              
      DO 112 I=1,N                                                       
         C(I) = CH(I)                                                    
  112 CONTINUE                                                           
      RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1)                                
      DIMENSION       CH(IDO,2,L1)           ,CC(IDO,L1,2)           ,   
     1                WA1(1)                                             
      DO 101 K=1,L1                                                      
         CH(1,1,K) = CC(1,K,1)+CC(1,K,2)                                 
         CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)                               
  101 CONTINUE                                                           
      IF (IDO-2) 107,105,102                                             
  102 IDP2 = IDO+2                                                       
      DO 104 K=1,L1                                                      
         DO 103 I=3,IDO,2                                                
            IC = IDP2-I                                                  
            TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)                
            TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)                
            CH(I,1,K) = CC(I,K,1)+TI2                                    
            CH(IC,2,K) = TI2-CC(I,K,1)                                   
            CH(I-1,1,K) = CC(I-1,K,1)+TR2                                
            CH(IC-1,2,K) = CC(I-1,K,1)-TR2                               
  103    CONTINUE                                                        
  104 CONTINUE                                                           
      IF (MOD(IDO,2) .EQ. 1) RETURN                                      
  105 DO 106 K=1,L1                                                      
         CH(1,2,K) = -CC(IDO,K,2)                                        
         CH(IDO,1,K) = CC(IDO,K,1)                                       
  106 CONTINUE                                                           
  107 RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)                            
      DIMENSION       CH(IDO,3,L1)           ,CC(IDO,L1,3)           ,   
     1                WA1(1)     ,WA2(1)                                 
      DATA TAUR,TAUI /-.5,.866025403784439/                              
      DO 101 K=1,L1                                                      
         CR2 = CC(1,K,2)+CC(1,K,3)                                       
         CH(1,1,K) = CC(1,K,1)+CR2                                       
         CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))                          
         CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2                                
  101 CONTINUE                                                           
      IF (IDO .EQ. 1) RETURN                                             
      IDP2 = IDO+2                                                       
      DO 103 K=1,L1                                                      
         DO 102 I=3,IDO,2                                                
            IC = IDP2-I                                                  
            DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)                
            DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)                
            DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)                
            DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)                
            CR2 = DR2+DR3                                                
            CI2 = DI2+DI3                                                
            CH(I-1,1,K) = CC(I-1,K,1)+CR2                                
            CH(I,1,K) = CC(I,K,1)+CI2                                    
            TR2 = CC(I-1,K,1)+TAUR*CR2                                   
            TI2 = CC(I,K,1)+TAUR*CI2                                     
            TR3 = TAUI*(DI2-DI3)                                         
            TI3 = TAUI*(DR3-DR2)                                         
            CH(I-1,3,K) = TR2+TR3                                        
            CH(IC-1,2,K) = TR2-TR3                                       
            CH(I,3,K) = TI2+TI3                                          
            CH(IC,2,K) = TI3-TI2                                         
  102    CONTINUE                                                        
  103 CONTINUE                                                           
      RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)                        
      DIMENSION       CC(IDO,L1,4)           ,CH(IDO,4,L1)           ,   
     1                WA1(1)     ,WA2(1)     ,WA3(1)                     
      DATA HSQT2 /.7071067811865475/                                     
      DO 101 K=1,L1                                                      
         TR1 = CC(1,K,2)+CC(1,K,4)                                       
         TR2 = CC(1,K,1)+CC(1,K,3)                                       
         CH(1,1,K) = TR1+TR2                                             
         CH(IDO,4,K) = TR2-TR1                                           
         CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)                               
         CH(1,3,K) = CC(1,K,4)-CC(1,K,2)                                 
  101 CONTINUE                                                           
      IF (IDO-2) 107,105,102                                             
  102 IDP2 = IDO+2                                                       
      DO 104 K=1,L1                                                      
         DO 103 I=3,IDO,2                                                
            IC = IDP2-I                                                  
            CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)                
            CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)                
            CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)                
            CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)                
            CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)                
            CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)                
            TR1 = CR2+CR4                                                
            TR4 = CR4-CR2                                                
            TI1 = CI2+CI4                                                
            TI4 = CI2-CI4                                                
            TI2 = CC(I,K,1)+CI3                                          
            TI3 = CC(I,K,1)-CI3                                          
            TR2 = CC(I-1,K,1)+CR3                                        
            TR3 = CC(I-1,K,1)-CR3                                        
            CH(I-1,1,K) = TR1+TR2                                        
            CH(IC-1,4,K) = TR2-TR1                                       
            CH(I,1,K) = TI1+TI2                                          
            CH(IC,4,K) = TI1-TI2                                         
            CH(I-1,3,K) = TI4+TR3                                        
            CH(IC-1,2,K) = TR3-TI4                                       
            CH(I,3,K) = TR4+TI3                                          
            CH(IC,2,K) = TR4-TI3                                         
  103    CONTINUE                                                        
  104 CONTINUE                                                           
      IF (MOD(IDO,2) .EQ. 1) RETURN                                      
  105 CONTINUE                                                           
      DO 106 K=1,L1                                                      
         TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))                          
         TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))                           
         CH(IDO,1,K) = TR1+CC(IDO,K,1)                                   
         CH(IDO,3,K) = CC(IDO,K,1)-TR1                                   
         CH(1,2,K) = TI1-CC(IDO,K,3)                                     
         CH(1,4,K) = TI1+CC(IDO,K,3)                                     
  106 CONTINUE                                                           
  107 RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)                    
      DIMENSION       CC(IDO,L1,5)           ,CH(IDO,5,L1)           ,   
     1                WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)         
      DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,       
     1-.809016994374947,.587785252292473/                                
      DO 101 K=1,L1                                                      
         CR2 = CC(1,K,5)+CC(1,K,2)                                       
         CI5 = CC(1,K,5)-CC(1,K,2)                                       
         CR3 = CC(1,K,4)+CC(1,K,3)                                       
         CI4 = CC(1,K,4)-CC(1,K,3)                                       
         CH(1,1,K) = CC(1,K,1)+CR2+CR3                                   
         CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3                       
         CH(1,3,K) = TI11*CI5+TI12*CI4                                   
         CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3                       
         CH(1,5,K) = TI12*CI5-TI11*CI4                                   
  101 CONTINUE                                                           
      IF (IDO .EQ. 1) RETURN                                             
      IDP2 = IDO+2                                                       
      DO 103 K=1,L1                                                      
         DO 102 I=3,IDO,2                                                
            IC = IDP2-I                                                  
            DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)                
            DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)                
            DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)                
            DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)                
            DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)                
            DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)                
            DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)                
            DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)                
            CR2 = DR2+DR5                                                
            CI5 = DR5-DR2                                                
            CR5 = DI2-DI5                                                
            CI2 = DI2+DI5                                                
            CR3 = DR3+DR4                                                
            CI4 = DR4-DR3                                                
            CR4 = DI3-DI4                                                
            CI3 = DI3+DI4                                                
            CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3                            
            CH(I,1,K) = CC(I,K,1)+CI2+CI3                                
            TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3                          
            TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3                            
            TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3                          
            TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3                            
            TR5 = TI11*CR5+TI12*CR4                                      
            TI5 = TI11*CI5+TI12*CI4                                      
            TR4 = TI12*CR5-TI11*CR4                                      
            TI4 = TI12*CI5-TI11*CI4                                      
            CH(I-1,3,K) = TR2+TR5                                        
            CH(IC-1,2,K) = TR2-TR5                                       
            CH(I,3,K) = TI2+TI5                                          
            CH(IC,2,K) = TI5-TI2                                         
            CH(I-1,5,K) = TR3+TR4                                        
            CH(IC-1,4,K) = TR3-TR4                                       
            CH(I,5,K) = TI3+TI4                                          
            CH(IC,4,K) = TI4-TI3                                         
  102    CONTINUE                                                        
  103 CONTINUE                                                           
      RETURN                                                             
      END                                                                
c
c
c
      SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)               
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,   
     1                C1(IDO,L1,IP)          ,C2(IDL1,IP),               
     2                CH2(IDL1,IP)           ,WA(1)                      
      DATA TPI/6.28318530717959/                                         
      ARG = TPI/FLOAT(IP)                                                
      DCP = COS(ARG)                                                     
      DSP = SIN(ARG)                                                     
      IPPH = (IP+1)/2                                                    
      IPP2 = IP+2                                                        
      IDP2 = IDO+2                                                       
      NBD = (IDO-1)/2                                                    
      IF (IDO .EQ. 1) GO TO 119                                          
      DO 101 IK=1,IDL1                                                   
         CH2(IK,1) = C2(IK,1)                                            
  101 CONTINUE                                                           
      DO 103 J=2,IP                                                      
         DO 102 K=1,L1                                                   
            CH(1,K,J) = C1(1,K,J)                                        
  102    CONTINUE                                                        
  103 CONTINUE                                                           
      IF (NBD .GT. L1) GO TO 107                                         
      IS = -IDO                                                          
      DO 106 J=2,IP                                                      
         IS = IS+IDO                                                     
         IDIJ = IS                                                       
         DO 105 I=3,IDO,2                                                
            IDIJ = IDIJ+2                                                
            DO 104 K=1,L1                                                
               CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)   
               CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)     
  104       CONTINUE                                                     
  105    CONTINUE                                                        
  106 CONTINUE                                                           
      GO TO 111                                                          
  107 IS = -IDO                                                          
      DO 110 J=2,IP                                                      
         IS = IS+IDO                                                     
         DO 109 K=1,L1                                                   
            IDIJ = IS                                                    
            DO 108 I=3,IDO,2                                             
               IDIJ = IDIJ+2                                             
               CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)   
               CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)     
  108       CONTINUE                                                     
  109    CONTINUE                                                        
  110 CONTINUE                                                           
  111 IF (NBD .LT. L1) GO TO 115                                         
      DO 114 J=2,IPPH                                                    
         JC = IPP2-J                                                     
         DO 113 K=1,L1                                                   
            DO 112 I=3,IDO,2                                             
               C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)                    
               C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)                       
               C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)                          
               C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)                     
  112       CONTINUE                                                     
  113    CONTINUE                                                        
  114 CONTINUE                                                           
      GO TO 121                                                          
  115 DO 118 J=2,IPPH                                                    
         JC = IPP2-J                                                     
         DO 117 I=3,IDO,2                                                
            DO 116 K=1,L1                                                
               C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)                    
               C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)                       
               C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)                          
               C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)                     
  116       CONTINUE                                                     
  117    CONTINUE                                                        
  118 CONTINUE                                                           
      GO TO 121                                                          
  119 DO 120 IK=1,IDL1                                                   
         C2(IK,1) = CH2(IK,1)                                            
  120 CONTINUE                                                           
  121 DO 123 J=2,IPPH                                                    
         JC = IPP2-J                                                     
         DO 122 K=1,L1                                                   
            C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)                             
            C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)                            
  122    CONTINUE                                                        
  123 CONTINUE                                                           
C                                                                        
      AR1 = 1.                                                           
      AI1 = 0.                                                           
      DO 127 L=2,IPPH                                                    
         LC = IPP2-L                                                     
         AR1H = DCP*AR1-DSP*AI1                                          
         AI1 = DCP*AI1+DSP*AR1                                           
         AR1 = AR1H                                                      
         DO 124 IK=1,IDL1                                                
            CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)                            
            CH2(IK,LC) = AI1*C2(IK,IP)                                   
  124    CONTINUE                                                        
         DC2 = AR1                                                       
         DS2 = AI1                                                       
         AR2 = AR1                                                       
         AI2 = AI1                                                       
         DO 126 J=3,IPPH                                                 
            JC = IPP2-J                                                  
            AR2H = DC2*AR2-DS2*AI2                                       
            AI2 = DC2*AI2+DS2*AR2                                        
            AR2 = AR2H                                                   
            DO 125 IK=1,IDL1                                             
               CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)                        
               CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)                     
  125       CONTINUE                                                     
  126    CONTINUE                                                        
  127 CONTINUE                                                           
      DO 129 J=2,IPPH                                                    
         DO 128 IK=1,IDL1                                                
            CH2(IK,1) = CH2(IK,1)+C2(IK,J)                               
  128    CONTINUE                                                        
  129 CONTINUE                                                           
C                                                                        
      IF (IDO .LT. L1) GO TO 132                                         
      DO 131 K=1,L1                                                      
         DO 130 I=1,IDO                                                  
            CC(I,1,K) = CH(I,K,1)                                        
  130    CONTINUE                                                        
  131 CONTINUE                                                           
      GO TO 135                                                          
  132 DO 134 I=1,IDO                                                     
         DO 133 K=1,L1                                                   
            CC(I,1,K) = CH(I,K,1)                                        
  133    CONTINUE                                                        
  134 CONTINUE                                                           
  135 DO 137 J=2,IPPH                                                    
         JC = IPP2-J                                                     
         J2 = J+J                                                               
       DO 136 K=1,L1    
            CC(IDO,J2-2,K) = CH(1,K,J)                                   
            CC(1,J2-1,K) = CH(1,K,JC)                                    
  136    CONTINUE                                                        
  137 CONTINUE                                                           
      IF (IDO .EQ. 1) RETURN                                             
      IF (NBD .LT. L1) GO TO 141                                         
      DO 140 J=2,IPPH                                                    
         JC = IPP2-J                                                     
         J2 = J+J                                                        
         DO 139 K=1,L1                                                   
            DO 138 I=3,IDO,2                                             
               IC = IDP2-I                                               
               CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)                 
               CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)                
               CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)                       
               CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)                      
  138       CONTINUE                                                     
  139    CONTINUE                                                        
  140 CONTINUE                                                           
      RETURN                                                             
  141 DO 144 J=2,IPPH                                                    
         JC = IPP2-J                                                     
         J2 = J+J                                                        
         DO 143 I=3,IDO,2                                                
            IC = IDP2-I                                                  
            DO 142 K=1,L1                                                
               CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)                 
               CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)                
               CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)                       
               CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)                      
  142       CONTINUE                                                     
  143    CONTINUE                                                        
  144 CONTINUE                                                           
      RETURN                                                             
      END                 
c
c
c=========================================================================
c=========================================================================
c		INPUT FILE NEEDED WITH PARAMETERS
c=========================================================================
c
c 0,5
c 1152        
c 5,5                                     out1,nstore
c 10.e7,50.e7                                 dp(1),dp(2)
c 1                          1 --> free slip, 2 --> no-slip
c 1.e6                        viscosity in cgs units
c
c=========================================================================
