!--external mode--
!  u->eastward velocity,v->northward velocity,
!  el->surface elevation
!--this program uses (mks) unit--
      program main
      include 'abiki.R8.h90'
      real*8  :: eps,days,horcon,smoth,prts,prtd,dte2,tck,&
                 cbc0,alpha,vel,elemax,elemin,velmax,velmax1,pai
      integer :: i,j,n,ispdif,iend,iprst,iprint,nck,&
                 imax,jmax,nfile,nout
!--
!*********************************************************************
!   open output files                                                 *
!*********************************************************************
      open(10,file='./ECS-TOPO.DATA',form='unformatted')
      open(50,file='./OUTPUT.DATA',  form='unformatted',access='direct',recl=4*im*jm)     
!*********************************************************************
!   input parameters(mks unit)                                       *
!*********************************************************************
!--external time step(sec)--
      dte=3.d0
!--calculation stop(days)--
      days=20./(24.*12.)
!--data output start time(days)-
      prts=0.
!--data output interval(days)--
      prtd=1./(24.*12.)
!--boundary condition--
      nbc=1
      if(abs(nbc).gt.2) call warning('abs(nbc) must be 1/')
!--smagorinsky parameterization coefficient--
      horcon=0.2
      if(horcon.gt.0.) then
        call notice('HORIZONTAL VISCOSITY: SMAGORINSKY SCHEME')
      elseif(horcon.lt.0.) then
        call notice('HORIZONTAL VISCOSITY: CONSTANT')
      else        
        call warning('horcon=0/')
      end if
!--nonlinear or linear calculation(eps=1:nonlinear, eps=0:linear)--
      eps=1.d0
      if(eps.ne.0.d0.and.eps.ne.1.d0) call warning('eps must be 1 or 0/')
      if(eps.eq.0.) then
        call notice('LINEAR CALCULATION')
      else
        call notice('NON-LINEAR CALCULATION')
      end if
!--diffusion+advection term calculation interval--
      ispdif=5
      if(ispdif.eq.0) call warning('ispdif=0/')
!--bottom friction coefficient--
      cbc0=0.0025d0
!--gravitaional accelaration--
      grav=9.806
!--constant in time smoother--
      smoth=0.10d0
!-----------------------------
      dte2=dte*2
      iend  = days*(24*3600)/dte+2
      iprint= prtd*(24*3600)/dte
      iprst = prts*(24*3600)/dte
!------------------------------
      write(6,7030) days,dte,prts,prtd
 7030 format(/'days(days)=',f15.7,/ &
              'dte(sec)    =',f15.7,/ &
              'prts (days)=',f15.3,'    prtd(days)=',f15.3)
      write(6,7031) grav,ispdif,smoth,horcon,cbc0
 7031 format(/'grav  =',f15.7,/ &
              'ispdif=',i7,/    &
              'smoth =',f15.7,/ &
              'horcon=',f15.7,'     cbc=',f13.10/)
!*********************************************************************
!   set grid,depth and coriolis parameter etc.                       *
!*********************************************************************
      call inform1(cbc0)
!*********************************************************************
!   set initial condition                                            *
!*********************************************************************
      call initia(eps)
!**********************************************************************
!                                                                     *
!   !!! begin numerical integration !!!                               *
!                                                                     *
!**********************************************************************
      tck=0
      nck=0
      nout=0
!
!
      do 8000 iext=1,iend
!***************************************************************************************
!***************************************************************************************
      time=dte*dble(iext-1)/(24.*3600.d0) !days
!-------------------------------------------------------------------     
      if(time.ge.tck-0.00001) then
       tck=tck+1./(24.*60.)
       nck=nck+1
       elemax=-1.d+10
       elemin= 1.d+10
       velmax=0.0
       do j=1,jm
       do i=1,im
        if(el(i,j).gt.elemax) elemax=el(i,j)
        if(el(i,j).lt.elemin) elemin=el(i,j)
        vel=dsqrt(u(i,j)**2+v(i,j)**2)
        if(vel.gt.velmax) velmax=vel
       end do
       end do
        write(6,667) iext,time,velmax,elemax,elemin
  667   format('iext=',I10,' time(days)=',f8.2,&
               ' vel-max=',e21.15,' ele-max=',e21.15,' ele-min=',e21.15)
      end if
