      subroutine setedt(ixo,iyo)
      common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      common /plot/ dxp(2),yp(4),xpix,ypix,xscale,yscale
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
      character label(10)*8, prompt*79, ifile*80
      dimension mpen(10)
      data nmenu/10/,label/' new mod',' editmod',' new body',
     &'del body',' new cor',' del cor','move cor','  prop  ',
     &'  calc  ',' return  '/
      data mpen/0,0,0,0,0,0,0,0,0,0/
      icolor=600
      call menu(nmenu,label,mpen)
c      ixo=600
c      iyo=300
c  100 prompt=' '
c      call tpromp(0,0,prompt)
  100 call menudvr(xdum,y,ifunc,ixo,iyo)
      go to (100,101,102,103,104,105,106,107,108,109,110) ifunc+1
c new model
  101 if(iopn.gt.0) then
        prompt='click on "newmod" again to write out and erase the curre
     1nt model'
        call tpromp(0,0,prompt)
        call menudvr(xdum,y,ifunc,ixo,iyo)
        go to (100,201,100,100,100,100,100,100,100,100,100) ifunc+1
  201   call writemod
        call tdraw(1)
        call plotmod
        call tdraw(0)
        iopn=0
      endif
      call newmod
      go to 100
c edit model
  102 if(iopn.gt.0) then
        prompt='click on "editmod" again to write out and erase the curr
     1ent model'
        call tpromp(0,0,prompt)
        call menudvr(xdum,y,ifunc,ixo,iyo)
        go to (100,100,202,100,100,100,100,100,100,100,100) ifunc+1
  202   call writemod
        call tdraw(1)
        call plotmod
        call tdraw(0)
        iopn=0
      endif
      call editmod(ixo,iyo)
      icolor=500
      if(iopn.gt.0) call plotmod
      go to 100
c new body
  103 call newbod(ixo,iyo)
      go to 100
c delete body
  104 call delbod(ixo,iyo)
      go to 100
c new corner
  105 call newcor(ixo,iyo)
      go to 100
c delete corner
  106 call delcor(ixo,iyo)
      go to 100
c move corner
  107 call movcor(ixo,iyo)
      go to 100
c change properties
  108 call prop(ixo,iyo)
      go to 100
c calculate profiles
  109 call calc
      go to 100
c return
  110 call tdraw(1)
      call menu(nmenu,label,mpen)
      call tdraw(0)
      prompt=' '
      call tpromp(0,0,prompt)
c      call tpromp(0,2,0,prompt)
      return
      end
c ****************************************************************************
      subroutine newmod
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
        common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      character prompt*79, ifile*80, gfile*80, sfile*80
      prompt='enter new model file name: '
      call tpromp(0,0,prompt)
      iopn=0
      do 7 i=1,50
      rho(i)=0.
      sus(i)=0.
      pl(i,1)=0.
      pl(i,2)=0.
      rem(i,1)=0.
      rem(i,2)=0.
      rem(i,3)=0.
    7 continue
      do 8 i=1,150
      ncrn(i)=0
      loca(i)=0
    8 continue
      ifile=' '
      read(5,200,end=9) ifile
  200 format(a)
    9 if(ifile(1:1).eq.' ') go to 55
      close(12)
      close(11,status='delete')
c      iedit=1
      open(12,file=ifile,status='new',form='formatted',
     &access='direct',recl=68,err=10)
      go to 11
   10 prompt='file exists, OK to overwrite? :'
      call tpromp(0,0,prompt)
      if(noyes().eq.0) go to 998
      prompt=' '
      call tpromp(0,0,prompt)
      open(12,file=ifile,status='old',form='formatted',
     &access='direct',recl=68,err=999)
      close(12,status='delete')
      open(12,file=ifile,status='new',form='formatted',
     &access='direct',recl=68)
   11 open(11,file='pdepth.bod',status='unknown',form='formatted')
      nbod=0
      ncor=0
c      ef(1)=0.0
c      if(igrv.gt.0) then
c        do 12 i=1,ng
c        xb(i)=gx(i)
c        zb(i)=gz(i)
c   12   write(12,203,rec=i) gx(i),gz(i)
c        ncor=ng
c      endif
      iopn=1
  203 format(2g16.6)
      go to 999
   55 prompt=' '
      call tpromp(0,0,prompt)
      go to 999
  998 prompt='error creating file'
      call tpromp(0,0,prompt)
  999 return
      end
c ****************************************************************************
      subroutine editmod
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
        common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      character prompt*79, ifile*80, aline*80, blank*80
      do 7 i=1,50
      rho(i)=0.
      sus(i)=0.
      pl(i,1)=0.
      pl(i,2)=0.
      rem(i,1)=0.
      rem(i,2)=0.
      rem(i,3)=0.
    7 continue
      do 8 i=1,150
      ncrn(i)=0
      loca(i)=0
    8 continue
c      ef(1)=0.0
      blank=' '
      if(ifile(1:1).eq.' ') then
        prompt='enter name of model file to edit: '
      else
        len=len_trim(ifile)
        write(prompt,101) ifile(1:len)
  101   format('enter name of model file to edit: [',a,'] ')
      endif
      call tpromp(0,0,prompt)
      iopn=0
      prompt=' '
      read(5,199,end=3) prompt
  199 format(a)
    3 if(prompt(1:1).eq.' ') then
        if(ifile(1:1).eq.' ') go to 55
      else
        ifile=prompt
      endif
      close(10)
      close(11,status='delete')
      close(12)
c      iedit=0
      open(10,file=ifile,status='old',form='formatted',err=998)
c      len=index(ifile,'.')
c      jfile=ifile(1:len)//'bak'
c      open(11,file=jfile,status='unknown',form='formatted',err=15)
c      go to 16
c   15 jfile=ifile(1:len)//'bk1'
c      open(11,file=jfile,status='unknown',form='formatted',err=998)
c   16 continue
c
        iunit=10
c        junit=11
        call countv(ncor,iunit,iend)
        if(iend.eq.1) then
          prompt='eof on model file'
          call tpromp(0,0,prompt)
          return
        endif
        nbod=0
4       nbod=nbod+1
5       read(iunit,'(a)',end=6) aline
       if(aline.eq.blank) go to 5
c       len=len_trim(aline)
c       write(junit,'(a)') aline(1:len)
        call countv(ncrn(nbod),iunit,iend)
        if(iend.eq.1) go to 6
        go to 4
6       nbod=nbod-1
        close(11)
        rewind 10
        ncor=ncor/2
c  read vertex coordinates
        read(10,*,err=66) (xb(i),zb(i),i=1,ncor)
c  skip blank lines, find body list delimiter
        do 10  i=1,5
          read(10,*,err=11) delimt
10      continue         
c  read physical property and vertex lookup lists
c  put in read for pc Microsoft Fortran
11     read(10,'(a)') aline
c        ivert=0
        is=1
        do 20 ibody=1,nbod
        loca(ibody)=is-1
        ie=loca(ibody)+ncrn(ibody)
c  blank line before next body spec?
        read(10,13) aline
13      format(a)
        if (aline.ne.blank) backspace (10)
        read(10,*,err=77,end=99) rho(ibody),sus(ibody),pl(ibody,1),
     1 pl(ibody,2), rem(ibody,1),rem(ibody,2),rem(ibody,3)
        read(10,*,err=88,end=99) (icrn(i),i=is,ie)
c        xm=xb(icrn(is))
c        ym=zb(icrn(is))
c        call line(xm,ym,1,0,500)
c        do 32 j=is+1,ie
c        xm=xb(icrn(j))
c        ym=zb(icrn(j))
c        call line(xm,ym,1,1,500)
c   32   continue
c        xm=xb(icrn(is))
c        ym=zb(icrn(is))
c        call line(xm,ym,1,1,500)
        is=ie+1
20      continue
c
c  get body tags
c
c      close(10,status='delete')
      close(10)
