c  pdepth - magnetic and gravity profile analysis and SAKI model building.
c           Input profile and model data are in SAKI format.
c
c           written by J.Phillips
c
      common /menu1/ menux,menudy
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
      common /plot/ dxp(2),dyp(2),xpix,ypix,xscale,yscale
      common /parms/ bfile,mfile,grfile,hfile,efield,einc,edec,jplotr,
     1               iave,jave
      common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      common /pscale/ xmin,xmax,zmin,zmax
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /pmagdata/ pmmin,pmmax,pmminsav,pmmaxsav,pminc,pm(512)
      common /pgrvdata/ pgmin,pgmax,pgminsav,pgmaxsav,pginc,pg(512)
      common /rmagdata/ rmin,rmax,rminsav,rmaxsav,ra(512)
      common /contdata/ nco,comin,comax,cominsav,comaxsav,coinc,zup,
     1co(512)
      common /angles/ pmaginc,pmagdec,continc,contdec
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      common /calc1/icalc,icalcd,gc(512),hc(512)
      common /save/ gminsav,gmaxsav,aminsav,amaxsav
      dimension mpen(7)
      character ifile*80, label(7)*8, prompt*79, gfile*80, ans*1
      character sfile*80
      character*56 bfile,mfile,grfile,hfile
      data nmenu/7/,label/' profile','  model ',' window ','  reset ',
     &' geology','  clear ','  exit  '/
      data mpen/0,0,0,0,0,0,0/
      data imag/0/,igrv/0/,irmag/0/,ipmag/0/,ipgrv/0/,icont/0/
      data icalc/0/,icalcd/0/,ef(1)/0.0/
c
        call namepd
        ifile=mfile
        iplotr=jplotr
        if(iplotr.ne.9.and.iplotr.ne.10) then
   21     iplotr= ask_int('Enter plot device code (9=EGA,10=VGA): ')
        endif
        if(iplotr.eq.9) then
          nhor=640
          tspace=40
          nvert=350
        else if(iplotr.eq.10) then
          nhor=640
          tspace=55
          nvert=480
        else
          print *,'try again'
          go to 21
        endif
      iopn=0
      icolor=500
      icsw=1
c      imag=0
      if(hfile(1:1).eq.' ') then
  211   write(6,201)
  201   format(' enter magnetic profile file name (or return): '\)
        read(5,200,end=8,err=211) hfile
  200   format(a)
      endif
    8 if(hfile(1:1).ne.' ') then
        open(10,file=hfile,status='old',form='formatted')
        imag=1
        amax=-1.e+38
        amin=1.e+38
        zmax=-1.e+38
        zmin=1.e+38
        n=0
   10   n=n+1
        if(n.gt.500) then
          print*,'Warning: using first 500 points of profile'
          go to 12
        endif
        read(10,*,end=12) x(n),z(n),a(n)
        if(a(n).gt.amax) amax=a(n)
        if(a(n).lt.amin) amin=a(n)
        if(z(n).gt.zmax) zmax=z(n)
        if(z(n).lt.zmin) zmin=z(n)
        go to 10
   12   n=n-1
        close(10)
        xmin=x(1)
        xmax=x(n)
        delx=(xmax-xmin)/float(n-1)
        xsign=1.0
        if(delx.lt.0.) then
          print*,'Warning: sign of x-coordinates has been changed'
          xmin=-xmin
          xmax=-xmax
          xsign=-1.0
          x(1)=-x(1)
          delx=-delx
        endif
        do 13 i=2,n
        x(i)=xsign*x(i)
   13   if(abs(x(i)-x(i-1)-delx).gt.delx/10.) imag=2
        if(amax.eq.amin) then
          amax=amax+10.
          amin=amin-10.
        endif
      endif
c      igrv=0
      if(grfile(1:1).eq.' ') then
  212   write(6,202)
  202   format(' enter gravity profile file name (or return): '\)
        read(5,200,end=14,err=212) grfile
      endif
   14 if(grfile(1:1).ne.' ') then
        open(10,file=grfile,status='old',form='formatted')
        igrv=1
        gmax=-1.e+38
        gmin=1.e+38
        gzmax=-1.e+38
        gzmin=1.e+38
        ng=0
   16   ng=ng+1
        if(ng.gt.500) then
          print*,'Warning: using first 500 points of profile'
          go to 20
        endif
        read(10,*,end=20) gx(ng),gz(ng),g(ng)
        if(g(ng).gt.gmax) gmax=g(ng)
        if(g(ng).lt.gmin) gmin=g(ng)
        if(gz(ng).gt.gzmax) gzmax=gz(ng)
        if(gz(ng).lt.gzmin) gzmin=gz(ng)
        go to 16
   20   ng=ng-1
        close(10)
        xmin=gx(1)
        xmax=gx(ng)
        delx=(xmax-xmin)/float(ng-1)
        xsign=1.0
        if(delx.lt.0.) then
          print*,'Warning: sign of x-coordinates has been changed'
          xmin=-xmin
          xmax=-xmax
          xsign=-1.0
          gx(1)=-gx(1)
          delx=-delx
        endif
        do 22 i=2,ng
        gx(i)=xsign*gx(i)
   22   if(abs(gx(i)-gx(i-1)-delx).gt.delx/10.) igrv=2
        if(gmax.eq.gmin) then
          gmax=gmax+10.
          gmin=gmin-10.
        endif
      endif
c      if(imag.eq.0.and.igrv.eq.0) stop 'at least one profile required'
      if(imag.eq.0.and.igrv.eq.0) then
  213   write(6,203)
  203 format(' enter startx, delx, nobs, magalt, grvalt: '\)
        read(5,*,err=213) xmin,delx,n,zmin,gzmin
  214   write(6,204)
  204 format(' enter m, g, or b (for mag, grav, or both): '\)
        read(5,200,err=214) ans
        if(ans.eq.'m'.or.ans.eq.'M'.or.ans.eq.'b'.or.ans.eq.'B') then
          imag=1
          do 215 i=1,n
          x(i)=xmin+(i-1)*delx
          a(i)=0.
          z(i)=zmin
  215     continue
          amin=-10.
          amax=10.
c          zmin=0.
          zmax=zmin
        else
          if(ans.ne.'g'.and.ans.ne.'G') go to 214
        endif
        if(ans.eq.'g'.or.ans.eq.'G'.or.ans.eq.'b'.or.ans.eq.'B') then
          igrv=1
          ng=n
          do 216 i=1,ng
          gx(i)=xmin+(i-1)*delx
          g(i)=0.
          gz(i)=gzmin
  216     continue
          gmax=10.
          gmin=-10
c          gzmin=0.
          gzmax=gzmin
        endif
      endif
      if(imag.eq.0) then
        zmin=gzmin
        zmax=gzmax
        xmin=gx(1)
        xmax=gx(ng)
        delx=(xmax-xmin)/float(ng-1)
        call dplot(0,xmin,xmax,gmin,gmax)
        ipen=300
        call curv(gx,g,ng,ipen)
      else if(igrv.gt.0) then
        zmin=amin1(zmin,gzmin)
        zmax=amax1(zmax,gzmax)
        xmin=amin1(x(1),gx(1))
        xmax=amax1(x(n),gx(ng))
        delx=(xmax-xmin)/float(ng-1)
        call dplot(0,xmin,xmax,gmin,gmax)
        ipen=300
        call curv(gx,g,ng,ipen)
        call dplot(2,xmin,xmax,amin,amax)
        ipen=400
        call curv(x,a,n,ipen)
      else
        xmin=x(1)
        xmax=x(n)
        delx=(xmax-xmin)/float(n-1)
        call dplot(0,xmin,xmax,amin,amax)
        ipen=400
        call curv(x,a,n,ipen)
      endif
      xminsav=xmin
      xmaxsav=xmax
      amaxsav=amax
      aminsav=amin
      gmaxsav=gmax
      gminsav=gmin
c
      call dplot(1,xmin,xmax,zmin,zmax)
      ipen=400
      if(imag.gt.0) call curv(x,z,n,ipen)
      ipen=300
c      if(igrv.gt.0) call curv1(gx,gz,ng,ipen,2)
      if(igrv.gt.0) call curv(gx,gz,ng,ipen)
      call rescale
      call menu(nmenu,label,mpen)
      ixo=600
      iyo=300
      gfile=' '
      sfile=' '
      ibod=0
      iexit=0
      prompt=' '
      if(imag.eq.2) prompt='Warning: magnetic profile is unevenly sample
     1d'
      if(igrv.eq.2) prompt='Warning: gravity profile is unevenly sampled
     1'
      if(imag.eq.2.and.igrv.eq.2) prompt='Warning: both profiles are une
     1venly sampled'
      call tpromp(0,0,prompt)
  100 call menudvr(xdum,y,ifunc,ixo,iyo)
      go to (100,103,104,101,102,105,106,107) ifunc+1
c window
  101 prompt='click on left limit'
      call tpromp(0,0,prompt)
      call menudvr(xmin,y,ifunc,ixo,iyo)
      if(ifunc.ne.0) go to 163
c      call tgin(ix,iy,ixo,iyo)
c      call tgin2(2,0,0,0)
c      xmin=dxp(1)+ix*xpix
      prompt='click on right limit'
      call tpromp(0,0,prompt)
      call menudvr(xmax,y,ifunc,ixo,iyo)
      if(ifunc.ne.0) go to 163
c      call tgin(ix,iy,ixo,iyo)
c      call tgin2(2,0,0,0)
c      xmax=dxp(1)+ix*xpix
      if(xmax.le.xmin) go to 101
c      print *, xmin, xmax
c      call endit
      comax=-1.e38
      comin=1.e38
      if(igrv.gt.0) then
        gmax=-1.e38
        gmin=1.e38
        pmmax=-1.e38
        pmmin=1.e38
        do 150 i=1,ng
        if(gx(i).lt.xmin) go to 150
        if(gx(i).gt.xmax) go to 151
        if(g(i).gt.gmax) gmax=g(i)
        if(g(i).lt.gmin) gmin=g(i)
        if(ipmag.eq.1) then
          if(pm(i).gt.pmmax) pmmax=pm(i)
          if(pm(i).lt.pmmin) pmmin=pm(i)
        endif
        if(icont.gt.0.and.nco.eq.ng) then
          if(co(i).gt.comax) comax=co(i)
          if(co(i).lt.comin) comin=co(i)
        endif
  150   continue
  151   if(gmin.ge.gmax) then
          gmin=-10.
          gmax=10.
          pmmin=-10.
          pmmax=10.
          comin=-10.
          comax=10.
        endif
      endif
      if(imag.gt.0) then
        amax=-1.e38
        amin=1.e38
        rmax=-1.e38
        rmin=1.e38
        pgmax=-1.e38
        pgmin=1.e38
        do 152 i=1,n
        if(x(i).lt.xmin) go to 152
        if(x(i).gt.xmax) go to 153
        if(a(i).gt.amax) amax=a(i)
        if(a(i).lt.amin) amin=a(i)
        if(irmag.eq.1) then
          if(ra(i).gt.rmax) rmax=ra(i)
          if(ra(i).lt.rmin) rmin=ra(i)
        endif
        if(ipgrv.eq.1) then
          if(pg(i).gt.pgmax) pgmax=pg(i)
          if(pg(i).lt.pgmin) pgmin=pg(i)
        endif
        if(icont.gt.0.and.nco.eq.n) then
          if(co(i).gt.comax) comax=co(i)
          if(co(i).lt.comin) comin=co(i)
        endif
  152   continue
  153   if(amin.ge.amax) then
          amin=-10.
          amax=10.
          rmin=-10.
          rmax=10.
          pgmin=-10.
          pgmax=10.
          comin=-10.
          comax=10.
        endif
      endif
  120  call setmod(3)
      icalcd=0
      if(igrv.gt.0) then
        ipen=300
        call dplot(0,xmin,xmax,gmin,gmax)
        call curv(gx,g,ng,ipen)
        if(imag.gt.0) then
          call dplot(2,xmin,xmax,amin,amax)
          ipen=400
          call curv(x,a,n,ipen)
        endif
      else
        call dplot(0,xmin,xmax,amin,amax)
        ipen=400
        call curv(x,a,n,ipen)
      endif
      if(irmag.eq.1) then
        call dplot(2,xmin,xmax,rmin,rmax)
        ipen=500
        call curv(x,ra,n,ipen)
      endif
      if(ipmag.eq.1) then
        call dplot(2,xmin,xmax,pmmin,pmmax)
        ipen=600
        call curv(gx,pm,ng,ipen)
      endif
      if(ipgrv.eq.1) then
        call dplot(2,xmin,xmax,pgmin,pgmax)
        ipen=100
        call curv(x,pg,n,ipen)
      endif
      if(icont.gt.0) then
        call dplot(2,xmin,xmax,comin,comax)
        ipen=200
c       if(nco.eq.n) then
       if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
         call curv(x,co,n,ipen)
       else
         call curv(gx,co,ng,ipen)
       endif
      endif
      call dplot(1,xmin,xmax,zmin,zmax)
      ipen=400
      if(imag.gt.0) call curv(x,z,n,ipen)
      ipen=300
c      if(igrv.gt.0) call curv1(gx,gz,ng,ipen,2)
      if(igrv.gt.0) call curv(gx,gz,ng,ipen)
      call rescale
      if(icont.gt.0) then
        ipen=200
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          call curv2(x,z,zup,n,ipen)
        else
          call curv2(gx,gz,zup,ng,ipen)
        endif
      endif
      call menu(nmenu,label,mpen)
      if(iopn.gt.0) call plotmod
      if(gfile(1:1).ne.' ') then
        open(10,file=gfile,form='formatted',status='old',err=163)
        call plotgeo
        close(10)
      endif
      if(sfile(1:1).ne.' ') then
        open(10,file=sfile,form='formatted',status='old',err=163)
        call plotsect
        close(10)
      endif
c      ixo=600
c      iyo=300
      go to 100
c reset window
  102 xmin=xminsav
      xmax=xmaxsav
      amin=aminsav
      amax=amaxsav
      gmin=gminsav
      gmax=gmaxsav
      pmmin=pmminsav
      pmmax=pmmaxsav
      rmin=rminsav
      rmax=rmaxsav
      pgmin=pgminsav
      pgmax=pgmaxsav
      comin=cominsav
      comax=comaxsav
      go to 120
