merge.f

Go to the documentation of this file.

C> @file grid2D/merge.f
C> Merge two 2D grids   
C GRIDMAN, grid managment library. Author: Vladislav Kotov, v.kotov@fz-juelich.de

!    Copyright (c) 2017  Forschungszentrum Juelich GmbH
!                        Vladislav Kotov 
!
!    This file is part of GRIDMAN.
!
!    GRIDMAN is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    GRIDMAN is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with GRIDMAN.  If not, see <http://www.gnu.org/licenses/>.

C> Merge two 2D grids by connecting their boundary edges
C>
C> GRID2 is merged into GRID1 by connecting their boundary edges.
C>
C> WARNING: the object GRID must not be equal to GRID1 or GRID2!
C>
C> WARNING: The object GRID will be overwritten if already exists!
C>
C> WARNING: the algorithm implemented at the moment assumes 
C> that edges of GRID1 and GRID2 which overlap with each over 
C> can be completely removed. That is, there are no 
C> "hanging tails" which do not belong to the intersection 
C>
C> Units of GRID are same as units of GRID1
C
      SUBROUTINE gridman_grid2d_merge(GRID,GRID1,GRID2,TOL,IERR)
      USE gridman,ONLY:gridman_grid,gridman_sp,gridman_dp,
     u                 gridman_unit,gridman_dbg,gridman_check
      USE gridman_lib,ONLY:gridman_grid_allocate,
     u                     gridman_index_transform, gridman_index_copy
      IMPLICIT NONE
      INTRINSIC sqrt,huge,max,trim,tiny
C> Resulting grid object
      TYPE(gridman_grid) :: GRID  
C> First grid which stays intact 
      TYPE(gridman_grid) :: GRID1 
C> Second grid which is merged into the first one
      TYPE(gridman_grid) :: GRID2       
C> Relative accuracy of the points coordinates
C>
C> - If TOL is too large than the code can mistakenly take two different points for one point
C> - If TOL is too small then the code cannot recognise a point which is same in both grids 
C> 
C> Recommended value TOL=1e-5.
C> It might be necessary to increase or decrease TOL in each individual case.
      REAL(GRIDMAN_DP),INTENT(IN) :: TOL 
C> Error code     
      INTEGER,INTENT(OUT) :: IERR
     
      INTEGER(GRIDMAN_SP) :: 
     i           iedge,iedge1,iedge2,nedges,enedges,nedges0,
     i           p(2),q(2),npoints,ipoint,ind(6),
     i           ipoint1,ipoint2, npoints1,ncells,tpq(2,2),              
     i           points(2,2*grid1%NEDGES+2*grid2%NEDGES),!TEMPORARY ARRAY FOR POINT INDICES
     i           cells(2,2*grid1%NEDGES+2*grid2%NEDGES), !TEMPORARY ARRAY FOR CELL INDICES
     i           mapp(grid2%NPOINTS),      !INDEX OF POINT IN COMBINED GRID FOR EACH POINT IN GRID2
     i           rpoints(grid2%NPOINTS)
      INTEGER(GRIDMAN_SP) ::
     i           imape1(2,grid1%NEDGES+grid2%NEDGES),   !MAPPING OF THE EDGES OF GRID TO GRID1  
     i           imape2(2,grid1%NEDGES+grid2%NEDGES),   !MAPPING OF THE EDGES OF GRID TO GRID2  
     i           lmap0(max(grid1%NEDGES,grid2%NEDGES))    !MAPPING OF GRID1,2 TO GRID  
      INTEGER :: RES,IERR0,NEDGEINDEX,NCELLINDEX,II,II2

      LOGICAL MEDGE1(grid1%nedges),MEDGE2(grid2%nedges)

      REAL(GRIDMAN_DP) :: LE1(grid1%nedges),    !LENGTH COVERED BY THE INTERSECTION FOR EACH EDGE:
     r                    le2(grid2%NEDGES)     !TO CHECK IF INTERSECTION COVERS THE WHOLE EDGE
      REAL(GRIDMAN_DP) :: XP(2),YP(2),XQ(2),YQ(2),
     r                    x1,y1,x2,y2,dlp,dlq,fscl,dlt,eps
    
      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_MERGE"

      ierr=0
     
      IF(grid1%TYPE.NE.2.OR.grid1%PDIM.NE.2.OR.grid1%EDIM.NE.2) THEN
       ierr=100
       WRITE(gridman_unit,*) "GRIDMAN_GRID2D_TRIANG: ",
     w                       "1st grid is not of 2D type"
       WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
     w                         grid1%TYPE,grid1%PDIM,grid1%EDIM
       RETURN 
      END IF 

      IF(grid2%TYPE.NE.2.OR.grid2%PDIM.NE.2.OR.grid2%EDIM.NE.2) THEN
       ierr=100
       WRITE(gridman_unit,*) "GRIDMAN_GRID2D_TRIANG: ",
     w                       "1st grid is not of 2D type"
       WRITE(gridman_unit,*) " TYPE, PDIM, EDIM ",
     w                         grid2%TYPE,grid2%PDIM,grid2%EDIM
       RETURN 
      END IF 

      IF(gridman_check) THEN
       CALL gridman_grid2d_check(grid1,res,ierr)
       IF(res.NE.0.OR.ierr.GT.0) THEN
        ierr=100
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                        "incorrect GRID1"
        RETURN
       END IF
       CALL gridman_grid2d_check(grid2,res,ierr)
       IF(res.NE.0.OR.ierr.GT.0) THEN
        ierr=100
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                        "incorrect GRID2"
        RETURN
       END IF
      END IF

      eps=10.*tiny(eps) !MACHINE PRECISION
