***********************************************************************
      program AddBound
***********************************************************************
c
c replaces all elements in a specified box or cylinder by a single element
c or removes them if element volume is zero
c
c S. Finsterle
c
      parameter (maxel=150000)
     
      character*5  cbound,elem,elem1,elem2,cmat,cbmat,cmaterial
      character*10 caht,cpermmod
      character*40 filein,fileout
      character*80 line

      dimension elem(maxel)
      dimension ibound(maxel)

      data ibound/maxel*0/

      write(*,*)
      write(*,*) ' Add Boundary Element'
      write(*,*) ' ********************'
      write(*,*)
      write(*,*) ' Input MESH file                    : ?  '
      read(*,'(a)')filein
      open(1,file=filein, status='old')
      write(*,*) ' Output MESH file                   : ?  '
      read(*,'(a)')fileout
      open(2,file=fileout,status='unknown')
      write(*,*) ' Boundary element name (AAAII)      : ?  '
      read(*,'(a)')cbound
      write(*,*) ' Boundary rock type name (A5)       : ?  '
      read(*,'(a)')cbmat
      write(*,*) ' Boundary element volume            : ?  '
      read(*,*)    volume
      write(*,*) ' Nodal distance to boundary element : ?  '
      read(*,*)    distance
      write(*,*) ' Domain shape: 1=cube                    '
      write(*,*) '               2=cylinder           : ?  '
      write(*,*) '               3=material name      : ?  '
      read(*,*)    ishape
      if (ishape.eq.3) then
         write(*,*) ' Material name                      : ?  '
         read(*,*)    cmaterial
      else
         write(*,*) ' Xmin                               : ?  '
         read(*,*)    xmin
         write(*,*) ' Xmax                               : ?  '
         read(*,*)    xmax
         write(*,*) ' Ymin                               : ?  '
         read(*,*)    ymin
         write(*,*) ' Ymax                               : ?  '
         read(*,*)    ymax
         write(*,*) ' Zmin                               : ?  '
         read(*,*)    zmin
         write(*,*) ' Zmax                               : ?  '
         read(*,*)    zmax
         if (ishape.eq.2) then
            write(*,*) ' Radius                             : ?  '
            read(*,*)    radius
         endif
      endif
      xave=(xmin+xmax)/2.
      yave=(ymin+ymax)/2.
      zave=(zmin+zmax)/2.
      write(*,*)
      if (volume.ne.0.0) then
         ibm=2
      else
         ibm=1
      endif
 2000 continue
      read(1,7000,end=9000) line
 7000 format(a)
      call lenos(line,ll)
      write(2,7000) line(:ll)
      if (line(1:5).ne.'ELEME') goto 2000
      iel=0
 2001 continue
c
c --- elements
      iel=iel+1
      if (mod(iel,1000).eq.1)
     &    write(*,*) ' Working on element ',max(1,iel-1)
      read(1,7001) elem(iel),cmat,vol,caht,cpermmod,x,y,z
 7001 format(a5,10x,a5,e10.4,2a10,3e10.4)
      if (elem(iel).eq.'    ') goto 2002
      call select(ishape,x,y,z,xmin,xmax,ymin,ymax,zmin,zmax,
     &            radius,cmat,cmaterial,ibound(iel))
      if (ibound(iel).eq.0) then
         read(cpermmod,'(e10.4)',iostat=ios) permmod
         if (abs(permmod).lt.1.0e-10) cpermmod='          '
         read(caht,'(e10.4)',iostat=ios) aht
         if (abs(aht).lt.1.0e-10) caht='          '
         write(2,7001) elem(iel),cmat,vol,caht,cpermmod,x,y,z
      endif
      goto 2001
 2002 continue
      if (abs(volume).gt.1.0e-10) then
         caht='          '
         cpermmod='          '
         write(2,7001) cbound,cbmat,volume,caht,cpermmod,
     &                 xave,yave,zave
      endif
      iel=iel-1
c
c --- connections
      read(1,7000) elem1
      write(2,7000) '     '
      write(2,7000) 'CONNE'
      icon=0
 2003 continue
      icon=icon+1
      if (mod(icon,1000).eq.1)
     &    write(*,*) ' Working on connection ',max(1,icon-1)
      read(1,7002) elem1,elem2,isot,d1,d2,areax,betax
 7002 format(2a5,15x,i5,4e10.4)
      if (elem1.eq.'     '.or.elem1.eq.'+++  ') goto 2006
      do 1001 i=1,iel
         if (elem1.eq.elem(i)) goto 2004
 1001 continue
 2004 continue
      do 1002 j=1,iel
         if (elem2.eq.elem(j)) goto 2005
 1002 continue
 2005 continue
      if (i.gt.iel.or.j.gt.iel) goto 2003
      if (ibound(i).eq.1) then
         elem1=cbound
         d1=distance
      endif
      if (ibound(j).eq.1) then
         elem2=cbound
         d2=distance
      endif
      if (ibound(i)+ibound(j).lt.ibm) then
         write(2,7002) elem1,elem2,isot,d1,d2,areax,betax
      endif
      goto 2003
 2006 continue
      write(2,7000) elem1
 2007 continue
      read(1,7000,end=9000) line
      call lenos(line,ll)
      write(2,7000) line(:ll)
      goto 2007
 9000 continue
      close(1)
      close(2)
      end

c --- end of AddBound

***********************************************************************
      subroutine select(ishape,x,y,z,xmin,xmax,ymin,ymax,zmin,zmax,
     &                  radius,cmat,cmaterial,ini)                    
***********************************************************************

      character*5 cmat,cmaterial

      ini=0
      if (ishape.eq.1) then
         if (x.ge.xmin.and.x.le.xmax.and.
     &       y.ge.ymin.and.y.le.ymax.and.
     &       z.ge.zmin.and.z.le.zmax) ini=1
      else if (ishape.eq.2) then
         dotp=(xmin-x)*(xmax-xmin)+
     &        (ymin-y)*(ymax-ymin)+
     &        (zmin-z)*(zmax-zmin)
         x2x1=(xmax-xmin)*(xmax-xmin)+
     &        (ymax-ymin)*(ymax-ymin)+
     &        (zmax-zmin)*(zmax-zmin)
         t=-dotp/x2x1
         if (t.ge.0. .and. t.le. 1.0) then
            x1x=(xmin-x)*(xmin-x)+
     &          (ymin-y)*(ymin-y)+
     &          (zmin-z)*(zmin-z)
            dist=(x1x*x2x1-dotp*dotp)/x2x1
            if (dist.le.radius*radius) ini=1
         endif
      else if (ishape.eq.3) then
         if (cmat.eq.cmaterial) ini=1
      endif
      end

c --- end of select


***********************************************************************
      subroutine lenos(string,lstr)
***********************************************************************
      character string*(*)

      lstr=1
      ilenstr=len(string)
      if (ilenstr.eq.0) return
      do 1000 i=ilenstr,1,-1
        if (string(i:i).ne.char(9).and.string(i:i).ne.' ') then
          lstr=i
          return
        endif
1000  continue
      end

c --- end of lenos