C                
      PROGRAM MAIN
C************************************************************************
C           TEST WITH SIMPLE FLAT BOTTOMED BASIN                        *
C                                                                       *
C THIS IS A VERSION OF THE THREE DIMENSIONAL, TIME DEPENDENT,           *
C PRIMITIVE EQUATION, OCEAN MODEL DEVELOPED BY ALLAN BLUMBERG AND ME    *
C WITH SUBSEQUENT CONTRIBUTIONS BY LEO OEY AND OTHERS. TWO REFERENCES   *
C ARE:                                                                  *
C                                                                       *
C     BLUMBERG, A.F. AND G.L. MELLOR, DIAGNOSTIC AND PROGNOSTIC         *
C     NUMERICAL CIRCULATION STUDIES OF THE SOUTH ATLANTIC BIGHT         *
C     J. GEOPHYS. RES. 88, 4579-4592, 1983.                             *
C                                                                       *
C     BLUMBERG, A.F. AND G.L. MELLOR, A DESCRIPTION OF A THREE          *
C     COASTAL OCEAN CIRCULATION MODEL, THREE DIMENSIONAL SHELF          *
C     MODELS, COASTAL AND ESTUARINE SCIENCES, 5, N.HEAPS, ED.,          *
C     AMERICAN GEOPHYSICAL UNION, 1987.                                 *
C                                                                       *
C IN SUBROUTINE PROFQ THE MODEL MAKES USE OF THE TURBULENCE CLOSURE     *
C SUB-MODEL MOST RECENTLY DESCRIBED IN:                                 *
C                                                                       *
C     MELLOR, G.L. AND T. YAMADA, DEVELOPMENT OF A TURBULENCE CLOSURE   *
C     MODEL FOR GEOPHYSICAL FLUID PROBLEMS, REV. GEOPHYS. SPACE PHYS.,  *
C     851-879, 1982.                                                    *
C                                                                       *
C                                          GEORGE MELLOR, 03/87         *
C                                                                       *
C IN THIS VERSION THE HORIZONTAL COORDINATES ARE GENERALIZED CURVI-     *
C LINEAR ORTHOGANAL COORDINATES AND THE MODEL HAS BEEN CONFIGURED       *
C TO RUN IN HALF PRECISION (32 BITS ) ON THE CYBER 205.                 *
C TO REDUCE ROUNDOFF ERROR, SALINTIES AND TEMPERATURES ARE REDUCED      *
C BY 35. AND 10. RESPECTIVELY. 64 BIT PRECISION IS USED LOCALLY         *
C IN SUBROUTINES DENS (WHEREIN THE T AND S OFFSET HAS BEEN TAKEN        *
C INTO ACCOUNT) AND PROFT.                                              *
C                                          GM, 09/87                    *
C
C THIS VERSION HAS BEN CONVERTED TO STANDARD FORTRAN 77 BY STEVE        *
C BRENNER WITH SOME ADDITIONAL CHANGES BY ME IN PROFQ WHICH IS NOW      *
C SOMEWHAT SIMPLER. HALF PRECISION HAS BEEN REMOVED.                    *
C A USER'S GUIDE IS AVAILABLE:       *
C    
C    MELLOR, G.L., A USER'S GUIDE FOR A THREE-DIMENSIONAL,              *
C    NUMERICAL OCEAN MODEL. PRINCETON UNIVERSITY REPORT, 1990.          *
C                                          GM, 04/90                    *
C                                                                       *
C Previous to this date the code was called pmod.f and all changes      *
C were recorded in pmod.change. Henceforth, the code will be called     * 
C pom.f and changes will be recorded in pom.change.                     *
C                                          GM, 11/92                    * 
C                                                                       *
C THE LAST CODE CHANGE, as recorded in pom.change, WAS ON NOV. 1, 1994  *
C************************************************************************
C
         INCLUDE 'comblk.h'    
         INCLUDE 'params'    
      DIMENSION UTB(IM,JM),VTB(IM,JM),UTF(IM,JM),VTF(IM,JM),
     1     ADVX(IM,JM,KB),ADVY(IM,JM,KB),ADVUA(IM,JM),ADVVA(IM,JM),
     2     TSURF(IM,JM),SSURF(IM,JM),
     2     DRHOX(IM,JM,KB),DRHOY(IM,JM,KB),TRNU(IM,JM),TRNV(IM,JM),
     3     TMEAN(IM,JM,KB),SMEAN(IM,JM,KB),ADVUU(IM,JM),ADVVV(IM,JM),
     4     SWRAD(IM,JM),RMEAN(IM,JM,KB),
     6     TEMPDP(1000),TEMPIS(1000),TEMPWC(1000),TEMPWA(1000)
c=====================================================================
c      DEFINE SEDIMENT TRANSPORT VARIABLES
C===================================================================
      DIMENSION TTT1(IM,JM,KB),TTT2(IM,JM,KB)
      DIMENSION srho1(im,jm),srho2(im,jm)
      DIMENSION TEMPT(1000),RB1(im,jm,kb),RB2(im,jm,kb)
      DIMENSION RRMEAN(IM,JM,KB),sink1(im,jm,kb)
c      DIMENSION RRMEAN(IM,JM,KB)
C--------------------------------------------------------------------
      REAL ISPI,ISP2I
C     NAMELIST/PRMTR/MODE,TIME,DTI,ISPLIT,NREAD,IPRTD1,ISWTCH,IPRTD2,
C    1       IDAYS,ISPADV,HORCON,TPRNU,SMOTH
      DATA PI/3.1416E0/,RAMP/1.E0/,SMALL/1.E-10/,TIME0/0.E0/
      DATA BETA/1.98E-11/,GRAV/9.806E0/,UMOL/2.E-5/
	
	OPEN(40,FORM='UNFORMATTED',STATUS='NEW')
c	OPEN(70,FORM='UNFORMATTED',STATUS='OLD')
	OPEN(80,FORM='UNFORMATTED',STATUS='NEW')
	OPEN(81,FORM='FORMATTED',STATUS='NEW')
	OPEN(85,FORM='UNFORMATTED',STATUS='NEW')
	OPEN(90,FORM='UNFORMATTED',STATUS='NEW')

      IINT=0
      TIME=0.
      DTI=30
      isplit=10
      ISPADV=5
      TPRNU=1.0  
      SMOTH=0.1
      HORCON=.2
      MODE=3
      NBC=1
      TBIAS=10.
      SBIAS=35.
C
      IPRTD1=100
      ISWTCH=22
      IPRTD2=1
C--------------------------------------------------------------------------
C     READ NAMELIST (DISCONNECTED HERE) FOR PROBLEM PARAMETERS: 
C       MODE = 2; 2-D CALCULATION (BOTTOM STRESS CALCULATED IN ADVAVE)
C              3; 3-D CALCULATION (BOTTOM STRESS CALCULATED IN PROFU,V)
C              4; 3-D CALCULATION WITH T AND S HELD FIXED
C       NBC = TYPE OF THERMO. B.C.; SEE SUBROUTINE PROFT
C       DTI = INTERNAL TIME STEP
C       DTE = EXTERNAL TIME STEP
C       ISPLIT = DTI/DTE
C       NREAD=0, NO RESTART INPUT FILE; NREAD=1, RESTART
C       IPRTD1 = PRINT INTERVAL IN DAYS; AFTER ISWTCH DAYS, PRINT
C                  INTERVAL = IPRTD2
C       IDAYS = LENGTH OF RUN IN DAYS
C       ISPADV = STEP INTERVAL WHERE EXTERNAL MODE ADVECTIVE TERMS ARE
C                  NOT UPDATED
C       HORCON = CONSTANT IN SMAGORINSKY HORIZONTAL DIFFUSIVITY
C       TPRNU = HORIZONTAL DIFFUSIVITY PRANDTL NUMBER
C       SMOTH = CONSTANT IN TIME SMOOTHER TO PREVENT SOLUTION SPLITTING
C
C     READ(5,PRMTR)
C--------------------------------------------------------------------------
      DTE=DTI/FLOAT(ISPLIT)
      DTE2=DTE*2
      DTI2=DTI*2
      IPRINT=IPRTD1*24*3600/INT(DTI)
      IPRINT2=IPRTD2*24*3600/INT(DTI)
      IPRINT3=1*3600/INT(DTI)
      ISWTCH=ISWTCH*24*3600/INT(DTI)
      IEND=IDAYS*24*3600/INT(DTI)+1     
      RBIAS=0.002
C:: Sediment entrainment velocity
        svel1=1.e-5

      WRITE(6,7030) MODE,DTE,DTI,IEND,ISPLIT,ISPADV,IPRINT,SMOTH,
     1   HORCON,TPRNU   
 7030 FORMAT(//,' MODE =',I3,
     1 '  DTE =',F7.2,'   DTI =',F9.1,'     IEND =',I6,
     2  '   ISPLIT =',I6,'    ISPADV =',I6,'   IPRINT =',I6,/,
     3  '   SMOTH =',F7.2,'   HORCON =',F7.3,'   TPRNU ='F7.3,/)
      WRITE(6,'('' NREAD ='',I5,//)') NREAD
C
C----------------------------------------------------------------------
C             ESTABLISH PROBLEM CHARACTERISTICS
C          ****** ALL UNITS IN M.K.S. SYSTEM ******
C      F,BLANK AND B REFERS TO FORWARD,CENTRAL AND BACKWARD TIME LEVELS.
C----------------------------------------------------------------------
C
      DAYI=1.E0/86400.E0
      ISPI=1.E0/FLOAT(ISPLIT)
      ISP2I=1.E0/(2.E0*FLOAT(ISPLIT))
C
C     SET UP  GRID AND INITIAL CONDITIONS FOR STAND-ALONE TEST RUN
      CALL DEPTH(Z,ZZ,DZ,DZZ,DZR,KB,KBM1)
      DO 11 K=1,KB
      DO 11 J=1,JM
      DO 11 I=1,IM
      Q2B(I,J,K)=1.E-8   
      Q2LB(I,J,K)=1.E-8   
      KH(I,J,K)=2.E-4
      KM(I,J,K)=KH(I,J,K)
C:: Settling velocity Ws
      sink1(i,j,k)=-1.1E-4
  11  AAM(I,J,K)=500.
      DO 12 J=1,JM
      DO 12 I=1,IM
      AAM2D(I,J)=AAM(I,J,1)    
	WRSURF(I,J)=0.0
      DX(I,J)=5.E2   
      DY(I,J)=5.E2
c      if (i.ge.40) dx(i,j)=1.0e3+float(i-40)*5.0e2
      ART(I,J)=DX(I,J)*DY(I,J)
  12  CONTINUE


      DO 13 J=2,JM
      DO 13 I=2,IM
      ARU(I,J)=.25*(DX(I,J)+DX(I-1,J))*(DY(I,J)+DY(I-1,J))
  13  ARV(I,J)=.25*(DX(I,J)+DX(I,J-1))*(DY(I,J)+DY(I,J-1))
      DO 121 J=1,JM
      ARU(1,J)=ARU(2,J)
121   ARV(1,J)=ARV(2,J)
      DO 122 I=1,IM
      ARU(I,1)=ARU(I,2)
122      ARV(I,1)=ARV(I,2)
      DO 14 I=1,IM
      DO 14 J=1,JM
      COR(I,J)=-0.8E-4
        h(i,j)=20.0
  14  CONTINUE

      DO 18 I=1,IM
      DO 18 J=1,JM
c      H(IM,J)=1.E0
      H(1,J)=1.E0
      H(I,1)=1.E0
      H(I,JM)=1.E0
  18  CONTINUE

C###############################


      DO 15 K=1,KB
      DO 15 J=1,JM
      DO 15 I=1,IM
       TINTER=25.0
      TB(I,J,K)=TINTER-TBIAS  
      SB(I,J,K)=35.0-SBIAS  
  15  CONTINUE
      DO 17 K=1,KB
      DO 17 J=1,JM
      DO 17 I=1,IM
C:: RB1 & 2 are sediment concentration C_fine and C_coarse
C:: In Wang JPO, 2002, only RB1 is used.
C::  XH Wang, March 2005
         RB1(I,J,K)=0.0
         RB2(I,J,K)=0.0
         RR1(I,J,K)=0.
         RR2(I,J,K)=0.
	RMEAN(I,J,K)=0.
      TMEAN(I,J,K)=TB(I,J,K)
  17  SMEAN(I,J,K)=SB(I,J,K)

      DO K=1,KB
      DO J=1,JM
      SBE(J,K)=0.0
      TBE(J,K)=TB(IM,4,K)
      END DO
      END DO

	DO K=1,KB
	DO J=1,JM
	DO I=1,IM
	RRMEAN(I,J,K)=0.0
	END DO
	END DO
	END DO

C
C------------------------------------------------------------------------
C     READ IN GRID DATA AND INITIAL CONDITIONS
C------------------------------------------------------------------------
C     REWIND(40)
C     REWIND(60)
C     READ(40) KBB,Z,ZZ,DZ,DZZ,IMM,JMM,ALON,ALAT,DX,DY,H,COR,
C    1      ART,ARU,ARV,TMEAN,SMEAN,RMEAN,TBN,TBE,TBS,SBN,SBE,SBS,
C    2      ELN,ELE,ELS,VABN,UABE,VABS
      TIME0=0.
C     READ(50) TIME0,UB,VB,UAB,VAB,TB,SB,ELB,Q2B,Q2LB,KM,KH
C*****************************************************************************
      DO 16 J=1,JM
      DO 16 I=1,IM
      srho1(i,j)=0.0
      srho2(i,j)=0.0
      TSURF(I,J)=TB(I,J,1)
  16  SSURF(I,J)=SB(I,J,1)
      DO 20 K=1,KBM1
  20  DZR(K)=1./DZ(K)
      PERIOD=DAYI*(2.E0*PI)/ABS(COR(IM/2,JM/2))
      DO 30 J=1,JM
      DO 30 I=1,IM
      FSM(I,J)=1.E0
      DUM(I,J)=1.E0
      DVM(I,J)=1.E0
      IF(H(I,J).GT.1.E0) GO TO 30
      FSM(I,J)=0.E0
      DUM(I,J)=0.E0
      DVM(I,J)=0.E0
 30   CONTINUE
      DO 35 J=1,JMM1
      DO 35 I=1,IM
      IF(FSM(I,J).EQ.0.E0.AND.FSM(I,J+1).NE.0.E0) DVM(I,J+1)=0.E0
 35   CONTINUE
      DO 40 J=1,JM
      DO 40 I=1,IMM1
      IF(FSM(I,J).EQ.0.E0.AND.FSM(I+1,J).NE.0.E0) DUM(I+1,J)=0.E0
 40   CONTINUE
C
      DO 41 J=1,JM
      DUM(1,J)=0.E0
  41  DVM(1,J)=0.E0
C A very empirical specification of the roughness parameter follows
      DO 45 J=1,JM
      DO 45 I=1,IM
      Z0B=.0010
      CBCMIN=small
  45  CBC(I,J)=MAX(CBCMIN,.16E0/LOG((ZZ(KBM1)-Z(KB))*H(I,J)
     1        /Z0B)**2)*FSM(I,J)
c  45  CBC(I,J)=.0025E0
C
      IP=IM
      ISK=1

       WRITE(40) KB,Z,ZZ,DZ,DZZ,IM,JM,ALON,ALAT,DX,DY,H,FSM,DUM,DVM

        CALL PRXY('  TOPOGRAPHY   ',TIME, H,IP,ISK,JM,1,1.E0)
c        CALL PRXY(' FSM            ',TIME,FSM,IP,ISK,JM,1,1.E0)
c        CALL PRXY(' DUM            ',TIME,DUM,IP,ISK,JM,1,1.E0)
c       CALL PRXY(' DVM            ',TIME,DVM,IP,ISK,JM,1,1.E0)
        CALL PRXY(' CBC            ',TIME,CBC,IP,ISK,JM,1,0.E0)

c      call prxyz('rb10',time,rb1,im,1,jm,1,kbm1,0.)
c      call prxyz('rb20',time,rb2,im,1,jm,1,kbm1,0.)

      VTOT=0.E0
      TTOT=0.E0
      STOT1=0.E0
      STOT2=0.E0
      RTOT1=0.E0
      RTOT2=0.E0
      ATOT=0.E0
      DO  J = 1, JM
      DO  I = 1, IM
        RTOT1=RTOT1+srho1(i,j)*ART(I,J)*FSM(I,J)
        RTOT2=RTOT2+srho2(i,j)*ART(I,J)*FSM(I,J)
c      ATOT=ATOT+ART(I,J)*FSM(I,J)
        end do
        end do
      DO  K=1,KBM1
      DO  J=1,JM
      DO  I=1,IM
      DVTOT=DX(I,J)*DY(I,J)*H(I,J)*DZ(K)*FSM(I,J)
      VTOT=VTOT+DVTOT
      STOT1=STOT1+(RB1(I,J,K)+RBIAS)*DVTOT
      STOT2=STOT2+(RB2(I,J,K)+RBIAS)*DVTOT
c      RTOT=RTOT+RB1(I,J,K)*DVTOT
        end do
        end do
        end do
c==========================================================================
c     CALCULATION OF THE BASIN AVERAGED TEMPERATURE AND SALINITY
c==========================================================================
      TTOT1=STOT1/VTOT
      TTOT2=STOT2/VTOT
c      RTOT=RTOT/ATOT
      ETOT1=STOT1+RTOT1
      ETOT2=STOT2+RTOT2
      WRITE(6,'('' ETOT1 =,'',E20.8,'' TTOT1 ='',E20.8,
     1 '' STOT1 ='',E20.8,'' RTOT1='',E20.8)') ETOT1,TTOT1,STOT1,RTOT1
      WRITE(6,'('' ETOT2 =,'',E20.8,'' TTOT2 ='',E20.8,
     1 '' STOT2 ='',E20.8,'' RTOT2='',E20.8)') ETOT2,TTOT2,STOT2,RTOT2

C
        DO 47 J=1,JM
        DO 47 I=1,IM
   47   TPS(I,J)=0.5E0/SQRT(1.E0/DX(I,J)**2+1.E0/DY(I,J)**2)
     1         /SQRT(GRAV*H(I,J))*FSM(I,J)