c      ifile=' '
c      prompt='enter name of model file to create [pdepth.cor]: '
c      call tpromp(0,0,prompt)
c      read(5,199,end=29) ifile
c   29 if(ifile(1:1).eq.' ') ifile='pdepth.cor'
c      open(12,file=ifile,status='new',form='formatted',
      open(12,file='pdepth.cor',status='unknown',form='formatted',
     &access='direct',recl=68)
c      go to 15
c   14 prompt='file exists, OK to overwrite? :'
c      call tpromp(0,0,prompt)
c      if(noyes().eq.0) go to 998
c      prompt=' '
c      call tpromp(0,0,prompt)
c      open(12,file=ifile,status='old',form='formatted',
c     &access='direct',recl=68,err=999)
c      close(12,status='delete')
c      open(12,file=ifile,status='new',form='formatted',
c     &access='direct',recl=68)
c   15 continue
      do 30 i=1,ncor
      call checkcor(xb(i),zb(i))
      write(12,203,rec=i) xb(i),zb(i)
   30 continue
      call writebod
  203 format(2g16.6)
c
      iopn=2
   55 prompt=' '
      call tpromp(0,0,prompt)
      go to 999
   66 prompt='vertex read error'
      call tpromp(0,0,prompt)
      go to 999
   77 prompt='phys-prop input error'
      call tpromp(0,0,prompt)
      go to 999
   88 write(prompt,200) ibody
  200 format('body corner read error, body ',i5)
c      prompt='body corner read error'
      call tpromp(0,0,prompt)
      go to 999
   99 prompt='eof in body definition'
      call tpromp(0,0,prompt)
      go to 999
  998 prompt='error opening model file'
      call tpromp(0,0,prompt)
  999 return
      end
c******************************************************************************
        subroutine countv(ival,iunit,iend)
c  count words delimited by spaces or commas
c  return when non-alpanumeric character begins word
        character a*200
        iend=0
        ival=0
10      a=' '
        read(iunit,20,end=99) a
c        len=len_trim(a)
c        write(junit,20) a(1:len)
20      format(a200)
30      na=leftj(a)
        if(na.eq.0) go to 10
        if(a(1:1).eq.',') then
          if(na.gt.1) a(1:na-1)=a(2:na)
          a(na:na)=' '
          go to 30
        endif
        ich=ichar(a(1:1))
        if(ich.le.42 .or. ich.ge.123) return
        if(ich.ge.58 .and. ich.le.64) return
        if(ich.ge.91 .and. ich.le.96) return
        mb=index(a,' ')
        mc=index(a,',')
        m=mb
        if(mc.gt.0 .and. mc.lt.mb) m=mc 
        ival=ival+1
        na1=na-m+1
        if(na1.eq.0) go to 10
        a(1:na1)=a(m:na)
        do 40 i=na1+1,200
40      a(i:i)=' '
        go to 30
99      iend=1
        return
        end
c******************************************************************************
        function leftj(a)
c  left justifies a string and returns the position
c  of the last nonblank character
        character a*(*)
        n=len(a)
        if(a(1:1).ne.' ') go to 15
        do 1 m=2,n
1       if(a(m:m).ne.' ') go to 5
        leftj=0
        return
c
5       i2=0
        do 10 i=m,n
        i2=i2+1
        a(i2:i2)=a(i:i)
10      continue
        do 11 i3=i2+1,n
11      a(i3:i3)=' '
        n=n-m+1
c
15      do 20 leftj=n,1,-1
20      if(a(leftj:leftj).ne.' ') go to 25
25      return
        end
c
c ****************************************************************************
      subroutine writebod
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      character fmt1*24,ifile*80
      open(11,file='pdepth.bod',status='unknown',form='formatted')
c
c
      do 31 i=1,nbod
      write(11,200) rho(i),sus(i),pl(i,1),pl(i,2),(rem(i,j),j=1,3),i
      nleft=ncrn(i)
      is=loca(i)+1
      ie=loca(i)+ncrn(i)
   50 if(nleft.le.10) go to 60
      write(11,201) (icrn(j),j=is,is+9)
      is=is+10
      nleft=nleft-10
      go to 50
   60 write(fmt1,202) nleft
      write(11,fmt1) (icrn(j),j=is,ie),i
   31 continue
      return
  200 format((1x,f6.3), (1x,1pg11.3), 0p2(1x,f8.1),
     1 (1x,1pg11.3), 0p2(1x,f6.1), ' >P',i2)
  201 format(10i4)
  202 format('(',i2,'i4,'' >B'',i2)')
      end
c ****************************************************************************
      subroutine plotmod
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
        character ifile*80
c        ivert=0
c        is=1
        do 10 i=1,ncor
        xm=xb(i)
        ym=zb(i)
        call vchar(xm,ym,2,1,icolor+1,.05,0.,0.,0.)
   10   continue
        do 20 ibody=1,nbod
c        loca(ibody)=is-1
        is=loca(ibody)+1
        ie=loca(ibody)+ncrn(ibody)
c  blank line before next body spec?
c        read(10,13) aline
c13      format(a)
c        if (aline.ne.blank) backspace (10)
c        read(10,*,err=77,end=99) rho(ibody),sus(ibody),pl(ibody,1),
c     1 pl(ibody,2), rem(ibody,1),rem(ibody,2),rem(ibody,3)
c        read(10,*,err=88,end=99) (icrn(i),i=is,ie)
        xm=xb(icrn(is))
        ym=zb(icrn(is))
c        print *,is,icrn(is)
c        call vchar(xm,ym,2,1,icolor+1,.05,0.,0.,0.)
        call line(xm,ym,1,0,icolor)
        do 32 j=is+1,ie
        xm=xb(icrn(j))
        ym=zb(icrn(j))
c        call vchar(xm,ym,2,1,icolor+1,.05,0.,0.,0.)
        call line(xm,ym,1,1,icolor)
   32   continue
        xm=xb(icrn(is))
        ym=zb(icrn(is))
c        call vchar(xm,ym,2,1,icolor+1,.05,0.,0.,0.)
        call line(xm,ym,1,1,icolor)
c        is=ie+1
20      continue
        if(ibod.eq.1) then
          is=loca(nbod+1)+1
          ie=loca(nbod+1)+ncrn(nbod+1)
c          print *,nbod,is,ie
          if(is.gt.500) then
            print*,"is = ",is
            return
          else if (icrn(is).gt.712) then
            print*,"icrn(is) = ",icrn(is)
            return
          endif
          xm=xb(icrn(is))
          ym=zb(icrn(is))
c          call vchar(xm,ym,2,1,icolor+1,.05,0.,0.,0.)
          call line(xm,ym,1,0,icolor)
          do 33 j=is+1,ie
          xm=xb(icrn(j))
          ym=zb(icrn(j))
c          call vchar(xm,ym,2,1,icolor+1,.05,0.,0.,0.)
          call line(xm,ym,1,1,icolor)
   33     continue
        endif
        return
        end
c ****************************************************************************
      subroutine newbod(ixo,iyo)
c        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
      dimension icor(200)
      character prompt*79, fmt1*24, ifile*80
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(iopn.eq.2) iopn=3
      prompt='click on points (clockwise); click "newbod" to close'
      call tpromp(0,0,prompt)
      if(ibod.eq.0) then
        j=0
        if(nbod.eq.0) then
          is=0
        else
          is=loca(nbod)+ncrn(nbod)
        endif
        loca(nbod+1)=is
      else
        call line(xmo,ymo,1,0,500)
      endif
      ibod=1
   10 j=j+1
      is=is+1
   11 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.eq.3) go to 99
      if (ifunc.eq.10) then
        is=is-1
        j=j-1
        prompt='input of corners suspended, click "newbod" to continue'
        call tpromp(0,0,prompt)
        go to 999
      endif
      xmo=xm
      ymo=ym
      if (ifunc.ne.0) go to 11
c      call line(xm,ym,1,j-1,500)
      if(ncor.eq.0) go to 30
c      rewind(10)
      eps=delx/2.
      do 20 i=1,ncor
      read(12,203,rec=i) xc,yc
      dist=sqrt((xc-xmo)**2+(yc-ymo)**2)
      if(dist.lt.eps) then
        eps=dist
        xm=xc
        ym=yc
        isav=i
c        go to 40
      endif
   20 continue
      if(eps.lt.delx/2.) go to 40
   30 ncor=ncor+1
      isav=ncor
      call checkcor(xm,ym)
      write(12,203,rec=ncor) xm,ym
      xb(ncor)=xm
      zb(ncor)=ym
   40 if(j.eq.1) then
        xm1=xm
        ym1=ym
      endif
      icor(j)=isav
      ncrn(nbod+1)=ncrn(nbod+1)+1
      icrn(is)=isav
      call line(xm,ym,1,j-1,500)
      go to 10
   99 prompt='body completed'
      call tpromp(0,0,prompt)
      ibod=0
c      iopn=2
      if(j.lt.4) go to 998
      nbod=nbod+1
      call line(xm1,ym1,1,j-1,500)
      write(11,200) nbod
      nleft=j-1
      is=1
   50 if(nleft.le.10) go to 60
      write(11,201) (icor(i),i=is,is+9)
      is=is+10
      nleft=nleft-10
      go to 50
   60 write(fmt1,202) nleft
      write(11,fmt1) (icor(i),i=is,j-1),nbod
      go to 999
  200 format('0.00 0.00e-03 0 0 0 0 0 >P',i2)
  201 format(10i4)
  202 format('(',i2,'i4,'' >B'',i2)')
  203 format(2g16.6)
c      call tgin(ix,iy,ixo,iyo)
c      x=dxp(1)+float(ix)/float(nhor)*(dxp(2)-dxp(1))
c      y=dyp(1)+float(iy)/float(nnvert)*(dyp(2)-dyp(1))
  998 prompt='not enough corners to define a body'
      call tpromp(0,0,prompt)
  999 return
      end
c ****************************************************************************
      subroutine checkcor(xm,ym)
      common /magdata/ n,amin,amax,x(512),z(512),a(512)
      common /grvdata/ ng,gmin,gmax,gx(512),gz(512),g(512)
      logical inside
