!----------------------------------------------------------------------!
!     This POM code has wetting and drying (WAD) capability            !
!                                                                      !
!     Set nwad=1 for WAD,set =0 to turn off WAD                        !
!     Set BOTH nwad=0 & nsmolar=0 to recover the original non-WAD POM  !
!                                                                      !
!     For details see:                                                 !
!                                                                      !
!     O2005 = Oey [2005; OceanModelling,  9, 133-150]                  !
!     O2006 = Oey [2006; OceanModelling, 13, 176-195]                  !
!     O2007 = Oey, Ezer et al [2007; Ocean Dynamics, 57, 205-221]      !
!                                                                      !
!     or, http://www.aos.princeton.edu/WWWPUBLIC/PROFS/                !
!     also, .../WWWPUBLIC/PROFS/pom98_wad_release_descriptions.html    !
!                                                                      !
!     For changes that I made, search for the following keyword:       !
!                                                                      !
      !lyo:!wad:  or simply "wad:"                                     !
!     !tne:!wad:  (changes by T. Ezer to run pom98 seamount case)      !
!                                                                      !
!     I recommed using "DOUBLE  PRECISION (i.e. REAL*8)" option to     !
!     compile, with or without WAD; e.g. w/Intel ifort, we use -r8:    !
!                                                                      !
!     ifort pom98_wad.f -o a.out -r8                                   !
!                                                                      !
!     The -r8 option is the same as specifying -real_size 64 or        !
!     -auto_double.  However, the code works with or without -r8       !
!                                                                      !
!     NOTE: for nsmolar=1, the PARAMETER (IM=???,JM=???) in subroutines!
!     advt2d & smol_adif must match that specified in comblk*.h        !
!                                                                      !
!     Send bugs and/or improvements to: lyo@princeton.edu              !
!                                                                      !
!                       ... l.oey (May/02/2007)                        !
!                                 (Jul/24,30/2007)                     !
!                                                                      !
!----------------------------------------------------------------------!
C***********************************************************************
C     pom98.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 ---------------------------------------------------------
      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 Sep 12, 2001  *
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
!lyo:!wad:Add WAD variables in comblk98_wad.h 
      INCLUDE 'comblk98_wad.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,KAPPA
      DATA PI/3.141592654/,SMALL/1.E-10/,KAPPA/0.4/
      GRAV=9.806
      TBIAS=0.
      SBIAS=0.
      RHOREF=1025.
      IINT=0
      TIME=0.
C
      SEAMT=.TRUE.
!----------------------------------------------------------------------!
!                                                                      !
!lyo:!wad:For 3-d runs w/rivers (strong stratification)                !
!     strong currents can (and physically they should) develop         !
!     in thin film of water over nearly dry cells, so dte and          !
!     isplit (below) may need to be smaller.  For dx~dy~1km,           !
!     I have used values as small as DTE=3.0 and ISPLIT=5.             !
!                                                                      !
!     Note#1: Since alpha=0 when nwad=1 (see below), the corresponding !
!             DTE needs to be smaller than the DTE used for nwad=0.    !
!     Note#2: ISPLIT should also be smaller than its value for no WAD; !
!             see Appendix A of Oey [2006; Ocean Modell. 13, 176-195]. !
!                                                                      !
!     For example, for the seamount problem that Tal Ezer has tested   !
!     (this version), he found that DTE=60 & ISPLIT=30 worked for      !
!     nwad=0 but blew up for nwad=1. However, by decreasing ISPLIT     !
!     =5->10, the model with nwad=1 works for DTE=20->60, except that  !
!     isolated cells w/thin fluid layers (~5cm) appear, surrounded by  !
!     cells which are mostly dry; this probably is to be expected for  !
!     the relatively large dx&dy >~4km used here.  Using smaller DTE=15!
!     & ISPLIT=5, or DTE=6 & ISPLIT=10, or DTE=10 & ISPLIT=5 (as used  !
!     below) can remove the isolated thin-fluid cells.                 !
!----------------------------------------------------------------------!
      DTE=10.
      ISPLIT=5
c     DTE=24.   ! works with nowad & HMAX=200m 
c     ISPLIT=30     
      DAYS=3.0
      PRTD1=0.125    ! 3h interval print  
      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'
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:set nwad=1 for WAD run, =0 to recover non-WAD run with min   !
!     depths set at say 10m                                            !
!                                                                      !
      nwad=1
!                                                                      !
!lyo:!wad:set nsmolar=1 if water depth is to be solved by Smolarkiewicz!
!     See O2005 or O2006.                                              !
!     NOTE: nwad=1 w/nsmolar=1 may require finer grid AND even a larger!
!     hc, below, say hc=0.1 instead of default hc=0.05                 !
!                                                                      !
      nsmolar=0; !nsmolar=nsmolar*nwad
!                                                                      !
      WRITE(6,'('' nwad,nsmolar ='',2i5,/)') nwad,nsmolar
!----------------------------------------------------------------------!
!lyo:!wad:                                                             !
!                                                                      !
!     HHI, HMIN & HC (all +ve) are defined as:                         !
!                                                                      !
!     HHI = HIghest water [above MSL] ever expected, this should       !
!           be the observed_highest + say 1 meter to make sure;        !
!     H = hmin ---> absolute land area [never flooded], FSM=0.0;       !
!     H >= hc  ---> water depth wrt High water level, FSM=1.0;         !
!                                                                      !
!     Choose hhi=20 (below) assuming that the highest tide or surge    !
!     is not expected to rise above 20m wrt MSL                        !
!                                                                      !
      hhi=20.0*float(nwad)
      hmin=0.01; hc=0.05; hco=hc       !hco is used to avoid round-
!     hc=hc*1.0001                     !off in initial elevation
      hc=hc*1.01                       !tne:!wad:- using a larger 
!                 multiplying factor, "1.01" instead of "1.0001",      !
!                 can help eliminate singular wet spots due to roundoff!
!                                                                      !
      WRITE(6,'('' hhi,hmin,hc,hco ='',1p4e13.4,/)') hhi,hmin,hc,hco