c        CALL PRXY(' CFL TIME STEP   ',TIME,TPS,IP,ISK,JM,1,1.E0)
cc      CALL PRXY(' COR           ',TIME,COR,IP,ISK,JM,1,0.E0)
C
      DO 48 I=1,IM
   48 COVRHN(I)=.1E0*SQRT(GRAV/H(I,JMM2))
      DO 49 J=1,JM
   49 COVRHE(J)=.1E0*SQRT(GRAV/H(IMM1,J))
      DO 50 J=1,JM
      DO 50 I=1,IM
      L(I,J,1)=0.E0
   50 L(I,J,KB)=0.E0
      DO 51 J=1,JM
      DO 51 I=1,IM
      UA(I,J)=UAB(I,J)
      VA(I,J)=VAB(I,J)
      EL(I,J)=ELB(I,J)
      ET(I,J)=ETB(I,J)
      ETF(I,J)=ET(I,J)
      D(I,J)=H(I,J)+EL(I,J)
   51 DT(I,J)=H(I,J)+ET(I,J)
      DO 52 K=1,KB
      DO 52 J=1,JM
      DO 52 I=1,IM
      L(I,J,K)=Q2LB(I,J,K)/(Q2B(I,J,K)+SMALL)
      KQ(I,J,K)=0.2E0*L(I,J,K)*SQRT(Q2B(I,J,K))
      Q2(I,J,K)=Q2B(I,J,K)
      Q2L(I,J,K)=Q2LB(I,J,K)
      T(I,J,K)=TB(I,J,K)
      S(I,J,K)=SB(I,J,K)
      U(I,J,K)=UB(I,J,K)
   52 V(I,J,K)=VB(I,J,K)
C
      CALL DENS(S,T,RR1,RR2,RHO)
C
C
      DO 80 K=1,KBM1
      DO 80 J=1,JM
      DO 80 I=1,IM
        if (fsm(i,j).gt.0.5) then
        rhosum=rhosum+rho(i,j,k)
        icount=icount+1
        end if
80      continue
        rhosum=rhosum/float(icount)

      do k=1,kb
      do j=1,jm
      do i=1,im
      RMEAN(I,J,K)=rhosum
      end do
      end do
      end do

C
      CALL BAROPG(DRHOX,DRHOY,RMEAN)
C
C  THE FOLLOWING DATA ARE NEEDED FOR A SEAMLESS RESTART.
      DO 10 N=1,NREAD
      READ(70) TIME0,
     1  WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB,
     2  UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADVUU,ADVVV,ADVUA,ADVVA,
     3  KM,KH,KQ,L,Q2,Q2B,AAM,Q2L,Q2LB
  10  CONTINUE

      DO 81 J=1,JM
      DO 81 I=1,IM
      D(I,J)=H(I,J)+EL(I,J)
   81 DT(I,J)=H(I,J)+ET(I,J)
C
      TIME=TIME0
c      CALL PRXYZ(' TB    ',TIME,TB,IP,ISK,JM,1,KB,0.E0)
c      CALL PRXZ('TB' ,TIME,TB,IM,ISK,JM,35,45,KB,.01,DT,ZZ)
c     CALL PRXYZ('RHO         ',TIME,RHO,IP,ISK,JM,1,KB,0.E0)
c     CALL PRXYZ('DRHOX     ',TIME,DRHOX,IP,ISK,JM,1,KB,0.E0)
c     CALL PRXYZ('DRHOY     ',TIME,DRHOY,IP,ISK,JM,1,KB,0.E0)
c     CALL PRXY(' TRNU      ',TIME,TRNU,IP,ISK,JM,1,0.E0)
c     CALL PRXY(' TRNV      ',TIME,TRNV,IP,ISK,JM,1,0.E0)
c     CALL PRXY('  FREE SURFACE ',TIME, ELB,IP,ISK,JM,1,1.E-3)
c     CALL FINDPSI
C     CALL PRXY('  VAB          ',TIME, VAB,IP,ISK,JM,1,0.E0)
C***********************************************************************
C                                                                      *
C                  BEGIN NUMERICAL INTEGRATION                         *
C                                                                      *
C***********************************************************************
C
C
                     DO 9000 IINT=1,IEND
C
      TIME=DAYI*DTI*FLOAT(IINT)+TIME0
      RAMP=TIME/PERIOD
c      IF(RAMP.GT.0.5E0) RAMP=1.E0
c      RAMP=1.E0
c      IF (TIME.GT.5.) MODE=3
c     write(6,'('' IINT,TIME ='',I5,F9.2)') IINT,TIME      
C------------------------------------------------------------------------
C     SET TIME DEPENDENT, SURFACE AND LATERAL BOUNDARY CONDITIONS.
C     THE LATTER WILL BE USED IN SUBROUTINE BCOND .
C------------------------------------------------------------------------
C    INTRODUCE SIMPLE WIND STRESS, VALUE IS NEGATIVE FOR WESTERLY WIND
C  The following wind stress has been tapered along the boundary to suppress 
C  numerically induced oscilations near the boundary (Jamart&Ozer,
C  J.G.R.,91, 10621-10631)

	time=time-time0

	J1=2
	J2=JMM1

      IF(MODE.NE.2) THEN
        CALL ADVCT(ADVX,ADVY)
        CALL BAROPG(DRHOX,DRHOY,RMEAN)
C**********************************************************************
C     HOR VISC = HORCON*DX*DY*SQRT((DU/DX)**2+(DV/DY)**2
C                                 +.5*(DU/DY+DV/DX)**2)
C**********************************************************************
C  If MODE.EQ.2 then initial values of AAM2D are used. If one wishes 
C  Smagorinsky lateral viscosity and diffusion for an external mode 
C  calculation, then appropiate code can be adapted from that below
C  and installed after s.n 102 and before s.n. 5000 in subroutine ADVAVE.
        DO 95 K=1,KBM1
        DO 95 J=2,JMM1
        DO 95 I=2,IMM1
        AAM(I,J,K)=HORCON*DX(I,J)*DY(I,J)
     1          *SQRT( ((U(I+1,J,K)-U(I,J,K))/DX(I,J))**2
     2                 +((V(I,J+1,K)-V(I,J,K))/DY(I,J))**2
     3 +.5E0*(.25E0*(U(I,J+1,K)+U(I+1,J+1,K)-U(I,J-1,K)-U(I+1,J-1,K))
     4              /DY(I,J)
     5       +.25E0*(V(I+1,J,K)+V(I+1,J+1,K)-V(I-1,J,K)-V(I-1,J+1,K))
     6              /DX(I,J)) **2)
       if (i.gt.110) aam(i,j,k)=500.00
   95   CONTINUE
C --  Form vertical averages of 3-D fields for use in external mode --
        DO 96 J=1,JM
        DO 96 I=1,IM
        ADVUU(I,J)=0.
        ADVVV(I,J)=0.
        TRNU(I,J)=0.
        TRNV(I,J)=0.
        AAM2D(I,J)=0.
   96   CONTINUE    
        DO 100 K=1,KBM1
        DO 100 J=1,JM
        DO 100 I=1,IM
        ADVUU(I,J)=ADVUU(I,J)+ADVX(I,J,K)*DZ(K)
        ADVVV(I,J)=ADVVV(I,J)+ADVY(I,J,K)*DZ(K)
        TRNU(I,J)=TRNU(I,J)+DRHOX(I,J,K)*DZ(K)
        TRNV(I,J)=TRNV(I,J)+DRHOY(I,J,K)*DZ(K)
        AAM2D(I,J)=AAM2D(I,J)+AAM(I,J,K)*DZ(K)
  100   CONTINUE
C ----------------------------------------------------------------------
        CALL ADVAVE(ADVUA,ADVVA,MODE)
        DO 87 J=1,JM
        DO 87 I=1,IM
        ADVUU(I,J)=ADVUU(I,J)-ADVUA(I,J)
   87   ADVVV(I,J)=ADVVV(I,J)-ADVVA(I,J)
C
      ENDIF
C------------------------------------------------------------------------
      DO 399 J=1,JM
      DO 399 I=1,IM
  399 EGF(I,J)=EL(I,J)*ISPI
      DO 400 J=2,JM
      DO 400 I=2,IM
      UTF(I,J)=UA(I,J)*(D(I,J)+D(I-1,J))*ISP2I
 400  VTF(I,J)=VA(I,J)*(D(I,J)+D(I,J-1))*ISP2I
C********** BEGIN EXTERNAL MODE ***************************************
                      DO 8000 IEXT=1,ISPLIT
      DO 405 J=2,JM
      DO 405 I=2,IM
      FLUXUA(I,J)=.25E0*(D(I,J)+D(I-1,J))*(DY(I,J)+DY(I-1,J))*UA(I,J)
 405  FLUXVA(I,J)=.25E0*(D(I,J)+D(I,J-1))*(DX(I,J)+DX(I,J-1))*VA(I,J)
C
      DO 410 J=2,JMM1
      DO 410 I=2,IMM1
  410 ELF(I,J)=ELB(I,J)
     1    -DTE2*(FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J))
     2                    /ART(I,J)
C
      CALL BCOND(1)
C
      IF(MOD(IEXT,ISPADV).EQ.0) CALL ADVAVE(ADVUA,ADVVA,MODE)
C  Note that ALPHA = 0. is perfectly acceptable. The value, ALPHA = .225
C  permits a longer time step.
      ALPHA=0.225     
      DO 420 J=2,JMM1
      DO 420 I=2,IMM1
  420 UAF(I,J)=ADVUU(I,J)+ADVUA(I,J)   
     1    -ARU(I,J)*.25*(  COR(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J))
     2              +COR(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) )
     3         +.25E0*GRAV*(DY(I,J)+DY(I-1,J))*(D(I,J)+D(I-1,J))
     4             *( (1.E0-2.E0*ALPHA)*(EL(I,J)-EL(I-1,J))
     4            +ALPHA*(ELB(I,J)-ELB(I-1,J)+ELF(I,J)-ELF(I-1,J)) )
     5               +TRNU(I,J)   
     6      +ARU(I,J)*( WUSURF(I,J)-WUBOT(I,J)   )
      DO 425 J=2,JMM1
      DO 425 I=2,IMM1
  425 UAF(I,J)=
     1         ((H(I,J)+ELB(I,J)+H(I-1,J)+ELB(I-1,J))*ARU(I,J)*UAB(I,J)
     2                -4.E0*DTE*UAF(I,J))
     3        /((H(I,J)+ELF(I,J)+H(I-1,J)+ELF(I-1,J))*ARU(I,J))
      DO 430 J=2,JMM1
      DO 430 I=2,IMM1
 430  VAF(I,J)=ADVVV(I,J)+ADVVA(I,J)  
     1    +ARV(I,J)*.25*(  COR(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J))
     2               +COR(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) )
     3         +.25E0*GRAV*(DX(I,J)+DX(I,J-1))*(D(I,J)+D(I,J-1))
     4             *( (1.E0-2.E0*ALPHA)*(EL(I,J)-EL(I,J-1))
     4            +ALPHA*(ELB(I,J)-ELB(I,J-1)+ELF(I,J)-ELF(I,J-1)) )
     5               +TRNV(I,J)  
     6    + ARV(I,J)*( WVSURF(I,J)-WVBOT(I,J)   )
      DO 435 J=2,JMM1
      DO 435 I=2,IMM1
  435 VAF(I,J)=
     1        ((H(I,J)+ELB(I,J)+H(I,J-1)+ELB(I,J-1))*VAB(I,J)*ARV(I,J)
     2              -4.E0*DTE*VAF(I,J))
     3       /((H(I,J)+ELF(I,J)+H(I,J-1)+ELF(I,J-1))*ARV(I,J))
      CALL BCOND(2)
C
      IF(IEXT.LT.(ISPLIT-2)) GO TO 440
      IF(IEXT.EQ.(ISPLIT-2))THEN
        DO 431 J=1,JM
        DO 431 I=1,IM
  431   ETF(I,J)=.25*SMOTH*ELF(I,J)
       ENDIF
      IF(IEXT.EQ.(ISPLIT-1)) THEN
        DO 432 J=1,JM
        DO 432 I=1,IM
  432   ETF(I,J)=ETF(I,J)+.5*(1.-.5*SMOTH)*ELF(I,J)
       ENDIF
      IF(IEXT.EQ.(ISPLIT-0)) THEN
        DO 433 J=1,JM
        DO 433 I=1,IM
  433   ETF(I,J)=(ETF(I,J)+.5*ELF(I,J))*FSM(I,J)
       ENDIF
 440  CONTINUE
C
C  TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP
C
      VAMAX=0.       
      DO 442 J=1,JM
      DO 442 I=1,IM
      IF(ABS(VAF(I,J)).GE.VAMAX)VAMAX=ABS(VAF(I,J))   
  442 CONTINUE
      IF(VAMAX.GT.100.) GO TO 9001
C    
C
C       APPLY FILTER TO REMOVE TIME SPLIT. RESET TIME SEQUENCE.
      DO 445 J=1,JM
      DO 445 I=1,IM
      UA(I,J)=UA(I,J)+.5E0*SMOTH*(UAB(I,J)-2.E0*UA(I,J)+UAF(I,J))
      VA(I,J)=VA(I,J)+.5E0*SMOTH*(VAB(I,J)-2.E0*VA(I,J)+VAF(I,J))
      EL(I,J)=EL(I,J)+.5E0*SMOTH*(ELB(I,J)-2.E0*EL(I,J)+ELF(I,J))
      ELB(I,J)=EL(I,J)
      EL(I,J)=ELF(I,J)
      D(I,J)=H(I,J)+EL(I,J)
      UAB(I,J)=UA(I,J)
      UA(I,J)=UAF(I,J)
      VAB(I,J)=VA(I,J)
      VA(I,J)=VAF(I,J)
  445 CONTINUE
C
      IF(IEXT.EQ.ISPLIT) GO TO 8000
      DO 450 J=2,JM
      DO 450 I=2,IM
      EGF(I,J)=EGF(I,J)+EL(I,J)*ISPI
      UTF(I,J)=UTF(I,J)+UA(I,J)*(D(I,J)+D(I-1,J))*ISP2I
 450  VTF(I,J)=VTF(I,J)+VA(I,J)*(D(I,J)+D(I,J-1))*ISP2I
C
C
 8000                    CONTINUE
C---------------------------------------------------------------------
C          END EXTERNAL (2-D) MODE CALCULATION
C     AND CONTINUE WITH INTERNAL (3-D) MODE CALCULATION
C---------------------------------------------------------------------
      IF(IINT.EQ.1.AND.TIME0.EQ.0.) GO TO 8200
      IF(MODE.EQ.2) GO TO 8200
C---------------------------------------------------------------------
C      ADJUST U(Z) AND V(Z) SUCH THAT
C      VERTICAL AVERAGE OF (U,V) = (UA,VA)
C---------------------------------------------------------------------
C
      DO 299 J=1,JM
      DO 299 I=1,IM
  299 TPS(I,J)=0.E0
      DO 300 K=1,KBM1
      DO 300 J=2,JMM1
      DO 300 I=2,IMM1
 300  TPS(I,J)=TPS(I,J)+U(I,J,K)*DZ(K)
      DO 302 K=1,KBM1
      DO 302 J=2,JMM1
      DO 302 I=2,IMM1
 302  U(I,J,K)=(U(I,J,K)-TPS(I,J))+
     1     (UTB(I,J)+UTF(I,J))/(DT(I,J)+DT(I-1,J))
      DO 303 J=1,JM
      DO 303 I=1,IM
 303  TPS(I,J)=0.
      DO 304 K=1,KBM1
      DO 304 J=2,JMM1
      DO 304 I=2,IMM1
 304  TPS(I,J)=TPS(I,J)+V(I,J,K)*DZ(K)
      DO 306 K=1,KBM1
      DO 306 J=2,JMM1
      DO 306 I=2,IMM1
 306  V(I,J,K)=(V(I,J,K)-TPS(I,J))+
     2     (VTB(I,J)+VTF(I,J))/(DT(I,J)+DT(I,J-1))
C
C----------------------------------------------------------------
C     VERTVL INPUT = U,V,DT(=H+ET),ETF,ETB; OUTPUT = W
C----------------------------------------------------------------
      CALL VERTVL(DTI2)
      CALL BCOND(5)
C
C
      DO 307 K=1,KB
      DO 307 J=1,JM
      DO 307 I=1,IM
      UF(I,J,K)=0.E0
  307 VF(I,J,K)=0.E0
C----------------------------------------------------------------
C     COMPUTE Q2F AN Q2LF USING UF AND VF AS TEMPORARY VARIABLES
C----------------------------------------------------------------
      CALL ADVQ(Q2B,Q2,DTI2,UF)
      CALL ADVQ(Q2LB,Q2L,DTI2,VF)
      CALL PROFQ(DTI2)
      CALL BCOND(6)
      DO 310 K=1,KB
      DO 310 J=1,JM
      DO 310 I=1,IM
      Q2(I,J,K)=Q2(I,J,K)
     1    +.5*SMOTH*(UF(I,J,K)+Q2B(I,J,K)-2.*Q2(I,J,K))
      Q2L(I,J,K)=Q2L(I,J,K)
     1    +.5*SMOTH*(VF(I,J,K)+Q2LB(I,J,K)-2.*Q2L(I,J,K))
      Q2B(I,J,K)=Q2(I,J,K)
      Q2(I,J,K)=UF(I,J,K)
      Q2LB(I,J,K)=Q2L(I,J,K)
      Q2L(I,J,K)=VF(I,J,K)
  310 CONTINUE
C---Re-calculate drag coeff according to Adams & Weatherly (1981)-
	do j=1,jm
	do i=1,im
	if (rff(i,j,kbm1).le.0.0) rff(i,j,kbm1)=0.0E0
	if (rff(i,j,kbm1).ge.0.21) rff(i,j,kbm1)=0.21E0
      CBC(I,J)=MAX(CBCMIN,.16E0/(1.E0+5.5*RFF(I,J,KBM1))**2
     1        /LOG((ZZ(KBM1)-Z(KB))*H(I,J)
     1        /Z0B)**2)*FSM(I,J)
	end do
	end do
C----------------------------------------------------------------
C     COMPUTE TF AN SF USING UF AND VF AS TEMPORARY VARIABLES
C----------------------------------------------------------------
      IF(MODE.EQ.4) GO TO 360
