!--Output Topo. Data Size--
      implicit none
      integer,parameter :: im=451,jm=301
!--Input Topo. Data Size--
      integer,parameter :: kq=2*60*360/60+60*360+1,lq=60*180+1
!--
      integer :: itopo0(kq,lq)
      real*8  :: xd0(kq),yd0(lq)
      integer :: ks(im),ls(jm)
      real*8 :: topo(im,jm),xd(0:im+1),yd(0:jm+1),dx(im,jm),dy(im,jm)
      real*8 :: pai,re,xds,xde,yds,yde,dxd,dyd,topomin, &
                dxd0,dyd0,dxd1,dyd1,x00,x11,y00,y11,t00,t01,t10,t11
      integer :: i,j,k,l,k1,l1,kd,ld,im1,jm1,nd,nc,nsmooth,mgx,mgy
      common /main1/itopo0
      common /main2/topo,dx,dy
      pai=2.*dasin(1.0d+0)
!--input file--
      open(10,file='./GEBCO.GLOBAL.BIN.DATA',form='unformatted')
!--output file--
      open(20,file='./ECS-TOPO.DATA'  ,form='unformatted'                             )      
      open(21,file='./ECS-TOPO.DIRECT',form='unformatted',access='direct',recl=4*im*jm)
!--NOTE:this program uses MKS unit--
!--Radius of Earth(m)--
      re=6370.*1000.
!--Range of the Output Topo. Data--
      xds= 125.0
      xde= 132.5
      yds=  30.0
      yde=  35.0
!---Latitudinal and Longitudinal Grid Interval (deg)---
      dxd=1.d0/60.d0
      dyd=1.d0/60.d0
!--------------------------------    
      im1=idnint((xde-xds)/dxd)+1
      jm1=idnint((yde-yds)/dyd)+1
        write(6,*) 'im=',im1,' jm=',jm1
        if(im1.ne.im) call warning('change parameter im')
        if(jm1.ne.jm) call warning('change parameter jm')
!--Input data--
      read(10) kd,ld
        write(6,*) kd,ld
        if(kd.ne.kq.or.ld.ne.lq) call warning('kq or lq mismatch')
      read(10) xd0,yd0
      read(10) itopo0
      dxd0=abs(xd0(kq/2+1)-xd0(kq/2))
      dyd0=abs(yd0(lq/2+1)-yd0(lq/2))
      do k=1,kq-1
        dxd1=abs(xd0(k+1)-xd0(k))
        if(abs(dxd1-dxd0)/abs(dxd0).gt.1.d-5) call warning('dxd0 is not constatnt')
      end do
      do l=1,lq-1
        dyd1=abs(yd0(l+1)-yd0(l))
        if(abs(dyd1-dyd0)/abs(dyd0).gt.1.d-5) call warning('dyd0 is not constatnt')
      end do
        write(6,*) '**** GEBCO DATA *****'
        write(6,*) 'long.=',xd0(1),'-',xd0(kq),' 1/dlong.=',1.d0/(xd0(kq/2+1)-xd0(kq/2))
        write(6,*) 'lat. =',yd0(1),'-',yd0(lq),' 1/dlat. =',1.d0/(yd0(lq/2+1)-yd0(lq/2))
!------------------------------------
      do i=1,im
        xd(i)=xds+(xde-xds)*dble(i-1)/dble(im-1)
      end do
      do j=1,jm
        yd(j)=yds+(yde-yds)*dble(j-1)/dble(jm-1)
      end do
!-------------------------------------------------------
      do j=1,jm
      do i=1,im
       dx(i,j)=re*(dxd*pai/180.)*dcos(yd(j)*pai/180.)
       dy(i,j)=re*(dyd*pai/180.)
      end do
      end do 
      write(6,*) 
      write(6,*) '**** OUTPUT  DATA *****'
      write(6,*) 'im=',im,',jm=',jm
      write(6,*) 'long.=',xd(1),'-',xd(im),' 1/dlong.=',1.d0/(xd(im/2+1)-xd(im/2))
      write(6,*) 'lat. =',yd(1),'-',yd(jm),' 1/dlat. =',1.d0/(yd(jm/2+1)-yd(jm/2))
      write(6,*) 'long. range=',real(xd(im)-xd(1))
      write(6,*) 'lat.  range=',real(yd(jm)-yd(1))
      write(6,*) 'dx,dy(im/2,jm/2)(km)=',real(dx(im/2,jm/2)/1000.),',',real(dy(im/2,jm/2)/1000.)
      write(6,*) 
