c
c                        program greduc
c
c  basic gravity reduction sans terrain calculation
c
c  Mike Webring
c  U.S. Geological Survey
c  po box 25046, mail stop 964
c  Denver Federal Center 
c  Denver, CO  80225 
c
c  The proper location for the data file GMETER.DAT needs to be
c  specified in the 'sfile = ' statement.
c
c  The third title line now contains the following information:
c  key, meter, project, elv units, grv units, sort key,
c  inner terr range, outer terr range
c
c  key is an integer representing the major steps of the reduction
c   0 null
c   1 meter reading converted to mgal units
c   2 drift corrected and tied to the regional net
c
c  The next seven specifiers are each 8 characters in length
c   meter - if the gravity field contains meter readings that need to be
c           converted to mgal units using some constant, then this is
c           matched with meter ids in a data file to obtain those constants.
c   project - any 8 character project identifier
c   elv units - either 'elev=f' , with a precision of .1 foot
c                   or 'elev=m' , with a precision of .01 metre
c   grv units - either 'gu=.01' , where 7 digits are reserved with
c                                 2 places to the right of the decimal 
c                   or 'gu=.001' , where 9 digits are stored with 3 places
c   sort key - 'srt=???' with any of: id,tme,lat,lon,elv,grv,itc,otc
c                and is also the suffix of the output filename.
c   inner terr - specifies the Hammer zone range for the inner correction
c                usually done by hand or high resolution methods.
c                'itc=???'    example 'itc=d-e' for zones d through e
c   outer terr - the range for the outer zones normally done with a digital
c                terrain data base.  example 'otc=f-x' for zones f through
c                166.7 kilometres.  The Hammer zones only go up to 'm'
c                (21.9 km) so 167 is arbitarily set to 'x'.  This leaves 'y'
c                and 'z' for custom calculations beyond the standard distance.
c                  
c  >>data format<<
c
c  example record
c  sta01   833651200  40 435999-1055999 52800co8000000  1599  1599-20000 -123
c  a       b c  d   e   f  g   h   i   j     k l      m     n     o     p    q
c
c  1234567890123456789012345678901234567890123456789012345678901234567890123456
c          10        20        30        40        50        60        70      
c
c  explanation
c  a	8 character station id
c  b	year
c  c	julian day
c  d	time, 24 hour clock
c  e	time zone, one decimal place
c  f	latitude degrees
c  g	latitude minutes, 2 decimals
c  h	longitude degrees
c  i	longitude minutes, 2 decimals
c  j	elevation, 1 decimal for feet, 2 for metres
c  k	state abbreviation, overwritten in a microgal survey
c  l	gravity (minus 9*10^5 mgal) or meter reading (meter units), 2 decimals,
c  	3 for microgal surveys
c  m	inner terrain correction, mgal units, 2 decimals
c  n	outer terrain correction, mgal units, 2 decimals
c  o	special field value like isostatic or regional anomalies,
c  	mgal units, 2 decimals
c  p	tide (deleted after the drift correction),
c  	units of microgal, no decimal
c  q	columns 75-80 unassigned
c       A 'd' in column 80 is a delete flag that program ganom recognizes.
c
c
c  fortran format when elevation is in feet and standard precision for gravity:
c  (a8, i2,i3,2i2,f4.1,  i3,f4.2,i4,f4.2, f6.1, a2, f7.2, 3f6.2, f5.0, a6)
c   id      time            lat/lon       elev  st  grav  corr   tide  codes
c
c  and for metres and/or microgal surveys       
c  (a8, i2,i3,2i2,f4.1,  i3,f4.2,i4,f4.2, f6.2,     f9.3, 3f6.2, f5.0, a6)
c                                            ^       ^ ^
c                                            |       | |
c
	common /gfiles/ icont, itty, jtty, f1,f2,f3,f4,f5
	character*56                       f1,f2,f3,f4,f5
	character       ans1*4, ans16*16, prog*8
	equivalence     (ans1,ans16)
	itty = 5
	jtty = 6
	write(jtty,5)
5	format(/' basic gravity processing system '/)
	prog='greduc'
c
	ifunc = 1
	icont = 0
10	write(jtty,15)
15	format(' function :'$)
	read(itty,20) ans16
20	format(a8)
	call cvc(ans1,2)
c
 	if(ans1.eq.'cont' .or. ans1.eq.'star') icont=1
	if(ans1.eq.'cont' .and. ifunc.eq.1) then
	  write(jtty,*,err=10) ' begin at'
	  go to 10
	endif
	if(ans1.eq.'cont' .and. ifunc.gt.1) then
	  ifunc=ifunc+1
	  go to (30,100,200,300,400,500,900,1000) ifunc
	endif	
	ifunc=1
	if(ans1.eq.'star') ifunc= 2
	if(ans1.eq.'conv') ifunc= 2
	if(ans1.eq.'tide') ifunc= 3
	if(ans1.eq.'mete') ifunc= 4
	if(ans1.eq.'sort') ifunc= 5
	if(ans1.eq.'drif') ifunc= 6
	if(ans1.eq.'anom') ifunc= 7
	if(ans1.eq.'exit') go to 1000
	if(ans1.eq.'stop') go to 1000
	go to (30,100,200,300,400,500,900) ifunc
30	write(jtty,40)
40	format(' functions :start, continue, stop',/,
     1 ' convert, tide, meter, sort, drift, anomaly')
	go to 10
c
100	call gconv(key)
	if(key.eq.1) go to 400
	if(key.eq.2) go to 510
	if(icont.eq.0) go to 10
200	call gtide
	if(icont.eq.0) go to 10
300	call gmeter
	if(icont.eq.0) go to 10
400	call gsort
	if(icont.eq.0) go to 10
500	call gdrift
	if(icont.eq.0) go to 10
510	write(jtty,550)
550	format(/' do you want an anomaly listing ?')
        if(noyes().eq.0) stop
900	call ganom
1000	stop
	end
c******************************************************************************
	subroutine gconv(key)

c  simple checks, reorder 3rd header, do elevation unit change.

	common/gfiles/ icont,itty,jtty,f1,f2,f3,f4,f5
	dimension elcon(2)
c
	character*56 f1,f2,f3,f4,f5
	character*8 meter,projct,feet,gunit,sort,terr1,terr2
        character suff*3,id*8,state*2
	character rec*80, title*80, hdr79*79, title3*56
c
	equivalence (title3(1:8),meter),(title3(9:16),projct),
     1 (title3(17:24),feet),(title3(25:32),gunit),(title3(33:40),sort),
     1 (title3(41:48),terr1),(title3(49:56),terr2)
c
	data suff/'cnv'/, xm2f/3.28084/
c
c  elevation conversion = 0 no change, = 1 is ft to mt, = 2 mt to ft.
	ielcon = 0
	f2m    = 1.0 / xm2f
	icnt   = 0
c1       format(a80)
1       format(a)
c
	write(jtty,100)
100	format(' Enter your input filename: '$)
	read(itty,110) f1
110	format(a56)

c  blank = 'zero' causes fortran to read values and place decimals strictly
c  by the format specification.  default is to pick up a value any place
c  in the field and assign decimal place in the nonzero string that it finds.
c  ie. if a user leaves off a trailing zero they've divided by ten...

	open (10, file=f1, form='formatted', status='old',
     1     blank='zero')
	call getfn(f1,suff,f3)
	open (11, file=f3, form='formatted', status='unknown')
c
	do 121 i = 1, 2
	  read(10, 1, end=380) title
	  write(11, 1)         title