c      if(igrv.eq.1) then
c        ipen=300
c        call dplot(0,xmin,xmax,gmin,gmax,2.)
c        call curv(gx,g,ng,ipen)
c        if(imag.eq.1) then
c          call dplot(2,xmin,xmax,amin,amax,50.)
c          ipen=400
c          call curv(x,a,n,ipen)
c        endif
c      else
c        call dplot(0,xmin,xmax,amin,amax,50.)
c        ipen=400
c        call curv(x,a,n,ipen)
c      endif
c      if(irmag.eq.1) then
c        call dplot(2,xmin,xmax,rmin,rmax,50.)
c        ipen=500
c        call curv(x,ra,n,ipen)
c      endif
c      call dplot(1,xmin,xmax,zmin,zmax,1.)
c      if(imag.eq.1) then
c        ipen=400
c        call curv(x,z,n,ipen)
c      endif
c      if(igrv.eq.1) then
c        ipen=300
c        call curv(gx,gz,ng,ipen)
c      endif
c      call rescale
c      call menu(nmenu,label)
c      ixo=600
c      iyo=300
c      if(iopn.eq.1) call plotmod
c      if(gfile(1:1).ne.' ') then
c        open(10,file=gfile,form='formatted',status='old',err=163)
c        call plotgeo
c        close(10)
c      endif
c      go to 100
c profile manipulation
  103 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
c      call setedt(ixo,iyo)
      call setpro(ixo,iyo)
      icolor=500
      call menu(nmenu,label,mpen)
      go to 100
c model manipulation
  104 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
c     call setdep(ixo,iyo)
      call setmodel(ixo,iyo)
      icolor=500
      call menu(nmenu,label,mpen)
      go to 100
c geology file
  105 prompt='enter name of geology file: '
      call tpromp(0,0,prompt)
      gfile=' '
      read(5,200,end=31) gfile
   31 if(gfile(1:1).eq.' ') go to 162
      open(10,file=gfile,form='formatted',status='old',err=163)
      call plotgeo
      close(10)
  162 prompt='enter name of cross section file: '
      call tpromp(0,0,prompt)
      sfile=' '
      read(5,200,end=31) sfile
      if(sfile(1:1).eq.' ') go to 163
      open(10,file=sfile,form='formatted',status='old',err=163)
      call plotsect
      close(10)
      call menu(nmenu,label,mpen)
  163 prompt=' '
      call tpromp(0,0,prompt)
c      print *,ifunc,iy,nvert,tspace,nmenu,menudy
      go to 100
c clear
  106 call clear
      call menu(nmenu,label,mpen)
      go to 100
c exit
  107 prompt='click on exit again to confirm exit'
      call tpromp(0,0,prompt)
      call menudvr(xdum,y,ifunc,ixo,iyo)
      prompt=' '
      call tpromp(0,0,prompt)
      go to (100,103,104,101,102,105,106,108) ifunc+1
c really exit
  108 if(iopn.eq.3) then
        prompt='save changes to model file? [y] '
        call tpromp(0,0,prompt)
        read(5,200) ans
        if(ans.eq.' ') ans='y'
        if(ans.eq.'y'.or.ans.eq.'Y') then
          call writemod
        else
          iopn=0
        endif
      else if(iopn.eq.1) then
        call writemod
      endif
      if(icalc.gt.0.or.igrv.gt.0.or.irmag.gt.0.or.ipmag.gt.0.or.ipgrv
     1.gt.0) then
        prompt='output calculated profiles? [n] '
        call tpromp(0,0,prompt)
        read(5,200) ans
        if(ans.eq.' ') ans='n'
        if(ans.eq.'y'.or.ans.eq.'Y') go to 117
        icalc=0
        igrv=0
        irmag=0
        ipmag=0
        ipgrv=0
        go to 116
      endif
  117 if(icalc.gt.0) then
        if(imag.gt.0) then
          open(10,file='hcalc.out',form='formatted',status='new',
     1      err=109)
          go to 110
  109     prompt='hcalc.out exists, OK to overwrite? [n] '
          call tpromp(0,0,prompt)
          read(5,200) ans
          if(ans.eq.' ') ans='n'
          if(ans.ne.'y'.and.ans.ne.'Y') then
            imag=0
            go to 112
          endif
          prompt=' '
          call tpromp(0,0,prompt)
          open(10,file='hcalc.out',form='formatted',status='unknown')
  110     continue
          do 111 i=1,n
          if(icalc.eq.1) then
            write(10,205) x(i),z(i),hc(i)
          else
            write(10,205) x(i),z(i),a(i)
          endif
  111     continue
          close(10)
        endif
  112   if(igrv.gt.0) then
          open(10,file='gcalc.out',form='formatted',status='new',
     1      err=113)
          go to 114
  113     prompt='gcalc.out exists, OK to overwrite? [n] '
          call tpromp(0,0,prompt)
          read(5,200) ans
          if(ans.eq.' ') ans='n'
          if(ans.ne.'y'.and.ans.ne.'Y') then
            igrv=0
            go to 408
          endif
          prompt=' '
          call tpromp(0,0,prompt)
          open(10,file='gcalc.out',form='formatted',status='unknown')
  114     continue
          do 115 i=1,ng
          if(icalc.eq.1) then
            write(10,205) gx(i),gz(i),gc(i)
          else
            write(10,205) gx(i),gz(i),g(i)
          endif
  115     continue
          close(10)
        endif
      endif
  408 if(irmag.gt.0) then
        open(10,file='rcalc.out',form='formatted',status='new',
     1    err=409)
        go to 410
  409   prompt='rcalc.out exists, OK to overwrite? [n] '
        call tpromp(0,0,prompt)
        read(5,200) ans
        if(ans.eq.' ') ans='n'
        if(ans.ne.'y'.and.ans.ne.'Y') then
          irmag=0
          go to 412
        endif
        prompt=' '
        call tpromp(0,0,prompt)
        open(10,file='rcalc.out',form='formatted',status='unknown')
  410   continue
        do 411 i=1,n
        write(10,205) x(i),z(i),ra(i)
  411   continue
        close(10)
      endif
  412 if(ipmag.gt.0) then
        open(10,file='pmcalc.out',form='formatted',status='new',
     1    err=413)
        go to 414
  413   prompt='pmcalc.out exists, OK to overwrite? [n] '
        call tpromp(0,0,prompt)
        read(5,200) ans
        if(ans.eq.' ') ans='n'
        if(ans.ne.'y'.and.ans.ne.'Y') then
          ipmag=0
          go to 416
        endif
        prompt=' '
        call tpromp(0,0,prompt)
        open(10,file='pmcalc.out',form='formatted',status='unknown')
  414   continue
        do 415 i=1,ng
        write(10,205) gx(i),gz(i),pm(i)
  415   continue
        close(10)
      endif
  416 if(ipgrv.gt.0) then
        open(10,file='pgcalc.out',form='formatted',status='new',
     1    err=309)
        go to 310
  309   prompt='pgcalc.out exists, OK to overwrite? [n] '
        call tpromp(0,0,prompt)
        read(5,200) ans
        if(ans.eq.' ') ans='n'
        if(ans.ne.'y'.and.ans.ne.'Y') then
          ipgrv=0
          go to 116
        endif
        prompt=' '
        call tpromp(0,0,prompt)
        open(10,file='pgcalc.out',form='formatted',status='unknown')
  310   continue
        do 311 i=1,n
        write(10,205) x(i),z(i),pg(i)
  311   continue
        close(10)
      endif
  116 call setmod(3)
      print*
      if(iopn.eq.1) then
        len=len_trim(ifile)
        print*,'new model written to ',ifile(1:len)
      endif
      if(iopn.eq.3) then
        len=len_trim(ifile)
        print*,'edited model written to ',ifile(1:len)
        lend=index(ifile,'.')
        if(ifile(lend:len).ne.'.bak') then
          print*,'   old model written to ',ifile(1:lend)//'bak'
c        else
c          print*,'old model written to ',ifile(1:lend)//'bk1'
        endif
        print*
      endif
      if(icalc.gt.0) then
        if(imag.gt.0) then
          print*,'calculated magnetic profile written to hcalc.out'
          print*
        endif
        if(igrv.gt.0) then
          print*,'calculated gravity  profile written to gcalc.out'
          print*
        endif
      endif
      if(irmag.gt.0) then
        print*,'reduced-to-pole magnetic profile written to rcalc.out'
        print*
      endif
      if(ipmag.gt.0) then
        print*,'pseudomagnetic profile written to pmcalc.out'
        print*
      endif
      if(ipgrv.gt.0) then
        print*,'pseudogravity profile written to pgcalc.out'
        print*
      endif
      stop
  205 format(3g16.6)
      end
c******************************************************************************
c      subroutine dplot(iview, xmin, xmax, delx, zmin, zmax, aint0)
      subroutine dplot(iview, xmin, xmax, zmin, zmax)
      common /menu1/ menux,menudy
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
      common /plot/ dxp(2),dyp(2),xpix,ypix,xscale,yscale
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      dimension xp(4), yp(4)
      dimension x(n), y(n), mpen(nbox)
      character label(nbox)*8, dfmt*6, afmt*6, gfile*80, sfile*80
      if(iview.eq.0) then
        call pltset(iplotr, xp(4), yp(4), 1)
      endif
      xp(4) = min(xp(4),10.)
      yp(4) = min(yp(4),8.)
      dxp(1) = xmin
      dxp(2) = (xmax + .001) - mod(xmax,.001)
      dyp(1) = zmin
      dyp(2) = zmax
      xp(1) = xp(4) - 2.1
      xp(2) = 0.
      xp(3) = 0.9
      zint=(xmax-xmin)*(yp(4)/2.-0.2)/xp(1)
      if(zint.lt.0.01) then
        dfmt='(f6.4)'
        xint=.0005
      else if(zint.lt.0.1) then
        dfmt='(f6.3)'
        xint=.005
      else if(zint.lt.1.) then
        dfmt='(f6.2)'
        xint=0.05
      else if(zint.lt.10.) then
        dfmt='(f6.1)'
        xint=0.5
      else if(zint.lt.100.) then
        dfmt='(f7.0)'
        xint=5.
      else
        dfmt='(f7.0)'
        xint=50.
      endif
      if(iview.ne.1) then
        yp(1) = yp(4)/2. - 1.2
        yp(2) = 0.
        yp(3) = yp(4)/2. + 0.9
      else
        yp(1) = yp(4)/2. - 0.2
        yp(2) = 0.
        yp(3) = 0.2
        dyp(1) = (xmax-xmin)*yp(1)/xp(1)+zmin
        zmax = dyp(1)
        dyp(2) = zmin
        xscale=(dxp(2)-dxp(1))/xp(1)
        yscale=(dyp(2)-dyp(1))/yp(1)
      endif
      zint=abs(zmax-zmin)
      if(zint.lt.0.01) then
        afmt='(f6.4)'
        aint=0.0005
      else if(zint.lt.0.1) then
        afmt='(f6.3)'
        aint=0.005
      else if(zint.lt.1.) then
        afmt='(f6.2)'
        aint=0.05
      else if(zint.lt.10.) then
        afmt='(f6.1)'
        aint=0.5
      else if(zint.lt.100.) then
        afmt='(f7.0)'
        aint=5.
      else
        afmt ='(f7.0)'
        aint=50.
      endif
      call scale(dxp, dyp, xp, yp, 4, ier)
      if (ier .ne. 0) goto 90
      if(iview.ne.2) then
c      print*,'dyp = ',dyp,' aint = ',aint,' afmt = ',afmt
c      print*,'dxp = ',dxp,' xint = ',xint,' dfmt = ',dfmt
        call yaxis(dyp, dxp, yp, aint, 2, .12, afmt, 6)
        if(iview.eq.0) then
          ipen=302
          if(igrv.eq.0) ipen=402
          call xaxis(dxp, dyp, xp, xint, 2, .12, dfmt, 6)
          call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'distance',8, 302,
     &    .12, 0., -.54, -.54)
          call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'amplitude', 9, ipen,
     &    .12, 1.5706, -.58, .54)
        else
          call xaxis(dxp, dyp, xp, xint, 2, .12, dfmt, 0)
          call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'depth', 5, 302,
     &    .12, 1.5706, -.58, .54)
        endif
        call neatl
      endif
   64 return 
   90 stop
c
      entry curv(x, y, n, m)
      call line(x, y, n, 0, m)
      return
c
      entry curv1(x, y, n, m, isym)
      call line(x, y, n, 0, m)
      do 91 i=1,n
      call vchar(x(i),y(i),isym,1,m+1,.08,0.,0.,0.)
   91 continue
      return
c
      entry curv2(x, y, zup, n, m)
      do 92 i=1,n
   92 y(i)=y(i)-zup
      call line(x, y, n, 0, m)
      do 93 i=1,n
   93 y(i)=y(i)+zup
      return
c
      entry rescale
      xscale=(dxp(2)-dxp(1))/xp(1)
      yscale=(dyp(2)-dyp(1))/yp(1)
      xpix=xscale*xp(4)/nhor
      ypix=yscale*yp(4)/nvert
      xp(1)=xp(4)
      yp(1)=yp(4)
      dxp(1)=dxp(1)-xp(3)*xscale
      dyp(1)=dyp(1)-yp(3)*yscale
      dxp(2)=xscale*xp(4)+dxp(1)
      dyp(2)=yscale*yp(4)+dyp(1)
      xp(3)=0.
      yp(3)=0.
      call scale(dxp,dyp,xp,yp,4,ier)
      if(ier.ne.0) stop 'rescaling error'
      return
c
      entry menu(nbox,label,mpen)
      xleft=dxp(2)-1.0*xscale
      xright=dxp(2)-0.1*xscale
      ybot=0.1*yscale
      call line(xleft,dyp(2),1,0,icolor)
      call line(xright,dyp(2),1,1,icolor)
      call line(xright,dyp(1)+ybot,1,1,icolor)
      call line(xleft,dyp(1)+ybot,1,1,icolor)
      call line(xleft,dyp(2),1,1,icolor)
      dy=(dyp(2)-dyp(1))/float(nbox)
      ybox=dyp(1)+dy
      do 10 i=1,nbox-1
      call line(xleft,ybox,1,0,icolor)
      call line(xright,ybox,1,1,icolor)
      ybox=ybox+dy
   10 continue
      xm=xleft+0.1*xscale
      ym=dyp(2)-dy/2.
      do 20 i=1,nbox
      if(icsw.eq.-1) then
        call vchar(xm,ym,label(i),8,icsw*2,.1,0.,0.,0.)
      else
        call vchar(xm,ym,label(i),8,mpen(i)+2,.1,0.,0.,0.)
      endif
      ym=ym-dy
   20 continue
      menux=(xleft-dxp(1))/xpix
      menudy=nvert/nbox
      return
c
      entry menudvr(xm,ym,ifunc,ixo,iyo)
        ifunc=0
        xm=0.
        ym=0.
        call tgin(ix,iy,ixo,iyo)
        call tgin2(2,0,0,0)