c
c   check for vertices above terrain (observations)
      if(n.le.0.or.n.gt.510) go to 15
        x(n+1)=x(n)
        x(n+2)=x(1)
        z(n+1)=-100.
        z(n+2)=-100.
10      if(inside(n+2,x,z,xm,ym)) then
          ym=ym+.001
          go to 10
        endif
15    if(ng.le.0.or.ng.gt.510) go to 30
        gx(ng+1)=gx(ng)
        gx(ng+2)=gx(1)
        gz(ng+1)=-100.
        gz(ng+2)=-100.
20        if(inside(ng+2,gx,gz,xm,ym)) then
          ym=ym+.001
          go to 20
        endif
30      return
        end
c ****************************************************************************
      subroutine delbod(ixo,iyo)
c deletes a body; does not delete the body corners yet.
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      dimension icor(200)
      common /swork/ x(512),y(512)
      character prompt*79, ifile*80
        logical inside
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(nbod.lt.1) then
        prompt='no bodies available'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(iopn.eq.2) iopn=3
      prompt='click on body to delete'
      call tpromp(0,0,prompt)
   11 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) go to 99
      do 12 jbod=1,nbod
      is=loca(jbod)+1
      n=ncrn(jbod)
      do 13 i=1,n
      j=icrn(is+i-1)
      x(i)=xb(j)
      y(i)=zb(j)
   13 continue
c      print *,is,n
      if(inside(n,x,y,xm,ym)) go to 14
c      print*,xm,(xb(i),i=is,is+n-1)
c      print*,ym,(zb(i),i=is,is+n-1)
c      if(incpt(xm,ym,x,y,n).eq.0) go to 14
   12 continue
      go to 11
   14   ie=loca(jbod)+ncrn(jbod)
        xm=xb(icrn(is))
        ym=zb(icrn(is))
        call line(xm,ym,1,0,1500)
        do 32 j=is+1,ie
        xm=xb(icrn(j))
        ym=zb(icrn(j))
        call line(xm,ym,1,1,1500)
   32   continue
        xm=xb(icrn(is))
        ym=zb(icrn(is))
        call line(xm,ym,1,1,1500)
      prompt='click on body again to confirm'
      call tpromp(0,0,prompt)
   21 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) then
         kbod=-1
         go to 24
      endif
      do 22 kbod=1,nbod
      is=loca(kbod)+1
      n=ncrn(kbod)
      do 23 i=1,n
      j=icrn(is+i-1)
      x(i)=xb(j)
      y(i)=zb(j)
   23 continue
c      print *,is,n
      if(inside(n,x,y,xm,ym)) go to 24
c      print*,xm,(xb(i),i=is,is+n-1)
c      print*,ym,(zb(i),i=is,is+n-1)
c      if(incpt(xm,ym,x,y,n).eq.0) go to 24
   22 continue
c      go to 21
      kbod=-1
   24 if(kbod.ne.jbod) then
        icolor=500
        call plotmod
        go to 99
      endif
       close(11,status='delete')
c      print*,'deleting body ',jbod
c   14 is=loca(jbod)+1
c      ie=loca(jbod)+ncrn(jbod)
c      iee=loca(nbod)+ncrn(nbod)
c      i=is-1
cc      do 17 i=is,ie
c   17 i=i+1
c      if(i.gt.ie) go to 18
c      do 15 j=1,is-1
c      if(icrn(j).eq.icrn(i)) go to 17
c   15 continue
c      do 16 j=ie+1,iee
c      if(icrn(j).eq.icrn(i)) go to 17
c   16 continue
c      print *,i
c      call delcor1(i)
c      i=i-1
c      ie=ie-1
c      iee=iee-1
c      go to 17
c   18 continue
      call tdraw(1)
      call plotmod
      call tdraw(0)
      icolor=500
c
c      ie=loca(jbod)+ncrn(jbod)
c      i=loca(jbod)
c   19 i=i+1
c      if(i.gt.ie) go to 30
c      do 21 j=1,nbod
c      if(j.eq.jbod) go to 21
c      iss=loca(j)+1
c      iee=loca(j)+ncrn(j)
c      do 20 ii=iss,iee
c      if(icrn(ii).eq.icrn(i)) go to 19
c   20 continue
c   21 continue
c      print*,'deleting corner ',i
c      call delcor1(i)
c      i=i-1
c      go to 19
   30 continue
c
      nbod=nbod-1
      do 40 i=jbod,nbod
      loca(i)=loca(i+1)
      ncrn(i)=ncrn(i+1)
      rho(i)=rho(i+1)
      sus(i)=sus(i+1)
      pl(i,1)=pl(i+1,1)
      pl(i,2)=pl(i+1,2)
      rem(i,1)=rem(i+1,1)
      rem(i,2)=rem(i+1,2)
      rem(i,3)=rem(i+1,3)
   40 continue
c      call plotmod
      ncrn(nbod+1)=0
c      i=0
c   18 i=i+1
c      if(i.gt.ncor) go to 999
c      do 19 j=1,nbod
c      is=loca(j)+1
c      ie=loca(j)+ncrn(j)
c      do 19 k=is,ie
c      if(icrn(k).eq.i) go to 18
c   19 continue
c      call delcor1(i)
c8      i=i-1
      call plotmod
      call writebod
c      iopn=2
      go to 99
c      go to 18
c      open(11,file='pdepth.bod',status='unknown',form='formatted')
c      do 19 j=1,nbod
c      write(11,200) rho(j),sus(j),pl(j,1),pl(j,2),(rem(j,i),i=1,3),j
c        nleft=ncrn(j)
c        is=loca(j)+1
c        ie=loca(j)+ncrn(j)
c   20   if(nleft.le.10) go to 21
c        write(11,201) (icrn(i),i=is,is+9)
c        is=is+10
c        nleft=nleft-10
c        go to 20
c   21   write(fmt1,202) nleft
c        write(11,fmt1) (icrn(i),i=is,ie),j
c   19 continue
  200 format((1x,f6.3), (1x,1pg11.3), 0p2(1x,f8.1),
     1 (1x,1pg11.3), 0p2(1x,f6.1), ' >P',i2)
  201 format(10i4)
  202 format('(',i2,'i4,'' >B'',i2)')
   99   prompt=' '
        call tpromp(0,0,prompt)
  999 return
      end
c******************************************************************************
        logical function inside(n,x,z,xp,zp)
c   is xp,zp inside polygon defined by x,z ?
        dimension x(1),z(1)
        logical z1,z2
        inside=.false.
        do 5 i=1,n
        i1=i+1
        if(i1.gt.n) i1=1
        dz=z(i1)-z(i)
        if(dz.eq.0.0) go to 5
        z1 = zp.le.z(i)
        z2 = zp.gt.z(i1)
        if((z1.and.z2) .or. (.not.z1 .and. .not.z2)) go to 3
        go to 5
c   points on boundary are outside with condition .lt.
c   for inclusion change to .le.
3       rslope=(x(i1)-x(i))/dz
        d=(xp-x(i))-(zp-z(i))*rslope
c        if(abs(d).lt.1.e-4) d=0.0
        if(abs(d).lt.1.e-6) d=0.0
c        print*,xp,zp,d
        if(d.lt.0.0) inside = .not.inside
5       continue
        return
        end