!--output field data--
      if(iext.ge.iprst+1.and.mod(iext,iprint).eq.1) then
        nout=nout+1
        call backup1(nout)
      end if
!--surface elevation--
      do j=1,jm+1
      do i=1,im+1
        xflux(i,j)=0.25d0*(d(i,j)+d(i-1,j))*(dy(i,j)+dy(i-1,j))*u(i,j)
        yflux(i,j)=0.25d0*(d(i,j)+d(i,j-1))*(dx(i,j)+dx(i,j-1))*v(i,j)
      end do
      end do
      do j=1,jm
      do i=1,im
       elf(i,j)=elb(i,j)-dte2*(xflux(i+1,j)-xflux(i,j)+yflux(i,j+1)-yflux(i,j))/art(i,j)
       elf(i,j)=elf(i,j)*spng(i,j)
      end do
      end do
      call bcond(1)
!============
      call edge(elf)
!============
      do j=0,jm+1
      do i=0,im+1
       df(i,j)=h(i,j)+eps*elf(i,j)
      end do
      end do
!--horizontal diffusion+advection and bottom friction--
      if(mod(iext,ispdif).eq.1) call advdiff(eps,horcon)
!--external u & v--
      alpha=0.225d0
      do j=1,jm
      do i=1,im
       uf(i,j)=-aru(i,j)*0.25d0*( cor(i  ,j)*d(i  ,j)*(v(i  ,j+1)+v(i  ,j)) &
                                 +cor(i-1,j)*d(i-1,j)*(v(i-1,j+1)+v(i-1,j)) ) &
               +0.25d0*grav*(dy(i,j)+dy(i-1,j))*(d(i,j)+d(i-1,j)) &
               *(  (1.d0-2.d0*alpha)*(  el(i,j)- el(i-1,j)                     )   &
                              +alpha*( elb(i,j)-elb(i-1,j)+elf(i,j)-elf(i-1,j) )  ) &
               +aru(i,j)*(wusurf(i,j)-wubot(i,j))+advx(i,j)
      end do
      end do
      do j=1,jm
      do i=1,im
       uf(i,j)=( (db(i,j)+db(i-1,j))*aru(i,j)*ub(i,j)-4.d0*dte*uf(i,j) )/( (df(i,j)+df(i-1,j))*aru(i,j) )
       uf(i,j)=uf(i,j)*spng(i,j)
      end do
      end do
      do j=1,jm
      do i=1,im
       vf(i,j)=+arv(i,j)*0.25d0*( cor(i,j  )*d(i,j  )*(u(i+1,j  )+u(i,j  )) &
                                 +cor(i,j-1)*d(i,j-1)*(u(i+1,j-1)+u(i,j-1)) ) &
               +0.25d0*grav*(dx(i,j)+dx(i,j-1))*(d(i,j)+d(i,j-1))  &
               *(   (1.d0-2.d0*alpha)*(  el(i,j)- el(i,j-1)                     )    &
                               +alpha*( elb(i,j)-elb(i,j-1)+elf(i,j)-elf(i,j-1) )  ) &
               +arv(i,j)*(wvsurf(i,j)-wvbot(i,j))+advy(i,j)
      end do
      end do
      do j=1,jm
      do i=1,im
       vf(i,j)=( (db(i,j)+db(i,j-1))*vb(i,j)*arv(i,j)-4.d0*dte*vf(i,j) )/( (df(i,j)+df(i,j-1))*arv(i,j) )
       vf(i,j)=vf(i,j)*spng(i,j)
      end do
      end do
      call bcond(2)
!============
      call edge(uf); call edge(vf)
!============
!--test for cfl violation. if so,print and stop--
      if(mod(iext,ispdif).eq.1) then
        velmax1=100.
        velmax=0.
        do j=1,jm
        do i=1,im
          vel=dsqrt(vf(i,j)**2+uf(i,j)**2)
          if(vel.ge.velmax) then
            velmax=vel
            imax=i
            jmax=j
          endif
        end do
        end do
        if(velmax.gt.velmax1) then
          write(6,'(''*********** abnormal job end *****************'')')
          write(6,'(''************ user terminated *****************'')')
          write(6,'(''   vamax ='',e12.3,'' imax,jmax ='',2i7)') velmax,imax,jmax
          write(6,'(''longitude='',f12.3,'' latitude='',f12.3)') xdeg(imax),ydeg(jmax)  
          write(6,'(''iext ='',i10)') iext
          goto 9999
        endif
      end if