121	continue
c
	read(10,160) key,hdr79
160	format(i1,a79)
	call gheadr(hdr79,title3)
	title = ' '
	if (key.eq.1)
     1   title=' key=1, relative gravity without drift correction'
	if (key.eq.2) 
     1   title=' key=2, observed gravity file'
	write(jtty,1) title
c
	write(*,*)' are you inputting west longitude values ?'
        ny = noyes()
	if(ny.eq.1) isign1 = -1.
	if(ny.eq.0) isign1 =  1.
c
200	if (feet.ne.'elev=f ' .and. feet.ne.'elev=m ') then
	  write(*,*)' are your input elevations in feet or metres :'
	  read(*, 1) feet
	  call cvc(feet,2)
	  if (feet.eq.'feet' .or. feet.eq.'metres') then
	    if (feet .eq. 'feet')   feet = 'elev=f'
	    if (feet .eq. 'metres') feet = 'elev=m'
	  else
		    write(*,*)'>>',feet, '...answer must be ''feet'' or ''metres'' '
		    go to 200
	  endif
	endif
c
	write(*,*) 'do you want to change elevation units ?'
        if (noyes() .eq. 1) then
	  if (feet(1:6) .eq. 'elev=f' ) then
c           increase decimals to hundredths and store conversion
	    ielcon   = 1
	    elcon(1) = 10.0 * f2m
	    feet     = 'elev=m'
	    write(*,*)' converting elevations to metres'
	    else if (feet(1:6) .eq. 'elev=m' ) then
c             decrease decimals to tenths and store conversion
	      ielcon   = 2
	      elcon(2) = 0.1 * xm2f
	      feet     = 'elev=f'
	      write(*,*)' converting elevations to feet'
	    else
	    write(*,*)' cannot decode input units'
	    go to 200
	  endif
	endif	
c
	if (gunit.ne.'gu=.01' .and. gunit.ne.'gu=.001') then
	  write(*,*)' is this a microgal survey ? '
	  gunit = 'gu=.01'
          if (noyes() .eq. 1) gunit = 'gu=.001'
	endif
	iu = 0
	if (gunit .eq. 'gu=.001') iu = 1
c
	write(11,299) key,meter,projct,feet,gunit,sort,terr1,terr2
299	format(i1,7a8)
c
c check data records 
c
300	read(10, 1, err=888, end=400) rec
	read(rec, 310, err=311) id, iyr, iday, ihour, min, izone,
     1    latd, latm, lond, lonm, iele, state, ig
310	format( a8, i2,i3,2i2,i4, i3,3i4, i6, a2,i7 )
	go to 312
311	  write(*,*)'error parsing station record: id=', rec(1:8)
	  iec = iec + 1
	  if (iec .gt. 10) stop
	  go to 350
312	continue
c
	lond = isign( lond, isign1 )
	write( rec(29:32), 315) lond
315	format(i4)
c
	if (ielcon.ne.0) then
	  iele = int( elcon(ielcon) * float(iele) + .5)
	  write( rec(37:42), 320) iele
320	  format( i6 )
	endif
c
c  for a microgal survey clear spaces for more precision
	if (iu .eq. 1) then
	  n = ichar( state(1:1) )
	  if (n.gt.57 .or. (n.lt.48 .and. n.ne.45)) state(1:1)=' '
	  n = ichar( state(2:2) )
	  if (n.gt.57 .or. (n.lt.48 .and. n.ne.45)) state(2:2)=' '
	endif
c
350	write(11,1) rec
	icnt = icnt + 1
c
c  warning messages
	if (iday.gt.366) write(*,*)' station ',id,' has a day error'
	if (ihour.ge.24) write(*,*)' station ',id,
     1                       ' has a time (hour) error'
	if (min.ge.60)   write(*,*)' station ',id,
     1                       ' has a time (minute) error'
	if(latd.gt.90)   write(*,*)' station ',id,' has latitude error'
	if(latm.ge.6000) write(*,*)' station ',id,' has latitude error'
	if(lonm.ge.6000) write(*,*)' station ',id,' has longitude error'
	go to 300
c  end data record loop
c
380	write(390)
390	format(' premature end of titles')
	icont = 0
	go to 999
400	write(jtty,410) icnt
410	format(' grcvt finished,', i7, ' stations checked')
	f1 = f3
	go to 999
888	write(*,*)' input conversion error at record',icnt
	iec = iec + 1
	if (iec.gt.10) stop
	go to 300
999	close(10)
	close(11)
	return
	end
c******************************************************************************
	subroutine gheadr(headr3,title3)
c  Parsing routine for the third header record of 
c  the Denver standard gravity file.
c  Returned is a record with the keywords reordered and formatted.
	character headr3*79, hdr79*79, title*56, title3*56
	character*8 meter, projct, feet, gunit, sort, terr1, terr2
	character*8 str(10), test
c
	equivalence (title(1:8),meter), (title(9:16),projct),
     1 (title(17:24),feet),(title(25:32),gunit), (title(33:40),sort),
     1 (title(41:48),terr1), (title(49:56),terr2)
c
	hdr79 = headr3
	title = ' '
	call cvc(hdr79,2)
	nstr = 0
	do 20 nword = 1, 10
	  if (leftj(hdr79) .eq. 0) go to 50
	  mblk = index(hdr79,' ')
	  if (mblk.ge.2) then
	    nstr = nstr + 1
	    str(nstr) = hdr79(1:mblk-1)
	    do 10 i = 1, mblk
	      hdr79(i:i) = ' '
10	    continue
	  endif
20	continue
c
50	nword = nstr
	mtr   = 0
	iprj  = 0
	do 100 iword = 1, nword
	test = str(iword)
	meq = index(test,'=')
	if (meq.gt.0) then
	  if (test(1:meq) .eq. 'elev=') feet  = str(iword)	
	  if (test(1:meq) .eq. 'gu=')   gunit = str(iword)	
	  if (test(1:meq) .eq. 'srt=')  sort  = str(iword)	
	  if (test(1:meq) .eq. 'itc=')  terr1 = str(iword)	
	  if (test(1:meq) .eq. 'otc=')  terr2 = str(iword)	
	  else
	  if (mtr .eq. 0) then
	    meter = str(iword)
	    mtr = 1
	    else if (iprj .eq. 0) then
	      projct = str(iword)
	      iprj = 1
	  endif
	endif
100	continue
c
	title3 = title
	return
	end
c******************************************************************************
	subroutine gtide
	common /gfiles/ icont,itty,jtty,f1,f2,f3,f4,f5
	character*56 f1,f2,f3,f4,f5
	character*8 meter,projct,feet,gunit,sort,terr1,terr2
        character suff*3,title*80,idat1*69,idat2*5,idat3*6
	double precision hour,min,gmt,year,elev,latd,latm,
     1 lond,lonm,atide,elevm
	data suff/'tid'/
	ier=1
	if(icont.gt.0) go to 120
	write(jtty,100)
100	format(' Enter your input filename :'$)
	read(itty,110)f1
110	format(a56)
	go to 140
120	write(jtty,130)
130	format(/' do you want tides calculated ?'$)
        if(noyes().eq.0) return

c
140	open(10,file=f1,form='formatted',status='old',
     1 blank='zero')
	call getfn(f1,suff,f3)
	open(11,file=f3,form='formatted',status='unknown')
c
	do 200 i=1,2
	read(10,210,end=340)title
200	write(11,210)title
210	format(a80)
	read(10,220) key,meter,projct,feet,gunit,sort,terr1,terr2