cC+
cC      INCPT   - finds intercepts of the horizontal line thru a body
cC
cC    Call INCPT (x,y,xbod,ybod,npt)
cC
cC    x,y= point to find nearest intercepts to
cC    xbod,ybod= two-dimensional arrays with corners of body (dim 50,50)
cC    npt= number of points in body
cC    nbd= body number (first subscript of xbod, ybod)
cC    x1,x2= return for closest intercepts to x,y
cC
cC    RETURNS:  0- x,y inside
cC              1- x,y outside
cC             -1- x,y on parallel edge of body
cC             -2- no intercept or just one intercept
cC-
cC  *********************************************
c      Integer Function incpt(x,y,xbod,ybod,npt)
c      Dimension xbod(npt),ybod(npt),xint(100)
cC
cC     statement function to find intersection of line with horizontal
cC
c      xsect(a1,b1,a2,b2,b)=(b-b1)*((a2-a1)/(b2-b1))+a1
cC
cC     int is number of intersections, ipar # of parallel lines found
cC
c      int=1
c      ipar=0
c      Do 50 i=1,npt-1
c      xp1=xbod(i)
c      yp1=ybod(i)
c      xp2=xbod(i+1)
c      yp2=ybod(i+1)
c      Call insect(x,y,xp1,yp1,xp2,yp2,xint,int,ipar)
c   50 Continue
c      xp1=xbod(npt)
c      yp1=ybod(npt)
c      xp2=xbod(1)
c      yp2=ybod(1)
c      Call insect(x,y,xp1,yp1,xp2,yp2,xint,int,ipar)
cC
cC     examine int and ipar and report results
cC
c      If (ipar.EQ.0)Go To 60
cC
cC     on parallel edge
Cc
cc      x1=xint(1)
cc      x2=xint(2)
c      incpt=-1
c      Return
c   60 Continue
c      If (int.GT.1)Go To 70
cC
cC     no intercept
cC
cc      x1=x
cc      x2=x
c      incpt=-2
c      Return
c   70 Continue
c      If (int.GT.2)Go To 80
cC
cC     only one intercept
cC
cc      x1=xint(1)
cc      x2=xint(2)
c      incpt=-2
c      Return
c   80 Continue
cC
cC     two or more intercepts
cC
c   85 isort=0
cC
cC     bubble sort of intercepts
cC
c      Do 100 j=1,int-2
c      If (xint(j).LE.xint(j+1))Go To 100
c      xtemp=xint(j+1)
c      xint(j+1)=xint(j)
c      xint(j)=xtemp
c  100 Continue
c      If (isort.NE.0)Go To 85
cC
cC     count # of intersections less than x
cC
cC
cC    check for duplicate intersections
cC
c      l=1
c  190 Continue
c      If (l.GE.int)Go To 195
c      If (abs(xint(l+1)-xint(l)).GT.1e-5)Go To 192
c      If (l+1.EQ.int)Go To 193
c      Do 191 l1=l+1,int-1
c      xint(l1)=xint(l1+1)
c  191 Continue
c  193 int=int-1
c  192 l=l+1
c      Go To 190
c  195 Continue
c      If (int.LE.1)Go To 70
c      Do 200 k=1,int-1
c      If (xint(k).GE.x)Go To 210
c  200 Continue
c  210 Continue
c      If (k.EQ.1)Go To 400
c      If (mod(k-1,2).EQ.0)Go To 350
cC
cC     odd number of intersections to the left- inside
cC
cc      x1=xint(k-1)
cc      x2=xint(k)
c      incpt=0
c      Return
c  350 Continue
cC
cC     even number to the left, outside
cC
cc      x1=xint(k-2)
cc      x2=xint(k-1)
c      incpt=1
c      Return
c  400 Continue
cC
cC     no intersections to the left, outside
cC
cc      x1=xint(k)
cc      x2=xint(k+1)
c      incpt=1
c      Return
c      End
cC
cC     INSECT  finds intersection with individual side
cC
cC  *********************************************
c      Subroutine insect(x,y,a1,b1,a2,b2,xint,int,ipar)
c      Dimension xint(100)
cC
cC     statement function to find intersection
cC
c      xsect(x1,y1,x2,y2,yloc)=(yloc-y1)*((x2-x1)/(y2-y1))+x1
cC
cC     check if y is between b1 and b2
cC
c      If ((b1.LE.y).AND.(b2.GE.y))Go To 100
c      If ((b1.GE.y).AND.(b2.LE.y))Go To 100
c      Return
c  100 Continue
c      If (abs(b2-b1).lt.1.e-20)Go To 105
c      xint(int)=xsect(a1,b1,a2,b2,y)
c      int=int+1
c      Return
c  105 Continue
cC
cC     parallel
cC
c      xint(1)=a1
c      xint(2)=a2
c      ipar=ipar+1
c      Return
c      End
c******************************************************************************
      subroutine newcor(ixo,iyo)
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      dimension icor(200)
      character prompt*79, ifile*80
c        logical inside
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(nbod.lt.1) then
        prompt='no bodies available'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(iopn.eq.2) iopn=3
      prompt='click on one endpoint of the line segment to split'
      call tpromp(0,0,prompt)
      eps=delx/2.
   11 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) go to 99
      x1=1.e38
      y1=1.e38
      do 12 i=1,ncor
      read(12,203,rec=i) xc,yc
      dist=sqrt((xc-xm)**2+(yc-ym)**2)
      if(dist.lt.eps) then
        eps=dist
        x1=xc
c      print *,dist,eps,icor
        y1=yc
        icor=i
c        go to 30
      endif
   12 continue
      if(x1.eq.1.e38) go to 11
   30 prompt='click on the other endpoint of the line segment to split'
      call tpromp(0,0,prompt)
      prompt=' '
      eps=delx/2.
   13 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) go to 99
      x2=1.e38
      y2=1.e38
      do 14 i=1,ncor
      read(12,203,rec=i) xc,yc
      dist=sqrt((xc-xm)**2+(yc-ym)**2)
      if(dist.lt.eps) then
        eps=dist
        x2=xc
c      print *,dist,eps,icor
        y2=yc
        jcor=i
c        go to 40
      endif
   14 continue
      if(x2.eq.1.e38) go to 13
   40 if(jcor.eq.icor) go to 99
      call tdraw(1)
      call plotmod
      call tdraw(0)
      icolor=500
      ncor=ncor+1
      xb(ncor)=(x2+x1)/2.
      zb(ncor)=(y2+y1)/2.
      write(12,203,rec=ncor) xb(ncor),zb(ncor)
c
      iee=0
      do 50 i=1,nbod
      is=loca(i)+1
      ie=loca(i)+ncrn(i)
      do 51 j=is,ie-1
      if(icrn(j).eq.icor.and.icrn(j+1).eq.jcor) go to 52
      if(icrn(j).eq.jcor.and.icrn(j+1).eq.icor) go to 52
   51 continue
      j=ie
      if(icrn(is).eq.icor.and.icrn(ie).eq.jcor) go to 52
      if(icrn(is).eq.jcor.and.icrn(ie).eq.icor) go to 52
      go to 50
   52 iee=loca(nbod)+ncrn(nbod)
      do 53 ii=iee+1,j+1,-1
   53 icrn(ii)=icrn(ii-1)
      icrn(j+1)=ncor
      ncrn(i)=ncrn(i)+1
      if(i.eq.nbod) go to 50
      do 54 ii=i+1,nbod
   54 loca(ii)=loca(ii)+1
   50 continue
      if(iee.eq.0) prompt='something wrong'
      call plotmod
c
c   61 close(11,status='delete')
c      open(11,file='pdepth.bod',status='unknown',form='formatted')
cc      ifirst=1
c      do 19 j=1,nbod
c      write(11,200) rho(j),sus(j),pl(j,1),pl(j,2),(rem(j,i),i=1,3),j
c        is=loca(j)+1
c        ie=loca(j)+ncrn(j)
c        nleft=ncrn(j)
c   20   if(nleft.le.10) go to 21
c        write(11,201) (icrn(i),i=is,is+9)
c        is=is+10
c        nleft=nleft-10
c        go to 20
c   21   write(fmt1,202) nleft
c        write(11,fmt1) (icrn(i),i=is,ie),j
c   19 continue
       close(11,status='delete')
       call writebod
c       iopn=2
cc      entry delcor1(icor)
cc      print *,icor,xb(icor),zb(icor)
   99  call tpromp(0,0,prompt)
c      print *,'at 40',icor
  999 return
c  200 format((1x,f6.3), (1x,1pg11.3), 0p2(1x,f8.1),
c     1 (1x,1pg11.3), 0p2(1x,f6.1), ' >P',i2)
c  201 format(10i4)
c  202 format('(',i2,'i4,'' >B'',i2)')
  203 format(2g16.6)
      end
c******************************************************************************
      subroutine delcor(ixo,iyo)
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      dimension icor(200)
      character prompt*79, ifile*80
c        logical inside
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(ncor.lt.1) then
        prompt='no corners available'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(iopn.eq.2) iopn=3
      prompt='click on corner to delete'
      call tpromp(0,0,prompt)
      eps=delx/2.
   11 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) go to 99
      dist1=1.e38
c      print*,ncor
      do 12 icor=1,ncor
      read(12,203,rec=icor) xc,yc
      dist=sqrt((xc-xm)**2+(yc-ym)**2)
      if(dist.lt.dist1) then
        dist1=dist
        icorsav=icor
        xn=xc
c      print *,dist,dist1,eps,icor
        yn=yc
      endif
   12 continue
      if(dist1.gt.eps) go to 11
      call vchar(xn,yn,2,1,1501,.05,0.,0.,0.)
      prompt='click on corner again to confirm'
      call tpromp(0,0,prompt)
   13 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) then
        call vchar(xn,yn,2,1,501,.05,0.,0.,0.)
        go to 99
      endif
      dist1=1.e38
c      print*,ncor
      do 14 icor=1,ncor
      read(12,203,rec=icor) xc,yc
      dist=sqrt((xc-xm)**2+(yc-ym)**2)
      if(dist.lt.dist1) then
        dist1=dist
        jcorsav=icor
c        xn=xc
c      print *,dist,dist1,eps,icor
c        yn=yc
      endif
   14 continue
c      if(dist1.gt.eps) go to 13
      if(jcorsav.ne.icorsav.or.dist1.gt.eps) then
        call vchar(xn,yn,2,1,501,.05,0.,0.,0.)
        go to 99
      endif
   40 call tdraw(1)
      call plotmod
      call tdraw(0)
      icolor=500
      call delcor1(icorsav)
      call plotmod
      call writebod
