***********************************************************************
      program Perm2Mesh
***********************************************************************
c
c Inserts log-permeability field values from coord-log(k) file as
c permeability modifiers in columns 41-50 of TOUGH2 mesh file.
c
c S. Finsterle
c
      parameter (maxel=150000,maxihist=1000)

      character*5  elem,cmat
      character*10 caht
      character*40 filein,fileout,fileperm,filehist
      character*80 line

      dimension elem(maxel),xx(maxel),yy(maxel),zz(maxel)
      dimension permmod(maxel),ihist(maxihist)
      dimension readin(4)
     
      data maxhist/1/,minhist/maxihist/

      write(*,*)
      write(*,*) ' Add Perms to MESH'
      write(*,*) ' *****************'
      write(*,*)
      write(*,*) ' Input file with coord-log(k) values: ?  '
      read(*,'(a)')fileperm
      open(3,file=fileperm, status='old')
      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(*,*) ' Number of header lines             : ?  '
      read(*,*)    nheader
      write(*,*) ' Dimension of log(k) field          : ?  '
      write(*,*) ' 2 - x y log(k)'
      write(*,*) ' 3 - x y z log(k)'
      read(*,*)    ndim   
      write(*,*) ' TOUGH2 mesh is                     : ?  '
      write(*,*) ' 0 - xy '
      write(*,*) ' 1 - x z'
      write(*,*) ' 2 -  yz'
      write(*,*) ' 3 - xyz'
      read(*,*)    ndimt2 
      write(*,*) ' Replace, add, or multiply          : ?  '
      write(*,*) ' 1 - replace'
      write(*,*) ' 2 - add'
      write(*,*) ' 3 - multiply'
      read(*,*)    nram
      write(*,*) ' Output histogram file              : ?  '
      read(*,'(a)')filehist
      open(4,file=filehist, status='unknown')
      write(*,*) ' Histogram class size               : ?  '
      read(*,*)    class
      write(*,*)
c
c --- read file with coord-log(k) values
      do 1000 i=1,nheader
         read(3,7000,end=9000) line
 1000 continue
      iel=0             
 1500 continue
      iel=iel+1
      read(3,*,end=1600) (readin(i),i=1,ndim+1)        
c
c --- Assign coordinates
      if (ndimt2.eq.1) then
c
c --- TOUGH2 mesh is 2d-xz
         xx(iel)=readin(1)
         yy(iel)=0.0
         if (ndim.eq.2) then
            zz(iel)=readin(2)
         else
            zz(iel)=readin(3)
         endif
      else if (ndimt2.eq.2) then
c
c --- TOUGH2 mesh is 2d-yz
         xx(iel)=0.0
         if (ndim.eq.2) then
            yy(iel)=readin(1)
            zz(iel)=readin(2)
         else
            yy(iel)=readin(2)
            zz(iel)=readin(3)
         endif
      else
c
c --- TOUGH2 mesh is 2d-xy or 3d
         xx(iel)=readin(1)
         yy(iel)=readin(2)
         if (ndimt2.eq.3) then
            zz(iel)=readin(3)
         else
            zz(iel)=0.0
         endif
      endif
c
c --- convert log(k) to permeability modifier
      permmod(iel)=-10**readin(ndim+1)
      goto 1500
 1600 continue
      nel=iel-1
 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,permin,x,y,z
 7001 format(a5,10x,a5,e10.4,a10,4e10.4)
      if (elem(iel).eq.'    ') goto 2002
      distmin=1.0e10
      do 1001 i=1,nel
         dist=(xx(i)-x)*(xx(i)-x)+
     &        (yy(i)-y)*(yy(i)-y)+
     &        (zz(i)-z)*(zz(i)-z)
         if (dist.lt.distmin) then
            distmin=dist
            ith=i
         endif      
 1001 continue
      read(caht,'(e10.4)',iostat=ios) aht
      if (abs(aht).lt.1.0e-10) caht='          '
      if (abs(permin).lt.1.0e-10) permin=-1.0
      if (nram.eq.1) then
         permfact=-abs(permmod(ith))
      else if (nram.eq.2) then
         permfact=-abs(permin)-abs(permmod(ith))
      else if (nram.eq.3) then
         permfact=-abs(permin*permmod(ith))
      endif
      write(2,7001) elem(iel),cmat,vol,caht,permfact,x,y,z
      nhist=int(log10(-permfact)/class)+maxihist/2
      maxhist=max(maxhist,nhist)
      minhist=min(minhist,nhist)
      ihist(nhist)=ihist(nhist)+1
      goto 2001
 2002 continue
c
c --- copy rest of file
      write(2,7000) '     '
 2007 continue
      read(1,7000,end=9000) line
      call lenos(line,ll)
      write(2,7000) line(:ll)
      goto 2007
 9000 continue
      write(4,7002)
 7002 format('variables = "log(permfact)" n',/,'zone')
      do 2008 i=minhist,maxhist
         write(4,'(f10.4,i6)') class*(i-maxihist/2)+0.5*class,ihist(i)
 2008 continue
      close(1)
      close(2)
      close(3)
      close(4)
      end

c --- end of Perm2Mesh


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