*DECK MAIN
c gomtide.f 
c this program is set for tides in the gulf of mexico 
c                      provided by Bill O'Connor (NOAA/NCEP) 
      PROGRAM MAIN
C
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 THIS VERSION HAS BEN CONVERTED TO STANDARD FORTRAN 77 BY STEVE        *
C BRENNER WITH SOME ADDITIONAL CHANGES BY ME IN PROFQ WHICH IS NOW      *
C SOMEWHAT SIMPLER. HALF PRECISION HAS BEEN REMOVED.                    *
C THE DIVISOR 4. HAS BEEN ELIMINATED FROM COR4, CURV4 AND CURV42D AND   *
C THEY HAVE BEEN RENAMED COR, ETC.  A SMALL DIFFERENCING CHANGE HAS     *
C HAS BEEN MADE IN SUROUTINE BAROPG. A USER'S GUIDE IS AVAILABLE:       *
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
      INCLUDE 'comdec'    
      include 'comdec.tides'
      include 'comdec.astro' 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      DIMENSION 
     1     ADVUA(IM,JM),ADVVA(IM,JM),
     3     ADVUU(IM,JM),ADVVV(IM,JM)
      REAL ALON2(IM,JM),ALAT2(IM,JM),H2(IM,JM)     
c set for tides
      real coran(im,jm),copha(im,jm),htide(im,jm)
      real*8 T,h0,s0,p0 
      real*8 pi
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      DATA PI/3.1415926535D0/,RAMP/1.E0/,SMALL/1.E-10/
      DATA BETA/1.98E-11/,GRAV/9.806E0/
c-------------------------------------------------------------------
c initialize arrays
      do j=1,jm
      do i=1,im
      uaf(i,j)=0.0
      ua(i,j)=0.0
      uab(i,j)=0.0
      vaf(i,j)=0.0
      va(i,j)=0.0
      vab(i,j)=0.0
      elf(i,j)=0.0
      el(i,j)=0.0
      elb(i,j)=0.0
      advua(i,j)=0.0
      advva(i,j)=0.0
      advuu(i,j)=0.0
      advvv(i,j)=0.0
      fluxua(i,j)=0.0
      fluxva(i,j)=0.0
      wusurf(i,j)=0.0
      wvsurf(i,j)=0.0  
      wubot(i,j)=0.0
      wvbot(i,j)=0.0 
      aam2d(i,j)=50.0
      enddo
      enddo
c
c set for tides 
      do j=1,jm 
      do i=1,im  
      coran(i,j)=0.0
      copha(i,j)=0
      htide(i,j)=0.0
      enddo
      enddo
c
      do i=1,im
      ampm2s(i)=0.0
      pham2s(i)=0.0
      ampn2s(i)=0.0
      phan2s(i)=0.0
      amps2s(i)=0.0
      phas2s(i)=0.0
      ampk2s(i)=0.0
      phak2s(i)=0.0 
      ampo1s(i)=0.0
      phao1s(i)=0.0 
      ampk1s(i)=0.0
      phak1s(i)=0.0 
      ampp1s(i)=0.0 
      phap1s(i)=0.0
      ampq1s(i)=0.0
      phaq1s(i)=0.0  
c
      ampm2n(i)=0.0
      pham2n(i)=0.0
      ampn2n(i)=0.0
      phan2n(i)=0.0
      amps2n(i)=0.0
      phas2n(i)=0.0
      ampk2n(i)=0.0
      phak2n(i)=0.0 
      ampo1n(i)=0.0
      phao1n(i)=0.0 
      ampk1n(i)=0.0
      phak1n(i)=0.0 
      ampp1n(i)=0.0 
      phap1n(i)=0.0
      ampq1n(i)=0.0
      phaq1n(i)=0.0  
c
      eln(i)=0.0
      els(i)=0.0 
      covrhn(i)=0.0
      covrhs(i)=0.0   
      enddo
c
      do j=1,jm
      ampm2e(j)=0.0
      pham2e(j)=0.0
      ampn2e(j)=0.0 
      phan2e(j)=0.0 
      amps2e(j)=0.0 
      phas2e(j)=0.0 
      ampk2e(j)=0.0
      phak2e(j)=0.0
      ampo1e(j)=0.0 
      phao1e(j)=0.0 
      ampk1e(j)=0.0 
      phak1e(j)=0.0 
      ampp1e(j)=0.0 
      phap1e(j)=0.0
      ampq1e(j)=0.0
      phaq1e(j)=0.0 
c
      ele(j)=0.0 
      covrhe(j)=0.0
      enddo
c-----------------------------------------------------------
c--------------------------------------------------------------------
c read in file with Schwiderski tidal BC 
c aplitude in meters, phase in degrees 
      open(unit=10,form='formatted',status='old',
     1     file='gal_bc.dat')  
  997 format(2f10.4)   
      do i=1,im
      read(10,997) ampm2s(i),pham2s(i) 
      enddo 
      do i=1,im
      read(10,997) amps2s(i),phas2s(i) 
      enddo 
      do i=1,im
      read(10,997) ampk1s(i),phak1s(i) 
      enddo 
      do i=1,im
      read(10,997) ampo1s(i),phao1s(i) 
      enddo 
      do i=1,im
      read(10,997) ampn2s(i),phan2s(i) 
      enddo 
      do i=1,im
      read(10,997) ampp1s(i),phap1s(i) 
      enddo 
      do i=1,im
      read(10,997) ampk2s(i),phak2s(i) 
      enddo 
      do i=1,im
      read(10,997) ampq1s(i),phaq1s(i) 
      enddo 