c      iopn=2
c      entry delcor1(icor)
c      print *,icor,xb(icor),zb(icor)
   99 prompt=' '
      call tpromp(0,0,prompt)
c      print *,'at 40',icor
  999 return
  203 format(2g16.6)
      end
c******************************************************************************
      subroutine delcor1(icor)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      dimension icor(200)
      character prompt*79, fmt1*24, ifile*80
      ncor=ncor-1
      do 50 i=icor,ncor
      xb(i)=xb(i+1)
      zb(i)=zb(i+1)
c      icrn(i)=icrn(i+1)
      read(12,203,rec=i+1) xc,yc
      write(12,203,rec=i) xc,yc
   50 continue
c      ind=0
c      ind1=0
      iee=loca(nbod)+ncrn(nbod)
cc      do 58 i=1,nbod
cc   58 print *,i,loca(i),ncrn(i)
cc      do 59 i=1,ie
cc   59 print *,i,icrn(i)
cc      print *
c      do 60 i=1,ie
c      if(i+ind.gt.ie) go to 61
c      if(icrn(i).eq.icor) ind1=1
cc        if(ind1.eq.0) ind1=1
cc        ind=ind+1
cc      endif
c      if(icrn(i+ind).eq.icor) ind=ind+1
c      icrn(i)=icrn(i+ind) - ind1
cc      print *,i,icrn(i)
c   60 continue
   61 close(11,status='delete')
c      open(11,file='pdepth.bod',status='unknown',form='formatted')
c      ifirst=1
      do 19 j=1,nbod
c      write(11,200) rho(j),sus(j),pl(j,1),pl(j,2),(rem(j,i),i=1,3),j
        is=loca(j)+1
        ie=loca(j)+ncrn(j)
        i=is
   51   if(icrn(i).eq.icor) then
c          ifirst=0
          ncrn(j)=ncrn(j)-1
          ie=ie-1
          if(ncrn(j).le.2) then
            prompt='Warning: zero area body '
            call tpromp(0,0,prompt)
          endif
c          print *
c          print *,i,icrn(i),j,ncrn(j)
          do 515 ii=j+1,nbod
  515     loca(ii)=loca(ii)-1
          iee=iee-1
          do 516 ii=i,iee
  516     icrn(ii)=icrn(ii+1)
          i=i-1
        endif
        i=i+1
        if(i.le.ie) go to 51
        do 52 i=is,ie
        if(icrn(i).gt.icor) icrn(i)=icrn(i)-1
   52   continue
        nleft=ncrn(j)
        ie=loca(j)+ncrn(j)
   20   if(nleft.le.10) go to 21
c        write(11,201) (icrn(i),i=is,is+9)
        is=is+10
        nleft=nleft-10
        go to 20
   21   write(fmt1,202) nleft
c        write(11,fmt1) (icrn(i),i=is,ie),j
   19 continue
  999 return
  200 format((1x,f6.3), (1x,1pg11.3), 0p2(1x,f8.1),
     1 (1x,1pg11.3), 0p2(1x,f6.1), ' >P',i2)
  201 format(10i4)
  202 format('(',i2,'i4,'' >B'',i2)')
  203 format(2g16.6)
      end

c******************************************************************************
      subroutine movcor(ixo,iyo)
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      dimension icor(200)
      character prompt*79, ifile*80
c        logical inside
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
c        call tpromp(0,0,prompt)
        go to 999
      endif
      if(ncor.lt.1) then
        prompt='no corners available'
c        call tpromp(0,0,prompt)
        go to 999
      endif
      if(iopn.eq.2) iopn=3
      prompt='click on corner to move'
      call tpromp(0,0,prompt)
      prompt=' '
      eps=delx/2.
   11 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) go to 999
      dist1=1.e38
      do 12 icor=1,ncor
      read(12,203,rec=icor) xc,yc
      dist=sqrt((xc-xm)**2+(yc-ym)**2)
      if(dist.lt.dist1) then
        dist1=dist
        icorsav=icor
        xn=xc
c      print *,dist,eps,icor
        yn=yc
      endif
   12 continue
      if(dist1.gt.eps) go to 11
      call vchar(xn,yn,2,1,1501,.05,0.,0.,0.)
   40 prompt='click on new position'
      call tpromp(0,0,prompt)
      prompt=' '
      call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) then
        call vchar(xn,yn,2,1,501,.05,0.,0.,0.)
        go to 999
      endif
      call tdraw(1)
      call plotmod
      call tdraw(0)
      icolor=500
      call checkcor(xm,ym)
      xb(icorsav)=xm
      zb(icorsav)=ym
      write(12,203,rec=icorsav) xm,ym
      call plotmod
c      iopn=2
c      call writebod
c      entry delcor1(icor)
c      print *,icor,xb(icor),zb(icor)
  999 call tpromp(0,0,prompt)
c      print *,'at 40',icor
      return
  203 format(2g16.6)
      end
c ****************************************************************************
      subroutine prop(ixo,iyo)
c sets physical property values
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c      common /data/ n,xmin,xmax,amin,amax,gmin,gmax,zmin,zmax,
c     &x(512),z(512),a(512),gx(512),gz(512),g(512)
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
c      dimension icor(200)
      dimension x(100),y(100)
      character prompt*79, fmt*10, ifile*80
        logical inside
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(nbod.lt.1) then
        prompt='no bodies available'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(iopn.eq.2) iopn=3
      prompt='click on body'
      call tpromp(0,0,prompt)
   11 call menudvr(xm,ym,ifunc,ixo,iyo)
      if (ifunc.ne.0) go to 999
      do 12 jbod=1,nbod
      is=loca(jbod)+1
      n=ncrn(jbod)
      do 13 i=1,n
      j=icrn(is+i-1)
      if(i.gt.100) then
        print*,"i = ",i
        go to 12
      else if (j.gt.712) then
        print*,"j = ",j
        go to 12
      endif
      x(i)=xb(j)
      y(i)=zb(j)
   13 continue
c      print *,is,n
      if(inside(n,x,y,xm,ym)) go to 14
c      print*,xm,(xb(i),i=is,is+n-1)
c      print*,ym,(zb(i),i=is,is+n-1)
c      if(incpt(xm,ym,x,y,n).eq.0) go to 14
   12 continue
      go to 11
   14 close(11,status='delete')
c      print*,'deleting body ',jbod
c   14 is=loca(jbod)+1
c      ie=loca(jbod)+ncrn(jbod)
c      iee=loca(nbod)+ncrn(nbod)
c      i=is-1
cc      do 17 i=is,ie
c   17 i=i+1
c      if(i.gt.ie) go to 18
c      do 15 j=1,is-1
c      if(icrn(j).eq.icrn(i)) go to 17
c   15 continue
c      do 16 j=ie+1,iee
c      if(icrn(j).eq.icrn(i)) go to 17
c   16 continue
c      print *,i
c      call delcor1(i)
c      i=i-1
c      ie=ie-1
c      iee=iee-1
c      go to 17
c   18 continue
c      call tdraw(1)
c      call plotmod
c      call tdraw(0)
c      icolor=500
c
c      ie=loca(jbod)+ncrn(jbod)
c      i=loca(jbod)
c   19 i=i+1
c      if(i.gt.ie) go to 30
c      do 21 j=1,nbod
c      if(j.eq.jbod) go to 21
c      iss=loca(j)+1
c      iee=loca(j)+ncrn(j)
c      do 20 ii=iss,iee
c      if(icrn(ii).eq.icrn(i)) go to 19
c   20 continue
c   21 continue
c      print*,'deleting corner ',i
c      call delcor1(i)
c      i=i-1
c      go to 19
c   30 continue
c
c      nbod=nbod-1
c      do 40 i=jbod,nbod
c      loca(i)=loca(i+1)
c      ncrn(i)=ncrn(i+1)
      write(prompt,300) rho(jbod)
      call tpromp(0,0,prompt)
c      read(5,*,end=40) rho(jbod)
      read (*,'(a)') prompt
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) rho(jbod)
      endif 
   40 write(prompt,301) sus(jbod)
      call tpromp(0,0,prompt)
c      read(5,*,end=41) sus(jbod)
      read (*,'(a)') prompt
      if(prompt(1:1).eq.'\') go to 46
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) sus(jbod)
      endif 
   41 write(prompt,302) pl(jbod,1)
      call tpromp(0,0,prompt)
c      read(5,*,end=41) sus(jbod)
      read (*,'(a)') prompt
      if(prompt(1:1).eq.'\') go to 46
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) pl(jbod,1)
      endif 
   42 write(prompt,303) pl(jbod,2)
      call tpromp(0,0,prompt)
