      program extract_xyz
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   program to extract subregions in GRAVSOFT ASCII format into
c   single line format having lat,lon,fieldv values on each line
c   Only data WITHIN user specified region will be output
c   Se text below  OA 04-2008   
c      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      implicit real*8 (a-h,o-z)
      character*56 fname,oname,ename
      real*4 grow(23000),erow(23000)
      igrow =     23000

      write(*,*)'Extract and Convert 1 minute DTU13 gravsoft file'
      write(*,*)'Into GMT ascii format or AsCII grid file format'
      write(*,*)'Extract two fields simultaneously (i.e. grav,err)'
      write(*,*)'-----------------------------------------------'
      write(*,*)'For interpolation and averating please use geoip'
      write(*,*)' ' 
      write(*,*)'Program uses ASCII files which takes long time '
      write(*,*)'to process but works across platforms ' 
      write(*,*)'Much faster execution time can be achieved by ' 
      write(*,*)'compiling gbin and executing gbin.job program ' 
      write(*,*)'Subsequently xtract use BinExtract program '
      write(*,*)'-----------------------------------------------'
      write(*,*)' ' 
 8000 write(*,*) 'Enter First Field to extract: '
      write(*,*) ' 1: DTU13 1 minute Gravity field '
      write(*,*) ' 2: DTU13 1 minute Mean Sea Surface '
      write(*,*) ' 3: DTU13 1 minute Bathymetry'
      write(*,*) ' 4: DTU13 1 minute Mean Dynamic Topography'
      write(*,*) ' 5: DTU13 1 minute Error file '
      read (*,*) ifil
      write(*,*) 'Enter Second Field to extract: '
      write(*,*) ' 1: DTU13 1 minute Gravity field '
      write(*,*) ' 2: DTU13 1 minute Mean Sea Surface '
      write(*,*) ' 3: DTU13 1 minute Bathymetry'
      write(*,*) ' 4: DTU13 1 minute Mean Dynamic Topography'
      write(*,*) ' 5: DTU13 1 minute Error file '
      read (*,*) jfil
      if (ifil.eq.jfil) write(*,*)'Choose different files'
      if (ifil.eq.jfil) goto 8000
      if (ifil.le.1) open(20,file='DTU13GRA_1min.gra')
      if (ifil.eq.2) open(20,file='DTU13MSS_1min.mss')
      if (ifil.eq.3) open(20,file='DTU13BAT_1min.bat')
      if (ifil.eq.4) open(20,file='DTU13MDT_1min.mdt') 
      if (ifil.ge.5) open(20,file='DTU13ERR_1min.err') 
      if (jfil.le.1) open(21,file='DTU13GRA_1min.gra')
      if (jfil.eq.2) open(21,file='DTU13MSS_1min.mss')
      if (jfil.eq.3) open(21,file='DTU13BAT_1min.bat')
      if (jfil.eq.4) open(21,file='DTU13MDT_1min.mdt') 
      if (jfil.ge.5) open(21,file='DTU13ERR_1min.err') 
      read(21,*) rla1,rla2,rlo1,rlo2,dla,dlo 
      read(20,*) rla1,rla2,rlo1,rlo2,dla,dlo
      if (dlat.lt.0.017) dlat = 1.0d0/60.0d0
      if (dlon.lt.0.017) dlon = 1.0d0/60.0d0
      nn = (rla2-rla1)/dla + 1.5
      ne = (rlo2-rlo1)/dlo + 1.5
      write(*,30)
      write(*,*) ' '
30    format(' Input grid file label:')
      write(*,*)' Min/max latitude, min/max long, Spacing(lat,long)'
      write(*,31) rla1,rla2,rlo1,rlo2,dla,dlo
31    format(' ',4f11.6,2f12.7)
      write(*,33) nn,ne,nn*ne
33    format(' Number of points in grid (north,east,total): ',2i7,i12)
      if (ne.gt.igrow) stop 'rows in grid too long - change igrow' 
      write(*,*) ' '

c     write(*,*)' Extraction of 1 minute files '
 9000 continue
      write(*,*) ' '
      write(*,*)' Point within the user specified region is extracted'
      write(*,*)' Enter area limits(lat min,lat max,lon min,lon max)'
      write(*,*)' Enter 0,0,0,0 if you want the entire grid output '
      read(*,*) ala1,ala2,alo1,alo2
       