220	format(i1,7a8)
	write(11,220) key,meter,projct,feet,gunit,sort,terr1,terr2
	if(key.eq.2) go to 400
	iconv=-1
	if(feet.eq.'elev=m') iconv=0
	if(feet.eq.'elev=f') iconv=1
        if(iconv.eq.-1) stop ' elevation units in header not correct'
c
230	irec=0
	inrec=0
240	read(10,310,end=290) year,iday,hour,min,gmt,latd,latm,
     1 lond,lonm,elev,idat1,idat2,idat3
310	format(8x,f2.0,i3,2f2.0,f4.1,f3.0,f4.2,f4.0,f4.2,f6.2,
     1 t1,a69,a5,a6)
	inrec=inrec+1
	if(year.le.0.d0 .or. iday.le.0) go to 280
	if(latd.eq.0.d0 .and. latm.eq.0.d0 .and.
     1 lond.eq.0.d0 .and. lonm.eq.0.d0) go to 280
	if(latd.gt.90.d0 .or. latd.lt.-90.d0) go to 280
	irec=irec+1
	lond=-lond
	lonm=dsign(lonm,lond)
	elevm=elev
	if(iconv.eq.1) elevm=elev*.3048
	zone=aint((lond+lonm/60.)/15. + .5)
	if(zone.gt.gmt+2.d0 .or. zone.lt.gmt-2.d0) write(jtty,260) inrec
260	format(' record number ',i4,' may have an incorrect time zone')
	call tides(hour,min,gmt,iday,1,year,elevm,latd,latm,
     1 lond,lonm,atide)
	itide=int(atide*1.d3+dsign(.5d0,atide))
270	write(11,320) idat1,itide,idat3
320	format(a69,i5,a6)
	go to 240
280	idat2=' '
	write(11,330)idat1,idat2,idat3
330	format(a69,a5,a6)
	go to 240
c
290	ier=0
	f1=f3
	write(jtty,300)inrec,irec
300	format(1x,i6,' Input',' stations,',i6,' Tides calculated')
	if(inrec.ne.irec) write(jtty,360)
360	format(' warning: Time and location fields (lat or lon nonzero)',/,
     1 ' must be valid for tide to be computed.',/,
     1 ' The loops without tide values should be short enough that',/,
     1 ' a linear trend is a good approximation of drift.')
	go to 999
340	write(jtty,350)
350	format(' premature end of titles')
	icont=0
	go to 999
400	write(jtty,410)
410	format(' key indicates this is an observed gravity file',/,
     1 ' no tides calculated')
	icont=0
999	close(10)
	close(11)
	return
	end
c******************************************************************************
      subroutine tides(hrsr,xmin,gmt,iday,month,yrs,elev,xlat,xlatm,
     &  xlon,xlonm,grav)
c computation routine for tides revised by r wahl, 12 sep 74
c  ref geophysics monograph no 16, r.a. brouche  (?)
c  positive west longitude, elev in meters, grav in mgals
      implicit real*8         (a-h,o-z)
	common /smtide/gs,gm
      dimension mthday(12)
      data mthday/0,31,59,90,120,151,181,212,243,273,304,334/
      idate=iday+mthday(month)
      iyears=yrs
      if (mod(iyears,4).eq.0.and.month.gt.2)  idate=idate+1
      daysr=dfloat(idate)
      daysq=yrs*365.d0+daysr-.5d0+dfloat(idint((yrs-1.d0)*.25d0+
     & 1.d-6))
      hrsd=xmin*.01666666666667d0+gmt+hrsr
      daysd=daysq+.041666666666667d0*hrsd
      ht=elev*1.d-9
      t=daysd*2.737850787d-5
      ta=dfloat(idint(t*100.d0+1.d-6))*1.d-2
      tc=(t-ta)*100.d0
      tb=ta*1336.d0
      tb=360.d0*(tb-dfloat(idint(tb)))
      al=(xlat+xlatm*.0166666666666667d0)*.0174532925199433d0
      cosal=dcos(al)
      sinal=dsin(al)
      cos2al=cosal*cosal
      r=ht+6.37816d-3*(.99664712d0+cos2al*(3.381028d-3
     1 -cos2al*2.8152d-5))
      a=13.36d0*tc
      b=dfloat(idint(a))
      c=(a-b)*360.d0+tb
  100 s=270.43435833d0+c+t*(307.883141667d0+t*(-1.13333d-3+t*
     1 1.888888d-6))
      s=dmod(s,360.d0)*.0174532925199433d0
      a=11.d0*t
      b=dfloat(idint(a))
      c=(a-b)*360.d0
  110 p=334.32965777d0+c+t*(109.03403d0+t*(-1.0325d-2-t*1.25d-5))
      p=dmod(p,360.d0)*.0174532925199433d0
      c=360.d0*tc
  120 h=279.69668d0+c+t*(.768925d0+3.03d-4*t)
      h=dmod(h,360.d0)*.0174532925199433d0
      a=5.d0*t
      b=dfloat(idint(a))
      c=(a-b)*360.d0
  130 oln=259.183275d0-c+t*(-134.14201d0+t*(2.07777d-3+t*2.d-6))
      oln=dmod(oln,360.d0)
      if (oln.lt.0.d0)  oln=oln+360.d0
      oln=oln*.0174532925199433d0
      sinoln=dsin(oln)
      cosoln=dcos(oln)
      ci=.91369d0-0.03569d0*cosoln
      si=dsqrt(1.d0-ci*ci)
      sinsi=sinoln/si
      sn=.08968d0*sinsi
      cn=dsqrt(1.d0-sn*sn)
      sit=.39798d0*sinsi
      cit=cosoln*cn+.91739d0*sinoln*sn
      tit=sit/(1.d0+cit)
      et=2.d0*datan(tit)
      if (et.lt.0.d0)  et=et+6.28318530717958d0
      cxi=oln-et
      sig=s-cxi
      smp=s-p
      sinsmp=dsin(smp)
      cossmp=dcos(smp)
      sm2hp=s-h-h+p
      smhx2=2.d0*(s-h)
      olm=sig+sinsmp*(.1098d0+7.535024d-3*cossmp)+.0154003d0*
     1 dsin(sm2hp)+9.75199d-3*dsin(smhx2)
      v=datan2(sn,cn)
      ha=(15.d0*(hrsd-12.d0)-xlon-xlonm*.016666666666667d0)*
     1 .0174532925199433d0
      chis=ha+h
      chi=chis-v
      sinolm=dsin(olm)
      ct=sinal*si*sinolm+0.5d0*cosal*((1.d0+ci)*dcos(olm-chi)
     1 +(1.d0-ci)*dcos(olm+chi))
      da=2.59356456d0+cossmp*(.14325025d0+.01572888d0*cossmp)
     1 +.0200919d0*dcos(sm2hp)+.01859186d0*dcos(smhx2)
      ct2=ct*ct
      gm=r*da*da*da*.49049d0*(3.d0*ct2-1.d0+1.5d0*r*da*ct*
     1 (5.d0*ct2-3.d0))
      es=.01675104d0-t*(4.18d-5+t*1.26d-7)
      ps=281.22083d0+t*(1.71902d0+t*(4.53d-4+t*3.d-6))
      ps=dmod(ps,360.d0)*.0174532925199433d0
      hmps=h-ps
      ols=h+2.d0*es*dsin(hmps)
      ds=6.684559d-3*(1.d0+es*dcos(hmps)/(1.d0-es*es))
      cf=.39798d0*sinal*dsin(ols)+cosal*(.95869d0*
     1 dcos(ols-chis)+.0413d0*dcos(ols+chis))
      gs=13.2916d6*r*(3.d0*cf*cf-1.d0)*ds*ds*ds
      grav=1.2d0*(gm+gs)
  140 return
      end