c      read(5,*,end=41) sus(jbod)
      read (*,'(a)') prompt
      if(prompt(1:1).eq.'\') go to 46
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) pl(jbod,2)
      endif 
   43 write(prompt,304) rem(jbod,1)
      call tpromp(0,0,prompt)
      read(*,'(a)') prompt
      if(prompt(1:1).eq.'\') go to 46
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) rem(jbod,1)
      endif 
   44 write(prompt,305) rem(jbod,2)
      call tpromp(0,0,prompt)
      read(*,'(a)') prompt
      if(prompt(1:1).eq.'\') go to 46
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) rem(jbod,2)
      endif 
   45 write(prompt,306) rem(jbod,3)
      call tpromp(0,0,prompt)
      read(*,'(a)') prompt
      if(prompt(1:1).eq.'\') go to 46
      leng=len_trim(prompt)
      if(leng.ne.0) then
        lenm=leng-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=999) rem(jbod,3)
      endif
   46 continue
c      rho(i)=rho(i+1)
c      sus(i)=sus(i+1)
c      pl(i,1)=pl(i+1,1)
c      pl(i,2)=pl(i+1,2)
c      rem(i,1)=rem(i+1,1)
c      rem(i,2)=rem(i+1,2)
c      rem(i,3)=rem(i+1,3)
c   40 continue
c      call plotmod
c      ncrn(nbod+1)=0
c      i=0
c   18 i=i+1
c      if(i.gt.ncor) go to 999
c      do 19 j=1,nbod
c      is=loca(j)+1
c      ie=loca(j)+ncrn(j)
c      do 19 k=is,ie
c      if(icrn(k).eq.i) go to 18
c   19 continue
c      call delcor1(i)
c8      i=i-1
c      call plotmod
      prompt=' '
      call tpromp(0,0,prompt)
      call writebod
c      iopn=2
      return
c      go to 18
c      open(11,file='pdepth.bod',status='unknown',form='formatted')
c      do 19 j=1,nbod
c      write(11,200) rho(j),sus(j),pl(j,1),pl(j,2),(rem(j,i),i=1,3),j
c        nleft=ncrn(j)
c        is=loca(j)+1
c        ie=loca(j)+ncrn(j)
c   20   if(nleft.le.10) go to 21
c        write(11,201) (icrn(i),i=is,is+9)
c        is=is+10
c        nleft=nleft-10
c        go to 20
c   21   write(fmt1,202) nleft
c        write(11,fmt1) (icrn(i),i=is,ie),j
c   19 continue
  100 format(a)
  200 format((1x,f6.3), (1x,1pg11.3), 0p2(1x,f8.1),
     1 (1x,1pg11.3), 0p2(1x,f6.1), ' >P',i2)
  201 format(10i4)
  202 format('(',i2,'i4,'' >B'',i2)')
c  300 format('change density from ',f10.5,' to :')
c  301 format('change susceptibility from ',f10.5,' to :')
  300 format('new density [',f10.5,']:')
  301 format('new susceptibility [',f10.5,']:')
  302 format('new near-side extent (positive) [',f10.5,']:')
  303 format('new far-side extent (negative) [',f10.5,']:')
c  304 format(f10.5)
  304 format('new remanent intensity [',f10.5,']:')
  305 format('new remanent inclination [',f10.5,']:')
  306 format('new remanent declination [',f10.5,']:')
  999 return
      end
c ****************************************************************************
      subroutine calc
      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 /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      common /switch/ imag,igrv,irmag,ipmag,ipgrv,gfile,sfile,icont
      common /model/nbod,ncor,iopn,delx,ibod,ifile
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
      common /calc1/icalc,icalcd,gc(512),hc(512)
      character prompt*79, gfile*80, ifile*80, sfile*80
      data isym/2/
        data d2r/1.745329e-2/
      if(iopn.eq.0) then
        prompt='open a model file first using "newmod" or "editmod"'
        call tpromp(0,0,prompt)
        go to 999
      endif
      if(ef(1).eq.0.0.and.imag.gt.0) then
5       prompt='enter field magnitude, inc, dec, and profile azimuth: '
        call tpromp(0,0,prompt)
        read(5,*,end=5,err=5) ef,azmuth
      endif
100     ei=ef(2)*d2r
	ed=(ef(3)-azmuth)*d2r
	eh=cos(ei)
	ev(1)=cos(ed)*eh
	ev(2)=sin(ed)*eh
	ev(3)=sin(ei)
	do 110 i=1,3
110     tdir(i)=ev(i)
      if(icalcd.eq.1) then
        call tdraw(1)
        if(igrv.gt.0) then
c          ipen=300
          call dplot(2,xmin,xmax,gmin,gmax)
          call curv1(gx,gc,ng,icolor,isym)
          if(imag.gt.0) then
            call dplot(2,xmin,xmax,amin,amax)
c            ipen=400
            call curv1(x,hc,n,icolor,isym)
          endif
        else
          call dplot(2,xmin,xmax,amin,amax)
c          ipen=400
          call curv1(x,hc,n,icolor,isym)
        endif
        call tdraw(0)
      endif
      prompt='calculating model field(s)...'
      call tpromp(0,0,prompt)
      call dplot(1,xmin,xmax,zmin,zmax)
      call rescale
      call taldvr(rmsg,rmsh,aveg,avem,gc,hc,ierr)
      if(ierr.eq.1) then
        icalc=0
        icalcd=0
        go to 999
      endif
      icalc=1
      icalcd=1
c      call clear
      if(igrv.gt.0) then
        ipen=300
        call dplot(2,xmin,xmax,gmin,gmax)
        call curv(gx,g,ng,ipen)
        call curv1(gx,gc,ng,ipen,isym)
        if(imag.gt.0) then
          call dplot(2,xmin,xmax,amin,amax)
          ipen=400
          call curv(x,a,n,ipen)
          call curv1(x,hc,n,ipen,isym)
        endif
      else
        call dplot(2,xmin,xmax,amin,amax)
        ipen=400
        call curv(x,a,n,ipen)
        call curv1(x,hc,n,ipen,isym)
      endif
      call dplot(1,xmin,xmax,zmin,zmax)
      call rescale
        if(igrv.gt.0.and.imag.gt.0) then
          write(prompt,205) rmsg,aveg,rmsh,avem
        else if(igrv.gt.0) then
          write(prompt,206) rmsg,aveg
        else if(imag.gt.0) then
           write(prompt,207) rmsh,avem
        endif
        call tpromp(0,0,prompt)
205     format('error, shift: grav =',1pg11.3,1pg11.3,
     1' mag =',1pg11.3,1pg11.3)
206     format('gravity rms error, shift = ',1pg13.3,1pg13.3)
207     format('magnetic rms error, shift = ',1pg13.3,1pg13.3)
999   return
      end
c ****************************************************************************
        subroutine taldvr(rmsg,rmsm,aveg,avem,gc,hc,ierr)
c   basic forward calculation, the talwani driver cycles through
c   the body array, removes a dc level, and figures rms error.
c        common /parm1/ nvert,nbody,nobs1,mobs1,grvswt,magswt,
c     1 iave,jave,vex,idepth,iparsh,indiv,idum
c        common /swork/ xb1(200),zb1(200),dum(8020)
c        common /model/ xb(200),zb(200),rho(50),sus(50),pl(50,2),
c     1 rem(50,3),dum1(150)
c        common /vertx/ icrn(500),ncrn(150),loca(150),nbord,ibord(200)
c        dimension xg(1),zg(1),go(1),gc(1),gerr(1),
c     1 xh(1),zh(1),ho(1),hc(1),herr(1)
        logical inside
      character prompt*79
      character*56 bfile,mfile,grfile,hfile
c
      common /model/nbod,ncor,iopn,delx,ibod,ifile
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
      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 /swork/ xb1(512),zb1(512)
      common /parms/ bfile,mfile,grfile,hfile,efield,einc,edec,jplotr,
     1               iave,jave
      character gfile*80, ifile*80, sfile*80
      dimension gc(512),hc(512)
        ierr=0
        if(igrv.gt.0) then
        do 1 i=1,ng
1       gc(i)=0.0
        rmsg=1.e30
        endif
        if(imag.gt.0) then
        do 2 i=1,n
2       hc(i)=0.0
        rmsm=1.e30
        endif
        call ckmodl(2,ng,gx,gz,ierr)
c       if(grvswt) call ckmodl(1,nobs,xg,zg,ierr)
        if(imag.gt.0) call ckmodl(1,n,x,z,ierr)
        if(ierr.eq.1) return
c
c  check for overlapping bodies
        ic=1
        do 100 j=1,nbod
        do 5 i=1,ncrn(j)
        xb1(i)=xb(icrn(loca(j)+i))
5       zb1(i)=zb(icrn(loca(j)+i))
        do 6 k=1,ncor
        do 9 l=1,ncrn(j)
        if(k.eq.icrn(loca(j)+l)) go to 6
9       continue
        if(.not.inside(ncrn(j),xb1,zb1,xb(k),zb(k))) go to 6
      call vchar(xb(k),zb(k),2,1,1501,.05,0.,0.,0.)
        prompt='Error -- highlighted vertex is inside a body'
        call tpromp(0,0,prompt)
        ierr=1
c        print  7,k,j
c7       format(' vertex',i3,' is inside body',i3)
        return