C     FACTOR TO RESCALE UNITS OF GRID2 INTO UNITS OF GRID1
      IF(grid1%UNIT2SI.LT.eps) THEN
       ierr=100
       WRITE(gridman_unit,*) "ERROR detected in GRIDMAN_GRID2D_MERGE: ",
     w                       "wrong unit scale"
       WRITE(gridman_unit,*) " GRID1%UNIT2SI ",grid1%UNIT2SI
       RETURN
      END IF
      fscl=grid2%UNIT2SI/grid1%UNIT2SI

C
C     COLLECT INTERSECTING EDGES AND THEIR NON-INTERSECTING PARTS
C      
      rpoints=0 !NEW INDICES OF TNE POINTS IN GRID2
      npoints=grid1%NPOINTS
      ncells=grid1%NCELLS
      nedges=0
      enedges=0     
      medge1=.true.
      medge2=.true.
      le1=0.0
      le2=0.0
      DO iedge1=1,grid1%NEDGES
       IF(grid1%CELLS(1,iedge1).GT.0.AND.
     f    grid1%CELLS(2,iedge1).GT.0) cycle !SKIP NON-BOUNDARY EDGES
C      TAKE EDGE OF GRID1
       p(1)=grid1%POINTS(1,iedge1)
       p(2)=grid1%POINTS(2,iedge1)
       xp(1)=grid1%X(1,p(1))
       yp(1)=grid1%X(2,p(1))
       xp(2)=grid1%X(1,p(2))
       yp(2)=grid1%X(2,p(2))
       DO iedge2=1,grid2%NEDGES
       IF(grid2%CELLS(1,iedge2).GT.0.AND.
     f    grid2%CELLS(2,iedge2).GT.0) cycle
C       TAKE EDGE OF GRID2
        q(1)=grid2%POINTS(1,iedge2)
        q(2)=grid2%POINTS(2,iedge2)
        xq(1)=grid2%X(1,q(1))*fscl
        yq(1)=grid2%X(2,q(1))*fscl
        xq(2)=grid2%X(1,q(2))*fscl
        yq(2)=grid2%X(2,q(2))*fscl
C       SHIFT FOR THE INDICES OF THE POINTS OF GRID2
        q(1)=q(1)+npoints 
        q(2)=q(2)+npoints
C       CHECK IF TWO EDGES OVERLAP
        CALL find_overlap(p,q,xp,yp,xq,yq,tol,ind,tpq)
C
C       IND(1) -- HANGING -- IND(2)-IND(3) -- OVERLAP --IND(4)-IND(5) -- HANGING  -- IND(6)
C
C       IND(3) ALWAYS CORRESPONDS TO Q(1), IND(4) ALWAYS CORRESPONDS TO Q(2) 
        IF(ind(3).GT.0.AND.ind(4).GT.0) THEN