c******************************************************************************
	subroutine gmeter
c  conversion from meter reading to mGal units
c  Mike Webring, USGS, VAX version 6/83
c
c      Format of the meter constant file where w-351 is a Worden meter
c      and g-159 is a LaCoste meter.
c  meter w-351  1       
c  0.0949     
c  meter g-159  2                                              
c  * 1.00028                  supplemental                     
c  0000,0,1.05333                                              
c  0100,105.33,1.05352                                         
c  0200,210.69,1.05351                                         
c  0300,316.04,1.05338                                         
c  ... and continuing for the rest of the LaCoste interval constants.
c
	common /gfiles/ icont,itty,jtty,f1,f2,f3,f4,f5
	dimension basm(100),factr(100)
	character*56 f1,f2,f3,f4,f5,sfile
	character*8 gmtr,projct,feet,gunit,sort,terr1,terr2
        character*6 meter,mname,name
	character suff*3,mtype*2,flag1*1,cfudge*12,
     1 idat*24,title*80,card1*44,card2*29
	logical normg
	data suff/'mtr'/,maxint/100/
	itty=5
	jtty=6
c
	write(jtty,99)
99	format(/,'   meter reading conversion to mgal')
	if(icont.gt.0) go to 140
	write(jtty,100)
100	format(' enter your input filename  :'$)
	read(itty,110)f1
110	format(a56)
140	write(jtty,130)
130	format(' enter your supplemental meter multiplier :'$)
	read(itty,*) fudge
	if (fudge.le.0.0) fudge = 1.0
c
	open(9,file=f1,form='formatted',status='old',
     1 blank='zero')
c
c  --  !! Change following directory to location of GMETER.DAT !!
c
        sfile = '\PF\BIN\gmeter.dat'
	open(10,file=sfile,form='formatted',status='old')
	call getfn(f1,suff,f3)
	open(11,file=f3,form='formatted',status='unknown')
c
c  read/write titles and check key
	do 160 i=1,2
	read(9,150,end=220) title
150	format(a80)
	write(11,150) title
160	continue
	read(9,170,end=220,err=540)
     1 key,gmtr,projct,feet,gunit,sort,terr1,terr2
	na=leftj(gmtr)
	mname(1:6)=gmtr(1:6)
	call cvc(mname,2)
c  character comparisons are all done in lowercase
	if(key.eq.2) go to 180
	if(key.eq.1) go to 200
	key=1
	write(11,170) key,gmtr,projct,feet,gunit,sort,terr1,terr2
170	format(i1,7a8)
	n=leftj(gunit)
	call cvc(gunit,2)
	normg=.true.
	if(gunit.eq.'gu=.001') normg=.false.
	go to 240
180	write(jtty,190)
190	format(' cannot process - already observed gravity.')
	go to 540
200	write(jtty,210)
210	format(' file already converted')
	go to 540
220	write(jtty,230)
230	format(' Premature end of titles',/,'please',
     1 ' check your data for the proper format.')
	go to 540
c
c  default meter scaling
240	if (mname.ne.'one   ') go to 260
	scalm=fudge
	itype=1
	go to 450
c
c  search for 'meter' and meter id
260	read(10,270,end=300) meter,name,mtype
270	format(2a6,a2)
	if(meter.ne.'meter ') go to 260
	if(name.ne.mname) go to 260
	read(mtype,280) itype
280	format(i2)
	if (itype.eq.1)  go to 320
	if (itype.eq.2)  go to 350
	write(jtty,290)itype,name
290	format(' erroneous type number (',i2,') for meter: ',a6)
	go to 540
300	write(jtty,310)mname
310	format(' meter ',a6,' not in calibration file')
	go to 540
c
c  worden meter
c
320	read(10,330,end=410) scalm
330	format(f7.4)
	scalm = scalm * fudge
	write(*, 333)name,scalm
333	format(' scaling factor for worden meter ',a6,' is ',1pg13.5)
	go to 450
c
c  la coste meter
c
c  read supplemental meter multiplier associated with lacoste meter.
350	read(10,354,end=415) flag1,cfudge
354	format(a1,a)
	if (flag1.ne.'*') stop ' missing supplemental multiplier'
	read(cfudge,355) fudge1
355	format(f10.6)
	if (fudge1.le.1.e-4) fudge1 = 1.0
	write(jtty,356) name,fudge1
356	format(' meter ',a6,' multiplier from table file = ',f9.6)
	fudge = fudge * fudge1
	write(jtty,357) fudge
357	format(' total meter multiplier for data = ',f9.6)
c
c  get interval constants for lacoste
	imtr=0
	kcc=0
	v=0.0
360	read(10,370,end=410) idat
370	format(a24)
	if(idat(1:5).eq.'meter' .or. imtr.gt.maxint) go to 410
	imtr=imtr+1
	call lacoste(idat,icc,basm(imtr),factr(imtr))
	if (icc.ne.kcc) then
	 write(jtty,400)name
400	 format(' sequence error for meter:',a6)
	 go to 540
	endif
	kcc=kcc+100
	if(abs(v-basm(imtr)).gt..02) write(*, 402)
	v=v+factr(imtr)*100.
402	format(' possible error in meter contant file')
	go to 360
c
410	if (imtr.eq.0) then
415	 write(jtty,420)name
420	 format(' no interval constants for meter ',a6)
	 go to 540
	endif
	write(*, 600)name
600	format(' do you want to see the interval constants for meter ',
     1 a6,' ?'$)
        if(noyes().eq.1) then
	icc=0
	v=0.0
	write(*, 601)
601	format(' reading  mgal   factor     check')
	do 610 i=1,imtr
	write(*, 605)icc,basm(i),factr(i),v
605	format(1x,i7,f8.2,f8.5,3x,f8.2)
	v=v+100.*factr(i)
610	icc=icc+100
	endif
c
c  process data file
450	if(normg) then
	read(9,460,end=550) card1,value,card2
460	format(a44,f7.2,a29)
	else
	read(9,461,end=550) card1,value,card2
461	format(a42,f9.3,a29)
	endif
	go to (470,480),itype
470	value=value*scalm
	go to 510
480	k=int(value*.01)+1
	if (k.gt.imtr) then
	ivalue=int(value*100.+.5)
	 write(*,  490)card1,ivalue,card2
490	 format(' meter value out of range.'/
     1 1x,a44,i7,a29)
	 go to 450
	endif
	value=(basm(k)+factr(k)*amod(value,100.))*fudge
c
510	if(normg) then
c  10,000.00 is a marker that indicates untied gravity readings.
	ivalue=int(sngl(dble(value)*100.0d0+.5d0))+1000000
	write(11,520) card1,ivalue,card2
520	format(a44,i7,a29)
	else
	ivalue=int(sngl(dble(value)*1000.0d0+.5d0))+10000000
	write(11,521) card1,ivalue,card2
521	format(a42,i9,a29)
	endif
	go to 450
c
540	icont=0
	go to 999
550	f1=f3
999	close(9)
	close(10)
	close(11)
	return
	end
c******************************************************************************
	subroutine lacoste(inrec,ia,b,c)
c  inrec has no embedded blanks
	character tmp*24,ca*8,cb*8
	character inrec*(*)
	m = index(inrec,' ') - 1
	if(m.le.10) go to 4
	ca = inrec(1:4)
	read(ca,1) ia
