program BinExtract ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c program to extract subregions in the GRAVSOFT BINARY file c single line format having lat,lon,field values on each line c This version has the option of specifying a small region c and only data within this region will be output c se text below - OA 04-2008 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) fname = 'DTU13GRA_1min.bin' if (ifil.eq.2) fname = 'DTU13MSS_1min.bin' if (ifil.eq.3) fname = 'DTU13BAT_1min.bin' if (ifil.eq.4) fname = 'DTU13MDT_1min.bin' if (ifil.ge.5) fname = 'DTU13ERR_1min.bin' if (jfil.le.1) ename = 'DTU13GRA_1min.bin' if (jfil.eq.2) ename = 'DTU13MSS_1min.bin' if (jfil.eq.3) ename = 'DTU13BAT_1min.bin' if (jfil.eq.4) ename = 'DTU13MDT_1min.bin' if (jfil.ge.5) ename = 'DTU13ERR_1min.bin' open(20,file=fname,form='unformatted',status='unknown', .access='direct',recl=64,err=900) open(21,file=ename,form='unformatted',status='unknown', .access='direct',recl=64,err=900) read(21,rec=1) icode,rla1,rla2,rlo1,rlo2,dla,dlo,lutm,iell,izone read(20,rec=1) icode,rla1,rla2,rlo1,rlo2,dla,dlo,lutm,iell,izone if (icode.ne.777) stop 'binary check code 777 missing' 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 write(*,*)' Your area will most likely not match the extracted' write(*,*)' region perfectly due to the 1 minute grid spacing' 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) do i = la2,la1,-1 call inrow(20,grow,i,nn,ne) call inrow(21,erow,i,nn,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 subroutine inrow(iunit,grow,i,nn,ne) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c i n r o w c c reads one row of the grid in binary direct access format. the row c records are starting from the north. 'iunit' must be opened with c a record length equal to 16 reals c c (c) rf dec 89 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real*4 grow(*) jrec = (ne-1)/16 irec = (nn-i)*(jrec+1)+2 do 10 j = 0, jrec 10 read(iunit,rec=(irec+j)) (grow(j*16+k),k=1,16) return end c