C        OVERLAP FOUND, ADD AS NEW EDGE
         nedges=nedges+1
         points(1,nedges)=ind(3)
         points(2,nedges)=ind(4)
         cells(1,nedges)=nonzer_cell(grid1%CELLS(:,iedge1))
         cells(2,nedges)=nonzer_cell(grid2%CELLS(:,iedge2))
         IF(cells(2,nedges).GT.0) cells(2,nedges)=cells(2,nedges)+ncells
         medge1(iedge1)=.false. !EDGES WHICH MUST BE EXCLUDED 
         medge2(iedge2)=.false.
C       FOR CHECKS
C        LENGTHS OF THE OVERLAP
         ipoint1=ind(3)
         IF(ipoint1.GT.npoints) THEN          
          x1=grid2%X(1,ipoint1-npoints)*fscl
          y1=grid2%X(2,ipoint1-npoints)*fscl
         ELSE
          x1=grid1%X(1,ipoint1)
          y1=grid1%X(2,ipoint1)
         END IF
         ipoint2=ind(4)
         IF(ipoint2.GT.npoints) THEN          
          x2=grid2%X(1,ipoint2-npoints)*fscl
          y2=grid2%X(2,ipoint2-npoints)*fscl
         ELSE
          x2=grid1%X(1,ipoint2)
          y2=grid1%X(2,ipoint2)
         END IF
         dlp=sqrt((x1-x2)**2+(y1-y2)**2)
         le1(iedge1)=le1(iedge1)+dlp
         le2(iedge2)=le2(iedge2)+dlp
C        FOR MAPPING
         imape1(1,nedges)=iedge1
         imape1(2,nedges)=nedges
         imape2(1,nedges)=iedge2
         imape2(2,nedges)=nedges
        END IF !IF(IND(3).GT.0.AND.IND(4).GT.0)
C       REPLACE POINTS IN GRID2 BY POINTS OF GRID1, IF NECESSARY
        ipoint1=tpq(2,1)
        ipoint2=tpq(2,2)
        IF(ipoint1.GT.0) rpoints(tpq(1,1)-npoints)=ipoint1
        IF(ipoint2.GT.0) rpoints(tpq(1,2)-npoints)=ipoint2 
C      "HANGING EDGES" ARE NOT SAVED IN THE CURRENT IMPLEMENTATION
       END DO !DO IEDGE2=1,GRID2%NEDGES
      END DO !DO IEDGE1=1,GRID1%NEDGES
  
      IF(nedges.EQ.0) THEN
       WRITE(gridman_unit,*) "WARNING from GRIDMAN_GRID2D_MERGE:",        
     w                      " grids are not attached (no common edges)"    
       ierr0=-400 
      ELSE 
       ierr0=0
      END IF

C     CHECK THAT OVERLAPED POINTS ARE CORRECT
      DO ipoint2=1,grid2%NPOINTS
       ipoint1=rpoints(ipoint2)
       IF(ipoint1.GT.0) THEN
        x1=grid1%X(1,ipoint1)
        x2=grid2%X(1,ipoint2)*fscl
        y1=grid1%X(2,ipoint1)
        y2=grid2%X(2,ipoint2)*fscl
        dlp=(x1-x2)**2+(y1-y2)**2
        dlt=abs(x2-x1)*(abs(x1)+abs(x2)+eps)+
     +      abs(y2-y1)*(abs(y1)+abs(y2)+eps)
        dlt=2.*tol*(dlt+eps)
        IF(dlp.GT.dlt) THEN
          ierr=400
          WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                          "internal error"
          WRITE(gridman_unit,*) 
     w          "points marked as overlapped are not same, "
          WRITE(gridman_unit,*) 
     w          " IPOINT1, IPOINT2, DLP, DLT, TOL", 
     w            ipoint1, ipoint2, dlp, dlt, tol    
          WRITE(gridman_unit,*) " X1, Y1 ",x1,y1
          WRITE(gridman_unit,*) " X2, Y2 ",x2,y2
          RETURN
        END IF
       END IF
      END DO