1	format(i4)
	tmp=inrec(6:m)
	n = index(tmp,',') - 1
	cb=tmp(1:n)
	read(cb,2) b
2	format(f8.2)
	n=n+2
	tmp=tmp(n:m)
	read(tmp,3) c
3	format(f8.5)
	return
4	tmp=inrec
	write(*, 5)inrec
5	format(' record ',a24,' is too short')
	return
	end
c*****************************************************************************
c******************************************************************************
	subroutine gdrift
c  meter drift correction routine
c  A linear drift is assumed between any two bases 
c  obtained from a base station file. 
c  The input data must be time sorted.
c  M Webring, 6/83
	common /gfiles/ icont,itty,jtty,f1,bases,f3,gtime,ofile
	double precision basg(200),rg,sread,sbase,shift,shift2,dgrav,og
	character*56 f1,bases,f3,gtime,ofile
	character*8 sid,basid(200),gmtr,projct,feet,gunit,sort,
     1           terr1,terr2
	character suff1*3,title*80,irec1*27,irec2*18,irec3*6,
     1         irec4*6,bunit*8
	character chog*15
	logical first,normg,normb
	data suff1/'og'/,irec3/' '/,izero/0/
	write(*, 1)
1	format(/,'   drift correction routine',/,
     1 ' do you want to see the drift constants in ''drift.tmp'' ?'$)
	title='scratch'
        if(noyes().eq.1) title='new'
	open(8,file='drift.tmp',form='formatted',status=title)
	write(8,2)
2	format(' rec #       date         base        shift',
     1     '          mgal/hr')
	if(icont.eq.0) then
	  write(*, 3)
3	  format(' enter your time sorted data filename :'$)
	  read(itty,37) gtime
37	  format(a56)
	endif
	open(10,file=gtime,form='formatted',status='old',
     1 blank='zero')
	write(*, 4)
4	format(' enter your base station filename :'$)
	read(itty,37) bases

	open(9,file=bases,form='formatted',status='old',
     1 blank='zero')
c
	call getfn(gtime,suff1,ofile)
	open(11,file=ofile,form='formatted',status='unknown')
	read(10,5) title
	read(10,5) title
5	format(a80)
	read(10,6) key,gmtr,projct,feet,gunit,sort,terr1,terr2
6	format(i1,7a8)
	if(key.ne.1) then
	  write(jtty,7)
7	  format(' input data has not been converted to mgals')
	  return
	endif
	if(key.eq.2) then
	  write(*, 8)
8	  format(' input data is already drift corrected')
	  return
	endif
	ju = -1
	n  = leftj(gunit)
	call cvc(gunit,2)
	if( gunit .eq. 'gu=.01'  ) ju = 0
	if( gunit .eq. 'gu=.001' ) ju = 1
	if( ju .eq. -1 ) then
10	  write(*, 11)
11	  format(' enter gravity data precision, .01 or .001 mgal :'$)
	  read(itty,*,err=10) prec
	  ju = 0
	  if( prec .lt. 0.005 ) ju = 1
	endif
	if( ju .eq. 0 ) normg = .true.
	if( ju .eq. 1 ) normg = .false.
c
c  read in base station data
c
	do 15 i = 1, 2
	  read(9,5)
15	continue
	read(9,16) bunit 
16	format(25x,a8)
	normb = .true.
	if( bunit .eq. 'gu=.001' ) then
	  normb = .false.
	  write(*,*)' base station file being read with .001 precision'
	  else
	  write(*,*)' base station file being read with .01 precision'
	endif	
	nbase = 0
20	nbase = nbase + 1
	if( normb ) then
	  read(9,30,end=40) basid(nbase), basg(nbase)
30	  format(a8,36x,f7.2)
	  else
	  read(9,31,end=40) basid(nbase), basg(nbase)
31	  format(a8,34x,f9.3)
	endif
	n = leftj( basid(nbase) )
	call cvc(  basid(nbase), 2 )
	go to 20
c
c  meter drift for each base loop
c
40	nbase = nbase - 1
	first = .true.
	drift = 0.0
	ista  = 0
100	ista  = ista + 1
	if(normg) then
	  read(10,110,end=200) sid,iyr,iday,ihr,min,irg,tide
110	  format(a8,i2,i3,2i2,27x,i7,18x,f5.3)
	  rg = dfloat(irg) * 1.0d-2
	else
	  read(10,111,end=200) sid,iyr,iday,ihr,min,irg,tide
111	  format(a8,i2,i3,2i2,25x,i9,18x,f5.3)
	  rg = dfloat(irg) * 1.0d-3
	endif
	hour = float(ihr) + float(min) / 60.
	n    = leftj( sid )
	call cvc( sid, 2 )
c
	do 150 ibase=1,nbase
	if(sid.ne.basid(ibase)) go to 150
	shift = basg(ibase)-(rg+dble(tide))
c
	if(.not.first) then
c  calculate drift for preceeding interval
	dread = sngl( (rg+dble(tide))-sread )
	dbase = sngl( basg(ibase)-sbase )
	drift=0.0
	nday=0
	if(iyr.gt.isyr) then
	 nday=365
	 if(mod(isyr,4).eq.0) nday=366
	endif
	dtime=float((iday+nday)-isday)*24.+(hour-shour)
	if(dtime.lt.0.0) write(*,*)' incorrect base loop at record',ista
	if(dtime.eq.0) dtime=1.e-10
	drift=(dread-dbase)/dtime
	endif
c
	if(first .and. ista.ne.1) write(8,120) izero,isyr,isday,shour,
     1 basg(ibase),shift,drift
	write(8,120) ista,iyr,iday,hour,basg(ibase),shift,drift
120	format(1x,3i4,f7.3,1p2d15.7,1pe15.5)
	first=.false.
	lastb=ista
	isyr=iyr
	isday=iday
	shour=hour
	sread=rg+dble(tide)
	sbase=basg(ibase)
	go to 100
150	continue
c  test for closed loop at beginning of data
	if(ista.ne.1) go to 100
	isyr=iyr
	isday=iday
	shour=hour
	go to 100
c
c  drift correction
c
200	nsta=ista-1
	if(first) stop ' no base stations matched with data'
	drift=0.0
c  test for closed loop at the end of data
	if(lastb.ne.nsta) write(8,120)
     1 nsta,iyr,iday,hour,sbase,shift,drift
	write(8,201)
201	format(/,' base and shift are for the station listed',/,
     1 ' the drift is for the preceeding interval',/,
     1 ' shift = base-relgrv   where relgrv = meter+tide',/,
     1 ' base loop drift rate = (drelgrv-dbase)/dtime of base loop',/,
     1 ' and obs gravity = shift + relgrv-drift*dtime from begin loop')
	rewind 8
	read(8,120)
	read(8,120) idumm,isyr,isday,shour,basgrv,shift,zero
	read(8,120) istop,isyr2,isday2,shour2,basgrv,shift2,drift
	rewind 10
	do 205 i=1,2
	read(10,5) title
	write(11,5) title
205	continue
	read(10,6) key,gmtr,projct,feet,gunit,sort,terr1,terr2
	key=2
	write(11,6) key,gmtr,projct,feet,gunit,sort,terr1,terr2
c
	do 290 ista=1,nsta
	if(normg) then
	  read(10,210) sid,iyr,iday,ihr,min,irec1,irg,irec2,tide,irec4
210	  format(a8,i2,i3,2i2,a27,i7,a18,f5.3,a6)
	  rg=dfloat(irg)*1.d-2
	  else
	  read(10,211) sid,iyr,iday,ihr,min,irec1,irg,irec2,tide,irec4