!--apply filter to remove time split--
!--reset time sequence.--
      do j=0,jm+1
      do i=0,im+1
         u(i,j)= u(i,j)+0.5d0*smoth*( ub(i,j)-2.d0* u(i,j)+ uf(i,j))
         v(i,j)= v(i,j)+0.5d0*smoth*( vb(i,j)-2.d0* v(i,j)+ vf(i,j))
        el(i,j)=el(i,j)+0.5d0*smoth*(elb(i,j)-2.d0*el(i,j)+elf(i,j))
      end do
      end do
      do j=0,jm+1
      do i=0,im+1
        elb(i,j)=el(i,j)
        db(i,j) =h(i,j)+eps*elb(i,j)
        ub(i,j) =u(i,j)
        vb(i,j) =v(i,j)
      end do
      end do
      do j=0,jm+1
      do i=0,im+1
        el(i,j)=elf(i,j)
        d(i,j) =h(i,j)+eps*el(i,j)
        u(i,j) =uf(i,j)
        v(i,j) =vf(i,j)
      end do
      end do
!***************************************************************************************
!***************************************************************************************
 8000 continue
      write(6,*) '***end of time loop***'
!***********************************************************************
!   end of numerical integration                                       *
!***********************************************************************
 9999 continue
      stop
      end
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine backup1(nout)
      include 'abiki.R8.h90'
      integer,intent(inout) :: nout
      integer :: i,j
      write(50,rec=nout*3-2) ((real(el(i,j)),i=1,im),j=1,jm)
      write(50,rec=nout*3-1) ((real( 0.5*(u(i,j)+u(i+1,j  )) ),i=1,im),j=1,jm)
      write(50,rec=nout*3  ) ((real( 0.5*(v(i,j)+v(i  ,j+1)) ),i=1,im),j=1,jm)
      write(6,*) '*********************************************************'
      write(6,*) '  output field data, nout=',nout,' time(days)=',time
      write(6,*) '*********************************************************'
      write(6,*) ' '      
      return
      end
!
!
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine bcond(idx)
      include 'abiki.R8.h90'
      integer,intent(in) :: idx  
      integer :: i,j,n
!-
!-- nbc < 0 : East- West- Boundaries => Cyclic Boundaries-- 
!--abs(nbc)=1:solid boundary condition--
      if(idx.eq.1.and.iext.eq.1) then
        write(6,*) ' '
        if(nbc.lt.0)      call notice('cyclic b.c. for east- west- boundaries') 
        if(abs(nbc).eq.1) call notice('solid wall boundary')
        if(abs(nbc).eq.2) call notice('open boundary condition')
        write(6,*) ' '
      end if
!------------------------------------
      go to (10,20), idx
!--external mode boundary condition--
!--elevation--
   10 continue
      if(abs(nbc).eq.1) goto 11
      if(abs(nbc).eq.2) goto 22
      call warning('check nbc/')
!--solid wall boundary--
   11 continue
!-- North- South- Boundaries--
      do i=1,im
        elf(i, 1)=elf(i, 2)
      end do
      do i=1,im
        elf(i,jm)=elf(i,jm-1)
      end do
!-- East- West- Boundaries--
      if(nbc.lt.0) goto 99
      do j=1,jm
        elf( 1,j)=elf(   2,j)
      end do
      do j=1,jm
        elf(im,j)=elf(im-1,j)
      end do
      goto 99
!--open bounadry condition--
   22 continue
!-- North- South- Boundaries--
      do i=1,im
        elf(i, 1)=0.
      end do
      do i=1,im
        elf(i,jm)=0.
      end do
!-- East- West- Boundaries--
      if(nbc.lt.0) goto 99
      do j=1,jm
        elf( 1,j)=0.
      end do
      do j=1,jm
        elf(im,j)=0.
      end do
      goto 99      