C  In this calculation, S = constant. However, we calculate
C  S to demonstrate that the same subroutines can be used to 
C  calculate multiple scalars.
      CALL ADVT(TB,T,TMEAN,DTI2,UF)
      CALL ADVT(SB,S,SMEAN,DTI2,VF)
	if (time.ge.5.0) then
      CALL ADVTRSML(RB1,RR1,RRMEAN,DTI2,RF1,3,sink1)
	end if
      CALL PROFT(UF,WTSURF,SWRAD,TSURF,NBC,DTI2)
      CALL PROFT(VF,WSSURF,SWRAD,SSURF,  1,DTI2)
	if (time.ge.5.0) then
      CALL PROFTR
     1 (RF1,WRSURF,DTI2,TTT1,srho1,sink1,svel1)
	end if
      CALL BCOND(4)
C
c==================================================================
c     SEDIMENT CONCENTRATION OPEN BOUNDARY CONDITIONS
c==================================================================

      CALL BCOND(7)
C
      DO 355 K=1,KB
      DO 355 J=1,JM
      DO 355 I=1,IM
c==================================================================
c     CHECK FOR NEGATIVE SALINITIES (IF NEEDED.....)
c==================================================================
      if((rf1(i,j,k)+RBIAS).lt.0.) rf1(i,j,k)=-RBIAS
      if((rf2(i,j,k)+RBIAS).lt.0.) rf2(i,j,k)=-RBIAS
c
c --------------------------------------------------------------
c   TIME SMOOTHING AND STEP FORWARD IN TIME FOR SED, T, S
c -------------------------------------------------------------

c   ------- sediment conc. --------------------

      RR1(I,J,K)=RR1(I,J,K)+.5*SMOTH*(RF1(I,J,K)
     1  +RB1(I,J,K)-2.E0*RR1(I,J,K))
      RR2(I,J,K)=RR2(I,J,K)+.5*SMOTH*(RF2(I,J,K)
     1  +RB2(I,J,K)-2.E0*RR2(I,J,K))

      RB1(I,J,K)=RR1(I,J,K)
      RB2(I,J,K)=RR2(I,J,K)

      RR1(I,J,K)=RF1(I,J,K)
      RR2(I,J,K)=RF2(I,J,K)

      T(I,J,K)=T(I,J,K)+.5*SMOTH*(UF(I,J,K)
     1  +TB(I,J,K)-2.E0*T(I,J,K))
      S(I,J,K)=S(I,J,K)+.5*SMOTH*(VF(I,J,K)
     1  +SB(I,J,K)-2.E0*S(I,J,K))
      TB(I,J,K)=T(I,J,K)
      T(I,J,K)=UF(I,J,K)
      SB(I,J,K)=S(I,J,K)
  355 S(I,J,K)=VF(I,J,K)

  360 CONTINUE
C
      CALL DENS(S,T,RR1,RR2,RHO)
C
C----------------------------------------------------------------
C     COMPUTE UF AND VF
C----------------------------------------------------------------
      CALL ADVU(DRHOX,ADVX,DTI2)
      CALL ADVV(DRHOY,ADVY,DTI2)
      CALL PROFU(DTI2)
      CALL PROFV(DTI2)
      CALL BCOND(3)
C
      DO 369 J=1,JM
      DO 369 I=1,IM
  369 TPS(I,J)=0.E0
      DO 370 K=1,KBM1
      DO 370 J=1,JM
      DO 370 I=1,IM
  370 TPS(I,J)=TPS(I,J)+(UF(I,J,K)+UB(I,J,K)-2.E0*U(I,J,K))*DZ(K)
      DO 372 K=1,KBM1
      DO 372 J=1,JM
      DO 372 I=1,IM
  372 U(I,J,K)=U(I,J,K)+.5*SMOTH*(UF(I,J,K)+UB(I,J,K)-2.E0*U(I,J,K)
     1        -TPS(I,J))
      DO 3721 J=1,JM
      DO 3721 I=1,IM
 3721 TPS(I,J)=0.E0
      DO 374 K=1,KBM1
      DO 374 J=1,JM
      DO 374 I=1,IM
  374 TPS(I,J)=TPS(I,J)+(VF(I,J,K)+VB(I,J,K)-2.E0*V(I,J,K))*DZ(K)
      DO 376 K=1,KBM1
      DO 376 J=1,JM
      DO 376 I=1,IM
  376 V(I,J,K)=V(I,J,K)+.5*SMOTH*(VF(I,J,K)+VB(I,J,K)-2.E0*V(I,J,K)
     1        -TPS(I,J))
      DO 377 K=1,KB
      DO 377 J=1,JM
      DO 377 I=1,IM
      UB(I,J,K)=U(I,J,K)
      U(I,J,K)=UF(I,J,K)
      VB(I,J,K)=V(I,J,K)
  377 V(I,J,K)=VF(I,J,K)
C
C
 8200 CONTINUE
C
      DO 8210 J=1,JM
      DO 8210 I=1,IM
      EGB(I,J)=EGF(I,J)
      ETB(I,J)=ET(I,J)
      ET(I,J)=ETF(I,J)
      DT(I,J)=H(I,J)+ET(I,J)
      UTB(I,J)=UTF(I,J)
 8210 VTB(I,J)=VTF(I,J)
C
C
C--------------------------------------------------------------------
C           BEGIN PRINT SECTION
C--------------------------------------------------------------------
      IF(MOD(IINT,IPRINT3).EQ.0) THEN
	j0=5
	t9=24.0*time
      DO 888 K=1,KBM1
      DO 888 I=1,IM
888      WRITE(85) t9,UB(I,j0,K),VB(I,j0,K),UAB(I,j0),VAB(I,j0)
     1  ,ELB(I,j0),RB2(I,j0,K),srho1(i,j0),RB1(i,j0,k),ttt1(i,j0,k)
     3  ,kh(i,j0,k),km(i,j0,k),rff(i,j0,k),cbc(i,j0),ustar2(i,j0,k)
     4  ,sink1(i,j0,k)
	ENGTOT=0.E0
      VTOT=0.E0
      TTOT=0.E0
      STOT=0.E0
      DO 8888 J=1,JM
      DO 8888 I=1,IM
	ENGTOT=ENGTOT+(h(i,j)+elb(i,j))*0.5*(UAB(I,J)**2+VAB(I,J)**2)
      DO 8888 K=1,KBM1
      DVTOT=DX(I,J)*DY(I,J)*DT(I,J)*DZ(K)*FSM(I,J)
      VTOT=VTOT+DVTOT
      TTOT=TTOT+TB(I,J,K)*DVTOT
      STOT=STOT+SB(I,J,K)*DVTOT
 8888 CONTINUE
      TTOT=TTOT/VTOT
      STOT=STOT/VTOT
	WRITE(81,*) TIME,ENGTOT
      END IF

        IF(IINT.GE.ISWTCH) IPRINT=iprint3
      IF(MOD(IINT,IPRINT).NE.0) GO TO 7000 
      WRITE(80) t9,UB,VB,UAB,VAB,ELB,SB,KH,
     *   RB1,ttt1,srho1,KM,RB2,
     *   RFF,CBC,USTAR2
9001  CONTINUE
        WRITE(6,'(/,'' TIME ='',F10.2,''    IINT ='',I8,
     1  ''   IEXT ='',I8,''   IPRINT ='',I5,//)') TIME,IINT,IEXT,IPRINT
c        CALL PRXY(' TRNU      ',TIME,TRNU,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' TRNV      ',TIME,TRNV,IP,ISK,JM,1,0.E0)
c       CALL PRXY(' ADVUA     ',TIME,ADVUA,IP,ISK,JM,1,0.0)
c       CALL PRXY(' ADVVA     ',TIME,ADVVA,IP,ISK,JM,1,0.0)
c       CALL PRXY(' ADVUU     ',TIME,ADVUU,IP,ISK,JM,1,0.0)
c       CALL PRXY(' ADVVV     ',TIME,ADVVV,IP,ISK,JM,1,0.0)
c        CALL PRXY(' WTSURF    ',TIME,WTSURF,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' WUSURF    ',TIME,WUSURF,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' WVSURF    ',TIME,WVSURF,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' WUBOT     ',TIME,WUBOT,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' WVBOT     ',TIME,WVBOT,IP,ISK,JM,1,0.E0)
        CALL PRXY(' CBC            ',TIME,CBC,IP,ISK,JM,1,0.E0)
c        CALL PRXY('  FREE SURFACE ',TIME, ELB,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' AVERAGE U COMP.',TIME,UAB,IP,ISK,JM,1,0.E0)
c        CALL PRXY(' AVERAGE V COMP.',TIME,VAB,IP,ISK,JM,1,0.E0)
c        CALL FINDPSI
c        IF(MODE.NE.2) THEN
          CALL PRXYZ('AAM       ',TIME,AAM,IP,ISK,JM,1,KB,0.E0)
c          CALL PRXYZ('UB ',TIME,UB,IP,ISK,JM,1,KB,0.E0)
c          CALL PRXYZ('VB ',TIME,VB,IP,ISK,JM,1,KB,0.E0)
         CALL PRXYZ(' RB1    ',TIME,RB1,IP,ISK,JM,1,KB,0.E0)
c         CALL PRXYZ(' RB2    ',TIME,RB2,IP,ISK,JM,1,KB,0.E0)
         CALL PRXYZ(' SB    ',TIME,SB,IP,ISK,JM,1,KB,0.E0)
        CALL PRXY('WSSURF    ',TIME,WSSURF,IP,ISK,JM,1,0.E0)
c         CALL PRXYZ(' TB    ',TIME,TB,IP,ISK,JM,1,KB,0.E0)
c          CALL PRXZ('TB' ,TIME,TB,IM,ISK,JM,35,45,KB,.01,DT,ZZ)
c          CALL PRXZ('SB' ,TIME,SB,IM,ISK,JM,5,6,KB,1.E0,DT,ZZ)
c          CALL PRXZ('RHO',TIME,RHO,IM,ISK,JM,25,45,KB,0.0001E0,DT,ZZ)
c          CALL PRXZ('Q2' ,TIME,Q2,IM,ISK,JM,35,45,KB,0.0,DT,Z)
c          CALL PRXZ('KH' ,TIME,KH,IM,ISK,JM,35,47,KB,0.0,DT,Z)
c          CALL PRXZ('KM' ,TIME,KM,IM,ISK,JM,35,47,KB,0.0,DT,Z)
c          CALL PRXZ('UB' ,TIME,UB,IM,ISK,JM,35,47,KB,0.0,DT,ZZ)
c          CALL PRXZ('VB' ,TIME,VB,IM,ISK,JM,35,47,KB,0.0,DT,ZZ)
c        ENDIF

C
C
c      WRITE(6,'(''  VTOT ='',E20.8,''    TTOT ='',E20.8,
c     1      ''   STOT ='',E20.8)') VTOT,TTOT,STOT

      VTOT=0.E0
      TTOT=0.E0
      STOT1=0.E0
      STOT2=0.E0
      RTOT1=0.E0
      RTOT2=0.E0
      ATOT=0.E0
      DO  J = 1, JM
      DO  I = 1, IM
        RTOT1=RTOT1+srho1(i,j)*ART(I,J)*FSM(I,J)
        RTOT2=RTOT2+srho2(i,j)*ART(I,J)*FSM(I,J)
c      ATOT=ATOT+ART(I,J)*FSM(I,J)
        end do
        end do
      DO  K=1,KBM1
      DO  J=1,JM
      DO  I=1,IM
      DVTOT=DX(I,J)*DY(I,J)*H(I,J)*DZ(K)*FSM(I,J)
      VTOT=VTOT+DVTOT
      STOT1=STOT1+(RB1(I,J,K)+RBIAS)*DVTOT
      STOT2=STOT2+(RB2(I,J,K)+RBIAS)*DVTOT
c      RTOT=RTOT+RB1(I,J,K)*DVTOT
        end do
        end do
        end do
c==========================================================================
c     CALCULATION OF THE BASIN AVERAGED TEMPERATURE AND SALINITY
c==========================================================================
      TTOT1=STOT1/VTOT
      TTOT2=STOT2/VTOT
c      RTOT=RTOT/ATOT
      ETOT1=STOT1+RTOT1
      ETOT2=STOT2+RTOT2
      WRITE(6,'('' ETOT1 =,'',E20.8,'' TTOT1 ='',E20.8,
     1 '' STOT1 ='',E20.8,'' RTOT1='',E20.8)') ETOT1,TTOT1,STOT1,RTOT1
      WRITE(6,'('' ETOT2 =,'',E20.8,'' TTOT2 ='',E20.8,
     1 '' STOT2 ='',E20.8,'' RTOT2='',E20.8)') ETOT2,TTOT2,STOT2,RTOT2

      IF(VAMAX.GT.100.) THEN
        CALL PRXY('  EL ',TIME, EL ,IP,ISK,JM,1,0.E0)
        CALL PRXY('  ELF',TIME, ELF,IP,ISK,JM,1,0.E0)
        CALL PRXY('  UA            ',TIME,UA ,IP,ISK,JM,1,0.E0)
        CALL PRXY('  VA            ',TIME,VA ,IP,ISK,JM,1,0.E0)
        CALL PRXY('  UAF           ',TIME,UAF,IP,ISK,JM,1,0.E0)
        CALL PRXY('  VAF           ',TIME,VAF,IP,ISK,JM,1,0.E0)
      WRITE(6,'(//////////////////)')
      WRITE(6,'(''************************************************'')')
      WRITE(6,'(''************ ABNORMAL JOB END ******************'')')
      WRITE(6,'(''************* USER TERMINATED ******************'')')
      WRITE(6,'(''************************************************'')')
      WRITE(6,'('' VAMAX ='',E12.3)') VAMAX
      WRITE(6,'('' IINT,IEXT ='',2I6)') IINT,IEXT
      STOP
      ENDIF
 7000 CONTINUE
C--------------------------------------------------------------------
C             END PRINT SECTION
C--------------------------------------------------------------------
 9000                     CONTINUE
C***********************************************************************
C                                                                      *
C                  END NUMERICAL INTEGRATION                           *
C                                                                      *
C***********************************************************************
c        CALL EPRXYZ('W VEL. COMP.',TIME,W,IM,1,JM,1,KB,.00001E0)
c        CALL EPRXYZ('Q2         ',TIME,Q2,IM,1,JM,1,KB,0.0)
C
      WRITE(90) TIME,
     1  WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB,
     2  UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADVUU,ADVVV,ADVUA,ADVVA,
     3  KM,KH,KQ,L,Q2,Q2B,AAM,Q2L,Q2LB

      STOP
      END
C
      SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE)
       INCLUDE 'comblk.h'
      DIMENSION ADVUA(IM,JM),ADVVA(IM,JM),CURV2D(IM,JM)
      EQUIVALENCE (TPS,CURV2D)
C---------------------------------------------------------------------
C      CALCULATE U-ADVECTION & DIFFUSION
C---------------------------------------------------------------------
C
C-------- ADVECTIVE FLUXES -------------------------------------------
C
      DO 200 J=1,JM
      DO 200 I=1,IM
 200  ADVUA(I,J)=0.
      DO 300 J=2,JM
      DO 300 I=2,IMM1
 300  FLUXUA(I,J)=.125E0*((D(I+1,J)+D(I,J))*UA(I+1,J)
     1                 +(D(I,J)+D(I-1,J))*UA(I,J))
     2                  *(UA(I+1,J)+UA(I,J))
      DO 400 J=2,JM
      DO 400 I=2,IM
 400  FLUXVA(I,J)=.125E0*((D(I,J)+D(I,J-1))*VA(I,J)
     1                 +(D(I-1,J)+D(I-1,J-1))*VA(I-1,J))
     2                    *(UA(I,J)+UA(I,J-1))
C----------- ADD VISCOUS FLUXES ---------------------------------
      DO  460 J=2,JM
      DO  460 I=2,IMM1
 460  FLUXUA(I,J)=FLUXUA(I,J)
     1         -D(I,J)*2.E0*AAM2D(I,J)*(UAB(I+1,J)-UAB(I,J))/DX(I,J)
      DO  470 J=2,JM
      DO  470 I=2,IM
      TPS(I,J)=.25E0*(D(I,J)+D(I-1,J)+D(I,J-1)+D(I-1,J-1))
     1            *(AAM2D(I,J)+AAM2D(I,J-1)+AAM2D(I-1,J)+AAM2D(I-1,J-1))
     2                *((UAB(I,J)-UAB(I,J-1))
     3                /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1))
     4                 +(VAB(I,J)-VAB(I-1,J))
     5                /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)) )
      FLUXUA(I,J)=FLUXUA(I,J)*DY(I,J)
      FLUXVA(I,J)=(FLUXVA(I,J)-TPS(I,J))
     1            *.25E0*(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1))
 470  CONTINUE
C----------------------------------------------------------------
      DO  480 J=2,JMM1
      DO  480 I=2,IMM1
 480  ADVUA(I,J)=FLUXUA(I,J)-FLUXUA(I-1,J)
     1           +FLUXVA(I,J+1)-FLUXVA(I,J)
C----------------------------------------------------------------
C       CALCULATE V-ADVECTION & DIFFUSION
C----------------------------------------------------------------
C
      DO 481 J=1,JM
      DO 481 I=1,IM
 481  ADVVA(I,J)=0.
C
C---------ADVECTIVE FLUXES ----------------------------
      DO 700 J=2,JM
      DO 700 I=2,IM
 700  FLUXUA(I,J)=.125E0*((D(I,J)+D(I-1,J))*UA(I,J)
     1         +(D(I,J-1)+D(I-1,J-1))*UA(I,J-1))*
     2                        (VA(I-1,J)+VA(I,J))
      DO 800 J=2,JMM1
      DO 800 I=2,IM
 800  FLUXVA(I,J)=.125E0*((D(I,J+1)+D(I,J))
     1       *VA(I,J+1)+(D(I,J)+D(I,J-1))*VA(I,J))
     2      *(VA(I,J+1)+VA(I,J))
C------- ADD VISCOUS FLUXES -----------------------------------
      DO  860 J=2,JMM1
      DO  860 I=2,IM
 860  FLUXVA(I,J)=FLUXVA(I,J)
     1        -D(I,J)*2.E0*AAM2D(I,J)*(VAB(I,J+1)-VAB(I,J))/DY(I,J)
      DO  870 J=2,JM
      DO  870 I=2,IM
      FLUXVA(I,J)=FLUXVA(I,J)*DX(I,J)
 870  FLUXUA(I,J)=(FLUXUA(I,J)-TPS(I,J))
     3             *.25E0*(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1))