6       continue
c
c  call forward routines
        j1=j
        if(igrv.gt.0) call pickmg(0,j1,rho(j),ncrn(j),xb1,zb1,
     1 ng,gx,gz,gc)
        if(imag.gt.0)  call pickmg(1,j1,sus(j),ncrn(j),xb1,zb1,
     1 n,x,z,hc)
100     continue
c
c  gravity errors
        if(igrv.eq.0) go to 30
        if(iave.eq.0) go to 19
        aveg=0.0
        do 16 i=1,ng
16      aveg=aveg+(g(i)-gc(i))
        aveg=aveg/float(ng)
c        print 17,ave
c17      format(' removing',1pg15.3,' mgals from calculated gravity')
        do 18 i=1,ng
18      gc(i)=gc(i)+aveg
19      errt=0.0
        do 20 i=1,ng
c        gerr(i)=g(i)-gc(i)
        gerr=g(i)-gc(i)
c20      errt=errt+gerr(i)*gerr(i)
20      errt=errt+gerr*gerr
        rmsg=sqrt(errt/float(ng))
c
c  magnetic errors
30      if(imag.eq.0) return
        if(jave.eq.0) go to 39
        avem=0.0
        do 36 i=1,n
36      avem=avem+(a(i)-hc(i))
        avem=avem/float(n)
c        print 37,ave
c37      format(' removing',1pg15.3,' nTesla from calculated magnetics')
        do 38 i=1,n
38      hc(i)=hc(i)+avem
39      errt=0.0
        do 40 i=1,n
c        herr(i)=a(i)-hc(i)
        herr=a(i)-hc(i)
c40      errt=errt+herr(i)*herr(i)
40      errt=errt+herr*herr
        rmsm=sqrt(errt/float(n))
        return
        end
c******************************************************************************
        subroutine pickmg(mgswt,ibody,rho,nvert,xb,zb,
     1 nfld,xf,zf,pf)
c  switching routine for the individual forward algorithms
c  rho is the passed linear parameter (density or susceptibility)
c  prism lengths and remanent magnetic values are passed thru common
c  xb,zb is one body coordinate array
c  xf,zf array of locations to be used
        common /model1/ dum(1524),pl(50,2),rf(50,3)
c        common /model/ dum(500),pl(50,2),rf(50,3),dum1(150)
        common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
        dimension yl(2),vm(3),xb(1),zb(1),xf(1),zf(1),pf(1)
        data d2r/1.745329e-2/
        yl(1)=pl(ibody,1)
        yl(2)=pl(ibody,2)
        if(mgswt.eq.0) go to 10
        rinc=d2r*rf(ibody,2)
        rdec=d2r*(rf(ibody,3)-azmuth)
        rl=cos(rdec)*cos(rinc)
        rm=sin(rdec)*cos(rinc)
        rn=sin(rinc)
c  convert emu/cc to gamma
        rem=rf(ibody,1)*1.e5
c  polarization vector
        vm(1)=rho*ef(1)*ev(1)+rem*rl
        vm(2)=rho*ef(1)*ev(2)+rem*rm
        vm(3)=rho*ef(1)*ev(3)+rem*rn
        ylen=abs(yl(1))
        if(ylen.eq.0.0) ylen=1000.
        call mag2hd(vm,tdir,ylen,nvert,xb,zb,nfld,xf,zf,pf)
        return
10      if(yl(1).eq.0. .and. yl(2).eq.0.) go to 20
        call grv2hd(rho,yl,nvert,xb,zb,nfld,xf,zf,pf)
        return
20      call talwon(rho,nvert,xb,zb,nfld,xf,zf,pf)
        return
        end
c******************************************************************************
        subroutine talwon(rho,nvert,xb,zb,nobs,
     & xf,zf,gc)
c  gravity profile using the Talwani 2-d prism algorithm
c  xb,zb,xf,zf in kilometers, rho in gm/cc
c  reference Talwani et al, JGR, v64, no1.
c  coded by Mike Webring 9/80
        dimension xb(1),zb(1),xf(1),zf(1),gc(1)
        if(rho.eq.0.0) return
        const=2.0*6.673
        pi=3.1415927
        hfpi=pi*.5
        twopi=2.0*pi
        tol=1.e-5
c
        do 200 j=1,nobs
        gtmp=0.0
c
        do 101 i=1,nvert
        capz=0.0
        i1=i+1
        if(i.eq.nvert) i1=1
        x1=xb(i)-xf(j)
        x2=xb(i1)-xf(j)
        z1=zb(i)-zf(j)
        z2=zb(i1)-zf(j)
        ax1=abs(x1)
        ax2=abs(x2)
        if(ax1.lt.tol .and. abs(z1).lt.tol) go to 100
        if(ax2.lt.tol .and. abs(z2).lt.tol) go to 100
        dx=x2-x1
        dz=z2-z1
        t1=atan2(z1,x1)
        t2=atan2(z2,x2)
        if(t1.lt.0.0) t1=t1+twopi
        if(t2.lt.0.0) t2=t2+twopi
        dt=t1-t2
        if(abs(dt).gt.pi) dt=dt-sign(twopi,dt)
        if(abs(dt).lt.tol) go to 100
c
        if(abs(dx).gt.tol) go to 10
        ct2=cos(t2)
        if(abs(ct2).lt.tol) go to 100
        f1=cos(t1)/ct2
        capz=x1*log(f1)
        go to 100
c
10      if(abs(dz).gt.tol) go to 20
        capz=-dt*z1
        go to 100
c
20      tphi=dz/dx
        dxs=dx*dx
        r2=dxs+dz*dz
        if(r2.lt.tol) go to 100
        acs=(dx*dz*x2-dxs*z2)/r2
        if(ax1.gt.tol) go to 30
        f2=sin(t2)-cos(t2)*tphi
        capz=-acs*(-dt+tphi*log(f2))
        go to 100
c
30      if(ax2.gt.tol) go to 40
        f1=sin(t1)-cos(t1)*tphi
        if(f1.le.0.0) then
        print *,j,xf(j),zf(j),xb(i),zb(i),xb(i1),zb(i1)
        stop
        endif
        capz=acs*(dt+tphi*log(f1))
        go to 100
c
40      r1s=x1*x1+z1*z1
        r2s=x2*x2+z2*z2
        capz=acs*(dt+tphi*0.5*log(r2s/r1s))
        go to 100
c
100     gtmp=gtmp+capz
101     continue
200     gc(j)=gc(j)+const*rho*gtmp
        return
        end
c******************************************************************************
        subroutine grv2hd(rho,yin,nvert,xb,zb,nobs,
     & xf,zf,gc)
c  profile gravity calculation of a 2.5d prism
c  yin1 > yin2, z positive down
c  reference Rasmussen and Pedersen, Geophysical Prospecting, v27 p749-760.
c  coded by Mike Webring, 11/81
        dimension x(2),y(2),yin(2),z(2),xzr(2),yr(2,2),f(4,2),ef(8)
        dimension xb(nvert),zb(nvert),xf(nobs),zf(nobs),gc(nobs)
        logical sym
        equivalence (f(1,1),ef(1))
c  rho in gm/cc, distance in kilometer
        const=6.673*rho
        tol=1.e-5
        pi=3.1415926
        twopi=2.*pi
        sym=.false.
        if(yin(1).eq.-yin(2)) sym=.true.
        it=2
        if(sym) it=1
        add=sign(1.,yin(1)*yin(2))
        y(1)=abs(yin(1))
        y(2)=abs(yin(2))
        if(y(1).gt.y(2)) go to 1
        y(1)=abs(yin(2))
        y(2)=abs(yin(1))
1       if(y(1).lt.tol .or. rho.eq.0.0) return
        if(y(2).lt.tol) it=1
c
        do 200 j=1,nobs
        x(1)=xb(1)-xf(j)
        z(1)=zb(1)-zf(j)
        r2=x(1)*x(1)+z(1)*z(1)
        xzr(1)=sqrt(r2)
        yr(1,1)=sqrt(r2+y(1)*y(1))
        yr(1,2)=sqrt(r2+y(2)*y(2))
c
        do 10 k=1,8
10      ef(k)=0.
        do 100 i=1,nvert
        i1=i+1
        if(i.eq.nvert) i1=1
        x(2)=xb(i1)-xf(j)
        z(2)=zb(i1)-zf(j)
        r2=x(2)*x(2)+z(2)*z(2)
        xzr(2)=sqrt(r2)
        dx=x(2)-x(1)
        dz=z(2)-z(1)
        dxz2=dx*dx+dz*dz
        dl=sqrt(dxz2)
        if(dl.lt.tol) go to 90
        cosp=dx/dl
        zn=-cosp
        sinp=dz/dl
        w=cosp*z(1)-sinp*x(1)
        u1=cosp*x(1)+sinp*z(1)
        u2=cosp*x(2)+sinp*z(2)
