C***********************************************************************
C     pom97.f for island/seamount problem with eastward flow 
C     (Set SEAMT=.TRUE.) or user provided problem (Set SEAMT=
C     .FALSE.).
C     This code contains two versions of S.R. BCOND. Check
C     the subroutines for details.
C ---------------------------------------------------------
C Keep record of subroutines that are "commented out".
C      none
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 BEEN CONVERTED TO STANDARD FORTRAN 77 BY STEVE       *
C BRENNER WITH SOME ADDITIONAL CHANGES BY ME IN PROFQ WHICH IS NOW      *
C SOMEWHAT SIMPLER. CODE RESTORED TO SINGLE PRECISION.                  *
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 April 5, 2000 *
C                                                                       *
C                                                                       *
C  Copyright 1996 The Trustees of Princeton University                  *
C                                                                       *
C This program is free software; you can redistribute it and/or         *
C modify it under the terms of the GNU General Public License           *
C as published by the Free Software Foundation, either version 2        *  
C of the License, or (at your option) any later version.                *
C                                                                       *
C This program is distributed in the hope that it will be useful,       *
C but WITHOUT ANY WARRANTY; without even the implied warranty of        *
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
C GNU General Public License for more details.                          *
C                                                                       *
C A copy of the GNU General Public License is available at              *
C        http://www.gnu.org/copyleft/gpl.html#SEC3                      *
C or by writing to The Free Software Foundation, Inc., 59 Temple Place  *
C Suite 330, Boston, MA  02111-1307, USA.                               *
C                                                                       *
C************************************************************************
C
         INCLUDE 'comblk97.h'    
      LOGICAL SEAMT
      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),DRX2D(IM,JM),DRY2D(IM,JM),
     3     TCLIM(IM,JM,KB),SCLIM(IM,JM,KB),ADX2D(IM,JM),ADY2D(IM,JM),
     4     SWRAD(IM,JM)
C--------------------------------------------------------------------
      REAL ISPI,ISP2I
      DATA PI/3.141592654/,SMALL/1.E-10/
C     OPEN(40,FILE='INCOND ',FORM='UNFORMATTED')
      GRAV=9.806
      TBIAS=0.
      SBIAS=0.
      IINT=0
      TIME=0.
C
      SEAMT=.TRUE.
      DTE=6.    
      ISPLIT=30     
      DAYS=5. 
      PRTD1=5.   
      ISPADV=5
      SMOTH=0.10
      HORCON=0.2
C  N.B.: TPRNI=0. yields zero  horizontal diffusivity !!
      TPRNI=.0   
      UMOL=2.E-5
      MODE=3
      NREAD=0
C
      ISWTCH=100000
      IPRTD2=1
C  One can create a 'params' file in a runscript to overwrite parameters.
C     INCLUDE 'params'
C
C--------------------------------------------------------------------------
C    Herewith is a list of parameters that are user discretionary.
C--------------------------------------------------------------------------
C                      In comblk97.h
C       IM, JM, KB = grid dimension
C
C                      In MAIN
C       SEAMT=.T. The ready-to-run seamount probelem will be executed.
C       SEAMT=.F. Code expects user prepared input for READ(40) which
C              will overwrite seamount problem.
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       DTE = EXTERNAL TIME STEP
C       DTI = INTERNAL TIME STEP
C       ISPLIT = DTI/DTE
C       DAYS = LENGTH OF RUN IN DAYS
C       PRTD1 = PRINT INTERVAL IN DAYS; AFTER ISWTCH DAYS, PRINT
C                  INTERVAL = IPRTD2
C       NREAD=0, NO RESTART INPUT FILE; NREAD=1, RESTART
C       ISPADV = STEP INTERVAL WHEREIN EXTERNAL MODE ADVECTIVE TERMS ARE
C                  NOT UPDATED
C       HORCON = CONSTANT IN SMAGORINSKY HORIZONTAL VISCOSITY  
C       TPRNI = HORIZONTAL DIFFUSIVITY/VISCOSITY, INVERSE PRANDTL NUMBER
C                (see S.R. ADVT)
C       UMOL = Background viscosity used in S.R. PROFT, PROFU and PROFV.       
C       SMOTH = CONSTANT IN TIME SMOOTHER TO PREVENT SOLUTION SPLITTING
C       TBIAS, SBIAS = Value set in data statement. For 32 bit arithmetic,
C                      setting S=35. and T=10. will reduce round-off error.
C       RA, HO, DELH = Island/seamount topographical parameters
C       RAMP = Device to slowly ramp up wind stress and baroclinic forcing.
C              It can be used to eliminate initial inertial oscillations.
C              However, RAMP is set to a constant 1.0 below and is 
C              therefor disabled.
C
C                       In BCOND
C       HMAX = maximum depth
C  
C                       In DEPTH
C       KL1, KL2  = # of logarithmic grid points in surface and bottom
C                   respectively
C
C                       In PROFT
C       NBC = choice of surface boundary conditions
C       NTP = choice of radiative penetration parameter.
C
C                       In SLPMIN
C       SLMIN = see subroutine  
C--------------------------------------------------------------------------
      DTI=DTE*FLOAT(ISPLIT)
      DTE2=DTE*2
      DTI2=DTI*2
      IEND=DAYS*24*3600/DTI+2
      IPRINT=PRTD1*24*3600/DTI
      ISWTCH=ISWTCH*24*3600/INT(DTI)
      write(6,'(''DAYS ='',F10.2)') DAYS
      WRITE(6,7030) MODE,DTE,DTI,IEND,ISPLIT,ISPADV,IPRINT,SMOTH,
     1   HORCON,TPRNI   
 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,'   TPRNI ='F7.3,/)
      WRITE(6,'('' NREAD ='',I5,//)') NREAD
C
      ISPI=1.E0/FLOAT(ISPLIT)
      ISP2I=1.E0/(2.E0*FLOAT(ISPLIT))
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----------------------------------------------------------------------
      IP=IM
      ISK=2
      JSK=2
C
C     READ IN GRID DATA AND INITIAL AND LATERAL BOUNDARY CONDITIONS
C------------------------------------------------------------------------
      IF(SEAMT) THEN         
        CALL SEAMOUNT(TSURF,SSURF,TCLIM,SCLIM)
      ELSE
      REWIND(40)
      READ(40) KBB,Z,ZZ,DZ,DZZ,IMM,JMM,ALON,ALAT,DX,DY,H,FSM,
     1    DUM,DVM,ART,ARU,ARV,COR,TB,SB,UAB,VAB,ELB,ETB,UB,VB,
     2    TCLIM,SCLIM,RMEAN,TBE,TBW,TBN,TBS,SBE,SBW,SBN,SBS,
     3    UABE,UABW,VABN,VABS,RFE,RFW,RFN,RFS
      ENDIF
      TIME0=0.
C
C********************************************************************
C  Note that lateral thermodynamic boundary conditions are often set equal
C  to the initial conditions and are held constant thereafter. Users can
C  of course create varable boundary conditions.
C********************************************************************
C
C  Additional initial conditions
      DO 52 K=1,KB
      DO 52 J=1,JM
      DO 52 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)
      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))
   52 CONTINUE
      DO 62 I=1,IM
      DO 62 J=1,JM
      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)
      DT(I,J)=H(I,J)+ET(I,J)
      DRX2D(I,J)=0.
      DRY2D(I,J)=0.
   62 CONTINUE
      DO 63 K=1,KBM1
      DO 63 I=1,IM
      DO 63 J=1,JM
      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)
      V(I,J,K)=VB(I,J,K)
   63 CONTINUE
C
      CALL DENS(S,T,RHO)
      RAMP=1.
      CALL BAROPG(DRHOX,DRHOY)
      DO K=1,KBM1
      DO J=1,JM
      DO I=1,IM 
      DRX2D(I,J)=DRX2D(I,J)+DRHOX(I,J,K)*DZ(K)
      DRY2D(I,J)=DRY2D(I,J)+DRHOY(I,J,K)*DZ(K)
      ENDDO
      ENDDO
      ENDDO
C A very empirical specification of the bottom roughness
C          parameter follows
      DO 45 J=1,JM
      DO 45 I=1,IM
      DT(I,J)=H(I,J)
      Z0B=.01
      CBCMIN=.0025E0
  45  CBC(I,J)=MAX(CBCMIN,.16E0/LOG((ZZ(KBM1)-Z(KB))*H(I,J)
     1        /Z0B)**2)*FSM(I,J)
