***********************************************************************
      program CutCavity
***********************************************************************
c
c Cuts cavity from TOUGH2 mesh.      
c Version 1.0
c
c S. Finsterle, February 9, 2000
c 
      parameter (maxel=150000)
     
      character*5  elem,elem1,elem2,cavit,cmat,cdum
      character*10 caht,cpermmod
      character*40 filein,fileout
      character*80 line

      dimension elem(maxel)
      dimension iniche(maxel)

      data cavit/'NIC99'/
      data iniche/maxel*0/

      write(*,*)
      write(*,*) ' Cut Ellipsoidal Cavity'
      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(*,*) ' Volume of cavity element           : ?  '
      read(*,*)    volume
      write(*,*) ' Nodal distance from cavity wall    : ?  '
      read(*,*)    distance
      write(*,*) ' Cosine multiplication factor       : ?  '
      read(*,*)    cosine
      write(*,*) ' Xcenter                            : ?  '
      read(*,*)    xcenter
      write(*,*) ' Ycenter                            : ?  '
      read(*,*)    ycenter
      write(*,*) ' Zcenter                            : ?  '
      read(*,*)    zcenter
      write(*,*) ' X semiaxis                         : ?  '
      read(*,*)    a     
      write(*,*) ' Y semiaxis                         : ?  '
      read(*,*)    b     
      write(*,*) ' Z semiaxis                         : ?  '
      read(*,*)    c     
      open(3,file='ellipsoid.conne',status='unknown')

 1000 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 1000
      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 cavity(x,y,z,xcenter,ycenter,zcenter,a,b,c,iniche(iel))
      if (iniche(iel).le.0) then
         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
      iel=iel-1
      cavit(4:5)='97'
      write(2,7001) cavit,'CAVIT',volume/3.0,' ',' ',0.0,0.0,zcenter
      cavit(4:5)='98'
      write(2,7001) cavit,'CAVIT',volume/3.0,' ',' ',0.0,0.0,zcenter
      cavit(4:5)='99'
      write(2,7001) cavit,'CAVIT',volume/3.0,' ',' ',0.0,0.0,zcenter
c
c --- connections
      write(2,7000) '     '
      write(2,7000) 'CONNE'
 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.'+++  ') then
         write(2,7000) elem1
         goto 2007
      endif
      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 (iniche(i).eq.1) then
         elem1=cavit
         d1=distance
         betax=betax*cosine
      endif
      if (iniche(j).eq.1) then
         elem2=cavit
         d2=distance
         betax=betax*cosine
      endif
      if (iniche(i)+iniche(j).lt.2) then
         if (iniche(j).eq.1) then
            cdum=elem1
            elem1=elem2
            elem2=cdum
            d=d1
            d1=d2
            d2=d
            betax=-betax
         endif
         if (iniche(i).eq.1.or.iniche(j).eq.1) then
            if (isot.eq.1) then
               elem1(4:5)='97'
            else if (isot.eq.2) then
               elem1(4:5)='98'
            else
               elem1(4:5)='99'
            endif
            cdum=elem2
            if (elem2(4:4).eq.' ') cdum(4:4)='_'
            write(3,7003)  elem1,cdum
 7003 format(9x,a5,1x,a5,2x'&')
         endif
         write(2,7002) elem1,elem2,isot,d1,d2,areax,betax
      endif
      goto 2003
c
c --- write rest of file
 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)
      close(3)
      end

c --- end of CutNiche


***********************************************************************
      subroutine cavity(x,y,z,xcenter,ycenter,zcenter,a,b,c,ini)
***********************************************************************

      ini=0
c 
c --- cut out ellipsoidal cavity
      dist=(x-xcenter)**2/a**2+
     &    (y-ycenter)**2/b**2+
     &    (z-zcenter)**2/c**2
      if (dist.le.1.0) then
         ini=1
         return
      endif
      end

c --- end of cavity


***********************************************************************
      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