c        print *, ix, iy
        if(ix.gt.menux) then
c          iy=nvert+tspace-iy
c          iy=(nbox+1)*menudy-iy
c          ifunc=nbox + 1 - iy/menudy
           ifunc=iy/menudy + 1
          if(ifunc.gt.nbox) ifunc=nbox
c      print *,ifunc,iy
        else
          iy=nvert-iy
          xm=float(ix)*xpix+dxp(1)
          ym=float(iy)*ypix+dyp(1)
        endif
        return
c
      entry endit()
      call endpt(ie)
      return 
      end
c******************************************************************************
      function ask_int(request)
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(1x,a\)
      read(unit=5, fmt=*, err=1) ask_int
      return 
    1 write(unit=6, fmt=101) 
  101 format(/,20h **ERROR - try again)
      call bell
      goto 2
      end
c******************************************************************************
      subroutine bell()
      character ding*1
      ding = char(7)
      write(unit=*, fmt=*) ding
      return 
      end
c******************************************************************************
      interface to subroutine
     +settextposition[far,c,alias:"__settextposition"](row,column)
      integer*2 row,column
      end
cc      interface to function
cc     +settextcolor[far,c,alias:"__settextcolor"](index)
cc      integer*2 index
cc      end
cc        INTERFACE TO FUNCTION settextcolor(index)
cc        INTEGER*2 settextcolor[FAR,C,ALIAS:"__settextcolor"],index
cc        END
cc        INTERFACE TO SUBROUTINE
cc     +   settextposition[FAR,C,ALIAS:"__f_settextposition"](row,col,s)
cc        INTEGER*2 row,col
cc        STRUCTURE/rccoord/
cc          INTEGER*2 row
cc          INTEGER*2 col
cc        END STRUCTURE
cc        RECORD/rccoord/s[FAR,REFERENCE]
cc        END
cc
c***********************************************************************
        subroutine tpromp(iwrt,icol,prompt)
c
c        include 'fgraph.fi'
c        include 'fgraph.fd'
c
        character*16 ifmt
c        integer*2 idum
        character bpromp*70
           character prompt*79
        data bpromp/'
     1                    '/
        iwrt=0
c        call setcur(line,icol)
        call settextposition(1,1)
c        idum=settextcolor(7)
        write(*,'(a\)') bpromp
c        print'(a\)',bpromp
c        call setcur(line-1,icol)
        call settextposition(1,1)
        if(icol.eq.0) icol=1
        num=len_trim(prompt)
        write(ifmt,10) icol,num
   10   format('(',i2,'x,a',i2,'\)')
        write(*,ifmt) prompt(1:num)
        return
        end
c******************************************************************************
cc        subroutine tpromp(iwrt,line,icol,prompt)
c        subroutine tpromp(iwrt,icol,prompt)
c        character*16 ifmt
c        character bpromp*70
c           character prompt*79
c        data bpromp/' '/
c       iwrt=0
cc        call setcur(line,icol)
c        call setcur(1,icol)
c        write(*,*) bpromp
cc        call setcur(line-1,icol)
c        call setcur(0,icol)
c        if(icol.eq.0) icol=1
cc        num=lench(prompt)
c       num=len_trim(prompt)
cc        write(ifmt,10) icol
c        write(ifmt,10) icol,num
cc   10   format('(',i2,'x,a','\)')
c   10   format('(',i2,'x,a',i2,'\)')
c        write(*,ifmt) prompt(1:num)
c        return
c        end
c***********************************************************************
c        function lench(data)
c        character data*79
c        do 10 i=79,1,-1
c        if(data(i:i).ne.' ') go to 20
c   10   continue
c   20   lench=i
c        return
c        end
c******************************************************************************
        subroutine tgin(ix,iy,ixo,iyo)
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c  gets cursor position
        m1=0
        call mouse(m1,m2,m3,m4)
        m1=4
        call mouse(m1,m2,ixo,iyo)
        m1=1
        call mouse(m1,m2,m3,m4)
        m2=0
        m1=3
c  loop on n added to fix Windows95 bug
   11   n=0
   10   call mouse(m1,m2,ix,iy)
        n=n+1
        if(n.lt.10000) go to 10
        if(m2.eq.0) go to 10
        if(n.eq.10000) go to 11
        iyo=iy
        ixo=ix
c        iy=nvert+tspace-iy
        return
      end
c******************************************************************************
      subroutine tgin2(m1,m2,m3,m4)
      call mouse(m1,m2,m3,m4)
      return
      end
c******************************************************************************
	subroutine tdraw(iop)
c  drawing mode:  iop= 0 pixel on, 1 off, 2 toggled
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
        iop1=iop+1
        go to (10,20,30),iop1
   10   icolor=0
        icsw=1
        go to 40
   20   icolor=-7
        icsw=-1
        go to 40
   30   if(icolor.eq.0) then
        icolor=-7
        icsw=-1
        else
        icolor=0
        icsw=1
        endif
   40   return
	end
c******************************************************************************
c      subroutine setctl(ixo,iyo)
c        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /plot/ dxp(2),xpix,ypix,xscale,yscale
c      common /data/ n,xmin,xmax,amin,amax,zmin,zmax,x(512),z(512),
c     &a(512)
c      character label(7)*8, prompt*79
c      data nmenu/7/,label/' window ',' vertex ','  reset ','  clear ',
c     &                    '         ','         ',' return '/
c      icolor=600
c      call menu(nmenu,label)
cc      ixo=600
cc      iyo=300
c  100 call menudvr(x,y,ifunc,ixo,iyo)
c      go to (101,102,103,104,105,106,107) ifunc
c  105 continue
c  106 continue
c      go to 100
c  107 call tdraw(1)
c      call menu(nmenu,label)
c      call tdraw(0)
cc      prompt=' '
cc      call tpromp(0,0,prompt)
cc      call tpromp(0,2,0,prompt)
c      return
c      end
c******************************************************************************
      subroutine setmodel(ixo,iyo)
c       common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c     common /plot/ dxp(2),dyp(2),xpix,ypix,xscale,yscale
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      common /pscale/ xmin,xmax,zmin,zmax
      common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /pmagdata/ pmmin,pmmax,pmminsav,pmmaxsav,pminc,pm(512)
      common /pgrvdata/ pgmin,pgmax,pgminsav,pgmaxsav,pginc,pg(512)
      common /rmagdata/ rmin,rmax,rminsav,rmaxsav,ra(512)
      common /contdata/ nco,comin,comax,cominsav,comaxsav,coinc,zup,
     1co(512)
      common /angles/ pmaginc,pmagdec,continc,contdec
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      common /work1/xc(512)
      common /calc1/icalc,icalcd,gc(512),hc(512)
      common /save/ gminsav,gmaxsav,aminsav,amaxsav
      complex xc,cmplx
      character label(7)*8, prompt*79, gfile*80, ifile*80, fmt*24
      character sfile*80
      dimension mpen(7)
      data nmenu/7/,label/' depth  ',' editmod',' upcont ',
     &         ' itgrav ','  swap  ','  clear ',' return  '/
      data mpen/0,0,200,0,0,0,0/
      icolor=600
      call menu(nmenu,label,mpen)
c      ixo=600
c      iyo=300
  100 call menudvr(xdum,y,ifunc,ixo,iyo)
      go to (100,101,103,104,102,105,106,107) ifunc+1
  101 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
c      call source(ixo,iyo,isor)
c      if(isor.eq.0) go to 111
      call setdep(ixo,iyo)
      icolor=600
  111 call menu(nmenu,label,mpen)
      go to 100
  102 if(igrv.eq.0.and.ipgrv.eq.0) go to 100
      call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
  121 call source(ixo,iyo,isor)
      if(isor.eq.0.or.isor.eq.7) go to 111
      if(isor.ne.2.and.isor.ne.5.and.isor.ne.6) go to 121
      if(isor.eq.6) then
        if(icont.ne.2.and.icont.ne.5) go to 121
      endif
      prompt='enter number of iterations and epsilon [5 .01]: '
      call tpromp(0,0,prompt)
      nmax=5
      eps=.01
c      read(5,*,err=106) nmax,eps
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        leni=index(prompt(1:leng),' ')-1
        lenm=leng-index(prompt(1:leng),'.')
        if(lenm.eq.leng) lenm=0
        lenf=leng-leni
        write(fmt,'(a,i3,a,i3,a,i3,a)') '(i',leni,',f',lenf,'.',lenm,')'
        read(prompt(1:leng),fmt,err=106) nmax,eps
      endif
  163 prompt=' '
      call tpromp(0,0,prompt)
      if(isor.eq.2) then
        sfact=1.
        soff=0.
        call itgrv(nmax,eps,soff,sfact,ng,g,gx,gz)
        prompt=' '
      else if(isor.eq.5) then
        if(igrv.gt.0) then
          sfact=(pgmax-pgmin)/(gmax-gmin)
          soff=pgmin-gmin*sfact
        else
          sfact=(pgmax-pgmin)/(amax-amin)
          soff=pgmin-amin*sfact
        endif
        call itgrv(nmax,eps,soff,sfact,n,pg,x,z)
        prompt=' '
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
      else if(isor.eq.6) then
        if(igrv.gt.0) then
          sfact=(comax-comin)/(gmax-gmin)
          soff=comin-gmin*sfact
        else
          sfact=(comax-comin)/(amax-amin)
          soff=comin-amin*sfact
        endif
        if(icont.eq.2) then
          do i=1,nco
          gz(i)=gz(i)-zup
          end do
          call itgrv(nmax,eps,soff,sfact,nco,co,gx,gz)
          do i=1,nco
          gz(i)=gz(i)+zup
          end do
        else if(icont.eq.5) then
          do i=1,nco
          z(i)=z(i)-zup
          end do
          call itgrv(nmax,eps,soff,sfact,nco,co,x,z)
          do i=1,nco
          z(i)=z(i)+zup
          end do
        endif
        prompt=' '
      endif
      call tpromp(0,0,prompt)
      go to 111
  103 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call setedt(ixo,iyo)
      icolor=600
      call menu(nmenu,label,mpen)
      go to 100
  104 continue
      call tdraw(1)
      if(icont.gt.0) then
        call dplot(2,xmin,xmax,comin,comax)
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          call curv(x,co,n,-7)
        else
          call curv(gx,co,ng,-7)
        endif
        call dplot(1,xmin,xmax,zmin,zmax)
        call rescale
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          call curv2(x,z,zup,n,-7)
        else
          call curv2(gx,gz,zup,ng,-7)
        endif
        icont=0
      endif
      call menu(nmenu,label,mpen)
      call tdraw(0)
  141 call source(ixo,iyo,isor)
      if(isor.eq.0.or.isor.ge.6) go to 111
c     if(isor.ne.2.and.isor.ne.5) go to 121
      if(isor.eq.1) then
        continc=ef(2)
        contdec=ef(3)
      else if(isor.eq.3) then
        continc=90.
        contdec=0.
      else if(isor.eq.4) then
        continc=pmaginc
        contdec=pmagdec
      else
        continc=1.e38
        contdec=0.
      endif
      icont=isor
      prompt='enter continuation distance (positive up): '
      call tpromp(0,0,prompt)
      read(5,*,err=104) zup
      pcd=1.e38
      if(zup.le.0.) then
        prompt='enter lowpass filter parameters startf (0 to 3.14) and d
     1ecay (>0): '
        call tpromp(0,0,prompt)
        read(5,*,err=104) pcd,alphad
      endif
  143 prompt=' '
      call tpromp(0,0,prompt)
      go to (151,152,153,154,155) isor
  151   call setfft(0,npad,n,a)
        coinc=50.
        nco=n
        go to 156
  152   call setfft(0,npad,ng,g)
        coinc=2.
        nco=ng
        go to 156
  153   call setfft(0,npad,n,rtp)
        coinc=50.
        nco=n
        go to 156
  154   call setfft(0,npad,ng,pm)
        coinc=pminc
        nco=ng
        go to 156
  155   call setfft(0,npad,n,pg)
        coinc=pginc
        nco=n
  156   pi=3.1415927
c      if(nco.eq.n) then
      if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
        up=zup/abs(x(2)-x(1))
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
      else
        up=zup/abs(gx(2)-gx(1))
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
      endif
      lim=npad/2+1
      fil=1.0
      do 157 i=1,lim
      j=i-1
c      p=j*2.*pi/(npad-1)
      p=j*2.*pi/float(npad)
c <<<<<<<<  low-pass filter & upward continuation   >>>>>>>>>>>>>
      if(p.gt.pcd) fil=exp(-alphad*(p-pcd)**2)
      if(fil.lt.0.0001)fil=0.0001
       fac=exp(-up*p)
c      cx(i)=cx(i)*fac*fil
       xc(i)=xc(i)*fac*fil
       xr=real(xc(i))
       xi=-1.0*aimag(xc(i))
       if(i.gt.1) xc((npad+2)-i)=cmplx(xr,xi)
  157  continue
       call fork(npad,xc,-1.0)
       comin=1.e+38
       comax=-1.e+38
       cominsav=1.e+38
       comaxsav=-1.e+38
c       if(nco.eq.n) then
       if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
         do 158 i=1,n
         co(i)=real(xc(i))
         if(co(i).gt.comaxsav) comaxsav=co(i)
         if(co(i).lt.cominsav) cominsav=co(i)
         if(x(i).lt.xmin) go to 158
         if(x(i).gt.xmax) go to 158
         if(co(i).gt.comax) comax=co(i)
         if(co(i).lt.comin) comin=co(i)
  158    continue
       else
         do 159 i=1,ng
         co(i)=real(xc(i))
         if(co(i).gt.comaxsav) comaxsav=co(i)
         if(co(i).lt.cominsav) cominsav=co(i)
         if(gx(i).lt.xmin) go to 159
         if(gx(i).gt.xmax) go to 159
         if(co(i).gt.comax) comax=co(i)
         if(co(i).lt.comin) comin=co(i)
  159    continue
       endif
       call dplot(2,xmin,xmax,comin,comax)
c       if(nco.eq.n) then
       if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
         call curv(x,co,n,200)
       else
         call curv(gx,co,ng,200)
       endif
      call dplot(1,xmin,xmax,zmin,zmax)
      call rescale
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          call curv2(x,z,zup,n,200)
        else
          call curv2(gx,gz,zup,ng,200)
        endif