211	  format(a8,i2,i3,2i2,a25,i9,a18,f5.3,a6)
	  rg=dfloat(irg)*1.0d-3
	endif
	hour=float(ihr)+float(min)/60.
	nday=0
	if(iyr.gt.isyr) then
	 nday=365
	 if(mod(isyr,4).eq.0) nday=366
	endif
	dtime=float((iday+nday)-isday)*24.+(hour-shour)
	if(dtime.lt.-0.01) stop ' stations not time sorted'
	dgrav=rg+dble(tide-drift*dtime)
	if(normg) then
	  iog=int(sngl(100.d0*(shift+dgrav)+.5d0))
	  write(11,220) sid,iyr,iday,ihr,min,irec1,iog,irec2,irec3,irec4
220	  format(a8,i2,i3,2i2,a27,i7,a18,a5,a6)
	else
c  coded to maintain 8 place output.
	  og=1.0d3*(shift+dgrav)
	  write(chog,222) og
222	  format(d15.8)
	  n=index(chog,'.')+1
	  n2=n+7
	  chog(2:9)=chog(n:n2)
	  chog(1:1)=' '
	  write(11,221) sid,iyr,iday,ihr,min,irec1,chog(1:9),irec2,irec3,irec4
221	  format(a8,i2,i3,2i2,a25,a9,a18,a5,a6)
	endif
	if(ista.eq.istop) then 
	 isyr=isyr2
	 isday=isday2
	 shour=shour2
	 shift=shift2
	 if(istop.lt.nsta) read(8,120) istop,isyr2,isday2,shour2,
     1 basgrv,shift2,drift
	endif
290	continue
	close(8)
	close(9)
	close(10)
	close(11)
	return
	end
c******************************************************************************
	subroutine ganom
c  Bouguer anomaly calculation, standard ascii file input and choice
c  of ascii listing file or binary 'posting' file on output.
c
c  Input/output elevation units: ascii files can be either feet or metres
c  with a unit specifier in the header, binary file elevations units are the 
c  same as the input file.
c
c  Observed gravity in the ascii standard is less 900,000 mgal (maintains
c  readability), but to ensure accuracy on 32 bit machines 980,000 mgal
c  is removed for binary output.
c
c  Terrain correction density is assumed to be 2.67 gm/cc.
c
c  A 'd' in column 80 will cause that record to be skipped.


	common /gfiles/  icont,itty,jtty,f1,f2,f3,f4,f5
	character*56                     f1,f2,f3,f4,f5

	dimension        xl(4)
	double precision ogd
	logical          bin, form, normg

	character   rec*80, title*80, suff1*3, suff2*3
	character*8 id,meter,projct,feet,gunit,sort,terr1,terr2
	character*5 junit(2)
	character   ofmt*80, ofmte*4, ofmtg*5

	data suff1/'lis'/,suff2/'pst'/, junit/'metre','feet '/

	in     = 0
	iout   = 0
	idelet = 0

	if( icont .eq. 0 ) then
	  write(*,  5)
5	  format(/' enter your standard gravity filename')
	  read(5,37) f5
	  f1=f5
37	  format(a56)
	endif
	open (10,file=f5,status='old',form='formatted',
     1 blank='zero')
	write(*,  10)
10	format(' output listing=1, post=2, both=3 :'$)
	read(5,*) ifunc
	bin  = .false.
	form = .false.
	if ( ifunc .eq. 1  .or.  ifunc .ge. 3 ) form = .true.
	if ( ifunc .ge. 2 ) bin = .true.
c
	call getfn( f1, suff1, f2 )
	call getfn( f1, suff2, f3 )

c  'list' is VAX convention to suppress fortran page control characters.

	if ( form ) open (11, file=f2, status='unknown', form='formatted')
c  optional - fixed length speeds i/o.
	if ( ifunc .ge. 2 ) open(12, file=f3, status='unknown',
     1  form='unformatted', recl=10)

15	write(jtty,20)
20	format(' enter Bouguer reduction density :'$)
	read(itty,*,err=15) den

	read(10,25) title
25	format(a80)
	if( form ) write(11,25) title
	read(10,25) title
	if( form ) write(11,25) title

	read(10,30) key,meter,projct,feet,gunit,sort,terr1,terr2
30	format(i1,7a8)
	if( form ) write(11,31) meter,projct,sort,terr1,terr2
31	format( 1x, 'meter=',a8,'  project=',a8,2x,a8,/,
     1 1x,'terrain correction zones ',2a8,'.')

	n     = leftj(gunit)
	normg = .true.
	ofmtg = 'f11.2'
	if( gunit .eq. 'gu=.001' ) then
	  write(*,*)'data header indicates gravity precision is .001 mgal' 
	  normg = .false.
	  ofmtg = 'f11.3'
	endif
c
	n  = leftj(feet)
	iu = -1
	if( feet .eq. 'elev=f  ' ) iu = 1
	if( feet .eq. 'elev=m  ' ) iu = 0
	if( iu .eq. -1 ) then
	  write(*,*)' data header does not contain a elevation specifier'
	  write(*, 35)
35	  format(' enter 0 for input elevations in metres, 1 for feet :'$)
	  read(itty,*) iu
	endif
c
	if( form ) then
40	  write(*, 41)
41	  format(' for the listing, enter 0 for metre or 1 for feet :'$)
	  read(itty,*,err=40) ju
	  if( .not. ( ju.eq.0 .or. ju.eq.1) ) go to 40
	  ju1   = ju + 1
	  ofmte = 'f8.2'
	  if( ju .eq. 1 ) ofmte = 'f8.1'
	endif
c
c  encode output format
	write(ofmt,42) ofmte,ofmtg 
42	format( '(a8,i3,f6.2,i5,f6.2,', a4,
     1 ',', a5,',2(1x,f5.2), 1x,f7.2,2x,f7.2)' )
	if(form) write(11,25)
	if(form) write(11,45)
45	format('station  latitude  longitude  elev.  observed',
     1 '    terrain   free-air   Bouguer')
	if(form) write(11,46) junit(ju1),den
46	format(37x,'gravity   inner outer       anomaly',/,
     1 8x,'deg  min   deg  min   ',a5,'    mGal     mGal  mGal',
     1 4x,'mGal',3x,f4.2,' gm/cc'/)
c
	dratio = den / 2.67
c
c  begin processing
c  og in the standard ascii file is -900,000.
c
50	read(10, 51, err=97, end=99) rec
51	format( a )
	in = in + 1
	if( rec(80:80) .eq. 'd' .or. rec(80:80) .eq. 'D' ) then
	  idelet = idelet + 1
	  go to 50
	endif
c
	if(normg) then
	  read(rec, 55, err=97, end=99) id, xl, elvin, ogd, 
     1      tc1, tc2, spec
55	  format(a8,14x,f2.0,f4.2,f4.0,f4.2,f6.2,2x,d7.2,3f6.2)
	else
	  read(rec, 56, err=97, end=99) id, xl, elvin, ogd,
     1      tc1,tc2,spec
56	  format(a8,14x,f2.0,f4.2,f4.0,f4.2,f6.2,   d9.3,3f6.2)
	endif
c
c  elevation precision .01 metre (default) or .1 foot
	em = elvin
	if( iu .eq. 1 ) em = elvin / .328084
	phi = xl(1) + sign( xl(2), xl(1) ) * 1.666667e-2
	call g67( phi, em, tg98, fa, boug, cc )