!--set elf=0 at land--
   99 continue
      do j=1,jm
      do i=1,im
        elf(i,j)=elf(i,j)*fsm(i,j)
      end do
      end do
      return
!--external mode velocities--
   20 continue
      if(abs(nbc).eq.1) goto 111
      if(abs(nbc).eq.2) goto 222
!--solid wall boundary--
  111 continue
!-- North- South- Boundaries--
      do i=1,im
        vf(i, 2)=0.d0
        vf(i, 1)=0.d0
        uf(i, 1)=0.d0
      end do
      do i=1,im
        vf(i,jm)=0.d0
        uf(i,jm)=0.d0
      end do
!-- East- West- Boundaries--
      if(nbc.lt.0) goto 999
      do j=1,jm
        uf( 2,j)=0.d0
        uf( 1,j)=0.d0
        vf( 1,j)=0.d0
      end do
      do j=1,jm
        uf(im,j)=0.d0
        vf(im,j)=0.d0
      end do
      goto 999
!--open boundary--
  222 continue
!-- North- South- Boundaries--
      do i=1,im
        vf(i, 2)=vf(i,3)
        vf(i, 1)=vf(i,3)
        uf(i, 1)=uf(i,2)
      end do
      do i=1,im
        vf(i,jm)=vf(i,jm-1)
        uf(i,jm)=uf(i,jm-1)
      end do
!-- East- West- Boundaries--
      if(nbc.lt.0) goto 999
      do j=1,jm
        uf( 2,j)=uf(3,j)
        uf( 1,j)=uf(3,j)
        vf( 1,j)=vf(2,j)
      end do
      do j=1,jm
        uf(im,j)=uf(im-1,j)
        vf(im,j)=vf(im-1,j)
      end do
      goto 999
!--set uaf,vaf=0 at land--
  999 continue
      do j=1,jm
      do i=1,im
        uf(i,j)=uf(i,j)*dum(i,j)
        vf(i,j)=vf(i,j)*dvm(i,j)
      end do
      end do
      return
      end
!
!
!
      subroutine edge(a)
      include 'abiki.R8.h90'
      real*8,dimension(0:im+1,0:jm+1),intent(inout) :: a
      integer :: i,j
!-- North- South- Edges--
      do i=1,im
        a(i,0)=a(i,1)
      end do
      do i=1,im
        a(i,jm+1)=a(i,jm)
      end do
!-- East- West- Edges --
      if(nbc.ge.0) then
        do j=0,jm+1
         a(0,j)=a(1,j)
        end do
        do j=0,jm+1
         a(im+1,j)=a(im,j)
        end do
      else  
        do j=0,jm+1
         a(0,j)=a(im,j)   ! Cyclic Bondary Condition    
        end do
        do j=0,jm+1
         a(im+1,j)=a(1,j) ! Cyclic Bondary Condition 
        end do
      end if
      return
      end
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine inform1(cbc0)
      include 'abiki.R8.h90'
      real*8,intent(in) :: cbc0
      real*8 :: pai
      real*8 :: dxdeg,dydeg,rearth,taux,rhoref
      integer :: i,j,imd,jmd
      integer :: idist,idist_spng
      pai=2.*dasin(1.d+0)
!-- input topo data--
      read(10) imd,jmd
        write(6,*) imd,jmd
        if(imd.ne.im) call warning('check im')
        if(jmd.ne.jm) call warning('check jm')
      read(10) (xdeg(i),i=1,im),(ydeg(j),j=1,jm)
      read(10) ((h(i,j),i=1,im),j=1,jm)
!----------------------
      dxdeg=(xdeg(im)-xdeg(1))/dble(im); dydeg=(ydeg(jm)-ydeg(1))/dble(jm)
      xdeg(0)=xdeg(1)-dxdeg; xdeg(im+1)=xdeg(im)+dxdeg
      ydeg(0)=ydeg(1)-dydeg; ydeg(jm+1)=ydeg(jm)+dydeg
--radius of the earth (m)------
      rearth=6400.d0*1000.d0