c      icont=1
c      icont=isor
      call tpromp(0,0,prompt)
      go to 111
  105 if(icalc.gt.0) then
        prompt='click on "swap" again to swap observed and calculated pr
     1ofiles'
        call tpromp(0,0,prompt)
      else
        prompt='no calculated profiles to swap with observed profiles'
        call tpromp(0,0,prompt)
        go to 100
      endif
      call menudvr(xdum,y,ifunc,ixo,iyo)
      prompt=' '
      call tpromp(0,0,prompt)
      go to (100,100,100,100,100,205,100,100) ifunc+1
  205 if(igrv.gt.0) then
        gmax=-1.e38
        gmin=1.e38
        gmaxsav=-1.e38
        gminsav=1.e38
        do 210 i=1,ng
        temp=gc(i)
        gc(i)=g(i)
        g(i)=temp
        if(temp.gt.gmaxsav) gmaxsav=temp
        if(temp.lt.gminsav) gminsav=temp
        if(gx(i).lt.xmin) go to 210
        if(gx(i).gt.xmax) go to 210
        if(temp.gt.gmax) gmax=temp
        if(temp.lt.gmin) gmin=temp
  210   continue
        if(gmin.eq.gmax) then
          gmin=gmin-10.
          gmax=gmax+10.
        endif
      endif
      if(imag.gt.0) then
        amax=-1.e38
        amin=1.e38
        amaxsav=-1.e38
        aminsav=1.e38
        do 215 i=1,n
        temp=hc(i)
        hc(i)=a(i)
        a(i)=temp
        if(temp.gt.amaxsav) amaxsav=temp
        if(temp.lt.aminsav) aminsav=temp
        if(x(i).lt.xmin) go to 215
        if(x(i).gt.xmax) go to 215
        if(temp.gt.amax) amax=temp
        if(temp.lt.amin) amin=temp
  215   continue
        if(amin.eq.amax) then
          amin=amin-10.
          amax=amax+10.
        endif
      endif
      irmag=0
      ipmag=0
      ipgrv=0
      icont=0
      icalc=3-icalc
  106 call clear
      call menu(nmenu,label,mpen)
      go to 100
  107 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
c      prompt=' '
c      call tpromp(0,0,prompt)
c      call tpromp(0,2,0,prompt)
      return
      end
c******************************************************************************
      subroutine setpro(ixo,iyo)
c       common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c     common /plot/ dxp(2),dyp(2),xpix,ypix,xscale,yscale
      common /pscale/ xmin,xmax,zmin,zmax
      common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /pmagdata/ pmmin,pmmax,pmminsav,pmmaxsav,pminc,pm(512)
      common /pgrvdata/ pgmin,pgmax,pgminsav,pgmaxsav,pginc,pg(512)
      common /rmagdata/ rmin,rmax,rminsav,rmaxsav,ra(512)
      common /contdata/ nco,comin,comax,cominsav,comaxsav,coinc,zup,
     1co(512)
      common /angles/ pmaginc,pmagdec,continc,contdec
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      character label(7)*8, prompt*79, gfile*80, sfile*80
      dimension mpen(7)
      data nmenu/7/,label/'terrace ',' redpol ','psudomag',
     &         'psudogrv','clrprof ',' clear  ',' return  '/
      data mpen/0,500,600,100,0,0,0/
      icolor=600
      call menu(nmenu,label,mpen)
c      ixo=600
c      iyo=300
  100 call menudvr(xdum,y,ifunc,ixo,iyo)
      go to (100,101,102,103,104,105,106,107) ifunc+1
  101 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call source(ixo,iyo,isor)
      call menu(nmenu,label,mpen)
      if(isor.eq.0) go to 100
c      prompt='enter number of iterations: '
c      call tpromp(0,0,prompt)
c      read(5,*) iter
      prompt=' '
      go to (121,122,123,124,125,126) isor
  121   call dplot(2,xmin,xmax,amin,amax)
        call terrace1(20,400,n,a,x)
        call terrace1(21,400,n,a,x)
        go to 127
  122   call dplot(2,xmin,xmax,gmin,gmax)
        call terrace1(20,300,ng,g,gx)
        call terrace1(21,300,ng,g,gx)
        go to 127
  123   call dplot(2,xmin,xmax,rmin,rmax)
        call terrace1(20,500,n,ra,x)
        call terrace1(21,500,n,ra,x)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 127
  124   continue
        call dplot(2,xmin,xmax,pmmin,pmmax)
        call terrace1(20,600,ng,pm,gx)
        call terrace1(21,600,ng,pm,gx)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 127
  125   continue
        call dplot(2,xmin,xmax,pgmin,pgmax)
        call terrace1(20,100,n,pg,x)
        call terrace1(21,100,n,pg,x)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 127
  126   continue
        call dplot(2,xmin,xmax,comin,comax)
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          call terrace1(20,200,nco,co,x)
          call terrace1(21,200,nco,co,x)
          if(imag.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled magnetic profile'
        else
          call terrace1(20,200,nco,co,gx)
          call terrace1(21,200,nco,co,gx)
          if(igrv.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled gravity profile'
        endif
  127 call dplot(1,xmin,xmax,zmin,zmax)
      call rescale
      call tpromp(0,0,prompt)
      go to 100
  102 if(imag.eq.0) go to 100
        if(ef(1).eq.0.0) then
          prompt='enter current inc, dec, strike of profile: '
          call tpromp(0,0,prompt)
          read(5,*,err=100) ef(2), ef(3), azmuth
        endif
        call redpol(ef(2),ef(3),azmuth,0)
        call dplot(1,xmin,xmax,zmin,zmax)
        call rescale
        prompt=' '
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        call tpromp(0,0,prompt)
      go to 100
  103 if(igrv.eq.0) go to 100
        if(ef(1).eq.0.0) then
          prompt='enter desired inc, dec, and current profile strike: '
          call tpromp(0,0,prompt)
          read(5,*,err=100) pmaginc, pmagdec, azmuth
        else
          prompt='enter desired inc, dec: '
          call tpromp(0,0,prompt)
          read(5,*,err=100) pmaginc, pmagdec
        endif
        call redpol(pmaginc,pmagdec,azmuth,1)
        call dplot(1,xmin,xmax,zmin,zmax)
        call rescale
        prompt=' '
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        call tpromp(0,0,prompt)
      go to 100
  104 go to 102
  105 irmag=0
      ipmag=0
      ipgrv=0
      icont=0
  106 call clear
      call menu(nmenu,label,mpen)
      go to 100
  107 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
c      prompt=' '
c      call tpromp(0,0,prompt)
c      call tpromp(0,2,0,prompt)
      return
      end
c******************************************************************************
      subroutine source(ixo,iyo,ifunc)
c       common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c     common /plot/ dxp(2),dyp(2),xpix,ypix,xscale,yscale
c      common /magdata/ n,xmin,xmax,amin,amax,zmin,zmax,
c     &x(512),z(512),a(512)
c      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      character label(7)*8, prompt*79, gfile*80, sfile*80
      dimension mpen(7)
      data nmenu/7/,label/'  mag   ','  grav  ',' rtpmag ',
     &         '  pmag  ',' pgrav  ',' upcont ','  none   '/
      data mpen/400,300,500,600,100,200,0/
      icolor=600
      call menu(nmenu,label,mpen)
c      ixo=600
c      iyo=300
      prompt='select source from menu'
      call tpromp(0,0,prompt)
  100 call menudvr(xdum,y,ifunc,ixo,iyo)
      go to (100,101,102,103,104,105,106,107) ifunc+1
  101 if(imag.eq.0) go to 100
      go to 200
  102 continue
      if(igrv.eq.0) go to 100
      go to 200
  103 continue
      if(irmag.ne.1) go to 100
      go to 200
  104 continue
      if(ipmag.ne.1) go to 100
      go to 200
  105 if(ipgrv.ne.1) go to 100
      go to 200
  106 if(icont.eq.0) go to 100
      go to 200
  107 ifunc=0
  200 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      prompt=' '
      call tpromp(0,0,prompt)
      return
      end
c******************************************************************************
      subroutine setdep(ixo,iyo)
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
      common /plot/ dxp(2),dyp(2),xpix,ypix,xscale,yscale
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      common /pscale/ xmin,xmax,zmin,zmax
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /pmagdata/ pmmin,pmmax,pmminsav,pmmaxsav,pminc,pm(512)
      common /pgrvdata/ pgmin,pgmax,pgminsav,pgmaxsav,pginc,pg(512)
      common /rmagdata/ rmin,rmax,rminsav,rmaxsav,ra(512)
      common /contdata/ nco,comin,comax,cominsav,comaxsav,coinc,zup,
     1co(512)
      common /angles/ pmaginc,pmagdec,continc,contdec
      common /work2/dummy(512),tmp(512)
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
        common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      character label(7)*8, prompt*79, ans*1, sym*1, gfile*80, sfile*80
      character ifile*80,fmt*24
      dimension mpen(7)
      data nmenu/7/,label/'  MSW   ','  euler ','  ansig ',
     &         ' autocor',' werner ','  clear ',' return  '/
      data mpen/0,0,0,0,0,0,0/
c        data d2r/1.745329e-2/
      icolor=600
      open(20,file='pdepth.dep',status='unknown',form='formatted')
      call menu(nmenu,label,mpen)
c      ixo=600
c      iyo=300
  100 call menudvr(xdum,y,ifunc,ixo,iyo)
      go to (100,105,102,103,104,101,106,107) ifunc+1
  101 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call source(ixo,iyo,iflg)
      call menu(nmenu,label,mpen)
      if(iflg.eq.0) go to 100
      prompt='enter degree of interference polynomial (0 to 5): [2] '
      call tpromp(0,0,prompt)
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        write(fmt,'(a,i1,a)') '(i',leng,')'
        read(prompt(1:leng),fmt,err=101) ndeg
      else
        ndeg=2
      endif
c      read(5,*,err=101) ndeg
      if(ndeg.lt.0.or.ndeg.gt.5) go to 101
      prompt='enter clustering radius as a percent of depth: [5] '
      call tpromp(0,0,prompt)
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=index(prompt(1:leng),'.')
        if(lenm.eq.leng) lenm=0
        write(fmt,'(a,i3,a,i3,a)') '(f',leng,'.',lenm,')'
        read(prompt(1:leng),fmt,err=101) eps
      else
        eps=5.
      endif
      prompt='plot individual solutions? [y]: '
      call tpromp(0,0,prompt)
      read(5,200,err=101) ans
      ipltall=1
      if(ans.eq.' ') ans='y'
      if(ans.eq.'n'.or.ans.eq.'N') ipltall=0
      prompt=' '
      call tpromp(0,0,prompt)
      go to (111,112,113,114,115,116) iflg
  111   ipen=400
      write(20,200) 'werner depths - observed magnetics'
        call werner(ipltall,ndeg,eps,ipen,zmin,zmax,n,x,z,a)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 117
  112   ipen=300
      write(20,200) 'werner depths - observed gravity'
        call werner(ipltall,ndeg,eps,ipen,zmin,zmax,ng,gx,gz,g)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 117
  113   ipen=500
      write(20,200) 'werner depths - RTP magnetics'
        call werner(ipltall,ndeg,eps,ipen,zmin,zmax,n,x,z,ra)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 117
  114   ipen=600
      write(20,200) 'werner depths - pseudomagnetics'
        call werner(ipltall,ndeg,eps,ipen,zmin,zmax,ng,gx,gz,pm)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 117
  115   ipen=100
      write(20,200) 'werner depths - pseudogravity'
        call werner(ipltall,ndeg,eps,ipen,zmin,zmax,n,x,z,pg)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 117
  116   ipen=200
      write(20,200) 'werner depths - continued field'
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          do i=1,nco
          z(i)=z(i)-zup
          end do
          call werner(ipltall,ndeg,eps,ipen,zmin,zmax,nco,x,z,co)
          do i=1,nco
          z(i)=z(i)+zup
          end do
          if(imag.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled magnetic profile'
        else
          do i=1,nco
          gz(i)=gz(i)-zup
          end do
          call werner(ipltall,ndeg,eps,ipen,zmin,zmax,nco,gx,gz,co)
          do i=1,nco
          gz(i)=gz(i)+zup
          end do
          if(igrv.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled gravity profile'
        endif
  117 call tpromp(0,0,prompt)
      go to 100
  102 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call source(ixo,iyo,iflg)
      call menu(nmenu,label,mpen)
      if(iflg.eq.0) go to 100
      prompt='do you want to specify the euler dimension?: '
      call tpromp(0,0,prompt)
      read(5,200,err=102) ans
  200 format(a)
      mdim=4
      if(ans.eq.'y'.or.ans.eq.'Y') then
        mdim=3
        prompt='enter euler dimension (0 to 3): '
        call tpromp(0,0,prompt)
        read(5,*,err=102) edim
        if(edim.lt.0.0.or.edim.gt.3.0) go to 102
        if(edim.eq.0.0) edim=.0001
      endif
      prompt='enter window length (5 to 20): '
      call tpromp(0,0,prompt)
      read(5,*,err=102) iwin
      if(iwin.lt.5.or.iwin.gt.20) go to 102
c      prompt='enter maximum mean square error: '
      prompt='enter tolerance (0 to plot all solutions): '
      call tpromp(0,0,prompt)
      read(5,*,err=102) err
      prompt=' '
      call tpromp(0,0,prompt)
      go to (121,122,123,124,125,126) iflg
  121   ipen=400
      write(20,200) 'euler depths - observed magnetics'
        call euler(edim,mdim,iwin,err,ipen,n,x,z,a)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 127
  122   ipen=300
      write(20,200) 'euler depths - observed gravity'
        call euler(edim,mdim,iwin,err,ipen,ng,gx,gz,g)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 127
  123   ipen=500
      write(20,200) 'euler depths - RTP magnetics'
        call euler(edim,mdim,iwin,err,ipen,n,x,z,ra)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 127
  124   ipen=600
      write(20,200) 'euler depths - pseudomagnetics'
        call euler(edim,mdim,iwin,err,ipen,ng,gx,gz,pm)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 127
  125   ipen=100
      write(20,200) 'euler depths - pseudogravity'
        call euler(edim,mdim,iwin,err,ipen,n,x,z,pg)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 127
  126   ipen=200
      write(20,200) 'euler depths - continued field'
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          do i=1,nco
          z(i)=z(i)-zup
          end do
          call euler(edim,mdim,iwin,err,ipen,nco,x,z,co)
          do i=1,nco
          z(i)=z(i)+zup
          end do
          if(imag.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled magnetic profile'
        else
          do i=1,nco
          gz(i)=gz(i)-zup
          end do
          call euler(edim,mdim,iwin,err,ipen,nco,gx,gz,co)
          do i=1,nco
          gz(i)=gz(i)+zup
          end do
          if(igrv.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled gravity profile'
        endif
  127  call tpromp(0,0,prompt)
       go to 100
c      if(iflg.eq.1) then
c        prompt='enter einc, edec, strike (in degrees): '
c        call tpromp(0,0,prompt)
c        read(5,*,err=102) einc, edec, strike
c        prompt=' '
c        call tpromp(0,0,prompt)
c        call euler(einc,edec,strike,zmin,zmax,n,x,z,a)
c      else
c        if(iflg.eq.2) call euler(90.,0.,0.,zmin,zmax,ng,gx,gz,g)
c      endif
c      go to 100
  103 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call source(ixo,iyo,iflg)
      call menu(nmenu,label,mpen)
      if(iflg.eq.0) go to 100
c  103 call tdraw(1)
c      call menu(nmenu,label)
c      call tdraw(0)
c      call source(ixo,iyo,isor)
c      call menu(nmenu,label)
c      if(isor.eq.0) go to 100
c      prompt='include derivative solutions?: '
c      call tpromp(0,0,prompt)
c      iopt=noyes()
      iopt=0
      prompt='plot curves?: [y] '
      call tpromp(0,0,prompt)
      read(5,200,err=101) ans
      ipen=1
      if(ans.eq.' ') ans='y'
      if(ans.eq.'n'.or.ans.eq.'N') ipen=0
c      ipen=noyes()
      ipen=2*ipen
      iopt=iopt+ipen
      prompt=' '
      call tpromp(0,0,prompt)
      go to (131,132,133,134,135,136) iflg
  131 continue
      if(ef(1).eq.0.0) then
  138   prompt='enter field magnitude, inc, dec, and profile azimuth: '
        call tpromp(0,0,prompt)
        read(5,*,end=138,err=138) ef,azmuth
c        ei=ef(2)*d2r
c        ed=(ef(3)-azmuth)*d2r
c        eh=cos(ei)
c        ev(1)=cos(ed)*eh
c        ev(2)=sin(ed)*eh
c        ev(3)=sin(ei)
c        do 139 i=1,3
c  139   tdir(i)=ev(i)
      endif
      write(20,200) 'analytic signal depths - observed magnetics'
      write(20,202)
      write(20,203)
      ipen=400
      sym='D'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1ef(2),ef(3),n,a,x,z)
      ipen=1200
      delx=abs(x(n)-x(1))/float(n-1)
      call spline(a,tmp,n,delx)
      write(20,204)
      sym='C'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1ef(2),ef(3),n,tmp,x,z)
      if(imag.eq.2) prompt='Warning: results unreliable due to unevenly
     1sampled magnetic profile'
      go to 137
  132 ipen=300
      write(20,200) 'analytic signal depths - observed gravity'
      write(20,202)
      write(20,203)
      sym='P'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     11.e38,0.,ng,g,gx,gz)
      ipen=1100
      delx=abs(gx(ng)-gx(1))/float(ng-1)
      call spline(g,tmp,ng,delx)
      write(20,204)
      sym='D'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     11.e38,0.,ng,tmp,gx,gz)
      if(igrv.eq.2) prompt='Warning: results unreliable due to unevenly
     1sampled gravity profile'
      go to 137
  133 ipen=500
      write(20,200) 'analytic signal depths - RTP magnetics'
      write(20,202)
      write(20,203)
      sym='D'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,90.,0.,n,ra,x,z)
      ipen=1300
      delx=abs(x(n)-x(1))/float(n-1)
      call spline(ra,tmp,n,delx)
      write(20,204)
      sym='C'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,90.,0.,n,tmp,x,z)
      if(imag.eq.2) prompt='Warning: results unreliable due to unevenly
     1sampled magnetic profile'
      go to 137
  134 ipen=600
      write(20,200) 'analytic signal depths - pseudomagnetics'
      write(20,202)
      write(20,203)
      sym='D'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1pmaginc,pmagdec,ng,pm,gx,gz)
      ipen=1400
      delx=abs(gx(ng)-gx(1))/float(ng-1)
      call spline(pm,tmp,ng,delx)
      write(20,204)
      sym='C'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1pmaginc,pmagdec,ng,tmp,gx,gz)
      if(igrv.eq.2) prompt='Warning: results unreliable due to unevenly
     1sampled gravity profile'
      go to 137
  135 ipen=100
      write(20,200) 'analytic signal depths - pseudogravity'
      write(20,202)
      write(20,203)
      sym='P'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,1.e38,0.,n,pg,x,z)
      ipen=900
      delx=abs(x(n)-x(1))/float(n-1)
      call spline(pg,tmp,n,delx)
      write(20,204)
      sym='D'
      call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,1.e38,0.,n,tmp,x,z)
      if(imag.eq.2) prompt='Warning: results unreliable due to unevenly
     1sammpled magnetic profile'
      go to 137
  136 ipen=200
      write(20,200) 'analytic signal depths - continued field'
      write(20,202)
      write(20,203)