!                                                                      !
!----------------------------------------------------------------------!
!                                                                      !
C--------------------------------------------------------------------------
C    Herewith is a list of parameters that are user discretionary.
C--------------------------------------------------------------------------
C                      In comblk98_wad.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 = Bias value subtracted from T and S. For 32 bit
C                      arithmetic, setting S=35. and T=10. will reduce
C                      round-off error.
C       RHOREF=Reference density. For ocean aplications recommended value
C              is 1025; for fresh water, 1000.
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       Z0B = bottom roughness parameter = (Nikuradse equivalent
C               roughness)/30. Should not be larger than the vertical grid
C               spacing at the bottom
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))
!
!lyo:!wad:dry-cell T/S relaxation using leapfrog-trapezoidal:
      dtrat=dti/(0.5*86400.)  !1/2 day relaxation time-scale
      cwetrlx1=(1.-dtrat)/(1.+dtrat); cwetrlx2=2.*dtrat/(1.+dtrat)
!
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)
       write(6,'(''DO SEAMOUNT'')')
      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
!lyo:!wad:                                                             !
c     hkeep(:,:)=h(:,:) !Keep original for plots etc
      do j=1,jm; do i=1,im
      hkeep(i,j)=h(i,j) !Keep original for plots etc
      enddo; enddo
!                                                                      !
      if (nwad.eq.1) then
!                                                                      !
!lyo:!wad:                                                             !
!     The above "H" is wrt MSL at which level H=0, and H<0 for depth   !
!     below the MSL, and H>0 for land above the MSL (i.e. for hills    !
!     and mountains etc.  We now first reverse the sign of H, then     !
!     change the reference from MSL to ALB, the Absolute Land Boundary,!
!     according to O2006, Fig.1:                                       !
!                                                                      !
      do j=1,jm; do i=1,im
       h(i,j)=-h(i,j)+hhi
c
       FSM(I,J)=1.; DUM(I,J)=1.; DVM(I,J)=1.
      IF(H(I,J).LT.0.0)THEN !Define absolute land areas (never flood):
                            ! temporarily set H=-1.0 (some -ve number)
       H(I,J)=-1.0; FSM(I,J)=0.; DUM(I,J)=0.; DVM(I,J)=0.
      ENDIF
c
      enddo; enddo
c
      DO J=1,JM-1; DO I=1,IM
       IF(FSM(I,J).EQ.0..AND.FSM(I,J+1).NE.0.)DVM(I,J+1)=0.
      enddo; enddo
      DO J=1,JM;   DO I=1,IM-1
       IF(FSM(I,J).EQ.0..AND.FSM(I+1,J).NE.0.)DUM(I+1,J)=0.
      enddo; enddo
!                                                                      !
!     The above yields the following definitions for H:                !
! 
!     H=-1, absolute land areas [never flooded], FSM=0.0, always       !
!     H>=0, wet but potentially dry areas,  FSM also=1.0, always       !
!                                                                      !
!     We now define a wetmask, assuming that at t=0, the free surface  !
!     is at the MSL (i.e. elevation=0 wrt MSL).  If this run does not  !
!     begin with time=0 (or if other special free-surface distribution !
!     is desired), then wetmask needs to be redefined or read in e.g.  !
!     from a RESTART file, as in below when NREAD.ne.0                 !
!                                                                      !
      do j=1,jm; do i=1,im
      wetmask(i,j)=fsm(i,j)
      if (h(i,j).lt.hhi) wetmask(i,j)=0.0
      enddo; enddo
!                                                                      !
!     Finalize definitions of "H":                                     !
!                                                                      !
!     H=hmin, absolute land areas [never flooded], FSM=0.0, always     !
!     H>=hc,  wet but potentially dry areas,  FSM also=1.0, always     !
!                                        but wetmask can be 0 or 1     !
!                                                                      !
      do j=1,jm; do i=1,im
      if (h(i,j).lt.0.0) then
       h(i,j)=hmin
      else
       h(i,j)=h(i,j)+hc
      endif
      enddo; enddo
!                                                                      !
!lyo:!wad:
      else    !nwad.eq.0
C
      DO J=1,JM
        DO I=1,IM
      H(I,J)=-H(I,J) !MAKE SEA-DEPTH POSITIVE VALUE
         IF(H(I,J).LE.0.) H(I,J)=0.1
        ENDDO
      ENDDO
C
      DO J=1,JM
       DO I=1,IM
        IF(H(I,J).GT.0.1.and.H(I,J).LT.10.0) then
         H(I,J)=10.   !In order to avoid exposure of intertidal flats
       ENDIF
      ENDDO
      ENDDO
c
C---- If the Depth H(I,J).GT.0.1, then Sea
      DO J=1,JM
      DO I=1,IM
C      DT(I,J)=H(I,J)
      FSM(I,J)=0.E0
      DUM(I,J)=0.E0
      DVM(I,J)=0.E0
c
      IF(H(I,J).GT.0.1) FSM(I,J)=1.
c
      ENDDO
      ENDDO
C
      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
C
C
!lyo:!wad:                                                             !
      do j=1,jm; do i=1,im
      wetmask(i,j)=fsm(i,j)
      enddo; enddo
!
      endif   !if (nwad.eq.1) then...
!                                                                      !
!lyo:!wad:                                                             !
!     Print to check:                                                  !
      CALL PRXY(' Read-in H ',TIME, Hkeep(1,1),IM,1,JM,1,0.0)
      CALL PRXY(' H ',TIME, H(1,1),IM,1,JM,1,0.0)
      CALL PRXY(' FSM ',TIME, FSM(1,1),IM,1,JM,1,1.0)
      CALL PRXY(' WETMASK ',TIME, WETMASK(1,1),IM,1,JM,1,1.0)
      do j=1,jm; do i=1,im
      tps(i,j)=FSM(i,j)-WETMASK(i,j)
      enddo; enddo
      CALL PRXY('FSM-WETMASK',TIME,tps(1,1),IM,1,JM,1,1.0)
c
      if (nwad.eq.0) then
      do j=1,jm; do i=1,im
      if (fsm(i,j).ne.wetmask(i,j)) then
      WRITE(6,'('' Stopped, nwad ='',i4)') nwad
      WRITE(6,'('' fsm.ne.wetmask @ (i,j) ='',2i4,/)') i,j
c
      stop
c
      endif
      enddo; enddo
      endif
!----------------------------------------------------------------------!
      DO 51 J=1,JM
      DO 51 I=1,IM
      ELB(I,J)=0.0
      ETB(I,J)=0.0
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Note elevation is always <=0, measured downward from ALB     !
!                                                                      !
!     Set initial elevation to be flat MSL in wet regions, but set     !
!     initial elevation s.t. total water depth = H+ELB=HC in regions   !
!     where wetmask indicates dry but not land                         !
!                                                                      !
      ELB(I,J)=-hhi*fsm(i,j)
!lyo:!wad:note:The following "if" is not satisfied if nwad=0 since     !
!     then fsm=wetmask at all cells                                    !
!     if (fsm(i,j).ne.0.0.and.wetmask(i,j).eq.0.0) ! dry but not land
      if (fsm(i,j).ne.wetmask(i,j))                ! dry but not land
     &elb(i,j)=-h(i,j)+hco !note slightly smaller hco<hc to ensure that
                           !initially, cell is wet with elb+h=hco < hc
      ETB(I,J)=ELB(I,J)
!                                                                      !
!----------------------------------------------------------------------!
! 
      AAM2D(I,J)=AAM(I,J,1)
   51 CONTINUE
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Make sure that initial wetmask remains compatible with       !
!     initial elb.  Also define marsh evaporation, +ve upward or       !
!     evaporate.  Also initialize wriv for rivers.                     !
!lyo:!wad:note:wriv is added here but the implementation following     !
!     Oey [1996; JPO, 26, 145-175; but see p.154 in particular] is not !
!     complete in this code.  E-mail me if want to implement it.       !
!                                                                      !
      if (nwad.eq.1) then
      do j=1,jm; do i=1,im
      wetmask(i,j)=fsm(i,j)
      if ((h(i,j)+elb(i,j)).le.hc) wetmask(i,j)=0.0
      enddo; enddo
      endif
!                                                                      !
      do j=1,jm; do i=1,im
      wmarsh(i,j)=0.0; wriv(i,j)=0.0
      if (fsm(i,j).ne.wetmask(i,j)) then           ! dry but not land
!     2 examples of wmarsh:
!--   wmarsh(i,j)=fsm(i,j)*10./art(i,j)
!--   wmarsh(i,j)=fsm(i,j)*(4.e-7)  !Gill's [1982] value
      endif
      enddo; enddo
!                                                                      !
!----------------------------------------------------------------------!
! 
C
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))
      AAM(I,J,K)=500.
   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 -----------------------------------------------------------------