!--grid size---------------------------------------
      do j=0,jm+1
      do i=1,im
        dx(i,j)=2.d0*pai*rearth*cos(ydeg(j)*pai/180.d0)*dxdeg/360.d0
        dy(i,j)=2.d0*pai*rearth*dydeg/360.d0
      end do
      end do
!----------------------
        write(6,7040) xdeg(1),xdeg(im),ydeg(1),ydeg(jm), &
                      dx(im/2,jm/2)/1000.,dy(im/2,jm/2)/1000.
 7040 format('longitude range=',f5.1,'-',f5.1,/ ' latitude range=',f5.1,'-',f5.1,/ &
             'dx(km)=',f6.2,',dy(km)=',f6.2/)
!============
     call edge(dx); call edge(dy)
!============
!--at land  (h= 0.),set h=1.0--
      do j=1,jm
      do i=1,im
        if(h(i,j).le.0.) h(i,j)=1.
      end do
      end do
!============
      call edge(h)
!============
!--grid area--
      do j=1,jm
      do i=1,im
        art(i,j)=dx(i,j)*dy(i,j)
        aru(i,j)=0.25d0*(dx(i,j)+dx(i-1,j))*(dy(i,j)+dy(i-1,j))
        arv(i,j)=0.25d0*(dx(i,j)+dx(i,j-1))*(dy(i,j)+dy(i,j-1))
      end do
      end do
!============
      call edge(art); call edge(aru); call edge(arv)
!============ 
!--fsm:mask for scalar variables;=0 over land,=1 over water--
!--dum,dvm:mask for u and v--
      do j=1,jm
      do i=1,im
        fsm(i,j)=0.
        if(h(i,j).gt.1.) fsm(i,j)=1.
      end do
      end do
!============
      call edge(fsm)
!============ 
      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)
      end do
      end do
!============
      call edge(dum); call edge(dvm)
!============ 
!--coriolis parameter--
      do j=1,jm
      do i=1,im
        cor(i,j)=2.*(2.*pai/(24.*3600.))*dsin(ydeg(j)*pai/180.)
      end do
      end do
!============
      call edge(cor)
!============ 
!--drag coefficient of bottom stress--
      do j=1,jm
      do i=1,im
        cbc(i,j)=cbc0*fsm(i,j)
      end do
      end do
!============
      call edge(cbc)
!============ 
!--top surface condition--
      do j=1,jm
      do i=1,im
        wusurf(i,j)=-0.3*(1.d0/1000.d0)*cos(pai*(40.d0-ydeg(j))/30.d0)
        wvsurf(i,j)=0.d0
      end do
      end do
!============
      call edge(wusurf); call edge(wvsurf)
!============ 
      idist_spng=20      
      do j=1,jm
      do i=1,im
        idist=min(abs(i-1),abs(im-i),abs(j-1),abs(jm-j))
        if(idist.le.idist_spng) then
          spng(i,j)=dble(idist)/dble(idist_spng)
        else
          spng(i,j)=1.d0
        end if
      end do
      end do
!============
      call edge(spng)
!============ 
      return
      end
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++x++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine initia(eps,cfile)
      include 'abiki.R8.h90'
      real*8,intent(in) :: eps
      integer :: i,j,i1,j1
      character,intent(in) :: cfile*(*)
      real*8 :: dist
!--external mode initial condition--
      do j=0,jm+1
      do i=0,im+1
        ub(i,j)=0.
        vb(i,j)=0.
        dist=sqrt( (xdeg(i)-xdeg(im/2))**2+(ydeg(j)-ydeg(jm/2))**2 )
        elb(i,j)=exp(-(dist/0.5d0)**2)
        aam(i,j)=500.
      end do
      end do 
      do j=0,jm+1
      do i=0,im+1
          u(i,j)=ub(i,j)
         uf(i,j)= u(i,j)
          v(i,j)=vb(i,j)
         vf(i,j)= v(i,j)
         el(i,j)=elb(i,j)
        elf(i,j)=el(i,j)
      end do
      end do
      do j=0,jm+1
      do i=0,im+1
         d(i,j)=h(i,j) +eps*el(i,j)
        db(i,j)=h(i,j)+eps*elb(i,j)
      end do
      end do
      return
      end