c
      do j=1,jm
      read(10,997) ampm2e(j),pham2e(j)
      enddo
      do j=1,jm
      read(10,997) amps2e(j),phas2e(j)
      enddo
      do j=1,jm
      read(10,997) ampk1e(j),phak1e(j)
      enddo
      do j=1,jm
      read(10,997) ampo1e(j),phao1e(j)
      enddo
      do j=1,jm
      read(10,997) ampn2e(j),phan2e(j)
      enddo
      do j=1,jm
      read(10,997) ampp1e(j),phap1e(j)
      enddo
      do j=1,jm
      read(10,997) ampk2e(j),phak2e(j)
      enddo
      do j=1,jm
      read(10,997) ampq1e(j),phaq1e(j)
      enddo
c
c      do i=1,im
c      read(10,997) ampm2n(i),pham2n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) amps2n(i),phas2n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) ampk1n(i),phak1n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) ampo1n(i),phao1n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) ampn2n(i),phan2n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) ampp1n(i),phap1n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) ampk2n(i),phak2n(i) 
c      enddo 
c      do i=1,im
c      read(10,997) ampq1n(i),phaq1n(i) 
c      enddo 
c
      close(10)  
C--------------------------------------------------------------------
      IINT=0
      MODE=2
      NREAD=0
      DTE=20.0  
      IPRTD1=10000
      ISWTCH=0
      IPRTD2=10000
      IDAYS=15
      ISAVE=10000
      ISPADV=1
c      HORCON=.10
      HORCON=0.0
      SMOTH=0.1
      EIM=IM
      EJM=JM
C--------------------------------------------------------------------------
C       MODE = 2; 2-D CALCULATION (BOTTOM STRESS CALCULATED IN ADVAVE)
C       DTE = EXTERNAL TIME STEP
C       NREAD=0, NO RESTART INPUT FILE; NREAD=1, RESTART
C       IPRTD1 = PRINT INTERVAL IN DAYS; AFTER ISWTCH DAYS, PRINT
C                  INTERVAL = IPRTD2
C       IDAYS = LENGTH OF RUN IN DAYS
C       ISPADV = STEP INTERVAL WHERE EXTERNAL MODE ADVECTIVE TERMS ARE
C                  NOT UPDATED
C       HORCON = CONSTANT IN SMAGORINSKY HORIZONTAL DIFFUSIVITY
C       TPNU = HORIZONTAL DIFFUSIVITY PRANDT NUMBER
C       SMOTH = CONSTANT IN TIME SMOOTHER TO PREVENT SOLUTION SPLITTING
C
C--------------------------------------------------------------------------
      DTE2=DTE*2
      IPRINT=IPRTD1*24*3600/IFIX(DTE)
      ISWTCH=ISWTCH*24*3600/IFIX(DTE)
      IEND=IDAYS*24*3600/IFIX(DTE)
c 
      ISAVE=ISAVE*24*3600/IFIX(DTE)   
      WRITE(6,7030) MODE,DTE,IEND,ISPADV,IPRINT,SMOTH,HORCON
 7030 FORMAT(//,' MODE =',I3,
     1 '  DTE =',F7.2,'     IEND =',I6,
     2  '    ISPADV =',I6,'   IPRINT =',I6,/,
     3  '   SMOTH =',F7.2,'   HORCON =',F7.3/)
      WRITE(6,'('' NREAD ='',I5,//)') NREAD
C
C----------------------------------------------------------------------
C CALCULATE ASTRONOMICAL FORCING TIME DEPENDENT TERMS 
c----------------------------------------------------
c jday0 = julian day when tide phases are initialized  
c         integer, =0 at 00 Z  01 Jan 
c set for tides 
      TIME=0.
      TIME0=0. 
      IYEAR=1995
      JDAY0=0 
c calculate astronomical forcing time dependent terms 
c this is set for special case of starting 00 Z on 1 Jan 95
c
      bjday=jday0 + 1 + 365*(iyear - 1975) + (iyear - 1973)/4 
      T=(27392.500528d0 + 1.0000000356d0*bjday)/36525.d0  
c
      h0=279.69668d0+36000.768930485d0*t+3.03d-4*t**2
      s0=270.434358d0+481267.88314137d0*t-0.001133d0*t**2+1.9d-6*t**3
      p0=334.329653d0+4069.0340329575d0*t-0.010325d0*t**2-1.2d-5*t**3
C-----astronomical arguments--------------------
      chim2 =2.0d0*pi*(h0-s0)/180.d0     
      chis2=0.d0                
      chin2=pi*(2.d0*h0-3.d0*s0+p0)/180.d0   
      chik1=pi*(h0+90.d0)/180.d0           
      chio1=pi*(h0-2.d0*s0-90.d0)/180.d0    
      chip1=-pi*(h0+90.d0)/180.d0         
      chik2=2.d0*pi*h0/180.d0 
      chiq1=pi*(h0-3.d0*s0+p0-90.d0)/180.d0 
c
      write(6,7033) jday0,chim2  
 7033 format(10x,'JDAY0 = ',i10,'  chim2 = ',f20.8) 
c end of specifying tidal constants. the phase angles are in 
c the common block /astro/ to be communicated to subroutine BCOND   
C       
C----------------------------------------------------------------------
C             ESTABLISH PROBLEM CHARACTERISTICS
C          ****** ALL UNITS IN M.K.S. SYSTEM ******
C      F,BLANK AND B REFERS TO FORWARD,CENTRAL AND BACKWARD TIME LEVELS.
C----------------------------------------------------------------------
C
      DAYI=1.E0/86400.E0
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC      OPEN BATHYMETRY FILE
      OPEN(UNIT=20,FORM='FORMATTED',STATUS='OLD',
     1     FILE='bathy_gal.dat') 