C
        CALL PRXY('  TOPOGRAPHY 2 ',TIME, H,IP,1,JM,1,10.)
        CALL PRXY(' FSM            ',TIME,FSM,IP,1,JM,1,1.)
        CALL PRXY(' DUM            ',TIME,DUM,IP,1,JM,1,1.)
        CALL PRXY(' DVM            ',TIME,DVM,IP,1,JM,1,1.)
c       CALL PRXY(' CBC            ',TIME,CBC,IP,1,JM,1,0.)
C   Evaluate external CFL time step
        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)
        CALL PRXY(' EXT. CFL TIME STEP ',TIME,TPS,IP,ISK,JM,JSK,1.E0)
c       CALL PRXY(' COR           ',TIME,COR,IP,ISK,JM,1,0.E0)
C
C
C  THE FOLLOWING DATA ARE NEEDED FOR A SEAMLESS RESTART. IF NREAD=1,
C  DATA HAD BEEN CREATED BY A PREVIOUS RUN (See WRITE(140) at end of
C  this program). NREAD=0 DENOTES A FIRST TIME RUN.      
      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,ADX2D,ADY2D,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
      CALL PRXYZ(' TEMP.      ',TIME,TB,IP,ISK,JM,JSK,KB,0.0E0)
      CALL PRXYZ(' SAL.       ',TIME,SB,IP,ISK,JM,JSK,KB,0.0E0)
