c GTOPO30データを切り出すプログラム
c 2001/01/29 by A.Tsubota
c コンパイル時 要-byteswapio

c 変数定義

       parameter(im=5000,jm=6000)
       parameter(nlo=4800,nla=6000,nle=nlo*2)

       integer*2 ivv(nlo),iv(nlo,nla)
       dimension map(im,jm)
       real*8 fblon,felon,fblat,felat
       character flnam(40)*80

c 必要緯度経度読み込み

       open(3,file='init.dat')
       read(3,*) idate,itime,fblon,felon,felat,fblat
       close(3)

c 

       open(4,file='src/fname.in',status='old')
       do jfile=3,1,-1
          do ifile=1,9
             mt=10*jfile+ifile
             read(4,'(a80)') flnam(mt)
          end do
       end do

c 読み込み

       dd=1.0/120.0
       ddd=1.0/240.0

       fw=fblon+ddd
       fe=felon-ddd
       fn=fblat-ddd
       fs=felat+ddd

       iw=int(fblon+180)*120
       ie=int(felon+180)*120
       in=int(90-fblat)*120
       is=int(90-felat)*120
       ifilew=int(fblon+180)/40+1
       ifilee=int(felon+180)/40+1
       jfilen=int(90-fblat)/60+1
       jfiles=int(90-felat)/60+1

       if (mod(int(felon+180),40).eq.0) then
          ifilee=ifilee-1
       end if

c       write(*,*)iw,ie,in,is
c       write(*,*)ifilew,ifilee,jfilen,jfiles

       do jfile=jfilen,jfiles
          do ifile=ifilew,ifilee
             mt=10*jfile+ifile
             open(mt,file=flnam(mt),access='direct',recl=nle)
c             write(*,*) flnam(mt)
          end do
       end do

       jo=0
       jfile=jfilen

 9000  continue

       do ifile=ifilew,ifilee
          mt=10*jfile+ifile
          do j=1,nla
          read(mt,rec=j) ivv
          do i=1,nlo
             iv(i,j)=ivv(i)
          end do
          end do
          if(jfile.gt.jfilen) then
             jn1=1
             else
                jn1=in-(jfilen-1)*nla+1
                jnshift=jn1
          end if
          if(jfile.ne.jfiles) then
             js1=nla
             else
                js1=is-(jfiles-1)*nla
          end if
          if(ifile.gt.ifilew) then
             iw1=1
             else
                iw1=iw-(ifilew-1)*nlo+1
                iwshift=iw1
          end if
          if(ifile.ne.ifilee) then
             ie1=nlo
             else
                ie1=ie-(ifilee-1)*nlo
          end if
c          write(*,*) iw1,ie1,jn1,js1
          do j=jn1,js1
             jo=(jfile-jfilen)*nla+j-jnshift+1
             do i=iw1,ie1
                io=(ifile-ifilew)*nlo+i-iwshift+1
                map(io,jo)=iv(i,j)
c                write(*,*) io,jo,i,j,iv(i,j)
             end do
          end do
       imo=imo+ie1-iw1+1
       end do

       if(jfile.ne.jfiles) then
       jfile=jfile+1
       go to 9000
       end if

       open(52,file='topo/gtopo30.xyz',status='new')
c       write(*,*) is-in,ie-iw
       do j=1,is-in
          do i=1,ie-iw
             alon=fblon+(i-1.)*dd
             alat=fblat-(j-1.)*dd
             write(52,'(2f10.5,2x,i7)') alon,alat,map(i,j)
c             write(*,'(2f10.5,2x,i7)') alon,alat,map(i,j)
          end do
       end do

       close(52)

       stop
       end
 