C     CHECK THAT INTERSECTIONS COVER THE WHOLE EDGES
      DO iedge1=1,grid1%NEDGES    
       IF(.NOT.medge1(iedge1)) THEN
        ipoint=grid1%POINTS(1,iedge1)
        x1=grid1%X(1,ipoint)
        y1=grid1%X(2,ipoint)
        ipoint=grid1%POINTS(2,iedge1)
        x2=grid1%X(1,ipoint)
        y2=grid1%X(2,ipoint)
        dlp=(x1-x2)**2+(y1-y2)**2
        dlq=le1(iedge1)**2
        dlt=abs(x2-x1)*(abs(x1)+abs(x2)+eps)+
     +      abs(y2-y1)*(abs(y1)+abs(y2)+eps)
        dlt=4.*tol*(dlt+eps)
C       IF LENGTH OF THE INTERSECTION IS NOT EQUAL TO THE LENGTH OF THE EDGE THEN ERROR
        IF(abs(dlp-dlq).GT.dlt) GOTO 30
       END IF
      END DO

      DO iedge2=1,grid2%NEDGES 
       IF(.NOT.medge2(iedge2)) THEN
        ipoint=grid2%POINTS(1,iedge2)
        x1=grid2%X(1,ipoint)*fscl
        y1=grid2%X(2,ipoint)*fscl
        ipoint=grid2%POINTS(2,iedge2)
        x2=grid2%X(1,ipoint)*fscl
        y2=grid2%X(2,ipoint)*fscl
        dlp=(x1-x2)**2+(y1-y2)**2
        dlq=le2(iedge2)**2
        dlt=abs(x2-x1)*(abs(x1)+abs(x2)+eps)+
     +      abs(y2-y1)*(abs(y1)+abs(y2)+eps)
        dlt=4.*tol*(dlt+eps)
C       IF LENGTH OF THE INTERSECTION IS NOT EQUAL TO THE LENGTH OF THE EDGE THEN ERROR
        IF(abs(dlp-dlq).GT.dlt) GOTO 30
       END IF
      END DO

      GOTO 31
 30   ierr=400
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "internal error"
      WRITE(gridman_unit,*) 
     w     " Intersection does not cover the whole edge." 
      WRITE(gridman_unit,*) 
     w     " Unfortunately, current version can not ",
     w     " treat this situation properly."       
      WRITE(gridman_unit,*) " DLP, DLQ, DLT, TOL ",dlp,dlq,dlt,tol
      RETURN  
 31   CONTINUE  
       
C     ADD NON-BOUNDARY EDGES OF GRID1
      lmap0=0
      nedges0=nedges
      DO iedge1=1,grid1%NEDGES      
       IF(medge1(iedge1)) THEN    
        nedges0=nedges0+1
C       EDGE INDEX OF THE COMBINED GRID AS FUNCTION OF THE EDGE INDEX OF GRID1
C       THIS TABLE IS REQUIRED TO RE-CALCULATE EDGE INDICES FOR TELEPORTATION
        lmap0(iedge1)=nedges0
       END IF 
      END DO   

      DO iedge1=1,grid1%NEDGES      
       IF(medge1(iedge1)) THEN  
        nedges=nedges+1
        IF(grid1%CELLS(1,iedge1).GE.0) THEN
         cells(1,nedges)=grid1%CELLS(1,iedge1)
        ELSE IF(grid1%CELLS(1,iedge1).LT.0) THEN
C        PRESERVING INDEX FOR TELEPORTATION
         iedge=-grid1%CELLS(1,iedge1)
         IF(iedge.GT.grid1%NEDGES) GOTO 200
         IF(iedge.GT.0) cells(1,nedges)=-lmap0(iedge)
        END IF
        IF(grid1%CELLS(2,iedge1).GE.0) THEN        
         cells(2,nedges)=grid1%CELLS(2,iedge1)
        ELSE IF(grid1%CELLS(2,iedge1).LT.0) THEN     
C        PRESERVING INDEX FOR TELEPORTATION
         iedge=-grid1%CELLS(2,iedge1)
         IF(iedge.GT.grid1%NEDGES) GOTO 200
         cells(2,nedges)=-lmap0(iedge)
        END IF
        points(1,nedges)=grid1%POINTS(1,iedge1)
        points(2,nedges)=grid1%POINTS(2,iedge1)
        imape1(1,nedges)=iedge1
        imape1(2,nedges)=nedges
       END IF
      END DO