C     CALL PRXYZ('RHO         ',TIME,RHO,IP,ISK,JM,1,KB,0.E0)
      CALL PRXZ('TB' ,TIME,TB,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
      CALL PRXZ('RHO  ' ,TIME,RHO  ,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
      CALL PRXZ('RMEAN' ,TIME,RMEAN,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c     CALL PRXZ('DRHOX' ,TIME,DRHOX,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c     CALL PRXZ('DRHOY' ,TIME,DRHOY,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
      CALL PRXY(' DRX2D      ',TIME,DRX2D,IP,ISK,JM,JSK,0.)
      CALL PRXY(' DRY2D      ',TIME,DRY2D,IP,ISK,JM,JSK,0.)
c     CALL PRXY('  FREE SURFACE ',TIME, ELB,IP,ISK,JM,JSK,1.E-3)
c     CALL FINDPSI
C     CALL PRXY('  VAB          ',TIME, VAB,IP,ISK,JM,JSK,0.)
      PERIOD=(2.E0*PI)/ABS(COR(IM/2,JM/2))/86400.
C***********************************************************************
C                                                                      *
C                  BEGIN NUMERICAL INTEGRATION                         *
C                                                                      *
C***********************************************************************
C
                     DO 9000 IINT=1,IEND
C     IF(TIME.GT.20) MODE=3
C
      TIME=DTI*FLOAT(IINT)/86400.+TIME0
      RAMP=TIME/PERIOD
      IF(RAMP.GT.1.E0) RAMP=1.
      RAMP=1.
C     write(6,'('' MODE,IINT,TIME ='',2I5,F9.2)') MODE,IINT,TIME      
C------------------------------------------------------------------------
C     SET TIME DEPENDENT, SURFACE AND LATERAL BOUNDARY CONDITIONS.
C     THE LATTER WILL BE USED IN SUBROUTINE BCOND. Users may
C     wish to create a subroutine to supply WUSURF, WVSURF, WTSURF,
C     WSSURF and SWRAD.
C------------------------------------------------------------------------
C    INTRODUCE SIMPLE WIND STRESS, VALUE IS NEGATIVE FOR WESTERLY OR 
C    SOUTHERLY WINDS
C  The following wind stress has been tapered along the boundary
C  to suppress numerically induced oscilations near the boundary (Jamart
C  &Ozer,J.G.R.,91, 10621-10631). To make a healthy surface Ekman layer
C  it would be well to set KL1=9 in S.R. DEPTH.
      DO 85 J=2,JMM1
      DO 85 I=2,IMM1
c     WUSURF(I,J)=RAMP*(1.E-4*COS(PI*(J-1)/JMM1))
c    1         *.25*(DVM(I,J+1)+DVM(I-1,J+1)+DVM(I-1,J)+DVM(I,J))
      WTSURF(I,J)=0.E0
      WUSURF(I,J)=0.
      WVSURF(I,J)=0.
   85 SWRAD(I,J)=0.E0
C------------------------------------------------------------------------
C
      IF(MODE.EQ.2) GOTO 8001
        CALL ADVCT(ADVX,ADVY)
        CALL BAROPG(DRHOX,DRHOY)
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)
   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
        ADX2D(I,J)=0.
        ADY2D(I,J)=0.
        DRX2D(I,J)=0.
        DRY2D(I,J)=0.
        AAM2D(I,J)=0.
   96   CONTINUE    
        DO 100 K=1,KBM1
        DO 100 J=1,JM
        DO 100 I=1,IM
        ADX2D(I,J)=ADX2D(I,J)+ADVX(I,J,K)*DZ(K)
        ADY2D(I,J)=ADY2D(I,J)+ADVY(I,J,K)*DZ(K)
        DRX2D(I,J)=DRX2D(I,J)+DRHOX(I,J,K)*DZ(K)
        DRY2D(I,J)=DRY2D(I,J)+DRHOY(I,J,K)*DZ(K)
        AAM2D(I,J)=AAM2D(I,J)+AAM(I,J,K)*DZ(K)
  100   CONTINUE
C       CALL PRXY(' DRX2D     ',TIME,DRX2D,IP,ISK,JM,JSK,0.0)
C ----------------------------------------------------------------------
        CALL ADVAVE(ADVUA,ADVVA,MODE)
        DO 87 J=1,JM
        DO 87 I=1,IM
        ADX2D(I,J)=ADX2D(I,J)-ADVUA(I,J)
   87   ADY2D(I,J)=ADY2D(I,J)-ADVVA(I,J)
C
 8001 CONTINUE
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
c     write(6,'('' IEXT,TIME ='',I5,F9.2)') IEXT,TIME      
      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,IM
  420 UAF(I,J)=ADX2D(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               +DRX2D(I,J)   
     6      +ARU(I,J)*( WUSURF(I,J)-WUBOT(I,J)   )
      DO 425 J=2,JMM1
      DO 425 I=2,IM
  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,JM
      DO 430 I=2,IMM1
      VAF(I,J)=ADY2D(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               +DRY2D(I,J)  
     6    + ARV(I,J)*( WVSURF(I,J)-WVBOT(I,J)   )
 430  CONTINUE
      DO 435 J=2,JM
      DO 435 I=2,IMM1
      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))
  435 CONTINUE
      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
      VMAXL=100.
      VAMAX=0.       
      DO 442 J=1,JM
      DO 442 I=1,IM
      IF(ABS(VAF(I,J)).GE.VAMAX) THEN
        VAMAX=ABS(VAF(I,J))   
	IMAX=I
	JMAX=J
      ENDIF
  442 CONTINUE
      IF(VAMAX.GT.VMAXL) GO TO 9001
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
      IF(VAMAX.GT.VMAXL) GO TO 9001
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=1,JM
      DO 300 I=1,IM
 300  TPS(I,J)=TPS(I,J)+U(I,J,K)*DZ(K)
      DO 302 K=1,KBM1
      DO 302 J=1,JM
      DO 302 I=2,IM
 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=1,JM
      DO 304 I=1,IM
 304  TPS(I,J)=TPS(I,J)+V(I,J,K)*DZ(K)
      DO 306 K=1,KBM1
      DO 306 J=2,JM
      DO 306 I=1,IM
 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----------------------------------------------------------------
C     COMPUTE TF AN SF USING UF AND VF AS TEMPORARY VARIABLES
C----------------------------------------------------------------
      IF(MODE.EQ.4) GO TO 360
      CALL ADVT(TB,T,TCLIM,DTI2,UF)
      CALL ADVT(SB,S,SCLIM,DTI2,VF)
      CALL PROFT(UF,WTSURF,SWRAD,TSURF,1,DTI2)
      CALL PROFT(VF,WSSURF,SWRAD,SSURF,1,DTI2)
      CALL BCOND(4)
C
      DO 355 K=1,KB
      DO 355 J=1,JM
      DO 355 I=1,IM
      T(I,J,K)=T(I,J,K)+.5*SMOTH*(UF(I,J,K)
     1  +TB(I,J,K)-2.*T(I,J,K))
      S(I,J,K)=S(I,J,K)+.5*SMOTH*(VF(I,J,K)
     1  +SB(I,J,K)-2.*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)
      S(I,J,K)=VF(I,J,K)
  355 CONTINUE
      CALL DENS(S,T,RHO)
C
  360 CONTINUE
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,IPRINT).NE.0) GO TO 7000
 9001 CONTINUE
        IF(IINT.GE.ISWTCH) IPRINT=IPRTD2*24*3600/INT(DTI)
        WRITE(6,'(/,'' TIME ='',F10.2,''    IINT ='',I8,
     1  ''   IEXT ='',I8,''   IPRINT ='',I5,//)') TIME,IINT,IEXT,IPRINT
c       CALL PRXY(' ADVUA     ',TIME,ADVUA,IP,ISK,JM,JSK,0.0)
c       CALL PRXY(' ADVVA     ',TIME,ADVVA,IP,ISK,JM,1,0.0)
c       CALL PRXY(' ADX2D     ',TIME,ADX2D,IP,ISK,JM,1,0.0)
c       CALL PRXY(' ADY2D     ',TIME,ADY2D,IP,ISK,JM,1,0.0)
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,JSK,0.E0)
c       CALL PRXY(' WVBOT     ',TIME,WVBOT,IP,ISK,JM,JSK,0.E0)
        CALL PRXY(' DRX2D      ',TIME,DRX2D,IP,ISK,JM,JSK,0.E0)
        CALL PRXY(' DRY2D      ',TIME,DRY2D,IP,ISK,JM,JSK,0.E0)
        CALL PRXY('  FREE SURFACE ',TIME, ELB,IP,ISK,JM,JSK,0.E0)
        CALL PRXY(' AVERAGE U COMP.',TIME,UAB,IP,ISK,JM,JSK,0.E0)
        CALL PRXY(' AVERAGE V COMP.',TIME,VAB,IP,ISK,JM,JSK,0.E0)
C       CALL FINDPSI
        IF(MODE.NE.2) THEN
c       CALL PRXZ('DRHOX' ,TIME,DRHOX,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c       CALL PRXZ('DRHOY' ,TIME,DRHOY,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c         CALL PRXYZ('AAM       ',TIME,AAM,IP,ISK,JM,1,KB,0.E0)
          CALL PRXYZ(' T ',TIME,T  ,IP,ISK,JM,JSK,KB,0.E0)
          CALL PRXYZ(' S ',TIME,S  ,IP,ISK,JM,JSK,KB,0.E0)
C         CALL PRXYZ(' RHO ',TIME,RHO,IP,ISK,JM,JSK,KB,0.E0)
          CALL PRXZ('TB' ,TIME,TB,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
          CALL PRXZ('SB' ,TIME,SB,IM,ISK,JM,10,25,KB,1.,DT,ZZ)
          CALL PRXZ('RHO',TIME,RHO,IM,ISK,JM,10,25,KB,0.0,DT,ZZ)
          CALL PRXZ('RMEAN' ,TIME,RMEAN,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c         CALL PRXZ('Q2' ,TIME,Q2,IM,ISK,JM,10,25,KB,0.0,DT,Z)
          CALL PRXZ('KM' ,TIME,KM,IM,ISK,JM,10,25,KB,0.0,DT,Z)
          CALL PRXZ('U' ,TIME,U,IM,ISK,JM,10,25,KB,0.0,DT,ZZ)
          CALL PRXZ('V' ,TIME,V,IM,ISK,JM,10,25,KB,0.0,DT,ZZ)
        ENDIF
C
      VTOT=0.
      ATOT=0.
      TAVER=0.
      SAVER=0.
      EAVER=0.
      DO 8888 K=1,KBM1
      DO 8888 J=1,JM
      DO 8888 I=1,IM
      DAREA=DX(I,J)*DY(I,J)*FSM(I,J)
      DVOL =DAREA          *DT(I,J)*DZ(K)          
      VTOT=VTOT+DVOL    
      TAVER=TAVER+TB(I,J,K)*DVOL    
      SAVER=SAVER+SB(I,J,K)*DVOL   
 8888 CONTINUE
      DO 8889 J=1,JM
      DO 8889 I=1,IM
      DAREA=DX(I,J)*DY(I,J)*FSM(I,J)
      ATOT=ATOT+DAREA
      EAVER=EAVER+ET(I,J)*DAREA
 8889 CONTINUE
      TAVER=TAVER/VTOT
      SAVER=SAVER/VTOT
      EAVER=EAVER/ATOT
      WRITE(6,'(''  VTOT ='',E20.8,''      ATOT ='',E20.8,
     1 ''      EAVER ='',E20.8 )') VTOT,ATOT,EAVER            
      WRITE(6,'(''  TAVER ='',E20.8,
     1      ''     SAVER ='',E20.8/)') TAVER,SAVER
      IF(VAMAX.GT.VMAXL) 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,JSK,0.E0)
        CALL PRXY('  UAF           ',TIME,UAF,IP,ISK,JM,JSK,0.E0)
        CALL PRXY('  VA            ',TIME,VA ,IP,ISK,JM,JSK,0.E0)
        CALL PRXY('  VAF           ',TIME,VAF,IP,ISK,JM,JSK,0.E0)
      WRITE(6,'(//////////////////)')
      WRITE(6,'(''************************************************'')')
      WRITE(6,'(''************ ABNORMAL JOB END ******************'')')
      WRITE(6,'(''************* USER TERMINATED ******************'')')
      WRITE(6,'(''************************************************'')')
      WRITE(6,'('' VAMAX ='',E12.3,''   IMAX,JMAX ='',2I5)')
     1          VAMAX,IMAX,JMAX
      WRITE(6,'('' IINT,IEXT ='',2I6)') IINT,IEXT
      STOP
      ENDIF
 7000 CONTINUE
c     IF(MOD(IINT,IPRINT).EQ.0)
c    1  WRITE(40) TIME,Z,ZZ,X,Y,DX,DY,H,FSM,U,V,UA,VA,T,S,EL
C
C--------------------------------------------------------------------
C             END PRINT SECTION
C--------------------------------------------------------------------
 9000                     CONTINUE
C***********************************************************************
C                                                                      *
C                  END NUMERICAL INTEGRATION                           *
C                                                                      *
C***********************************************************************
        CALL EPRXYZ('W VEL. COMP.',TIME,W,IM,ISK,JM,JSK,KB,.00001E0)
c       CALL EPRXYZ('Q2         ',TIME,Q2,IM,ISK,JM,JSK,KB,0.0)
C             SAVE THIS DATA FOR A SEAMLESS RESTART
C     WRITE(140) TIME,
C    1  WUBOT,WVBOT,AAM2D,UA,UAB,VA,VAB,EL,ELB,ET,ETB,EGB,
C    2  UTB,VTB,U,UB,W,V,VB,T,TB,S,SB,RHO,ADX2D,ADY2D,ADVUA,ADVVA,
C    3  KM,KH,KQ,L,Q2,Q2B,AAM,Q2L,Q2LB
C
      STOP
      END
C
      SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE)
       INCLUDE 'comblk97.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))
  120 CONTINUE  
      DO 130 J=3,JMM2
      DO 130 I=3,IMM2
      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)) )
  130 CONTINUE  
 5000 CONTINUE
C 
      RETURN
      END
C               
      SUBROUTINE ADVQ(QB,Q,DTI2,QF)
       INCLUDE 'comblk97.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 ADVT(FB,F,FCLIM,DTI2,FF)
C
C     THIS SUBROUTINE INTEGRATES CONSERVATIVE SCALAR EQUATIONS
C
       INCLUDE 'comblk97.h'   
      DIMENSION FB(IM,JM,KB),F(IM,JM,KB),FF(IM,JM,KB),FCLIM(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******  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)-FCLIM(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))*TPRNI
     2    *(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J)/(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))*TPRNI
     2    *(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J)/(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)+FCLIM(I,J,K)