c     An empirical specification of the roughness parameter
c   * If the height of the nearest bottom grid < Z0B
c     then set CBC to CBCmax

         Z0B    = .01              ! Roughness parameter (in m)
         CBCMIN = .0025            ! min. Drag Coef. Bottom log layer
                                   ! not resolved.
         CBCmax = 1.000            ! max. Drag Coef
c -----------------------------------------------------------------
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:                                                             !
!     Define CBC [see O2006, Appendix A]:                              !
!     Also see Mellor,JPO,2002:Osc.Blayer...                           !
!                                                                      !
      vonk=0.4; vonk2=vonk*vonk
      ZSH=Z0B*float(nwad)    !ZSH=0 is used by original POM w/o WAD
      WRITE(6,'('' cbcmin,vonk,vonk2,z0b,zsh ='',1p5e13.4,/)')
     &             cbcmin,vonk,vonk2,z0b,zsh
         do j = 1, JM
           do i = 1, IM
c--- original pom98
c           CBC(I,J) = (KAPPA/LOG((1.+ZZ(KBM1))*H(I,J)/Z0B))**2 
c           CBC(I,J) = MAX(CBCmin,CBC(I,J))  
c           CBC(I,J)=  MIN(CBCmax,CBC(I,J))
c--- wad -----------
            CBC(I,J)=  CBCMIN
      CBC(I,J)=MAX(CBCMIN,vonk2/LOG((ZSH+(ZZ(KBM1)-Z(KB))*H(I,J))
     1        /Z0B)**2)*FSM(I,J)
          end do                          
         end do                          
C
        CALL PRXY('  TOPOGRAPHY 2 ',TIME, H,IP, 1 ,JM, 1 ,10.)
        CALL PRXY(' DX             ',TIME,DX ,IP,ISK,JM,JSK,1.)
        CALL PRXY(' DY             ',TIME,DY ,IP,ISK,JM,JSK,1.)
        CALL PRXY(' FSM            ',TIME,FSM,IP,ISK,JM,JSK,1.)
        CALL PRXY(' DUM            ',TIME,DUM,IP,ISK,JM,JSK,1.)
        CALL PRXY(' DVM            ',TIME,DVM,IP,ISK,JM,JSK,1.)