C READ LON, LAT, AND BATHYMETRY 
c      read(20,995) alon,alat
  995 format(8f9.5) 
      read(20,996) h
  996 format(10f7.1)  
      close(20) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c OPEN THE FILE FOR MODEL OUTPUT 
      OPEN(UNIT=30,FORM='FORMATTED',STATUS='UNKNOWN',
     1     FILE='gom2d.dat')
c open the  file for station water level data
      open(unit=40,form='formatted',status='unknown',
     1     file='stations.dat') 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SET UP  GRID AND INITIAL CONDITIONS FOR STAND-ALONE TEST RUN
      DO 12 J=1,JM
      DO 12 I=1,IM
      EI=I 
      EJ=J
      ALON(I,J)=262.42+(EI-1.0)/12.0 
      ALON2(I,J)=ABS(ALON(I,J)-360.0) 
      ALAT(I,J)=25.91+(EJ-1.0)/12.0  
      AL1=25.91+(EJ-1.0)/12.0
      AL1=AL1*PI/180.0
      DX(I,J)=2.0*PI*6.371E6*COS(AL1)/(360.0*12.0)  
      DY(I,J)=2.0*PI*6.371E6/(360.0*12.0) 
      COR(I,J)=2.0*(7.292E-5)*SIN(AL1) 
      ART(I,J)=DX(I,J)*DY(I,J)
   12 CONTINUE
      write(6,4443) dx(im/2,jm/2),dy(im/2,jm/2),cor(im/2,jm/2)
 4443 format(5x,'midpt dx = ',e10.3,' dy = ',e10.3,' cor = ',e10.3)
      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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC       READ IN BATHYMETRY
Cwpo      READ(20,950) H
C
c      DO 14 J=1,JM
c      DO 14 I=1,IM
c      H(I,J)=2000.E0
c      EJ=J
c      S=ALOG(5.0)/EJM
c      H(I,J)=20.0*EXP(S*(EJ-1.0))
c   14 CONTINUE 
c      DO 1121 J=1,JM
c      H(1,J)=1.0
c      H(IM,J)=1.0
c 1121 CONTINUE
c      DO 1122 I=1,IM
c      H(I,1)=1.0
C      H(I,JM)=1.0
c 1122 CONTINUE  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc
C------------------------------------------------------------------------
C     READ IN GRID DATA AND INITIAL CONDITIONS
C------------------------------------------------------------------------
C     REWIND(40)
C     REWIND(60)
C     READ(40) ALON,ALAT,DX,DY,H,COR,
C    1      ART,ARU,ARV
C    2      ELN,
C*****************************************************************************
c      PERIOD=DAYI*(2.E0*PI)/COR(IM/2,JM/2)
      PERIOD=5. 
      WRITE(6,990) PERIOD
  990 FORMAT(10X,'PERIOD = ',F10.4) 
      DO 21 J=1,JM
      DO 21 I=1,IM
  21  CURV2D(I,J)=COR(I,J)
      DO 30 J=1,JM
      DO 30 I=1,IM
      FSM(I,J)=1.E0
      DUM(I,J)=1.E0
      DVM(I,J)=1.E0
C     IF(H(I,J).GT.1.E0.AND.H(I,J).LT.50.E0) H(I,J)=50.E0
cwpo      IF(H(I,J).GT.1.E0) GO TO 30
cwpo      IF(H(I,J).GE.1.E0) GO TO 30
      IF(H(I,J).GE.1.5E0) GO TO 30
      FSM(I,J)=0.E0
      DUM(I,J)=0.E0
      DVM(I,J)=0.E0
 30   CONTINUE
      DO 35 J=1,JMM1
      DO 35 I=1,IM
      IF(FSM(I,J).EQ.0.E0.AND.FSM(I,J+1).NE.0.E0) DVM(I,J+1)=0.E0
 35   CONTINUE
      DO 40 J=1,JM
      DO 40 I=1,IMM1
      IF(FSM(I,J).EQ.0.E0.AND.FSM(I+1,J).NE.0.E0) DUM(I+1,J)=0.E0
 40   CONTINUE
C
      DO 41 J=1,JM
      DUM(1,J)=0.E0
  41  DVM(1,J)=0.E0
      DO 45 J=1,JM
      DO 45 I=1,IM
      CBCMIN=.0025E0
      CBC(I,J)=0.0025
   45 CONTINUE
C
      DO 46 J=1,JM
      DO 46 I=1,IM
      ALAT2(I,J)=ALAT(I,J)*FSM(I,J) 
      ALON2(I,J)=ALON2(I,J)*FSM(I,J)   
      H2(I,J)=H(I,J)*FSM(I,J)   
   46 CBC(I,J)=CBC(I,J)*FSM(I,J)
      IP=IM
      ISK=1
c      CALL PRXY(' DX ',TIME,DX,IM,10,JM,10,1.0)
c      CALL PRXY(' DY ',TIME,DY,IM,10,JM,10,1.0)
c      CALL PRXY('COR ',TIME,COR,IM,10,JM,10,0.0)
        CALL PRXY('  TOPOGRAPHY   ',TIME, H2(1,1),IP,ISK,JM,1,1.E0)