C---------------------------------------------------------------
      DO  880 J=2,JMM1
      DO  880 I=2,IMM1
 880  ADVVA(I,J)=FLUXUA(I+1,J)-FLUXUA(I,J)
     1          +FLUXVA(I,J)-FLUXVA(I,J-1)
C
C---------------------------------------------------------------
      IF(MODE.NE.2) GO TO 5000
      DO 100 J=2,JMM1
      DO 100 I=2,IMM1
      WUBOT(I,J)=-0.5E0*(CBC(I,J)+CBC(I-1,J))
     1     *SQRT(UAB(I,J)**2+(.25E0*(VAB(I,J)
     2     +VAB(I,J+1)+VAB(I-1,J)+VAB(I-1,J+1)))**2)*UAB(I,J)
  100 CONTINUE
      DO 102 J=2,JMM1
      DO 102 I=2,IMM1
      WVBOT(I,J)=-0.5E0*(CBC(I,J)+CBC(I,J-1))
     1     *SQRT((.25E0*(UAB(I,J)+UAB(I+1,J)
     2     +UAB(I,J-1)+UAB(I+1,J-1)))**2+VAB(I,J)**2)*VAB(I,J)
  102 CONTINUE
      DO 120 J=2,JMM1
      DO 120 I=2,IMM1
      CURV2D(I,J)=.25*((VA(I,J+1)+VA(I,J))*(DY(I+1,J)-DY(I-1,J))
     1               -(UA(I+1,J)+UA(I,J))*(DX(I,J+1)-DX(I,J-1)) )
     2                /(DX(I,J)*DY(I,J))
      ADVUA(I,J)=ADVUA(I,J)   
     1    -ARU(I,J)*.25*(  CURV2D(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J))
     2          +CURV2D(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) )
      ADVVA(I,J)=ADVVA(I,J)
     1    +ARV(I,J)*.25*(  CURV2D(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J))
     2          +CURV2D(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) )
  120 CONTINUE  
 5000 CONTINUE
C 
      RETURN
      END
C               
      SUBROUTINE ADVQ(QB,Q,DTI2,QF)
       INCLUDE 'comblk.h'   
      DIMENSION QB(IM,JM,KB),Q(IM,JM,KB),QF(IM,JM,KB)
      DIMENSION XFLUX(IM,JM,KB),YFLUX(IM,JM,KB)
      EQUIVALENCE (XFLUX,A),(YFLUX,C)
C
C******* HORIZONTAL ADVECTION ************************************
      DO 110 K=2,KBM1
      DO 110 J=2,JM
      DO 110 I=2,IM
      XFLUX(I,J,K)=.125E0*(Q(I,J,K)+Q(I-1,J,K))
     1     *(DT(I,J)+DT(I-1,J))*(U(I,J,K)+U(I,J,K-1))
 110  YFLUX(I,J,K)=.125E0*(Q(I,J,K)+Q(I,J-1,K))*(DT(I,J)+DT(I,J-1))
     1        *(V(I,J,K)+V(I,J,K-1))
C******* HORIZONTAL DIFFUSION ************************************
      DO 315 K=2,KBM1
      DO 315 J=2,JM
      DO 315 I=2,IM
      XFLUX(I,J,K)=XFLUX(I,J,K)
     1    -.25*(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J,K-1)+AAM(I-1,J,K-1))
     2     *(H(I,J)+H(I-1,J))*(QB(I,J,K)-QB(I-1,J,K))*DUM(I,J)
     3           /(DX(I,J)+DX(I-1,J))
      YFLUX(I,J,K)=YFLUX(I,J,K)
     1    -.25*(AAM(I,J,K)+AAM(I,J-1,K)+AAM(I,J,K-1)+AAM(I,J-1,K-1))
     2     *(H(I,J)+H(I,J-1))*(QB(I,J,K)-QB(I,J-1,K))*DVM(I,J)
     3           /(DY(I,J)+DY(I,J-1))
      XFLUX(I,J,K)=.5E0*(DY(I,J)+DY(I-1,J))*XFLUX(I,J,K)
      YFLUX(I,J,K)=.5E0*(DX(I,J)+DX(I,J-1))*YFLUX(I,J,K)
  315 CONTINUE
C****** VERTICAL ADVECTION; ADD FLUX TERMS ;THEN STEP FORWARD IN TIME **
      DO 230 K=2,KBM1
      DO 230 J=2,JMM1
      DO 230 I=2,IMM1
      QF(I,J,K)=(W(I,J,K-1)*Q(I,J,K-1)-W(I,J,K+1)*Q(I,J,K+1))
     1                     /(DZ(K)+DZ(K-1))*ART(I,J)
     2                      +XFLUX(I+1,J,K)-XFLUX(I,J,K)
     3                      +YFLUX(I,J+1,K)-YFLUX(I,J,K)
 230  QF(I,J,K)=((H(I,J)+ETB(I,J))*ART(I,J)*QB(I,J,K)-DTI2*QF(I,J,K))
     1             /((H(I,J)+ETF(I,J))*ART(I,J))
C
      RETURN
      END
C
      SUBROUTINE ADVCT(ADVX,ADVY)        
C======================================================================
C  This routine calculates the horizontal portions of momentum advection
C  well in advance of their use in ADVU and ADVV so that their vertical 
C  integrals (created in MAIN) may be used in the external mode calculation.
C======================================================================
       INCLUDE 'comblk.h'   
      DIMENSION ADVX(IM,JM,KB),ADVY(IM,JM,KB),
     1          XFLUX(IM,JM,KB),YFLUX(IM,JM,KB),
     1          CURV(IM,JM,KB)
      EQUIVALENCE (A,XFLUX),(C,YFLUX),(EE,CURV)
C
C
      DO 60 K=1,KBM1
      DO 60 J=2,JMM1
      DO 60 I=2,IMM1
  60  CURV(I,J,K)=
     1      +.25E0*((V(I,J+1,K)+V(I,J,K))*(DY(I+1,J)-DY(I-1,J))
     2               -(U(I+1,J,K)+U(I,J,K))*(DX(I,J+1)-DX(I,J-1)) )
     3                /(DX(I,J)*DY(I,J))
C-----------------------------------------------------------------
C      Calculate x-component of velocity advection                 
C-----------------------------------------------------------------
      DO 59 K=1,KB
      DO 59 J=1,JM
      DO 59 I=1,IM
      ADVX(I,J,K)=0.E0
      XFLUX(I,J,K)=0.E0
  59  YFLUX(I,J,K)=0.E0
C
C******** HORIZONTAL ADVECTION FLUXES *****************************
      DO 100 K=1,KBM1
      DO 100 J=1,JM
      DO 100 I=2,IMM1
 100  XFLUX(I,J,K)=.125E0*((DT(I+1,J)+DT(I,J))*
     1              U(I+1,J,K)+(DT(I,J)+DT(I-1,J))
     2             *U(I,J,K))*(U(I+1,J,K)+U(I,J,K))
      DO 120 K=1,KBM1
      DO 120 J=2,JM
      DO 120 I=2,IM
 120  YFLUX(I,J,K)=.125E0*((DT(I,J)+DT(I,J-1))
     1        *V(I,J,K)+(DT(I-1,J)+DT(I-1,J-1))
     2        *V(I-1,J,K))*(U(I,J,K)+U(I,J-1,K))
C****** ADD HORIZONTAL DIFFUSION FLUXES ****************************
      DO  130 K=1,KBM1
      DO  130 J=2,JM
      DO  130 I=2,IMM1
      XFLUX(I,J,K)=XFLUX(I,J,K)
     1        -DT(I,J)*AAM(I,J,K)*2.E0*(UB(I+1,J,K)-UB(I,J,K))/DX(I,J)
      DTAAM=.25E0*(DT(I,J)+DT(I-1,J)+DT(I,J-1)+DT(I-1,J-1))
     1           *(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J-1,K)+AAM(I-1,J-1,K))
      YFLUX(I,J,K)=YFLUX(I,J,K)
     1          -DTAAM*((UB(I,J,K)-UB(I,J-1,K))
     2                    /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1))
     3                    +(VB(I,J,K)-VB(I-1,J,K))
     4                    /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)))
C
      XFLUX(I,J,K)=DY(I,J)*XFLUX(I,J,K)
      YFLUX(I,J,K)=
     1     .25E0*(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1))*YFLUX(I,J,K)
  130 CONTINUE 
C
C
C******** DO HORIZ. ADVECTION *******
      DO 146 K=1,KBM1
      DO 146 J=2,JMM1
      DO 146 I=2,IMM1
 146  ADVX(I,J,K)= 
     1            +XFLUX(I,J,K)-XFLUX(I-1,J,K)
     2            +YFLUX(I,J+1,K)-YFLUX(I,J,K)
     3   -ARU(I,J)*.25*(CURV(I,J,K)*DT(I,J)*(V(I,J+1,K)+V(I,J,K))
     4             +CURV(I-1,J,K)*DT(I-1,J)*(V(I-1,J+1,K)+V(I-1,J,K)))
C
C-----------------------------------------------------------------
C      Calculate y-component of velocity advection                 
C-----------------------------------------------------------------
      DO 299 K=1,KB
      DO 299 J=1,JM
      DO 299 I=1,IM
      ADVY(I,J,K)=0.E0
      XFLUX(I,J,K)=0.E0
 299  YFLUX(I,J,K)=0.E0
C
C********** HORIZONTAL ADVECTION FLUXES **************************
      DO 300 K=1,KBM1
      DO 300 J=2,JM
      DO 300 I=2,IM
 300  XFLUX(I,J,K)=.125E0*((DT(I,J)+DT(I-1,J))*U(I,J,K)
     1                  +(DT(I,J-1)+DT(I-1,J-1))*U(I,J-1,K))
     2                        *(V(I,J,K)+V(I-1,J,K))
      DO 320 K=1,KBM1
      DO 320 J=2,JMM1
      DO 320 I=1,IM
 320  YFLUX(I,J,K)=.125E0*((DT(I,J+1)+DT(I,J))*V(I,J+1,K)
     1                  +(DT(I,J)+DT(I,J-1))*V(I,J,K))
     2                        *(V(I,J+1,K)+V(I,J,K))
C******* ADD HORIZONTAL DIFFUSION FLUXES **************************
      DO  700 K=1,KBM1
      DO  700 J=2,JMM1
      DO  700 I=2,IM
      DTAAM=.25E0*(DT(I,J)+DT(I-1,J)+DT(I,J-1)+DT(I-1,J-1))
     1           *(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J-1,K)+AAM(I-1,J-1,K))
      XFLUX(I,J,K)=XFLUX(I,J,K)
     1     -DTAAM*((UB(I,J,K)-UB(I,J-1,K))
     2         /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1))
     3               +(VB(I,J,K)-VB(I-1,J,K))
     4         /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)))
      YFLUX(I,J,K)=YFLUX(I,J,K)
     1        -DT(I,J)*AAM(I,J,K)*2.E0*(VB(I,J+1,K)-VB(I,J,K))/DY(I,J)
C
      XFLUX(I,J,K)
     1    =.25E0*(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1))*XFLUX(I,J,K)
      YFLUX(I,J,K)=DX(I,J)*YFLUX(I,J,K)
  700 CONTINUE
C
C********** DO HORIZ. ADVECTION ************
      DO 400 K=1,KBM1
      DO 400 J=2,JMM1
      DO 400 I=2,IMM1
 400  ADVY(I,J,K)=
     1            +XFLUX(I+1,J,K)-XFLUX(I,J,K)
     2            +YFLUX(I,J,K)-YFLUX(I,J-1,K)
     3   +ARV(I,J)*.25*(CURV(I,J,K)*DT(I,J)*(U(I+1,J,K)+U(I,J,K))
     4             +CURV(I,J-1,K)*DT(I,J-1)*(U(I+1,J-1,K)+U(I,J-1,K)))
C
      RETURN
      END
C               
      SUBROUTINE ADVT(FB,F,FMEAN,DTI2,FF)
C
C     THIS SUBROUTINE INTEGRATES CONSERVATIVE SCALAR EQUATIONS
C
       INCLUDE 'comblk.h'
      DIMENSION FB(IM,JM,KB),F(IM,JM,KB),FF(IM,JM,KB),FMEAN(IM,JM,KB)
      DIMENSION XFLUX(IM,JM,KB),YFLUX(IM,JM,KB)
      EQUIVALENCE (XFLUX,A),(YFLUX,C)
C
      DO 529 J=1,JM
      DO 529 I=1,IM
      F(I,J,KB)=F(I,J,KBM1)
 529  FB(I,J,KB)=FB(I,J,KBM1)
C
C******* DO ADVECTION FLUXES **************************************
      DO 530 K=1,KBM1
      DO 530 J=2,JM
      DO 530 I=2,IM
      XFLUX(I,J,K)=.25E0*((DT(I,J)+DT(I-1,J))
     2            *(F(I,J,K)+F(I-1,J,K))*U(I,J,K))
 530  YFLUX(I,J,K)=.25E0*((DT(I,J)+DT(I,J-1))
     2            *(F(I,J,K)+F(I,J-1,K))*V(I,J,K))
c      XFLUX(I,J,K)=0.0
c 530  YFLUX(I,J,K)=0.0
C******  ADD DIFFUSIVE FLUXES *************************************
C
      DO 99 K=1,KB
      DO 99 J=1,JM
      DO 99 I=1,IM
  99  FB(I,J,K)=FB(I,J,K)-FMEAN(I,J,K)
C
      DO 100 K=1,KBM1
      DO 100 J=2,JM
      DO 100 I=2,IM
      XFLUX(I,J,K)=XFLUX(I,J,K)
     1    -.5E0*(AAM(I,J,K)+AAM(I-1,J,K))*(H(I,J)+H(I-1,J))
     2    *(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J)/(TPRNU*(DX(I,J)+DX(I-1,J)))
      YFLUX(I,J,K)=YFLUX(I,J,K)
     1    -.5E0*(AAM(I,J,K)+AAM(I,J-1,K))*(H(I,J)+H(I,J-1))
     2    *(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J)/(TPRNU*(DY(I,J)+DY(I,J-1)))
      XFLUX(I,J,K)=.5E0*(DY(I,J)+DY(I-1,J))*XFLUX(I,J,K)
      YFLUX(I,J,K)=.5E0*(DX(I,J)+DX(I,J-1))*YFLUX(I,J,K)
  100 CONTINUE
C
      DO 101 K=1,KB
      DO 101 J=1,JM
      DO 101 I=1,IM
 101  FB(I,J,K)=FB(I,J,K)+FMEAN(I,J,K)
C
C****** DO VERTICAL ADVECTION *************************************
      DO 505 J=2,JMM1
      DO 505 I=2,IMM1
 505  FF(I,J,1)=-.5E0*DZR(1)*(F(I,J,1)+F(I,J,2))*W(I,J,2)*ART(I,J)
      DO 520 K=2,KBM1
      DO 520 J=2,JMM1
      DO 520 I=1,IM
 520  FF(I,J,K)=.5E0*DZR(K)*((F(I,J,K-1)+F(I,J,K))*W(I,J,K)
     1                  -(F(I,J,K)+F(I,J,K+1))*W(I,J,K+1))*ART(I,J)
C****** ADD NET HORIZONTAL FLUXES; THEN STEP FORWARD IN TIME **********
      DO 120 K=1,KBM1
      DO 120 J=2,JMM1
      DO 120 I=2,IMM1
      FF(I,J,K)=FF(I,J,K)
     1              +XFLUX(I+1,J,K)-XFLUX(I,J,K)
     2              +YFLUX(I,J+1,K)-YFLUX(I,J,K)
      FF(I,J,K)=(FB(I,J,K)*(H(I,J)+ETB(I,J))*ART(I,J)-DTI2*FF(I,J,K))
     1                    /((H(I,J)+ETF(I,J))*ART(I,J))
 120  CONTINUE
      RETURN
      END
C

      SUBROUTINE ADVU(DRHOX,ADVX,DTI2)
       INCLUDE 'comblk.h'   
      DIMENSION DRHOX(IM,JM,KB),ADVX(IM,JM,KB)                        
C
C  Do vertical advection
      DO 100 K=1,KB  
      DO 100 J=1,JM
      DO 100 I=1,IM
  100 UF(I,J,K)=0.
      DO 140 K=2,KBM1
      DO 140 J=1,JM
      DO 140 I=2,IM
 140  UF(I,J,K)=.25E0*(W(I,J,K)+W(I-1,J,K))*(U(I,J,K)+U(I,J,K-1))
C****COMBINE HOR. and VERT. ADVECTION with
C           -FVD + GDEG/DX + BAROCLINIC TERM **********************
      DO 150 K=1,KBM1
      DO 150 J=2,JMM1
      DO 150 I=2,IMM1
 150  UF(I,J,K)=ADVX(I,J,K)+DZR(K)*(UF(I,J,K)-UF(I,J,K+1))*ARU(I,J)
     1   -ARU(I,J)*.25*(COR(I,J)*DT(I,J)*(V(I,J+1,K)+V(I,J,K))
     2             +COR(I-1,J)*DT(I-1,J)*(V(I-1,J+1,K)+V(I-1,J,K)))
     3        +GRAV*.25E0*(DT(I,J)+DT(I-1,J))
     4        *.5*(EGF(I,J)-EGF(I-1,J)+EGB(I,J)-EGB(I-1,J))
     5        *(DY(I,J)+DY(I-1,J))
     6        +DRHOX(I,J,K)
C******* STEP FORWARD IN TIME ***********************************
      DO 190 K=1,KBM1
      DO 190 J=2,JMM1
      DO 190 I=2,IMM1
 190  UF(I,J,K)=
     1      ((H(I,J)+ETB(I,J)+H(I-1,J)+ETB(I-1,J))*ARU(I,J)*UB(I,J,K)
     2         -2.E0*DTI2*UF(I,J,K))
     3     /((H(I,J)+ETF(I,J)+H(I-1,J)+ETF(I-1,J))*ARU(I,J))
      RETURN
      END