C
C****** DO VERTICAL ADVECTION *************************************
      DO 505 J=2,JMM1
      DO 505 I=2,IMM1
 505  FF(I,J,1)=-.5*(F(I,J,1)+F(I,J,2))*W(I,J,2)*ART(I,J)/DZ(1)
      DO 520 K=2,KBM1
      DO 520 J=2,JMM1
      DO 520 I=2,IMM1
 520  FF(I,J,K)=.5*((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)/DZ(K)
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 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 'comblk97.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
      DO 59 K=1,KB
      DO 59 J=1,JM
      DO 59 I=1,IM
      CURV(I,J,K)=0.
      ADVX(I,J,K)=0.
      XFLUX(I,J,K)=0.
  59  YFLUX(I,J,K)=0.
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-----------------------------------------------------------------
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
      ADVX(I,J,K)= 
     1            +XFLUX(I,J,K)-XFLUX(I-1,J,K)
     2            +YFLUX(I,J+1,K)-YFLUX(I,J,K)
  146 CONTINUE
      DO 150 K=1,KBM1
      DO 150 J=3,JMM2
      DO 150 I=3,IMM2
      ADVX(I,J,K)=ADVX(I,J,K)
     1   -ARU(I,J)*.25*(CURV(I,J,K)*DT(I,J)*(V(I,J+1,K)+V(I,J,K))
     2             +CURV(I-1,J,K)*DT(I-1,J)*(V(I-1,J+1,K)+V(I-1,J,K)))
  150 CONTINUE
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
      ADVY(I,J,K)=
     1            +XFLUX(I+1,J,K)-XFLUX(I,J,K)
     2            +YFLUX(I,J,K)-YFLUX(I,J-1,K)
 400  CONTINUE
      DO 410 K=1,KBM1
      DO 410 J=2,JMM1
      DO 410 I=2,IMM1
      ADVY(I,J,K)=ADVY(I,J,K)
     1   +ARV(I,J)*.25*(CURV(I,J,K)*DT(I,J)*(U(I+1,J,K)+U(I,J,K))
     2             +CURV(I,J-1,K)*DT(I,J-1)*(U(I+1,J-1,K)+U(I,J-1,K)))
 410  CONTINUE
C
      RETURN
      END
C               
      SUBROUTINE ADVU(DRHOX,ADVX,DTI2)
       INCLUDE 'comblk97.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)+(UF(I,J,K)-UF(I,J,K+1))*ARU(I,J)/DZ(K)
     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 'comblk97.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)+(VF(I,J,K)-VF(I,J,K+1))*ARV(I,J)/DZ(K)
     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)
       INCLUDE 'comblk97.h'   
      DIMENSION DRHOX(IM,JM,KB),DRHOY(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,JMM1
      DO 300 I=2,IMM1
  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,JMM1
      DO 310 I=2,IMM1
 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,JMM1
      DO 360 I=2,IMM1
 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,JMM1
      DO 500 I=2,IMM1
 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,JMM1
      DO 510 I=2,IMM1
 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,JMM1
      DO 560 I=2,IMM1
 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,JMM1
      DO 561 I=2,IMM1
      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 'comblk97.h'   
C
C  Closed boundary conditions are automatically enabled through
C  specification of the masks, DUM, DVM and FSM in which case open
C  boundary condition, included below, will be overwritten.
C
C
C*******************************************************************************
C                            POM (C-Grid)   
C                          ================
C
C    The diagram below and some of the text was provided by D.-S. Ko. 
C    It is for the case where U and V are the primary boundary conditions
C    together with T and S (co-located with E). E = EL itself is rather
C    unimportant and is substituted from an adjacent interior point.
C    
C    Inside ....... indicate the interior (non-boundary) grid points.
C    In general only those variables in the interior are computed and
C    variables at open boundary have to be specified.
C    All interpolations are centered in space except those at lateral
C    open boundary where an upstream scheme is usually used.
C
C    Horizontal locations of E(EL), T and S (etc.) are coincident.
C    NU = Not Used is indicated by *. However, for attractive output
C    adjacent interior values may be filled in at these points.   
C
C    I have been asked many times what kind of B.C. along the side wall
C    POM uses from people not acquainted with sigma coordinates. Although
C    the issue is not important as it might be for z-level grids, a direct
C    answer is "half slip" which, of course, is between free slip and
C    non-slip B.C.s.
C
C-------------------------------------------------------------------------------
C     |                               N O R T H
C     |
C     |    1     2     3           I-1   I    I+1         IM-2  IM-1   IM 
C-----+-------------------------------------------------------------------------
C     |   NU BC BC                                                    BC BC
C     |    v  v  v                                                     v  v
C     |
C     |BC> U* E  U  E  U  E  .  .  U  E  U  E  U  E  .  .  U  E  U  E  U  E  <BC
C     |    |     |     |           |     |     |           |     |     |      
C  JM |BC> +--V--+--V--+--V--   .  +--V--+--V--+--V--   .  +--V--+--V--+--V- <BC
C     |    |     | ....|...........|.....|.....|...........|.....|.... |      
C     |    U* E  U :E  U  E  .  .  U  E  U  E  U  E  .  .  U  E  U  E: U  E 
C     |    |     | :   |           |     |     |           |     |   : |      
C JM-1|    +--V--+--V--+--V--   .  +--V--+--V--+--V--   .  +--V--+--V--+--V-
C     |    |     | :   |           |     |     |           |     |   : |
C     |    U* E  U :E  U  E  .  .  U  E  U  E  U  E  .  .  U  E  U  E: U  E
C     |    |     | :   |           |     |     |           |     |   : |
C JM-2|    +--V--+--V--+--V--   .  +--V--+--V--+--V--   .  +--V--+--V--+--V-
C     |            :                                                 :          
C W   |       .  . :.  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .: .  .    E
C E   |            :                    Interior                     :         A
C S   |       .    :.     .     .     .     .     .     .     .     .:    .    S
C T   |            :                                                 :         T
C     |    |     | :   |           |     |     |           |     |   : |
C     |    U* E  U :E  U  E  .  .  U  E  U  E  U  E  .  .  U  E  U  E: U  E
C     |    |     | :   |           |     |     |           |     |   : |
C   3 |    +--V--+--V--+--V--   .  +--V--+--V--+--V--   .  +--V--+--V--+--V-
C     |    |     | :   |           |     |     |           |     |   : |      
C     |    U* E  U :E  U  E  .  .  U  E  U  E  U  E  .  .  U  E  U  E: U  E
C     |    |     | ....|...........|.....|.....|...........|.....|.... |      
C   2 |BC> +--V--+--V--+--V--   .  +--V--+--V--+--V--   .  +--V--+--V--+--V- <BC
C     |    |     |     |           |     |     |           |     |     |      
C     |BC> U* E  U  E  U  E  .  .  U  E  U  E  U  E  .  .  U  E  U  E  U  E  <BC
C     |    |     |     |           |     |     |           |     |     |      --
C   1 |NU> +--V*-+--V*-+--V*-   .  +--V*-+--V*-+--V*-   .  +--V*-+--V*-+--V* <NU
C     |
C     |    ^  ^  ^                                                     ^  ^
C     |   NU BC BC                                                    BC BC
C-----+-------------------------------------------------------------------------
C     |    1     2     3           I-1   I     I+1         IM-2  IM-1  IM
C     |
C     |                                S O U T H
C-------------------------------------------------------------------------------
C
C
C
      DATA PI/3.141592654/,GEE/9.807/
C
      GO TO (10,20,30,40,50,60), IDX
C
C-----------------------------------------------------------------------
C                   EXTERNAL B.C.'S
C  In this example the governing boundary conditions are a radiation   
C  condition on UAF in the east and in the west and VAF in the north
C  and south. The tangential velocities are set to zero on both boundaries.
C  These are only one set of possibilities and may not represent a choice
C  which yields the most physically realistic result.
C-----------------------------------------------------------------------
C
C
 10   CONTINUE
C---------ELEVATION --------------------
C  In this application, elevation is not a primary B.C.     
      DO 120 J=1,JM
      ELF(1,J)=ELF(2,J)
      ELF(IM,J)=ELF(IMM1,J)
120   CONTINUE
      DO 130 I=1,IM
      ELF(I,1)=ELF(I,2)
      ELF(I,JM)=ELF(I,JMM1)
130   CONTINUE
C
      DO 140 J=1,JM
      DO 140 I=1,IM
 140  ELF(I,J)=ELF(I,J)*FSM(I,J)
      RETURN
C
C
 20   CONTINUE
C---------- VELOCITY --------------        
C--- Governing external B.C.s -------------
      DO 210 J=2,JMM1
C            EAST
      UAF(IM,J)=UABE(J)+RFE*SQRT(GRAV/H(IMM1,J))*(EL(IMM1,J)-ELE(J))
      UAF(IM,J)=RAMP*UAF(IM,J)
      VAF(IM,J)=0.0
C            WEST
      UAF(2,J)=UABW(J)-RFW*SQRT(GRAV/H(2,J))*(EL(2,J)-ELW(J))
      UAF(2,J)=RAMP*UAF(2,J)
      UAF(1,J)=UAF(2,J)
      VAF(1,J)=0.0
 210  CONTINUE  
      DO 220 I=2,IMM1
C            NORTH
      VAF(I,JM)=VABN(I)+RFN*SQRT(GRAV/H(I,JMM1))*(EL(I,JMM1)-ELN(I))
      VAF(I,JM)=RAMP*VAF(I,JM)
      UAF(I,JM)=0.0
C            SOUTH
      VAF(I,2)=VABS(I)-RFS*SQRT(GRAV/H(I,2))*(EL(I,2)-ELS(I))
      VAF(I,2)=RAMP*VAF(I,2)
      VAF(I,1)=VAF(I,2)
      UAF(I,1)=0.0
 220  CONTINUE
C---------------------------
C
      DO 240 J=1,JM
      DO 240 I=1,IM
      UAF(I,J)=UAF(I,J)*DUM(I,J)
 240  VAF(I,J)=VAF(I,J)*DVM(I,J)
      RETURN
C
C
 30   CONTINUE
C-----------------------------------------------------------------------
C                   INTERNAL VEL B.C.'S
C-----------------------------------------------------------------------
C         Radiation boundary conditions. Smoothing is used
C         in the direction tangential to the boundaries.
C
      HMAX=4500.
      DO 340 K=1,KBM1
      DO 340 J=2,JMM1
C            EAST
      GA=SQRT(H(IM,J)/HMAX)    
      UF(IM,J,K)
     1  =GA*(.25*U(IMM1,J-1,K)+.5*U(IMM1,J,K)+.25*U(IMM1,J+1,K))
     2   +(1.E0-GA)*(.25*U(IM,J-1,K)+.5*U(IM,J,K)+.25*U(IM,J+1,K))
      VF(IM,J,K)=0.                                        
C            WEST
      GA=SQRT(H(1,J)/HMAX)    
      UF(2,J,K)
     1  =GA*(.25*U(3,J-1,K)+.5*U(3,J,K)+.25*U(3,J+1,K))
     2   +(1.-GA)*(.25*U(2,J-1,K)+.5*U(2,J,K)+.25*U(2,J+1,K))
      UF(1,J,K)=UF(2,J,K)
      VF(1,J,K)=0.
  340 CONTINUE
      DO 360 K=1,KBM1
      DO 360 I=2,IMM1
C            NORTH
      GA=SQRT(H(I,JM)/HMAX)    
      VF(I,JM,K)
     1  =GA*(.25*V(I-1,JMM1,K)+.5*V(I,JMM1,K)+.25*V(I+1,JMM1,K))
     2   +(1.-GA)*(.25*V(I-1,JM,K)+.5*V(I,JM,K)+.25*V(I+1,JM,K))
      UF(I,JM,K)=0.
C            SOUTH
      GA=SQRT(H(I,1)/HMAX)    
      VF(I,2,K)
     1  =GA*(.25*V(I-1,3,K)+.5*V(I,3,K)+.25*V(I+1,3,K))
     2   +(1.-GA)*(.25*V(I-1,2,K)+.5*V(I,2,K)+.25*V(I+1,2,K))
      VF(I,1,K)=VF(I,2,K)
      UF(I,1,K)=0.
  360 CONTINUE
C
      DO 380 K=1,KBM1
      DO 380 J=1,JM
      DO 380 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)
  380 CONTINUE