c        CALL PRXY(' ALON2',TIME,ALON2,IM,1,JM,1,0.0)
c        CALL PRXY(' ALAT2',TIME,ALAT2,IM,1,JM,1,0.0)   
c        CALL PRXY(' FSM            ',TIME,FSM(1,1),IP,ISK,JM,1,1.E0)
c        CALL PRXY(' DUM            ',TIME,DUM(1,1),IP,ISK,JM,1,1.E0)
c        CALL PRXY(' DVM            ',TIME,DVM(1,1),IP,ISK,JM,1,1.E0)
c        CALL PRXY(' CBC            ',TIME,CBC(1,1),IP,ISK,JM,1,0.E0)
C
        DO 47 J=1,JM
        DO 47 I=1,IM
   47   TPS(I,J)=0.5E0/SQRT(1.E0/DX(I,J)**2+1.E0/DY(I,J)**2)
     1         /SQRT(GRAV*H(I,J))*FSM(I,J)
c        CALL PRXY(' CFL TIME STEP   ',TIME,TPS(1,1),IP,ISK,JM,1,1.E0)
c        CALL PRXY(' COR           ',TIME,COR(1,1),IP,ISK,JM,1,0.E0)
C------------------------------------------------------------------------
C             DEFINE LATERAL B.C.'S
C------------------------------------------------------------------------
C
      DO 48 I=1,IM
c      covrhs(i)=0.1*sqrt(grav/h(i,1)) 
      covrhs(i)=sqrt(grav/h(i,1)) 
   48 COVRHN(I)=.1E0*SQRT(GRAV/H(I,JMM2))
c
      do j=1,jm
c      covrhe(j)=0.1*sqrt(grav/h(im,j))
      covrhe(j)=sqrt(grav/h(im,j))
      enddo  
c 
      DO 51 J=1,JM
      DO 51 I=1,IM
      UA(I,J)=UAB(I,J)
      VA(I,J)=VAB(I,J)
      EL(I,J)=ELB(I,J)
      D(I,J)=H(I,J)+EL(I,J)
   51 CONTINUE
C
C
      DO 81 J=1,JM
      DO 81 I=1,IM
      D(I,J)=H(I,J)+EL(I,J)
   81 CONTINUE
C
C    INTRODUCE SIMPLE WIND STRESS, VALUE IS NEGATIVE FOR WESTERLY WIND
      DO 85 J=1,JM
      DO 85 I=1,IM
      WUSURF(I,J)=0.00E-4 
      WVSURF(I,J)=0.00E-4 
   85 CONTINUE 
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('  FREE SURFACE ',TIME, ELB(1,1),IP,ISK,JM,1,1.E-3)
c      CALL FINDPSI
c      CALL PRXY('  VAB          ',TIME, VAB(1,1),IP,ISK,JM,1,0.E0)
c***********
c      STOP
C***********************************************************************
C                                                                      *
C                  BEGIN NUMERICAL INTEGRATION                         *
C                                                                      *
C***********************************************************************
C
                     DO 9000 IINT=1,IEND
C
      TIME=DAYI*DTE*FLOAT(IINT)+TIME0
      RAMP=TIME/PERIOD
      IF(RAMP.GT.1.E0) RAMP=1.E0
cc      RAMP=1.E0
c calculate the astronomical tide generating potential
      call tidebody(htide)  
C**********************************************************************
C     HOR VISC = HORCON*DX*DY*SQRT((DU/DX)**2+(DV/DY)**2
C                                 +.5*(DU/DY+DV/DX)**2)
C**********************************************************************
      DO 95 J=2,JMM1
      DO 95 I=2,IMM1
      AAM2D(I,J)=HORCON*DX(I,J)*DY(I,J)
     1          *SQRT( ((UAB(I+1,J)-UAB(I,J))/DX(I,J))**2
     2                 +((VAB(I,J+1)-VAB(I,J))/DY(I,J))**2
     3 +.5E0*(.25E0*(UAB(I,J+1)+UAB(I+1,J+1)-UAB(I,J-1)-UAB(I+1,J-1))
     4              /DY(I,J)
     5       +.25E0*(VAB(I+1,J)+VAB(I+1,J+1)-VAB(I-1,J)-VAB(I-1,J+1))
     6              /DX(I,J)) **2)
  95  CONTINUE
      DO 405 J=2,JM
      DO 405 I=1,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=1,IM
      if(fsm(i,j).eq.0.e0) go to 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)
  410 continue 
C
      CALL BCOND(1)
C
      IF(MOD(IINT,ISPADV).EQ.0) CALL ADVAVE(ADVUA,ADVVA,MODE)
C
      ALPHA=0.225E0
      DO 420 J=2,JMM1
      DO 420 I=1,IM
      if(dum(i,j).eq.0.e0) go to 420 
      UAF(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)) )
     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)) )
     6      +ARU(I,J)*( RAMP*WUSURF(I,J)-WUBOT(I,J)   )
  420 continue 
c
c apply the tidal forcing 
      do j=2,jmm1
      do i=1,im
      if(dum(i,j).eq.0.e0) go to 421
      uaf(i,j)=uaf(i,j)
     1        -0.250*grav*(dy(i,j)+dy(i-1,j))*(d(i,j)+d(i-1,j))* 
     2        ramp*(htide(i,j)-htide(i-1,j)) 
  421 continue 
      enddo
      enddo  
c 
      DO 425 J=2,JMM1
      DO 425 I=1,IM
      if(dum(i,j).eq.0.e0) go to 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))
  425 continue 
      DO 430 J=3,JMM1
      DO 430 I=1,IM
      if(dvm(i,j).eq.0.e0) go to 430 
      VAF(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)) )
     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)) )
     6    + ARV(I,J)*( RAMP*WVSURF(I,J)-WVBOT(I,J)   )
  430 continue  