C               
      SUBROUTINE ADVV(DRHOY,ADVY,DTI2)
       INCLUDE 'comblk.h'   
      DIMENSION DRHOY(IM,JM,KB),ADVY(IM,JM,KB)                    
C
C  Do vertical advection
      DO 100 K=1,KB  
      DO 100 J=1,JM
      DO 100 I=1,IM
  100 VF(I,J,K)=0.
      DO 140 K=2,KBM1
      DO 140 J=2,JM
      DO 140 I=1,IM
 140  VF(I,J,K)=.25*(W(I,J,K)+W(I,J-1,K))*(V(I,J,K)+V(I,J,K-1))
C****COMBINE HOR. and VERT. ADVECTION with
C           +FUD + GDEG/DY + BAROCLINIC TERM **********************
      DO 340 K=1,KBM1
      DO 340 J=2,JMM1
      DO 340 I=2,IMM1
 340  VF(I,J,K)=ADVY(I,J,K)+DZR(K)*(VF(I,J,K)-VF(I,J,K+1))*ARV(I,J)
     1   +ARV(I,J)*.25*(COR(I,J)*DT(I,J)*(U(I+1,J,K)+U(I,J,K))
     2             +COR(I,J-1)*DT(I,J-1)*(U(I+1,J-1,K)+U(I,J-1,K)))
     3       +GRAV*.25E0*(DT(I,J)+DT(I,J-1))
     4       *.5*(EGF(I,J)-EGF(I,J-1)+EGB(I,J)-EGB(I,J-1))
     5       *(DX(I,J)+DX(I,J-1))
     6       +DRHOY(I,J,K)
C******* STEP FORWARD IN TIME ***************************************
      DO 390 K=1,KBM1
      DO 390 J=2,JMM1
      DO 390 I=2,IMM1
  390 VF(I,J,K)=
     1     ((H(I,J)+ETB(I,J)+H(I,J-1)+ETB(I,J-1))*ARV(I,J)*VB(I,J,K)
     2          -2.E0*DTI2*VF(I,J,K))
     3    /((H(I,J)+ETF(I,J)+H(I,J-1)+ETF(I,J-1))*ARV(I,J))
      RETURN
      END
C               
      SUBROUTINE BAROPG(DRHOX,DRHOY,RMEAN)
       INCLUDE 'comblk.h'   
      DIMENSION DRHOX(IM,JM,KB),DRHOY(IM,JM,KB),
     & RMEAN(IM,JM,KB)
C
      DO 299 K=1,KB
      DO 299 J=1,JM
      DO 299 I=1,IM
 299  RHO(I,J,K)=RHO(I,J,K)-RMEAN(I,J,K)
C
C               X COMPONENT OF BAROCLINIC PRESSURE GRADIENT
C
      DO 300 J=2,JM
      DO 300 I=2,IM
  300 DRHOX(I,J,1)=.5E0*GRAV*(-ZZ(1))*(DT(I,J)+DT(I-1,J))
     2     *(RHO(I,J,1)-RHO(I-1,J,1))
      DO 310 K=2,KBM1
      DO 310 J=2,JM
      DO 310 I=2,IM
 310  DRHOX(I,J,K)=DRHOX(I,J,K-1)
     1      +GRAV*.25E0*(ZZ(K-1)-ZZ(K))*(DT(I,J)+DT(I-1,J))
     2      *(RHO(I,J,K)-RHO(I-1,J,K)+RHO(I,J,K-1)-RHO(I-1,J,K-1))
     3      +GRAV*.25E0*(ZZ(K-1)+ZZ(K))*(DT(I,J)-DT(I-1,J))
     4      *(RHO(I,J,K)+RHO(I-1,J,K)-RHO(I,J,K-1)-RHO(I-1,J,K-1))
C
      DO 360 K=1,KBM1
      DO 360 J=2,JM
      DO 360 I=2,IM
 360  DRHOX(I,J,K)=.25E0*(DT(I,J)+DT(I-1,J))*DRHOX(I,J,K)*DUM(I,J)
     1     *(DY(I,J)+DY(I-1,J))
C
C               Y COMPONENT OF BAROCLINIC PRESSURE GRADIENT
C
      DO 500 J=2,JM
      DO 500 I=2,IM
 500  DRHOY(I,J,1)=.5E0*GRAV*(-ZZ(1))*(DT(I,J)+DT(I,J-1))
     1     *(RHO(I,J,1)-RHO(I,J-1,1))
      DO 510 K=2,KBM1
      DO 510 J=2,JM
      DO 510 I=2,IM
 510  DRHOY(I,J,K)=DRHOY(I,J,K-1)
     1        +GRAV*.25E0*(ZZ(K-1)-ZZ(K))*(DT(I,J)+DT(I,J-1))
     2        *(RHO(I,J,K)-RHO(I,J-1,K)+RHO(I,J,K-1)-RHO(I,J-1,K-1))
     3        +GRAV*.25E0*(ZZ(K-1)+ZZ(K))*(DT(I,J)-DT(I,J-1))
     4        *(RHO(I,J,K)+RHO(I,J-1,K)-RHO(I,J,K-1)-RHO(I,J-1,K-1))
C
      DO 560 K=1,KBM1
      DO 560 J=2,JM
      DO 560 I=2,IM
 560  DRHOY(I,J,K)=.25E0*(DT(I,J)+DT(I,J-1))*DRHOY(I,J,K)*DVM(I,J)
     1       *(DX(I,J)+DX(I,J-1))
C
      DO 561 K=1,KB
      DO 561 J=2,JM
      DO 561 I=2,IM
      DRHOX(I,J,K)=RAMP*DRHOX(I,J,K)
 561  DRHOY(I,J,K)=RAMP*DRHOY(I,J,K)
C
      DO 571 K=1,KB
      DO 571 J=1,JM
      DO 571 I=1,IM
 571  RHO(I,J,K)=RHO(I,J,K)+RMEAN(I,J,K)
C
      RETURN
      END
C               
      SUBROUTINE BCOND(IDX)
       INCLUDE 'comblk.h'   
C
C  Closed boundary conditions are automatically enabled through
C  specification of the masks, DUM, DVM and FSM. For problems with
C  open boundaries, dependent variables must be specified on the
C  the boundary grid points. Sample boundary conditions, left over
C  from a problem with three open boundaries are to be seen below
C  but are commented out for the closed basin application.
C
C              SPECIFICATION OF OPEN BOUNDARY CONDITIONS
C        BRACKETS OF *  INDICATE INTERIOR (NON-BOUNDARY) GRID POINTS.
C        HORIZONTAL LOCATIONS OF EL, T AND S ARE COINCIDENT.
C
C        I-1           I         I+1
C
C
C      U(JM)       EL(JM)  =.   U(JM)           BC
C                           .
C                   V(JM)   .                   BC
C                           .
C     *U(JMM1)*   *EL(JMM1)*   *U(JMM1)*
C
C                 *V(JMM1)*
C                       .                         *V(IMM1)*           V(IM)
C                       .
C                       .              *U(IMM1)* *EL(IMM1)*   U(IM)   EL(IM)
C                       .                              ". . . . . . . ."
C                       .
C                       .                         *V(IMM1)*           V(IM)
C                       .
C                       .                                     BC        BC
C
C                    *V(3)*
C
C       *U(2)*       *EL(2)*   *U(2)*
C                           .
C                     V(2)  .             BC
C                           .
C         U(1)        EL(1)=.   U(1)      BC
C --------------------------------------------------------------------
C                     V(1) IS NOT USED
C
C
      DATA PI/3.14167E0/,GEE/9.807E0/,SMALL/1.0E-10/
C
      GO TO (10,20,30,40,50,60,70), IDX
C
 10   CONTINUE
C-----------------------------------------------------------------------
C                   EXTERNAL ELEV. B.C.'S
C-----------------------------------------------------------------------
C                   ***** SEAWARD ********
C
      DO 100 J=1,JM
C::External reading of tidal heights
        ELF(IMM1,J)=0.4*sin(2.*PI*TIME/(12.42/24.))
     1  +0.1*sin(2.*pi*TIME/(12./24.))
c        ELF(IMM1,J)=0.5*sin(2.*PI*TIME/(12.42/24.))
c     1  +0.1*sin(2.*pi*TIME/(12./24.))
ct        ELF(IMM1,J)=.5*cos(2.*PI*TIME/(12.42/24.)-PI*235./180.)
ct     &  +0.09*cos(2.*PI*TIME/(12./24.)-PI*272./180.)
ct     &  +0.22*cos(2.*PI*TIME/(23.93/24.)-PI*111./180.)
ct     &  +0.1*cos(2.*PI*TIME/(25.83/24.)-PI*67./180.)
       ELF(IM,J)=ELF(IMM1,J)
 100  CONTINUE

C                   **** NORTH AND SOUTH ****

       DO 110 I=1,IM
c       ELF(I,JMM1)=.5*sin(2.*PI*TIME/(12.42/24.))
c       ELF(I,2)=.5*sin(2.*PI*TIME/(12.42/24.))
c       ELF(I,1)=ELF(I,2)
c       ELF(I,JM)=ELF(I,JMM1)
C:: Cyclic Cond
c       ELF(I,1)=ELF(I,JMM2)
c       ELF(I,2)=ELF(I,JMM1)
c       ELF(I,JM)=ELF(I,3)
c       CPHNOTH = (ELB(I,JMM1)-ELF(I,JMM1))/
c     &          (ELF(I,JMM1)+ELB(I,JMM1)-2.0*EL(I,JMM2)+SMALL)
c       CPHNOTH= AMAX1(0.0,CPHNOTH)
c       CPHNOTH=AMIN1(1.0,CPHNOTH)
c       ELF(I,JM) =  (ELB(I,JM) - CPHNOTH*(ELB(I,JM)-2.0*EL(I,JMM1)))/
c     &              (1.0+CPHNOTH)
c       CPHSOTH = (ELB(I,2)-ELF(I,2))/
c     &          (ELF(I,2)+ELB(I,2)-2.0*EL(I,3)+SMALL)
c       CPHSOTH= AMAX1(0.0,CPHSOTH)
c       CPHSOTH=AMIN1(1.0,CPHSOTH)
c       ELF(I,1) =  (ELB(I,1) - CPHSOTH*(ELB(I,1)-2.0*EL(I,2)))/
c     &              (1.0+CPHSOTH)
110   CONTINUE

    
       DO 111 J=1,JM
       DO 111 I=1,IM
 111   ELF(I,J)=ELF(I,J)*FSM(I,J)
       RETURN
C
 20   CONTINUE
C-----------------------------------------------------------------------
C                   EXTERNAL VEL B.C.'S
C-----------------------------------------------------------------------
C                   ***** SEAWARD ********
C    
C::  FORCING WITH A M2 TIDE
      DO 120 J=1,JM
      UAF(IM,J)=UAF(IMM1,J)
c      UAF(IM,J)=.1*sin(2.*PI*TIME/(12.42/24.))
120   VAF(IM,J)=0.E0

C                   **** NORTH AND SOUTH ****
c      DO 130 I=1,IMM1
C:: Cyclic Cond
c      VAF(I,2)=VAF(I,JMM1)
c      VAF(I,1)=VAF(I,2)
c      VAF(I,JM)=VAF(I,3)
c      UAF(I,1)=UAF(I,JMM1)
C::GRADIENT COND
c      VAF(I,JM)=VAF(I,JMM1)
c      VAF(I,1)=VAF(I,2)
c      VMID=.5E0*(VA(I,JM)+VA(I-1,JM))
c      UAF(I,JM)=UA(I,JM)-DTE/(DY(I,JM)+DY(I,JMM1))
cc     1    *(   (VMID+ABS(VMID))*(UA(I,JM)-UA(I,JMM1))
c     2        +(VMID-ABS(VMID))*(0.E0-UA(I,JM))  )
c      VMID=.5E0*(VA(I,2)+VA(I-1,2))
c      UAF(I,1)=UA(I,1)-DTE/(DY(I,1)+DY(I,2))
c     1    *(   (VMID+ABS(VMID))*(UA(I,1)-0.E0)
c     2        +(VMID-ABS(VMID))*(UA(I,2)-UA(I,1)) )
c130      UAF(I,JM)=UAF(I,2)

      DO 131 J=1,JM
      DO 131 I=1,IM
      UAF(I,J)=UAF(I,J)*DUM(I,J)
 131  VAF(I,J)=VAF(I,J)*DVM(I,J)
      RETURN
 30   CONTINUE
C-----------------------------------------------------------------------
C                   INTERNAL VEL B.C.'S
C-----------------------------------------------------------------------
C                         **** SEAWARD *******
      DO 142 K=1,KBM1
      DO 141 J=2,JMM1
      GA=SQRT(H(IM,J)/120.E0)
      UF(IM,J,K)
     1  =GA*(.25E0*U(IMM1,J-1,K)+.5E0*U(IMM1,J,K)+.25E0*U(IMM1,J+1,K))
     2   +(1.E0-GA)*(.25E0*U(IM,J-1,K)+.5E0*U(IM,J,K)+.25E0*U(IM,J+1,K))
      UMID=.5E0*(U(IM,J,K)+U(IM,J-1,K))
      VF(IM,J,K)=V(IM,J,K)-DTI/(DX(IM,J)+DX(IMM1,J))
     1     *(  (UMID+ABS(UMID))*(V(IM,J,K)-V(IMM1,J,K))
     2        +(UMID-ABS(UMID))*(0.E0 -V(IM,J,K))  )
 141  CONTINUE
 142  CONTINUE

C                        **** NORTH AND SOUTH *****

c      DO 156 I=1,IMM1
c      DO 156 K=1,KBM1
C:: Cyclic Cond
c      VF(I,2,K)=VF(I,JMM1,K)
c      VF(I,1,K)=VF(I,2,K)
c      VF(I,JM,K)=VF(I,3,K)
c      UF(I,1,K)=UF(I,JMM1,K)
c      UF(I,JM,K)=UF(I,2,K)
C:: GRAD COND.
c      VF(I,JM,K) =  VF(I,JMM1,K)
c      UF(I,1,K) = UF(I,2,K)
c      VF(I,JM,K) =  VF(I,JMM1,K)
c      UF(I,1,K) = UF(I,2,K)
c 156  CONTINUE

      DO 160 K=1,KBM1
C    
      DO 160 J=1,JM
      DO 160 I=1,IM
      UF(I,J,K)=UF(I,J,K)*DUM(I,J)
      VF(I,J,K)=VF(I,J,K)*DVM(I,J)
 160  CONTINUE
C      
      RETURN
C
 40   CONTINUE
C-----------------------------------------------------------------------
C                   TEMP & SAL B.C.'S
C-----------------------------------------------------------------------
C                           *** SEAWARD *****
C    
      DO 220 K=1,KBM1
      DO 220 J=1,JM
      UF(IM,J,K)=T(IM,J,K)-DTI/(DX(IM,J)+DX(IMM1,J))
     1     *(  (U(IM,J,K)+ABS(U(IM,J,K)))*(T(IM,J,K)-T(IMM1,J,K))
     2        +(U(IM,J,K)-ABS(U(IM,J,K)))*(TBE(J,K) -T(IM,J,K))  )
 220  CONTINUE

C                        **** NORTH AND SOUTH *****
c      DO 230 K=1,KBM1
c      DO 230 I=1,IM
C-------GRAD COND-------------
c      UF(I,JM,K)=UF(I,JMM1,K)
c      UF(I,1,K)=UF(I,2,K)
c      VF(I,JM,K)=VF(I,JMM1,K)
c      VF(I,1,K)=VF(I,2,K)
C-------CYCLIC COND----------
c      UF(I,1,K)=UF(I,JMM2,K)
c      UF(I,2,K)=UF(I,JMM1,K)
c      UF(I,JM,K)=UF(I,3,K)
c      VF(I,1,K)=VF(I,JMM2,K)
c      VF(I,2,K)=VF(I,JMM1,K)
c      VF(I,JM,K)=VF(I,3,K)
C-------RADIATION COND---------
c      CPHSOTH=(TB(I,2,K)-T(I,2,K))/
c     1          (TB(I,3,K)-TB(I,2,K)+SMALL)
c       CPHSOTH=AMIN1(0.0,CPHSOTH)
c       CPHSOTH=AMAX1(-1.0,CPHSOTH)
c       UF(I,1,K)=T(I,1,K)-CPHSOTH*(T(I,2,K)-T(I,1,K))
c       CPHNOTH=(TB(I,JMM1,K)-T(I,JMM1,K))/
c     1          (TB(I,JMM1,K)-TB(I,JMM2,K)+SMALL)
c       CPHNOTH=AMIN1(1.0,CPHNOTH)
c       CPHNOTH=AMAX1(0.0,CPHNOTH)
c      UF(I,JM,K)=T(I,JM,K)-CPHNOTH*(T(I,JM,K)-T(I,JMM1,K))
230	CONTINUE

      DO 240 K=1,KBM1
      DO 240 J=1,JM
      DO 240 I=1,IM
      UF(I,J,K)=UF(I,J,K)*FSM(I,J)
      VF(I,J,K)=VF(I,J,K)*FSM(I,J)
 
 240  CONTINUE
      RETURN
C
C
 50   CONTINUE
C---------------VERTICAL VEL. B. C.'S --------------------------------
      DO 250 K=1,KBM1
      DO 250 J=1,JM
      DO 250 I=1,IM
      W(I,J,K)=W(I,J,K)*FSM(I,J)
 250  CONTINUE
      DO 251 K=1,KB
      DO 251 J=1,JM
 251  W(IM,J,K)=0.E0
C    
      RETURN
C
 60   CONTINUE
C---------------- Q2 AND Q2L B.C.'S -----------------------------------
      DO 300 K=1,KB
C     
      DO 295 J=1,JM
      UF(IM,J,K)=1.E-10
 295  VF(IM,J,K)=1.E-10
      DO 296 I=1,IM
      UF(I,JM,K)=1.E-10
      VF(I,JM,K)=1.E-10
      UF(I,1,K)=1.E-10
 296  VF(I,1,K)=1.E-10
      DO 297 J=1,JM
      DO 297 I=1,IM
      UF(I,J,K)=UF(I,J,K)*FSM(I,J)
 297  VF(I,J,K)=VF(I,J,K)*FSM(I,J)
 300  CONTINUE
      RETURN