c       CALL PRXY(' CBC            ',TIME,CBC,IP,ISK,JM,JSK,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
!lyo:!wad:add wetmask & wmarsh (in case it's time-dependent)           !
     & ,wetmask,wmarsh
! 
  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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Update CBC:                                                  !
!                                                                      !
      DO J=1,JM; DO I=1,IM
      CBC(I,J)=MAX(CBCMIN,vonk2/LOG((ZSH+(ZZ(KBM1)-Z(KB))*D(I,J))
     1        /Z0B)**2)*FSM(I,J)
      enddo; enddo
      CALL PRXY(' CBC           ',TIME,CBC,IM,1,JM,1,0.E0)
!                                                                      !
!----------------------------------------------------------------------!
! 
      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)
c     CALL PRXZ('TB' ,TIME,TB,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c     CALL PRXZ('RHO  ' ,TIME,RHO  ,IM,ISK,JM,10,25,KB,0.,DT,ZZ)
c     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)
c     CALL PRXY(' DRX2D      ',TIME,DRX2D,IP,ISK,JM,JSK,0.)
c     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
      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
      WUSURF(I,J)=RAMP*2.E-4
     1         *.25*(DVM(I,J+1)+DVM(I-1,J+1)+DVM(I-1,J)+DVM(I,J))
      WVSURF(I,J)=0.
      WTSURF(I,J)=0.
      WSSURF(I,J)=0.
   85 SWRAD(I,J)=0.
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      
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Water depth can be solved by Smolarkiewicz to prevent oscil. !
!                                                                      !
      if (nsmolar.eq.1) then
      dd0(:,:)=h(:,:)+ el(:,:)*fsm(:,:)
      ddb(:,:)=h(:,:)+elb(:,:)*fsm(:,:)
      nitera=3; sw=1.0; tps(:,:)=0.0           !aam2d(:,:)
      call advt2d(ddb,dd0,ddf,nitera,sw,dx,dy,ua,va,art,aru,arv,
     1    fsm,dum,dvm,tps,dte2)
      elf(:,:)=ddf(:,:)-h(:,:)
      else
      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)
      endif
!                                                                      !
!----------------------------------------------------------------------!
C
      CALL BCOND(1)
C
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Check wet/dry on ELF:                                        !
!                                                                      !
      if (nwad.eq.1) then
      do j=1,jm; do i=1,im
      DWET(I,J)=H(I,J)+ELF(I,J)*fsm(i,j)
      enddo; enddo
      do j=1,jm; do i=1,im
      wetmask(i,j)=fsm(i,j)
      if (dwet(i,j).le.hc) wetmask(i,j)=0.0
      enddo; enddo
      endif
!                                                                      !
!----------------------------------------------------------------------!
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     
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:alpha=0 for WAD because elf may be bad at this point!!       !
!     A better way may be found??                                      !
!                                                                      !
      ALPHA=0.225*float(1-nwad)
!                                                                      !
!----------------------------------------------------------------------!
      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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Check wet/dry on UAF & VAF:                                  !
!                                                                      !
      if (nwad.eq.1) then
      do j=1,jm; do i=2,im
      if (0.5*(dwet(i,j)+dwet(i-1,j)).le.hc) then
      uaf(i,j)=0.0
      else                      !Set flux=0 if coming from dry cell:
      if (wetmask(i-1,j).eq.0.0 .and. uaf(i,j).gt.0.0) uaf(i,j)=0.0
      if (wetmask(i,j)  .eq.0.0 .and. uaf(i,j).lt.0.0) uaf(i,j)=0.0
      endif
      enddo; enddo
!
      do j=2,jm; do i=1,im
      if (0.5*(dwet(i,j)+dwet(i,j-1)).le.hc) then
      vaf(i,j)=0.0
      else                      !Set flux=0 if coming from dry cell:
      if (wetmask(i,j-1).eq.0.0 .and. vaf(i,j).gt.0.0) vaf(i,j)=0.0
      if (wetmask(i,j)  .eq.0.0 .and. vaf(i,j).lt.0.0) vaf(i,j)=0.0
      endif
      enddo; enddo
      endif
!                                                                      !
!----------------------------------------------------------------------!
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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Update CBC:                                                  !
!                                                                      !
      DO J=1,JM; DO I=1,IM
      CBC(I,J)=MAX(CBCMIN,vonk2/LOG((ZSH+(ZZ(KBM1)-Z(KB))*D(I,J))
     1        /Z0B)**2)*FSM(I,J)
      enddo; enddo
!                                                                      !
!----------------------------------------------------------------------!
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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:The followings include Oey's [1996] and Oey & Chen's [1992]  !
!     implementations of river and temperature surface conditions      !
!     They are not directly for WAD, but ssurfwet is required below    !
!                                                                      !
      do j=1,jm
      do i=1,im
      ssurfwet(i,j)=SSURF(i,j)*(1.-WRIV(i,j)/(WRIV(i,j)+1.E-28))
      enddo
      enddo
c
      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)
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Prepare dry-cell T/S "initial conditions" for next time-step !
!                                                                      !
      if (nwad.eq.1) then
      do k=1,kb-1; do j=1,jm; do i=1,im
      if (wetmask(i,j).eq.0.0) then
      uf(i,j,k)=cwetrlx1*tb(i,j,k)+cwetrlx2*tsurf(i,j)
      vf(i,j,k)=cwetrlx1*sb(i,j,k)+cwetrlx2*ssurfwet(i,j)
      endif
      enddo; enddo; enddo
      endif
!----------------------------------------------------------------------!
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)
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Check wet/dry on UF & VF:                                    !
!                                                                      !
      if (nwad.eq.1) then
      do k=1,kb-1
      do j=1,jm; do i=2,im
      if (wetmask(i,j)*wetmask(i-1,j).eq.0.0) uf(i,j,k)=0.0
      enddo; enddo
      do j=2,jm; do i=1,im
      if (wetmask(i,j)*wetmask(i,j-1).eq.0.0) vf(i,j,k)=0.0
      enddo; enddo
      enddo
      endif
!----------------------------------------------------------------------!
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
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)
c       CALL PRXY(' DRX2D      ',TIME,DRX2D,IP,ISK,JM,JSK,0.E0)
c       CALL PRXY(' DRY2D      ',TIME,DRY2D,IP,ISK,JM,JSK,0.E0)
        CALL PRXY(' WETMASK ',TIME, WETMASK(1,1),IM,1,JM,1,1.0)
!tne:!wad: to print wet cells & surf. elev. relative to MSL
      do j=1,jm; do i=1,im
         tps(i,j)=fsm(i,j)-wetmask(i,j)
      enddo; enddo
        CALL PRXY('FSM-WETMASK',TIME,tps(1,1),IM,1,JM,1,1.0)
      do j=1,jm; do i=1,im
         tps(i,j)=(elb(i,j)+hhi)*wetmask(i,j)
      enddo; enddo
        CALL PRXY(' ELB+HHI (MSL)',TIME, tps ,IM,1,JM,1,0.E0)
