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 5 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(*,*)' ' open(20,file='DTU13MDT_5min.mdt') open(21,file='DTU13ERR_5min.err') read(21,*) rla1,rla2,rlo1,rlo2,dla,dlo read(20,*) rla1,rla2,rlo1,rlo2,dla,dlo if (dlat.gt.0.017) dlat = 5.0d0/60.0d0 if (dlon.gt.0.017) dlon = 5.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