c 
c apply the tidal forcing 
      do j=3,jmm1
      do i=1,im
      if(dvm(i,j).eq.0.e0) go to 431
      vaf(i,j)=vaf(i,j)
     1         -0.250*grav*(dx(i,j)+dx(i,j-1))*(d(i,j)+d(i,j-1))*
     2          ramp*(htide(i,j)-htide(i,j-1))
  431 continue 
      enddo
      enddo 
c 
      DO 435 J=3,JMM1
      DO 435 I=1,IM
      if(dvm(i,j).eq.0.e0) go to 435 
      VAF(I,J)=
     1        ((H(I,J)+ELB(I,J)+H(I,J-1)+ELB(I,J-1))*VAB(I,J)*ARV(I,J)
     2              -4.E0*DTE*VAF(I,J))
     3       /((H(I,J)+ELF(I,J)+H(I,J-1)+ELF(I,J-1))*ARV(I,J))
  435 continue 
      CALL BCOND(2)
C
C  TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP
C
      VAMAX=-1.E10
      VAMIN=1.E10
      DO 442 J=1,JM
      DO 442 I=1,IM
C      VAMAX=Q8SMAX(VAF(1,1;LIJ))
      IF(VAF(I,J).GE.VAMAX) THEN 
         VAMAX=VAF(I,J)
         IMAX=I
         JMAX=J 
      ENDIF
C      VAMIN=Q8SMIN(VAF(1,1;LIJ))
      IF(VAF(I,J).LE.VAMIN) THEN 
         VAMIN=VAF(I,J)
         IMIN=I
         JMIN=J   
      ENDIF
  442 CONTINUE
      IF(ABS(VAMAX).GT.100.E0.OR.ABS(VAMIN).GT.100.E0) GO TO 9001
C
C       APPLY FILTER TO REMOVE TIME SPLIT AND 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---------------------------------------------------------------
c write amphidromic system data if last tide cycle
c      if(iint.gt.((ncyc-1)*nstep)) then
c       if(iint.gt.62560) then
       if(time.gt.14.0) then
        do j=1,jm
        do i=1,im
            if(elf(i,j).gt.coran(i,j)) then
            coran(i,j)=elf(i,j) 
c            copha(i,j)=(float(iint-(ncyc-1)*nstep)/enstep)*360.0  
            copha(i,j)=(float(iint-62560)/2236.)*360.0  
            endif 
        enddo
        enddo 
      endif 
C---------------------------------------------------------------
C           FILE SAVE SECTION
C----------------------------------------------------------------
      IF(MOD(IINT,ISAVE).EQ.0) THEN
        WRITE(30,950) TIME,H,ELB,UAB,VAB
      ENDIF
  950 FORMAT(8(1PE10.3))
c
c write station data to a file 
c grid  (4,23)  = Corpus Christi, TX
c grid (28,37)  = Freeport, TX  
c grid (35,42)  = Galveston Pleasure Pier, TX
c grid (46,46)  = Sabine Pass, TX
c
      if(mod(iint,180).eq.0) then
       write(40,4328) time,
c     1       elb(9,119),elb(33,132),elb(40,137),elb(51,141), 
c     2       elb(98,136),elb(122,148),elb(147,147),elb(163,143),
c     3       elb(180,133),elb(182,121),elb(185,117),elb(193,102),
c     4       elb(195,97),elb(205,81),elb(203,81),elb(195,80)      
     1   elb(4,23),elb(28,37),elb(35,42),elb(46,46)   
      endif 
 4328 format(4(1pe10.3)) 