C     ADD NON-BOUNDARY EDGES OF GRID2, CONVERT INDICES OF POINTS
      lmap0=0
      nedges0=nedges
      DO iedge2=1,grid2%NEDGES      
       IF(medge2(iedge2)) THEN    
        nedges0=nedges0+1
C       EDGE INDEX OF THE COMBINED GRID AS FUNCTION OF THE  EDGE INDEX OF GRID1
C       THIS TABLE IS REQUIRED TO RE-CALCULATE EDGE INDICES FOR TELEPORTATION
        lmap0(iedge2)=nedges0
       END IF 
      END DO   

      npoints=grid1%NPOINTS
      nedges0=nedges
      ncells=grid1%NCELLS
      DO iedge2=1,grid2%NEDGES      
       IF(medge2(iedge2)) THEN  
        nedges=nedges+1
        IF(grid2%CELLS(1,iedge2).GT.0) THEN
         cells(1,nedges)=grid2%CELLS(1,iedge2)+ncells
        ELSE IF(grid2%CELLS(1,iedge2).LT.0) THEN
C        PRESERVING INDEX FOR TELEPORTATION
         iedge=-grid2%CELLS(1,iedge2)
         IF(iedge.GT.grid2%NEDGES) GOTO 300
         cells(1,nedges)=-lmap0(iedge)
        ELSE 
         cells(1,nedges)=0
        END IF
        IF(grid2%CELLS(2,iedge2).GT.0) THEN
         cells(2,nedges)=grid2%CELLS(2,iedge2)+ncells
        ELSE IF(grid2%CELLS(2,iedge2).LT.0) THEN
C        INDEX FOR TELEPORTATION
         iedge=-grid2%CELLS(2,iedge2)
         IF(iedge.GT.grid2%NEDGES) GOTO 300
         cells(2,nedges)=-lmap0(iedge)
        ELSE
         cells(2,nedges)=0
        END IF
        ipoint=grid2%POINTS(1,iedge2)
        IF(rpoints(ipoint).GT.0) THEN
C        THE POINT WAS REPLACED BY A POINT FROM GRID1
         points(1,nedges)=rpoints(ipoint) 
        ELSE         
C        NO REPLACEMENT, JUST SHIFT
         points(1,nedges)=ipoint+npoints
        END IF
        ipoint=grid2%POINTS(2,iedge2)
        IF(rpoints(ipoint).GT.0) THEN
         points(2,nedges)=rpoints(ipoint)
        ELSE         
         points(2,nedges)=ipoint+npoints
        END IF
        imape2(1,nedges)=iedge2
        imape2(2,nedges)=nedges
       END IF 
      END DO

C     SHIFT POINT INDICES (TO SKIP REPLACED POINTS)
      npoints=grid1%NPOINTS
      mapp=0
      DO ipoint=1,grid2%NPOINTS
       IF(rpoints(ipoint).GT.0) cycle
C      ADD ONLY POINTS WHICH WERE NOT REPLACED
       npoints=npoints+1
       mapp(ipoint)=npoints 
      END DO
        
C     CORRECT POINT INDICES
      npoints1=grid1%NPOINTS     
      DO iedge=1,nedges
       ipoint=points(1,iedge)
       IF(ipoint.GT.npoints1) THEN 
        IF(mapp(ipoint-npoints1).GT.0) THEN 
         points(1,iedge)=mapp(ipoint-npoints1)
        ELSE
         GOTO 50
        END IF      
       END IF
       ipoint=points(2,iedge)
       IF(ipoint.GT.npoints1) THEN
        IF(mapp(ipoint-npoints1).GT.0) THEN 
         points(2,iedge)=mapp(ipoint-npoints1)
        ELSE
         GOTO 50
        END IF   
       END IF
      END DO 
      GOTO 51
 50   ierr=400
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "internal error"
      WRITE(gridman_unit,*) "no shift information for a point, ",
     w                      " IEDGE, IPOINT ",iedge, ipoint
      RETURN
 51   CONTINUE