C      
      RETURN
C
C
 40   CONTINUE
C-----------------------------------------------------------------------
C                   TEMP(UF) & SAL(VF) B.C.'S
C-----------------------------------------------------------------------
C
      DO 420 K=1,KBM1
      DO 420 J=1,JM
C                     EAST
      U1=2.*U(IM,J,K)*DTI/(DX(IM,J)+DX(IMM1,J))
      IF(U1.LE.0.) THEN
        UF(IM,J,K)=T(IM,J,K)-U1*(TBE(J,K) -T(IM,J,K)) 
        VF(IM,J,K)=S(IM,J,K)-U1*(SBE(J,K) -S(IM,J,K)) 
      ELSE
        UF(IM,J,K)=T(IM,J,K)-U1*(T(IM,J,K)-T(IMM1,J,K))
        VF(IM,J,K)=S(IM,J,K)-U1*(S(IM,J,K)-S(IMM1,J,K))
        IF(K.NE.1.AND.K.NE.KBM1) THEN
          WM=.5*(W(IMM1,J,K)+W(IMM1,J,K+1))
     1               *DTI/((ZZ(K-1)-ZZ(K+1))*DT(IMM1,J))
          UF(IM,J,K)=UF(IM,J,K)-WM*(T(IMM1,J,K-1)-T(IMM1,J,K+1))
          VF(IM,J,K)=VF(IM,J,K)-WM*(S(IMM1,J,K-1)-S(IMM1,J,K+1))
        ENDIF 
      ENDIF
C                      WEST
      U1=2.*U(2,J,K)*DTI/(DX(1,J)+DX(2,J))
      IF(U1.GE.0.) THEN
        UF(1,J,K)=T(1,J,K)-U1*(T(1,J,K)-TBW(J,K)) 
        VF(1,J,K)=S(1,J,K)-U1*(S(1,J,K)-SBW(J,K)) 
      ELSE
        UF(1,J,K)=T(1,J,K)-U1*(T(2,J,K)-T(1,J,K))
        VF(1,J,K)=S(1,J,K)-U1*(S(2,J,K)-S(1,J,K))
        IF(K.NE.1.AND.K.NE.KBM1) THEN
          WM=.5*(W(2,J,K)+W(2,J,K+1))*DTI/((ZZ(K-1)-ZZ(K+1))*DT(2,J))
          UF(1,J,K)=UF(1,J,K)-WM*(T(2,J,K-1)-T(2,J,K+1))
          VF(1,J,K)=VF(1,J,K)-WM*(S(2,J,K-1)-S(2,J,K+1))
        ENDIF
      ENDIF
  420 CONTINUE
      DO 440 K=1,KBM1
      DO 440 I=1,IM
C                     NORTH
      U1=2.*V(I,JM,K)*DTI/(DY(I,JM)+DY(I,JMM1))
      IF(U1.LE.0.) THEN
        UF(I,JM,K)=T(I,JM,K)-U1*(TBN(I,K)-T(I,JM,K)) 
        VF(I,JM,K)=S(I,JM,K)-U1*(SBN(I,K)-S(I,JM,K)) 
      ELSE
        UF(I,JM,K)=T(I,JM,K)-U1*(T(I,JM,K)-T(I,JMM1,K))
        VF(I,JM,K)=S(I,JM,K)-U1*(S(I,JM,K)-S(I,JMM1,K))
        IF(K.NE.1.AND.K.NE.KBM1) THEN
          WM=.5*(W(I,JMM1,K)+W(I,JMM1,K+1))
     1                *DTI/((ZZ(K-1)-ZZ(K+1))*DT(I,JMM1))
          UF(I,JM,K)=UF(I,JM,K)-WM*(T(I,JMM1,K-1)-T(I,JMM1,K+1))
          VF(I,JM,K)=VF(I,JM,K)-WM*(S(I,JMM1,K-1)-S(I,JMM1,K+1))
        ENDIF
      ENDIF
C                    SOUTH
      U1=2.*V(I,2,K)*DTI/(DY(I,1)+DY(I,2))
      IF(U1.GE.0.) THEN
        UF(I,1,K)=T(I,1,K)-U1*(T(I,1,K)-TBS(I,K)) 
        VF(I,1,K)=S(I,1,K)-U1*(S(I,1,K)-SBS(I,K)) 
      ELSE
        UF(I,1,K)=T(I,1,K)-U1*(T(I,2,K)-T(I,1,K))
        VF(I,1,K)=S(I,1,K)-U1*(S(I,2,K)-S(I,1,K))
        IF(K.NE.1.AND.K.NE.KBM1) THEN
          WM=.5*(W(I,2,K)+W(I,2,K+1))*DTI/((ZZ(K-1)-ZZ(K+1))*DT(I,2))
          UF(I,1,K)=UF(I,1,K)-WM*(T(I,2,K-1)-T(I,2,K+1))
          VF(I,1,K)=VF(I,1,K)-WM*(S(I,2,K-1)-S(I,2,K+1))
        ENDIF
      ENDIF
  440 CONTINUE
C
      DO 460 K=1,KBM1
      DO 460 J=1,JM
      DO 460 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)
 460  CONTINUE
      RETURN
C
C
 50   CONTINUE
C---------------VERTICAL VEL. B. C.'S --------------------------------
      DO 550 K=1,KBM1
      DO 550 J=1,JM
      DO 550 I=1,IM
      W(I,J,K)=W(I,J,K)*FSM(I,J)
 550  CONTINUE
C    
      RETURN
C
C
 60   CONTINUE
C---------------- Q2 AND Q2L B.C.'S -----------------------------------
      DO 650 K=1,KB