c
      if(mod(iint,180).eq.0) then
      write(6,4338) htide(100,50)
 4338 format(10x,' htide(100,50) = ',e20.8) 
      endif
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/IFIX(DTE)
        WRITE(6,'(/,'' TIME ='',F10.2,''    IINT ='',I8,
     1  ''   IEXT ='',I8,''   IPRINT ='',I5,//)') TIME,IINT,IEXT,IPRINT
C       CALL PRXY(' ADVUA     ',TIME,ADVUA(1,1),IP,ISK,JM,1,.0001S9)
C       CALL PRXY(' ADVVA     ',TIME,ADVVA(1,1),IP,ISK,JM,1,.0001S9)
C       CALL PRXY(' ADVUU     ',TIME,ADVUU(1,1),IP,ISK,JM,1,.0001S9)
C       CALL PRXY(' ADVVV     ',TIME,ADVVV(1,1),IP,ISK,JM,1,.0001S9)
c        CALL PRXY(' WUBOT     ',TIME,WUBOT(1,1),IP,ISK,JM,1,0.E0)
c        CALL PRXY(' WVBOT     ',TIME,WVBOT(1,1),IP,ISK,JM,1,0.E0)
c        CALL PRXY('  FREE SURFACE ',TIME, ELB(1,1),IP,ISK,JM,1,0.E0)
c        CALL PRXY(' AVERAGE U COMP.',TIME,UAB(1,1),IP,ISK,JM,1,0.E0)
c        CALL PRXY(' AVERAGE V COMP.',TIME,VAB(1,1),IP,ISK,JM,1,0.E0)
c        CALL FINDPSI
C
C     WRITE(80) TIME,UAB,VAB,ELB
C
      VTOT=0.E0
      DO 8888 J=1,JM
      DO 8888 I=1,IM
      DVTOT=DX(I,J)*DY(I,J)*D(I,J)
      VTOT=VTOT+DVTOT
 8888 CONTINUE
      WRITE(6,'(''  VTOT =,'',E20.8,''    TTOT ='',E20.8,
     1      ''STOT ='',E20.8)') VTOT
      IF(ABS(VAMAX).GT.100.E0.OR.ABS(VAMIN).GT.100.E0) THEN
        CALL PRXY('  EL ',TIME, EL (1,1),IP,ISK,JM,1,0.E0)
        CALL PRXY('  ELF',TIME, ELF(1,1),IP,ISK,JM,1,0.E0)
        CALL PRXY('  UA            ',TIME,UA (1,1),IP,ISK,JM,1,0.E0)
        CALL PRXY('  VA            ',TIME,VA (1,1),IP,ISK,JM,1,0.E0)
        CALL PRXY('  UAF           ',TIME,UAF(1,1),IP,ISK,JM,1,0.E0)
        CALL PRXY('  VAF           ',TIME,VAF(1,1),IP,ISK,JM,1,0.E0)
      WRITE(6,'(//////////////////)')
      WRITE(6,'(''************************************************'')')
      WRITE(6,'(''************ ABNORMAL JOB END ******************'')')
      WRITE(6,'(''************* USER TERMINATED ******************'')')
      WRITE(6,'(''************************************************'')')
      WRITE(6,'(1X,2E12.3)') VAMAX,VAMIN
      write(6,1918) imax,jmax,imin,jmin
 1918 format(5x,'imax=',i5,' jmax=',i5,' imin=',i5,' jmin=',i5) 
      STOP
      ENDIF
 7000 CONTINUE
C--------------------------------------------------------------------
C             END PRINT SECTION
C--------------------------------------------------------------------
c write time step and model simulation time to file
      open(unit=27,form='formatted',status='unknown',
     1     file='gom2d.ck')
      write(27,1919) iint,time
 1919 format(10x,'IINT = ',i10,' TIME = ',f10.6)  
      close(27)
c
 9000                     CONTINUE
c--------------------------------------------------
ctide
c end of simulation. write amphidromic data to file to save 
C
c tides-open file for amphidromic system data
      open(unit=50,form='formatted',status='unknown',
     1     file='amphidrom.dat')       
      write(50,950) h,coran,copha 
      close(50) 
c--------------------------------------------
      STOP
      END
*DECK BCOND
      SUBROUTINE BCOND(IDX)
      INCLUDE 'comdec'
      include 'comdec.tides'
      include 'comdec.astro'  
C --------------------------------------------------------------------
      DATA PI/3.1415926535/,GEE/9.807E0/
c----------------------------------------------------
c specify tidal constants in s**(-1) 
      sigm2=1.40519e-4
      sigs2=1.45444e-4  
      sigp1=0.72523e-4 
      sigk1=0.72921e-4
      sigo1=0.67598e-4 
      sigk2=1.45842e-4 
      sigq1=0.64959e-4 
      sign2=1.37880e-4 
c
      perm2=(2.0*pi/sigm2)/86400.
      pers2=(2.0*pi/sigs2)/86400.
      perp1=(2.0*pi/sigp1)/86400.
      perk1=(2.0*pi/sigk1)/86400.
      pero1=(2.0*pi/sigo1)/86400.
      perk2=(2.0*pi/sigk2)/86400.
      perq1=(2.0*pi/sigq1)/86400. 
      pern2=(2.0*pi/sign2)/86400.      
c-----------------------------------------------------------------
C
      GO TO (10,20), IDX
C
 10   CONTINUE
C-----------------------------------------------------------------------
C                   EXTERNAL ELEV. B.C.'S
C-----------------------------------------------------------------------
C EASTERN BNDRY 
      DO 100 J=1,JM
c      ELF(IM,J)=ELF(IMM1,J)
c tidal bndry conditions on east boudnary  
c compute Schwiderski tidal BC on eastern bndry
      ele(j)= 
     1       0.6*ramp*ampm2e(j)*cos(2.0*pi*(time-jday0)/perm2  
     1                      +chim2 -(pham2e(j)-30.)*pi/180.0)    
c     2      +1.15*ramp*ampk1e(j)*cos(2.0*pi*(time-jday0)/perk1  
c     2                      +chik1 -(phak1e(j)-10.)*pi/180.0)    
c     3      +0.6*ramp*amps2e(j)*cos(2.0*pi*(time-jday0)/pers2 
c     3                      +chis2 -(phas2e(j)+20.)*pi/180.0)    
c     4      +0.85*ramp*ampo1e(j)*cos(2.0*pi*(time-jday0)/pero1  
c     4                      +chio1 -(phao1e(j)-25.)*pi/180.0)    
c     5      +0.66*ramp*ampn2e(j)*cos(2.0*pi*(time-jday0)/pern2  
c     5                      +chin2 -(phan2e(j)-40.)*pi/180.0)    
c     6      +0.8*ramp*ampp1e(j)*cos(2.0*pi*(time-jday0)/perp1  
c     6                      +chip1 -(phap1e(j)-36.)*pi/180.0)    
c     7      +ramp*ampk2e(j)*cos(2.0*pi*(time-jday0)/perk2  
c     7                      +chik2 -phak2e(j)*pi/180.0)    
c     8      +0.8*ramp*ampq1e(j)*cos(2.0*pi*(time-jday0)/perq1  
c     8                      +chiq1 -(phaq1e(j)-45.)*pi/180.0)    
      elf(im,j)=ele(j) 
c      elf(im,j)=elf(imm1,j)
      elf(imm1,j)=elf(im,j) 
  100 CONTINUE
C NORTHERN AND SOUTHERN BNDRY    
      DO 110 I=1,IM
c      ELF(I,JM)=ELF(I,JMM1)
c      ELF(I,1)=ELF(I,2)
c tidal boundary condition on south bndry
      els(i)=
     1       0.6*ramp*ampm2s(i)*cos(2.0*pi*(time-jday0)/perm2 
     1                     +chim2 -(pham2s(i)-30.)*pi/180.0)  
c     2      +1.15*ramp*ampk1s(i)*cos(2.0*pi*(time-jday0)/perk1  
c     2                     +chik1 -(phak1s(i)-10.)*pi/180.0)  
c     3      +0.6*ramp*amps2s(i)*cos(2.0*pi*(time-jday0)/pers2  
c     3                     +chis2 -(phas2s(i)+20.)*pi/180.0)  
c     4      +0.85*ramp*ampo1s(i)*cos(2.0*pi*(time-jday0)/pero1  
c     4                     +chio1 -(phao1s(i)-25.)*pi/180.0)  
c     5      +0.66*ramp*ampn2s(i)*cos(2.0*pi*(time-jday0)/pern2  
c     5                     +chin2 -(phan2s(i)-40.)*pi/180.0)  
c     6      +0.8*ramp*ampp1s(i)*cos(2.0*pi*(time-jday0)/perp1 
c     6                     +chip1 -(phap1s(i)-36.)*pi/180.0)  
c     7      +ramp*ampk2s(i)*cos(2.0*pi*(time-jday0)/perk2  
c     7                     +chik2 -phak2s(i)*pi/180.0)  
c     8      +0.8*ramp*ampq1s(i)*cos(2.0*pi*(time-jday0)/perq1   
c     8                     +chiq1 -(phaq1s(i)-45.)*pi/180.0)  
      elf(i,1)=els(i) 
c      elf(i,1)=elf(i,2) 
      elf(i,2)=elf(i,1)
  110 CONTINUE
C    
check  
      if(mod(iint,360).eq.0) then
      amplit=cos(2.0*pi*(time-jday0)/0.99727 + chim2) 
c      write(6,9943) time,chim2,amplit 
      write(6,9943) time,chik1,amplit 
 9943 format(5x,'time = ',f9.4,' chim2= ',f14.5,' amp = ',f7.4) 
      endif 
check
      DO 111 J=1,JM
      DO 111 I=1,IM
 111  ELF(I,J)=ELF(I,J)*FSM(I,J)
      RETURN
C
 20   CONTINUE
C-----------------------------------------------------------------------
C                   EXTERNAL VEL B.C.'S
C-----------------------------------------------------------------------
C                   ***** SEAWARD ********
c east 
      do j=1,jm
c      vaf(im,j)=va(im,j)-sqrt(9.8*h(im,j))*(dte/dx(im,j))*
c     1          (va(im,j)-va(imm1,j)) 
      vaf(im,j)=vaf(imm1,j) 
c      vaf(im,j)=vab(im,j)-sqrt(9.8*h(im,j))*(dte/dx(im,j))*
c     1          (vaf(im,j)+vab(im,j)-2.0*va(imm1,j))  
c
c      uaf(im,j)=ua(im,j)-sqrt(9.8*h(im,j))*(dte/dx(im,j))*
c     1          (ua(im,j)-ua(imm1,j)) 
c      uaf(im,j)=uab(im,j)-sqrt(9.8*h(im,j))*(dte/dx(im,j))*
c     1          (uaf(im,j)+uab(im,j)-2.0*ua(imm1,j))  
      uaf(im,j)=covrhe(j)*(el(imm1,j)-ele(j))  
      enddo
C     DO 120 J=1,JM
C     UAF(IM,J)=RAMP*UABE(J)+COVRHE(J)*(EL(IMM1,J)-ELE(J))
C     UMID=.5E0*(UA(IM,J)+UA(IM,J-1))
C     VAF(IM,J)=VA(IM,J)-DTI/(DX(IM,J)+DX(IMM1,J))
C    1     *(  (UMID+ABS(UMID))*(VA(IM,J)-VA(IMM1,J))
C    2        +(UMID-ABS(UMID))*(0.E0 -VA(IM,J))  )
C*120 VAF(IM,J)=0.E0
C                   **** NORTH AND SOUTH ****
C--------- SPONGE B.C. ----------------------------------------------
C     DO 124 I=1,IMM1
C124  VAF(I,JMM1)=RAMP*VABN(I)
C    1        *(H(I,JMM1)+H(I,JMM2))
C    2        /(H(I,JMM1)+ELF(I,JMM1)+H(I,JMM2)+ELF(I,JMM2))
C     JSPG=JMM2-4
C     DO 125 J=JMM2-4,JMM2
C     DO 125 I=1,IMM1
C     ALPH2DT=.0001*2.*DTE
C     UAF(I,J)=(1.-ALPH2DT)*UAF(I,J)
C 125 VAF(I,J)=(1.-ALPH2DT)*VAF(I,J)+ALPH2DT*VAF(I,JMM1)
C--------------------------------------------------------------------
C     DO 124 I=1,IMM1
C 124 VAF(I,JM)=RAMP*VABN(I)+COVRHN(I)*(EL(I,JMM1)-ELN(I))
C     DO 130 I=1,IMM1
C     VAF(I,2)=RAMP*VABS(I)
C    1        *(H(I,2)+H(I,1))/(H(I,2)+ELF(I,2)+H(I,1)+ELF(I,1))
C     VAF(I,1)=VAF(I,2)
C     VMID=.5E0*(VA(I,JM)+VA(I-1,JM))
C     UAF(I,JM)=UA(I,JM)-DTI/(DY(I,JM)+DY(I,JMM1))
C    1    *(   (VMID+ABS(VMID))*(UA(I,JM)-UA(I,JMM1))
C    2        +(VMID-ABS(VMID))*(0.E0-UA(I,JM))  )
C     VMID=.5E0*(VA(I,2)+VA(I-1,2))
C     UAF(I,1)=UA(I,1)-DTI/(DY(I,1)+DY(I,2))
C    1    *(   (VMID+ABS(VMID))*(UA(I,1)-0.E0)
C    2        +(VMID-ABS(VMID))*(UA(I,2)-UA(I,1)) )
C 130 CONTINUE
ccccccccccccccccccccccccccccccccc
c south 
      DO 135 I=1,IM
c      UAF(I,JM)=uaf(i,jmm1) 
c      uaf(i,1)=uaf(i,2) 
c      VAF(I,JM)=SQRT(GEE/H(I,JM))*EL(I,JM)
c      vaf(i,1)=-sqrt(gee/H(i,1))*el(i,1)  
c      vaf(i,1)=-covrhs(i)*(el(i,2)-els(i))  
      vaf(i,2)=-covrhs(i)*(el(i,2)-els(i))
      vaf(i,1)=vaf(i,2) 
c      VAF(I,JM)=WUSURF(I,JM)/(H(I,JM)*COR(I,JM))+
c     1          SQRT(GEE/H(I,JM))*EL(I,JM)
  135 CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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
      END
*DECK ADVAVE
      SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE)
      INCLUDE 'comdec'
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      DIMENSION ADVUA(IM,JM),ADVVA(IM,JM)
C---------------------------------------------------------------------
C      CALCULATE U-ADVECTION & DIFFUSION
C---------------------------------------------------------------------
C
C-------- ADVECTIVE FLUXES -------------------------------------------
C
      DO 299 J=1,JM
      DO 299 I=1,IM
 299  TPS(I,J)=0.E0
      DO 300 J=2,JM
      DO 300 I=1,IM
 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=1,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=1,IM
 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=1,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=1,IM
 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.0E0