C     ALLOCATE NEW GRID    
      ncells=grid1%NCELLS+grid2%NCELLS    
      nedgeindex=grid1%NEDGEINDEX+grid2%NEDGEINDEX
      ncellindex=grid1%NCELLINDEX+grid2%NCELLINDEX
      CALL gridman_grid_allocate(grid,2,nedges,npoints,ncells,ierr,
     c                                  nedgeindex,ncellindex)
      IF(ierr.GT.0) GOTO 400 

C     COPY EDGES AND POINTS
      grid%CELLS(1,1:nedges)=cells(1,1:nedges)
      grid%CELLS(2,1:nedges)=cells(2,1:nedges)
      grid%POINTS(1,1:nedges)=points(1,1:nedges)
      grid%POINTS(2,1:nedges)=points(2,1:nedges)    
      DO ipoint=1,grid1%NPOINTS   
       grid%X(1,ipoint)=grid1%X(1,ipoint)
       grid%X(2,ipoint)=grid1%X(2,ipoint)
      END DO
      DO ipoint=1,grid2%NPOINTS
       IF(mapp(ipoint).GT.0) THEN       
        grid%X(1,mapp(ipoint))=grid2%X(1,ipoint)*fscl
        grid%X(2,mapp(ipoint))=grid2%X(2,ipoint)*fscl
       END IF
      END DO
      grid%UNIT2SI=grid1%UNIT2SI !TAKE UNITS OF GRID1
      grid%UNITS=grid1%UNITS
      grid%DESCRIPTION=trim(grid1%DESCRIPTION)//' <merged with> '//
     /                 trim(grid2%DESCRIPTION)

      DO ii=1,grid1%NEDGEINDEX
       CALL gridman_index_transform(grid%EDGEINDEX(ii),
     c          grid1%EDGEINDEX(ii),imape1,nedges,ierr)
       IF(ierr.NE.0) GOTO 410  
      END DO
      DO ii2=1,grid2%NEDGEINDEX
       ii=ii2+grid1%NEDGEINDEX
       CALL gridman_index_transform(grid%EDGEINDEX(ii),
     c        grid2%EDGEINDEX(ii2),imape2,nedges,ierr)
       IF(ierr.NE.0) GOTO 420  
      END DO
      DO ii=1,grid1%NCELLINDEX
       CALL gridman_index_copy(grid%CELLINDEX(ii),
     c                         grid1%CELLINDEX(ii),ierr)
       IF(ierr.NE.0) GOTO 430  
      END DO
      DO ii2=1,grid2%NCELLINDEX
       ii=ii2+grid1%NCELLINDEX
       CALL gridman_index_copy(grid%CELLINDEX(ii),
     c                         grid2%CELLINDEX(ii2),ierr)
       grid%CELLINDEX(ii)%INDEXES(0,:)=grid%CELLINDEX(ii)%INDEXES(0,:)+
     +                                 grid1%NCELLS
       IF(ierr.NE.0) GOTO 440  
      END DO

      CALL gridman_grid2d_check(grid,res,ierr) 
      IF(res.NE.0.OR.ierr.GT.0) GOTO 100 
 
      IF(gridman_dbg) 
     w   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_MERGE finished"

      RETURN

 100  ierr=400
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "internal error - wrong resulting GRID"
      WRITE(gridman_unit,*) " RES ",res
      RETURN  
      
 200  ierr=100
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w      "index for teleportation on GRID1 is out of range"
      WRITE(gridman_unit,*) " IEDGE1, IEDGE ",iedge1, iedge
      RETURN

 300  ierr=100
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w      "index for teleportation on GRID2 is out of range"
      WRITE(gridman_unit,*) " IEDGE2, IEDGE ",iedge1, iedge
      RETURN

 400  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_MERGE terminated"
      RETURN

 410  ierr=420
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "cannot create edge index II1 ",ii
      RETURN
 420  ierr=420
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "cannot create edge index II2 ",ii2
      RETURN
 430  ierr=430
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "cannot create cell index II1 ",ii
      RETURN
 440  ierr=440
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_MERGE: ",
     w                      "cannot create cell index II2 ",ii2
      RETURN

      CONTAINS