c
ctne: -------------- save for matlab plot
c
       write(8,'(2i5,4f8.3)')im-2,jm-2,time*24.,99.9,99.9,99.9
      do j=2,jm-1; do i=2,im-1
       write(8,'(2i5,4f8.3)')i,j,tps(i,j),uab(i,j),vab(i,j),d(i,j)
      enddo; enddo
c
      do j=1,jm; do i=1,im
         tps(i,j)=d(i,j)*fsm(i,j)
      enddo; enddo
        CALL PRXY(' WATER DEPTH (D)',TIME, tps ,IM,1,JM,1,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)
c         CALL PRXZ('SB' ,TIME,SB,IM,ISK,JM,10,25,KB,1.,DT,ZZ)
c         CALL PRXZ('RHO',TIME,RHO,IM,ISK,JM,10,25,KB,0.0,DT,ZZ)
c         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)
c         CALL PRXZ('L ' ,TIME,L ,IM,ISK,JM,10,25,KB,0.0,DT,Z)
c         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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Sum over wet areas only:                                     !
!                                                                      !
      DAREA=DX(I,J)*DY(I,J)*WETMASK(I,J) !lyo:!wad:
!     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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Sum over wet areas only:                                     !
!                                                                      !
      DAREA=DX(I,J)*DY(I,J)*WETMASK(I,J) !lyo:!wad:
!     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
!lyo:!wad:add wetmask & wmarsh (in case it's time-dependent)           !
c    & ,wetmask,wmarsh
! 
C
!lyo:!wad:print final successful completion:
      WRITE(6,602) TIME
  602 FORMAT(/2X,'JOB SUCCESSFULLY COMPLT.; TIME = ',1P1E13.5,' DAYS')
!
      STOP
      END
C
      SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE)
       INCLUDE 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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/,SMALL/1.E-8/
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
 20   CONTINUE
C
!tne:!wad:
C Very large (9m) tide-like BC (imposed by tidal inflow vel. not elev.) 
C      note that ELW already include hhi (see SEAMOUNT)
!lyo:!wad:Specify time-dependent tidal tidal amplitude w/period=1day,
!         & amplitude=9m; note that tide is to be added to the geostro-
!         phically balanced ELW below.  So ETIDE is referenced wrt MSL:
!         [but perhaps it will be more direct to specify elevation??]
      ETIDE=-9.*SIN(2.*PI*TIME/1.0)
C
C---------- VELOCITY --------------        
C--- Governing external B.C.s -------------
      DO 210 J=2,JMM1
C            EAST
!tne:!wad:- for wad replace H with D
      UAF(IM,J)=UABE(J)+RFE*SQRT(GRAV/D(IMM1,J))*(EL(IMM1,J)-ELE(J))
      UAF(IM,J)=RAMP*UAF(IM,J)
      VAF(IM,J)=0.0
C            WEST
!tne:!wad:- for wad replace H with D add tide
      UAF(2,J)=UABW(J)-RFW*SQRT(GRAV/D(2,J))*(EL(2,J)-ELW(J)-ETIDE)
      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
!tne:!wad:- for wad replace H with D
      VAF(I,JM)=VABN(I)+RFN*SQRT(GRAV/D(I,JMM1))*(EL(I,JMM1)-ELN(I))
      VAF(I,JM)=RAMP*VAF(I,JM)
      UAF(I,JM)=0.0
C            SOUTH
!tne:!wad:- for wad replace H with D
      VAF(I,2)=VABS(I)-RFS*SQRT(GRAV/D(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
!tne:!wad:: test wad with shallower domain
c     HMAX=4500.
      HMAX= 50.
      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 605 K=1,KB
      DO 605 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)=Q2(IM,J,K)-U1*(SMALL-Q2(IM,J,K)) 
        VF(IM,J,K)=Q2L(IM,J,K)-U1*(SMALL-Q2L(IM,J,K)) 
      ELSE
        UF(IM,J,K)=Q2(IM,J,K)-U1*(Q2(IM,J,K)-Q2(IMM1,J,K))
        VF(IM,J,K)=Q2L(IM,J,K)-U1*(Q2L(IM,J,K)-Q2L(IMM1,J,K))
      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)=Q2(1,J,K)-U1*(Q2(1,J,K)-SMALL) 
        VF(1,J,K)=Q2L(1,J,K)-U1*(Q2L(1,J,K)-SMALL)       
      ELSE
        UF(1,J,K)=Q2(1,J,K)-U1*(Q2(2,J,K)-Q2(1,J,K))
        VF(1,J,K)=Q2L(1,J,K)-U1*(Q2L(2,J,K)-Q2L(1,J,K))
      ENDIF
 605  CONTINUE
      DO 610 K=1,KB
      DO 610 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)=Q2(I,JM,K)-U1*(SMALL-Q2(I,JM,K)) 
        VF(I,JM,K)=Q2L(I,JM,K)-U1*(SMALL-Q2L(I,JM,K)) 
      ELSE
        UF(I,JM,K)=Q2(I,JM,K)-U1*(Q2(I,JM,K)-Q2(I,JMM1,K))
        VF(I,JM,K)=Q2L(I,JM,K)-U1*(Q2L(I,JM,K)-Q2L(I,JMM1,K))
      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)=Q2(I,1,K)-U1*(Q2(I,1,K)-SMALL)    
        VF(I,1,K)=Q2L(I,1,K)-U1*(Q2L(I,1,K)-SMALL)        
      ELSE
        UF(I,1,K)=Q2(I,1,K)-U1*(Q2(I,2,K)-Q2(I,1,K))
        VF(I,1,K)=Q2L(I,1,K)-U1*(Q2L(I,2,K)-Q2L(I,1,K))
      ENDIF
 610  CONTINUE
      DO 615 K=1,KB
      DO 615 J=1,JM
      DO 615 I=1,IM
      UF(I,J,K)=UF(I,J,K)*FSM(I,J)+1.E-10
      VF(I,J,K)=VF(I,J,K)*FSM(I,J)+1.E-10
 615  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 'comblk98_wad.h'   