c  note og has 980000. subtracted for 32 bit machines
	og98  = sngl( ogd - 80000.d0 )
	faye  = og98 - ( tg98 - fa )
	terr  = boug + cc - tc1 - tc2
	banom = faye - terr*dratio - spec
c  bouguer anomaly = og - (tg-fa+(boug-tc)+spec)
c
	if(bin) then
	  xlon=xl(3)+sign(xl(4),xl(3))*1.666667e-2
	  tc=tc1+tc2
	  elvout=em
	  if(iu.eq.1) elvout=10.0*elvin
	  write(12) id,xlon,phi,faye,banom,elvout,tc1,tc2,og98
	endif
c
	if( form ) then
	  elev = em
	  if( ju .eq. 1 ) elev = em * 3.28084
	  lat = int( xl(1) )
	  lon = int( xl(3) )
	  write(11,ofmt) id,lat,xl(2),lon,xl(4),elev,ogd,
     1 tc1,tc2,faye,banom
	endif
c
	iout = iout + 1
	go to 50
c
97	write(*,  98)id
98	format(' error at station ',a8)
	if(in-iout.lt.10) go to 50
	write(*,  96)
96	format(' error count exceeded')

99	close (10)
	if( form )         close( 11 )
	if( ifunc .ge. 2 ) close( 12 )

	write(*,  101)in,iout
101	format(i7,' stations input and ',i7,' stations output')
	write(*, 102) idelet
102	format(i7,' stations deleted')
	return
	end
c******************************************************************************
	subroutine g67( phi, elev, tg98, fa, boug, cc )
c       calculate factors appearing in the bouguer anomaly
c       phi in decimal degree, elev in metre.
c       tg98, fa, boug, cc in mGal.
c  theoretic gravity reference:
c  Systeme Geodesique de Reference 1967, International Assoc of Geodesy
c    special publication #3.  
c  curvature correction reference:
c  W.D. Lambert, 1930, The reduction of observed values of gravity to
c    sea level, Bulletin Geodesique, vol 26, p 107-181.
c  fa and cc coded by G.I. Evenden, USGS.

	double precision p, p2, f1, f2, f3, tg
c        data gcon/6.67e-3/, bcon/.1118966/
        data bcon/.1118966/
	data f1,f2,f3 /978031.85d0, 5.278895d-3, 2.3462d-5/

	p    = dsin( dble( phi ) * 1.74532 92519 94329 d-2 )
	p2   = p * p
	tg   = f1 + f1 * p2 * ( f2 + f3 * p2 )
	phi2 = phi * phi * .0001
	fa   = elev * ( .30877 + phi2 * ( -.0013398 + phi2 * ( 
     1             .0013553 + phi2 * ( -.0005329 + phi2 *
     1              .0000911 ))) - elev * .072e-6 )
	cc   = elev * ( 1.4639108e-3 + elev * (4.449648e-14 * 
     1               elev - 3.532715e-7 ) )
	boug = bcon * elev
	tg98 = sngl( tg - 980000.d0 )
	return
	end
c******************************************************************************
	subroutine getfn(ifile,suffix,ofile)
c  strip off old extension and specified device/directory and add new suffix.
c  VAX filename conventions:    device:[directory]filename.extension
c  directory is level(n).level(n+1).level(n+2)...
	character ifile*(*),suffix*(*),ofile*(*)
c
	lenin = len(ifile)
	lenou = len(ofile)
	ofile = ' '
c  strip off device/directory and default to the working directory
	m = index(ifile,']') + 1
	n = lenin - m + 1
	ofile(1:n) = ifile(m:lenin)
c
c  find nonblank length
10	lennb=0
	do 20 i = n, 1, -1
	  if(ofile(i:i).ne.' ' .and. lennb.eq.0) lennb=i
20	continue
	if(lennb.eq.0) then
	  write(*,*)' getfn: cannot parse filename, enter new prefix'
	  write(*,*)'        to go with the suffix ',suffix
	  read(5,*,err=10) ofile
	  go to 10
	endif
c
c  terminate prefix with period.
	n = index(ofile,'.')
	if (n .eq. 0) then
	  n = lennb + 1
	  ofile(n:n)='.'
	endif	
c
	m = n + 1
	lens = len(suffix)
	n = n + lens
	ofile(m:n) = suffix(1:lens)
c
	do 30 i = n+1, lenou
	  ofile(i:i) = ' '
30	continue
	return
	end
c******************************************************************************
        function noyes()
	character inp*20
	ic=0
	noyes=-1
1	continue
	read(5,2,err=5) inp
2	format(a20)
	if(ic.ge.3) go to 9
4	if(inp(1:1).eq.'y' .or. inp(1:1).eq.'Y') noyes=1
	if(inp(1:1).eq.'n' .or. inp(1:1).eq.'N') noyes=0
	if(noyes.gt.-1) return
5	write(*,6)
6	format(' y or n:'$)
	ic=ic+1
	go to 1
9	write(*,*)' count exceeded, answering no'
	noyes=0
	return
	end
c******************************************************************************
	function leftj(a)
c  left justifies a string and returns the position
c  of the last nonblank character
	character a*(*)
	n=len(a)
	if(a(1:1).ne.' ') go to 15
	do 1 m=2,n
1	if(a(m:m).ne.' ') go to 5
	leftj=0
	return
c
5	i2=0
	do 10 i=m,n
	i2=i2+1
	a(i2:i2)=a(i:i)
10	continue
	do 11 i3=i2+1,n
11	a(i3:i3)=' '
	n=n-m+1
c
15	do 20 leftj=n,1,-1
20	if(a(leftj:leftj).ne.' ') go to 25
25	return
	end
c******************************************************************************
	subroutine cvc(a,i2)
c  i2=1 for lower case to upper
	dimension ilb(2),iub(2),icv(2)
	character a*(*)
	data ilb/97,65/, iub/122,90/, icv/-32,+32/
	n=len(a)
	if(i2.eq.1 .or. i2.eq.2) then
	do 10 i=1,n
	  j=ichar(a(i:i))
	  if(j.lt.ilb(i2) .or. j.gt.iub(i2)) go to 10
	  j=j+icv(i2)
	  a(i:i)=char(j)
10	continue
	else
	write(*,*)' cvc: invalid convert parameter'
	endif
	return
	end
	subroutine gsort
c  sorting routine for the gravity reduction system
c  m webring, 6/83
	common /gfiles/ icont,itty,jtty,f1,f2,f3,f4,f5
        dimension is(11), iw(11), is1(2), iw1(2)
	character*56 f1,f2,f3,f4,f5,tmp1,tmp2
	character*4 suff(11)
	character*8 sort, srtkey(11)
	character idata*80, title1*80, title2*80,
     1 title3*33, title4*16
c       added 5-9-88 for rosort
        integer*2 ipoint(5000),n
        character*12 ikey(5000)
c
	data suff/'id', 'tme', 'lat',  'lon', 'elv', 'grv',
     1 'itc', 'otc', 'spc', 'cde1', 'cde2' /
c
	data srtkey/ 'srt=id',  'srt=tme', 'srt=lat', 'srt=lon',
     1 'srt=elv', 'srt=grv', 'srt=itc', 'srt=otc', 
     1 'srt=spc', 'srt=c1',  'srt=c2'/
c
        data tmp1/'sort1.tmp'/,tmp2/'sort2.tmp'/
c        , iacend,ichrec/0,0/
	data is/1, 9, 22, 29, 37, 45, 52, 58, 64, 70, 75/
	data iw/8, 9,  7,  8,  6,  7,  6,  6,  6,  5,  6/