C     
      DO 605 J=1,JM
      UF(1,J,K)=1.E-10
      VF(1,J,K)=1.E-10
      UF(IM,J,K)=1.E-10
      VF(IM,J,K)=1.E-10
 605  CONTINUE
      DO 610 I=1,IM
      UF(I,1,K)=1.E-10
      VF(I,1,K)=1.E-10
      UF(I,JM,K)=1.E-10
      VF(I,JM,K)=1.E-10
 610  CONTINUE
      DO 697 J=1,JM
      DO 697 I=1,IM
      UF(I,J,K)=UF(I,J,K)*FSM(I,J)
 697  VF(I,J,K)=VF(I,J,K)*FSM(I,J)
 650  CONTINUE
      RETURN
      END
C
      SUBROUTINE BCONDorl(IDX)
C
C  Herewith is an optional SUBROUTINE BCOND using Orlanski's
C  scheme (J. Comp. Phys. 21, 251-269 (1976)specialized for
C  the seamount problem. To make it work
C  for the seamount problem, I have had to add an extra condition
C  on an if statement in the T & S open boundary conditions 
C  which involves the sign of the normal velocity. But,the results 
C  are improved for the seamount problem if the following statement,   
C     if(cl.eq.0..and.ubw(j,k).ge.0.) uf(1,j,k)=tbw(j,k),
C  plus 3 others of the same kind are commented out. However, then
C  there is no connection to T & S open boundary values which would
C  be disaster for many other problems.
C
       INCLUDE 'comblk97.h'   
C
      DATA PI/3.14167/,GEE/9.807/
C
      GO TO (10,20,30,40,50,60), IDX
C
C-----------------------------------------------------------------------
C                   EXTERNAL B.C.'S
C  In this example the governing boundary conditions are a radition   
C  condition on UAF(IM,J) in the east and an inflow UAF(2,J) in the
C  west. The tangential velocities are set to zero on both boundaries.
C  These are only one set of possibilities and may not represent a choice
C  which yields the most physically realistic result.
C-----------------------------------------------------------------------
C
C
 10   CONTINUE
C---------ELEVATION --------------------
C  In this application, elevation is not a primary B.C.     
      DO 110 J=1,JM
      ELF(1,J)=ELF(2,J)
      ELF(IM,J)=ELF(IMM1,J)
110   CONTINUE
C
      DO 111 J=1,JM
      DO 111 I=1,IM
 111  ELF(I,J)=ELF(I,J)*FSM(I,J)
      RETURN
C
C
 20   CONTINUE
C---------- VELOCITY --------------        
      DO 210 J=2,JMM1
C--- Governing B.C.s -------------
C            WEST
      UAF(2,J)=RAMP*UABW(J)-SQRT(GRAV/H(2,J))*(EL(2,J)-ELW(J)) 
      UAF(1,J)=UAF(2,J)
      VAF(1,J)=0.0
C            EAST
      UAF(IM,J)=RAMP*UABE(J)+SQRT(GRAV/H(IMM1,J))*(EL(IMM1,J)-ELE(J)) 
      VAF(IM,J)=0.0
C---------------------------
 210  CONTINUE  
C
      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
C
C
 30   CONTINUE
C-----------------------------------------------------------------------
C                   INTERNAL VEL B.C.'S
C-----------------------------------------------------------------------
C         Eastern and western radiation boundary conditions according
C         to Orlanski's explicit scheme. 
C
      DO 140 K=1,KBM1
      DO 140 J=2,JMM1
c
c                 WEST 
c
      denom=(uf(3,j,k)+ub(3,j,k)-2.*u(4,j,k))
      if(denom.eq.0.)denom=0.01
      cl=(ub(3,j,k)-uf(3,j,k))/denom
      if(cl.gt.1.)cl=1.
      if(cl.lt.0.) then
         cl=0.
      endif
      uf(2,j,k)=(ub(2,j,k)*(1.-cl)+2.*cl*u(3,j,k))/(1.+cl)
      uf(1,j,k)=uf(2,j,k)
      vf(1,j,k)=0.
C
C                 EAST
C
      denom=(uf(im-1,j,k)+ub(im-1,j,k)-2.*u(im-2,j,k))
      if(denom.eq.0.)denom=0.01
      cl=(ub(im-1,j,k)-uf(im-1,j,k))/denom
      if(cl.gt.1.)cl=1.
      if(cl.lt.0.) then
         cl=0.
      endif
      uf(im,j,k)=(ub(im,j,k)*(1.-cl)+2.*cl*u(im-1,j,k))/(1.+cl)
      vf(im,j,k)=0.
  140 CONTINUE
C
C
      DO 160 K=1,KBM1
      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
C
 40   CONTINUE
C-----------------------------------------------------------------------
C                   TEMP(UF) & SAL(VF) B.C.'S
C-----------------------------------------------------------------------
      DO 220 K=1,KBM1
      DO 220 J=1,JM
C
C                 WEST 
C
      UBW(J,K)=UB(2,J,K)
      denom=(uf(2,j,k)+tb(2,j,k)-2.*t(3,j,k))
      if(denom.eq.0.)denom=0.01
      cl=(tb(2,j,k)-uf(2,j,k))/denom
      if(cl.gt.1.) cl=1.
      if(cl.lt.0.) then
         cl=0.
      endif
         uf(1,j,k)=(tb(1,j,k)*(1.-cl)+2.*cl*t(2,j,k))/(1.+cl)
      if(cl.eq.0..and.ubw(j,k).ge.0.) uf(1,j,k)=tbw(j,k)
C
      denom=(vf(2,j,k)+sb(2,j,k)-2.*s(3,j,k))
      if(denom.eq.0.)denom=0.01
      cl=(sb(2,j,k)-vf(2,j,k))/denom
      if(cl.gt.1.) cl=1.
      if(cl.lt.0.) then
         cl=0.
      endif
         vf(1,j,k)=(sb(1,j,k)*(1.-cl)+2.*cl*s(2,j,k))/(1.+cl)
      if(cl.eq.0..and.ubw(j,k).ge.0.) vf(1,j,k)=sbw(j,k)
C
C                 EAST
C
      UBE(J,K)=UB(IM,J,K)
      denom=(uf(im-1,j,k)+tb(im-1,j,k)-2.*t(im-2,j,k))
      if(denom.eq.0.)denom=0.01
      cl=(tb(im-1,j,k)-uf(im-1,j,k))/denom
      if(cl.gt.1.) cl=1.
      if(cl.lt.0.) then
         cl=0.
      endif
         uf(im,j,k)=(tb(im,j,k)*(1.-cl)+2.*cl*t(im-1,j,k))/(1.+cl)
      if(cl.eq.0..and.ube(j,k).le.0.) uf(im,j,k)=tbe(j,k)
C
      denom=(vf(im-1,j,k)+sb(im-1,j,k)-2.*s(im-2,j,k))
      if(denom.eq.0.)denom=0.01
      cl=(sb(im-1,j,k)-vf(im-1,j,k))/denom
      if(cl.gt.1.) cl=1.
      if(cl.lt.0.) then
         cl=0.
      endif
         vf(im,j,k)=(sb(im,j,k)*(1.-cl)+2.*cl*s(im-1,j,k))/(1.+cl)
      if(cl.eq.0..and.ube(j,k).le.0.) vf(im,j,k)=sbe(j,k)
C
  220 CONTINUE
C
      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
C    
      RETURN
C
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
      VF(IM,J,K)=1.E-10
      UF(1,J,K)=1.E-10
      VF(1,J,K)=1.E-10
 295  CONTINUE
      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
      END
C
      SUBROUTINE DENS(SI,TI,RHOO)
       INCLUDE 'comblk97.h'
      DIMENSION SI(IM,JM,KB),TI(IM,JM,KB),RHOO(IM,JM,KB)
C       If using 32 bit precision, it is recommended that
C       TR, SR, P, RHOR , CR be made double precision.
C       and the E's in the constants should be changed 
C       to D's.
C
C         THIS SUBROUTINE COMPUTES DENSITY- 1.000
C         T = POTENTIAL TEMPERATURE
C    ( See: Mellor, 1991, J. Atmos. Oceanic Tech., 609-611)
C
      DO 1 K=1,KBM1
      DO 1 J=1,JM
      DO 1 I=1,IM
      TR=TI(I,J,K)+TBIAS
      SR=SI(I,J,K)+SBIAS
      TR2=TR*TR
      TR3=TR2*TR
      TR4=TR3*TR
C            Approximate pressure in units of bars
      P=-GRAV*1.025*ZZ(K)*DT(I,J)*0.01
C
      RHOR = -0.157406  + 6.793952E-2*TR
     $        - 9.095290E-3*TR2 + 1.001685E-4*TR3
     $        - 1.120083E-6*TR4 + 6.536332E-9*TR4*TR
C
      RHOR = RHOR + (0.824493 - 4.0899E-3*TR
     $       + 7.6438E-5*TR2 - 8.2467E-7*TR3
     $       + 5.3875E-9*TR4) * SR
     $       + (-5.72466E-3 + 1.0227E-4*TR
     $       - 1.6546E-6*TR2) * ABS(SR)**1.5
     $       + 4.8314E-4 * SR*SR
C
      CR=1449.2+.0821*P+15.E-7*P**2+4.55*TR-.045*TR2
     $                              +1.34*(SR-35.)
      RHOR=RHOR + 1.E5*P/(CR*CR)*(1.-2.*P/(CR*CR))
C
      RHOO(I,J,K)=RHOR*1.E-3*FSM(I,J)
    1 CONTINUE
C
      RETURN
      END
C               
      SUBROUTINE DEPTH(Z,ZZ,DZ,DZZ,KB)
      DIMENSION Z(KB),ZZ(KB),DZ(KB),DZZ(KB)
      DIMENSION KDZ(12)
      DATA KDZ/1,1,2,4,8,16,32,64,128,256,512,1024/
C ***
C    THIS SUBROUTINE ESTABLISHES THE VERTICAL SIGMA GRID WITH LOG
C    DISTRIBUTIONS  AT THE TOP AND BOTTOM AND A LINEAR DISTRIBUTION
C    BETWEEN. CHOOSE KL1 AND KL2. THE VALUES, KL1 = 2 AND KL2 = KB-2,
C    YIELD AN EQUALLY SPACED DISTRIBUTION WITH NO LOG PORTIONS.   
C ***
      KL1=6
      KL2=KB-2
C
      Z(1)=0.
      DO K=2,KL1
      Z(K)=Z(K-1)+KDZ(K-1)
      ENDDO
      DELZ=Z(KL1)-Z(KL1-1)
      DO K=KL1+1,KL2
      Z(K)=Z(K-1)+DELZ
      ENDDO
      DO K=KL2+1,KB
      DZ(K)=FLOAT(KDZ(KB-K+1))*DELZ/FLOAT(KDZ(KB-KL2))
      Z(K)=Z(K-1)+DZ(K)                                         
      ENDDO
      DO K=1,KB
      Z(K)=-Z(K)/Z(KB)
      ENDDO
      DO K=1,KB-1
      ZZ(K)=0.5*(Z(K)+Z(K+1))
      ENDDO
      ZZ(KB)=2.*ZZ(KB-1)-ZZ(KB-2)
      DO 60 K=1,KB-1
      DZ(K)=Z(K)-Z(K+1)
   60 DZZ(K)=ZZ(K)-ZZ(K+1)
      DZ(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)
   90 CONTINUE
   80 FORMAT(' ',I5,4F10.3)
      RETURN
      END
C               
      SUBROUTINE EPRXYZ(LABEL,TIME,A,IM,ISKP,JM,JSKP,KB,DUMY)
      DIMENSION A(IM,JM,KB),NUM(300)
      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        WITH FLOATING POINT NUMBERS.
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
C  This subroutine calculates the stream function first assuming 
C  zero on the southern boundary and then, using the values on the 
C  the western boundary, the stream function is calculated again.
C  If the elevation field is near steady state, the two calculations
C  should agree; otherwise not.
         INCLUDE 'comblk97.h'    
        DO 9004 J=1,JM
        DO 9004 I=1,IM
 9004   PSI(I,J)=0.E0
C    Sweep northward.
        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,2,JM,2,1.E6)