C
      DATA PI/3.14167/
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 'comblk98_wad.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- 1000.)/RHOREF
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*RHOREF*ZZ(K)*H(I,J)*1.E-5 
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.1+.0821*P+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/RHOREF*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=9
      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 'comblk98_wad.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 'comblk98_wad.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),STF(IM,JM,KB)
      EQUIVALENCE (A,SM),(C,SH),(PROD,KN),(TPS,DH)
      EQUIVALENCE (DTEF,CC)
      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*RHOREF*ZZ(K)*H(I,J)*1.E-4
      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)*H(I,J)))
     1 +GRAV**2*2./(CC(I,J,K-1)**2+CC(I,J,K)**2)
  102 CONTINUE
      DO 220 K=2,KBM1
      DO 220 J=1,JM
      DO 220 I=1,IM
      L(I,J,K)=Q2LB(I,J,K)/Q2B(I,J,K)
      GH(I,J,K)=L(I,J,K)**2/Q2B(I,J,K)*BOYGR(I,J,K)              
      GH(I,J,K)=MIN(GH(I,J,K),.028)
  220 CONTINUE
      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.  
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)
      GHC=-6.0
      DO 110 K=1,KB
      DO 110 J=1,JM
      DO 110 I=1,IM
      STF(I,J,K)=1.
      IF(GH(I,J,K).LT.0.) STF(I,J,K)=1.0-0.9*(GH(I,J,K)/GHC)**1.5
      IF(GH(I,J,K).LT.GHC) STF(I,J,K)=0.1
  110 DTEF(I,J,K)=Q2B(I,J,K)*SQRT(Q2B(I,J,K))/(B1*Q2LB(I,J,K)+SMALL)
     1           *STF(I,J,K)
      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***********************************************************************
      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 230 K=1,KB
      DO 230 J=1,JM
      DO 230 I=1,IM
      COEF1=A2*(1.-6.*A1/B1*STF(I,J,K))
      COEF2=3.*A2*B2/STF(I,J,K)+18.*A1*A2
      COEF3=A1*(1.-3.*C1-6.*A1/B1*STF(I,J,K))
      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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:See O2006, Appendix A:                                       !
!                                                                      !
      KN(I,J,K)=(vonk*zsh+L(I,J,K))*SQRT(ABS(Q2(I,J,K)))
!     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 'comblk98_wad.h'
      DIMENSION F(IM,JM,KB),WFSURF(IM,JM),FSURF(IM,JM),DH(IM,JM)
      DIMENSION SWRAD(IM,JM),RAD(IM,JM,KB),R(5),AD1(5),AD2(5)  
      EQUIVALENCE (TPS,DH)
C
C Irradiance parameters after Paulson and Simpson, JPO, 1977, 952-956.
C
      NTP=2
C       NTP         =     1      2       3       4       5
C   JERLOV TYPE     =     I      IA      IB      II      III
      DATA R   /        .58 ,   .62  ,  .67  ,  .77  ,  .78   /
      DATA AD1 /        .35 ,   .60  ,  1.0  ,  1.5  ,  1.4   /
      DATA AD2 /        23. ,   20.  ,  17.  ,  14.  ,  7.9   /
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 PENETRATES WATER COLUMN
C   NBC=3; SURF. B.C. IS TSURF. NO SW RADIATIVE PENETRATION.
C   NBC=4; SURF. B.C. IS TSURF. 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
        RAD(I,J,K)=SWRAD(I,J)
     1     *(   R(NTP)*EXP(Z(K)*DH(I,J)/AD1(NTP))
     2             +(1.-R(NTP))*EXP(Z(K)*DH(I,J)/AD2(NTP))   )
  511   CONTINUE
      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) 
     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 'comblk98_wad.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 'comblk98_wad.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 'comblk98_wad.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.
!lyo:!wad:Use finer grid for nsmolar=1
      DELX=8.0E3  !=4.E3
      nvar=1  !=1 for variable grid; =0 otherwise
      if (nvar.eq.0) delx=delx*0.5
!
c     DELH=0.9      ! seamount
      DELH=1.15     ! island (for wad)
c     RA=25.E3      ! narrow
      RA=50.E3      ! wideer (for wad)
c     H0=4500.      ! deep open ocean
      H0= 50.       ! shallow (for wad) - update HMAX in BCOND
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     
C  Variable DX, DY.
      DX(I,J)=DELX-DELX/2.*SIN(PI*FLOAT(I)/FLOAT(IM))*float(nvar)
      DY(I,J)=DELX-DELX/2.*SIN(PI*FLOAT(J)/FLOAT(JM))*float(nvar)
      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
C
C     Define depth:
C
!tne:!wad:!lyo:!wad:
!     Note that unlike original seamount, H<0 is now water below MSL, and 
!     H>0 is above the MSL and is either land (if h>hland, see below)
!     or is potential wet-and-dry region
!     (MAIN will reverse these to make "h" comparible w/POM)
!
      hland=5.e0 !note all pnts are potential WAD if we set
                 !hland >= -h0*(1.e0-delh) (=7.5m)
!
      DO I=1,IM
      DO J=1,JM
      H(I,J)=(-1.)*H0*(1.-DELH*EXP(-
     1  ((X(I)-X((IM+1)/2))**2+(Y(J)-Y((JM+1)/2))**2)/RA**2) )
c
!lyo:!wad:!slpmin: CALL SLPMIN, here before WAD redefinitions
!
!lyo:!wad:Make region near top of seamount absolute land region:
      if(h(i,j).gt.hland) h(i,j)=hhi+1.e0
!
      ENDDO
!tne:!wad: Close the north and south boundaries to form a channel.
      H(I,1)=1.+hhi
      H(I,JM)=1.+hhi
      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
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
!lyo:wad:DT is not used, so delete/comment-out the following line:
!lyo:wad:     DT(I,J)=abs(H(I,J))  !tne:!wad:
      FSM(I,J)=0.
      DUM(I,J)=0.
      DVM(I,J)=0.