C
C     FIND INTERSECTION OF TWO EDGES: OVERLAP AND "HANGING" EDGES
C
C     THIS SUBROUTINE LOOKS SO STRANGE BECAUSE ORIGINALY I WANTED TO 
C     WRITE A GENERAL ALGORITHM WHICH CAN ALSO PROCESS "HANGING TAILS"
C     OF THE INTERSECTIONS
C 
      SUBROUTINE find_overlap(P,Q,XP,YP,XQ,YQ,TOL,IND,TPQ)
      IMPLICIT NONE
      INTRINSIC abs,tiny

      INTEGER(GRIDMAN_SP),INTENT(IN) :: P(2),Q(2) !INDICES OF POINTS OF P AND Q EDGES
      REAL(GRIDMAN_DP),INTENT(IN) :: XP(2),YP(2), !COORDINATES OF THE P AND Q EDGES
     r                               xq(2),yq(2)
      REAL(GRIDMAN_DP),INTENT(IN) :: TOL          !ACCURACY OF GEOMETRICAL OPERATIONS
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: IND(6),  !SEQUENCE OF INDICES IN RESULTING INTERSECTION
     i                                   tpq(2,2) !INDEX P WHICH CORRESPONDS TO Q IF REPLACED

      REAL(GRIDMAN_DP) :: T,DEL,ZP1,ZP2,ZQ1,ZQ2,EPS
      INTEGER(GRIDMAN_SP) :: P1,P2,Q1,Q2  
      INTEGER(GRIDMAN_SP),PARAMETER :: NULL=0
      eps=10.*tiny(eps) !MACHINE PRECISION

      tpq=0
C     POINT Q(1)
      del=abs(yp(2)-yp(1))*(abs(xq(1))+abs(xp(1))+eps)+
     +   abs(xq(1)-xp(1))*(abs(yp(2))+abs(yp(1))+eps)+
     +   abs(xp(2)-xp(1))*(abs(yq(1))+abs(yp(1))+eps)+
     +   abs(yq(1)-yp(1))*(abs(xp(2))+abs(xp(1))+eps)        
      del=tol*(del+eps)
      t=(xq(1)-xp(1))*(yp(2)-yp(1))
      t=t-(yq(1)-yp(1))*(xp(2)-xp(1))
C     DEL IS THE ESTIMATE OF THE ERROR OF T CALCULATED
C     ON ASSUMPTION THAT RELATIVE ERRORS OF ALL X, Y
C     ARE EQUAL TO TOL
      IF(abs(t).GT.del) THEN
C       POINT Q(1) DOES NOT LIE ON THE LINE P, RETURN
       ind=0
       RETURN
      END IF     

C     POINT Q(2)
      del=abs(yp(2)-yp(1))*(abs(xq(2))+abs(xp(1)))+
     +   abs(xq(2)-xp(1))*(abs(yp(2))+abs(yp(1)))+
     +   abs(xp(2)-xp(1))*(abs(yq(2))+abs(yp(1)))+
     +   abs(yq(2)-yp(1))*(abs(xp(2))+abs(xp(1)))        
      del=tol*del
      t=(xq(2)-xp(1))*(yp(2)-yp(1))
      t=t-(yq(2)-yp(1))*(xp(2)-xp(1))
      IF(abs(t).GT.del) THEN
C      POINT Q(2) DOES NOT LIE ON THE LINE P, RETURN
       ind=0
       RETURN
      END IF   

      IF(abs(xp(2)-xp(1)).GT.abs(yp(2)-yp(1))) THEN
C      X PROJECTION OF THE INTERVAL P IS LARGER
C      USE X COORDINATE TO ARRANGE POINTS
       IF(xp(2).GT.xp(1)) THEN
C       PRESCRIBED ORDER
        zp1=xp(1)
        zp2=xp(2)
        p1=p(1)
        p2=p(2)
       ELSE
C       INVERSE ORDER
        zp1=xp(2)
        zp2=xp(1)
        p1=p(2)
        p2=p(1)
       END IF
       IF(xq(2).GT.xq(1)) THEN
C       PRESCRIBED ORDER
        zq1=xq(1)
        zq2=xq(2)
        q1=q(1)
        q2=q(2)
       ELSE
C       INVERSE ORDER
        zq1=xq(2)
        zq2=xq(1)
        q1=q(2)
        q2=q(1)
       END IF
      ELSE