c
	maxkey = 11
	ier    = 1
	ifld   = 2
	if( icont .gt. 0 ) go to 135
	write(jtty,100)
100	format(' Enter your gravity filename  :'$)
	read(5,37) f1
37	format(a56)
109	write(jtty,110)
110	format('      sort field number',/,
     1 ' 1=id,  2=tim, 3=lat,  4=lon,   5=elv, 6=grv',/,
     1 ' 7=itc, 8=otc, 9=spc, 10=cde1, 11=cde2 :'$)
	read(itty,*) ifld
	if( ifld.lt.1 .or. ifld.gt.maxkey ) go to 109
c
135	continue
	call getfn(f1,suff(ifld),f4)
	if(icont.eq.0) write(*,*)'output file is ',f4
	nkey = 2
	if( ifld .eq. 3 ) then
	is1(1) = is(3)
	is1(2) = is(4)
	iw1(1) = iw(3)
	iw1(2) = 4
	else if( ifld .eq. 4 ) then
	is1(1) = is(3)
	is1(2) = is(4)
	iw1(1) = 3
	iw1(2) = iw(4)
	else
	nkey   = 1
	is1(1) = is(ifld)
	iw1(1) = iw(ifld)
	endif
c
c  generate temp file without title records
	open(10, file=f1, form='formatted', status='old')
c 5-9-88 Changed unit 11 to d.a. file for use in sort
        lrec=80
	open(11, file=tmp1, form='unformatted', status='unknown',
     1 access='direct',recl=4*lrec)
	icnt=1
	read(10,179,end=300) title1
179	format(a80)
	read(10,179,end=300) title2
	read(10,180,end=300) title3, sort, title4
180	format(a33, a8, a16)
	icnt = 1
140	read(10,179,end=160) idata
	write(11,rec=icnt) idata 
        if (icnt .gt. 5000) then
        write (*,*) ' Must increase arrays ipoint & ikey beyond 5000'
        stop
        endif 
        ipoint(icnt)=icnt
c       store sorting keys
c       each character key must be 12 of fewer chars
        ikey(icnt)=idata(is1(1):is1(1)+iw1(1)-1)
        if (nkey .eq.2) then
        ikey(icnt)=ikey(icnt)(1:iw1(1))
     1  //idata(is1(2):is1(2)+iw1(2)-1)
15      continue
        endif
	icnt = icnt + 1
	go to 140
160	close(10)
	close(11)        
        icnt=icnt-1
	if( icnt .eq. 0 ) go to 260
	if( icnt .gt. 1 ) go to 170
	tmp2 = tmp1
	go to 175
c
c170	call vaxsort( tmp1, tmp2, nkey, is1, iw1, iacend, ichrec )
c       substitute Grundy's quicker sort
170     n=icnt
        call rosort(ikey,ipoint,n)
c
c  generate sorted data file
175	open(11,file=tmp1,form='unformatted',status='old',
     1  access='direct',recl=4*lrec)
	open(13,file=f4,form='formatted',status='unknown')
	write(13,179) title1
	write(13,179) title2
	write(13,201) title3,srtkey(ifld),title4
201	format(a33,a8,a16)
	do 210 i = 1, icnt
	read(11,rec=ipoint(i)) idata
	write(13,179) idata
210	continue
220	close(13)
c	open(10,file=tmp1,status='old')
c	close(10)
	write(jtty,230) icnt
230	format(/,1x,i7,' stations sorted')
	return
c
300	write(jtty,350) i
350	format(' grsort: eof at record ',i7)
	icont=0
	go to 999
260	write(jtty,270)
270	format(' grsort: no data records')
	icont=0
999	close(10)
	close(11)
	return
	end
      SUBROUTINE ROSORT(CHAR8,IPOINT,N)
      implicit  integer *2 (i-n)
      CHARACTER*12 CHAR8
C***  THIS IS A FAST ROUTINE FOR SORTING A SET OF POINTERS
C***  CONTAINED IN ARRAY IPOINT BASED UPON INCREASING VALUES
C***  IN CHARACTER ARRAY CHAR8.  IT IS SIMILAR TO FSORT.FOR.
C***  ROUTINE ORIGINALLY WRITTEN BY CALVIN SMITH, U. S. ATOMIC
C***  ENERGY COMMISSION, GRAND JUNCTION OFFICE, FOR A CDC 3100
C***  COMPUTER.  THIS SORT IS NOT STABLE.
C***  RECODED FOR IBM PERSONAL COMPUTER BY W. D. GRUNDY APRIL 1984
C***
C***
C***  DIMENSION REQUIREMENTS: 
C***  IF CHAR8 IS DIMENSIONED AS CHAR8(M) IN THE MAIN PROGRAM, SET
C***  DIMENSIONS IU(N) AND IL(N) SUCH THAT 2**(N+1)-1 IS 
C***  GREATER THAN OR EQUAL TO M.
C***
      DIMENSION CHAR8(1),IPOINT(1),IU(15),IL(15)
      I=1
      M=1
      J=N
   10 IF(I.GE.J) GO TO 110
   20 K=I
      L=J
      IJ=(I+J)/2
      IF(CHAR8(IPOINT(I)).LE.CHAR8(IPOINT(IJ))) GO TO 30
      ITEMP=IPOINT(I)
      IPOINT(I)=IPOINT(IJ)
      IPOINT(IJ)=ITEMP
   30 IF(CHAR8(IPOINT(J)).GE.CHAR8(IPOINT(IJ))) GO TO 60
      ITEMP=IPOINT(J)
      IPOINT(J)=IPOINT(IJ)
      IPOINT(IJ)=ITEMP
      IF(CHAR8(IPOINT(IJ)).GE.CHAR8(IPOINT(I))) GO TO 60
      ITEMP=IPOINT(I)
      IPOINT(I)=IPOINT(IJ)
      IPOINT(IJ)=ITEMP
      GO TO 60
   40 ITEMP=IPOINT(L)
      IPOINT(L)=IPOINT(K)
      IPOINT(K)=ITEMP
      IF(IJ.NE.L) GO TO 50
      IJ=K
      GO TO 60
   50 IF(IJ.NE.K) GO TO 60
      IJ=L
   60 L=L-1
      IF(CHAR8(IPOINT(IJ)).LT.CHAR8(IPOINT(L))) GO TO 60
   70 K=K+1
      IF(CHAR8(IPOINT(K)).LT.CHAR8(IPOINT(IJ))) GO TO 70
      IF(K-L) 40,80,90
   80 K=K+1
      L=L-1
   90 IF((L-I).LE.(J-K)) GO TO 100
      IL(M)=I
      IU(M)=L
      I=K
      M=M+1
      GO TO 120
  100 IL(M)=K
      IU(M)=J
      J=L
      M=M+1
      GO TO 120
  110 M=M-1
      IF(M.EQ.0) RETURN
      I=IL(M)
      J=IU(M)
  120 IF((J-I).GE.11) GO TO 20
      IF(I.EQ.1) GO TO 10
      GO TO 140
  130 I=I+1
  140 IF(I.EQ.J) GO TO 110
      IF(CHAR8(IPOINT(I)).LE.CHAR8(IPOINT(I+1))) GO TO 130
      K=I
  150 ITEMP=IPOINT(K)
      IPOINT(K)=IPOINT(K+1)
      IPOINT(K+1)=ITEMP
      K=K-1
      IF(CHAR8(IPOINT(K+1)).LT.CHAR8(IPOINT(K))) GO TO 150
      GO TO 130
      END
