program meshadd2

     integer, parameter:: np = 10001
     integer :: i,max
     real(8), dimension(0:np) :: xmesh1, ymesh1
     real(8), dimension(0:np) :: xmesh2, ymesh2
     real(8), dimension(np) :: x, y
     character(len=20) file1,file2
! ------------------------------------------------------ !
     write(6,*) 'First infile'
     read(5,*) file1
     open(10,file=file1)
     i=1
     do 
        read(10,*,end=10) xmesh1(i), ymesh1(i)
        i=i+1
     enddo
10   continue
     close(10)
     i1max=i
     print*,'i1max=',i1max
! ------------------------------------------------------ !
     write(6,*) 'Second infile'
     read(5,*) file2
     open(11,file=file2)
     i=1
     do 
        read(11,*,end=11) xmesh2(i), ymesh2(i)
        i=i+1
     enddo
11   continue
     close(11)
     i2max=i
     print*,'i2max=',i2max

     call grid_interleave(xmesh1,ymesh1,i1max,xmesh2,ymesh2,i2max,x,y,max)
     print*,'max=',max
     print*,'x(1)=',x(1)
     do i=1,max
        write(14,*) x(i),' ',y(i)
     enddo

    stop 'Code Complete'
    end program meshadd2
! ----------------------------------------------------------------------- !
     subroutine grid_interleave(xmesh1,ymesh1,i1max,xmesh2,ymesh2,i2max,&
     & x,y,nsz)
! ----------------------------------------------------------------------- !
     integer, parameter:: np = 10001
     integer :: i,max
     integer :: i1  = 1
     integer :: i2  = 1
     integer :: i1max
     integer :: i2max
     integer :: ip1 = 1
     integer :: ip2 = 1
     integer :: flg = 0
     integer :: ict 
     integer :: nsz
     integer :: ixm ,maxixm
     integer :: ixpt,maxixpt 
     real(8) :: eps = 1.0e-10
     real(8) :: diff
     integer, dimension(np) :: ilower1
     integer, dimension(np) :: ilower2
     real(8), dimension(np) :: x,y 
     real(8), dimension(0:np) :: x1, y1, x2, y2
     real(8), dimension(0:np) :: xmesh1, ymesh1
     real(8), dimension(0:np) :: xmesh2, ymesh2
     real(8), dimension(np) :: xpts1, ypts1 
     real(8), dimension(np) :: xpts2, ypts2 
     real(8), external :: func
     do 
        diff=abs(xmesh1(i1)-xmesh2(i2))
        if(xmesh1(i1).gt.xmesh2(i2))then
            xpts1(ip1)=xmesh2(i2)
            ilower1(ip1)=i1
            ip1=ip1+1
          if(i2.le.i2max)then
            i2=i2+1
          else
            exit
          endif
          flg=1
        elseif(xmesh1(i1).lt.xmesh2(i2))then
            xpts2(ip2)=xmesh1(i1)
            ilower2(ip2)=i2
            ip2=ip2+1
          if(i1.le.i1max)then
            i1=i1+1
          else
            exit            
          endif
        elseif(diff.le.eps)then
          if(i1.le.i1max)then
            i1=i1+1
          else
            exit            
          endif
          if(i2.le.i2max)then
            i2=i2+1
          else
            exit 
          endif
        else
          Print*,'ERROR ELSE TRAP'
        endif
     enddo 
     
     do i=1,ip1-1
        ypts1(i)=func(ymesh1(0),xmesh1(0),xpts1(i),ilower1(i))
        print*,'ypts1(',i,')=',ypts1(i)
     enddo
     do i=1,ip2-1
        ypts2(i)=func(ymesh2(0),xmesh2(0),xpts2(i),ilower2(i))
        print*,'ypts2(',i,')=',ypts2(i)
     enddo

     do i=1,ip1-1
        write(40,*) xpts1(i), ypts1(i)
     enddo
     do i=1,ip2-1
        write(41,*) xpts2(i), ypts2(i)
     enddo
     call combine(np,i1max,ip1-1,xmesh1,ymesh1,xpts1,ypts1,x1,y1)
     max1=aint(x1(0))
     do i=1,max1
        write(10,*) x1(i),' ',y1(i)
     enddo
     call combine(np,i2max,ip2-1,xmesh2,ymesh2,xpts2,ypts2,x2,y2)
     max2=aint(x2(0))
     if(max2.ne.max1) print*,'ERROR: xmesh2 .ne. xmesh1'
     do i=1,max2
        write(11,*) x2(i),' ',y2(i)
     enddo
     do i=1,max2
        x(i)=x1(i)
        y(i)=y1(i)+y2(i)
        write(12,*) x1(i),y(i)
     enddo
     nsz=x2(0)
     end subroutine grid_interleave

! -------------------------------------------------------------------------
     subroutine combine(np,maxixm,maxixpt,xmesh,ymesh,xpts,ypts,x,y)
! -------------------------------------------------------------------------
     integer :: ict
     integer :: ixm ,maxixm
     integer :: ixpt,maxixpt
     real(8), dimension(0:np) :: xmesh, ymesh
     real(8), dimension(np)   :: xpts, ypts
     real(8), dimension(0:np) :: x, y 
     ict  = 0
     ixm  = 1
     ixpt = 1 
     do 
     ict = ict+1
     if(xmesh(ixm).lt.xpts(ixpt)) then
       x(ict)=xmesh(ixm)
       y(ict)=ymesh(ixm)
       ixm=ixm+1
       if(ixm.gt.maxixm) goto 200
     elseif(xmesh(ixm).gt.xpts(ixpt)) then
       x(ict)=xpts(ixpt)
       y(ict)=ypts(ixpt)
       ixpt = ixpt+1
       if(ixpt.gt.maxixpt)then
          do           ! copy the rest of xmesh into x
            ict=ict+1
            x(ict)=xmesh(ixm)
            y(ict)=ymesh(ixm) 
            ixm=ixm+1 
            if(ixm.gt.maxixm) goto 200
          enddo
       endif
     endif
     enddo
200  continue
     x(0)=ict-1
     y(0)=ict-1
     return
     end subroutine combine 
! -------------------------------------------------------
      double precision function func(y,x,x0,ilower)
! -------------------------------------------------------
! Region: Profile
! Purpose: Linear Interpolation
!          Using bracket to interpolate at grid point x0
     
      integer,parameter :: np=10001
      integer :: iupper, ilower
      real(8),dimension(np) :: y
      real(8),dimension(np) :: x
      real(8) :: x0
      real(8) :: slope
      
      iupper=ilower+1
      slope=(y(iupper)-y(ilower))/(x(iupper)-x(ilower))
      func=y(iupper)+slope*(x0-x(iupper))
      end function func