70    CONTINUE
C-----------------------------------------------------------------------
C                   Sediment concentration BC
C-----------------------------------------------------------------------
C                           *** SEAWARD *****

C-------RADIATION COND---------
      DO 930 K=1,KBM1
      DO 930 J=1,JM
      RF1(IM,J,K)=RR1(IM,J,K)-DTI/(DX(IM,J)+DX(IMM1,J))
     1     *(  (U(IM,J,K)+ABS(U(IM,J,K)))*(RR1(IM,J,K)-RR1(IMM1,J,K))
     2        +(U(IM,J,K)-ABS(U(IM,J,K)))*(0.0E0 -RR1(IM,J,K))  )
      RF2(IM,J,K)=RR2(IM,J,K)-DTI/(DX(IM,J)+DX(IMM1,J))
     1     *(  (U(IM,J,K)+ABS(U(IM,J,K)))*(RR2(IM,J,K)-RR2(IMM1,J,K))
     2        +(U(IM,J,K)-ABS(U(IM,J,K)))*(0.0E0 -RR2(IM,J,K))  )
 930  CONTINUE
c
C                        **** NORTH AND SOUTH *****
c      DO 935 K=1,KBM1
c      DO 935 I=1,IM

C-------GRAD COND----------
c      RF1(I,1,K)=RF1(I,2,K)
c      RF2(I,1,K)=RF2(I,2,K)
c      RF2(I,JM,K)=RF2(I,JMM1,K)
c      RF2(I,JM,K)=RF2(I,JMM1,K)
C-------CYCLIC COND----------
c      RF1(I,1,K)=RF1(I,JMM2,K)
c      RF1(I,2,K)=RF1(I,JMM1,K)
c      RF1(I,JM,K)=RF1(I,3,K)
c      RF2(I,1,K)=RF2(I,JMM2,K)
c      RF2(I,2,K)=RF2(I,JMM1,K)
c      RF2(I,JM,K)=RF2(I,3,K)
c935     CONTINUE

      DO 940 K=1,KBM1
      DO 940 J=1,JM
      DO 940 I=1,IM

      RF1(I,J,K)=RF1(I,J,K)*FSM(I,J)
      RF2(I,J,K)=RF2(I,J,K)*FSM(I,J)

 940  CONTINUE
      RETURN

      END


      SUBROUTINE DENS(SI,TI,RR1I,RR2I,RHOO)
       INCLUDE 'comblk.h'
      DIMENSION SI(IM,JM,KB),TI(IM,JM,KB),RHOO(IM,JM,KB)
      DIMENSION RR1I(IM,JM,KB),RR2I(IM,JM,KB)
C       If using 32 bit precision, it is recommended that
C       TR, SR, P, RHOR , CR be made double precision.
C
C         THIS SUBROUTINE COMPUTES DENSITY- 1.025
C         T = POTENTIAL TEMPERATURE
C
      DO 1 K=1,KBM1
      DO 1 J=1,JM
      DO 1 I=1,IM
      TR=TI(I,J,K)+10.0E0
      SR=SI(I,J,K)+35.0E0
      CRR1=RR1I(I,J,K)+2.0E-3
      CRR2=RR2I(I,J,K)+2.0E-3
C            Approximate pressure in units of bars
      P=-GRAV*1.025*ZZ(K)*DT(I,J)*0.01
C
      RHOR = 999.842594 + 6.793952E-2*TR
     $        - 9.095290E-3*TR**2 + 1.001685E-4*TR**3
     $        - 1.120083E-6*TR**4 + 6.536332E-9*TR**5
C
      RHOR = RHOR + (0.824493 - 4.0899E-3*TR
     $       + 7.6438E-5*TR**2 - 8.2467E-7*TR**3
     $       + 5.3875E-9*TR**4) * SR
     $       + (-5.72466E-3 + 1.0227E-4*TR
     $       - 1.6546E-6*TR**2) * SR**1.5
     $       + 4.8314E-4 * SR**2
C
      CR=1449.1+.0821*P+4.55*TR-.045*TR**2
     $                              +1.34*(SR-35.)
      RHOR=RHOR + 1.E5*P/CR**2
     $     *(1.-2.*P/CR**2)
C
      RHOR=RHOR
     $     + (1.-RHOR/1100.0)*CRR1
C
      RHOO(I,J,K)=(RHOR-1025.)*1.E-3*FSM(I,J)
    1 CONTINUE
C
      DO 3 J=1,JM
      DO 3 I=1,IM
    3 RHOO(I,J,KB)=1.
      RETURN
      END
C               
               
      SUBROUTINE DEPTH(Z,ZZ,DZ,DZZ,DZR,KB,KBM1)
C >>>
      DIMENSION Z(KB),ZZ(KB),DZ(KB),DZZ(KB),DZR(KB)
      KL1=5
      KL2=KB-5
C ***
C    THIS SUBROUTINE ESTABLISHES THE VERTICAL RESOLUTION WITH LOG
C    DISTRIBUTIONS  AT THE TOP AND BOTTOM AND A LINEAR DISTRIBUTION
C    BETWEEN. CHOOSE KL1 AND KL2. THE DEFAULT KL1 = .3*KB AND KL2 = KB-2
C    YIELDS A LOG DISTRIBUTION AT THE TOP AND NONE AT THE BOTTOM.
C ***
      BB=FLOAT(KL2-KL1)+4.
      CC=FLOAT(KL1)-2.
      DEL1=2./BB/EXP(.693147*FLOAT(KL1-2))
      DEL2=2./BB/EXP(.693147*FLOAT(KB-KL2-1))
      Z(1)=0.
      ZZ(1)=-DEL1/2.
      DO 30 K=2,KL1-2
      Z(K)=-DEL1*EXP(.693147*FLOAT(K-2))
   30 ZZ(K)=-DEL1*EXP(.693147*(FLOAT(K)-1.5))
      DO 40 K=KL1-1,KL2+1
      Z(K)=-(FLOAT(K)-CC)/BB
   40 ZZ(K)=-(FLOAT(K)-CC+0.5)/BB
      DO 50 K=KL2+1,KBM1
      Z(K)=(1.0-DEL2*EXP(.693147*FLOAT(KB-K-1)))*(-1.)
   50 ZZ(K)=(1.0-DEL2*EXP(.693147*(FLOAT(KB-K)-1.5)))*(-1.)
      Z(KB)=-1.0
      ZZ(KBM1)=-1.*(1.-DEL2/2.)
      ZZ(KB)=-1.*(1.+DEL2/2.)
C ***
c      DO 55 K=1,KB
c      DENOM=FLOAT(KB-1)
c      Z(K)=-FLOAT(K-1)/DENOM
c      ZZ(K)=Z(K)-0.5/DENOM
c  55  CONTINUE
C ***
      DO 60 K=1,KBM1
      DZ(K)=Z(K)-Z(K+1)
      DZR(K)=1./DZ(K)
   60 DZZ(K)=ZZ(K)-ZZ(K+1)
      DZ(KB)=0.
      DZR(KB)=0.
      DZZ(KB)=0.
      WRITE(6,70)
   70 FORMAT(/2X,'K',7X,'Z',9X,'ZZ',9X,'DZ',9X,'DZZ',/)
      DO 90 K=1,KB
      WRITE(6,80) K,Z(K),ZZ(K),DZ(K),DZZ(K)
   80 FORMAT(' ',I5,4F10.3)
   90 CONTINUE
      RETURN
      END
C               
      SUBROUTINE EPRXYZ(LABEL,TIME,A,IM,ISKP,JM,JSKP,KB,DUMY)
      DIMENSION A(IM,JM,KB),NUM(48)
      DIMENSION KP(4)
      CHARACTER LABEL*(*)
      DATA KP/1,2,11,12/
      KP(3)=KB-1
      KP(4)=KB
      KEND=4
      IF(KB.LE.4) KEND=KB
C
C        THIS SUBROUTINE WRITES HORIZONTAL LAYERS OF A 3-D FIELD
C
C      TIME=TIME IN DAYS
C      A= ARRAY(IM,JM,KB) TO BE PRINTED
C      ISPL=PRINT SKIP FOR I
C      JSPL=PRINT SKIP FOR J
C
      SCALE=1.
      WRITE(6,9980) LABEL
 9980 FORMAT(' ',A30)
      WRITE(6,9981) TIME,SCALE
 9981 FORMAT(' TIME = ',F8.2,'    MULTIPLY ALL VALUES BY',F8.4)
      DO 10 KM=1,KEND 
      K=KP(KM)
      WRITE(6,30) K
 30   FORMAT(3X,/7H LAYER ,I2)
      IB=1
      IE=IB+11*ISKP
 
  50  CONTINUE
      IF(IE.GT.IM) IE=IM
      DO 100 I=IB,IE,ISKP
 100  NUM(I)=I
      WRITE(6,9990) (NUM(I),I=IB,IE,ISKP)
 9990 FORMAT(/,12I10,/)
      DO 200 J=1,JM,JSKP
      JWR=JM+1-J
      WRITE(6,9966) JWR,(A(I,JWR,K),I=IB,IE,ISKP)
 200  CONTINUE
      WRITE(6,9984)
      IF(IE.EQ.IM) GO TO 10
      IB=IB+12*ISKP
      IE=IB+11*ISKP
      GO TO 50
 9966 FORMAT(1X,I2,12(E10.2))
 9984 FORMAT(//)
 10   CONTINUE
      RETURN
      END
C                  
      SUBROUTINE FINDPSI
         INCLUDE 'comblk.h'    
        DO 9004 J=1,JM
        DO 9004 I=1,IM
 9004   PSI(I,J)=0.E0
        DO 9005 J=2,JM
        DO 9005 I=1,IMM1
 9005   PSI(I+1,J)=PSI(I,J)-VAB(I,J)*.25E0*(D(I,J)+D(I,J-1))
     1           *(DX(I,J)+DX(I,J-1))
        CALL PRXY('   PSI FROM V       ',TIME,PSI,IM,1,JM,1,1.E3)
        DO 9006 J=2,JMM1
        DO 9006 I=2,IM
 9006   PSI(I,J+1)=PSI(I,J)+.25E0*UAB(I,J)*(D(I,J)+D(I-1,J))
     1           *(DY(I,J)+DY(I-1,J))
        CALL PRXY('   PSI FROM U       ',TIME,PSI,IM,1,JM,1,1.E3)
      RETURN
      END
C          
C               
      SUBROUTINE PROFQ(DT2)
C
       INCLUDE 'comblk.h'    
      REAL KN,KAPPA        
      DIMENSION GH(IM,JM,KB),SM(IM,JM,KB),SH(IM,JM,KB)
c      DIMENSION PROD(IM,JM,KB),KN(IM,JM,KB),BOYGR(IM,JM,KB)
      DIMENSION KN(IM,JM,KB),BOYGR(IM,JM,KB)
      DIMENSION DH(IM,JM),CC(IM,JM,KB)
      EQUIVALENCE (A,SM),(C,SH),(PROD,KN),(TPS,DH)
      EQUIVALENCE (DTEF,CC,GH)
      DATA A1,B1,A2,B2,C1/0.92,16.6,0.74,10.1,0.08/
      DATA E1/1.8/,E2/1.33/,E3/1.0/
      DATA KAPPA/0.40/,SQ/0.20/,SEF/1./
      DATA GEE/9.806/
      DATA SMALL/1.E-8/
C
C
      DO 50 J=1,JM
      DO 50 I=1,IM
   50 DH(I,J)=H(I,J)+ETF(I,J)
      DO 100 K=2,KBM1
      DO 100 J=1,JM
      DO 100 I=1,IM
      A(I,J,K)=-DT2*(KQ(I,J,K+1)+KQ(I,J,K)+2.*UMOL)*.5
     1    /(DZZ(K-1)*DZ(K)*DH(I,J)*DH(I,J))
      C(I,J,K)=-DT2*(KQ(I,J,K-1)+KQ(I,J,K)+2.*UMOL)*.5
     1    /(DZZ(K-1)*DZ(K-1)*DH(I,J)*DH(I,J))
  100 CONTINUE
C***********************************************************************
C                                                                      *
C        THE FOLLOWING SECTION SOLVES THE EQUATION                     *
C        DT2*(KQ*Q2')' - Q2*(2.*DT2*DTEF+1.) = -Q2B                    *
C                                                                      *
C***********************************************************************
C------  SURFACE AND BOTTOM B.C.S ------------
      CONST1=16.6**.6666667*SEF
      DO 90 J=1,JMM1
      DO 90 I=1,IMM1
      EE(I,J,1)=0.
      GG(I,J,1)=SQRT( (.5*(WUSURF(I,J)+WUSURF(I+1,J)))**2
     1                +(.5*(WVSURF(I,J)+WVSURF(I,J+1)))**2 )*CONST1
      UF(I,J,KB)=SQRT( (.5*(WUBOT(I,J)+WUBOT(I+1,J)))**2
     1                +(.5*(WVBOT(I,J)+WVBOT(I,J+1)))**2 )*CONST1
   90 CONTINUE
C---------------------------------------------------------------------
C   RESTORE 10 DEG, 35 PPT AND 1.025 BIAS IN T, S AND RHO
C---------------------------------------------------------------------
      DO 101 K=1,KBM1
      DO 101 J=1,JM
      DO 101 I=1,IM
      TP=T(I,J,K)+10.0E0
      SP=S(I,J,K)+35.0E0
C      Calculate pressure in units of decibars
      P=-GEE*1.025*ZZ(K)*DH(I,J)*.1 
      CC(I,J,K)=1449.1+.00821*P+4.55*TP -.045*TP**2
     $                              +1.34*(SP-35.)
      CC(I,J,K)=CC(I,J,K)/SQRT((1.-.01642*P/CC(I,J,K))
     $       *(1.-0.40*P/CC(I,J,K)**2))
  101 CONTINUE
      DO 102 K=2,KBM1
      DO 102 J=1,JM
      DO 102 I=1,IM
      Q2B(I,J,K)=ABS(Q2B(I,J,K))
      Q2LB(I,J,K)=ABS(Q2LB(I,J,K))
      BOYGR(I,J,K)=GEE*((RHO(I,J,K-1)-RHO(I,J,K))/(DZZ(K-1)*DH(I,J)))      
     1 +GEE**2*2.*1.025/(CC(I,J,K-1)**2+CC(I,J,K)**2)
  102 CONTINUE
C------ CALC. T.K.E. PRODUCTION + Richardson No (rii & rff)---------------
      DO 120 K=2,KBM1
      DO 120 J=1,JMM1
      DO 120 I=1,IMM1
      shear=.25*SEF
     1       *( (U(I,J,K)-U(I,J,K-1)+U(I+1,J,K)-U(I+1,J,K-1))**2
     2         +(V(I,J,K)-V(I,J,K-1)+V(I,J+1,K)-V(I,J+1,K-1))**2 )
     3              /(DZZ(K-1)*DH(I,J))**2
      PROD(I,J,K)=KM(I,J,K)*shear
	ustar2(i,j,k)=sqrt(km(i,j,k)*prod(i,j,k))
	rff(i,j,k)=-kh(i,j,k)*boygr(i,j,k)/(prod(i,j,k)+small)
  120 PROD(I,J,K)=PROD(I,J,K)+KH(I,J,K)*BOYGR(I,J,K)
C----Cosmetic work--------------
	do j=1,jmm1
	do i=1,imm1
        rff(i,j,1)=rff(i,j,2)
        ustar2(i,j,1)=ustar2(i,j,2)
        rff(i,j,kb)=rff(i,j,kbm1)
        ustar2(i,j,kb)=0.0
	end do
	end do
      DO 110 K=2,KBM1
      DO 110 J=1,JM
      DO 110 I=1,IM
  110 DTEF(I,J,K)=Q2B(I,J,K)*SQRT(Q2B(I,J,K))/(B1*Q2LB(I,J,K)+SMALL)
      DO 140 K=2,KBM1
      DO 140 J=1,JM
      DO 140 I=1,IM
      GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))
     1    -(2.*DT2*DTEF(I,J,K)+1.) )
      EE(I,J,K)=A(I,J,K)*GG(I,J,K)
      GG(I,J,K)=(-2.*DT2*PROD(I,J,K)
     1  +C(I,J,K)*GG(I,J,K-1)-UF(I,J,K))*GG(I,J,K)
  140 CONTINUE
      DO 150 K=1,KBM1
      KI=KB-K
      DO 150 J=1,JM
      DO 150 I=1,IM
      UF(I,J,KI)=EE(I,J,KI)*UF(I,J,KI+1)+GG(I,J,KI)
  150 CONTINUE
C***********************************************************************
C                                                                      *
C        THE FOLLOWING SECTION SOLVES THE EQUATION                     *
C        DT2(KQ*Q2L')' - Q2L*(DT2*DTEF+1.) = -Q2LB                     *
C                                                                      *
C***********************************************************************
      DO 155 J=1,JM
      DO 155 I=1,IM
      EE(I,J,1)=0.0
      GG(I,J,1)=0.0
      VF(I,J,KB)=0.     
  155 CONTINUE   
      DO 160 K=2,KBM1
      DO 160 J=1,JM
      DO 160 I=1,IM
      DTEF(I,J,K) =DTEF(I,J,K)*(1.+E2*((1./ABS(Z(K)-Z(1))+
     1    1./ABS(Z(K)-Z(KB))) *L(I,J,K)/(DH(I,J)*KAPPA))**2)
      GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))
     1    -(DT2*DTEF(I,J,K)+1.))
      EE(I,J,K)=A(I,J,K)*GG(I,J,K)
      GG(I,J,K)=(DT2*(-PROD(I,J,K)
     1   *L(I,J,K)*E1)+C(I,J,K)*GG(I,J,K-1)-VF(I,J,K))*GG(I,J,K)
  160 CONTINUE
      DO 170 K=1,KBM1
      KI=KB-K
      DO 170 J=1,JM
      DO 170 I=1,IM
      VF(I,J,KI)=EE(I,J,KI)*VF(I,J,KI+1)+GG(I,J,KI)
  170 CONTINUE
      DO 180 K=2,KBM1
      DO 180 J=1,JM
      DO 180 I=1,IM
      IF(UF(I,J,K).LT.SMALL.OR.VF(I,J,K).LT.SMALL) THEN
        UF(I,J,K)=SMALL
        VF(I,J,K)=SMALL
      ENDIF
  180 CONTINUE