!
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine advdiff(eps,horcon)
      include 'abiki.R8.h90'
      real*8,intent(in) :: eps,horcon
      integer :: i,j
!--horizontal diffusivity:smagorinski parameterization--
      if(horcon.gt.0.) then
        do j=1,jm
        do i=1,im
          aam(i,j)=horcon*dx(i,j)*dy(i,j)  &
                  *dsqrt(  ( (u(i+1,j)-u(i,j))/dx(i,j) )**2 & 
                          +( (v(i,j+1)-v(i,j))/dy(i,j) )**2  &
               +0.5d0*( 0.25d0*( (u(i,j+1)+u(i+1,j+1))-(u(i,j-1)+u(i+1,j-1)) )/dy(i,j)        &
                       +0.25d0*( (v(i+1,j)+v(i+1,j+1))-(v(i-1,j)+v(i-1,j+1)) )/dx(i,j) )**2  )
        end do
        end do
      else
        do j=1,jm
        do i=1,im
          aam(i,j)=abs(horcon)
        end do
        end do
      end if        
!============
      call edge(aam)
!============
!--calculate u-advection & diffusion for external mode--
!--advective fluxes--
      do j=0,jm+1
      do i=0,im+1
        advx(i,j)=0.0
        xflux(i,j)=0.0
        yflux(i,j)=0.0
      end do
      end do
!--
      if(eps.ne.0.) then
        do j=1,jm
        do i=1,im
          xflux(i,j)=0.125d0*(  (d(i+1,j)+d(i  ,j))*u(i+1,j) &
                               +(d(i  ,j)+d(i-1,j))*u(i  ,j) )*(u(i+1,j)+u(i,j))
        end do
        end do
        do j=1,jm
        do i=1,im
          yflux(i,j)=0.125d0*(  (d(i  ,j)+d(i  ,j-1))*v(i  ,j) &
                               +(d(i-1,j)+d(i-1,j-1))*v(i-1,j) )*(u(i,j)+u(i,j-1))
        end do
        end do
      end if
!--add viscous fluxes--
      do j=1,jm
      do i=1,im
        xflux(i,j)=xflux(i,j)-d(i,j)*2.d0*aam(i,j)*(ub(i+1,j)-ub(i,j))/dx(i,j)
        xflux(i,j)=xflux(i,j)*dy(i,j)
      end do
      end do 
      do j=1,jm
      do i=1,im
        tps(i,j)=0.25d0*(  d(i,j)+  d(i,j-1)+  d(i-1,j)  +d(i-1,j-1))   &
                       *(aam(i,j)+aam(i,j-1)+aam(i-1,j)+aam(i-1,j-1))   &
                       *( (ub(i,j)-ub(i,j-1))/(dy(i,j)+dy(i-1,j)+dy(i,j-1)+dy(i-1,j-1))   &
                         +(vb(i,j)-vb(i-1,j))/(dx(i,j)+dx(i-1,j)+dx(i,j-1)+dx(i-1,j-1)) )
      end do
      end do
      do j=1,jm
      do i=1,im
        yflux(i,j)=(yflux(i,j)-tps(i,j))*0.25d0*(dx(i,j)+dx(i-1,j)+dx(i,j-1)+dx(i-1,j-1))
      end do
      end do
!============
      call edge(xflux); call edge(yflux)
!============
      do  j=1,jm
      do  i=1,im
        advx(i,j)=(xflux(i,j)-xflux(i-1,j))+(yflux(i,j+1)-yflux(i,j))
      end do
      end do
!--calculate v-advection & diffusion for external mode--
!--advective fluxes--
      do j=1,jm
      do i=1,im
        advy(i,j)=0.0
        xflux(i,j)=0.0
        yflux(i,j)=0.0
      end do
      end do
!--   
      if(eps.ne.0.) then
        do j=1,jm
        do i=1,im
          xflux(i,j)=0.125d0*( (d(i,j  )+d(i-1,j  ))*u(i,j  )        &
                              +(d(i,j-1)+d(i-1,j-1))*u(i,j-1))*(v(i-1,j)+v(i,j))
        end do
        end do
        do j=1,jm
        do i=1,im
          yflux(i,j)=0.125d0*( (d(i,j+1)+d(i,j  ))*v(i,j+1)     &
                              +(d(i,j  )+d(i,j-1))*v(i,j  ))*(v(i,j+1)+v(i,j))
        end do
        end do
      end if