!lyo:wad:we define absolute land grid-cells fsm=0.0 above and also 
!        potential WAD cells fsm=1.0 below to "go-along" with 
!        this old seamount subroutine; the correct masks together 
!        with wet/dry mask will be defined later in MAIN
!lyo:wad:      IF(abs(H(I,J)).GT.1.) FSM(I,J)=1.  !tne:!wad:
      IF(H(I,J).LT.hhi) 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
!                                                                      !
!----------------------------------------------------------------------!
!lyo:!wad:Note that if it is necessary to adjust "H" with SLPMIN       !
!     then better do the adjustments above right after the bathymetry  !
!     is defined (search backward !slpmin:), BUT before the subsequent !
!     WAD definitions                                                  !
!----------------------------------------------------------------------!
!                                                                      !
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
!lyo:wad:
!     TB(I,J,K)=5.+15.*EXP(ZZ(K)*H(I,J)/1000.)-TBIAS !original code
      if (h(i,j).le.0.0) then !lyo:wad:stratified water below MSL:
      TB(I,J,K)=5.+15.*EXP(ZZ(K)*(-H(I,J))/1000.)
      else                    !lyo:wad:well-mixed water above MSL:
      TB(I,J,K)=5.+15.
      endif
      TB(I,J,K)=TB(I,J,K)-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
      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.  ! =0 for no radiation BC only UA
      RFW=1.
      RFN=0.
      RFS=0.
      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
!tne:!wad:- for wad
      DO 62 I=1,IM
       ELN(I)=0.-hhi
       ELS(I)=0.-hhi
  62  CONTINUE
!tne:!wad:- for wad
      DO 63 J=1,JM
       ELE(J)=(ELE(J)-hhi)*FSM(IM,J)
       ELW(J)=(ELW(J)-hhi)*FSM(2,J)
  63  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 'comblk98_wad.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
!----------------------------------------------------------------------!
!lyo:!wad:Additional subroutines used by WAD:                          !
!                                                                      !
cwetdry391d_dam02:2-d version of Smolarkiewicz from /home/lyo/pom/pom2k/
c     pom2k.f's subroutine advt2
c     This version gets rid of common block [i.e. the routine is now
c     self-contained] and also 2-d calculation [solve for D]:
c
c     d(D)/dt + d(UD)/dx + d(VD)/dy = Diffusion_of_(D)
c
cwetdry391d_dam02:
c--   subroutine advt2(fb,f,fclim,ff,xflux,yflux,nitera,sw)
      subroutine advt2d(fb,f,ff,nitera,sw,dx,dy,u,v,art,aru,arv,
     1    fsm,dum,dvm,aam,dti2)
c    1    fsm,dum,dvm,aam,dti2,im,jm)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Integrates conservative scalar equations.           *
C *                                                                    *
C *                This is a first-order upstream scheme, which        *
C *                reduces implicit diffusion using the Smolarkiewicz  *
C *                iterative upstream scheme with an antidiffusive     *
C *                velocity.                                           *
C *                                                                    *
C *                It is based on the subroutines of Gianmaria Sannino *
C *                (Inter-university Computing Consortium, Rome, Italy)*
C *                and Vincenzo Artale (Italian National Agency for    *
C *                New Technology and Environment, Rome, Italy),       *
C *                downloaded from the POM FTP site on 1 Nov. 2001.    *
C *                The calculations have been simplified by removing   *
C *                the shock switch option. It should be noted that    *
C *                this implementation does not include cross-terms    *
C *                which are in the original formulation.              *
C *                                                                    *
C *                fb,f,fclim,ff . as used in subroutine advt1         *
C *                xflux,yflux ... working arrays used to save memory  *
C *                nitera ........ number of iterations. This should   *
C *                                be in the range 1 - 4. 1 is         *
C *                                standard upstream differencing;     *
C *                                3 adds 50% CPU time to POM.         *
C *                sw ............ smoothing parameter. This should    *
C *                                preferably be 1, but 0 < sw < 1     *
C *                                gives smoother solutions with less  *
C *                                overshoot when nitera > 1.          *
C *                                                                    *
C *                Reference:                                          *
C *                                                                    *
C *                Smolarkiewicz, P.K.; A fully multidimensional       *
C *                  positive definite advection transport algorithm   *
C *                  with small implicit diffusion, Journal of         *
C *                  Computational Physics, 54, 325-362, 1984.         *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      real sw,dti2
      integer nitera,im,jm
c
!lyo:!wad:Use finer grid for nsmolar=1
      PARAMETER (IM=65,JM=49)         !dx=0.1
!     PARAMETER (IM=131,JM=99)                                    
c
      real fb(im,jm),f(im,jm),ff(im,jm)
      real dx(im,jm),dy(im,jm),u(im,jm),v(im,jm),aam(im,jm)
      real fsm(im,jm),dum(im,jm),dvm(im,jm)
      real art(im,jm),aru(im,jm),arv(im,jm)

c     integer ix,jx
c     PARAMETER (IX=62,JX=18)
      real xflux(im,jm),yflux(im,jm)
      real fbmem(im,jm)
      real xmassflux(im,jm),ymassflux(im,jm)
      integer i,j,k,itera,imm1,jmm1
C
C     Calculate horizontal mass fluxes:
C
      imm1=im-1; jmm1=jm-1
c
        do j=1,jm
          do i=1,im
            xmassflux(i,j)=0.e0
            ymassflux(i,j)=0.e0
          end do
        end do
C
        do j=2,jmm1
          do i=2,im
            xmassflux(i,j)=0.5e0*(dy(i-1,j)+dy(i,j))*u(i,j)
          end do
        end do
C
        do j=2,jm
          do i=2,imm1
            ymassflux(i,j)=0.5e0*(dx(i,j-1)+dx(i,j))*v(i,j)
          end do
        end do
C
        do j=1,jm
          do i=1,im
            fbmem(i,j)=fb(i,j)
          end do
        end do
C
C     Start Smolarkiewicz scheme:
C
      do itera=1,nitera
