program main
      include 'wav.c'   
      dimension cbc(im,jm),wubot(im,jm),wvbot(im,jm)
      dimension elips(im,jm)
      data pi/3.141592654/,gee/9.807/
        rewind 10
        open(10, file='specavs',form='formatted')
C---------------------------------------------------------------
C   This wave model is described in the paper:
C
C  A Wave Model for Coupling with Numerical Ocean Circulation Models
C             by George Mellor and Mark Donelan
C
C   Also see references to M03 and M05 cited in this paper.
C
C     This main program is a surogate for pom2k.f and can be
C   discarded when coupling with pom2k.f or another circulation
C   model. When combined with a circulation model one must include
C   the following subroutines: 
C   wave, botdrag, cur2wave, specavs, speeds, stress, stepmpd,
C   wvfcns, wave2cur and the functions: summ, summ1, sumk,fsprm
C   and mxloc. The file specavs needs to be opened.
C
C     The common variables in wav.c will have to be melded with 
C   pom2k.c. 
C     After the above is achieved, simply call wave and call botdrag
C   for wave affected bottom drag (cbc).
C     The wave propagation time step, dtw, must be supplied.
C   Hopefully it can be the same as the internal time step in
C   pom2k.f; otherwise, for example, subroutine wave will have to
C   be called twice with half the time step.
C     The print routines and subroutine depth can be expunged as 
C   they are duplicated in pom2k.f. For minumum wave info, print the 
C   significant wave height, Hs, and mean wave propagation angle,
C   thtav. However, other variables such as peak frequency, phase 
C   and group speeds and wave number may be of interest.
C     The surface wind stress is solved in subroutine stress
C   given the 10 meter wind speed u10x and u10y. For now, use the
C   relations, wusurf=-(tpx0+ttx0) and wvsurf=-(tpy0+tty0).
C   At a later time, we will differentiate between the two kinds
C   of surface stress as discribed in the aforementioned paper.
C     See further comments in subroutine wave2cur. 
C     In stand-alone mode, this program main solves for wave
C   properties in an elliptical closed basin.
C---------------------------------------------------------------
C Setup
      time=0.
      iend=4000               ! run to steady state 
      iprint=2000
      iend=40 
      iprint=20
      dxm=200
      dym=500
      dtw=dxm/15.     ! for U10=20
      dtw=dxm/10.     ! for U10=10
      dth=2.*pi/float(kseg)        ! wave angle increment
c     dtw=2.5
      write(6,'('' im,jm,mm ='', 3i5)') im,jm,mm
      write(6,'('' iprint ='',i5,''  iend ='',i5)') iprint,iend
      write(6,'('' time step (dtw) ='',f8.3)') dtw
      write(6,'('' grid step (dx) ='',f10.2)') dxm       
      write(6,'('' angle step (dth) ='',f9.4,''   No. of angle segments
     & (kseg) ='',i5)') dth,kseg
      call depth
C
C Set up elliptical basin
C
      a=155000.
      b=35000.
C do dx
      do j=1,jm
        dx(1,j)=200.
      enddo
      ico=40
      do j=1,jm
        do i=1,ico
          dx(i+1,j)=dx(i,j)*1.1
        enddo
        do i=ico,im-1
          dx(i+1,j)=dx(i,j)
        enddo
      enddo
C do dy
c       dy(1,jm/2)=447.9
        dy(1,jm/2)=203.9    
      do j=jm/2,jm-1  
        dy(1,j+1)=dy(1,j)*1.1
c       dy(1,j+1)=dy(1,j)*1.2
      enddo
      do j=1,jm/2 
        dy(1,j)=dy(1,jm+1-j)
      enddo
      do j=1,jm
        do i=1,im
          dy(i,j)=dy(1,j)
        enddo
      enddo
C end dx,dy
      call prxy('dx   ',time,dx   ,im,1,jm,1  ,0.)
      call prxy('dy   ',time,dy   ,im,1,jm,1  ,0.)
      write(6,'(''dx'')')
      write(6,'(10f10.1)') (dx(i,1),i=1,im)
      write(6,'(''dy'')')
      write(6,'(10f10.1)') (dy(1,j),j=1,jm)
      x(1)=0.
      x(2)=dx(2,1)/2.
      do i=3,im
        x(i)=x(i-1)+(dx(i-1,1)+dx(i,1))/2.
      enddo
      write(6,'(''x'')')
      write(6,'(10f10.1)') (x(i),i=1,im)
      y(1)=-b
      do j=2,jm
        y(j)=y(j-1)+(dy(1,j)+dy(1,j-1))*0.5
      enddo
      write(6,'(''y'')')
      write(6,'(10f10.1)') (y(j),j=1,jm)

C Caution: small d and subsequently kpD < ~3 can result in
C     large cth and instability unless dtw is suitably small

      do j=1,jm
        do i=1,im
          d(i,j)=10.
          elips(i,j)=((x(i)-a)/a)**2+(y(j)/b)**2
          fsm(i,j)=0.0
          if(elips(i,j).lt.1.0) fsm(i,j)=1.0
          enddo
        enddo
      call prxy('elips',time,elips,im,1,jm,1  ,0.)
c     call prxy('d ',time,d  ,im,1,jm,1  ,0.)
C
C  End basin setup
C

      do j=1,jm
        do i=1,im
          ageinv(i,j)=1.
          thtav(i,j)=angle*pi/180.
          wubot(i,j)=0.01
          wvbot(i,j)=0.01
          cbc(i,j)=0.
          do k=1,kb-1
            u(i,j,k)=0. 
            v(i,j,k)=0.      
          enddo
        enddo
      enddo
      do j=1,jm
        do i=1,im
          dum(i,j)=fsm(i,j)*fsm(i-1,j)
          dvm(i,j)=fsm(i,j)*fsm(i,j-1)
        enddo
      enddo
C Set parameters for bottom drag
      cbcmin=.0025
      cbcmax=1.000
      z0b=0.01
C
      call prxy('fsm  ',time,fsm  ,im,1,jm,1  ,0.)
c     call prxy('dum  ',time,dum  ,im,1,jm,1  ,0.)
c     call prxy('dvm  ',time,dvm  ,im,1,jm,1  ,0.)
C-----------------------------------------------------------------
C                   Begin time stepping                                  
C-----------------------------------------------------------------
c     iend=10
c     iprint=1
                   do 1000 iint=1,iend
      time=time+dtw/86400.
c     write(6,'(''time,iint,iend ='',f12.5,2i5)') time,iint,iend    
      do j=1,jm
        do i=1,im
          angle=120.
          u10(i,j)=10.
          u10x(i,j)=u10(i,j)*cos(pi*angle/180.)
          u10y(i,j)=u10(i,j)*sin(pi*angle/180.)
        enddo
      enddo

      call wave(dtw,dth,time)        
c     call botdrag (im,jm,fsm,
c    &     kp,cp,ent,wubot,wvbot,d,zz(kb-1),z0b,cbcmin,cbcmax,cbc)

 1000              continue
C--------------------------------------------------------------------
C                    End time stepping
C--------------------------------------------------------------------
      stop
      end
C-------------------------------------------------------------------

C---------------------------------------------------------------
      subroutine wave(dtw,dth,time)        
C------------------------------------------------------------
C Subroutine to calculate wave energy density as a function of
C x, y and theta. Grid is orthogonal curvilinear. For rectilinear
C grid, simply set dx(i,j) or dy(i,j) = constant
C
C  Input:
C    z(k),dz(k) = vertical coordinate and spacing (non-d.)
C    zz(k),dzz(k) = middle cell coordinate and spacing
C    dx(i,j), dy = spatial grid cell increments located at cell center
C    d(i,j) = water column depth (m)
C    fsm(i,j) = 1 for water cell, = 0 for land cell
C    u10x(i,j),u10y(i,j) = wind vector at 10 m height (m s-1)
C    u(i,j,k),v(i,j,k) = current vector (m s-1)
c    nitera = # of iterations in s.r. stepmpd

C  Output wave variables:
C    enb(i,j,m), en, enf = past, present,future energy densities
C       (m+3 s-2)
C    ent(i,j) en averaged over all wave directions
C    Hs(i,j) = sugnificant wave height = 4*sqrt(ent/gee) (m)
C    thtav(i,j) = average wave propagation direction (rad)