!----------------------------------------------------------------------------------------
      do i=1,im
        ks(i)=0
        do k=1,kq-1
          if(xd(i).ge.xd0(k).and.xd(i).lt.xd0(k+1)) then
             ks(i)=k
          end if
        end do
        if(ks(i).eq.0) call warning('check ks')
      end do
      do j=1,jm
        ls(j)=0
        do l=1,lq-1
           if(yd(j).ge.yd0(l).and.yd(j).lt.yd0(l+1)) then
             ls(j)=l
           end if
        end do
        if(ls(j).eq.0) call warning('check ls')
      end do
!---------------------------------------------------------------------------------------- 
      do j=1,jm
      do i=1,im
         x00=xd0(ks(i)  )
         x11=xd0(ks(i)+1)
         y00=yd0(ls(j)  )
         y11=yd0(ls(j)+1)
         t00=itopo0(ks(i)  ,ls(j)  )
         t10=itopo0(ks(i)+1,ls(j)  )
         t01=itopo0(ks(i)  ,ls(j)+1)
         t11=itopo0(ks(i)+1,ls(j)+1)
         topo(i,j)=( ( t00*(y11-yd(j))+t01*(yd(j)-y00) )*(  x11-xd(i) )    &
                    +( t10*(y11-yd(j))+t11*(yd(j)-y00) )*(xd(i)-x00   ) )  &
                  /( (x11-x00)*(y11-y00) )
         topo(i,j)=-topo(i,j)
      end do
      end do     
      write(6,*) 'End of Topography Data Interpolation'
!--Smoothing of the Sea Topography--
      nsmooth=1
      call smooth(im,jm,topo,nsmooth)
!---At land( topo=<0) => topo=0.--
      do j=1,jm
      do i=1,im
        if(topo(i,j).le.0.) topo(i,j)=0.
      end do
      end do
!--set no gradient boundary depth--
      do i=1,im
       if(topo(i, 1).ne.0.) topo(i, 1)=topo(i,   2)
       if(topo(i,jm).ne.0.) topo(i,jm)=topo(i,jm-1)
      end do 
      do j=1,jm
       if(topo( 1,j).ne.0.) topo( 1,j)=topo(   2,j)
       if(topo(im,j).ne.0.) topo(im,j)=topo(im-1,j)
      end do 
       write(6,*) 'No gradient boundary depth' 
!--Minimum depth--
      topomin=5.
      nc=0
      do j=1,jm
      do i=1,im
        if(topo(i,j).gt.0.and.topo(i,j).lt.topomin) then
          topo(i,j)=topomin
          nc=nc+1
        end if
      end do
      end do
       write(6,*) ' '  
       write(6,*) 'Minimum depth(m)=',topomin
       write(6,*) 'Change points=',nc,'/',im*jm
       write(6,*) ' '
!--Output data--
!--double precision data output--
      write(20) im,jm
      write(20) (xd(i),i=1,im),(yd(j),j=1,jm)
      write(20) ((topo(i,j),i=1,im),j=1,jm)
!------------------------      
      write(21,rec=1) ((real(topo(i,j)),i=1,im),j=1,jm)
!--
      write(6,*) 'Topo -> (m)unit'
      write(6,*) 'AT Ocean>0., AT Land=0., AT Masked Ocean=-1.'
      write(6,*) 'end of output data'
      stop
      end
!
!*******************************************************************************
      subroutine nearest1(x,it,is,ie,x0)
      implicit none
      integer,intent(in)  :: is,ie
      real*8 ,intent(in)  :: x0(is:ie)
      real*8 ,intent(in)  :: x
      integer,intent(out) :: it
      real*8  :: x0min,x0max,rmin,r
      integer :: i 
      x0min=dmin1(x0(is),x0(ie))
      x0max=dmax1(x0(is),x0(ie))
      if(x.lt.x0min.or.x.gt.x0max) call warning('stop at sub. nearest/')
      rmin=abs(x-x0(1))
      it=1
      do 10 i=is,ie      
       r=abs(x-x0(i))
       if(r.lt.rmin) then
        rmin=r
        it=i
       end if
   10 continue    
      return
      end 
!********************************************************************* 
      subroutine smooth(iq,jq,s,nsmoth)
      implicit none
      integer,intent(in)    :: iq,jq,nsmoth
      real*8 ,intent(inout) :: s(iq,jq)
      real*8 :: sd(6000,3000),eps,s1,s2,s3,s4
      integer :: i,j,n
      common /sub_smooth/ sd
      if(iq.gt.6000.or.jq.gt.3000) call warning('iq>6000 or jq> 3000/')