C
C     Upwind advection scheme:
C
          do j=2,jm
            do i=2,im
              xflux(i,j)=0.5e0
     $                      *((xmassflux(i,j)+abs(xmassflux(i,j)))
     $                        *fbmem(i-1,j)+
     $                        (xmassflux(i,j)-abs(xmassflux(i,j)))
     $                        *fbmem(i,j))
C
              yflux(i,j)=0.5e0
     $                      *((ymassflux(i,j)+abs(ymassflux(i,j)))
     $                        *fbmem(i,j-1)+
     $                        (ymassflux(i,j)-abs(ymassflux(i,j)))
     $                        *fbmem(i,j))
            end do
          end do
C
C     Add net advective fluxes and step forward in time:
C
          do j=2,jmm1
            do i=2,imm1
              ff(i,j)=xflux(i+1,j)-xflux(i,j)
     $                 +yflux(i,j+1)-yflux(i,j)
              ff(i,j)=(fbmem(i,j)*art(i,j)
     $                   -dti2*ff(i,j))/(art(i,j))
            end do
          end do
C
C     Calculate antidiffusion velocity:
C
      call smol_adif(xmassflux,ymassflux,ff,sw,fsm,aru,arv,dti2)
c     call smol_adif(xmassflux,ymassflux,ff,sw,fsm,aru,arv,dti2,im,jm)
C
        do j=1,jm
          do i=1,im
              fbmem(i,j)=ff(i,j)
          end do

        end do
C
C     End of Smolarkiewicz scheme
C
      end do
C
C     Add horizontal diffusive fluxes:
C
        do j=2,jm
          do i=2,im
            xmassflux(i,j)=0.5e0*(aam(i,j)+aam(i-1,j))
            ymassflux(i,j)=0.5e0*(aam(i,j)+aam(i,j-1))
          end do
        end do
C
        do j=2,jm
          do i=2,im
           xflux(i,j)=-xmassflux(i,j)
     $                   *(fb(i,j)-fb(i-1,j))*dum(i,j)
     $                   *(dy(i,j)+dy(i-1,j))/(dx(i,j)+dx(i-1,j))
           yflux(i,j)=-ymassflux(i,j)
     $                   *(fb(i,j)-fb(i,j-1))*dvm(i,j)
     $                   *(dx(i,j)+dx(i,j-1))/(dy(i,j)+dy(i,j-1))
          end do
        end do
C
C     Add net horizontal fluxes and step forward in time:
C
        do j=2,jmm1
          do i=2,imm1
            ff(i,j)=ff(i,j)-dti2*(xflux(i+1,j)-xflux(i,j)
     $                               +yflux(i,j+1)-yflux(i,j))
     $                           /(art(i,j))
          end do
        end do
C
      return
C
      end
C
C
      subroutine smol_adif(xmassflux,ymassflux,ff,sw,fsm,aru,arv,dti2)
c     subroutine smol_adif(xmassflux,ymassflux,ff,sw,fsm,aru,arv,dti2,
c    1                     im,jm)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Calculates the antidiffusive velocity used to       *
C *                reduce the numerical diffusion associated with the  *
C *                upstream differencing scheme.                       *
C *                                                                    *
C *                This is based on a subroutine of Gianmaria Sannino  *
C *                (Inter-university Computing Consortium, Rome, Italy)*
C *                and Vincenzo Artale (Italian National Agency for    *
C *                New Technology and Environment, Rome, Italy),       *
C *                downloaded from the POM FTP site on 1 Nov. 2001.    *
C *                The calculations have been simplified by removing   *
C *                the shock switch option.                            *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      integer im,jm
c
!lyo:!wad:Use finer grid for nsmolar=1
      PARAMETER (IM=65,JM=49)         !dx=0.1
!     PARAMETER (IM=131,JM=99)                                    
c
      real ff(im,jm)
      real xmassflux(im,jm),ymassflux(im,jm)
      real fsm(im,jm),aru(im,jm),arv(im,jm)
      real sw,dti2
      real mol,abs_1,abs_2
      real value_min,epsilon
      real udx,u2dt,vdy,v2dt,wdz,w2dt
      integer i,j,k,imm1,jmm1
C
      parameter (value_min=1.e-9,epsilon=1.0e-14)
C
c
      imm1=im-1; jmm1=jm-1
c
C     Apply temperature and salinity mask:
C
        do i=1,im
          do j=1,jm
            ff(i,j)=ff(i,j)*fsm(i,j)
          end do
        end do
C
C     Recalculate mass fluxes with antidiffusion velocity:
C
        do j=2,jmm1
          do i=2,im
            if(ff(i,j).lt.value_min.or.
     $         ff(i-1,j).lt.value_min) then
              xmassflux(i,j)=0.e0
            else
              udx=abs(xmassflux(i,j))
              u2dt=dti2*xmassflux(i,j)*xmassflux(i,j)
     $              /(aru(i,j))
              mol=(ff(i,j)-ff(i-1,j))
     $             /(ff(i-1,j)+ff(i,j)+epsilon)
              xmassflux(i,j)=(udx-u2dt)*mol*sw

              abs_1=abs(udx)
              abs_2=abs(u2dt)
              if(abs_1.lt.abs_2) xmassflux(i,j)=0.e0
            end if
          end do
        end do
C
        do j=2,jm
          do i=2,imm1
            if(ff(i,j).lt.value_min.or.
     $         ff(i,j-1).lt.value_min) then
              ymassflux(i,j)=0.e0
            else
             vdy=abs(ymassflux(i,j))
             v2dt=dti2*ymassflux(i,j)*ymassflux(i,j)
     $             /(arv(i,j))
             mol=(ff(i,j)-ff(i,j-1))
     $            /(ff(i,j-1)+ff(i,j)+epsilon)
             ymassflux(i,j)=(vdy-v2dt)*mol*sw
             abs_1=abs(vdy)
             abs_2=abs(v2dt)
             if(abs_1.lt.abs_2) ymassflux(i,j)=0.e0
            end if
          end do
        end do
C
      return
C
      end
!                                                                      !
!lyo:!wad:END of Additional subroutines used by WAD.                   !
!----------------------------------------------------------------------!
!                                                                      !