c     defining subregion for extraction
      
      if (abs(ala1).lt.0.01.and.abs(ala2).lt.0.01.and.
     . abs(alo1).lt.0.01.and.abs(alo2).lt.0.01)then 
        ala1 = rla1
        ala2 = rla2
        alo1 = rlo1
        alo2 = rlo2
        write(*,*) ala1,ala2,alo1,alo2
      end if
        la1 = (ala1-rla1)/dla + 1.5
        la2 = (ala2-rla1)/dla + 1.5
        lo1 = (alo1-rlo1)/dlo + 1.5
        lo2 = (alo2-rlo1)/dlo + 1.5
     
      if (la2.lt.la1) goto 9001 
      if (lo2.lt.lo1) goto 9001 
      if (la2.lt.1) goto 9001 
      if (la1.gt.nn) goto 9001

      if (la1.lt.1) la1 = 1
      if (la2.gt.nn) la2 = nn
      if (lo1.lt.1) lo1 = 1
      if (lo2.gt.ne) lo2 = ne
      if (lo2.lt.1) goto 9001 
      if (lo1.gt.ne) goto 9001
      
      if (rla1+(la1-1)*dla.lt.ala1-1.d-5) la1 = la1 + 1
      if (rla1+(la2-1)*dla.gt.ala2+1.d-5) la2 = la2 - 1
      if (rlo1+(lo1-1)*dlo.lt.alo1-1.d-5) lo1 = lo1 + 1
      if (rlo1+(lo2-1)*dlo.gt.alo2+1d-5) lo2 = lo2 - 1
c     correct for profile output 
      if (abs(la1-la2).eq.1) la2 = la1
      if (abs(lo1-lo2).eq.1) lo2 = lo1
      
      goto 9002
 9001 continue     
      write(*,*) '--------------------------------------------'
      write(*,*) 'User specified wrong input area parameters:'
      write(*,*) 'User specified (min,max latitude)  ',ala1,ala2
      write(*,*) 'User specified (min,max longitude) ',alo1,alo2
      write(*,*) 'Input MSS grid (min,max latitude)  ',rla1,rla2
      write(*,*) 'Input MSS grid (min,max longitude) ',rlo1,rlo2
      write(*,*) '--------------------------------------------'
      goto 9000
9002  continue
        
      write(*,*)'Corrected area to be within user and grid limits' 
      write(*,*)'Min/max latitude, min/max longitude '
      write(*,500) rla1+(la1-1)*dla,rla1+(la2-1)*dla,
     . rlo1+(lo1-1)*dlo,rlo1+(lo2-1)*dlo
      write(*,*) la1,la2,lo1,lo2
    
      write(*,*) ' '
      write(*,*) ' Enter name for output file '
      read(*,'(A)') oname
      write(*,*) ' '
      open (30,file = oname,form = 'formatted',err=901)
      write(*,*) 'Enter output file format '
      write(*,*) '1 = GMT ASCII format (long,lat,field1,field2)'
      write(*,*) '2 = GMT ASCII format (lat,long,field,field2)'
      write(*,*) '3 = GRAVSOFT ASCII(ista,lat,long,field1,field2)' 
      write(*,*) '4 = GRAVSOFT ASCII grid format (field1 as grid)' 
      read(*,*) iform
      if (iform.lt.1.or.iform.gt.4) iform = 1
      if (iform.eq.4)  write(30,501) 
     . rla1+(la1-1)*dla,rla1+(la2-1)*dla,
     . rlo1+(lo1-1)*dlo,rlo1+(lo2-1)*dlo,dla,dlo
 500  format(6f11.4)        
 501  format(6f13.7)
      write(*,*)' Extraction of very large files '
      write(*,*)' Please wait '
      itot = (nn-la1)
      do i = nn,la2+1,-1
	if (mod(nn-i,500).eq.1) 
     .     write(*,*)nn-i,' of ',itot,' processed'
        read (20,*) (grow(j),j=1,ne)
        read (21,*) (erow(j),j=1,ne)
	end do
      do i = la2,la1,-1
	if (mod(nn-i,500).eq.1) 
     .    write(*,*)nn-i,' of ',itot,' processed'
        read (20,*) (grow(j),j=1,ne)
        read (21,*) (erow(j),j=1,ne)
        a1 =  rla1 + (i-1)*dla 
	do j = lo1,lo2
         a2 =  rlo1 + (j-1)*dlo
	 if (iform.eq.1) write(30,555) a2,a1,grow(j),erow(j)
	 if (iform.eq.2) write(30,555) a1,a2,grow(j),erow(j)
	 if (iform.eq.3) write(30,556) i,a2,a1,grow(j),erow(j)
	 if (iform.eq.4) write(30,557) grow(j)
	 
	 end do
       end do
 555   format(2f11.5,3f11.3)
 556   format(i6,2f11.5,3f11.3)
 557   format(f11.3)
       goto 902
 900   write(*,*) 'Ensure inputfile is in current directory'
       write(*,*) fname
       goto 902
 901   write(*,*) 'Output file could not be opened - try again'
       write(*,*) oname
 902   continue
       close(20)
       close(21)
       close(30)
       end