C      Y PROJECTION OF THE INTERVAL P IS LARGER
C      USE Y COORDINATE TO ARRANGE POINTS 
       IF(yp(2).GT.yp(1)) THEN
C       PRESCRIBED ORDER
        zp1=yp(1)
        zp2=yp(2)         
        p1=p(1)
        p2=p(2)
       ELSE
C       INVERSE ORDER
        zp1=yp(2)
        zp2=yp(1)
        p1=p(2)
        p2=p(1)
       END IF      
       IF(yq(2).GT.yq(1)) THEN
C       PRESCRIBED ORDER     
        zq1=yq(1)
        zq2=yq(2)
        q1=q(1)
        q2=q(2)
       ELSE
C       INVERSE ORDER
        zq1=yq(2)
        zq2=yq(1)
        q1=q(2)
        q2=q(1)
       END IF    
      END IF    
  
      IF(abs(zq1-zp1).LT.tol*(abs(zq1)+abs(zp1)+eps).AND.
     f   abs(zq2-zp2).LT.tol*(abs(zq2)+abs(zp2)+eps)) THEN
       ind=(/null,null,p1,p2,null,null/)    !1
       tpq(1,:)=(/q1,q2/)
       tpq(2,:)=(/p1,p2/)
      ELSE IF(abs(zq1-zp1).LT.tol*(abs(zq1)+abs(zp1)+eps)) THEN
       IF(zq2.GT.zp2) THEN  
        ind=(/null,null,p1,p2,p2,q2/) !2      
       ELSE
        ind=(/null,null,p1,q2,q2,p1/) !3
       END IF
       tpq(1,1)=q1
       tpq(2,1)=p1  
      ELSE IF(abs(zq2-zp2).LT.tol*(abs(zq2)+abs(zp2)+eps)) THEN
       IF(zq1.LT.zp1) THEN
        ind=(/q1,p1,p1,p2,null,null/) !4
       ELSE
        ind=(/p1,q1,q1,p2,null,null/) !5
       END IF    
       tpq(1,1)=q2
       tpq(2,1)=p2  
      ELSE IF(abs(zq2-zp1).LT.tol*(abs(zq2)+abs(zp1)+eps)) THEN
       ind=(/q1,p1,p1,null,p1,p2/) !6
       tpq(1,1)=q2
       tpq(2,1)=p1
      ELSE IF(abs(zq1-zp2).LT.tol*(abs(zq1)+abs(zp2)+eps)) THEN
       ind=(/p1,p2,null,p2,p2,q2/) !7
       tpq(1,1)=q1
       tpq(2,1)=p2
      ELSE IF(zq1.LT.zp1) THEN
       IF(zq2.LT.zp1) THEN
        ind=(/q1,q2,null,null,p1,p2/) !8
       ELSE IF(zq2.GT.zp2) THEN 
        ind=(/q1,p1,p1,p2,p2,q2/) !9
       ELSE
        ind=(/q1,p1,p1,q2,q2,p1/) !10
       END IF
      ELSE IF(zq1.GT.zp2) THEN
        ind=(/p1,p2,null,null,q1,q2/)   !11
      ELSE ! ZP1<ZQ1<ZP2
       IF(zq2.LT.zp2) THEN
        ind=(/p1,q1,q1,q2,q2,p2/) !12
       ELSE
        ind=(/p1,q1,q1,p2,p2,q2/) !13
       END IF 
      END IF

      END SUBROUTINE find_overlap

C
C     RETURN ONE OF TWO CELLS WHICH IS LARGER THAN 0
C     IF BOTH LE 0 THEN RETURNS 0
C
      FUNCTION nonzer_cell(CELLS)
       INTEGER(GRIDMAN_SP),INTENT(IN) :: CELLS(2)
       INTEGER(GRIDMAN_SP) :: NONZER_CELL
       IF(cells(1).GT.0) THEN
        nonzer_cell=cells(1)
       ELSE IF(cells(2).GT.0) THEN
        nonzer_cell=cells(2)
       ELSE
        nonzer_cell=0         
       END IF
      END FUNCTION nonzer_cell      

      END SUBROUTINE gridman_grid2d_merge