c
        do 80 ip=1,it
        yr(2,ip)=sqrt(r2+y(ip)*y(ip))
        f(1,ip)=f(1,ip)+zn*y(ip)*alog((u2+yr(2,ip))/(u1+yr(1,ip)))
        f(2,ip)=f(2,ip)+zn*u2*alog((y(ip)+yr(2,ip))/xzr(2))
        f(3,ip)=f(3,ip)+zn*u1*alog((y(ip)+yr(1,ip))/xzr(1))
        xn=u2*y(ip)
        xd=w*yr(2,ip)
        t2=atan2(xn,xd)
        if(t2.lt.0.0) t2=t2+twopi
        xn2=u1*y(ip)
        xd2=w*yr(1,ip)
        t1=atan2(xn2,xd2)
        if(t1.lt.0.0) t1=t1+twopi
        dt=t2-t1
        if(abs(dt).gt.pi) dt=dt-sign(twopi,dt)
        f(4,ip)=f(4,ip)+zn*w*dt
80      continue
90      xzr(1)=xzr(2)
        yr(1,1)=yr(2,1)
        yr(1,2)=yr(2,2)
        x(1)=x(2)
        z(1)=z(2)
100     continue
c
        sum1=f(1,1)-f(4,1)+f(2,1)-f(3,1)
        sum2=sum1
        if(.not.sym) sum2=f(1,2)-f(4,2)+f(2,2)-f(3,2)
        gf=sum1-add*sum2
        gc(j)=gc(j)-const*gf
200     continue
        return
        end
c******************************************************************************
        subroutine mag2hd(vm,tf,y,nvert,xb,zb,nobs,
     & xf,zf,t)
c  computes the total field profile for a 2.5d prism
c  ref: shuey and pasquale, geophysics v38 n3
c  vm is the magnetization vector and tf is measurement unit vector
c  coded by mike webring, 10/81
        complex cdxz,cxz,ct1,ct2,cf,cfln
        dimension vm(3),tf(3),cf(2),x(2),z(2),r(2),
     1 xb(nvert),zb(nvert),xf(nobs),zf(nobs),t(nobs)
        if(y.gt.0.0) go to 1
        print *,' mag2hd: prism length <= zero, setting to 1000'
        y=1000.
1       ys=y*y
        do 50 iob=1,nobs
        px=0.
        pz=0.
        q=0.
        x(1)=xb(1)-xf(iob)
        z(1)=zb(1)-zf(iob)
        r(1)=sqrt(x(1)*x(1)+ys+z(1)*z(1))
c
        do 20 j=1,nvert
        j1=j+1
        if(j.eq.nvert) j1=1
        x(2)=xb(j1)-xf(iob)
        z(2)=zb(j1)-zf(iob)
          if(x(2).eq.0.0 .and. z(2).eq.0.0) then
          print *,' mag2hd: observation coincident with prism corner'
          return
        endif
        r(2)=sqrt(x(2)*x(2)+ys+z(2)*z(2))
        dx=x(2)-x(1)
        dz=z(2)-z(1)
        if(dx.eq.0.0 .and. dz.eq.0.0) then
          print *,' mag2hd: zero length segment, vertex ',j
          return
        endif
        cdxz=cmplx(dx,dz)
         do 10 k=1,2
         cxz=cmplx(x(k),z(k))
         t1=(1.0+r(k)/y)
         ct1=cmplx(t1,0.)
         t2=(x(k)*dz-z(k)*dx)/ys
         ct2=cmplx(0.,t2)
10       cf(k)=(cdxz*ct1)/cxz+ct2
        cfln=(clog(cf(2)/cf(1)))/cdxz
        px=px+dz*real(cfln)
        pz=pz-dx*aimag(cfln)
        q=q-dx*real(cfln)
        x(1)=x(2)
        z(1)=z(2)
20      r(1)=r(2)
c
        hx=2.0*(vm(1)*px+vm(3)*q)
        hy=2.0*vm(2)*(pz-px)
        hz=2.0*(vm(1)*q-vm(3)*pz)
50      t(iob)=t(iob)+hx*tf(1)+hy*tf(2)+hz*tf(3)
        return
        end
c******************************************************************************
	subroutine ckmodl(ifun,npt,xf,zf,ierr)
c  ifun=0 check both terrain and prisms
c  ifun=1 check terrain, ifun=2 check prism twists
c        common /parm1/ nvert,nbody,nobs,mobs,grvswt,magswt,dum(7)
c        common /model/ xb(200),zb(200),dum1(500)
c        common /vertx/icrn(500),ncrn(150),loca(150),nbord,ibord(200)
c        common /swork/ xf1(400),zf1(400),x(100),z(100),dum2(7420)
      common /model/nbod,ncor,iopn,delx,ibod,ifile
        common /model1/ xb(712),zb(712),rho(50),sus(50),
     1 pl(50,2),rem(50,3)
        common /vertx/ icrn(500),ncrn(150),loca(150)
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
      common /swork/ xf1(512),zf1(512)
      character prompt*79, gfile*80, ifile*80, sfile*80
        dimension lint(2),xf(npt),zf(npt)
        logical inside,member,planar
	ierr=0
        nn=loca(nbod)+ncrn(nbod)
	if(ifun.eq.2) go to 35
c
c   check for vertices above terrain (observations)
	do 20 i=1,npt
	xf1(i)=xf(i)
20	zf1(i)=zf(i)
	m=npt+2
	xf1(npt+1)=xf1(npt)
	xf1(m)=xf1(1)
	zf1(npt+1)=-100.
	zf1(m)=-100.
        do 30 i=1,ncor
	if(.not.member(i,icrn,nn)) go to 30
c        if(inside(m,xf1,zf1,xb(i),zb(i))) print  25,i
c25      format(' vertex',i3,' is above observations')
        if(inside(m,xf1,zf1,xb(i),zb(i))) then
          prompt='warning: vertices above observations'
          call tpromp(0,0,prompt)
        endif
30      continue
	if(ifun.eq.1) return
c
35      do 50 j=1,nbod
	if(planar(ncrn(j),xb,zb,icrn(loca(j)+1),lint)) go to 50
	ierr=1
          prompt='error: invalid body'
          call tpromp(0,0,prompt)
c        print *,' ckmodl>> body number ',j,' has an intersection'
c        print *,'       between line segments', lint(1), 'and', lint(2)
50	continue
c        if(planar(nbord,xb,zb,ibord,lint)) return
c        ierr=1
c        print *,' ckmodl>> the main model outer boundary has an'
c        print *,'          intersection between line segments'
c        print *,lint(1), 'and', lint(2)
	return
        end
c******************************************************************************
	logical function planar(n,x,z,lookup,lint)
c  check polygon for intersecting line segments
	dimension lint(2),lookup(n),x(1),z(1)
	lint(1)=0
	lint(2)=0
	planar=.true.
	if(n.eq.3) return
c        tolfac=.0001
        tolfac=.000001
c  first line segment
	do 10 line1=1,n-2
	i2=line1+1
	if(i2.gt.n) i2=1
	i=lookup(line1)
	i2=lookup(i2)
	dz=z(i2)-z(i)
	dx=x(i2)-x(i)
	adx=abs(dx)
	tol=tolfac*(adx+abs(dz)) + 1.e-10
	if(adx.lt.tol) dx=tol
	am=dz/dx
	amn=amin1(x(i),x(i2))-tol
	amx=amax1(x(i),x(i2))+tol
	aymn=amin1(z(i),z(i2))-tol
	aymx=amax1(z(i),z(i2))+tol
	k=n
	if(line1.eq.1) k=n-1
c  compared to segments: +2 to n segment-1
	do 10 line2=line1+2,k
	j2=line2+1
	if(j2.gt.n) j2=1
	j=lookup(line2)
	j2=lookup(j2)
	dz=z(j2)-z(j)
	dx=x(j2)-x(j)
	adx=abs(dx)
	tol=tolfac*(adx+abs(dz)) + 1.e-10
	if(adx.lt.tol) dx=tol
	bm=dz/dx
	dslope=am-bm
	if(abs(dslope).lt.tolfac) go to 10
	rs=1.0/dslope
	xi=(rs*am)*x(i)-(rs*bm)*x(j)+rs*(z(j)-z(i))
	yi=am*xi-am*x(i)+z(i)
	bmn=amin1(x(j),x(j2))-tol
	bmx=amax1(x(j),x(j2))+tol
	bymn=amin1(z(j),z(j2))-tol
	bymx=amax1(z(j),z(j2))+tol
c  test for intersection within segments
	if(xi.lt.amn .or. xi.lt.bmn) go to 10
	if(xi.gt.amx .or. xi.gt.bmx) go to 10
	if(yi.lt.aymn .or. yi.lt.bymn) go to 10
	if(yi.gt.aymx .or. yi.gt.bymx) go to 10
c  intersection between line segments
	lint(1)=line1
	lint(2)=line2
	planar=.false.
	go to 30
10	continue
30	return
        end
c******************************************************************************
	logical function member(itest,m,n)
	dimension m(n)
	member=.true.
	do 10 i=1,n
	if(m(i).eq.itest) return
10    continue
	member=.false.
	return
        end