!--Laplacian smoothing--
!     s(i,j)=s(i,j)+eps*( (s(i+1,j)+s(i-1,j)-2*s(i,j))
!                        +(s(i,j+1)+s(i,j-1)-2*s(i,j)) )
!-----------------------
      eps=0.125      
      do n=1,nsmoth
        do j=1,jq
        do i=1,iq
          s1=0.d0
          s2=0.d0
          s3=0.d0
          s4=0.d0
          if(i+1.le.iq) then 
            if(s(i+1,j).gt.0.) s1=s(i+1,j)-s( i,j)
          end if
          if(i-1.ge.1)  then
            if(s(i-1,j).gt.0.) s2=s(i-1,j)-s(i,j)
          end if
          if(j+1.le.jq) then
            if(s(i,j+1).gt.0.) s3=s(i,j+1)-s(i,j)
          end if
          if(j-1.ge.1) then
            if(s(i,j-1).gt.0.) s4=s(i,j-1)-s(i,j)
          end if
          if(s(i,j).gt.0.) sd(i,j)=s(i,j)+eps*(s1+s2+s3+s4)
        end do
        end do
        do j=1,jq
        do i=1,iq
          s(i,j)=sd(i,j)
        end do
        end do    
      end do 
      write(6,*) '------------------------------------------------'
      write(6,*) 'End of Laplacian Smoothing to Sea Topo., nsmooth=',nsmoth
      write(6,*) '------------------------------------------------'
      return
      end
!********************************************************************************
      subroutine smooth3(xs,xe,ys,ye,iq,jq,x,y,s,nsmoth)
      implicit none
      real*8, intent(in)    :: xs,xe,ys,ye,x(iq),y(jq)
      integer,intent(in)    :: iq,jq,nsmoth
      real*8, intent(inout) :: s(iq,jq)
      real*8 :: sd(6000,3000),s1,s2,s3,s4,eps
      integer :: i,j,n
      common /sub_smooth3/ sd
      if(iq.gt.6000.or.jq.gt.3000) call warning('iq > 6000 or jq> 3000/')
!--Laplacian smoothing--
!     s(i,j)=s(i,j)+eps*( (s(i+1,j)+s(i-1,j)-2*s(i,j))
!                        +(s(i,j+1)+s(i,j-1)-2*s(i,j)) )
!-----------------------
      if(xs.lt.x(1).or.xs.gt.x(iq)) call warning('xs is out of range')
      if(xe.lt.x(1).or.xe.gt.x(iq)) call warning('xe is out of range')
      if(ys.lt.y(1).or.ys.gt.y(jq)) call warning('ys is out of range')
      if(ye.lt.y(1).or.ye.gt.y(jq)) call warning('ye is out of range')
      eps=0.125d0      
      do n=1,nsmoth
        do j=1,jq
        do i=1,iq
          s1=0.d0
          s2=0.d0
          s3=0.d0
          s4=0.d0
          if(x(i).ge.xs.and.x(i).le.xe.and.&
             y(j).ge.ys.and.y(j).le.ye) then
            if(i+1.le.iq) then 
              if(s(i+1,j).gt.0.) s1=s(i+1,j)-s(i,j)
            end if
            if(i-1.ge.1)  then
              if(s(i-1,j).gt.0.) s2=s(i-1,j)-s(i,j)
            end if
            if(j+1.le.jq) then
              if(s(i,j+1).gt.0.) s3=s(i,j+1)-s(i,j)
            end if
            if(j-1.ge.1) then
              if(s(i,j-1).gt.0.) s4=s(i,j-1)-s(i,j)
            end if
          end if
          if(s(i,j).gt.0.) sd(i,j)=s(i,j)+eps*(s1+s2+s3+s4)
        end do
        end do
        do j=1,jq
        do i=1,iq
          s(i,j)=sd(i,j)
        end do
        end do     
      end do
      write(6,*) '-------------------------------------------------'
      write(6,*) 'End of Laplacian Smoothing to Sea Topo., nsmooth=',nsmoth
      write(6,'(A3,F10.5,A3,F10.5)') ' x=',xs,' - ',xe
      write(6,'(A3,F10.5,A3,F10.5)') ' y=',ys,' - ',ye
      write(6,*) '--------------------------------------------------'
      return
      end
!******************************************************
      subroutine warning(message)
      character,intent(in) :: message*(*)
      write(6,*) message
      stop
      end