C***********************************************************************
C                                                                      *
C               THE FOLLOWING SECTION SOLVES FOR KM AND KH             *
C                                                                      *
C***********************************************************************
      COEF1=A2*(1.-6.*A1/B1)
      COEF2=3.*A2*B2+18.*A1*A2
      COEF3=A1*(1.-3.*C1-6.*A1/B1)
      COEF4=18.*A1*A1+9.*A1*A2
      COEF5=9.*A1*A2
C NOTE THAT SM,SH LIMIT TO INFINITY WHEN GH APPROACHES 0.0288
      DO 220 K=2,KBM1
      DO 220 J=1,JM
      DO 220 I=1,IM
      L(I,J,K)=VF(I,J,K)/UF(I,J,K)
 220  GH(I,J,K)=L(I,J,K)**2/UF(I,J,K)*BOYGR(I,J,K)              
      DO 225 J=1,JM
      DO 225 I=1,IM
      L(I,J,1)=0.   
      L(I,J,KB)=0.   
      GH(I,J,1)=0.  
 225  GH(I,J,KB)=0.  
      DO 230 K=1,KB
      DO 230 J=1,JM
      DO 230 I=1,IM
      GH(I,J,K)=MIN(GH(I,J,K),.028)
      SH(I,J,K)=COEF1/(1.-COEF2*GH(I,J,K))
      SM(I,J,K)=COEF3+SH(I,J,K)*COEF4*GH(I,J,K)
      SM(I,J,K)=SM(I,J,K)/(1.-COEF5*GH(I,J,K))
  230 CONTINUE
C
      DO 280 K=1,KB
      DO 280 J=1,JM
      DO 280 I=1,IM
      KN(I,J,K)=L(I,J,K)*SQRT(ABS(Q2(I,J,K)))
      KQ(I,J,K)=(KN(I,J,K)*.41*SM(I,J,K)+KQ(I,J,K))*.5
      KM(I,J,K)=(KN(I,J,K)*SM(I,J,K)+KM(I,J,K))*.5
      KH(I,J,K)=(KN(I,J,K)*SH(I,J,K)+KH(I,J,K))*.5
c	if(TIME.gt.5.75) then
c      KM(I,J,K)=5.0E-5
c      KH(I,J,K)=0.002
c	end if
  280 CONTINUE
      RETURN   
      END
C
      SUBROUTINE PROFT(F,WFSURF,SWRAD,FSURF,NBC,DT2)
C
	INCLUDE 'comblk.h'
      DIMENSION F(IM,JM,KB),WFSURF(IM,JM),FSURF(IM,JM),DH(IM,JM)
      DIMENSION SWRAD(IM,JM),RAD(IM,JM,KB),TR(5),EXTC(5)
      EQUIVALENCE (TPS,DH)
C
      NTP=2
C       NTP         =     1      2       3       4       5
C   JERLOV TYPE     =     I      IA      IB      II      III
      DATA TR/          .32  ,  .31  ,  .29  ,  .26  ,  .24   /
      DATA EXTC/       .037  , .042  , .056  , .073  , .127   /
C
      UMOLPR=UMOL
C***********************************************************************
C                                                                      *
C        THE FOLLOWING SECTION SOLVES THE EQUATION                     *
C         DTI2*(KH*F')'-F=-FB                                          *
C                                                                      *
C***********************************************************************
      DO 10 J=1,JM
      DO 10 I=1,IM
      DH(I,J)=H(I,J)+ETF(I,J)
   10 CONTINUE
      DO 20 K=2,KBM1
      DO 20 J=1,JM
      DO 20 I=1,IM
      A(I,J,K-1)=-DT2*(KH(I,J,K)+UMOLPR)/(DZ(K-1)*DZZ(K-1)*DH(I,J)
     1     *DH(I,J))
      C(I,J,K)=-DT2*(KH(I,J,K)+UMOLPR)/(DZ(K)*DZZ(K-1)*DH(I,J)
     1     *DH(I,J))
   20 CONTINUE
C-----------------------------------------------------------------------
C   NBC=1: SURF. B.C. IS WFSURF. NO SW RADIATIVE PENETRATION.
C   NBC=2; SURF. B.C. IS WFSURF+SWRAD*(1.-TR). TR*SWRAD PENETRATES WATER COLUMN
C   NBC=3; SURF. B.C. IS TSURF. NO SW RADIATIVE PENETRATION.
C   NBC=4; SURF. B.C. IS TSURF. TR*SWRAD PENETRATES WATER COLUMN
C
C NOTE THAT WTSURF (=WFSURF) AND SWRAD ARE NEG. VALUES WHEN WATER COLUMN IS
C            WARMING.
C-----------------------------------------------------------------------
C------------------------------------------------------------------
C     Penetrative Radiation Calculation. At the bottom any 
C     unattenuated radiation is deposited in the bottom layer.
C------------------------------------------------------------------
        DO 512 K=1,KB
        DO 512 J=1,JM
        DO 512 I=1,IM
  512   RAD(I,J,K)=0.
      IF(NBC.EQ.2.OR.NBC.EQ.4.) THEN   
        DO 511 K=1,KBM1
        DO 511 J=1,JM
        DO 511 I=1,IM
  511   RAD(I,J,K)=SWRAD(I,J)*TR(NTP)*EXP(EXTC(NTP)*Z(K)*DH(I,J))
      ENDIF
      GO TO (50,51,52,52), NBC
   50 CONTINUE
      DO 500 J=1,JM
      DO 500 I=1,IM
      EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.0)
      GG(I,J,1)=-DT2*WFSURF(I,J)/(-DZ(1)*DH(I,J))-F(I,J,1)
  500 GG(I,J,1)=GG(I,J,1)/(A(I,J,1)-1.0)
      GO TO 53
C
   51 CONTINUE
      DO 510 J=1,JM
      DO 510 I=1,IM
      EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.0)
      GG(I,J,1)=DT2*(WFSURF(I,J)+SWRAD(I,J)*(1.0-TR(NTP))
     2      +RAD(I,J,1)-RAD(I,J,2))
     1     /(DZ(1)*DH(I,J))-F(I,J,1)
  510 GG(I,J,1)=GG(I,J,1)/(A(I,J,1)-1.0)
      GO TO 53
C
   52 CONTINUE
      DO 520 J=1,JM
      DO 520 I=1,IM
      EE(I,J,1)=0.
  520 GG(I,J,1)=FSURF(I,J)
C----------------------------------------------------------------------
   53 CONTINUE
C
C----------------------------------------------------------------------
      DO 101 K=2,KBM2
      DO 101 J=1,JM
      DO 101 I=1,IM
      GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))-1.)
      EE(I,J,K)=A(I,J,K)*GG(I,J,K)
      GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-F(I,J,K)
     1     +DT2*(RAD(I,J,K)-RAD(I,J,K+1))/(DH(I,J)*DZ(K)))*GG(I,J,K)
  101 CONTINUE
C-----  BOTTOM ADIABATIC B.C. ------------------------------------------
      DO 102 J=1,JM
      DO 102 I=1,IM
  102 F(I,J,KBM1)=((C(I,J,KBM1)*GG(I,J,KBM2)-F(I,J,KBM1)
     1         +DT2*(RAD(I,J,KBM1)-RAD(I,J,KB))/(DH(I,J)*DZ(KBM1)))
     2          /(C(I,J,KBM1)*(1.-EE(I,J,KBM2))-1.))
C----------------------------------------------------------------------
      DO 105 K=2,KBM1
      KI=KB-K
      DO 105 J=1,JM
      DO 105 I=1,IM
      F(I,J,KI)=(EE(I,J,KI)*F(I,J,KI+1)+GG(I,J,KI))
  105 CONTINUE
C
      RETURN
      END
C               
      SUBROUTINE PROFU(DT2)
       INCLUDE 'comblk.h'    
      DIMENSION DH(IM,JM)
      DATA UMOLB/1.36E-6 /
C***********************************************************************
C                                                                      *
C        THE FOLLOWING SECTION SOLVES THE EQUATION                     *
C         DT2*(KM*U')'-U=-UB                                           *
C                                                                      *
C***********************************************************************
C
      DO 84 J=1,JM
      DO 84 I=1,IM
84    DH(I,J)=1.0
 
      DO 85 J=2,JM
      DO 85 I=2,IM
   85 DH(I,J)=.5E0*(H(I,J)+ETF(I,J)+H(I-1,J)+ETF(I-1,J))
      DO 90 K=1,KB
      DO 90 J=2,JM
      DO 90 I=2,IM
   90 C(I,J,K)=(KM(I,J,K)+KM(I-1,J,K))*.5E0
      DO 100 K=2,KBM1
      DO 100 J=1,JM
      DO 100 I=1,IM
      A(I,J,K-1)=-DT2*(C(I,J,K)+UMOL  )/(DZ(K-1)*DZZ(K-1)*DH(I,J)
     1     *DH(I,J))
      C(I,J,K)=-DT2*(C(I,J,K)+UMOL  )/(DZ(K)*DZZ(K-1)*DH(I,J)
     1     *DH(I,J))
  100 CONTINUE
C
      DO 1011 J=1,JM
      DO 1011 I=1,IM
      EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.E0)
      GG(I,J,1)=(-DT2*WUSURF(I,J)/(-DZ(1)*DH(I,J))-UF(I,J,1))
     1   /(A(I,J,1)-1.E0)
1011  CONTINUE
      DO 101 K=2,KBM2
      DO 101 J=1,JM
      DO 101 I=1,IM
      GG(I,J,K)=1.E0/(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))-1.)
      EE(I,J,K)=A(I,J,K)*GG(I,J,K)
      GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-UF(I,J,K))*GG(I,J,K)
  101 CONTINUE
      DO 102 J=2,JMM1
      DO 102 I=2,IM
      TPS(I,J)=0.5E0*(CBC(I,J)+CBC(I-1,J))
     1     *SQRT(UB(I,J,KBM1)**2+(.25E0*(VB(I,J,KBM1)
     2     +VB(I,J+1,KBM1)+VB(I-1,J,KBM1)+VB(I-1,J+1,KBM1)))**2)
      UF(I,J,KBM1)=(C(I,J,KBM1)*GG(I,J,KBM2)-UF(I,J,KBM1))/(TPS(I,J)
     1 *DT2/(-DZ(KBM1)*DH(I,J))-1.E0-(EE(I,J,KBM2)-1.E0)*C(I,J,KBM1))
  102 UF(I,J,KBM1)=UF(I,J,KBM1)*DUM(I,J)
      DO 103 K=2,KBM1
      KI=KB-K
      DO 103 J=2,JMM1
      DO 103 I=2,IMM1
      UF(I,J,KI)=(EE(I,J,KI)*UF(I,J,KI+1)+GG(I,J,KI))*DUM(I,J)
  103 CONTINUE
C
      DO 104 J=2,JMM1
      DO 104 I=2,IMM1
104   WUBOT(I,J)=-TPS(I,J)*UF(I,J,KBM1)
      RETURN
      END
C            
      SUBROUTINE PROFV(DT2)
       INCLUDE 'comblk.h'    
      DIMENSION DH(IM,JM)
      DATA UMOLB/1.36E-6/
C***********************************************************************
C                                                                      *
C        THE FOLLOWING SECTION SOLVES THE EQUATION                     *
C         DT2*(KM*U')'-U=-UB                                           *
C                                                                      *
C***********************************************************************
C
      DO 84 J=1,JM
      DO 84 I=1,IM
84    DH(I,J)=1.
      DO 85 J=2,JM
      DO 85 I=2,IM
   85 DH(I,J)=.5E0*(H(I,J)+ETF(I,J)+H(I,J-1)+ETF(I,J-1))
      DO 90 K=1,KB
      DO 90 J=2,JM
      DO 90 I=2,IM
   90 C(I,J,K)=(KM(I,J,K)+KM(I,J-1,K))*.5E0
      DO 100 K=2,KBM1
C
      DO 100 J=1,JM
      DO 100 I=1,IM
      A(I,J,K-1)=-DT2*(C(I,J,K)+UMOL  )/(DZ(K-1)*DZZ(K-1)*DH(I,J)
     1     *DH(I,J))
      C(I,J,K)=-DT2*(C(I,J,K)+UMOL  )/(DZ(K)*DZZ(K-1)*DH(I,J)
     1     *DH(I,J))
  100 CONTINUE
C
      DO 1001 J=1,JM
      DO 1001 I=1,IM
      EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.E0)
      GG(I,J,1)=(-DT2*WVSURF(I,J)/(-DZ(1)*DH(I,J))-VF(I,J,1))
     1   /(A(I,J,1)-1.E0)
1001  CONTINUE
      DO 101 K=2,KBM2
      DO 101 J=1,JM
      DO 101 I=1,IM
      GG(I,J,K)=1.E0/(A(I,J,K)+C(I,J,K)*(1.E0-EE(I,J,K-1))-1.E0)
      EE(I,J,K)=A(I,J,K)*GG(I,J,K)
      GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-VF(I,J,K))*GG(I,J,K)
  101 CONTINUE
      DO 102 J=2,JMM1
      DO 102 I=2,IMM1
      TPS(I,J)=0.5E0*(CBC(I,J)+CBC(I,J-1))
     1     *SQRT((.25E0*(UB(I,J,KBM1)+UB(I+1,J,KBM1)
     2     +UB(I,J-1,KBM1)+UB(I+1,J-1,KBM1)))**2+VB(I,J,KBM1)**2)
      VF(I,J,KBM1)=(C(I,J,KBM1)*GG(I,J,KBM2)-VF(I,J,KBM1))/(TPS(I,J)
     1  *DT2/(-DZ(KBM1)*DH(I,J))-1.E0-(EE(I,J,KBM2)-1.E0)*C(I,J,KBM1))
  102 VF(I,J,KBM1)=VF(I,J,KBM1)*DVM(I,J)
      DO 103 K=2,KBM1
      KI=KB-K
      DO 103 J=2,JMM1
      DO 103 I=2,IMM1
      VF(I,J,KI)=(EE(I,J,KI)*VF(I,J,KI+1)+GG(I,J,KI))*DVM(I,J)
  103 CONTINUE
C
      DO 104 J=2,JMM1
      DO 104 I=2,IMM1
104   WVBOT(I,J)=-TPS(I,J)*VF(I,J,KBM1)
      RETURN
      END
C                  
      SUBROUTINE PRXY (LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA)
C >>>
C     THIS WRITES A 2-D FIELD
C     TIME=TIME IN DAYS
C     A = ARRAY(IM,JM) TO BE PRINTED
C     ISKP=PRINT SKIP FOR I
C     JSKP=PRINT SKIP FOR J
C     SCALE=DIVISOR FOR VALUES OF A
C
C   IMPLICIT HALF PRECISION (A-H,O-Z)
      DIMENSION A(IM,JM),NUM(150),LINE(150)
      CHARACTER LABEL*(*)
      DATA ZERO /1.E-12/
C
      SCALE=SCALA
      IF (SCALE.GT.ZERO) GO TO 160
      AMX=ZERO
      DO 150 J=1,JM,JSKP
      DO 150 I=1,IM,ISKP
      AMX=MAX(ABS(A(I,J)),AMX)
  150 CONTINUE
      IF(AMX.EQ.0.) THEN   
       SCALEI=0.   
       GOTO 165  
      ENDIF    
      SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103)
  160 CONTINUE
      SCALEI=1.E0/SCALE
  165 CONTINUE
      WRITE(6,170) LABEL
  170 FORMAT(1X,A40)
      WRITE(6,180) TIME,SCALE
  180 FORMAT(' TIME =',F9.4,' DAYS     MULTIPLY ALL VALUES BY',1PE10.3)
      DO 190 I=1,IM
  190 NUM(I)=I
      IB=1
C
  200 CONTINUE
      IE=IB+23*ISKP
      IF(IE.GT.IM) IE=IM
      WRITE(6,210) (NUM(I),I=IB,IE,ISKP)
  210 FORMAT(/,2X,24I5,/)
      DO 260 J=1,JM,JSKP
      JWR=JM+1-J
      DO 220 I=IB,IE,ISKP
  220 LINE(I)=INT(SCALEI*A(I,JWR))
      WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP)
 
  240 FORMAT(1X,I3,24I5)
  260 CONTINUE
      WRITE(6,280)
  280 FORMAT(//)
      IF(IE.GE.IM) RETURN
      IB=IB+24*ISKP
      GO TO 200
      END
      SUBROUTINE PRXYZ (LABEL,TIME,A,IM,ISKP,JM,JSKP,KB,SCALA)
C >>>
C
C     THIS WRITES HORIZONTAL LAYERS OF A 3-D FIELD
C     TIME=TIME IN DAYS
C     A = ARRAY(IM,JM,KB) TO BE PRINTED
C     ISKP=PRINT SKIP FOR I
C     JSKP=PRINT SKIP FOR J
C     SCALE=DIVISOR FOR VALUES OF A
C
      DIMENSION A(IM,JM,KB),NUM(150),LINE(150),KP(4)
      CHARACTER LABEL*(*)
      DATA KE,KP /4, 1,2,3,4/
      DATA ZERO /1.E-10/
C
      SCALE=SCALA
      IF (SCALE.GT.ZERO) GO TO 150
      AMX=ZERO
      DO 140  KM=1,KE
      K=KP(KM)
      DO 140 J=1,JM,JSKP
      DO 140 I=1,IM,ISKP
      AMX=MAX(ABS(A(I,J,K)),AMX)
  140 CONTINUE
      IF(AMX.EQ.0.) THEN   
       SCALEI=0.   
       GOTO 165  
      ENDIF    
      SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103)
  150 CONTINUE
      SCALEI=1.E0/SCALE
  165 CONTINUE
      WRITE(6,160) LABEL
  160 FORMAT(1X,A40)
      WRITE(6,170) TIME,SCALE
  170 FORMAT(' TIME = ',F9.4,' DAYS     MULTIPLY ALL VALUES BY',1PE10.3)
      DO 180 I=1,IM
  180 NUM(I)=I
C
      DO 500 KM=1,KE
      K=KP(KM)
      WRITE(6,190) K
  190 FORMAT(3X,/7H LAYER ,I2)
      IB=1