c      if(nco.eq.n) then
      if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
        do i=1,nco
        z(i)=z(i)-zup
        end do
        sym='D'
        call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1  continc,contdec,nco,co,x,z)
        ipen=1000
      delx=abs(x(n)-x(1))/float(n-1)
        call spline(co,tmp,nco,delx)
      write(20,204)
        sym='C'
        call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1  continc,contdec,nco,tmp,x,z)
        do i=1,nco
        z(i)=z(i)+zup
        end do
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
      else
        do i=1,nco
        gz(i)=gz(i)-zup
        end do
        sym='D'
        call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1  continc,contdec,nco,co,gx,gz)
        ipen=1000
      delx=abs(gx(ng)-gx(1))/float(ng-1)
        call spline(co,tmp,nco,delx)
      write(20,204)
        sym='C'
        call hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1  continc,contdec,nco,tmp,gx,gz)
        do i=1,nco
        gz(i)=gz(i)+zup
        end do
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
      endif
  137 call dplot(1,xmin,xmax,zmin,zmax)
      call rescale
      call tpromp(0,0,prompt)
      go to 100
  104 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call source(ixo,iyo,iflg)
      call menu(nmenu,label,mpen)
      if(iflg.eq.0) go to 100
      prompt='enter window length range (two values between 7 and 29): '
      call tpromp(0,0,prompt)
      read(5,*,err=104) ndeg1,ndeg2
      if(ndeg1.gt.ndeg2) go to 104
      if(ndeg1.lt.7.or.ndeg2.gt.29) go to 104
      prompt='enter required convergence for primary depths: '
      call tpromp(0,0,prompt)
      read(5,*,err=104) eps1
      prompt='enter required convergence for secondary depths: '
      call tpromp(0,0,prompt)
      read(5,*,err=104) eps2
      prompt=' '
      call tpromp(0,0,prompt)
      do 147 ndeg=ndeg1,ndeg2,2
      go to (141,142,143,144,145,146) iflg
  141   ipen=400
        ipen2=1200
      write(20,200) 'autocorrelation depths - observed magnetics'
      write(20,210) ndeg
      write(20,205)
        call adept(ndeg,ipen2,ipen,n,x,z,a,eps1,eps2)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 147
  142   ipen=300
        ipen2=1100
      write(20,200) 'autocorrelation depths - observed gravity'
      write(20,210) ndeg
      write(20,205)
        call adept(ndeg,ipen2,ipen,ng,gx,gz,g,eps1,eps2)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 147
  143   ipen=500
        ipen2=1300
      write(20,200) 'autocorrelation depths - RTP magnetics'
      write(20,210) ndeg
      write(20,205)
        call adept(ndeg,ipen2,ipen,n,x,z,ra,eps1,eps2)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 147
  144   ipen=600
        ipen2=1400
      write(20,200) 'autocorrelation depths - pseudomagnetics'
      write(20,210) ndeg
      write(20,205)
        call adept(ndeg,ipen2,ipen,ng,gx,gz,pm,eps1,eps2)
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 147
  145   ipen=100
        ipen2=900
      write(20,200) 'autocorrelation depths - pseudogravity'
      write(20,210) ndeg
      write(20,205)
        call adept(ndeg,ipen2,ipen,n,x,z,pg,eps1,eps2)
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 147
  146   ipen=200
        ipen2=1000
      write(20,200) 'autocorrelation depths - continued field'
      write(20,210) ndeg
      write(20,205)
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          do i=1,nco
          z(i)=z(i)-zup
          end do
          call adept(ndeg,ipen2,ipen,nco,x,z,co,eps1,eps2)
          do i=1,nco
          z(i)=z(i)+zup
          end do
          if(imag.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled magnetic profile'
        else
          do i=1,nco
          gz(i)=gz(i)-zup
          end do
          call adept(ndeg,ipen2,ipen,nco,gx,gz,co,eps2,eps2)
          do i=1,nco
          gz(i)=gz(i)+zup
          end do
          if(igrv.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled gravity profile'
        endif
  147 continue
      call tpromp(0,0,prompt)
      go to 100
  105 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      call source(ixo,iyo,iflg)
      call menu(nmenu,label,mpen)
      if(iflg.eq.0) go to 100
c  150 prompt='click on left limit'
c      call tpromp(0,0,prompt)
c      call menudvr(xmin1,y,ifunc,ixo,iyo)
c      if(ifunc.ne.0) go to 100
c      prompt='click on right limit'
c      call tpromp(0,0,prompt)
c      call menudvr(xmax1,y,ifunc,ixo,iyo)
c      if(ifunc.ne.0) go to 100
c      if(xmax1.le.xmin1) go to 150
C
      prompt='enter window length code (0=22, 1=44, 2=88, 3=176) [0] '
      call tpromp(0,0,prompt)
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        write(fmt,'(a,i3,a)') '(i',leng,')'
        read(prompt(1:leng),fmt,err=105) iwin
      else
        iwin=0
      endif
      iwin=22*2**iwin
c      print*, iwin
c      read(5,*,err=105) iwin,eps
      if(iwin.eq.0) go to 157
      if(iwin.lt.1) iwin=1
c      if(iwin.gt.21) iwin=21
      iwin=iwin-1
      prompt='enter clustering radius as a percent of depth: [5] '
      call tpromp(0,0,prompt)
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=index(prompt(1:leng),'.')
        if(lenm.eq.leng) lenm=0
        write(fmt,'(a,i3,a,i3,a)') '(f',leng,'.',lenm,')'
        read(prompt(1:leng),fmt,err=101) eps
      else
        eps=5.
      endif
      prompt='enter minimum number of solutions in cluster: [2] '
      call tpromp(0,0,prompt)
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        write(fmt,'(a,i3,a)') '(i',leng,')'
        read(prompt(1:leng),fmt,err=101) minclus
      else
        minclus=2
      endif
      prompt='enter min, max number of edges (two values <11): [3 3] '
      call tpromp(0,0,prompt)
c      read(5,*,err=105) ndeg1,ndeg2
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        leni=index(prompt(1:leng),' ')-1
        lenf=leng-leni
        write(fmt,'(a,i3,a,i3,a,i3,a)') '(i',leni,',i',lenf,')'
        read(prompt(1:leng),fmt,err=105) ndeg1,ndeg2
      else
        ndeg1=3
        ndeg2=3
      endif
      prompt='sheet or contact solutions (s or c)?: [c] '
      call tpromp(0,0,prompt)
      read(5,200,err=105) ans
      isheet=0
c      if(ans.eq.' ') ans='c'
      if(ans.eq.'s'.or.ans.eq.'S') isheet=1
      ipltall=1
      if(ndeg1.eq.ndeg2) then
        prompt='plot individual solutions?: [n] '
        call tpromp(0,0,prompt)
        read(5,200,err=105) ans
        ipltall=0
        if(ans.eq.'y'.or.ans.eq.'Y') ipltall=1
      endif
c
c      prompt='enter separation factor: [1] '
c      call tpromp(0,0,prompt)
c      read (*,'(a)') prompt
c      leng=len_trim(prompt)
c      if(leng.ne.0) then
c        lenm=index(prompt(1:leng),'.')
c        if(lenm.eq.leng) lenm=0
c        write(fmt,'(a,i3,a,i3,a)') '(f',leng,'.',lenm,')'
c        read(prompt(1:leng),fmt,err=101) cc
c      else
c        cc=1.
c      endif
c
      prompt='computing solutions...'
      call tpromp(0,0,prompt)
      go to (151,152,153,154,155,156) iflg
  151 if(ef(1).eq.0.0) then
  158   prompt='enter field magnitude, inc, dec, and profile azimuth: '
        call tpromp(0,0,prompt)
        read(5,*,end=158,err=158) ef,azmuth
        prompt='computing solutions...'
        call tpromp(0,0,prompt)
      endif
      write(20,200) 'multi-source werner depths - observed magnetics'
      write(20,211) iwin,eps,minclus,ndeg1,ndeg2
      if(isheet.eq.1) then
      write(20,203)
      write(20,206)
        ipen=400
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,n,a,x,z,eps,
     1  ef(2),ef(3),minclus)
      else
      write(20,204)
      write(20,206)
        ipen=1200
      delx=abs(x(n)-x(1))/float(n-1)
        call spline(a,tmp,n,delx)
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,n,tmp,x,z,eps,
     1  ef(2),ef(3),minclus)
      endif
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 157
  152 write(20,200) 'multi-source werner depths - observed gravity'
      write(20,211) iwin,eps,minclus,ndeg1,ndeg2
      if(isheet.eq.1) then
      write(20,203)
      write(20,206)
        ipen=300
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,ng,g,gx,gz,eps,
     1  1.e38,0.,minclus)
      else
      write(20,204)
      write(20,206)
        ipen=1100
      delx=abs(gx(ng)-gx(1))/float(ng-1)
        call spline(g,tmp,ng,delx)
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,ng,tmp,gx,gz,eps,
     1  1.e38,0.,minclus)
      endif
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 157
  153 write(20,200) 'multi-source werner depths - RTP magnetics'
      write(20,211) iwin,eps,minclus,ndeg1,ndeg2
      if(isheet.eq.1) then
      write(20,203)
      write(20,206)
        ipen=500
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,n,ra,x,z,eps,
     1  90.,0.,minclus)
      else
      write(20,204)
      write(20,206)
        ipen=1300
      delx=abs(x(n)-x(1))/float(n-1)
        call spline(ra,tmp,n,delx)
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,n,tmp,x,z,eps,
     1  90.,0.,minclus)
      endif
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 157
  154 write(20,200) 'multi-source werner depths - pseudomagnetics'
      write(20,211) iwin,eps,minclus,ndeg1,ndeg2
      if(isheet.eq.1) then
      write(20,203)
      write(20,206)
        ipen=600
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,ng,pm,gx,gz,eps,
     1  pmaginc,pmagdec,minclus)
      else
      write(20,204)
      write(20,206)
        ipen=1400
      delx=abs(gx(ng)-gx(1))/float(ng-1)
        call spline(pm,tmp,ng,delx)
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,ng,tmp,gx,gz,eps,
     1  pmaginc,pmagdec,minclus)
      endif
        if(igrv.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled gravity profile'
        go to 157
  155 write(20,200) 'multi-source werner depths - pseudogravity'
      write(20,211) iwin,eps,minclus,ndeg1,ndeg2
      if(isheet.eq.1) then
      write(20,203)
      write(20,206)
        ipen=100
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,n,pg,x,z,eps,
     1  1.e38,0.,minclus)
      else
      write(20,204)
      write(20,206)
        ipen=900
      delx=abs(x(n)-x(1))/float(n-1)
        call spline(pg,tmp,n,delx)
        call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,n,tmp,x,z,eps,
     1  1.e38,0.,minclus)
      endif
        if(imag.eq.2) prompt='Warning: results unreliable due to unevenl
     1y sampled magnetic profile'
        go to 157
  156 continue
      write(20,200) 'multi-source werner depths - continued field'
      write(20,211) iwin,eps,minclus,ndeg1,ndeg2