!--add viscous fluxes--
      do j=1,jm
      do i=1,im
        yflux(i,j)=yflux(i,j)-d(i,j)*2.d0*aam(i,j)*(vb(i,j+1)-vb(i,j))/dy(i,j)
        yflux(i,j)=yflux(i,j)*dx(i,j)
      end do
      end do
      do j=1,jm
      do i=1,im
        xflux(i,j)=(xflux(i,j)-tps(i,j))*0.25d0*(dy(i,j)+dy(i-1,j)+dy(i,j-1)+dy(i-1,j-1))
      end do
      end do
!============
      call edge(xflux); call edge(yflux)
!============
      do j=1,jm
      do i=1,im
        advy(i,j)=(xflux(i+1,j)-xflux(i,j))+(yflux(i,j)-yflux(i,j-1))
      end do
      end do
!--add curvature term--
      if(eps.ne.0.) then
        do j=1,jm
        do i=1,im
          curv2d(i,j)=0.25d0*( (v(i,j+1)+v(i,j))*(dy(i+1,j)-dy(i-1,j))     &
                              -(u(i+1,j)+u(i,j))*(dx(i,j+1)-dx(i,j-1)) )/(dx(i,j)*dy(i,j))
        end do
        end do
!============
        call edge(curv2d)
!============
        do j=1,jm
        do i=1,im
          advx(i,j)=advx(i,j)-aru(i,j)*0.25d0*( curv2d(i  ,j)*d(i  ,j)*(v(i  ,j+1)+v(i  ,j))     &
                                               +curv2d(i-1,j)*d(i-1,j)*(v(i-1,j+1)+v(i-1,j)) )
        end do
        end do
        do j=1,jm
        do i=1,im
          advy(i,j)=advy(i,j)+arv(i,j)*0.25d0*( curv2d(i,j  )*d(i,j  )*(u(i+1,j  )+u(i,j  ))     &
                                               +curv2d(i,j-1)*d(i,j-1)*(u(i+1,j-1)+u(i,j-1)) )
        end do
        end do
      end if
!--bottom stress terms and add curvature term--
      do j=1,jm
      do i=1,im
        wubot(i,j)=-0.5d0*(cbc(i,j)+cbc(i-1,j))   &
                  *dsqrt(  ub(i,j)**2 + ( 0.25d0*( vb(i  ,j)+vb(i  ,j+1)     &
                                                  +vb(i-1,j)+vb(i-1,j+1) ) )**2  )*ub(i,j)
      end do
      end do
      do j=1,jm
      do i=1,im
        wvbot(i,j)=-0.5d0*(cbc(i,j)+cbc(i,j-1))     &
                    *dsqrt(  ( 0.25d0*( ub(i,j  )+ub(i+1,j  )    &
                                       +ub(i,j-1)+ub(i+1,j-1) ) )**2+vb(i,j)**2  )*vb(i,j)
      end do
      end do
      return
      end
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine warning(message)
      character :: message*(*)
      write(6,*) message
      stop
      end
!
!
!
      subroutine notice(message)
      implicit none
      character,intent(in) :: message*(*)
      character :: message2*51,blank*1
      integer   :: length,nsp
      blank=' '
      length=len_trim(message)
      if(length.ge.40) length=40
      nsp=(51-4-length)
      if(nsp/2*2.eq.nsp) then
        message2='**'//repeat(blank,nsp/2)//trim(message)//&
                       repeat(blank,nsp/2)//'**'
      else
        message2='**'//repeat(blank,nsp/2)//trim(message)//&
                       repeat(blank,nsp/2)//' **'
      end if        
      write(6,*) ' '
      write(6,'(A51)') '***************************************************'
      write(6,'(A51)') '***************************************************'
      write(6,'(A51)') '**                                               **'
      write(6,'(A51)') message2
      write(6,'(A51)') '**                                               **'
      write(6,'(A51)') '***************************************************'
      write(6,'(A51)') '***************************************************'
      write(6,*) ' '
      return
      end