C  Internal wave variables (some could also be output)
C    theta(m) = wave directional angle located at cell edge (rad) 
C    cs(m),sn(m) = cosine and sine of theta
C    thetam(m) = wave directional angle located at cell center 
C    csm(m),snm(m) = cosine and sine of thetam
C    sigth(i,j,m)=theta dependent frequency
C    cth(i,j,m)=theta dependent phase speed
C    kth(i,j,m)= theta dependent wave number
C    cg(i,j,m) = spectrally averaged group speeds located at
C         cell center (m s-1)
C    sigp(i,j)=peak frequency
C    cp(i,j) = peak frequency phase speed (m s-1)
C    kp(i,j) = peak frequency wave number (m-1)
C    cgx(i,j,m), cgy(i,j,m) = propagation vector in x, y
C         directions located at cell edge (m s-1)
C    ud(i,j),vd(i,j) = Dopler velocities (m s-1)
C    cthth(i,j,m) = propagation speed in the theta direction 
C        located at cell edge (rad s-1)
C    ageinv(i,j)=u10*sigp/gee (= 0.83 for fully developed waves)
C    fspr(i,j)=spreading function for swin
C    beta = constant in fspr; higher beta = narrower spread
C        but recommended value = 2,2
C    swin(i,j,m) = wave energy source term
C    sdis(i,j,m) = surface wave energy dissipation 
C    bdis1(i,j,m) = bottom friction dissipation
C    bdis2(i,j,m) = bottom depth-induced dissipation
C    gama = empirical constant in formualtion of bdis2
C    Sradx(i,j),Srady(i,j),Sradxy(i,j),Srad(i,j) = wave radiation
C        stress related to wave velocity moments
C    Sradpx(i,j),Sradpy(i,j),Sradp(i,j) = wave radiation stress
C        related to wave pressure
C    Fcg(i,jm), Fct(i,j,m) = correction for cg and cth to  
C        account for spectral averaging 
C    Cin, swcon = constants in swin curve fit.
C    adis, bdis = empirically determined constants in sdis term
C    swfct = factor to reduce dissipation for swells.
C    dif = horizontal advective diffusion coeficient (see note
C        in x-y step.
C------------------------------------------------------------
      include 'wav.c' 
      integer init
      dimension               !  2-d variables
     &    xfluxtot(im,jm),yfluxtot(im,jm),swintot(im,jm),
     &    thtavdeg(im,jm), ustp(im,jm)
      dimension               !  3-d variables
     &    xflux(im,jm,mm),yflux(im,jm,mm),tflux(im,jm,mm),
     &    fspr(im,jm,mm),fsprout(im,jm,mm),
     &    uak(im,jm,mm)
      dimension usdif(im,jm),vsdif(im,jm),ud(im,jm),vd(im,jm)
      dimension Fcg(im,jm,mm),Fct(im,jm,mm)
      dimension tps(im,jm,mm),F1(kb),F2(kb),F3(kb),F4(kb),F5(kb),
     &      Sradx(im,jm),Sradxy(im,jm),Srady(im,jm),Srad(im,jm),
     &      Sradpx(im,jm),Sradpy(im,jm),Sradp(im,jm)
      dimension jo(2),mo(2)
      data nio/1/,io/12/,njo/2/,jo/10,30/,nmo/2/,mo/1,21/      
      data pi/3.141592654/,gee/9.807/,smoth/.05/
      data beta/2.2/                ! spreading const. 
      data Cin/370./                ! Sin const.
      data swcon/0.33/              ! Sin const.
      data adis/0.925/              ! Sdis const.
      data bdis/0.18e-4/            ! Sdis const.
      data swfct/0.2/               ! factor to decrease swell diss.
      data nitera/3/                ! # of iterations in MPDATA
      data dif/0.2/                 ! diffusion coeficient 
                                    ! in x-y step section
      data gama/0.70/               ! depth-induced breaking const.
      data init/0/
C------Initialize------------------------------------------------ 
                      if(init.eq.0) then
C----------------------------------------------------------------                                         
      iskp=1
      jskp=1
      dtw2=2.*dtw
      dthinv=1./dth                ! inverse angle increment
      write(6,'('' Cin ='',f8.1,''  a ='',f8.3,'' b ='',e12.4)')
     &      Cin,a,b
      write(6,'('' beta ='',f10.3)') beta
      write(6,'('' nitera ='',i5 )') nitera     
      write(6,'('' dif ='',f8.3)') dif
      do j=1,jm
        do i=1,im-1
          aren(i,j)=dx(i,j)*dy(i,j)
          areninv(i,j)=1./aren(i,j)
          ud(i,j)=0.
          vd(i,j)=0.
        enddo
      enddo
      d0=100.
      write(6,'(5x,''   m      theta       cs       sn      thetam'',
     &''      csm       snm'')')
      do m=1,mm
        thetam(m)=dth*float(m-mm/2)
        theta(m)=dth*float(m-mm/2)-0.5*dth
        cs(m)=cos(theta(m))
        sn(m)=sin(theta(m))
        csm(m)=cos(thetam(m))
        snm(m)=sin(thetam(m))
        write (6,'(5x,i5,6f10.3)') m, theta(m),cs(m),sn(m),
     &        thetam(m),csm(m),snm(m)
      enddo
      do j=1,jm    
        do i=1,im
          do m=1,mm
            cth(i,j,m)=10.
            kth(i,j,m)=.1
            kthD(i,j,m)=1.
          enddo
          ageinv(i,j)=1.
          sigp(i,j)=10.
          cp(i,j)=gee/sigp(i,j)
          kaveD(i,j)=sigp(i,j)/gee
          Sradx(i,j)=0.
          Sradxy(i,j)=0.
          Srady(i,j)=0.
          Srad(i,j)=0.
          Sradpx(i,j)=0.
          Sradpy(i,j)=0.
          Sradp(i,j)=0.
          ud(i,j)=0.
          vd(i,j)=0.
          tht_wnd(i,j)=0.
        enddo
      enddo
C
      do m=1,mm
        do j=1,jm
          do i=1,im
            sigthb(i,j,m)=sigp(i,j)
            sigth(i,j,m)=sigthb(i,j,m)
            sigthf(i,j,m)=sigth(i,j,m)
            cth(i,j,m)=gee/sigth(i,j,m)
            kth(i,j,m)=sigth(i,j,m)**2/gee
            kthD(i,j,m)=kth(i,j,m)*d(i,j)
            Hs(i,j)=1.95    
            Hs(i,j)=0.05   
            enb(i,j,m)=gee*(Hs(i,j)/4.)**2*fsprm(beta,m,mm,0.,thetam(m))
            enb(i,j,m)=max(.0001,enb(i,j,m))*fsm(i,j)
            en(i,j,m)=enb(i,j,m)
            enf(i,j,m)=en(i,j,m)
            xflux(i,j,m)=0.
            yflux(i,j,m)=0.
            tflux(i,j,m)=0.
            fspr(i,j,m)=fsprm(beta,m,mm,0.,thetam(m))
          enddo
        enddo
      enddo
      do j=1,jm
        do i=1,im
          ent(i,j)=summ(enb,dth,i,j,im,jm,mm)
        enddo
      enddo
C
C Open boundary conditions. (For a closed basin, not used)
C Add north south conditions as needed.
      Hsbc=2.          !Sig. wave height at boundaries
      Tp=10            !wave period at boundaries
      westang=0.       !western wave propagation angle; incoming
      eastang=pi       !eastern wave propagation angle; incoming
      do m=1,mm
        do j=1,jm
          enw(j,m)=gee*(Hsbc/4.)**2*fsprm(beta,m,mm,westang,thetam(m))
          enw(j,m)=max(.0001,enw(j,m))*fsm(1,j)
          ene(j,m)=gee*(Hsbc/4.)**2*fsprm(beta,m,mm,eastang,thetam(m))
          ene(j,m)=max(.0001,enw(j,m))*fsm(im,j) 
          sigthe(j,m)=2.*pi/Tp
          sigthw(j,m)=2.*pi/Tp
        enddo
      enddo
C
C
      call prxz('fspr',0.,fspr,im,iskp,jm,mm,jo,njo,0.,d,thetam)
      call prxz('enf',0.,enf,im,iskp,jm,mm,jo,njo,0.,d,thetam)
      call prxy('ent',0.,ent,im,1,jm,1  ,0.)
C---------------------------------------------------------------- 
                   init=1
                   endif
C---------------------------------------------------------------- 
      iskp=1
      jskp=1
      do j=1,jm
        do i=1,im
          u10(i,j)=sqrt(u10x(i,j)**2+u10y(i,j)**2)
          usdif(i,j)=U10x(i,j)-u(i,j,1)
          vsdif(i,j)=U10y(i,j)-v(i,j,1)
          Udif=sqrt(usdif(i,j)**2+vsdif(i,j)**2)
          tht_wnd(i,j)=0.
          if(Udif.ne.0.)
     &      tht_wnd(i,j)=acos(usdif(i,j)/(Udif))*sign(1.,vsdif(i,j))
        enddo
      enddo
      call stress (Hs,ageinv,usdif,vsdif,cd,tpx0,tpy0,
     &                          ttx0,tty0,ustw,im,jm)
C  Calculate sigp for wind driven portion of enf
      do i=2,im-1
        do j=2,jm-1
          if (fsm(i,j).gt.0.) then
            if(u10(i,j).gt.0.0) then
              mbig=mxloc(enf,fspr,i,j,im,jm,mm)
              sigp(i,j) =(gee/u10(i,j))
     &          *(gee*enf(i,j,mbig) /u10(i,j)**4/.0022)**(-.303)
            call wvfcns(sigp(i,j),d(i,j),kp(i,j),cp(i,j))   
            endif
            ageinv(i,j)=u10(i,j)*sigp(i,j)/gee       
          endif
        enddo
      enddo
C-----Tranfer sigp to sigth only for wind-driven portions of
C     en. Swell portions of en left unchanged.
      do m=1,mm
        do i=2,im-1
          do j=2,jm-1
            if(fspr(i,j,m).gt.0.10) sigth(i,j,m)=sigp(i,j)   
          enddo
        enddo
      enddo
C     call prxz('sigth',time,sigth,im,iskp,jm,mm,jo,njo,0.,d,thetam)
C     call prxy('u10      ',time,u10  ,im,1,jm,1  ,0. )
C---------------- Calculate propagation speeds-----------------------
      do j=1,jm
        do i=1,im  
          ent(i,j)=summ(en,dth,i,j,im,jm,mm)
          Hs(i,j)=0.5*(4.*sqrt(ent(i,j)/gee)+Hs(i,j))
C  Find kth and cth 
          do m=2,mm-1
            if(fsm(i,j).eq.1.) then
              call wvfcns(sigth(i,j,m),d(i,j),kth(i,j,m),cth(i,j,m))   
              kthD(i,j,m)=kth(i,j,m)*d(i,j)
C  Find spectrally averaged corrections to cg and cth
              call speeds(kthD(i,j,m),Fcg(i,j,m),Fct(i,j,m))
              cg(i,j,m)=Fcg(i,j,m)*cth(i,j,m)
            endif
          enddo
        enddo
      enddo
      do m=1,mm
        do j=2,jm-1  
          do i=2,im-1  
            if(fsm(i,j).eq.0.) then
              if(fsm(i+1,j).eq.1.) cg(i,j,m)=cg(i+1,j,m)
              if(fsm(i-1,j).eq.1.) cg(i,j,m)=cg(i-1,j,m)
              if(fsm(i,j+1).eq.1.) cg(i,j,m)=cg(i,j+1,m)
              if(fsm(i,j-1).eq.1.) cg(i,j,m)=cg(i,j-1,m)
            endif
          enddo
        enddo
      enddo
c     call prxz('enf ',time,enf ,im,iskp,jm,mm,jo,njo, 0.,d,thetam)
c     call prxy('ent',time,ent,im,1,jm,1  ,0. )
c     call prxy('sigp ',time,sigp ,im,1,jm,1  ,0. )
c     call prxy('kp ',time,kp ,im,1,jm,1  ,0. )
c     call prxy('cp ',time,cp ,im,1,jm,1  ,0. )
c     call prxy('cg ',time,cg ,im,1,jm,1  ,0. )
c     call prxz('Fcg ',time,Fcg ,im,iskp,jm,mm,jo,njo, 0.,d,thetam)
C ---------------------------------------------------------------
C   The solution is time split between x-y steps, theta steps and
C   source steps.
C ---------------------------------------------------------------
C                     x-y step fo wave energy
C ---------------------------------------------------------------
C  If dif=0.2, then diff = 0.2 everywhere, but, on open water
C  cells imediately next to land cells, diff = 0.5 for local upwind
C  B.C. If dif=0.5, then diff=0.5 everywhere.
C  N.B. not a great difference if diff = 0.5 everywhere.

      do m=1,mm
        do j=2,jm
          do i=2,im
            cgx(i,j,m)=0.5*(cg(i,j,m)+cg(i-1,j,m))*csm(m)    
     &                +0.5*(ud(i,j)+ud(i-1,j))
            cgy(i,j,m)=0.5*(cg(i,j,m)+cg(i,j-1,m))*snm(m)   
     &                +0.5*(vd(i,j)+vd(i,j-1))
            diff=0.5*(1.-dum(i,j))+dif*dum(i,j)
            xflux(i,j,m)= 0.5*cgx(i,j,m)*(en(i,j,m)+en(i-1,j,m))
     &           -diff*abs(cgx(i,j,m))*(enb(i,j,m)-enb(i-1,j,m))
            diff=0.5*(1.-dvm(i,j))+dif*dvm(i,j)
            yflux(i,j,m)= 0.5*cgy(i,j,m)*(en(i,j,m)+en(i,j-1,m))
     &           -diff*abs(cgy(i,j,m))*(enb(i,j,m)-enb(i,j-1,m))
            xflux(i,j,m)=0.5*(dy(i,j)+dy(i-1,j))*xflux(i,j,m)
            yflux(i,j,m)=0.5*(dx(i,j)+dx(i,j-1))*yflux(i,j,m)
          enddo
        enddo
      enddo
c     call pryz('cgy  ',time,cgy  ,im,jm,jskp,mm,12,1,0.,d,thetam)
c     call pryz('en   ',time,en   ,im,jm,jskp,mm,12,1,0.,d,thetam)
c     call pryz('enb  ',time,enb  ,im,jm,jskp,mm,12,1,0.,d,thetam)
c     call pryz('yflux',time,yflux,im,jm,jskp,mm,12,1,0.,d,thetam)
      do m=1,mm-1
        do j=2,jm-1
          do i=2,im-1
            enf(i,j,m)=enb(i,j,m)
     &        -(xflux(i+1,j,m)-xflux(i,j,m)
     &          +yflux(i,j+1,m)-yflux(i,j,m))*areninv(i,j)*dtw2
          enddo
        enddo
      enddo
C -- -------------------------------------------------------------
C                  theta step (MDATA) for wave energy
C ---------------------------------------------------------------
C  Cyclic B.C.s on m
      do j=1,jm
        do i=1,im
        enf(i,j,1)=enf(i,j,mm-1)
        enf(i,j,mm)=enf(i,j,2)
        enddo
      enddo
C
      do m=1,mm
        do j=1,jm
          do i=1,im 
            uak(i,j,m)=csm(m)*ud(i,j)+snm(m)*vd(i,j)
          enddo
        enddo
      enddo
      do j=2,jm-1
        do i=2,im-1 
          dx2=dx(i+1,j)+dx(i-1,j)
          dy2=dy(i,j+1)+dy(i,j-1)
          do m=1,mm
            cthth(i,j,m)=                    ! speed in theta coordinate
     &            Fct(i,j,m)*gee/(2.*cth(i,j,m)*cosh(kthD(i,j,m))*2)
     &         *( sn(m)*(d(i+1,j)-d(i-1,j))/dx2
     &           -cs(m)*(d(i,j+1)-d(i,j-1))/dy2 )
     &      +sn(m)*(uak(i+1,j,m)-uak(i-1,j,m))/dx2*fsm(i+1,j)*fsm(i-1,j)
     &      -cs(m)*(uak(i,j+1,m)-uak(i,j-1,m))/dy2*fsm(i,j+1)*fsm(i,j-1)
          enddo
        enddo
      enddo
      do m=1,mm
        do j=1,jm
          cthth(im,j,m)=cthth(im-1,j,m)
          cthth(1,j,m)=cthth(2,j,m)
         enddo   
         do i=1,im
          cthth(i,jm,m)=cthth(i,jm-1,m)
          cthth(i,1,m)=cthth(i,2,m)
         enddo
       enddo
      call stepmpd(enb,en,enf,cthth,im,jm,mm,dtw2,dthinv,nitera,
     &       thetam,iint,iprint)
C ---------------------------------------------------------------
C                  Source step for wave energy
C ---------------------------------------------------------------
      do i=2,im-1
        do j=2,jm-1
C         Calculate Hm and qdis for later use in depth-induced breaking
          Hm(i,j)=gama*d(i,j)
          fb=8.0*ent(i,j)/gee/(Hm(i,j))**2
          fb=min(fb,1.0)
          do loop=1,3
            qdis(i,j)=exp((qdis(i,j)-1)/fb)
          enddo
C
          ustp(i,j)=sqrt(sqrt(tpx0(i,j)**2+tpy0(i,j)**2))
          do m=2,mm-1
C Input source term
            fspr(i,j,m)=0.
            if(ustp(i,j).gt.0.001) then 
              tht_ops=tht_wnd(i,j)+pi
C Create neg. "tail" to oppose swell. Const. 0.4 from Donelan 1999
              fspr(i,j,m)=fsprm(beta,m,mm,tht_wnd(i,j),thetam(m))
     &               -0.4*fsprm(beta,m,mm,tht_ops     ,thetam(m))
            endif
            swin(i,j,m)=Cin *exp(-swcon *ageinv(i,j))*fsm(i,j)
            swin(i,j,m)=max(swin(i,j,m),10. )*ustp(i,j)**3*fspr(i,j,m)
C Surface dissipation. bdis is used for wind driven portion of en.
C Reduced swfct*bdis is used for swell portion.
            bb=bdis                   
            if(abs(fspr(i,j,m)).lt.0.10) bb=swfct*bdis  
            sdis(i,j,m)=adis*swin(i,j,m)+bb*sigth(i,j,m)*en(i,j,m)
C  sdis is incorporated into the following semi-implicitly.
            enf(i,j,m)=(enf(i,j,m)
     &        +(1.-adis)*swin(i,j,m)*dtw2)/(1.+bb*sigth(i,j,m) *dtw2)
C Bottom dissipation, first due to friction, second to depth-
C induced breaking a la Battjes and Janssen 1978.
           bdis1(i,j,m)=.003*(sigth(i,j,m)
     &             *sqrt(2.*en(i,j,m)/gee)/sinh(kthD(i,j,m)))**3
           bdis2(i,j,m)=(en(i,j,m)/ent(i,j))*gee*sigth(i,j,m)
     &                /(8.*pi)*qdis(i,j)*(Hm(i,j))**2
           enf(i,j,m)=enf(i,j,m)-(bdis1(i,j,m)+bdis2(i,j,m))*dtw2
           enf(i,j,m)=max(.0001,enf(i,j,m))*fsm(i,j)
          enddo
         enddo
       enddo
C ---------------------------------------------------------------
C  Include radiation/current stress terms
      do i=2,im-1
        do j=2,jm-1
          do m=1,mm
            tps(i,j,m)=en(i,j,m)*kth(i,j,m)
          enddo
          kaveD(i,j)=summ(kth,dth,i,j,im,jm,mm)*d(i,j)
        enddo
      enddo
       call cur2wave(im,jm,kb,mm,dx,dy,z,dz,u,v,kaveD,
     &       Sradx,Sradxy,Srady,Srad,Sradpx,Sradpy,Sradp,ud,vd)
      do i=2,im-1
        do j=2,jm-1
          dx2iv=1./(dx(i+1,j)+dx(i-1,j))
          dy2iv=1./(dy(i,j+1)+dy(i,j-1))
          do m=2,mm-1
            sorc=  (csm(m)*csm(m)*Sradx(i,j)
     &            +csm(m)*snm(m)*Sradxy(i,j)
     &            +snm(m)*snm(m)*Srady(i,j) 
     &            +              Srad(i,j)  )*en(i,j,m)
     &         +Sradpx(i,j)*(en(i+1,j,m)-en(i-1,j,m))*dx2iv
     &         +Sradpy(i,j)*(en(i,j+1,m)-en(i,j-1,m))*dy2iv
     &         +Sradp(i,j)*en(i,j,m)
            enf(i,j,m)=enf(i,j,m)+sorc*dtw2
          enddo
        enddo
      enddo 
C ---------------------------------------------------------------
C                     x-y step for frequency
C ---------------------------------------------------------------
      do m=1,mm-1
        do j=2,jm-1
          do i=2,im-1
            cgxm=cg(i,j,m)*csm(m)
            cgym=cg(i,j,m)*snm(m)
            adv=(0.5*cgxm*(sigth(i+1,j,m)-sigth(i-1,j,m))
     &           -0.5*abs(cgxm)
     &            *(sigthb(i+1,j,m)-2.*sigthb(i,j,m)+sigthb(i-1,j,m)) )
     &                  /dx(i,j)
     &         +(0.5*cgym*(sigth(i,j+1,m)-sigth(i,j-1,m))
     &           -0.5*abs(cgym)
     &            *(sigthb(i,j+1,m)-2.*sigthb(i,j,m)+sigthb(i,j-1,m)) )
     &                  /dy(i,j)
            sigthf(i,j,m)=sigthb(i,j,m)-adv*fsm(i,j)*dtw2
          enddo
        enddo
      enddo
C ---------------------------------------------------------------
C                  Source step for frequency     
C ---------------------------------------------------------------
      do i=2,im-1
        do j=2,jm-1
          do m=2,mm-1
            sor=-cg(i,j,m)*kth(i,j,m)*0.5*(
     &          (csm(m)*csm(m))*(ud(i+1,j)-ud(i-1,j))/dx(i,j)
     &         +(csm(m)*snm(m))*(vd(i+1,j)-vd(i-1,j))/dx(i,j)
     &         +(csm(m)*snm(m))*(ud(i,j+1)-ud(i,j-1))/dy(i,j)
     +         +(snm(m)*snm(m))*(vd(i,j+1)-vd(i,j-1))/dy(i,j) )
     &         +cg(i,j,m)*0.5*(sigth(i,j,m)/d(i,j))*(Fcg(i,j,m)-0.5)
     &             *(ud(i,j)*(d(i+1,j)-d(i-1,j))/dx(i,j)
     &              +vd(i,j)*(d(i,j+1)-d(i,j-1))/dy(i,j) )
           sigthf(i,j,m)=sigthf(i,j,m)+sor*dtw2
          enddo
        enddo
      enddo
C--------------------------------------------------------------
C                  Boundary Condtions
C--------------------------------------------------------------
C--- For cyclic B.C.s on the north-south boundaries ----------
c      do m=2,mm
c      do i=1,im
c        enf(i,1,m)=enf(i,jm-1,m)        !cyclic B.C.s
c        enf(i,jm,m)=enf(i,2,m)
c       enddo
c     enddo
C--- For open B.C.s on the eastern boundary -------------------
C     do j=1,jm
c        do m=1,mm
c          cgnd=cgx(im,j,m)/(dx(im,j)+dx(im-1,j))*dtw
c          enf(im,j,m)=en(im,j,m)
c    &       -(cgnd+abs(cgnd))*(en(im,j,m)-en(im-1,j,m))
c    &       -(cgnd-abs(cgnd))*(ene(j,m)-en(im,j,m)    )
c        enddo
c      enddo

C ---- Mask out land values and apply Asselin filter----------
      do m=1,mm
        do j=1,jm
          do i=1,im
            enf(i,j,m)=enf(i,j,m)*fsm(i,j)
            en(i,j,m)=en(i,j,m)
     &        +0.5*smoth*(enb(i,j,m)-2.*en(i,j,m)+enf(i,j,m))
            sigth(i,j,m)=sigth(i,j,m)
     &        +0.5*smoth*(sigthb(i,j,m)-2.*sigth(i,j,m)+sigthf(i,j,m))
          enddo
        enddo
      enddo
C------ Reset time sequence -----------------------------------
      do m=1,mm
        do j=1,jm
          do i=1,im
            enb(i,j,m)=en(i,j,m) 
            en(i,j,m)=enf(i,j,m)
            sigthb(i,j,m)=sigth(i,j,m) 
            sigth(i,j,m)=sigthf(i,j,m)
          enddo
        enddo
      enddo
      do j=1,jm
        do i=1,im
C  Determine average wave propagation angle
          do m=1,mm
            tps(i,j,m)=en(i,j,m)*csm(m)  
          enddo
          csmav=summ(tps,dth,i,j,im,jm,mm)
          do m=1,mm
            tps(i,j,m)=en(i,j,m)*snm(m)  
          enddo
          snmav=summ(tps,dth,i,j,im,jm,mm)
          if(csmav.eq.0.) csmav=csmav+1.e-7
          thtav(i,j)=atan(snmav/csmav)
          if (snmav.gt.0..and.csmav.lt.0.) thtav(i,j)=thtav(i,j)+pi
          if (snmav.lt.0..and.csmav.lt.0.) thtav(i,j)=thtav(i,j)-pi
C Calculate total dissipation for use in tke equation
          sdistot(i,j)=summ(sdis,dth,i,j,im,jm,mm)*fsm(i,j)
          do m=1,mm
            tps(i,j,m)=bdis1(i,j,m)+bdis2(i,j,m)
          enddo
          bdistot(i,j)=summ(tps,dth,i,j,im,jm,mm)           
        enddo
      enddo
C-----------------------------------------------------------
C           Begin Diagnostic and Print Section
C-----------------------------------------------------------
      if(mod(iint,iprint).eq.0) then
      write(6,'(''time,iend,iprint ='',f12.5,2i6)') time,iend,iprint
      do i=1,im
        do j=1,jm
          xfluxtot(i,j)=summ(xflux,dth,i,j,im,jm,mm)   
          yfluxtot(i,j)=summ(yflux,dth,i,j,im,jm,mm)   
          thtavdeg(i,j)=180.*thtav(i,j)/pi 
          swintot(i,j)=summ(swin,dth,i,j,im,jm,mm)
        enddo
      enddo
      iskp=1
      jskp=1
      call prxy('ageinv  ',time,ageinv  ,im,1,jm,1  ,0.)
      call prxy('sigp    ',time,sigp    ,im,1,jm,1  ,0.)
      call prxz('sigth',time,sigth,im,iskp,jm,mm,jo,njo,0.,d,thetam)
c     call prxy('cp      ',time,cp      ,im,1,jm,1  ,0.)
c     call prxy('Fcg     ',time,Fcg     ,im,1,jm,1  ,0.)
c     call prxy('cg      ',time,cg      ,im,1,jm,1  ,0.)
c     call prxy('vd      ',time,vd      ,im,1,jm,1  ,0.)
c     call prxy('Sradxy  ',time,Sradxy  ,im,1,jm,1  ,0.)
c     call prxy('xfluxtot',time,xfluxtot,im,1,jm,1  ,0.)
c     call prxy('yfluxtot',time,yfluxtot,im,1,jm,1  ,0.)
c     call prxy('swintot ',time,swintot ,im,1,jm,1  ,0.)
c     call prxy('bdistot ',time,bdistot ,im,1,jm,1  ,0.)
      call prxy('ent',time,ent,im,1,jm,1  ,0.)
c     call prxy('Hs ',time,Hs ,im,1,jm,1  ,0.)
c     call prxy('thtav   ',time,thtav   ,im,1,jm,1  ,0.   )
      call prxz('fspr',time,fspr,im,iskp,jm,mm,jo,njo,0.,d,thetam)
c     call prxz('cth ',time,cth ,im,iskp,jm,mm,jo,njo,0.,d,theta )
c     call prxz('swin',time,swin,im,iskp,jm,mm,jo,njo,0.,d,thetam)
c     call prxz('sdis',time,sdis,im,iskp,jm,mm,jo,njo,0.,d,thetam)
c     call prxz('cgx ',time,cgx ,im,iskp,jm,mm,jo,njo,0.,d,thetam)
c     call pryz('cgy',time,cgy,im,jm,jskp,mm,io,nio,0.,d,thetam)
      call prxz('enf ',time,enf ,im,iskp,jm,mm,jo,njo,-1.,d,thetam)
c     call pryz('xflux',time,xflux ,im,jm,1,mm,12,1,0.,dt,thetam)
c     call pryz('yflux',time,yflux ,im,jm,1,mm,12,1,0.,dt,thetam)
c     call pryz('cgy  ',time,cgy   ,im,jm,1,mm,12,1,0.,dt,thetam)
c     call pryz('enf  ',time,enf   ,im,jm,1,mm,12,1,0.,dt,thetam)
C Print out a j-slice 
      j=(jm+1)/2
      write(6,'(''Print properties for j ='',i5)') j
      thdegwd=tht_wnd(im/2,j)*180./pi
      write(6,'(''time ='',f9.5,''      u10x,u10y ='',2f8.2,
     & ''    thdegwd ='',f10.2 )')
     &         time, u10x(im/2,j),u10y(im/2,j),thdegwd
      write(6,'(''    i     x       d      cp    cg(21)  ud    '',
     &  '' vd     Tp  '',
     &''  ent theta/deg  Hs   ageinv      Cd        xnd        ennd'')')
      do i=1,im,1
        tp=2.*pi/sigp(i,j)
        xnd=x(i)*gee/u10(i,j)**2
        ennd=ent(i,j)*gee/u10(i,j)**4
        write(6,'(i5,f9.0,f7.0,9f7.2,3e12.3)')
     &     i,x(i),d(i,j),cp(i,j),cg(i,j,21),ud(i,j),vd(i,j),
     &     tp,ent(i,j),
     &     thtavdeg(i,j),Hs(i,j), ageinv(i,j),Cd(i,j),xnd,ennd
      enddo
      endif
      return  
      end
C---------------------------------------------------------------
       subroutine botdrag (im,jm,fsm,
     &     kp,cp,ent,wubot,wvbot,d,zzkbm1,z0b,cbcmin,cbcmax,cbc)
C--------------------------------------------------------------
C
C  Input:
C    fsm(I,j) = 1 for water cell, = 0 fro land cell
C    cp(i,j) = peak frequency phase speed (m s-1)
C    kp(i,j) = peak frequency wave number
C    ent(i,j) = wave energy averaged over all wave angles
C    wubot(i,j),wvbot(i,j) = bottom stress vector as determned
C      by bottom current (at zz(kb-1)) and cbc found here
C    zzkbn1 = zz(kb-1) (non-d.)
C    d(i,j) = water column depth (m)
C    zob = roughness parameter (m)
C    cbcmin,cbcmax = limits ob cbs (non-d.)    
C
C  Output:
C    cbc(i,j) = bottom drag coefficient (non-d.)
C--------------------------------------------------------------
      real kappa,kp
      dimension kp(im,jm),cp(im,jm),ent(im,jm),fsm(im,jm),
     &     wubot(im,jm),wvbot(im,jm),d(im,jm),cbc(im,jm)
      data kappa/0.4/,grav/9.807/
      do j=1,jm
        do i=1,im
          if (fsm(i,j).eq.1.) then
            uboscil=cp(i,j)*kp(i,j)*sqrt(2.*ent(i,j)/grav)
     &               /sinh(kp(i,j)*d(i,j))
            utau2=sqrt(wubot(i,j)**2+wvbot(i,j)**2)
            z0a=z0b*(1.+0.05*uboscil**2/utau2)
            cbc(i,j)=(kappa/log(1.+(1.0+zzkbm1)*d(i,j)/z0a))**2
            cbc(i,j)=max(cbcmin,cbc(i,j))
C
C     If the following is invoked, then it is probable that the wrong
C     choice of z0b or vertical spacing has been made:
C
            cbc(i,j)=min(cbcmax,cbc(i,j))
          endif
        enddo
      enddo
      return
      end

C--------------------------------------------------------------------
      subroutine depth
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Establishes the vertical sigma grid with log        *
C *                distributions at the top and bottom and a linear    *
C *                distribution in between. The number of layers of    *
C *                reduced thickness are kl1-2 at the surface and      *
C *                kb-kl2-1 at the bottom. kl1 and kl2 are defined in  *
C *                the main program. For no log portions, set kl1=2    *
C *                and kl2=kb-1.                                       *
C *                                                                    *
C **********************************************************************
C
C
      include 'wav.c'
C     implicit none
C
      real delz
      integer kdz(12)
      integer k
      data kdz/1,1,2,4,8,16,32,64,128,256,512,1024/
C
      kl1=5
      kl2=kb-5
      z(1)=0.e0
      do k=2,kl1
        z(k)=z(k-1)+kdz(k-1)
      end do
      delz=z(kl1)-z(kl1-1)
      do k=kl1+1,kl2
        z(k)=z(k-1)+delz
      end do
      do k=kl2+1,kb
        dz(k)=float(kdz(kb-k+1))*delz/float(kdz(kb-kl2))
        z(k)=z(k-1)+dz(k)
      end do
      do k=1,kb
        z(k)=-z(k)/z(kb)
      end do
      do k=1,kb-1
        zz(k)=0.5e0*(z(k)+z(k+1))
      end do
      zz(kb)=2.e0*zz(kb-1)-zz(kb-2)
      do k=1,kb-1
        dz(k)=z(k)-z(k+1)
        dzz(k)=zz(k)-zz(k+1)
      end do
      dz(kb)=0.e0
      dzz(kb)=0.e0
C
      write(6,1)
    1 format(/2x,'k',7x,'z',9x,'zz',9x,'dz',9x,'dzz',/)
      do k=1,kb
        write(6,2) k,z(k),zz(k),dz(k),dzz(k)
    2   format((' ',i5,4f10.3))
      end do
      write(6,3)
    3 format(//)
C
      return
      end
C--------------------------------------------------------------

      subroutine cur2wave(im,jm,kb,mm,dx,dy,z,dz,u,v,kpD,
     &      Sradx,Sradxy,Srady,Srad,Sradpx,Sradpy,Sradp,ud,vd)
C--------------------------------------------------------------------
C                  Current to wave interaction terms
C     Calculates vertical averages using spectrally averaged functions
C     provided by subroutine specavs.
C--------------------------------------------------------------------
C   Input       
      real kpD(im,jm)
      dimension dx(im,jm),dy(im,jm),z(kb),dz(kb),u(im,jm,kb),v(im,jm,kb)
C   Output
      dimension  Sradx(im,jm),Sradxy(im,jm),Srady(im,jm),Srad(im,jm),
     &    Sradpx(im,jm),Sradpy(im,jm),Sradp(im,jm),ud(im,jm),vd(im,jm)
C   Internal
      dimension tps(im,jm,mm),F1(kb),F2(kb),F3(kb),F4(kb),F5(kb),
     &    tp1(kb),tp2(kb),tp3(kb),tp4(kb),tp5(kb),
     &    gradxu(kb),gradxv(kb),gradyu(kb),gradyv(kb)
      do j=1,jm
        do i=1,im  
          call specavs (kpD(i,j),F1,F2,F3,F4,F5,z,kb,i,j)
        enddo
      enddo
      do j=3,jm-2
        do i=3,im-2  
          dx2iv=1./(dx(i+1,j)+dx(i-1,j))   
          dy2iv=1./(dy(i,j+1)+dy(i,j-1))   
          do k=1,kb-1
            tp1(k)=u(i,j,k)*(F5(k)-F5(k+1))
            tp2(k)=v(i,j,k)*(F5(k)-F5(k+1))
            gradxu(k)=(u(i+1,j,k)-u(i-1,j,k))*dx2iv
            gradxv(k)=0.5*(v(i+1,j,k)-v(i-1,j,k)
     &             +v(i+1,j+1,k)-v(i-1,j+1,k))*dx2iv
            gradyu(k)=0.5*(u(i,j+1,k)-u(i,j-1,k)
     &             +u(i+1,j+1,k)-u(i+1,j-1,k))*dy2iv
            gradyv(k)=(v(i,j+1,k)-v(i,j-1,k))*dy2iv
          enddo
          ud(i,j)=sumk(tp1,1.,kb)    ! Dopler velocity
          vd(i,j)=sumk(tp2,1.,kb)    !   ''     ''
C
          do k=1,kb-1
            tp1(k)=0.5*(F1(k)+F1(k+1))*gradxu(k)*dz(k)
            tp2(k)=0.5*(F1(k)+F1(k+1))*(gradxv(k)+gradyu(k))*dz(k)
            tp3(k)=0.5*(F1(k)+F1(k+1))*gradyv(k)*dz(k)
          enddo
          Sradx(i,j)=sumk(tp1,1.,kb)
          Sradxy(i,j)=sumk(tp2,1.,kb)
          Srady(i,j)=sumk(tp3,1.,kb)
          if(kpD(i,j).lt.3.) then 
            do k=1,kb-1
              tp1(k)=0.5*(F2(k)+F2(k+1))*gradxu(k)*dz(k)
              tp2(k)=0.5*(F2(k)+F2(k+1))*gradyv(k)*dz(k)
            enddo
            Srad(i,j)=sumk(tp1,1.,kb)+sumk(tp2,1.,kb)
          else
            Srad(i,j)=0.
          endif
          if(kpD(i,j).lt.3.) then
            do k=1,kb-2
              delu=(u(i,j,k)+u(i+1,j,k)-u(i,j,k+1)-u(i+1,j,k+1))
              tp1(k)=0.5*F3(k+1)*delu
              tp2(k)=0.5*F4(k+1)*delu
              delv=(v(i,j,k)+v(i,j+1,k)-v(i,j,k+1)-v(i,j+1,k+1))
              tp3(k)=0.5*F3(k+1)*delv
              tp4(k)=0.5*F4(k+1)*delv
            enddo
            Sradpx(i,j)=sumk(tp1,1.,kb)
            Sradpy(i,j)=sumk(tp3,1.,kb)
            Sradp(i,j)
     &              =(kpD(i+1,j)-kpD(i-1,j))*sumk(tp2,1.,kb)*dx2iv
     &              +(kpD(i,j+1)-kpD(i,j-1))*sumk(tp4,1.,kb)*dy2iv
          else
            Sradpx(i,j)=0.
            Sradpy(i,j)=0.
            Sradp(i,j)=0.
          endif
C For top vertical spacing too large to resolve F1
          if((kpD(i,j)*z(2)).lt.-0.5) then     
            Sradx(i,j)=0.520*gradxu(1)
            Sradxy(i,j)=0.520*(gradxv(1)+gradyu(1))
            Srady(i,j)=0.520*gradyv(1)
          endif
        enddo
      enddo
c     call prxy('Srad    ',time,Srad    ,im,1,jm,1  ,0.)
      return
      end
C----------------------------------------------------------------

      subroutine specavs (kpD,F1,F2,F3,F4,F5,z,kb,i,j)
C       This is a subroutine to deliver spectrally averaged F1 - F5.
C       averaged. Fcg and Fct are functions of kpD.
C       F1-F5 are spectrally averaged components of wave radiation
C       stresses. F1-F5 are functions of kpD and z.
C       N.B. This subroutine uses the file, specavs.
      parameter (me=17,ne=21)
      real kpD,kpDd(me)
      dimension F1d(me,ne),F2d(me,ne),F3d(me,ne),F4d(me,ne),F5d(me,ne) 
      dimension F1i(ne),F2i(ne),F3i(ne),F4i(ne),F5i(ne)
      dimension F1(kb),F2(kb),F3(kb),F4(kb),F5(kb)
      dimension zeta(ne),zetf(ne),z(kb) 
      integer init
      save kpDd,F1d,F2d,F3d,F4d,F5d,zeta,zetf                   
      data init/0/
      if(init.eq.0) then
        write(6,'(''Read spectrally averaged functions from specavs'')')
        read(10,'(2x)')
        write(6,'(/,''   F1'')')
        do n=1,ne 
          read(10,'(f6.2,15(f6.3),2f8.3)')
     &              zeta(n),(F1d(m,n),m=2,me),zetf(n)
          F1d(1,n)=0.
C         write(6,'(f6.2,16(f6.3),2f8.3)')
C    &              zeta(n),(F1d(m,n),m=1,me),zetf(n)
        enddo
        read(10,'(2x)')
        write(6,'(/,''   F2'')')
        do n=1,ne 
          read(10,'(f6.2,15(f6.3),2f8.3)')
     &              zeta(n),(F2d(m,n),m=2,me),zetf(n)
          F2d(1,n)=-zeta(n)
c         write(6,'(f6.2,16(f6.3),2f8.3)')
c    &              zeta(n),(F2d(m,n),m=1,me),zetf(n)
        enddo
        read(10,'(2x)')
        write(6,'(/,''   F3'')')
        do n=1,ne 
          read(10,'(f6.2,15(f6.3),2f8.3)')
     &              zeta(n),(F3d(m,n),m=2,me),zetf(n)
          F3d(1,n)=-zeta(n)*(1.+zeta(n))*0.5
c         write(6,'(f6.2,16(f6.3),2f8.3)')
c    &              zeta(n),(F3d(m,n),m=1,me),zetf(n)
        enddo
        read(10,'(2x)')
        write(6,'(/,''   F4'')')
        do n=1,ne 
          read(10,'(f6.2,15(f6.3),2f8.3)')
     &              zeta(n),(F4d(m,n),m=2,me),zetf(n)
          F4d(1,n)=-zeta(n)*(1.+zeta(n))/.01
c         write(6,'(f6.2,16(f6.3),2f8.3)')
c    &              zeta(n),(F4d(m,n),m=1,me),zetf(n)
        enddo
        read(10,'(2x)')
        write(6,'(/,''   F5'')')
        do n=1,ne 
          read(10,'(f6.2,15(f6.3),2f8.3)')
     &              zeta(n),(F5d(m,n),m=2,me),zetf(n)
          F5d(1,n)=1.+zeta(n)
C         write(6,'(f6.2,16(f6.3),2f8.3)')
C    &              zeta(n),(F5d(m,n),m=1,me),zetf(n)

       enddo 
       init=1
      endif  
C----------------------------------------------------------------
C  Solve for F1, F2, F3, F4, F5
C----------------------------------------------------------------
      m=kpD/0.2+1
c      write(6,'(''i,j,m = '',3i5,''kpD'',e12.3)') i,j,m,kpD
      if (m.lt.16) then
C Interpolate w.r.t. kpD =< 3. (0 < m< 17)
        do n=1,ne
          F1i(n)=F1d(m,n)+(F1d(m+1,n)-F1d(m,n))*(kpD-kpDd(m))/0.2
          F2i(n)=F2d(m,n)+(F2d(m+1,n)-F2d(m,n))*(kpD-kpDd(m))/0.2
          F3i(n)=F3d(m,n)+(F3d(m+1,n)-F3d(m,n))*(kpD-kpDd(m))/0.2
          F4i(n)=F4d(m,n)+(F4d(m+1,n)-F4d(m,n))*(kpD-kpDd(m))/0.2
          F5i(n)=F5d(m,n)+(F5d(m+1,n)-F5d(m,n))*(kpD-kpDd(m))/0.2
        enddo
C Interpolate w.r.t. z
        do k=1,kb-1
          n=-z(k)/.05+1
          F1(k)=F1i(n)+(F1i(n+1)-F1i(n))*(-z(k)+zeta(n))/.05
          F2(k)=F2i(n)+(F2i(n+1)-F2i(n))*(-z(k)+zeta(n))/.05
          F3(k)=F3i(n)+(F3i(n+1)-F3i(n))*(-z(k)+zeta(n))/.05
          F4(k)=F4i(n)+(F4i(n+1)-F4i(n))*(-z(k)+zeta(n))/.05
          F5(k)=F5i(n)+(F5i(n+1)-F5i(n))*(-z(k)+zeta(n))/.05
        enddo
        F1(kb)=F1i(ne)
        F2(kb)=F2i(ne)
        F3(kb)=F3i(ne)
        F4(kb)=F4i(ne)
        F5(kb)=F5i(ne)
      else
C Use similarity for large kpD > 3. (m = 17)
C e.g., F1(kpD, z) = (kpD/100)*F1(100, kpD*z/100)
        do k=1,kb
          F1(k)=0.
          F2(k)=0.
          F3(k)=0.
          F4(k)=0.
          F5(k)=0.
          n=1-z(k)*kpD/.1
          if(n.lt.ne) then
            F1(k)=F1d(17,n)+(F1d(17,n+1)-F1d(17,n))/.001
     &              *(-z(k)*kpD/100.+zetf(n))
            F1(k)=F1(k)*kpD*0.01
            F5(k)=F5d(17,n)+(F5d(17,n+1)-F5d(17,n))/.001
     &              *(-z(k)*kpD/100.+zetf(n))
          endif
        enddo
      endif
      return  
      end
C----------------------------------------------------------------
       subroutine speeds(kpD,Fcg,Fct)
C       This is a subroutine to deliver spectrally averaged Fcg, Fct
C       Fcg = ratio of cg to cp where cp is spectral peak phase speed.
C       Fct = is factor in cthth (theta speed) to make it a spectrally
C       averaged. Fcg and Fct are functions of kpD.
       parameter (me=17)
       real kpD,kpDd(me)
       dimension Fcgd(me),Fctd(me),c_cp(me),dum(me)
       data kpDd/
     &    0.0,   0.2,   0.4,   0.6,   0.8,   1.0,   1.2,   1.4,   1.6,
     &    1.8,   2.0,   2.2,   2.4,   2.6,   2.8,   3.0, 100.0/
       data c_cp/
     &  1.000, 1.005, 1.001, 0.968, 0.951, 0.936, 0.924, 0.915, 0.910,
     &  0.906, 0.905, 0.904, 0.904, 0.904, 0.905, 0.905, 0.910/
       data Fcgd/
     &  1.000, 0.984, 0.928, 0.834, 0.756, 0.689, 0.635, 0.592, 0.560,
     &  0.535, 0.516, 0.502, 0.491, 0.483, 0.477, 0.472, 0.472/
       data dum/
     &  1.000, 0.987, 0.950, 0.897, 0.837, 0.776, 0.720, 0.671, 0.631,
     &  0.598, 0.573, 0.554, 0.539, 0.529, 0.521, 0.515, 0.500/
       data Fctd/
     &  1.000, 0.984, 0.954, 0.931, 0.892, 0.858, 0.830, 0.809, 0.795,
     &  0.787, 0.784, 0.786, 0.791, 0.801, 0.814, 0.829, 0.829/

        k=kpD/0.2+1
      if (k.lt.16) then
        Fcg=Fcgd(k)+(Fcgd(k+1)-Fcgd(k))*(kpD-kpDd(k))/0.2
        Fct=Fctd(k)+(Fctd(k+1)-Fctd(k))*(kpD-kpDd(k))/0.2
      else
        k=16
        Fcg=Fcgd(k)+(Fcgd(k+1)-Fcgd(k))*(kpD-kpDd(k))/97.
        Fct=Fctd(k)+(Fctd(k+1)-Fctd(k))*(kpD-kpDd(k))/97.
      endif
      return
      end
C------------------------------------------------------------------

      subroutine stress (Hs,ageinv,usdif,vsdif,cd,tpx0,tpy0,
     &                          ttx0,tty0,ustw,im,jm)
      dimension Hs(im,jm),ageinv(im,jm),usdif(im,jm),vsdif(im,jm)  
      dimension cd(im,jm),tpx0(im,jm),tpy0(im,jm),ustw(im,jm) 
      dimension ttx0(im,jm),tty0(im,jm)
      real kappa,nu
      data nu/1.8e-6/,r/0.0011630/,rhalf/.0034100/,kappa/0.41/
      data z0w/1.e-6/,z0t/1.e-6/,gee/9.807/,pi/3.1415926/
      do j=2,jm-1
        do i=2,im-1
          u2=usdif(i,j)**2+vsdif(i,j)**2
          u1=sqrt(u2)
          z0w=1.38e-4*Hs(i,j)*(ageinv(i,j))**2.66+1.e-5   !Donelan
          z0t=0.18*nu/(ustw(i,j) +.000001)
          cdw=(kappa/log(10./z0w))**2
          cdt=(kappa/log(10./z0t))**2
          if(cdw.gt.cdt) then
            cdt=0.
          else
            cdw=0.
          endif
          cd(i,j)=cdw+cdt
          tpx0(i,j)=r*cdw*u1*usdif(i,j)
          tpy0(i,j)=r*cdw*u1*vsdif(i,j)
          ttx0(i,j)=r*cdt*u1*usdif(i,j)    
          tty0(i,j)=r*cdt*u1*vsdif(i,j)   
          ustw(i,j)=sqrt(r*(cdw+cdt)*u2)  
c        if(i.eq.2.and. j.eq.3) then
c        write(6,'(''n,usdif,vsdif ='',i5,2f10.3,3e12.3)')
c    &              n,usdif(i,j),vsdif(i,j),z0,Hs(i,j),sigp(i,j)
c        endif
        enddo
      enddo
      return
      end
C--------------------------------------------------------------------

      function fsprm(beta,m,mm,tht_wnd,thetam)
      save tht,fdat,dtht
      dimension tht(35),fdat(35) 
      data pi/3.141592654/
      data ifirst/0/
      if(ifirst.eq.0) then
        write(6,'(''spreading function'',/,
     &            ''    m     tht       fdat'')')
       dtht=pi/float(32)
        do k=1,35 
          fdat(k)=0.
          tht(k)=dtht*(k-1)
          if(k.lt.17) fdat(k)=0.5*beta/cosh(beta*tht(k))**2
          write(6,'(i5,2f10.4)') k,tht(k),fdat(k)
        enddo
      endif
      ifirst=1
C Create spreading function around the wind direction
      fsprm=0.
      thtrel=abs(thetam-tht_wnd)
      if(thtrel.gt.pi) thtrel=abs(thtrel-2.*pi)  
      kth=thtrel/dtht+1  
      fsprm=fdat(kth)
     &    +(fdat(kth+1)-fdat(kth))*(thtrel-tht(kth))/dtht   
      return
      end
C----------------------------------------------------------------
      function mxloc(a,b,i,j,im,jm,mm)
      dimension a(im,jm,mm),b(im,jm,mm)
      amx=0.
      mmx=1
      do m=2,mm
        if(b(i,j,m).gt.0.10) then
          if(a(i,j,m).gt.amx) then
            amx=a(i,j,m)
            mmx=m
          endif
        endif
c       if(i.eq.16.and.j.eq.31)then
c      write(6,'(''mxloc'',4i5,3e12.3)') i,j,m,mmx,amx,a(i,j,m),b(i,j,m)
c       endif
      enddo
      mxloc=mmx
      return
      end
C----------------------------------------------------------------


      subroutine stepmpd(fb,f,ff,c,im,jm,mm,dtw2,dthinv,nitera,
     &         thetam, iint,iprint)
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 *                                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
      parameter (value_min=1.e-9,epsilon=1.0e-14)
      real sw
      integer nitera
      real fbmem(im,jm,mm),ffmem(im,jm,mm),c(im,jm,mm)
      real flux(im,jm,mm),speed(im,jm,mm),thetam(mm)
      integer i,j,k,itera
      real fb(im,jm,mm),f(im,jm,mm),ff(im,jm,mm)
      real mol,abs_1,abs_2
      real value_min,epsilon
      real udx,u2dt,vdy,v2dt,wdz,w2dt 
      sw=1.0  
C
      
      do m=1,mm
        do j=1,jm
          do i=1,im
            speed(i,j,m)=c(i,j,m)
            fbmem(i,j,m)=fb(i,j,m)
            ffmem(i,j,m)=ff(i,j,m)
          end do
        end do
      end do
C
      do itera=1,nitera
C
C     Upwind advection scheme:
C
c       do j=1,jm-1
c         do i=2,im-1
c           flux(i,j,1)=0.e0
c           if(itera.eq.1) flux(i,j,1)=c(i,j,1)*f(i,j,1)
c           flux(i,j,mm)=0.e0
c         end do
c       end do
C
        do m=2,mm-1
          do j=1,jm
            do i=1,im
              flux(i,j,m)=0.5e0
     $                      *((speed(i,j,m)+abs(speed(i,j,m)))
     $                       *fbmem(i,j,m-1)+
     $                        (speed(i,j,m)-abs(speed(i,j,m)))
     $                       *fbmem(i,j,m))
            end do
          end do
        end do
        do j=1,jm
          do i=1,im
            flux(i,j,1)=flux(i,j,mm-1)
            flux(i,j,mm)=flux(i,j,2)
          end do
        end do
C
C     Add net advective fluxes and step forward in time:
C
        do m=1,mm-1
          do j=1,jm
            do i=1,im
              ff(i,j,m)=ffmem(i,j,m)
     $             +(flux(i,j,m)-flux(i,j,m+1))*dthinv*dtw2
            end do
          end do
        end do
c     if(mod(iint,iprint).eq.0) then
c     write(6,'(''iint ='',i5)') iint
c     call pryz('ff    ',time,ff    ,im,jm,1,mm,5,1,0.,dt,thetam)
c     endif
C
C *                Calculate the antidiffusive speed used to        *
C *                reduce the numerical diffusion associated with the  *
C *                upstream differencing scheme.                       *
C *                                                                    *
        do m=2,mm-1
          do j=1,jm
            do i=1,im
            if(ff(i,j,m).lt.value_min.or.
     $         ff(i,j,m-1).lt.value_min) then
              speed(i,j,m)=0.e0
            else
              wdz=abs(speed(i,j,m))
              w2dt=speed(i,j,m)*speed(i,j,m)*dtw2*dthinv              
              mol=(ff(i,j,m)-ff(i,j,m-1))
     $             /(ff(i,j,m)+ff(i,j,m-1)+epsilon)
              speed(i,j,m)=(wdz-w2dt)*mol*sw
              abs_1=abs(wdz)
              abs_2=abs(w2dt)
              if(abs_1.lt.abs_2)speed(i,j,m)=0.e0
            end if
          end do
        end do
      end do
C
        do m=1,mm
          do j=1,jm
            do i=1,im
              ffmem(i,j,m)=ff(i,j,m)
              fbmem(i,j,m)=ff(i,j,m)
            end do
          end do
        end do
C
C     End of Smolarkiewicz scheme
      end do
      return
      end
C----------------------------------------------------------------

      function summ(a,dpi,i,j,im,jm,mm)
      dimension a(im,jm,mm)
      summ=0.
      do m=2,mm-1 
         summ=summ+a(i,j,m)*dpi  
      enddo
      return
      end
C----------------------------------------------------------------
      function summ1(a,dpi,mm)
      dimension a(mm)
      summ1=0.
      do m=2,mm-1
         summ1=summ1+a(m)*dpi
      enddo
      return
      end
C----------------------------------------------------------------
      function sumk(a,dpi,kb)
      dimension a(kb)
      sumk=0.
      do k=2,kb-1 
         sumk=sumk+a(k)*dpi  
      enddo
      return
      end
C-----------------------------------------------------------------
      subroutine wave2cur (dth,Sxx,Sxy,Syy,Spx,Spy)
C  This routine supplies Sxx,Sxy.Syy,Spx,Spy for insertion
C  into the momentum equation when waves are coupled with 
C  pom2k.f

      include 'wav.c'  
      dimension Sxx(im,jm,kb),Sxy(im,jm,kb),Syy(im,jm,kb)
      dimension Spx(im,jm,kb),Spy(im,jm,kb)
      dimension F1(kb),F2(kb),F3(kb),F4(kb),F5(kb),
     &    tpxx(mm),tpxy(mm),tpyy(mm)
      do j=1,jm
        do i=1,im  
          do m=1,mm
            tpxx(m)=en(i,j,m)*csm(m)*csm(m)         
            tpxy(m)=en(i,j,m)*csm(m)*snm(m)
            tpyy(m)=en(i,j,m)*snm(m)*snm(m)
          enddo
            stpxx=summ1(tpxx,dth,mm)
            stpxy=summ1(tpxy,dth,mm)
            stpyy=summ1(tpyy,dth,mm)
            call specavs (kaveD(i,j),F1,F2,F3,F4,F5,z,kb,i,j)
          do k=1,kb-1
            F1m=0.5*(F1(k)+F1(k+1))
            F2m=0.5*(F2(k)+F2(k+1))
            F3m=0.5*(F3(k)+F3(k+1))
            F4m=0.5*(F4(k)+F4(k+1))
            Sxx(i,j,k)=stpxx*F1m+ent(i,j)*F2m
            Sxy(i,j,k)=stpxy*F1m
            Syy(i,j,k)=stpyy*F1m+ent(i,j)*F2m
            dxinv=1./(dx(i+1,j)+dx(i-1,j))
            Spx(i,j,k)=(ent(i+1,j)-ent(i-1,j))*dxinv*F3m
     &         +ent(i,j)*(kaveD(i+1,j)-kaveD(i-1,j))*dxinv*F4m      
            Spx(i,j,k)=Spx(i,j,k)*fsm(i+1,j)*fsm(i-1,j)
            dyinv=1./(dy(i,j+1)+dy(i,j-1))
            Spy(i,j,k)=(ent(i,j+1)-ent(i,j-1))*dyinv*F3m
     &         +ent(i,j)*(kaveD(i,j+1)-kaveD(i,j-1))*dyinv*F4m
            Spy(i,j,k)=Spy(i,j,k)*fsm(i,j+1)*fsm(i,j-1)
          enddo
        enddo
      enddo
      return
      end
C --------------------------------------------------------------
      subroutine wvfcns(sig,d,wvnu,c)   
      parameter(geeinv=1./9.807)
C
C This lookup table inputs sig and d and outputs wvnu, c and cg  
C using the dispersion relation
C
C sig = frequency (s-1)
C d = water column depth (m) 
C wvnu = wave number (m-1)
C c = phase speed (m s-1)
C cg = group speed (m s-1)  
C f1 = sig**2*d/gee
C f2 = wvnu*d
C f3 = cg/c
C
C ------------------------------------------------------------------
      dimension f1(81),f2(81),f3(81)
      data f1/ 
     &  0.0000,  0.0500,  0.1000,  0.1500,  0.2000,  0.2500,  0.3000,
     &  0.3500,  0.4000,  0.4500,  0.5000,  0.5500,  0.6000,  0.6500,
     &  0.7000,  0.7500,  0.8000,  0.8500,  0.9000,  0.9500,  1.0000,
     &  1.0500,  1.1000,  1.1500,  1.2000,  1.2500,  1.3000,  1.3500,
     &  1.4000,  1.4500,  1.5000,  1.5500,  1.6000,  1.6500,  1.7000,
     &  1.7500,  1.8000,  1.8500,  1.9000,  1.9500,  2.0000,  2.0500,
     &  2.1000,  2.1500,  2.2000,  2.2500,  2.3000,  2.3500,  2.4000,
     &  2.4500,  2.5000,  2.5500,  2.6000,  2.6500,  2.7000,  2.7500,
     &  2.8000,  2.8500,  2.9000,  2.9500,  3.0000,  3.0500,  3.1000,
     &  3.1500,  3.2000,  3.2500,  3.3000,  3.3500,  3.4000,  3.4500,
     &  3.5000,  3.5500,  3.6000,  3.6500,  3.7000,  3.7500,  3.8000,
     &  3.8500,  3.9000,  3.9500,  4.0000/
      data f2/
     &  0.0000,  0.2255,  0.3216,  0.3973,  0.4627,  0.5218,  0.5767,
     &  0.6284,  0.6778,  0.7255,  0.7717,  0.8168,  0.8611,  0.9046,
     &  0.9476,  0.9902,  1.0324,  1.0744,  1.1163,  1.1580,  1.1997,
     &  1.2414,  1.2831,  1.3249,  1.3668,  1.4088,  1.4511,  1.4934,
     &  1.5360,  1.5788,  1.6218,  1.6651,  1.7085,  1.7523,  1.7962,
     &  1.8405,  1.8849,  1.9297,  1.9746,  2.0199,  2.0653,  2.1110,
     &  2.1569,  2.2031,  2.2495,  2.2960,  2.3428,  2.3898,  2.4369,
     &  2.4843,  2.5318,  2.5795,  2.6273,  2.6752,  2.7233,  2.7716,
     &  2.8199,  2.8684,  2.9170,  2.9657,  3.0144,  3.0633,  3.1123,
     &  3.1613,  3.2104,  3.2596,  3.3088,  3.3581,  3.4074,  3.4568,
     &  3.5063,  3.5558,  3.6053,  3.6548,  3.7044,  3.7541,  3.8037,
     &  3.8500,  3.9000,  3.9500,  4.0000/
      data f3/
     &  1.0000,  0.9834,  0.9671,  0.9510,  0.9352,  0.9196,  0.9042,
     &  0.8891,  0.8743,  0.8598,  0.8455,  0.8315,  0.8179,  0.8045,
     &  0.7914,  0.7786,  0.7662,  0.7541,  0.7422,  0.7308,  0.7196,
     &  0.7088,  0.6983,  0.6882,  0.6784,  0.6689,  0.6598,  0.6511,
     &  0.6426,  0.6345,  0.6268,  0.6193,  0.6122,  0.6054,  0.5990,
     &  0.5928,  0.5870,  0.5814,  0.5761,  0.5711,  0.5664,  0.5619,
     &  0.5577,  0.5538,  0.5500,  0.5465,  0.5432,  0.5401,  0.5373,
     &  0.5345,  0.5320,  0.5297,  0.5274,  0.5254,  0.5235,  0.5217,
     &  0.5200,  0.5185,  0.5171,  0.5157,  0.5145,  0.5134,  0.5123,
     &  0.5114,  0.5104,  0.5096,  0.5088,  0.5081,  0.5075,  0.5069,
     &  0.5063,  0.5058,  0.5053,  0.5049,  0.5045,  0.5041,  0.5038,
     &  0.5000,  0.5000,  0.5000,  0.5000/
      fqnd=sig*sig*d*geeinv
      k=min(1+fqnd*20.,80.)  
        dknd=f2(k)+(f2(k+1)-f2(k))*(fqnd-f1(k))*20.
        cgnd=f3(k)+(f3(k+1)-f3(k))*(fqnd-f1(k))*20.
        if(fqnd.lt.0.2) dknd=sqrt(fqnd)
      wvnu=dknd/d
      c=sig/wvnu
      cg=cgnd*c
c     if(j.eq.5) then
c     if(i.eq.4.or.i.eq.8) then
c     write(6,'(2i5,3e10.3,/)')i,sig,d,fqnd   
c     write(6,'(i5,3e10.3,/)')k,f1(k),f2(k),f3(k)
c     write(6,'(5e10.3,/)')dknd,cgnd,wvnu,c,cg 
c     endif
c     endif
      return
      end
C----------------------------------------------------------------------
      subroutine prxy(label,time,a,im,iskp,jm,jskp,scala)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Writes a horizontal 2-D field.                      *
C *                                                                    *
C *                label ....... label for output                      *
C *                time ........ time (days)                           *
C *                a(im,jm,kb).. array to be printed                   *
C *                iskp ........ skipping interval for i               *
C *                jskp ........ skipping interval for j               *
C *                scala ....... < 0 for floating point numbers output *
C *                              0 for integer output, divisor for a   *
C *                                based on magnitudes of |a| values   *
C *                              > 0 for integer output, divisor for a *
C *                                given by scala                      *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      integer im,jm
      real a(im,jm)
      real time,scala
      integer iskp,jskp
      character label*(*)
      real amx,scale
      integer i,ib,ie,j,jwr,cols
C
      if(scala.ge.0.e0) then
        cols=24
      else
        cols=12
      endif
C
      if (scala.lt.0.e0) scale = 1.e0
      if (scala.eq.0.e0) then
        amx=1.e-12
        do j=1,jm,jskp
          do i=1,im,iskp
            amx=max(abs(a(i,j)),amx)
          end do
        end do
          scale=10.e0**(int(log10(amx)+100.e0)-103)
        endif
      if(scala.gt.0.e0) scale=scala
C
      write(6,1) label
    1 format(1x,a40/)
      write(6,2) time,scale
    2 format(' Time = ',f9.5,' days    multiply all values by ',1pe8.2)
C
      do ib=1,im,cols*iskp
C
        ie=ib+(cols-1)*iskp
        if(ie.gt.im) ie=im
C
        if(scala.ge.0.e0) then
          write(6,3) (i,i=ib,ie,iskp)
    3     format(/,2x,24i5,/)
        else
          write(6,4) (i,i=ib,ie,iskp)
    4     format(/,12i10,/)
        endif
C
        do j=1,jm,jskp
          jwr=jm+1-j
          if(scala.ge.0.e0) then
            write(6,5) jwr,(nint(a(i,jwr)/scale),i=ib,ie,iskp)
    5       format(1x,i3,24i5)
          else
            write(6,6) jwr,(a(i,jwr),i=ib,ie,iskp)
    6       format(1x,i2,12(e10.2))
          endif
        end do
C
        write(6,7)
    7   format(//)
C
      end do
C
      return
      end
C---------------------------------------------------------------------

      subroutine prxyz(label,time,a,im,iskp,jm,jskp,kb,ko,nko,scala)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Writes horizontal layers of a 3-D field with        *
C *                integers or floating point numbers.                 *
C *                                                                    *
C *                label ....... label for output                      *
C *                time ........ time (days)                           *
C *                a(im,jm,kb).. array to be printed                   *
C *                iskp ........ skipping interval for i               *
C *                jskp ........ skipping interval for j               *
C *                ko .......... 1-D array of k-indices for output     *
C *                nko ......... number of elements in ko              *
C *                scala ....... < 0 for floating point numbers output *
C *                              0 for integer output, divisor for a   *
C *                                based on magnitudes of |a| values   *
C *                              > 0 for integer output, divisor for a *
C *                                given by scala                      *
C *                                                                    *
C *                (NOTE that this combines functions of old prxyz and *
C *                 eprxyz)                                            *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      integer im,jm,kb
      real a(im,jm,kb)
      real time,scala
      integer ko(*)
      integer iskp,jskp,nko
      character label*(*)
      real amx,scale
      integer i,ib,ie,j,jwr,k,iko,cols
C
      if(scala.ge.0.e0) then
        cols=24
      else
        cols=12
      endif
C
      if (scala.lt.0.e0) scale = 1.e0
      if (scala.eq.0.e0) then
        amx=1.e-12
        do iko=1,nko
          k=ko(iko)
          do j=1,jm,jskp
            do i=1,im,iskp
              amx=max(abs(a(i,j,k)),amx)
            end do
          end do
        end do
          scale=10.e0**(int(log10(amx)+100.e0)-103)
        endif
      if(scala.gt.0.e0) scale=scala
C
      write(6,1) label
    1 format(1x,a40/)
      write(6,2) time,scale
    2 format(' Time = ',f9.4,' days    multiply all values by ',1pe8.2)
C
      do iko=1,nko
C
        k=ko(iko)
C
        write(6,3) k
    3   format(3x,/' Layer k = ',i2)
C
        do ib=1,im,cols*iskp
C
          ie=ib+(cols-1)*iskp
          if(ie.gt.im) ie=im
C
          if(scala.ge.0.e0) then
            write(6,4) (i,i=ib,ie,iskp)
    4       format(/,2x,24i5,/)
          else
            write(6,5) (i,i=ib,ie,iskp)
    5       format(/,12i10,/)
          endif
C
          do j=1,jm,jskp
            jwr=jm+1-j
            if(scala.ge.0.e0) then
              write(6,6) jwr,(nint(a(i,jwr,k)/scale),i=ib,ie,iskp)
    6         format(1x,i3,24i5)
            else
              write(6,7) jwr,(a(i,jwr,k),i=ib,ie,iskp)
    7         format(1x,i2,12(e10.2))
            endif
          end do
C
          write(6,8)
    8     format(//)
C
        end do
C
      end do
C
      return
      end
C--------------------------------------------------------------------

      subroutine prxz(label,time,a,im,iskp,jm,kb,jo,njo,scala,dt,zz)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Writes vertical section of a 3-D field, in the      *
C *                x- or i-direction .                                 *
C *                                                                    *
C *                label ....... label for output                      *
C *                time ........ time (days)                           *
C *                a(im,jm,kb).. array to be printed                   *
C *                iskp ........ skipping interval for i               *
C *                jo .......... 1-D array of j-indices for output     *
C *                njo ......... number of elements in jo              *
C *                scala ....... < 0 for floating point numbers output *
C *                              0 for integer output, divisor for a   *
C *                                based on magnitudes of |a| values   *
C *                              > 0 for integer output, divisor for a *
C *                                given by scala                      *
C *                dt(im,jm) ... total depth                           *
C *                zz(kb) ...... sigma coordinate                      *
C *                                                                    *
C **********************************************************************
C
      implicit none
C
      integer im,jm,kb
      real a(im,jm,kb),dt(im,jm),zz(kb)
      real time,scala
      integer jo(*)
      integer iskp,njo
      character label*(*)
      real amx,scale
      integer i,ib,ie,j,k,ijo,cols
C
      if(scala.ge.0.e0) then
        cols=24
      else
        cols=12
      endif
C
      if (scala.lt.0.e0) scale = 1.e0
      if (scala.eq.0.e0) then
        amx=1.e-12
        do  k=1,kb
          do ijo=1,njo
            j=jo(ijo)
            do i=1,im,iskp
              amx=max(abs(a(i,j,k)),amx)
            end do
          end do
        end do
          scale=10.e0**(int(log10(amx)+100.e0)-103)
        endif
      if(scala.gt.0.e0) scale=scala
C
      write(6,1) label
    1 format(1x,a40/)
      write(6,2) time,scale
    2 format(' Time = ',f9.5,' days    multiply all values by ',1pe8.2)
C
      do ijo=1,njo
C
        j=jo(ijo)
C
        write(6,3) j
    3   format(3x,/' Section j =',i3)
C
        do ib=1,im,cols*iskp
C
          ie=ib+(cols-1)*iskp
          if(ie.gt.im) ie=im
C
          if(scala.ge.0.e0) then
            write(6,4) (i,i=ib,ie,iskp)
    4       format(/,'    i =  ',24i5,/)
          else
            write(6,5) (i,i=ib,ie,iskp)
    5       format(/,'    i =  ',12i10,/)
          endif
C
          write(6,6) (nint(dt(i,j)),i=ib,ie,iskp)
    6     format(8x,'d =',24i5.0,/,'     z or zz')
C
          do k=1,kb
            if(scala.ge.0.e0) then
              write(6,7) k,zz(k),(nint(a(i,j,k)/scale),i=ib,ie,iskp)
    7         format(1x,i2,2x,f6.3,24i5)
            else
              write(6,8) k,zz(k),(a(i,j,k),i=ib,ie,iskp)
    8         format(1x,i2,2x,f6.3,12(e10.2))
            endif
          end do
C
          write(6,9)
    9     format(//)
C
        end do
C
      end do
C
      return
      end
C---------------------------------------------------------------------

      subroutine pryz(label,time,a,im,jm,jskp,kb,io,nio,scala,dt,zz)
C **********************************************************************
C *                                                                    *
C * FUNCTION    :  Writes vertical section of a 3-D field, in the      *
C *                y- or j-direction.                                  *
C *                                                                    *
C *                label ....... label for output                      *
C *                time ........ time (days)                           *
C *                a(im,jm,kb).. array to be printed                   *
C *                jskp ........ skipping interval for j               *
C *                io .......... 1-D array of i-indices for output     *
C *                nio ......... number of elements in io              *
C *                scala ....... < 0 for floating point numbers output *
C *                              0 for integer output, divisor for a   *
C *                                based on magnitudes of |a| values   *
C *                              > 0 for integer output, divisor for a *
C *                                given by scala                      *
C *                dt(im,jm) ... total depth                           *
C *                zz(kb) ...... sigma coordinate                      *
C *                                                                    *
C **********************************************************************
C
      implicit none
      integer im,jm,kb
      real a(im,jm,kb),dt(im,jm),zz(kb)
      real time,scala
      integer io(*)
      integer jskp,nio
      character label*(*)
      real amx,scale
      integer i,j,jb,je,k,iio,cols
C
      if(scala.ge.0.e0) then
        cols=24
      else
        cols=12
      endif
C
      if (scala.lt.0.e0) scale = 1.e0
      if (scala.eq.0.e0) then
        amx=1.e-12
        do  k=1,kb
          do j=1,jm,jskp
            do iio=1,nio
              i=io(iio)
              amx=max(abs(a(i,j,k)),amx)
            end do
          end do
        end do
          scale=10.e0**(int(log10(amx)+100.e0)-103)
        endif
      if(scala.gt.0.e0) scale=scala
C
      write(6,1) label
    1 format(1x,a40/)
      write(6,2) time,scale
    2 format(' Time = ',f9.5,' days    multiply all values by ',1pe8.2)
C
      do iio=1,nio
C
        i=io(iio)
C
        write(6,3) i
    3   format(3x,/' Section i =',i3)
C
        do jb=1,jm,cols*jskp
C
          je=jb+(cols-1)*jskp
          if(je.gt.jm) je=jm
C
          if(scala.ge.0.e0) then
            write(6,4) (j,j=jb,je,jskp)
    4       format(/,'    j =  ',24i5,/)
          else
            write(6,5) (j,j=jb,je,jskp)
    5       format(/,'    j =  ',12i10,/)
          endif
C
          write(6,6) (nint(dt(i,j)),j=jb,je,jskp)
    6     format(8x,'d =',24i5.0,/,'     z or zz')
C
          do k=1,kb
            if(scala.ge.0.e0) then
              write(6,7) k,zz(k),(nint(a(i,j,k)/scale),j=jb,je,jskp)
    7         format(1x,i2,2x,f6.3,24i5)
            else
              write(6,8) k,zz(k),(a(i,j,k),j=jb,je,jskp)
    8         format(1x,i2,2x,f6.3,12(e10.2))
            endif
          end do
C
          write(6,9)
    9     format(//)
C
        end do
C
      end do
C
      return
      end
C   END OF PROGRAM