c      if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
        do i=1,nco
        z(i)=z(i)-zup
        end do
        if(isheet.eq.1) then
      write(20,203)
      write(20,206)
          ipen=200
          call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,nco,co,x,z,eps,
     1    continc,contdec,minclus)
        else
      write(20,204)
      write(20,206)
      delx=abs(x(n)-x(1))/float(n-1)
          call spline(co,tmp,nco,delx)
          ipen=1000
          call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,nco,tmp,x,z,eps,
     1    continc,contdec,minclus)
        endif
        do i=1,nco
        z(i)=z(i)+zup
        end do
          if(imag.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled magnetic profile'
      else
        do i=1,nco
        gz(i)=gz(i)-zup
        end do
        if(isheet.eq.1) then
      write(20,203)
      write(20,206)
          ipen=200
          call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,nco,co,gx,gz,eps,
     1    continc,contdec,minclus)
        else
      write(20,204)
      write(20,206)
      delx=abs(gx(ng)-gx(1))/float(ng-1)
          call spline(co,tmp,nco,delx)
          ipen=1000
          call msw(ipltall,ierr,ipen,ndeg1,ndeg2,iwin,nco,tmp,gx,gz,eps,
     1    continc,contdec,minclus)
        endif
        do i=1,nco
        gz(i)=gz(i)+zup
        end do
          if(igrv.eq.2) prompt='Warning: results unreliable due to uneve
     1nly sampled gravity profile'
      endif
  157 continue
      if(ierr.eq.-1) then
        prompt='Error: inappropriate number of edges'
      else if(ierr.eq.1) then
        prompt='Error: more unknowns than data points'
      else if(ierr.eq.2) then
        prompt='Error: window too long'
      else
        prompt=' '
      endif
      call dplot(1,xmin,xmax,zmin,zmax)
      call rescale
      call tpromp(0,0,prompt)
      go to 100
  106 call clear
      call menu(nmenu,label,mpen)
      go to 100
  107 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      close(20)