C    Sweep eastward.
        DO 9005 J=2,JM
        DO 9005 I=2,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,2,JM,2,1.E6)
      RETURN
      END
C          
      SUBROUTINE PROFQ(DT2)
C
       INCLUDE 'comblk97.h'    
      REAL KN,KAPPA        
      DIMENSION GH(IM,JM,KB),SM(IM,JM,KB),SH(IM,JM,KB)
      DIMENSION PROD(IM,JM,KB),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 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----- Calculate speed of sound squared ----------------------------
      DO 101 K=1,KBM1
      DO 101 J=1,JM
      DO 101 I=1,IM
      TP=T(I,J,K)+TBIAS
      SP=S(I,J,K)+SBIAS
C      Calculate pressure in units of decibars
      P=-GRAV*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.0)  
      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
C----- Calculate buoyancy gradient ---------------------------------
      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)=GRAV*((RHO(I,J,K-1)-RHO(I,J,K))/(DZZ(K-1)*DH(I,J)))
     1 +GRAV**2*2.*1.025/(CC(I,J,K-1)**2+CC(I,J,K)**2)
  102 CONTINUE
C------ CALC. T.K.E. PRODUCTION -----------------------------------
      DO 120 K=2,KBM1
      DO 120 J=2,JMM1
      DO 120 I=1,IMM1
      PROD(I,J,K)=KM(I,J,K)*.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
  120 PROD(I,J,K)=PROD(I,J,K)+KH(I,J,K)*BOYGR(I,J,K)
      DO 110 K=2,KBM1
      DO 110 J=1,JM
      DO 110 I=1,IM
      STF=1.
C     GHC=-2.5
C     IF(GH(I,J,K).LT.0.) STF=1.0-0.9*(GH(I,J,K)/GHC)**1.5
C     IF(GH(I,J,K).LT.GHC) STF=0.1
  110 DTEF(I,J,K)=Q2B(I,J,K)*SQRT(Q2B(I,J,K))/(B1*Q2LB(I,J,K)+SMALL)*STF
      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
      IF(UF(I,J,2).LT.SMALL) UF(I,J,2)=SMALL
      EE(I,J,2)=0.
      GG(I,J,2)=-KAPPA*Z(2)*DH(I,J)*UF(I,J,2)
      VF(I,J,KB)=0.                                         
  155 CONTINUE   
      DO 160 K=3,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,KB-2
      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).GT.SMALL.AND.VF(I,J,K).GT.SMALL) GOTO 180
          UF(I,J,K)=SMALL
          VF(I,J,K)=SMALL
  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*SH(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
  280 CONTINUE
      RETURN   
      END
C
C--------------------------------------------------------------------
      SUBROUTINE PROFT(F,WFSURF,SWRAD,FSURF,NBC,DT2)
C
	INCLUDE 'comblk97.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 'comblk97.h'    
      DIMENSION DH(IM,JM)
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,IMM1
      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 'comblk97.h'    
      DIMENSION DH(IM,JM)
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',1PE9.2)
      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
      PARAMETER (KE=4)
      DIMENSION A(IM,JM,KB),NUM(150),LINE(150),KP(KE)
      CHARACTER LABEL*(*)
      DATA KP /1,2,7,14/
      DATA ZERO /1.E-10/
      KP(3)=KB/2
      KP(4)=KB-1
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 SLPMIN(H,IM,JM,FSM,SL)
      DIMENSION H(IM,JM),FSM(IM,JM)
      DIMENSION SL(IM,JM)
C     This subroutine limits the maximum difference in H divided
C     by twice the average of H of adjacent cells.
C     The maximum possible value is unity.
C     
      SLMIN=0.2
C
      DO 100 LOOP=1,10
C    sweep right
      DO 3 J=2,JM-1
      DO 1 I=2,IM-1
      IF(FSM(I,J).EQ.0..OR.FSM(I+1,J).EQ.0.) GOTO 1
      SL(I,J)=ABS(H(I+1,J)-H(I,J))/(H(I,J)+H(I+1,J))
      IF(SL(I,J).LT.SLMIN) GOTO 1
      DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I+1,J))
      SN=-1.
      IF(H(I+1,J).GT.H(I,J)) SN=1.
      H(I+1,J)=H(I+1,J)-SN*DH
      H(I,J)=H(I,J)+SN*DH
   1  CONTINUE
C    sweep left
      DO 2 I=IM-1,2,-1
      IF(FSM(I,J).EQ.0..OR.FSM(I+1,J).EQ.0.) GOTO 2
      SL(I,J)=ABS(H(I+1,J)-H(I,J))/(H(I,J)+H(I+1,J))
      IF(SL(I,J).LT.SLMIN) GOTO 2
      DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I+1,J))
      SN=-1.
      IF(H(I+1,J).GT.H(I,J)) SN=1.
      H(I+1,J)=H(I+1,J)-SN*DH
      H(I,J)=H(I,J)+SN*DH
   2  CONTINUE
   3  CONTINUE