C
C---------ADVECTIVE FLUXES ----------------------------
      DO 700 J=2,JM
      DO 700 I=1,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=1,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=1,JMM1
      DO  860 I=1,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,JMM1
      DO  870 I=1,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=3,JMM1
      DO  880 I=1,IM
 880  ADVVA(I,J)=FLUXUA(I+1,J)-FLUXUA(I,J)
     1          +FLUXVA(I,J)-FLUXVA(I,J-1)
C
      DO 881 J=1,JM
      ADVVA(1,J)=0.E0
 881  ADVVA(IM,J)=0.E0
C
C---------------------------------------------------------------
C     COMPUTE BOTTOM FRICTION IF RUN IS EXTERNAL MODE ONLY
C---------------------------------------------------------------
cwpo      IF(MODE.NE.2) GO TO 5000
      DO 100 J=2,JMM1
      DO 100 I=1,IM
      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,JMM2
      DO 102 I=1,IM
      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
 5000 CONTINUE
      RETURN
      END
*DECK PRXY
      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(600),LINE(600)
      CHARACTER LABEL*(*)
      DATA ZERO /1.E-10/
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
      SCALE=10.E0**(INT(ALOG10(AMX)+1.E2)-103)
  160 CONTINUE
      SCALEI=1.E0/SCALE
      WRITE(6,170) LABEL
  170 FORMAT(1X,A40)
      WRITE(6,180) TIME,SCALE
  180 FORMAT(' TIME =',F9.4,' DAYS     MULTIPLY ALL VALUES BY',1PE10.3)
      DO 190 I=1,IM
  190 NUM(I)=I
      IB=1