c      prompt=' '
c      call tpromp(0,0,prompt)
c      call tpromp(0,2,0,prompt)
      return
  202 format('         x              z             dips / contrasts...'
     1)
  203 format('  sheet solutions')
  204 format('  contact solutions')
  205 format('         x              z             err')
  206 format('         x              z             dip           contra
     1st         #')
  210 format(i10,' point window')
  211 format('  window  = ',i3,'  radius = ',f5.1,'%   minimum# = ',i3,
     1'  #edges = ',2i3)
      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	print  6
6	format(' y or n:'$)
	ic=ic+1
	go to 1
9	print *,' count exceeded, answering no'
	noyes=0
	return
        end
      subroutine fork(lx,cx,signi)
c
c                      lx
c  cx(k) = sqrt(1/lx) sum (cx(j)*exp(2*pi*signi*i*(j-1)*(k-1)/lx))
c                     j=1
c
c                              for k=1,2,...,(lx=2**integer)
c  written by j.f.claerbout.
      complex cx(lx),carg,cexp,cw,ctemp
      j=1
      sc=sqrt(1./lx)
      do 5 i=1,lx
      if(i.gt.j)go to 2
      ctemp=cx(j)*sc
      cx(j)=cx(i)*sc
      cx(i)=ctemp
    2 m=lx/2
    3 if(j.le.m)go to 5
      j=j-m
      m=m/2
      if(m.ge.1)go to 3
    5 j=j+m
      l=1
    6 istep=2*l
      do 8 m=1,l
      carg=(0.,1.)*(3.14159265*signi*(m-1))/l
      cw=cexp(carg)
      do 8 i=m,lx,istep
      ctemp=cw*cx(i+l)
      cx(i+l)=cx(i)-ctemp
    8 cx(i)=cx(i)+ctemp
      l=istep
      if(l.lt.lx)go to 6
    9 return
      end
c ****************************************************************************
c      program itgravplt
c**********************************************************************
c
c  L. Cordell's two-dimensional iterative gravity modelling program
c  with graphic output (set for Tektronix 4114.
c        ncol = number of points in input data traverse
c        x(i) = field point, not necessarily equi-spaced, i=1,1,..ncol
c        g(i) = gravity at x(i)
c        t(i) = thickness, from dzi positive downwards
c        dzi  = reference horizon-- can be zero (prisms at surface)
c               (prism top for ridge, bottom for basin model),
c               dz1 >or = 0, in kilometers.
c        rho  = density contrast
c        nmax = number of iterations
c        eps  = terminates when rms error <or = eps
c        dx   = grid interval in kilometers
c**********************************************************************
      subroutine itgrv(nmax,eps,soff,sfact,ncol,gsav,x,z)
c      character*50 infile,outfile1,outfile2
c      character*5  ans,ridge,basin
      character*1  type,answer,sym,sym1
      character*79 prompt
      dimension gsav(ncol),x(ncol)
      dimension g(512),gn(512),t(512),th(512),refhoriz(2),z(512)
c      common /grvdata/ ncol,gsmin,gsmax,x(512),gz(512),gsav(512)
c      data dval/0.5e38/,
      data size/.08/
      sym='#'
      sym1='.'
c
c101   format (a5)
c102   format (a50)
c      print 100
c100   format (' enter input file name:')
c      read 102, infile
c      open (10,file=infile,form='formatted',status='old',iostat=ios)
cc        if(ios.ne.0) then
cc        open (10,file=infile,form='unformatted',status='old')
cc        endif
c      open (11,file='itgrvplt.out',form='formatted',status='unknown')
c      print 301
c301   format (' enter number of iterations and epsilon')
c      read *, nmax,eps
      eps=abs(eps)
      if(nmax.lt.11) go to 307
      nmax=11
      prompt='setting iterations to 11'
      call tpromp(0,0,prompt)
c      print 305, nmax
c305   format (' normally only a few iterations are needed and after',/,
c     & ' many iterations the solution may oscillate. do you really',/,
c     & ' want to do ',i2,' iterations?')
c      read 103, answer
c      if(answer.eq.'y') go to 307
c      print 306
c306   format (' enter number of iterations:')
c      read *, nmax
307   continue
c      print *,'enter iplotr (5=hp, 8=cga, 9=ega, 10=vga)'
c      read *,iplotr
c      print *,'enter vertical exaggeration factor'
c      read *, vert
3     kk=0
2     kk=kk+1
      icode=100*kk
      if(kk.gt.14) go to 3
300   prompt='enter density contrast: '
      call tpromp(0,0,prompt)
c      print 300
c300   format (' enter density contrast')
      read *, rho
      if(rho.eq.0.0) go to 999
      rho=rho*sfact
362   prompt='enter reference plane depth: '
      call tpromp(0,0,prompt)
c362   print 302
c302   format (' enter reference plane depth')
      read *, dz1
      prompt='enter mode (ridge or basin): '
      call tpromp(0,0,prompt)
c      print 303
c303   format (' enter mode (ridge or basin)')
      read 103, type
      if(type.eq.'B') type='b'
      if(type.eq.'R') type='r'
      if(type.eq.'b'.or.type.eq.'r') go to 360
c363   print 361
c361   format (' try again')
      go to 362
360   if(dz1.ge.0.0e0) go to 304
c      print 364
c364   format (' reference plane depth must be > or = 0.')
c      go to 363
       go to 362
304   iout=1
c      read (10) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
c      i=0
c10    i=i+1
c      read(10,*,end=20) x(i),y,g(i)
c      go to 10
c20    ncol=i-1
c      gmax=g(1)
c      gmin=g(1)
c      gmax=gsmax
c      gmin=gsmin
c      do i=2,ncol
c      if(g(i).gt.gmax) gmax=g(i)
c      if(g(i).lt.gmin) gmin=g(i)
c      end do
c      print *,'gmin = ',gmin,'  gmax = ',gmax
c      do i=1,ncol
c      g(i)=gsav(i)
c      end do
      prompt='enter constant to add to the profile values: '
      call tpromp(0,0,prompt)
c      print *, 'enter constant to add to the profile values'
      read(5,*) const
      const=const*sfact-soff
      prompt='calculating interface...'
      call tpromp(0,0,prompt)
      zave=0.
      do i=1,ncol
      zave=zave+z(i)
      g(i)=gsav(i)+const
c      gsav(i)=g(i)
      end do
      zave=zave/float(ncol)
      dz1=dz1-zave
c      gmin=gsmin+const
c      gmax=gsmax+const
cc      if(dx.eq.0.0e0) read (10) (x(i),i=1,ncol)
cc      read (10) zilch,(g(i),i=1,ncol)
c      rewind 10
cc      if(nz.eq.1.and.nrow.eq.1) go to 309
cc      print 308
c308   format (' fatal error, program accepts only single precision
c     & real-valued, profile data. check the nrow, nz parameters of
c     & the input standard file.')
cc      go to 999
c309   continue
cc      print 310, (id(i),i=1,14),pgm(1),pgm(2),rho,dz1,type,dx,nmax,eps
c      print 310, rho,dz1,type,dx,nmax,eps
cc310   format (///,1x,14a4,/,1x,2a4,/,1x,'density contrast =',f6.2,
c310   format (///,1x,'density contrast =',f6.2,
c     & 5x,'reference plane depth =',f8.2,' km. (',a5,' mode.)',/
c     & 1x,'grid interval = ',f8.3,' km.',5x,'number of iterations
c     & =',i3,5x,'epsilon = ',f5.2,/)
      pi=3.1415926536
c
c  gravity constant for kilometer grid:
c
      c1=41.908847e0*rho
      c2=6.67e0*rho
      imax=ncol
      rimax=ncol
      refhoriz(1)=dz1
      refhoriz(2)=refhoriz(1)
      n=1
c
c  calculate initial thicknesses and x(i), if dx.ne.0.
c
      do 1 i=1,imax
c      if(dx.ne.0.0e0) x(i)=x0+(i-1)*dx
      if (rho*g(i).lt.0.) g(i)=0.
      t(i)=g(i)/c1
      th(i)=0.0
      if(type.eq.'r'.and.t(i).gt.refhoriz(1)) t(i)=refhoriz(1)
      if(t(i).ge.0.0e0) go to 1
c      print 200
c200   format (' fatal error, rho and g must have same algebraic sign.')
      go to 999
1     continue
c     type 23
c23   format (//,'  iter     max delta    loc  total rms delta',/)
c      z1sq=dz1*dz1
      dz2=dz1
c      z2sq=z1sq
      mode=1
      if(type.eq.'r') mode=2
c lock keyboard and initialize array processor
c     write(*,311) 27,1
c311  format('+',a1,'KL',i1,$)
c     call apinit(0,0,key)
      call itgrav(ncol,mode,nmax,eps,c2,dz1,dz2,x,t,g,gn,
     1            n,emax,iemax,rms)
c     call aprlse
c unlock keyboard
c     write(*,311) 27,0
c      print 6, n,emax,iemax,rms
c6     format(2(i5,e16.6))
      do 12 i=1,imax
c      if(abs(gn(i)).ge.dval) gn(i)=0.0
c      if(abs(t(i)).ge.dval) go to 12
c     e=g(i)-gn(i,kk)
      if(type.eq.'b') go to 11
      th(i)=refhoriz(1)-t(i)
      go to 12
   11 th(i)=refhoriz(1)+t(i)
   12 continue
      do 22 k=1,ncol
      if(abs(gn(k)-g(k)).le.sfact) then
        if(g(k).ne.0.0) then
          call vchar(x(k),z(k)+th(k),sym,1,icode,size,0.,0.,0.)
        endif
      else
        call vchar(x(k),z(k)+th(k),sym1,1,icode,size,0.,0.,0.)
      endif
   22 continue
      prompt='try a different model? [n] '
      call tpromp(0,0,prompt)
c      print *, 'try a different model?'
      read 103, answer
      if(answer.eq.'y'.or.answer.eq.'Y') go to 2
   23 continue
c      do 25 j=1,ncol
c      diff=dval
c      do 21 k=1,kk
c      delta=abs(gsav(j)-gn(j,k))
c      if(delta.lt.diff) then
c        diff=delta
c        ksav=k
c      endif
c   21 continue
c      write(11,405) x(j),gsav(j),gn(j,ksav),ksav,th(j,ksav)
c   25 continue
c      do 26 j=1,ncol
c   26 write(11,406) x(j),gsav(j),(gn(j,k),k=1,kk)
c      do 27 j=1,ncol
c   27 write(11,406) x(j),(th(j,k),k=1,kk)
c      do 21 j=1,ncol
c      write(11,405) x(j),(th(j,k),k=1,kk)
c   21 continue
c  405 format(3(1x,f9.3),i5,1x,f9.3)
c  406 format(8(1x,f9.3))
c      call pltset(iplotr,xboard,yboard,1)
c      xp(4)=xboard
c      yp(4)=yboard
c      xp(3)=1.
c      yp(3)=1.
c      xp(2)=0.
c      yp(2)=0.
c      xp(1)=xp(4)-1.5
c      yp(1)=yp(4)-1.5
cc     xncol=ncol
c      dxp(1)=x0
c      dxp(2)=x0+(ncol-1)*dx
c      if(dx.eq.0.) dxp(2)=x(ncol)-x(1)
c      dyp(2)=0.0
c      dyp(1)=dxp(2)*yp(1)/xp(1)/vert
c      call scale (dxp,dyp,xp,yp,4,ier)
c      if (ier.eq.0) go to 85
c      print *,' Error in scaling parameters'
c      go to 999
c85    call xaxis (dxp,dyp,xp,1.,5,.12,'(f4.0)',4)
c      call vchar ((dxp(1)+dxp(2))*.5,dyp(1),'DISTANCE-KM',
c     & 11,2,.12,0.,-.36,-.7)
c      call yaxis (dyp,dxp,yp,1.,5,.12,'(f4.0)',4)
c      call vchar (dxp,(dyp(1)+dyp(2))*.5, 'KILOMETERS',
c     & 10,2,.12,1.5706,-.4,.75)
c      call vchar (1.,.1,'PRISM MODEL',11,3,.15,0.,0.,0.)
c      call neatl
c      do 22 k=1,kk
c      call vchar (x,th(1,k),k,ncol,1,.04,0.,0.,0.)
c      call line  (x,th(1,k),ncol,0,k*100)
c      do 22 k=1,ncol
c      call curv(x,th(1,k),ncol,k*100)
c      call line (dxp,refhoriz(1,k),2,0,k*100)
c   22 continue
c      call endpt (ie)
c      print*, ' plot gravity? answer y/n'
c      read (5,103) answer
103   format (a1)
c      if (answer.eq.'n') go to 990
c      call pltset(iplotr,xboard,yboard,1)
c      xp(4)=xboard
c      yp(4)=yboard
c      xp(3)=1.
c      yp(3)=1.
c      xp(2)=0.
c      yp(2)=0.
c      xp(1)=xp(4)-1.5
c      yp(1)=yp(4)-1.5
cc     xncol=ncol
c      dxp(1)=x0
c      dxp(2)=x0+(ncol-1)*dx
c      if(dx.eq.0.) dxp(2)=x(ncol)-x(1)
cc      print *,' Enter estimated minimum & maximum gravity values-mgals'
cc      read *, dyp(1),dyp(2)
c      dyp(1)=gmin
c      dyp(2)=gmax
c      call scale (dxp,dyp,xp,yp,4,ier)
c      if (ier.eq.0) go to 86
c      print *, ' error in scaling parameters'
c      go to 999
c86    call xaxis (dxp,dyp,xp,1.,5,.12,'(f4.0)',4)
c      call yaxis (dyp,dxp,yp,1.,5,.12,'(f5.0)',5)
c      call vchar ((dxp(1)+dxp(2))*.5,dyp(1),'DISTANCE-KM',
c     & 11,2,.12,0.,-.36,-.7)
c      call vchar (dxp,(dyp(1)+dyp(2))*.5,'MILLIGALS',
c     & 9,2,.12,1.5706,-.4,.75)
c      call vchar (1.,.1,'GRAVITY PROFILE',
c     & 15,3,.15,0.,0.,0.)
c      call neatl
c      call line (x,gsav,ncol,0,7)
c      do 24 k=1,kk
c      call vchar (x,gn(1,k),k,ncol,1,.08,0.,0.,0.)
c      call line  (x,gn(1,k),ncol,0,k*100)
cc      call curv(x,gn(1,k),ncol,k*100)
c   24 continue
c      call endpt (ie)
cc   print input parameters
c990   write (6,130) infile
c130   format (20x,' Input file=',a50)
c      write (6,132) type
c132   format (20x,' Mode=',a5)
c      write (6,134) refhoriz(1,kk)
c134   format (20x,' Reference Horizon=',f4.1)
c      write (6,136) rho
c136   format (20x,' Density contrast=',f4.2)
c      write (6,138) n,rms
c138   format (20x,' Iterations=',i3,' RMS error=',f4.2)
c      write (6,140) emax,iemax
c140   format (20x,' Max error=',f4.2,' Location=',i3)
999   return
      end
c**********************************************************************
      subroutine itgrav(ncol,mode,nmax,eps,c2,dz1,dz2,x,t,g,gn,
     1                  n,emax,iemax,rms)
      dimension x(ncol),t(ncol),g(ncol),gn(ncol)
      pi=3.14159265
      imax=ncol
      rimax=ncol
      n=1
      lim=ncol-1
2     continue
      do 4 i=1,ncol
      gn(i)=0.0
c
c  first point, note xi1 = -infinity.
c
      xi2=0.5*(x(2)+x(1))
      if(t(1).eq.0.0) go to 400
      go to (401,402),mode
401   dz2=dz1+t(1)
c      z2sq=dz2*dz2
      go to 403
402   dz1=dz2-t(1)
c      z1sq=dz1*dz1
403   x2=xi2-x(i)
      x2sq=x2*x2
      z2sq=dz2*dz2
      z1sq=dz1*dz1
      r2sq=x2sq+z1sq
      r3sq=x2sq+z2sq
      ang1=atan2(dz1,x2)-pi
      ang2=pi-atan2(dz2,x2)
      gn(i)=gn(i)+c2*(2.0*(dz1*ang1+dz2*ang2)-x2*alog(r2sq/r3sq))
c
c  2 to ncol-1 th body elements
c
400   do 3 j=2,lim
      xi1=xi2
      xi2=0.5*(x(j+1)+x(j))
      if(t(j).eq.0.0) go to 3
      go to (350,351),mode
350   dz2=dz1+t(j)
c      z2sq=dz2*dz2
      go to 352
351   dz1=dz2-t(j)
c      z1sq=dz1*dz1
352   x1=xi1-x(i)
      x2=xi2-x(i)
      z2sq=dz2*dz2
      z1sq=dz1*dz1
      ddx=xi2-xi1
      x3=x(i)-0.5*(xi1+xi2)
c
c  approximate formula
c  not used in this version of the program
c
c  exact formula
c
371   continue
      x1sq=x1*x1
      x2sq=x2*x2
      r1sq=x1sq+z1sq
      r2sq=x2sq+z1sq
      r3sq=x2sq+z2sq
      r3sq=x2sq+z2sq
      r4sq=x1sq+z2sq
      ang1=atan2(dz1,x2)-atan2(dz1,x1)
      ang2=atan2(dz2,x1)-atan2(dz2,x2)
      grav=2.0*(dz1*ang1+dz2*ang2)+x1*alog(r1sq/r4sq)-x2*alog
     &(r2sq/r3sq)
      gn(i)=gn(i)+c2*grav
    3 continue
c
c  ncol th body element
c
      xi1=xi2
      if(t(ncol).eq.0.0) go to 4
      go to (404,405),mode
404   dz2=dz1+t(ncol)
c      z2sq=dz2*dz2
      go to406
405   dz1=dz2-t(ncol)
c      z1sq=dz1*dz1
406   x1=xi1-x(i)
      x1sq=x1*x1
      z2sq=dz2*dz2
      z1sq=dz1*dz1
      r1sq=x1sq+z1sq
      r4sq=x1sq+z2sq
      ang1=-atan2(dz1,x1)
      ang2=atan2(dz2,x1)
      gn(i)=gn(i)+c2*(2.0*(dz1*ang1+dz2*ang2)+x1*alog(r1sq/r4sq))
    4 continue
      emax=0.0
      iemax=1
      rms=0.0
      do 5 i=1,imax
      e=abs(g(i)-gn(i))
      if(e.lt.emax) go to 5
      emax=e
      iemax=i
    5 rms=rms+e*e
      rms=sqrt(rms/rimax)
c     type 6, n,emax,iemax,rms
c   6 format (2(i5,e16.6))
      if(n.eq.nmax.or.eps.ge.emax) go to 10
    8 n=n+1
      do 9 i=1,imax
      if (gn(i).eq.0.) go to 9
      t(i)=(t(i)*g(i))/gn(i)
9     if(mode.eq.2.and.t(i).gt.dz2) t(i)=dz2
c   9 gn(i)=0.0
      go to 2
   10 continue
      return
      end
c**********************************************************************
c  terrace1.for
c
c   To steepen gradients on basis of local curvature of profile
c    field.
c   No dvals!
c   This version is limited to x,y,z profiles only.
c   10 July 87
c
      subroutine terrace1(iterations,ipen,nobs,g,x)
c      common /grvdata/ nobs,gmin,gmax,gx(512),gz(512),g(512)
c      character*50 file1, file2
c        character*1 ans
c      dimension x(1300),z(1300),a(1300),id(14),pgm(2)
      dimension x(nobs),g(nobs),a(512)
c      write(6,1)
c1       format(1x,'Enter input file name:'/1x'*'$)
c        read(5,101)file1
c101     format(a50)
c        write(6,2)
c2       format(1x,'Enter output file name:'/1x'*'$)
c        read(5,101)file2
c        write(6,7)
c7       format(1x,'Enter number of iterations;'/1x'*'$)
c        read(5,*)iterations
c        open(11,file=file1,form='formatted',status='old')
c        open (12,file=file2,form='formatted',status='new')
cc        write(6,40)
c40      format(1x,'Enter x0, dx, and nobs (three values)'/' *'$)
cc        read(5,*)x0,dx,nobs
c        icount=0
c        nobs=0
c10      nobs=nobs+1
c        read(11,*,end=11) x(nobs),z(nobs),a(nobs)
c        go to 10
c41      format(3g15.5)
c11      nobs=nobs-1
        do 10 i=1,nobs
        a(i)=g(i)
10      continue
	do 32 iter=1,iterations
	sum=0.0
	a0=a(1)
	do 20 i=2,nobs-1
	ap=a(i+1)
	am=a0
	a0=a(i)
	fx=.5*(ap-am)
	fxx=ap+am-a0-a0
	if(fx)17,17,18
17	if(fxx)22,23,21
18	if(fxx)21,23,22
21	a(i)=ap
	sum=sum+abs(a0-ap)
	go to 20
22	a(i)=am
	sum=sum+abs(a0-am)
	go to 20
c23      if(iter.eq.iterations)icount=icount+1
23      continue
20	continue
c        write(6,33)sum
c33      format(1x,'sum = ',f10.2)
	a(1)=a(2)
	a(nobs)=a(nobs-1)
32	continue
c        do 34 i=1,nobs
c        xi=i-1
c        x=x0+xi*dx
c34      print*,i,x(i),a(i)
      call curv(x,a,nobs,ipen)
      go to 999
99    write(6,98)
98    format(1x,'Can`t handle it.')
999	continue
c        icount=(icount*100)/(nobs)
c        write(6,30)icount
c30      format(1x,'Percent of flat slopes =',i5)
c        close (11)
c        close (12)
      return
      end
C=====================================================================
C     SUBROUTINE TRANS CARRIES OUT LINEAR TRANSFORMATIONS :
C       * UPWARD CONTINUATION
C       * LOW-PASS FILTER        FCD, ALPHAD
C       * HORIZONTAL DERIVATIVE    DT/DX
C       * VERTICAL DERIVATIVE      DT/DZ
C       * REDUCTION TO THE POLE    RPOP
c      SUBROUTINE TRANS(NPAD,ZUP,PI,CX,HX,VX,RPX,FCD,ALPHAD,DCOSX,DCOSZ)
cC=====================================================================
cc      IMPLICIT REAL*8(A-H,O-Z)
cc      COMPLEX*16 CX(1),HX(1),VX(1),RPX(1),PCOMP,RPOP,RPM
c      COMPLEX CX(1),HX(1),VX(1),RPX(1),RPOP
c      RPOP1=(DCOSX**2-DCOSZ**2)
c      RPOP2=4*(DCOSX**2)*(DCOSZ**2)
c      RPOP3=RPOP1/(RPOP1**2+RPOP2)
c      RPOP4=-1.0*(2.0*DCOSX*DCOSZ)/(RPOP1**2+RPOP2)
c      RPOP=-1.0*CMPLX(RPOP3,RPOP4)
c      LIM=NPAD/2+1
c      PCD=2.*PI*FCD
c      FIL=1.0
c      DO 100 I=1,LIM
c      J=I-1
c      P=J*2.*PI/(NPAD-1)
cC <<<<<<<<  LOW-PASS FILTER & UPWARD CONTINUATION   >>>>>>>>>>>>>
c      IF(P.GT.PCD) FIL=EXP(-ALPHAD*(P-PCD)**2)
c      IF(FIL.LT.0.0001)FIL=0.0001
c      FAC=EXP(-ZUP*P)
c      CX(I)=CX(I)*FAC*FIL
c      XR=REAL(CX(I))
c      XI=-1.0*IMAG(CX(I))
cC <<<<<<<<<<<<<<<<<<<   REDUCTION TO THE POLE  >>>>>>>>>>>>>>>>>>>
c      RPX(I)=RPOP*(CX(I))
c      RPR=REAL(RPX(I))
c      RPI=-1.0*IMAG(RPX(I))
cC<<<<<<<<<<<<<<<<<<<    HORIZONTAL DERIVATIVE  >>>>>>>>>>>>>>>
c      HR=P*IMAG(RPX(I))
c      HI=-1.0*P*REAL(RPX(I))
c      HX(I)=CMPLX(HR,HI)
c      HI=-1.0*HI
cC <<<<<<<<<<<<<<<<<<<   VERTICAL DERIVATIVE   >>>>>>>>>>>>>>>>>
c      VR=P*REAL(RPX(I))
c      VI=P*IMAG(RPX(I))
c      VX(I)=CMPLX(VR,VI)
c      VI=-1.0*VI
cC <<<<<<<<<<<<<<<<<<<                          >>>>>>>>>>>>>>>>>>>
c      IF(I.GT.1) CX((NPAD+2)-I)=CMPLX(XR,XI)
c      IF(I.GT.1) VX((NPAD+2)-I)=CMPLX(VR,VI)
c      IF(I.GT.1) RPX((NPAD+2)-I)=CMPLX(RPR,RPI)
c  100 IF(I.GT.1) HX((NPAD+2)-I)=CMPLX(HR,HI)
c      CALL FORK(NPAD,CX,-1.0)
c      CALL FORK(NPAD,HX,-1.0)
c      CALL FORK(NPAD,VX,-1.0)
c      CALL FORK(NPAD,RPX,-1.0)
c      RETURN
c      END
c******************************************************************************
      subroutine redpol(ainc,adec,strike,mode)
c
c  ainc, adec:  magnetic field inclination and declination in degrees
c  strike:      strike of the profile in degrees
c  mode:        0 for reduction to pole and pseudogravity
c               1 for pseudomag
c
      common /pscale/ xmin,xmax,zmin,zmax
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /pmagdata/ pmmin,pmmax,pmminsav,pmmaxsav,pminc,pm(512)
      common /pgrvdata/ pgmin,pgmax,pgminsav,pgmaxsav,pginc,pg(512)
      common /rmagdata/ rmin,rmax,rminsav,rmaxsav,rtp(512)
      common /contdata/ nco,comin,comax,cominsav,comaxsav,coinc,zup,
     1co(512)
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      common /work1/xc(512)
      common /work2/pgx(512)
      common /work3/pmx(512)
      common /work4/rpx(512)
      complex xc,pgx,pmx,rpx,cmplx,rpop
      character*80 gfile, sfile
c
      einc=ainc/57.2957795
      sd=(strike-adec)/57.2957795
      dcosx=cos(einc)*cos(sd)
      dcosz=sin(einc)
      rpop1=(dcosx**2-dcosz**2)
      rpop2=4*(dcosx**2)*(dcosz**2)
      rpop3=rpop1/(rpop1**2+rpop2)
      rpop4=-1.0*(2.0*dcosx*dcosz)/(rpop1**2+rpop2)
      rpop=-1.0*cmplx(rpop3,rpop4)
      if(mode.eq.0.) then
        call setfft(1,npad,n,a)
      else
        call setfft(1,npad,ng,g)
      endif
      pi=3.1415927
      lim=npad/2+1
      fil=1.0
      do 100 i=1,lim
      j=i-1
c      p=j*2.*pi/(npad-1)
      p=j*2.*pi/float(npad)
c <<<<<<<<  low-pass filter & upward continuation   >>>>>>>>>>>>>
c      if(p.gt.pcd) fil=exp(-alphad*(p-pcd)**2)
c      if(fil.lt.0.0001)fil=0.0001
c      fac=exp(-zup*p)
c      cx(i)=cx(i)*fac*fil
c      xr=real(cx(i))
c      xi=-1.0*imag(cx(i))
c <<<<<<<<<<<<<<<<<<<   reduction to the pole  >>>>>>>>>>>>>>>>>>>
      if(mode.eq.0) then
        rpx(i)=rpop*(xc(i))
        rpr=real(rpx(i))
        rpi=-1.0*aimag(rpx(i))
        if(i.gt.1) rpx((npad+2)-i)=cmplx(rpr,rpi)
c <<<<<<<<<<<<<<<<<<<   pseudogravity  >>>>>>>>>>>>>>>>>>>
        if(p.ne.0.) then
          pgx(i)=rpop*(xc(i))/(p*750.)
        else
          pgx(i)=0.
        endif
        pgr=real(pgx(i))
        pgi=-1.0*aimag(pgx(i))
        if(i.gt.1) pgx((npad+2)-i)=(cmplx(pgr,pgi))
      else
c <<<<<<<<<<<<<<<<<<<   pseudomag  >>>>>>>>>>>>>>>>>>>
        pmx(i)=750.*p*(xc(i))/rpop
        pmr=real(pmx(i))
        pmi=-1.0*aimag(pmx(i))
        if(i.gt.1) pmx((npad+2)-i)=cmplx(pmr,pmi)
      endif
cc<<<<<<<<<<<<<<<<<<<    horizontal derivative  >>>>>>>>>>>>>>>
c      hr=p*imag(rpx(i))
c      hi=-1.0*p*real(rpx(i))
c      hx(i)=cmplx(hr,hi)
c      hi=-1.0*hi
cc <<<<<<<<<<<<<<<<<<<   vertical derivative   >>>>>>>>>>>>>>>>>
c      vr=p*real(rpx(i))
c      vi=p*imag(rpx(i))
c      vx(i)=cmplx(vr,vi)
c      vi=-1.0*vi
cc <<<<<<<<<<<<<<<<<<<                          >>>>>>>>>>>>>>>>>>>
c  100 if(i.gt.1) hx((npad+2)-i)=cmplx(hr,hi)
  100 continue
      if(mode.eq.0) then
        call fork(npad,pgx,-1.0)
        call fork(npad,rpx,-1.0)
        rmin=1.e+38
        rmax=-1.e+38
        pgmin=1.e+38
        pgmax=-1.e+38
        rminsav=1.e+38
        rmaxsav=-1.e+38
        pgminsav=1.e+38
        pgmaxsav=-1.e+38
      else
        pmmin=1.e+38
        pmmax=-1.e+38
        pmminsav=1.e+38
        pmmaxsav=-1.e+38
        call fork(npad,pmx,-1.0)
      endif
c      call fork(npad,vx,-1.0)
c      print *,rpop
c      zup=0.0
c      fcd=0.0
c      alphad=0.0
c  994 CALL TRANS(NPAD,ZUP,PI,xc,HX,VX,RPX,FCD,ALPHAD,DCOSX,DCOSZ)
      if(mode.eq.0) then
        npad=n
      else
       npad=ng
      endif
      DO 107 I=1,NPAD
c      DATA(I)=REAL(CX(I))
c      HDATA(I)=REAL(HX(I))
c      VDATA(I)=REAL(VX(I))
      if(mode.eq.0) then
        rtp(I)=REAL(RPX(I))
        if(rtp(i).gt.rmaxsav) rmaxsav=rtp(i)
        if(rtp(i).lt.rminsav) rminsav=rtp(i)
        pg(i)=real(pgx(i))
        if(pg(i).gt.pgmaxsav) pgmaxsav=pg(i)
        if(pg(i).lt.pgminsav) pgminsav=pg(i)
        if(x(i).lt.xmin) go to 107
        if(x(i).gt.xmax) go to 107
        if(rtp(i).gt.rmax) rmax=rtp(i)
        if(rtp(i).lt.rmin) rmin=rtp(i)
        if(pg(i).gt.pgmax) pgmax=pg(i)
        if(pg(i).lt.pgmin) pgmin=pg(i)
      else
        pm(i)=real(pmx(i))
        if(pm(i).gt.pmmaxsav) pmmaxsav=pm(i)
        if(pm(i).lt.pmminsav) pmminsav=pm(i)
        if(gx(i).lt.xmin) go to 107
        if(gx(i).gt.xmax) go to 107
        if(pm(i).gt.pmmax) pmmax=pm(i)
        if(pm(i).lt.pmmin) pmmin=pm(i)
      endif
  107 continue
c  107 CX(I)=CMPLX(XREAL(I),XIMAG(I))
c  999 CONTINUE
c      NDIF=(NPAD-N)/2
c      DO 1151 I=1,N
c      HDATA(I)=HDATA(I+NDIF)
c      VDATA(I)=VDATA(I+NDIF)
c      rtp(I)=rtp(I+NDIF)
c 1151 continue
c 1151 DATA(I)=DATA(I+NDIF)
c      k=lc/2
c      do 30 j=2,k
cc      xc(j)=xc(j)+xc(j)
cc   30 xc(lc+2-j)=0.
c      xc(i)=xc(i)*rpop
c      rpr=real(xc(i))
c      rpi=-1.0*aimag(xc(i))
c   30 xc(lc+2-j)=cmplx(rpr,rpi)
c      call fork(lc,xc,+1.)
c      do 31 i=1,n
c      rtp(i)=real(xc(i))
c   31 continue
c      print *, rmin,rmax
c      print *, rmin,rmax
      if(mode.eq.0) then
        call dplot(2,xmin,xmax,rmin,rmax)
        call curv(x,rtp,n,500)
        irmag=1
        if(igrv.gt.0) then
          pginc=2.*(pgmax-pgmin)/(gmax-gmin)
        else
          pginc=2.
        endif
        call dplot(2,xmin,xmax,pgmin,pgmax)
        call curv(x,pg,n,100)
        ipgrv=1
      else
        if(imag.gt.0) then
          pminc=50.*(pmmax-pmmin)/(amax-amin)
        else
          pminc=2.
        endif
        call dplot(2,xmin,xmax,pmmin,pmmax)
        call curv(gx,pm,ng,600)
        ipmag=1
      endif
c
      return
      end
c******************************************************************************
c******************************************************************************
      subroutine setfft(mflg,lc,n,a)
c      common /magdata/ n,xmin,xmax,amin,amax,zmin,zmax,
c     &x(512),z(512),a(512)
      common /work1/xc(512)
      common /work3/v(512),aa(512)
      dimension we(30),a(n)
      complex xc
      real*8 ap,xind
      data ncoef/16/
c remove mean
      ap=0.
      if(mflg.eq.1) then
        do 5 i=1,n
    5   ap=ap+a(i)
        ap=ap/float(n)
      endif
      do 11 i=1,n
c      if(x(i).gt.amx) amx=x(i)
c      if(x(i).lt.amn) amn=x(i)
   11 xc(i)=a(i)-ap
c  determine prediction filter
c      pi=3.1415927
      do 20 i=1,n
      aa(i)=a(i)
   20 v(i)=a(i)
      we(1)=1.
      do 23 j=2,ncoef
      ap=0.
      xind=0.
      do 21 i=j,n
      ap=ap+aa(i)*aa(i)+v(i-j+1)*v(i-j+1)
   21 xind=xind+aa(i)*v(i-j+1)
      rc=-2.*xind/ap
      do 22 i=j,n
      temp=aa(i)
      aa(i)=aa(i)+rc*v(i-j+1)
   22 v(i-j+1)=v(i-j+1)+rc*temp
      we(j)=0.
      jh=(j+1)/2
      do 23 i=1,jh
      k=j-i+1
      temp=we(k)+rc*we(i)
      we(i)=we(i)+rc*we(k)
   23 we(k)=temp
c  extend profile
      k=n+n/2
      do 24 i=1,11
      lc=2**i
      if(lc.gt.k) go to 25
   24 continue
   25 if(lc.gt.512) then
        k=n
        do 26 i=1,11
        lc=2**i
        if(lc.gt.k) go to 27
   26   continue
   27   if(lc.gt.512) stop 'profile too long'
      endif
      lp=n+1
      jh=(lp+lc)/2
      do 28 i=lp,jh
      xc(i)=0.
      xc(lc-i+lp)=0.
      do 28 j=2,ncoef
      xc(i)=xc(i)-we(j)*xc(i-j+1)
   28 xc(lc-i+lp)=xc(lc-i+lp)-we(j)*xc(mod(lc-i+lp+j-2,lc)+1)
c  fourier transform the complex profile
      call fork(lc,xc,+1.)
      return
      end
c******************************************************************************
      subroutine clear
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
      common /pscale/ xmin,xmax,zmin,zmax
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      common /pmagdata/ pmmin,pmmax,pmminsav,pmmaxsav,pminc,pm(512)
      common /pgrvdata/ pgmin,pgmax,pgminsav,pgmaxsav,pginc,pg(512)
      common /rmagdata/ rmin,rmax,rminsav,rmaxsav,ra(512)
      common /contdata/ nco,comin,comax,cominsav,comaxsav,coinc,zup,
     1co(512)
      common /calc1/icalc,icalcd,gc(512),hc(512)
      character gfile*80, ifile*80, sfile*80
      call setmod(3)
      icalcd=0
      if(igrv.gt.0) then
        ipen=300
        call dplot(0,xmin,xmax,gmin,gmax)
        call curv(gx,g,ng,ipen)
        if(imag.gt.0) then
          call dplot(2,xmin,xmax,amin,amax)
          ipen=400
          call curv(x,a,n,ipen)
        endif
      else
        call dplot(0,xmin,xmax,amin,amax)
        ipen=400
        call curv(x,a,n,ipen)
      endif
      if(irmag.eq.1) then
        call dplot(2,xmin,xmax,rmin,rmax)
        ipen=500
        call curv(x,ra,n,ipen)
      endif
      if(ipmag.eq.1) then
        call dplot(2,xmin,xmax,pmmin,pmmax)
        ipen=600
        call curv(gx,pm,ng,ipen)
      endif
      if(ipgrv.eq.1) then
        call dplot(2,xmin,xmax,pgmin,pgmax)
        ipen=100
        call curv(x,pg,n,ipen)
      endif
      if(icont.gt.0) then
        call dplot(2,xmin,xmax,comin,comax)
        ipen=200
c       if(nco.eq.n) then
       if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
         call curv(x,co,n,ipen)
       else
         call curv(gx,co,ng,ipen)
       endif
      endif
      call dplot(1,xmin,xmax,zmin,zmax)
      if(imag.gt.0) then
        ipen=400
        call curv(x,z,n,ipen)
      endif
      if(igrv.gt.0) then
        ipen=300
c        call curv1(gx,gz,ng,ipen,2)
        call curv(gx,gz,ng,ipen)
      endif
      call rescale
      if(icont.gt.0) then
        ipen=200
c        if(nco.eq.n) then
        if(icont.eq.1.or.icont.eq.3.or.icont.eq.5) then
          call curv2(x,z,zup,n,ipen)
        else
          call curv2(gx,gz,zup,ng,ipen)
        endif
      endif
c      call menu(nmenu,label,mpen)
      if(iopn.gt.0) then
        icolor=500
        call plotmod
      endif
      if(gfile(1:1).ne.' ') then
        open(10,file=gfile,form='formatted',status='old',err=163)
        call plotgeo
        close(10)
      endif
      if(sfile(1:1).ne.' ') then
        open(10,file=sfile,form='formatted',status='old',err=163)
        call plotsect
        close(10)
      endif
  163 return
      end
c******************************************************************************
      subroutine spline(f,p,n,dx)
      common /work3/c(512),d(512)
      dimension f(n),p(n)
      c(1)=.5
      d(1)=1.5*(f(2)-f(1))/dx
      n1=n-1
      do 10 i=2,n1
      dnm=(4.0-c(i-1))*dx
c      print *,'dnm = ',dnm,' d(i-1) = ',d(i-1)
      c(i)=dx/dnm
   10 d(i)=(3.0*(f(i+1)-f(i-1)) - dx*d(i-1))/dnm
c      print *,'c(n1) = ',c(n1)
      p(n)=(3.*(f(n)-f(n1))/dx - d(n1))/(2.-c(n1))
      do 9 i=1,n1
      k=n-i
    9 p(k)=d(k)-c(k)*p(k+1)
      return
      end
c ****************************************************************************
      subroutine writemod
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      character string*68, ifile*80, jfile*80
      if(iopn.eq.0) return
c      if(iedit.ne.1) return
c      return
      rewind(11)
      nrec=ncor+1
      write(12,201,rec=nrec)
   10 read(11,200,end=90) string
      nrec=nrec+1
      write(12,200,rec=nrec) string
      go to 10
   90 if(iopn.eq.3) then
        close(10)
        open(10,file=ifile,form='formatted',status='old')
        close(11)
        len=index(ifile,'.')
        jfile=ifile(1:len)//'bak'
        open(11,file=jfile,status='unknown',form='formatted',err=93)
   91   read(10,200,end=92) string
        len=len_trim(string)
        write(11,200) string(1:len)
        go to 91
   92   close(11)
   93   close(10,status='delete')
        open(10,file=ifile,form='formatted',status='new')
        do 94 i=1,nrec
        read(12,200,rec=i,end=95) string
        len=len_trim(string)
        write(10,200) string(1:len)
   94   continue
   95   close(10)
      endif
   99 close(12)
c      close(11,status='delete')
c      iopn=0
  200 format(a)
  201 format('>> BODY SPECS')
      return
      end