C
  200 CONTINUE
      IE=IB+23*ISKP
      IF(IE.GT.IM) IE=IM
      WRITE(6,220) (NUM(I),I=IB,IE,ISKP)
  220 FORMAT(/,2X,24I5,/)
      DO 260 J=1,JM,JSKP
      JWR=JM+1-J
      DO 230 I=IB,IE,ISKP
  230 LINE(I)=INT(SCALEI*A(I,JWR,K))
      WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP)
  240 FORMAT(1X,I3,24I5)
  260 CONTINUE
      WRITE(6,280)
  280 FORMAT(//)
      IF(IE.GE.IM) GO TO 500
      IB=IB+24*ISKP
      GO TO 200
  500 CONTINUE
      RETURN
      END
C               
      SUBROUTINE PRXZ(LABEL,TIME,A,IM,ISKP,JM,J1,J2,KB,SCALA,DT,ZZ)
      DIMENSION A(IM,JM,KB),NUM(100),LINE(100),IDT(100),
     1        ZZ(KB),DT(IM,JM)
      CHARACTER LABEL*(*)
      DATA ZERO/1.E-10/
C
C        THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD
C
C      TIME=TIME IN DAYS
C      A= ARRAY(IM,JM,KB) TO BE PRINTED
C      ISPL=PRINT SKIP FOR I
C      JSPL=PRINT SKIP FOR J
C      SCALE=DIVISOR FOR VALUES OF A
C
      SCALE=SCALA
      IF (SCALE.GT.ZERO) GO TO 150
      AMX=ZERO
      DO 140  K=1,KB
      DO 140 J=1,JM
      DO 140 I=1,IM,ISKP
      AMX=MAX(ABS(A(I,J,K)),AMX)
  140 CONTINUE
      IF(AMX.EQ.0.) THEN   
       SCALEI=0.   
       GOTO 165  
      ENDIF    
      SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103)
  150 CONTINUE
      SCALEI=1./SCALE
  165 CONTINUE
      WRITE(6,9980) LABEL
 9980 FORMAT(' ',A30)
      WRITE(6,9981) TIME,SCALE
 9981 FORMAT(' TIME = ',F8.2,'    MULTIPLY ALL VALUES BY',E8.2)
      DO 10 JPP=1,2
      J=J1
      IF(JPP.EQ.2) J=J2
      IF(J.EQ.0) RETURN
      WRITE(6,30) J
 30   FORMAT(3X,/' SECTION J =',I3)
      IB=1
      IE=IB+23*ISKP
  50  CONTINUE
      IF(IE.GT.IM) IE=IM
      DO 100 I=IB,IE,ISKP
      IDT(I)=DT(I,J)
 100  NUM(I)=I
      WRITE(6,9990) (NUM(I),I=IB,IE,ISKP)
 9990 FORMAT(/,'    I =  ',24I5,/)
      WRITE(6,9991) (IDT(I),I=IB,IE,ISKP)
 9991 FORMAT(9X,'D =',24I5.0,/,'       ZZ ')
      DO 200 K=1,KB
      DO 190 I=IB,IE,ISKP
  190 LINE(I)=SCALEI*A(I,J,K)
      WRITE(6,9966) K,ZZ(K),(LINE(I),I=IB,IE,ISKP)
 9966 FORMAT(1X,I2,2X,F6.3,24I5)
 200  CONTINUE
      WRITE(6,9984)
      IF(IE.EQ.IM) GO TO 10
      IB=IB+24*ISKP
      IE=IB+23*ISKP
      GO TO 50
 9984 FORMAT(//)
 10   CONTINUE
      RETURN
      END
C               
      SUBROUTINE PRYZ(LABEL,TIME,A,IM,ISKP,JM,JSKP,KB,SCALA,DT,ZZ)
      DIMENSION A(IM,JM,KB),NUM(200),PLINE(200),DT(IM,JM),ZZ(KB)
      DIMENSION IP(2)
      CHARACTER LABEL*(*)
      DATA ZERO/1.E-10/
      IP(1)=IM/2
      IP(2)=IM-3
C
C        THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD
C
C      TIME=TIME IN DAYS
C      A= ARRAY(IM,JM,KB) TO BE PRINTED
C      ISPL=PRINT SKIP FOR I
C      JSPL=PRINT SKIP FOR J
C      SCALE=DIVISOR FOR VALUES OF A
C
      SCALE=SCALA
      IF (SCALE.GT.ZERO) GO TO 150
      AMX=ZERO
      DO 140  K=1,KB
      DO 140 J=1,JM,JSKP
      DO 140 I=1,IM,ISKP
      AMX=MAX(ABS(A(I,J,K)),AMX)
  140 CONTINUE
      IF(AMX.EQ.0.) THEN   
       SCALEI=0.   
       GOTO 165  
      ENDIF    
      SCALE=10.E0**(INT( LOG10(AMX)+1.E2)-103)
  150 CONTINUE
      SCALEI=1./SCALE            
  165 CONTINUE
      WRITE(6,9980) LABEL
 9980 FORMAT(' ',A30)
      WRITE(6,9981) TIME,SCALE
 9981 FORMAT(' TIME = ',F8.2,'    MULTIPLY ALL VALUES BY',E8.2)
      DO 10 IPP=1,2
      I=IP(IPP)
      WRITE(6,30) I
 30   FORMAT(3X,/' SECTION I =',I3)
      JB=1
      JE=JB+23*JSKP
  50  CONTINUE
      IF(JE.GT.JM) JE=JM
      DO 100 J=JB,JE,JSKP
 100  NUM(J)=J
      WRITE(6,9990) (NUM(J),J=JB,JE,JSKP)
 9990 FORMAT(/,'    J =  ',24I5,/)
      WRITE(6,9991) (DT(I,J),J=JB,JE,JSKP)
 9991 FORMAT(7X,'D =',24F5.0,/,'       ZZ ')
      DO 200 K=1,KB
      DO 190 J=JB,JE,JSKP
  190 PLINE(J)=SCALEI*A(I,J,K)
      WRITE(6,9966) K,ZZ(K),(PLINE(J),J=JB,JE,JSKP)
 9966 FORMAT(1X,I2,2X,F6.3,24(F5.1))
 200  CONTINUE
      WRITE(6,9984)
      IF(JE.EQ.JM) GO TO 10
      JB=JB+24*JSKP
      JE=JB+23*JSKP
      GO TO 50
 9984 FORMAT(//)
 10   CONTINUE
      RETURN
      END
c             
      SUBROUTINE VERTVL(DTI2)
       INCLUDE 'comblk.h'    
      DIMENSION XFLUX(IM,JM,KB),YFLUX(IM,JM,KB)
      EQUIVALENCE (XFLUX,A),(YFLUX,C)
C
C
C CALCULATE NEW VERTICAL VELOCITY
C
C REESTABLISH BOUNDARY CONDITIONS
      DO 100 K=1,KBM1
      DO 100 J=2,JM
      DO 100 I=2,IM
 100  XFLUX(I,J,K)
     1 =.25E0*(DY(I,J)+DY(I-1,J))*(DT(I,J)+DT(I-1,J))*U(I,J,K)
      DO 120 K=1,KBM1
      DO 120 J=2,JM
      DO 120 I=1,IM
  120 YFLUX(I,J,K)
     1  =.25E0*(DX(I,J)+DX(I,J-1))*(DT(I,J)+DT(I,J-1))*V(I,J,K)
C
      DO 125 J=1,JM
      DO 125 I=1,IM
125   W(I,J,1)=0.E0
      DO 710 K=1,KBM1
      DO 710 J=2,JMM1
      DO 710 I=1,IMM1
 710  W(I,J,K+1)=W(I,J,K)
     1    +DZ(K)*((XFLUX(I+1,J,K)-XFLUX(I,J,K)
     2            +YFLUX(I,J+1,K)-YFLUX(I,J,K))/(DX(I,J)*DY(I,J))
     3                        +(ETF(I,J)-ETB(I,J))/DTI2 )
      RETURN
      END

      SUBROUTINE LINTERP(XA,YA,N,X,Y)

C Linear interpolation
C Given arrays XA(N) and YA(N), LINTERP finds the value of Y
C   for a given X

      DIMENSION XA(N),YA(N),XT(1000),YT(1000)

C Order the data
      if(XA(1).gt.XA(N)) then
       do 20 i=1,N
       XT(i)=XA(i)
20     YT(i)=YA(i)
       do 21 i=1,N
       ii=N+1-i
       XA(i)=XT(ii)
21     YA(i)=YT(ii)
      endif
 
C Find the array points that X lies between
      jl=0
10    jl=jl+1
      if(X.le.XA(jl)) goto 11
      if(jl.ge.(N+1)) then
        jl=N+1
        goto 11
      endif
      goto 10
 
11    J=jl-1
      if(J.eq.N) then
       Y=YA(N)
       RETURN
      endif
      if(J.eq.0) then
       Y=YA(1)
       RETURN
      endif
 
      Y=YA(J)+(X-XA(J))*(YA(J)-YA(J+1))/(XA(J)-XA(J+1))
      RETURN
      END


	SUBROUTINE ADVTSML(FB,F,FCLIM,DTI2,FF,NItera)
C**********************************************************************
C   THIS SUBROUTINE CALCULATES THE 3D ADVECTION OF PASSIVE SCALAR     *
C   TRACERS BY A FIRST-ORDER UPSTREAM SCHEME. TO REDUCE IMPLICIT      *
C   DIFFUSION THE SMOLARKIEWICZ ITERATIVE UPSTREAM SCHEME WITH A      *
C   SPECIAL DEFINED ANTIDIFFUSIVE VELOCITY IS USED.                   *
C                     			                              *
C       Gianmaria Sannino'          Vincenzo Artale''                 *
C                                                                     *
C  '  Inter-University Computing Consortium (CASPUR)                  *
C        Universita' "La Sapienza", Rome, ITALY                       *
C                                                                     *
C ''  Italian National Agency for New Technology and Environment      *
C     (ENEA C.R. Casaccia), Rome, ITALY                               *
C                                                                     *
C   Input:	 		       			              *
C    FB -----> As Subroutine ADVT in POM         	              *
C    F ------> As Subroutine ADVT in POM          	              *
C    FCLIM --> As Subroutine ADVT in POM        	              *
C    DTI2 ---> As Subroutine ADVT in POM        	              *
C    NItera -> Number of iteration		                      *
C   Output:	 		       		                      *
C    FF -----> As Subroutine ADVT in POM source code	              *
C                                                                     *
C  Smolarkiewicz, P.K. "A Fully Multidimensional Positive Definite    *
C   Advection Transport Algorithm with Small Implicit Diffusion"      *
C   Journal of Computational Physics 54, 325-362 (1984)               *
C                                                                     *
C  Notes for POM users:                                               *
C  1. replaces subroutine ADVT                                        *
C  2. suggested number of iterations: 1 to 4 (1 is standard upstream; *
C     more iterations are less diffusive; 3 itr. adds 50% to POM)     *
C  3. can run with no additional diff. (TPRNU=0.)                     *
C  4. the shock switch option may not work for most cases, to ignore  * 
C     it set SW=Constant (should be 1, but smaller value gives        * 
C     smoother solution with less overshoot for NItera>1 ).           * 
C                                                                     * 
C  (Try on your own risk, no advice or support will be provided) T.E. *
C**********************************************************************
C
      INCLUDE 'comblk.h'
C
      DIMENSION FB(IM,JM,KB),F(IM,JM,KB),FF(IM,JM,KB)
      DIMENSION FCLIM(IM,JM,KB),FBMem(IM,JM,KB),Eta(IM,JM)
      DIMENSION XMassFlux(IM,JM,KB),YMassFlux(IM,JM,KB),ZWFlux(IM,JM,KB)
      DIMENSION XFlux(IM,JM,KB),YFlux(IM,JM,KB),ZFlux(IM,JM,KB)
C
C-----------------------------------------------------------------------
C                 HORIZONTAL MASS FLUXES: DY*U*D  &  DX*V*D
C-----------------------------------------------------------------------
      DO 100 K=1,KB
      DO 100 J=1,JM
      DO 100 I=1,IM
      XMassFlux(I,J,K)=0.
      YMassFlux(I,J,K)=0.
 100  CONTINUE
C
      DO 110 K=1,KBM1
      DO 111 J=2,JMM1
      DO 111 I=2,IM
      XMassFlux(I,J,K)=0.25*(DY(I-1,J)+DY(I,J))*(DT(I-1,J)+DT(I,J))*
     &                  U(I,J,K)
  111 CONTINUE
      DO 112 J=2,JM
      DO 112 I=2,IMM1
      YMassFlux(I,J,K)=0.25*(DX(I,J-1)+DX(I,J))*(DT(I,J-1)+DT(I,J))*
     &                   V(I,J,K)
  112 CONTINUE
  110 CONTINUE
C
      DO 120 J=1,JM
      DO 120 I=1,IM
	FB(I,J,KB)=FB(I,J,KBM1)
  120 Eta(I,J)=ETB(I,J)
C
      DO 130 K=1,KB
      DO 130 J=1,JM
      DO 130 I=1,IM
	ZWFlux(I,J,K)=W(I,J,K)
  130 FBMem(I,J,K)=FB(I,J,K)
C-----------------------------------------------------------------------
C   		       START
C	   	 SMOLARKIEWICZ SCHEME 
C-----------------------------------------------------------------------
	DO 999 ITERA=1,NItera
C-----------------------------------------------------------------------
C           UPWIND ADVECTION SCHEME 
C-----------------------------------------------------------------------
      DO 140 K=1,KBM1
      DO 140 J=2,JM
      DO 140 I=2,IM
      XFlux(I,J,K)=0.5*(
     &             (XMassFlux(I,J,K)+ABS(XMassFlux(I,J,K)))
     &		   *FBMem(I-1,J,K)+
     &             (XMassFlux(I,J,K)-ABS(XMassFlux(I,J,K)))
     &	         *FBMem(I,J,K)) 
C
      YFlux(I,J,K)=0.5*(
     &             (YMassFlux(I,J,K)+ABS(YMassFlux(I,J,K)))
     &		    *FBMem(I,J-1,K)+
     &	       (YMassFlux(I,J,K)-ABS(YMassFlux(I,J,K)))
     &		    *FBMem(I,J,K)) 
  140 CONTINUE
C
      DO 150 J=2,JMM1
      DO 150 I=2,IMM1
      ZFlux(I,J,1)=0.0
  150 ZFlux(I,J,KB)=0.0
C
      DO 160 K=2,KBM1
      DO 160 J=2,JMM1
      DO 160 I=2,IMM1
      ZFlux(I,J,K)=0.5*(
     &		 (ZWFlux(I,J,K)+ABS(ZWFlux(I,J,K)))
     &                  *FBMem(I,J,K)+
     &             (ZWFlux(I,J,K)-ABS(ZWFlux(I,J,K)))
     &                 *FBMem(I,J,K-1))
  160 ZFlux(I,J,K)=ZFlux(I,J,K)*ART(I,J)
C-----------------------------------------------------------------------
C        ADD NET ADVECTIVE FLUXES; THEN STEP FORWARD IN TIME 
C-----------------------------------------------------------------------
      DO 170 K=1,KBM1
      DO 170 J=2,JMM1
      DO 170 I=2,IMM1
      FF(I,J,K)=     XFlux(I+1,J,K)-XFlux(I,J,K)
     &              +YFlux(I,J+1,K)-YFlux(I,J,K)
     &             +(ZFlux(I,J,K)-ZFlux(I,J,K+1))/DZ(K)
  170 FF(I,J,K)=(FBMem(I,J,K)*(H(I,J)+Eta(I,J))*ART(I,J)-
     &   DTI2*FF(I,J,K))/((H(I,J)+ETF(I,J))*ART(I,J))
C-----------------------------------------------------------------------
C                ANTIDIFFUSION VELOCITY 
C-----------------------------------------------------------------------
      CALL SMOL_ADIF(XMassFlux,YMassFlux,ZWFlux,FB,FF,DTI2)
C
      DO 180 J=1,JM
      DO 180 I=1,IM
      Eta(I,J)=ETF(I,J)
      DO 180 K=1,KB
      FBMem(I,J,K)=FF(I,J,K)
  180 CONTINUE
C-----------------------------------------------------------------------
C   		          END
C	   	 SMOLARKIEWICZ SCHEME 
C-----------------------------------------------------------------------
  999 CONTINUE
C-----------------------------------------------------------------------
C         ADD HORIZONTAL DIFFUSIVE FLUXES 
C-----------------------------------------------------------------------
      DO 190 K=1,KB   
      DO 190 J=1,JM
      DO 190 I=1,IM
 190  FB(I,J,K)=FB(I,J,K)-FCLIM(I,J,K)
C
      DO 200 K=1,KBM1
      DO 200 J=2,JM
      DO 200 I=2,IM
      XMassFlux(I,J,K)=0.5*(AAM(I,J,K)+AAM(I-1,J,K))
      YMassFlux(I,J,K)=0.5*(AAM(I,J,K)+AAM(I,J-1,K))
  200 CONTINUE
C 
      DO 210 K=1,KBM1
      DO 210 J=2,JM
      DO 210 I=2,IM
      XFlux(I,J,K)=
     &    -XMassFlux(I,J,K)*(H(I,J)+H(I-1,J))*TPRNU
     &    *(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J)/(DX(I,J)+DX(I-1,J))
     &    *0.5*(DY(I,J)+DY(I-1,J))
      YFlux(I,J,K)=
     &    -YMassFlux(I,J,K)*(H(I,J)+H(I,J-1))*TPRNU
     &    *(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J)/(DY(I,J)+DY(I,J-1))
     &    *0.5*(DX(I,J)+DX(I,J-1))
  210 CONTINUE
C
      DO 220 K=1,KB   
      DO 220 J=1,JM
      DO 220 I=1,IM
 220     FB(I,J,K)=FB(I,J,K)+FCLIM(I,J,K)
C
C----- ADD NET HORIZONTAL FLUXES; THEN STEP FORWARD IN TIME --------
      DO 230 K=1,KBM1
      DO 230 J=2,JMM1
      DO 230 I=2,IMM1
  230 FF(I,J,K)=FF(I,J,K)-DTI2*(
     &              +XFlux(I+1,J,K)-XFlux(I,J,K)
     &              +YFlux(I,J+1,K)-YFlux(I,J,K))
     &                    /((H(I,J)+ETF(I,J))*ART(I,J))
C--------------------------------------------------------------------
      RETURN
      END