C   sweep up     
      DO 13 I=2,IM-1
      DO 11 J=2,JM-1
      IF(FSM(I,J).EQ.0..OR.FSM(I,J+1).EQ.0.) GOTO 11
      SL(I,J)=ABS(H(I,J+1)-H(I,J))/(H(I,J)+H(I,J+1))
      IF(SL(I,J).LT.SLMIN) GOTO 11
      DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I,J+1))
      SN=-1.
      IF(H(I,J+1).GT.H(I,J)) SN=1.
      H(I,J+1)=H(I,J+1)-SN*DH
      H(I,J)=H(I,J)+SN*DH
   11 CONTINUE
C   sweep down
      DO 12 J=JM-1,2,-1
      IF(FSM(I,J).EQ.0..OR.FSM(I,J+1).EQ.0.) GOTO 12 
      SL(I,J)=ABS(H(I,J+1)-H(I,J))/(H(I,J)+H(I,J+1))
      IF(SL(I,J).LT.SLMIN) GOTO 12
      DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I,J+1))
      SN=-1.
      IF(H(I,J+1).GT.H(I,J)) SN=1.
      H(I,J+1)=H(I,J+1)-SN*DH
      H(I,J)=H(I,J)+SN*DH
  12  CONTINUE
  13  CONTINUE
C
  100 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 = ',F9.3,' DAYS    MULTIPLY ALL VALUES BY ',1PE8.2)
C23456   |         |         |         |         |         |         |  xxxxxxx|
      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(8X,'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

      SUBROUTINE SEAMOUNT(TSURF,SSURF,TCLIM,SCLIM)
       INCLUDE 'comblk97.h'
      DIMENSION X(IM),Y(JM)
      DIMENSION TSURF(IM,JM),SSURF(IM,JM),
     1          TCLIM(IM,JM,KB),SCLIM(IM,JM,KB)
      DATA PI/3.141592654/,SMALL/1.E-10/

C     SET UP  GRID AND INITIAL CONDITIONS FOR STAND-ALONE TEST RUN
C For an island, DELH > 1.0; for a seamount DELH < 1.0.
      DELX=8.E3
      DELH=0.9
      RA=25.E3
      H0=4500
C Flow magnitude
      VEL=0.2
      CALL DEPTH(Z,ZZ,DZ,DZZ,KB)
      DO 12 J=1,JM
      DO 12 I=1,IM
C  Constant DX, DY.
C     DX(I,J)=DELX   
C     DY(I,J)=DELX   
C  Variable DX, DY.
      DX(I,J)=DELX-DELX/2.*SIN(PI*FLOAT(I)/FLOAT(IM))
      DY(I,J)=DELX-DELX/2.*SIN(PI*FLOAT(J)/FLOAT(JM))
      ART(I,J)=DX(I,J)*DY(I,J)
      COR(I,J)=1.E-4
  12  CONTINUE     
      X(1)=0.
      DO I=2,IM
      X(I)=X(I-1)+DX(I-1,1)
      ENDDO
      Y(1)=0.
      DO J=2,JM
      Y(J)=Y(J-1)+DY(1,J-1)  
      ENDDO
      DO I=1,IM
      DO J=1,JM
      H(I,J)= 4500.*(1.-DELH*EXP(-
     1  ((X(I)-X((IM+1)/2))**2+(Y(J)-Y((JM+1)/2))**2)/RA**2) )
      IF(H(I,J).LT.1.0) H(I,J)=1.0
      ENDDO
C  Close the north and south boundaries to form a channel.
      H(I,1)=1.
      H(I,JM)=1.
      ENDDO
      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)
      ARV(1,J)=ARV(2,J)
 121  CONTINUE        
      DO 122 I=1,IM
      ARU(I,1)=ARU(I,2)
      ARV(I,1)=ARV(I,2)
 122  CONTINUE       
C
      write(6,'('' got here 1'')')
C
      IP=IM
      ISK=2
      JSK=2
C     CALL PRXY('  TOPOGRAPHY 1 ',TIME, H,IP,ISK,JM,JSK,10.)
C
      DO J=1,JM
      DO I=1,IM
      DT(I,J)=H(I,J)
      FSM(I,J)=0.
      DUM(I,J)=0.
      DVM(I,J)=0.
      IF(H(I,J).GT.1.) FSM(I,J)=1.
      ENDDO
      ENDDO
      DO J=2,JM
      DO I=2,IM
      DUM(I,J)=FSM(I,J)*FSM(I-1,J)                                     
      DVM(I,J)=FSM(I,J)*FSM(I,J-1)                                     
      ENDDO
      ENDDO
      write(6,'('' got here 2'')')
C
C  Adjust bottom topography so that cell to cell variations
C  in H do not exceed parameter that is set in SLPMIN
C
C     CALL SLPMIN (H,IM,JM,FSM,TPS)
C
C---------------------------------------------------------------------
C  Set initial conditions for seamount problem.
C---------------------------------------------------------------------
      DO 15 K=1,KBM1
      DO 15 J=1,JM
      DO 15 I=1,IM
      TB(I,J,K)=5.+15.*EXP(ZZ(K)*H(I,J)/1000.)-TBIAS          
      SB(I,J,K)=35.0-SBIAS   
      TCLIM(I,J,K)=TB(I,J,K)
      SCLIM(I,J,K)=SB(I,J,K)
      UB(I,J,K)=VEL*DUM(I,J)  
      VB(I,J,K)=0.            
      AAM(I,J,K)=500.
  15  CONTINUE
C
      write(6,'('' got here 3'')')
      CALL DENS(SB,TB,RHO)
C
C  In this application, the initial condition for density is 
C  a function of z (cartesian). When this is not so, make sure
C  that RMEAN has been area averaged before transfer to sigma 
C  coordinates.
      DO 50 K=1,KBM1
      DO 50 J=1,JM
      DO 50 I=1,IM
  50  RMEAN(I,J,K)=RHO(I,J,K)
      DO 16 J=1,JM
      DO 16 I=1,IM
      TSURF(I,J)=TB(I,J,1)
  16  SSURF(I,J)=SB(I,J,1)
      DO 51 J=1,JM
      DO 51 I=1,IM
      UAB(I,J)=VEL*DUM(I,J)   
      VAB(I,J)=0.
      ELB(I,J)=0.0           
      ETB(I,J)=0.0           
      AAM2D(I,J)=AAM(I,J,1)
   51 CONTINUE
C------------------------------------------------
C  Set lateral boundary conditions for use in S.R.BCOND. 
C  In the seamount problem the east and west BCs are open 
C  whereas the south and north BCs are closed through specification
C  of the masks, FSM, DUM and DVM.
C------------------------------------------------
      RFE=1.
      RFW=1.
      RFN=1.
      RFS=1.
      ELE(1)=0.
      ELW(1)=0.
      DO 60 J=2,JMM1
      UABW(J)=UAB(2,J)
      UABE(J)=UAB(IMM1,J)
C  Set geostrophically conditioned elevations at the boundaries.
      ELE(J)=ELE(J-1)-COR(IMM1,J)*UAB(IMM1,J)/GRAV*DY(IMM1,J-1)
      ELW(J)=ELW(J-1)-COR(2,J)*UAB(2,J)/GRAV*DY(2,J-1)
  60  CONTINUE
C  Adjust boundary elevations so they are zero in the middle 
C  of the channel.
      ELEJMID=ELE(JMM1/2)
      ELWJMID=ELW(JMM1/2)
      DO 61 J=2,JMM1
      ELE(J)=(ELE(J)-ELEJMID)*FSM(IM,J)
      ELW(J)=(ELW(J)-ELWJMID)*FSM(2,J)
  61  CONTINUE
C
C********************************************************************
C  For the seamount and  (other possible applications) lateral
C  thermodynamic  boundary conditions are set equal to the initial conditions
C  and are held constant thereafter. Users can of course create varable 
C  boundary conditions.
C********************************************************************
      DO K=1,KBM1
      DO J=1,JM
      TBE(J,K)=TB(IM,J,K)
      TBW(J,K)=TB(1,J,K)
      SBE(J,K)=SB(IM,J,K)
      SBW(J,K)=SB(1,J,K)
      ENDDO
      DO I=1,IM
      TBN(I,K)=TB(I,JM,K)
      TBS(I,K)=TB(I,1 ,K)
      SBN(I,K)=SB(I,JM,K)
      SBS(I,K)=SB(I,1 ,K)
      ENDDO
      ENDDO
      RETURN
      END
c             
      SUBROUTINE VERTVL(DTI2)
       INCLUDE 'comblk97.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=2,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