C
  200 CONTINUE
      IE=IB+23*ISKP
      IF(IE.GT.IM) IE=IM
      WRITE(6,210) (NUM(I),I=IB,IE,ISKP)
  210 FORMAT(/,2X,24I5,/)
      DO 260 J=1,JM,JSKP
      JWR=JM+1-J
      DO 220 I=IB,IE,ISKP
  220 LINE(I)=INT(SCALEI*A(I,JWR))
      WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP)
 
  240 FORMAT(1X,I3,24I5)
  260 CONTINUE
      WRITE(6,280)
  280 FORMAT(//)
      IF(IE.GE.IM) RETURN
      IB=IB+24*ISKP
      GO TO 200
      END
      subroutine tidebody(htide)    
c this subroutine computes the astronomical tide body forcing
c The time and astronomical argument of the tides is
c computed in the MAIN program at the start of the time
c integration, and passed to this subroutine (and BCOND)
c in common block file 'comdec.astro' 
c
      include 'comdec'   
      include 'comdec.astro' 
      real htide(im,jm) 
      data rad/0.0174533/,pi/3.1415926/ 
c
c set the amplitude of the astronomical tidal forcing
      ampm2=0.24233 
c 
c calculate the astronomical body forcing
      do 100 j=1,jm
      do 100 i=1,im
      htide(i,j)=
c m2 tide
     1      ampm2*cos(rad*alat(i,j))*cos(rad*alat(i,j))* 
     2 cos(2.0*pi*(time-jday0)/0.51754 + chim2 + 2.0*rad*alon(i,j))      
  100 continue   
c allow for tidal loading and self-attraction
      do 200 j=1,jm
      do 200 i=1,im   
      htide(i,j)=0.69*htide(i,j)-0.1*elb(i,j)
  200 continue 
c
      return 
      end