cut.f

Go to the documentation of this file.

C> @file grid2D/cut.f
C> Select part of 2D grid cut by polygon  
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**********************************************************************************************
C> Select part of a 2D grid cut by polygon  
C>
C> Take cells which lie inside or outside of a closed polygon. 
C> 
C> Current implementation cannot handle the case where polygon edges 
C> intersect the grid nodes, or grid edges intersect the polygon vertices
C>
C> Normally those conditions, as well as the codition that the polygon 
C> has no self intersections, are not checked because it may require 
C> too long time. The checks are enforced when GRIDMAN_CHECK=.TRUE. 
C>
C> If polygon segments are added to the grid, then an edge index with
C> indexes of the segments is added
      SUBROUTINE gridman_grid2d_cut(CUTGRID,GRID,XP,YP,NP,LEX,IERR)
      USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
     u                 gridman_dbg,gridman_check,gridman_indlist,
     u                 gridman_tol
      USE gridman_lib,ONLY:gridman_grid2d_check
      IMPLICIT NONE
      INTRINSIC tiny,count,abs,maxval,minval
C>    Resulting grid with eliminated and cut cells 
      TYPE(gridman_grid) :: CUTGRID
C>    Initial grid
      TYPE(gridman_grid) :: GRID
C>    Number of points in the polygon  
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NP
C>    X-coordinates of the polygon vertices, in METER
C>
C>    First and last points must be same - closed polygon 
      REAL(GRIDMAN_DP),INTENT(IN) :: XP(np)
C>    Y-coordinates of the polygon vertices, in METER
C>
C>    First and last points must be same - closed polygon 
      REAL(GRIDMAN_DP),INTENT(IN) :: YP(np)
C>    Defines if points inside the polygon are included or excluded 
C>
C>    If .FALSE. then points are included
C>    If .TRUE. then points are excluded 
      LOGICAL,INTENT(IN) :: LEX
C>    Error code
      INTEGER,INTENT(OUT) :: IERR   

      INTERFACE
       SUBROUTINE gridman_grid2d_polygon_trajectory(NPOL,NPEDGE,
     s            icell_pol,iedge_pol,x_pol,y_pol,indvert,indseg,
     s            xp,yp,np,grid,ierr)
       USE gridman,ONLY:gridman_grid,gridman_sp,gridman_dp
       TYPE(gridman_grid) :: GRID        
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: NPOL,NPEDGE(grid%nedges) 
       INTEGER(GRIDMAN_SP),ALLOCATABLE :: ICELL_POL(:),
     i                                    iedge_pol(:),indseg(:)                     
       REAL(GRIDMAN_DP),ALLOCATABLE :: X_POL(:),Y_POL(:)
       INTEGER(GRIDMAN_SP),INTENT(IN) :: NP 
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: INDVERT(np)
       REAL(GRIDMAN_DP) :: XP(np),YP(np)
       INTEGER,INTENT(OUT) :: IERR       
       END SUBROUTINE gridman_grid2d_polygon_trajectory
      END INTERFACE

C     LOCAL VARIABLES
      LOGICAL :: LFIRST_POINT_INSIDE_POLYGON
      LOGICAL :: LPOINT(grid%npoints),LCELL(grid%ncells) 
      INTEGER(GRIDMAN_SP) :: NPOL,NPMAX,NEMAX,NWORK,IEP,
     i                       npoints,nedges,ncells,ncells0,nseg,
     i                       ip,icell,icell1,iedge,iedge1,iedge0,
     i                       ipoint,ipoint1,ipoint2,
     i                       ie,ic,ne,n,ncount,st
      INTEGER :: RES,IERR0,IWARN
      INTEGER (GRIDMAN_SP) :: NPEDGE(grid%nedges),       
     i                        npcell(grid%NCELLS),
     i                        inewcell(grid%NCELLS),  !INDMAP OF CELLS, OLD->NEW
     i                        inewedge(grid%NEDGES),  !INDMAP OF EDGES, OLD->NEW
     i                        indpoint(grid%NPOINTS), !INDMAP OF POINTS, OLD->NEW
     i                        indvert(np)               !INDMAP OF POLYGON VERTICES, OLD->NEW
!    INDVERT is at the moment not used, is calculated in GRIDMAN_GRID2D_POLYGON_TRAJECTORY
!    for tje case if indices of points will be added in the future
      REAL(GRIDMAN_DP) :: X0,Y0
C     ALLOCATABLE ARRAYS
      REAL(GRIDMAN_DP),ALLOCATABLE :: X1(:),X2(:) 
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: IOLDEDGE(:), !INDMAP OF EDGES, NEW->OLD
     i                                   ioldcell(:), !INDMAP OF CELLS, NEW->OLD 
     i  icell_pol(:),iedge_pol(:),indseg_pol(:),  !see  GRIDMAN_GRID2D_POLYGON_TRAJECTORY for description
     i  cells(:,:),points(:,:) !temporary arrays for cells and points of the new grid
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: INDSEG(:,:) !INDSEG(1,:) - index of edge of the new grid, 
                                                     !INDSEG(2,:) - index of the polygon segment
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: IPWORK(:),TPWORK(:),ICOUNT(:)  
      REAL(GRIDMAN_DP),ALLOCATABLE :: X_POL(:),Y_POL(:),DWORK(:)
      LOGICAL,ALLOCATABLE :: LWORK(:)

      TYPE(gridman_indlist) :: CELLEDGES

C     CHECK INPUT GRID
      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "Starting GRIDMAN_GRID2D_CUT"

      ierr=0

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

      IF(gridman_check) THEN
C      CHECK INPUT GRID
       CALL gridman_grid2d_check(grid,res,ierr)
       IF(res.NE.0.OR.ierr.GT.0) THEN
        ierr=100
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "incorrect input grid"
        RETURN
       END IF
      END IF   

C     CHECK POLYGON
      CALL check_polygon(xp,yp,np,ierr)
      IF(ierr.GT.0) GOTO 100
      iwarn=ierr

C     1. CALCULATE TRAJECTORY OF THE POLYGON
      CALL gridman_grid2d_polygon_trajectory(npol,npedge,
     s                                       icell_pol,iedge_pol,
     s                                       x_pol,y_pol, 
     s                                       indvert,indseg_pol,
     s                                       xp,yp,np,grid,ierr)
      IF(ierr.NE.0) GOTO 100

C     2. FIND WHICH GRID POINTS (NODES) ARE INSIDE/OUTSIDE OF THE POLYGON 
      lfirst_point_inside_polygon=.true.
      DO ipoint=1,grid%NPOINTS
       x0=grid%X(1,ipoint)*grid%UNIT2SI
       y0=grid%X(2,ipoint)*grid%UNIT2SI
       lpoint(ipoint)=point_inside_polygon(xp,yp,np,x0,y0,lex)
      END DO   
C     NUMBER OF POINTS IN NEW GRID
      npoints=count(lpoint(1:grid%NPOINTS)) 
      npoints=npoints+npol       

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT: NPOINTS ",npoints

C     ALLOCATE TEMPORARY ARRAYS 
      npmax=maxval(npedge)
      CALL allocate_temporary(npoints,npol,npmax,ierr) 
      IF(ierr.GT.0) GOTO 100

C     3. TAKE COORDINATES OF THE POINTS TO THE NEW GRID      
      npoints=npol 
      IF(npol.GT.0) THEN
       x1(1:npol)=x_pol
       x2(1:npol)=y_pol
      END IF
      indpoint=0
      DO ipoint=1,grid%NPOINTS
       IF(lpoint(ipoint)) THEN
        npoints=npoints+1
        indpoint(ipoint)=npoints 
        x1(npoints)=grid%X(1,ipoint)
        x2(npoints)=grid%X(2,ipoint)
       END IF
      END DO

C     4. BUILD ARRAY OF EDGES
      nedges=0
      cells=0
      points=0
      ioldedge=0
      inewedge=0
      DO iedge=1,grid%NEDGES
C      BUILD TABLE OF POINTS FOR THE EDGE
C       IPWORK - index of point in the array of points
C       TPWORK - type of the point: in, out or boundary
       CALL edge_points(ipwork,tpwork,nwork,iedge,ierr) 
       IF(ierr.GT.0) GOTO 100
       DO iep=1,nwork-1 
        IF( (tpwork(iep).GT.0.AND.tpwork(iep+1).GE.0).OR.
     f      (tpwork(iep+1).GT.0.AND.tpwork(iep).GE.0)) THEN
C        IF BOTH IEP AND IEP+1 ARE IN OR 
C        IF ONE POINT IS IN AND ONE ON THE BOUNDARY THEN ADD EDGE
         nedges=nedges+1
         cells(:,nedges)=grid%CELLS(:,iedge)
         points(1,nedges)=ipwork(iep)
         points(2,nedges)=ipwork(iep+1)
         ioldedge(nedges)=iedge
         inewedge(iedge)=nedges
        ELSE IF( (tpwork(iep).EQ.0.AND.tpwork(iep+1).EQ.0) ) THEN 
C        IF BOTH ON THE BOUNDARY THEN TAKE MIDDLE POINT AND
C        CHECK IF THIS POINT INSIDE THE POLYGON
         ipoint1=ipwork(iep)
         ipoint2=ipwork(iep+1)
         x0=0.5*(x1(ipoint1)+x1(ipoint2))*grid%UNIT2SI
         y0=0.5*(x2(ipoint1)+x2(ipoint2))*grid%UNIT2SI
         IF(point_inside_polygon(xp,yp,np,x0,y0,lex)) THEN
          nedges=nedges+1
          cells(:,nedges)=grid%CELLS(:,iedge)
          points(1,nedges)=ipoint1
          points(2,nedges)=ipoint2
          ioldedge(nedges)=iedge
          inewedge(iedge)=nedges
         END IF 
        END IF    
       END DO 
      END DO !DO IEDGE=1,GRID%NEDGES

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT: NEDGES (*) ",nedges

C     5. ADD EDGES FORMED BY THE POLYGON SEGMENTS
      nseg=0 !number of edges which point to the polygon segments
      IF(npol.GT.0) THEN
       ALLOCATE(indseg(2,npol),stat=st)
       IF(st.NE.0) GOTO 400  
       DO ip=1,npol
        icell=icell_pol(ip)
        IF(icell.GT.0) THEN
C        ADD ONLY IF THE INTERVAL OF THE TRAJECTORY
C        INTERSECTS A CELL
         nedges=nedges+1
         cells(1,nedges)=icell
         cells(2,nedges)=0
         points(1,nedges)=ip
         IF(ip.LT.npol) THEN
          points(2,nedges)=ip+1
         ELSE
          points(2,nedges)=1
         END IF
         nseg=nseg+1
         indseg(1,nseg)=nedges
         indseg(2,nseg)=indseg_pol(ip)
        END IF
       END DO
      END IF

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT: NEDGES ",nedges

      IF(nedges.EQ.0) GOTO 200

C     6. INDEX MAPPING OF CELLS
      lcell=.false. 
      DO iedge=1,nedges
       icell=cells(1,iedge)
       IF(icell.GT.0) lcell(icell)=.true.
       icell=cells(2,iedge)
       IF(icell.GT.0) lcell(icell)=.true.
      END DO
      ncells=0 
      inewcell=0
      ioldcell=0
      DO icell=1,grid%NCELLS
       IF(lcell(icell)) THEN
        ncells=ncells+1
        inewcell(icell)=ncells
        ioldcell(ncells)=icell
       END IF
      END DO       

C     7. TRANSLATE CELL INDICES INTO NEW INDEXING    
      npcell(1:ncells)=0 !NUMBER OF INTERSECTIONS  WITH THE CELL
      DO iedge=1,nedges 
       DO ic=1,2
        icell=cells(ic,iedge)            
        IF(icell.GT.0) THEN
         icell1=inewcell(icell)
         IF(icell1.LT.1) GOTO 300
         cells(ic,iedge)=icell1
         iedge0=ioldedge(iedge) 
         IF(iedge0.GT.0)          
     f      npcell(icell1)=npcell(icell1)+npedge(iedge0)
        ELSEIF(icell.LT.0) THEN
C        EDGE INDEX FOR TELEPORTATION
         iedge1=inewedge(-icell)
         cells(ic,iedge)=-iedge1
        END IF
       END DO
      END DO

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT: NCELLS (*) ",ncells

C     8. CORRECT CELL INDICES 
C     TABLE OF EDGES WHICH BELONG TO EACH CELL    
      IF(maxval(npcell).GT.2) THEN 
       CALL list_of_edges(celledges,ierr)     
       IF(ierr.GT.0) GOTO 100 
      END IF
      ncells0=ncells
      DO icell=1,ncells0
C      SKIP CELLS WITH LESS THAN 2 INTERSECTIONS WHICH DO NOT REQUIRE EXTRA CORRECTION 
       IF(npcell(icell).LE.2) cycle 
C      TOTAL NUMBER OF ELEMENTS
       ne=celledges%ILAST(icell)-celledges%IFIRST(icell)+1
       n=0
C      UNMARK ALL ELEMENTS
       lwork(1:ne)=.true. 
       ic=1
 10    CONTINUE
C      FIND CLOUSED COUNTOUR CONSISTING OF UNMARKED ELEMENTS
       CALL find_contour(icount,ncount,lwork,icell,ierr)        
       IF(ierr.NE.0) RETURN
       IF(ic.GT.1) THEN
C       IF MORE THAN ONE CONTOUR, THEN ADD NEW CELL
        ncells=ncells+1
        ioldcell(ncells)=ioldcell(icell)
C       CORRECT CELL INDEX FOR ALL EDGES BELONGING TO THE CONTOUR
        DO ie=1,ncount
         iedge=icount(ie)
         IF(cells(1,iedge).EQ.icell) THEN
          cells(1,iedge)=ncells
         ELSE
          cells(2,iedge)=ncells
         END IF          
        END DO
       END IF
       n=n+ncount
       IF(n.LT.ne) THEN   
C       IF NOT ALL ELEMENTS OF THE CELL ARE IN THE CONTOUR,
C       THEN REPEAT ITERATION
        ic=ic+1
        GOTO 10
       END IF
      END DO !DO ICELL=1,GRID%NCELLS

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT: NCELLS ",ncells

      CALL create_cutgrid(ierr)
      IF(ierr.GT.0) GOTO 100

      CALL deallocate_temporary(ierr)

C     TO PRESERVE WARNING INDEX FROM CHECK_POLYGON
      IF(ierr.EQ.0) ierr=iwarn 

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT finished"

      RETURN

 100  WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT terminated"
      CALL deallocate_all(ierr0) 
      RETURN
 200  ierr=100 
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                      "resulting grid is empty"
      CALL deallocate_all(ierr0) 
      RETURN
 300  ierr=400 
      WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                      "internal error"
      WRITE(gridman_unit,*) 
     w     "Cell is taken into new grid but the index is <1"
      WRITE(gridman_unit,*) " ICELL, ICELL1 ",icell,icell1
      CALL deallocate_all(ierr0) 
 400  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                      "cannot allocate temporary array INDSEG"    
      WRITE(gridman_unit,*) " NPOL ",npol
      CALL deallocate_all(ierr0) 
      RETURN

      CONTAINS

C --------------------------------------------------------------------
C     FIND CHAIN OF CONNECTED EDGES BELONGING TO ICELL 
C     WHICH FORM CLOSED CONTOUR
C     TAKE ONLY EDGES WHICH WERE NOT MARGED IN LWORK AS EXCLUDED 
C     LWORK IS UPDATED
C --------------------------------------------------------------------
      SUBROUTINE find_contour(ICOUNT,NCOUNT,LWORK,ICELL,IERR)
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: ICOUNT(nemax), !LIST OF EDGES IN THE CONTOUR
     i                                    ncount         !NUMBER OF EDGES IN THE CONTOUR
       INTEGER(GRIDMAN_SP),INTENT(IN) :: ICELL
       LOGICAL,INTENT(INOUT) :: LWORK(nemax) 
       INTEGER,INTENT(OUT) :: IERR      

       INTEGER(GRIDMAN_SP) :: IE,I,IEDGE,IPOINT,IPOINT0,IPOINT00,NMAX

       ierr=0

C      FIND FIRST EDGE WHICH WAS NOT YET EXCLUDED     
       DO ie=celledges%IFIRST(icell),celledges%ILAST(icell)
        iedge=celledges%IND(ie)
        i=ie-celledges%IFIRST(icell)+1
        IF(lwork(i)) EXIT
       END DO 
C      TAKE ONE POINT OF THIS EDGE  
       ipoint0=points(1,iedge)
       ipoint00=ipoint0

       nmax=celledges%ILAST(icell)-celledges%IFIRST(icell)+1

       ncount=0       
       DO 
C       SCAN OVER ALL EDGES WHICH WERE NOT YET EXLCUDED
C       LOOK FOR THE POINT WHICH IS SAME AS IPOINT0 
        DO ie=celledges%IFIRST(icell),celledges%ILAST(icell)
         i=ie-celledges%IFIRST(icell)+1
         IF(lwork(i)) THEN
          iedge=celledges%IND(ie)
          IF(iedge.LT.1) THEN
           ierr=400
           WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "intenal error" 
           WRITE(gridman_unit,*) "No IEDGE in CELLEDGES"
           WRITE(gridman_unit,*) " ICELL, IE, IEDGE ",icell,ie,iedge
           RETURN 
          END IF
          ipoint=points(1,iedge)
          IF(ipoint.EQ.ipoint0) THEN
C          1st POINT = IPOINT0, THEN 2nd POINT IS THE NEXT IPOINT0
           ncount=ncount+1
           IF(ncount.GT.nmax) GOTO 100
           icount(ncount)=iedge
           ipoint0=points(2,iedge)
           lwork(i)=.false. 
           GOTO 20
          END IF 
          ipoint=points(2,iedge)
          IF(ipoint.EQ.ipoint0) THEN
C          2nd POINT = IPOINT0, THEN 1st POINT IS THE NEXT IPOINT0
           ncount=ncount+1
           IF(ncount.GT.nmax) GOTO 100
           icount(ncount)=iedge
           ipoint0=points(1,iedge)
           lwork(i)=.false. 
           GOTO 20
          END IF 
         END IF !IF(LWORK(I)) THEN
        END DO   
C       NO NEXT POINT IS FOUND - AN ERROR
        ierr=400
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "intenal error" 
        WRITE(gridman_unit,*) "Cannot find closed contour"
        WRITE(gridman_unit,*) " ICELL ",icell
        RETURN       
C       IF CURRENT POINT POINT0 == VERY FIRST POINT POINT00,
C       THEN THE CONTOUR IS CLOSED  
 20     IF(ipoint0.EQ.ipoint00) EXIT 
       END DO !DO
    
       RETURN

 100   ierr=400
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                       "intenal error" 
       WRITE(gridman_unit,*) "Index out of range in FIND_CONTOUR"
       WRITE(gridman_unit,*) " ICELL ",icell
       RETURN    

      END SUBROUTINE find_contour

C --------------------------------------------------------------------
C     CREATE TABLE (INDLIST) OF EDGES WHICH BELONG TO EACH CELL
C     ACCORDING TO ARRAY CELLS  
C --------------------------------------------------------------------
      SUBROUTINE list_of_edges(CELLEDGES,IERR)
       TYPE(gridman_indlist) :: CELLEDGES          
       INTEGER,INTENT(OUT) :: IERR
       INTEGER(GRIDMAN_SP) :: IEDGE,IC,ICELL,L,ST,IE,NLAST

       ierr=0
 
C      Here memory usage can be further optimized  because the list of 
C      edges is required not for all cells but only for cells
C      with >2 intersections
       ALLOCATE(celledges%IFIRST(ncells),
     a          celledges%ILAST(ncells),stat=st)      
       IF(st.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w      "cannot allocate temporary data structure CELLEDGES" 
        WRITE(gridman_unit,*) " NCELLS ",grid%NCELLS
        ierr=200
        RETURN
       END IF
       celledges%ILAST=0
       l=0
       DO iedge=1,nedges
        DO ic=1,2
         icell=cells(ic,iedge)
         IF(icell.GT.0) THEN
          celledges%ILAST(icell)=celledges%ILAST(icell)+1
          l=l+1
         END IF
        END DO
       END DO
C      CELLEDGES%ILAST NOW IS THE NUMBER OF EDGES WHICH BELONG TO ICELL
       nemax=maxval(celledges%ILAST)
       celledges%N=ncells
       celledges%L=l
       ALLOCATE(celledges%IND(l),stat=st)
       IF(st.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w      "cannot allocate temporary data structure CELLEDGES" 
        WRITE(gridman_unit,*) " L ",l
        ierr=200
        RETURN
       END IF
       celledges%IFIRST(1)=1
       nlast=celledges%ILAST(1)       
       celledges%ILAST(1)=0
       DO icell=2,celledges%N
        celledges%IFIRST(icell)=nlast+1
        nlast=nlast+celledges%ILAST(icell)
        celledges%ILAST(icell)=celledges%IFIRST(icell)-1
       END DO
       DO iedge=1,nedges
        DO ic=1,2
         icell=cells(ic,iedge)
         IF(icell.GT.0) THEN
          celledges%ILAST(icell)=celledges%ILAST(icell)+1         
          ie=celledges%ILAST(icell)
          celledges%IND(ie)=iedge
         END IF
        END DO
       END DO
       ALLOCATE(lwork(nemax),icount(nemax),stat=st)
       IF(st.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "cannot allocate temporary arrays" 
        WRITE(gridman_unit,*) " NEMAX ",nemax
        ierr=200
        RETURN
       END IF

      END SUBROUTINE list_of_edges
      
C --------------------------------------------------------------------
C     BUILD CHAIN OF POINTS WHICH REPRESENTS IEDGE
C --------------------------------------------------------------------
      SUBROUTINE edge_points(IPEDGE,TPEDGE,NP,IEDGE,IERR)
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: IPEDGE(npmax+2), !INDEX OF POINT IN X1,X2
     i                                                    !FIRST AND LAST POINTS
     i                                                    !ARE ALWAYS EDGE ENDS  
     i                                   tpedge(npmax+2), !TYPE OF POINT:
     i                                                    ! 1: INSIDE POLYGON
     i                                                    !-1: OUTSIDE POLYGON
     i                                                    ! 0: ON THE BOUNDARY 
     i                                    np              !NUMBER OF POINTS, >=2
      INTEGER(GRIDMAN_SP),INTENT(IN) :: IEDGE
      INTEGER,INTENT(OUT) :: IERR
      
      INTEGER(GRIDMAN_SP) :: IPOINT,IP,IPOINT1,ITMP
      REAL(GRIDMAN_DP) :: Xa,Ya,Xs,Ys,DTMP
      LOGICAL :: LSWAP

      IF(npedge(iedge).LT.1) THEN       
C      NO INTERSECTION WITH POLYGON
C      WHOLE EDGE IS IN OR OUT
       DO ip=1,2 
        ipoint=grid%POINTS(ip,iedge)
        IF(lpoint(ipoint)) THEN
         tpedge(ip)=1
        ELSE
         tpedge(ip)=-1  
        END IF
        ipedge(ip)=indpoint(ipoint)
       END DO
       IF(tpedge(1).NE.tpedge(2)) THEN
        ierr=400
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "internal error"
        WRITE(gridman_unit,*) 
     w      "If edge is not intersected by the polygon"
        WRITE(gridman_unit,*) 
     w      "then both its ends must be either inside ",
     w      "or outside the grid"
        WRITE(gridman_unit,*) " IEDGE, TPEDGE ",iedge,tpedge(1:2)
       END IF
       np=2
      ELSE !IF(NPEDGE(IEDGE).EQ.0) THEN   
C      EDGE INTERSECTS WITH POLYGON
C      BUILD CHAINS OF INTERSECTION POINTS
       np=1
C      TAKE ONE END OF THE EDGE
       ipoint=grid%POINTS(1,iedge)
       ipoint1=indpoint(ipoint)
       ipedge(np)=ipoint1
       xa=grid%X(1,ipoint)
       ya=grid%X(2,ipoint)
       IF(lpoint(ipoint)) THEN
        tpedge(np)=1
       ELSE
        tpedge(np)=-1  
       END IF
       DO ip=1,npol !CYCLE OVER WHOLE TRAJECTORY
        IF(iedge_pol(ip).EQ.iedge) THEN
C        TRAJECTORY INTERSECTS IEDGE
C        ADD INTERSECTION POIN TO IPEDGE, TPEDGE
         np=np+1
         ipoint1=ip
         ipedge(np)=ipoint1
         tpedge(np)=0
         xs=x1(ipoint1)
         ys=x2(ipoint1)
C        DISTANCE BETWEEN (Xs,Ys) AND (Xa,Ya)
         dwork(np)=(xa-xs)**2+(ya-ys)**2 
        END IF
       END DO
C      ADD ANOVER END OF THE EDGE
       np=np+1
       ipoint=grid%POINTS(2,iedge)
       ipoint1=indpoint(ipoint)
       ipedge(np)=ipoint1
       IF(lpoint(ipoint)) THEN
        tpedge(np)=1
       ELSE
        tpedge(np)=-1  
       END IF
       IF(np.GT.3) THEN
C       SORT IN INCREASING DISTANCE FROM (Xa,Ya) ORDER
C       SIMPLE BUBBLE SORT ALGORITHM
 10     lswap=.false.
        DO ip=2,np-2
         IF(dwork(ip).GT.dwork(ip+1)) THEN
          dtmp=dwork(ip)
          dwork(ip)=dwork(ip+1)
          dwork(ip+1)=dtmp  
          itmp=ipedge(ip)
          ipedge(ip)=ipedge(ip+1)
          ipedge(ip+1)=itmp
          itmp=tpedge(ip)
          tpedge(ip)=tpedge(ip+1)
          tpedge(ip+1)=itmp
          lswap=.true.
         END IF           
        END DO
        IF(lswap) GOTO 10
       END IF
      END IF !IF(NPEDGE(IEDGE).EQ.0) THEN   

      END SUBROUTINE edge_points

C --------------------------------------------------------------------
C     FIND IF THE POINT (X0,Y0) LIES
C     INSIDE THE POLYGON DEFINED BY XP, YP
C     THE RESULT IS INVERTED WHEN LEX=.TRUE,
C --------------------------------------------------------------------
      FUNCTION point_inside_polygon(XP,YP,NP,X0,Y0,LEX)
       USE gridman_lib,ONLY:gridman_point_in_polygon
       INTEGER(GRIDMAN_SP),INTENT(IN) :: NP
       REAL(GRIDMAN_DP),INTENT(IN) :: XP(np),YP(np),X0,Y0
       LOGICAL,INTENT(IN) :: LEX
       LOGICAL :: POINT_INSIDE_POLYGON

       REAL(GRIDMAN_DP),SAVE :: XMIN,XMAX,YMIN,YMAX
    
       IF(lfirst_point_inside_polygon) THEN
C       FOR OPTIMIZATION
        xmin=minval(xp)
        xmax=maxval(xp)
        ymin=minval(yp)
        ymax=maxval(yp)
        lfirst_point_inside_polygon=.false.
       END IF

       IF(x0.GT.xmax.OR.x0.LT.xmin.OR.
     f    y0.GT.ymax.OR.y0.LT.ymin) THEN
        IF(lex) THEN
         point_inside_polygon=.true.
        ELSE
         point_inside_polygon=.false.
        END IF
        RETURN
       END IF 
       
       point_inside_polygon=gridman_point_in_polygon(xp,yp,np,x0,y0) 
       IF(lex) point_inside_polygon=.NOT.point_inside_polygon

      END FUNCTION point_inside_polygon

C --------------------------------------------------------------------
C     CREATE NEW GRID OBJECT OUT OF TEMPORARY ARRAYS     
C --------------------------------------------------------------------
      SUBROUTINE create_cutgrid(IERR)
      USE gridman_lib,ONLY:gridman_grid_allocate
      INTRINSIC trim
      INTEGER,INTENT(OUT) :: IERR
      INTEGER(GRIDMAN_SP) :: ICELL,ICELL0,IEDGE,IEDGE0
      INTEGER :: II,NEI
      INTEGER (GRIDMAN_SP) :: M,IEMAP(2,nedges),ICMAP(2,ncells)

      IF(nseg.GT.0) THEN
       nei=grid%NEDGEINDEX+1
      ELSE
       nei=grid%NEDGEINDEX
      END IF
      CALL gridman_grid_allocate(cutgrid,grid%TYPE,nedges,npoints, 
     s                           ncells,ierr,nei,grid%NCELLINDEX)
      IF(ierr.NE.0) THEN
       WRITE(gridman_unit,*) "GRIDMAN_GRID2D_CUT terminated"
       RETURN
      END IF  
      cutgrid%CELLS=cells(1:2,1:nedges)
      cutgrid%POINTS(1:2,:)=points(1:2,1:nedges) 
      cutgrid%X(1,:)=x1(1:npoints)
      cutgrid%X(2,:)=x2(1:npoints)
      cutgrid%UNIT2SI=grid%UNIT2SI     
      cutgrid%UNITS=grid%UNITS
      cutgrid%DESCRIPTION=trim(grid%DESCRIPTION)//
     /                   " << processed by GRIDMAN_GRID2D_CUT"

C     TRANSFORM EDGE INDICES   
      DO ii=1,grid%NEDGEINDEX
       m=0
       DO iedge=1,nedges
        iedge0=ioldedge(iedge)
        IF(iedge0.GT.0) THEN
         m=m+1
         iemap(1,m)=iedge0 
         iemap(2,m)=iedge
        END IF
       END DO 
       CALL gridman_index_transform(cutgrid%EDGEINDEX(ii),
     c                              grid%EDGEINDEX(ii),iemap,m,ierr)
       IF(ierr.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "cannot create edge index II  ",ii
        RETURN
       END IF
      END DO
C     ADD INDEX MAPPING TO THE POLYGON SEGMENTS
      IF(nseg.GT.0) THEN
       CALL gridman_index_allocate(cutgrid%EDGEINDEX(nei),1,nseg,ierr)
       IF(ierr.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w       "cannot create the polygon segment index, NSEG ",nseg
        RETURN
       END IF
       cutgrid%EDGEINDEX(nei)%INDEXES(0:1,1:nseg)=indseg(1:2,1:nseg)
       cutgrid%EDGEINDEX(nei)%DESCRIPTION=
     w                       'Indices of the polygon segments'
       cutgrid%EDGEINDEX(nei)%COLUMNS(1)='CUT_ISEG'
      END IF

C     TRANSFORM CELL INDICES   
      DO ii=1,grid%NCELLINDEX
       m=0
       DO icell=1,ncells
        icell0=ioldcell(icell)
        IF(icell0.GT.0) THEN
         m=m+1
         icmap(1,m)=icell0 
         icmap(2,m)=icell
        END IF
       END DO 
       CALL gridman_index_transform(cutgrid%CELLINDEX(ii),
     c                              grid%CELLINDEX(ii),icmap,m,ierr)
       IF(ierr.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                        "cannot create edge index II  ",ii
        RETURN
       END IF
      END DO

      CALL gridman_grid2d_check(cutgrid,res,ierr) 
      IF(res.NE.0.OR.ierr.GT.0) THEN
       ierr=400
       WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                       "wrong resulting GRID object"
       RETURN       
      END IF

      END SUBROUTINE create_cutgrid

C --------------------------------------------------------------------
C     CHECK IF THE POLYGON CAN BE USED
C --------------------------------------------------------------------
      SUBROUTINE check_polygon(XP,YP,NP,IERR)
      USE gridman_lib,ONLY:gridman_intersect2d,gridman_close2interval2d
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NP
      REAL(GRIDMAN_DP),INTENT(IN) :: XP(np),YP(np)
      INTEGER,INTENT(OUT) :: IERR                
      INTEGER(GRIDMAN_SP) :: IP,IP1,NP1,IPOINT
      REAL(GRIDMAN_DP) :: DX,EX,DY,EY,EPS,X1,Y1,X2,Y2,X0,Y0,
     r                    x11,y11,x12,y12,x21,y21,x22,y22,xi,yi
       ierr=0

       IF(np.LT.4) THEN
        ierr=100
        WRITE(gridman_unit,*) "ERROR in  GRIDMAN_GRID2D_CUT: ",
     w                    "polygon must have at least 3 vertices" 
        WRITE(gridman_unit,*) " NP-1 ",np-1
        RETURN
       END IF

       DO ip=1,np-1
        dx=abs(xp(ip+1)-xp(ip))
        dy=abs(yp(ip+1)-yp(ip))
        ex=tiny(dx)*10
        IF(dx.LT.ex.AND.dy.LT.ex) THEN
         ierr=100
         WRITE(gridman_unit,*) "ERROR in  GRIDMAN_GRID2D_CUT: ",
     w                         "duplicate points in the polygon"
         WRITE(gridman_unit,*) " IP1, X1, Y1 ",ip,xp(ip),yp(ip)
         WRITE(gridman_unit,*) " IP2, X2, Y2 ",
     w                           ip+1,xp(ip+1),yp(ip+1)
         RETURN
        END IF
       END DO 

       ex=gridman_tol*(abs(xp(1))+abs(xp(np)))+tiny(dx)*10
       dx=abs(xp(1)-xp(np))
       ey=gridman_tol*(abs(yp(1))+abs(yp(np)))+tiny(dy)*10
       dy=abs(yp(1)-yp(np))
       IF(dx.GT.ex.OR.dy.GT.ey) THEN
        ierr=100
        WRITE(gridman_unit,*) "ERROR in  GRIDMAN_GRID2D_CUT: ",
     w                        "no closed polygon"
       END IF          

       IF(gridman_check) THEN
C       TESTS WHICH REQUIRE MORE TIME 
        eps=10.*tiny(x1)
C       1. SELF-INTERSECTIONS OF THE POLYGONE SEGMENTS
        np1=np-2
        DO ip=1,np-3                 
         x11=xp(ip)
         y11=yp(ip)
         x12=xp(ip+1)
         y12=yp(ip+1)
         DO ip1=ip+2,np1
          x21=xp(ip1)
          y21=yp(ip1)
          x22=xp(ip1+1)
          y22=yp(ip1+1)
          IF(gridman_intersect2d(x11,y11,x12,y12,
     f                           x21,y21,x22,y22,xi,yi)) THEN
            ierr=100
            WRITE(gridman_unit,*) "ERROR in  GRIDMAN_GRID2D_CUT: ",
     w                            "polygon have self intersection"
            WRITE(gridman_unit,*) " IP, IP1 ",ip,ip1
            RETURN          
          END IF 
         END DO
         np1=np-1
        END DO
C> \todo More efficient algorithm should be invented to check 
C>       for nodes on segments and vertices on edges
C       2. POLYGONE SEGMENTS COME TOO CLOSE TO GRID NODES
        DO ip=1,np-1
         x1=xp(ip)
         y1=yp(ip)
         x2=xp(ip+1)
         y2=yp(ip+1)
         DO ipoint=1,grid%NPOINTS
          x0=grid%X(1,ipoint)*grid%UNIT2SI
          y0=grid%X(2,ipoint)*grid%UNIT2SI
          IF(gridman_close2interval2d(x1,y1,
     f                                x2,y2,x0,y0,gridman_tol)) THEN
           ierr=-100 
           WRITE(gridman_unit,*) "WARNING from GRIDMAN_GRID2D_CUT: ",
     w             "polygone segment comes too close to a grid node"
           WRITE(gridman_unit,*) "     IP, X1, Y1, X2, Y2 ",
     w                                 ip, x1, y1, x2, y2
           WRITE(gridman_unit,*) " IPOINT, X, Y ",ipoint,x0,y0                      
          END IF 
         END DO
        END DO
C       3. GRID EDGES COME TOO CLOSE TO POLYGONE SEGMENTS
        DO iedge=1,grid%NEDGES
         ipoint=grid%POINTS(1,iedge)
         x1=grid%X(1,ipoint)*grid%UNIT2SI
         y1=grid%X(2,ipoint)*grid%UNIT2SI
         ipoint=grid%POINTS(2,iedge)
         x2=grid%X(1,ipoint)*grid%UNIT2SI
         y2=grid%X(2,ipoint)*grid%UNIT2SI
         DO ip=1,np-1
          IF(gridman_close2interval2d(x1,y1,x2,y2,
     w                                xp(ip),yp(ip),gridman_tol)) THEN
           ierr=-100 
           WRITE(gridman_unit,*) "WARNING from GRIDMAN_GRID2D_CUT: ",
     w            "grid edge comes too close to a polygone vertex"
           WRITE(gridman_unit,*) " IEDGE, X1, Y1, X2, Y2 ",
     w                             iedge, x1, y1, x2, y2
           WRITE(gridman_unit,*) "    IP, X, Y ",ip,xp(ip),yp(ip) 
          END IF 
         END DO
        END DO
       
       END IF !IF(GRIDMAN_CHECK) THEN
 
      END SUBROUTINE check_polygon   

C --------------------------------------------------------------------
C     ALLOCATE TEMPORARY ARRAYS
C --------------------------------------------------------------------
      SUBROUTINE allocate_temporary(NPOINTS,NPOL,NPMAX,IERR)        
       INTEGER(GRIDMAN_SP),INTENT(IN) :: NPOINTS,NPOL,NPMAX
       INTEGER,INTENT(OUT) :: IERR
       INTEGER(GRIDMAN_SP) :: NE,NC,ST

       ierr=0

       ne=grid%NEDGES+2*npol 
       nc=grid%NCELLS+npol
 
       IF(gridman_dbg)
     f   WRITE(gridman_unit,*) 
     w        "GRIDMAN_GRID2D_CUT: NPOL, NPOINTS, NPMAX ",
     w                             npol, npoints, npmax
 
       ALLOCATE(x1(npoints),x2(npoints),
     a          ioldedge(ne),ioldcell(nc),
     a          cells(2,ne),points(2,ne),
     a          ipwork(npmax+2),tpwork(npmax+2),dwork(npmax+2),stat=st)
       IF(st.NE.0) THEN
        WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w                       "cannot allocate temporary arrays"
        WRITE(gridman_unit,*) " NPOINTS, NPMAX, NPOL ",
     w                          npoints, npmax, npol
        WRITE(gridman_unit,*) " NCELLS, NEDGES ",grid%NCELLS,grid%NEDGES
        ierr=200
        RETURN
       END IF   

      END SUBROUTINE allocate_temporary 

C --------------------------------------------------------------------
C     DEALLOCATE TEMPORARY ARRAYS
C --------------------------------------------------------------------
      SUBROUTINE deallocate_temporary(IERR)
       INTEGER,INTENT(OUT) :: IERR
       INTEGER :: ST(17)
 
       st=0

       IF(ALLOCATED(x1)) DEALLOCATE(x1,stat=st(1))
       IF(ALLOCATED(x2)) DEALLOCATE(x2,stat=st(2))
       IF(ALLOCATED(ioldedge)) DEALLOCATE(ioldedge,stat=st(3))
       IF(ALLOCATED(ioldcell)) DEALLOCATE(ioldcell,stat=st(4))
       IF(ALLOCATED(cells)) DEALLOCATE(cells,stat=st(5))
       IF(ALLOCATED(points)) DEALLOCATE(points,stat=st(6))
       IF(ALLOCATED(icell_pol)) DEALLOCATE(icell_pol,stat=st(7))
       IF(ALLOCATED(iedge_pol)) DEALLOCATE(iedge_pol,stat=st(8))
       IF(ALLOCATED(ipwork)) DEALLOCATE(ipwork,stat=st(9))
       IF(ALLOCATED(tpwork)) DEALLOCATE(tpwork,stat=st(10))
       IF(ALLOCATED(icount)) DEALLOCATE(icount,stat=st(11))
       IF(ALLOCATED(x_pol)) DEALLOCATE(x_pol,stat=st(12))
       IF(ALLOCATED(y_pol)) DEALLOCATE(y_pol,stat=st(13))
       IF(ALLOCATED(dwork)) DEALLOCATE(dwork,stat=st(14))
       IF(ALLOCATED(lwork)) DEALLOCATE(lwork,stat=st(15))
       IF(ALLOCATED(indseg)) DEALLOCATE(indseg,stat=st(16))

       CALL gridman_indlist_deallocate(celledges,st(17))

       IF(any(st.NE.0)) THEN
         ierr=-200
         WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID2D_CUT: ",
     w         "problems with deallocation of temporary arrays"
       END IF
      

      END SUBROUTINE deallocate_temporary

C --------------------------------------------------------------------
      SUBROUTINE deallocate_all(IERR)
       INTEGER,INTENT(OUT) :: IERR
       CALL deallocate_temporary(ierr)
       CALL gridman_grid_deallocate(cutgrid,ierr)
      END  SUBROUTINE deallocate_all

      END SUBROUTINE gridman_grid2d_cut

C****************************************************************************
C     BUILD TRAJECTORY OF THE POLYGON ON THE GRID -
C     FIND EDGES AND CELL INTERSECTED BY THE POYLGON
C     auxiliary subroutine called from GRIDMAN_GRID2D_CUT
C****************************************************************************
      SUBROUTINE gridman_grid2d_polygon_trajectory(NPOL,NPEDGE,
     s                                             icell_pol,iedge_pol,
     s                                             x_pol,y_pol, 
     s                                             indvert,indseg,
     s                                             xp,yp,np,
     s                                             grid,ierr)
      USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,gridman_unit,
     u                 gridman_dbg,gridman_indlist
      USE gridman_lib,ONLY:gridman_grid2d_chains,gridman_grid_cells
      IMPLICIT NONE
C     INPUT 
      INTEGER(GRIDMAN_SP),INTENT(IN) :: NP !NUMBER OF POLYGON VERTICES 
      REAL(GRIDMAN_DP) :: XP(np),YP(np)    !COORDINATES OF THE POLYGON VERTICES
      TYPE(gridman_grid) :: GRID           !GRID OBJECT  
C     OUTPUT
      INTEGER(GRIDMAN_SP),INTENT(OUT) :: 
     i         npol,               !NUMBER OF POINTS IN THE RESULTING TRAJECTORY
     i         npedge(grid%NEDGES) !NUMBER OF INTERSECTIONS OF THE POLYGON WITH EACH EDGE 
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: 
     i                       icell_pol(:), !INDEX OF CELL WHERE THE POINT LIES
     i                                     !(OR THERE THE TRAJECTORY GOES IN)
     i                       iedge_pol(:), !INDEX OF EDGE ON WHICH THE POINT LIES
     i                                     
     i                       indseg(:)     !INDSEG(IP) IS THE INDEX OF THE POLYGON
     i                                     !SEGMENT ON WHICH THE INTERVAL IP,IP+1 LIES               
      INTEGER(GRIDMAN_SP),INTENT(OUT) ::  INDVERT(np)  !INDEX MAPPING OF THE POLYGON VERTICES
     i                                     !TO THE POINTS X_POL, Y_POL  (OLD->NEW)      
     i                                     !INDEX OF THIS VERTEX, 0 OTHERWISE 
      REAL(GRIDMAN_DP),ALLOCATABLE :: X_POL(:),Y_POL(:) !COORDINATES OF TRAJECTORY POINTS       
      INTEGER,INTENT(OUT) :: IERR       !ERROR CODE
C     LOCAL
      INTEGER(GRIDMAN_SP) :: NWORK,IEDGE0,IEDGE1,ICELL0,ICELL1,IP,
     i                       icell_first,iedge_first,ip_first,ipol
      REAL(GRIDMAN_DP) :: X0,Y0,X1,Y1,X_FIRST,Y_FIRST
C     LOCAL, ALLOCATABLE
      TYPE(gridman_indlist) :: CELLS,CHAINS
      REAL(GRIDMAN_DP),ALLOCATABLE :: XWORK(:),YWORK(:)

C     CHECK INPUT GRID
      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) 
     w        "Starting GRIDMAN_GRID2D_POLYGON_TRAJECTORY"

      ierr=0
   
C     GENERATE CLOSED CHAINS OF POINTS
      CALL gridman_grid2d_chains(grid,chains,ierr)
      IF(ierr.GT.0) RETURN

C     GENERATE LISTS OF EDGED
      CALL gridman_grid_cells(cells,grid,ierr)
      IF(ierr.GT.0) RETURN

C     ALLOCATE BUFFER ARRAYS XWORK, YWORK
      CALL allocate_work(chains,ierr)
      IF(ierr.GT.0) RETURN

C     FIND INITIAL POINT OF THE TRAJECTORY
      CALL first_point(x_first,y_first,
     s                 icell_first,iedge_first,ip_first)
      IF(ierr.GT.0) RETURN

      IF(gridman_dbg) THEN
       WRITE(gridman_unit,*) "GRIDMAN_GRID2D_POLYGON_TRAJECTORY:"
       WRITE(gridman_unit,*) " ICELL_FIRST, IEDGE_FIRST, IP_FIRST ",
     w                         icell_first, iedge_first, ip_first
       WRITE(gridman_unit,*) " X_FIRST, Y_FIRST ", x_first,y_first
      END IF

      npedge=0  
      
      IF(iedge_first.LT.0.AND.icell_first.LT.0) THEN
       npol=0
       RETURN
      END IF

C     Here instead of linked lists a primitive and old fashioned way
C     of dynamic memory allocation is applied. 
C     First the trajectory is run through to calculate its length.
C     After that arrays are allocated, and practically same 
C     calculations is repeated. This time to record the trajectory
C     in the allocated arrays

C     CALCULATE NUMBER OF POINTS NPOL IN THE RESULTING TRAJECTORY
      x0=x_first
      y0=y_first
      ip=ip_first
      icell0=icell_first
      iedge0=iedge_first
      npedge=0  
      IF(iedge0.GT.0) npedge(iedge0)=npedge(iedge0)+1 
      npol=1
      DO !INFINITE CYCLE 
       CALL find_next_point(iedge1,icell1,x1,y1,iedge0,icell0,
     c                      x0,y0,xp(ip),yp(ip))
       IF(iedge1.LT.1.AND.icell1.LT.1) GOTO 300
       x0=x1
       y0=y1
       icell0=icell1
       iedge0=iedge1      
       IF(iedge0.GT.0) npedge(iedge0)=npedge(iedge0)+1 
       IF(iedge0.LT.1) THEN    
C       IF POINT NOT ON THE EDGE THEN THIS IS THE 
C       POLYGON VERTEX, UPDATE INDEX OF THE NEXT VERTEX
        ip=ip+1 
C       IF THIS IS VERTEX NP THEN STOP THE CYCLE
        IF(ip.GT.np) EXIT
       END IF
       npol=npol+1
       IF(icell1.LT.1) THEN
C       TRAJECTORY HAS LEFT THE GRID, LOOK FOR THE POINT
C       WHERE TRAJECTORY ENTERS THE GRID AGAIN 
   11   CONTINUE  
        CALL find_closest_boundary(iedge1,x1,y1,iedge0,
     c                             x0,y0,xp(ip),yp(ip)) 
        IF(iedge1.GT.0) THEN
C        INTERSECTION FOUND, ADD NEXT POINT TO THE POLYGON         
         x0=x1
         y0=y1
         iedge0=iedge1
         npedge(iedge0)=npedge(iedge0)+1 
         icell0=grid%CELLS(1,iedge0)
         IF(icell0.LT.1) icell0=grid%CELLS(2,iedge0)
         npol=npol+1
        ELSE
C        NO INTERSECTION, SHIFT TO THE NEXT VERTEX OF THE POLYGON
         x0=xp(ip)
         y0=yp(ip)  
         iedge0=-1 
         icell0=0
         ip=ip+1
         IF(ip.GT.np) EXIT
         GOTO 11
        END IF !IF(IEDGE1.GT.1) THEN  
       END IF !IF(ICELL1.LT.1) THEN
      END DO !INFINITE CYCLE

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) 
     w        "GRIDMAN_GRID2D_POLYGON_TRAJECTORY: NPOL ",npol

C     ALLOCATE OUTPUT ARRAYS
      CALL allocate_arrays(npol,ierr)        
      IF(ierr.GT.0) RETURN

C     CALCULATE TRAJECTORY
      indvert=0
      x0=x_first
      y0=y_first
      ip=ip_first      
      icell0=icell_first
      iedge0=iedge_first
      ipol=1 
      x_pol(ipol)=x0
      y_pol(ipol)=y0
      icell_pol(ipol)=icell0
      iedge_pol(ipol)=iedge0 
      IF(iedge0.LT.0) THEN
       IF(ip.LT.2) GOTO 100
       indvert(ip-1)=ipol
      END IF   
      indseg(ipol)=ip-1
      DO !INFINITE CYCLE 
       CALL find_next_point(iedge1,icell1,x1,y1,iedge0,icell0,
     c                       x0,y0,xp(ip),yp(ip)) 
       IF(iedge1.LT.1.AND.icell1.LT.1) GOTO 300
       x0=x1
       y0=y1
       icell0=icell1
       iedge0=iedge1
       IF(iedge0.LT.1) THEN    
C       IF POINT NOT ON THE EDGE THEN THIS IS THE 
C       POLYGON VERTEX, UPDATE INDEX OF THE NEXT VERTEX
        indvert(ip)=ipol+1
        ip=ip+1 
C       IF THIS IS VERTEX NP THEN STOP THE CYCLE
        IF(ip.GT.np) EXIT
       END IF
       ipol=ipol+1
       x_pol(ipol)=x0
       y_pol(ipol)=y0
       icell_pol(ipol)=icell0
       iedge_pol(ipol)=iedge0
       indseg(ipol)=ip-1
       IF(icell1.LT.1) THEN
 22     CALL find_closest_boundary(iedge1,x1,y1,iedge0,
     c                                    x0,y0,xp(ip),yp(ip)) 
C       TRAJECTORY HAS LEFT THE GRID, LOOK FOR THE POINT
C       WHERE TRAJECTORY ENTERS THE GRID AGAIN 
        IF(iedge1.GT.0) THEN
         x0=x1
         y0=y1
         iedge0=iedge1
         icell0=grid%CELLS(1,iedge0)
         IF(icell0.LT.1) icell0=grid%CELLS(2,iedge0)
         ipol=ipol+1
         x_pol(ipol)=x0
         y_pol(ipol)=y0
         icell_pol(ipol)=icell0
         iedge_pol(ipol)=iedge0     
         indseg(ipol)=ip-1
        ELSE
C        NO INTERSECTION, SHIFT TO THE NEXT VERTEX OF THE POLYGON
         x0=xp(ip)
         y0=yp(ip)
         iedge0=-1 
         icell0=0
         ip=ip+1
         IF(ip.GT.np) EXIT
         GOTO 22
        END IF       
       END IF !IF(ICELL1.LT.1) THEN
      END DO

      IF(ipol.NE.npol) GOTO 200

      indvert(np)=indvert(1) !FIRST AND LAST POINTS OF THE CONTOUR ARE THE SAME

C     TRANSLATE COORDINATES FROM METER INTO UNITS OF GRID
      x_pol(1:npol)=x_pol(1:npol)/grid%UNIT2SI
      y_pol(1:npol)=y_pol(1:npol)/grid%UNIT2SI

      CALL deallocate_local(ierr)

      IF(gridman_dbg) 
     f   WRITE(gridman_unit,*) 
     w        "GRIDMAN_GRID2D_POLYGON_TRAJECTORY finished"

      RETURN       

 100  WRITE(gridman_unit,*) 
     w     "ERROR in GRIDMAN_GRID2D_POLYGON_TRAJECTORY: ",
     w     "internal error (1)"
      ierr=400
      CALL deallocate_local(ierr)
      RETURN
 200  WRITE(gridman_unit,*) 
     w     "ERROR in GRIDMAN_GRID2D_POLYGON_TRAJECTORY: ",
     w     "internal error (2)"
      ierr=400
      CALL deallocate_local(ierr)
      RETURN
 300  WRITE(gridman_unit,*) 
     w     "ERROR in GRIDMAN_GRID2D_POLYGON_TRAJECTORY: ",
     w     "internal error (3)"
      WRITE(gridman_unit,*) " IEDGE0, ICELL0 ",iedge0,icell0
      WRITE(gridman_unit,*) " X0, Y0 ",x0,y0
      WRITE(gridman_unit,*) " IEDGE1, ICELL1 ",iedge1,icell1
      WRITE(gridman_unit,*) " X1, Y1 ",x1,y1
      WRITE(gridman_unit,*) " IP ",ip
      WRITE(gridman_unit,*) " XP, YP ",xp(ip),yp(ip)
      ierr=400
      CALL deallocate_local(ierr)
      RETURN

      CONTAINS

C --------------------------------------------------------------------  
C     FIND FIRST POINT OF THE TRAJECTORY
C --------------------------------------------------------------------
      SUBROUTINE first_point(X_FIRST,Y_FIRST,
     s                       icell_first,iedge_first,ip_first)
       REAL(GRIDMAN_DP),INTENT(OUT) :: X_FIRST, !COORDINATES
     s                                 y_first 
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: ICELL_FIRST, !INDEX OF CELL 
     s                                    iedge_first, !INDEX OF EDGE (<0 IF INSIDE CELL)
     s                                    ip_first     !INDEX OF THE NEXT POLYGON VERTEX            
C
       INTEGER(GRIDMAN_SP) :: IP,ICELL1,ICELL2,ICELL0,IEDGE,IEDGE0
       REAL(GRIDMAN_DP) :: X0,Y0,X1,Y1              

C      SCALING FACTORS, XP,YP ARE  TRANSLATED INTO UNITS OF GRID
       DO ip=1,np-1
        CALL find_closest_edge(iedge,x1,y1,xp(ip),yp(ip),
     c                                     xp(ip+1),yp(ip+1))
        IF(iedge.GT.0) THEN
         ip_first=ip 
         EXIT
        END IF
       END DO
       IF(iedge.LT.1) THEN
        icell_first=-1
        iedge_first=-1
        RETURN
       END IF  

C      NOW WE HAVE TO DECIDE WHICH POINT HAS TO BE TAKEN AS THE FIRST ONE:
C      VERTEX IP OF THE POLYGON OR THE INTERSECTION POINT
       ip=ip_first
       x0=xp(ip)
       y0=yp(ip)
       iedge0=iedge
       icell1=grid%CELLS(1,iedge) 
       IF(icell1.GT.0) THEN
        IF(point_in_cell(x0,y0,icell1)) THEN
C        IF VERTEX IP IS IN ICELL1, 
C        THEN THE FIRST POINT IS IP, AND FIRST CELL IS ICELL1
         iedge0=-1
         icell0=icell1 
        ELSE
C        IF IP IS NOT IN ICELL1, THEN CHECK ICELL2
         icell2=grid%CELLS(2,iedge) 
         IF(icell2.GT.0) THEN
C         IF VERTEX IP IS IN ICELL2, 
C         THEN THE FIRST POINT IS IP, AND FIRST CELL IS ICELL2
          iedge0=-1
          icell0=icell2 
         ELSE
C         OTHERWISE: ICELL2<=0, ICELL1>0, BUT IP IS NOT IN ICELL1
C         THUS, IP IS OUTSIE THE GRID
C         TAKE INTERSECTION POINT AS THE FIRST POINT,
C         THAT STARTS INTO ICELL1
          x0=x1
          y0=y1
          icell0=icell1
         END IF 
        END IF              
       ELSE
C       ICELL1<=0, CHECK ICELL2 
        icell2=grid%CELLS(2,iedge)           
        IF(icell2.GT.0) THEN
         icell0=icell2         
         IF(point_in_cell(x0,y0,icell2)) THEN        
C         VERTEX IP IS IN ICELL2
C         FIRST POINT IS IP AND FIRST CELL IS ICELL2
          iedge0=-1
         ELSE
C         OTHERWISE, ICELL1<=0, ICELL2>0, BUT IP IS NOT IN ICELL2
C         THUS, IP IS OUTSIDE THE GRID,
C         TAKE INTERSECTION POINT AS THE FIRST POINT,
C         THAT STARTS INTO ICELL2
          x0=x1
          y0=y1
         END IF
        ELSE
C        BOTH ICELL1<=0, ICELL2<=0
C        TAKE INTERSECTION POINT AND SET FURST CELL TO ZERO 
         x0=x1
         y0=y1
         icell0=0
        END IF
       END IF
       ip_first=ip+1
       x_first=x0       
       y_first=y0       
       icell_first=icell0       
       iedge_first=iedge0       
       
      END SUBROUTINE first_point

C --------------------------------------------------------------------
C     FIND EDGE INTERSECTED BY THE INTERVAL
C     [(X1,Y1),(X2,Y2)] CLOSEST TO THE POINT (X1,Y1)
C --------------------------------------------------------------------
      SUBROUTINE find_closest_edge(IEDGE_OUT,X_OUT,Y_OUT,X1,Y1,X2,Y2)
       USE gridman_lib,ONLY:gridman_intersect2d
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: IEDGE_OUT !INDEX OF INTERESECTED EDGE,
                                                    
       REAL(GRIDMAN_DP),INTENT(OUT) :: X_OUT,Y_OUT  !COORDINATES OF THE INTERSECTION POINT
       REAL(GRIDMAN_DP),INTENT(IN) :: X1,Y1,X2,Y2  !COORDINATES OF THE ENDS OF THE INTERVAL

       REAL(GRIDMAN_DP) :: XE1,YE1,XE2,YE2,XI,YI,D,DMIN
       INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT

       ierr=0

       dmin=(x1-x2)**2+(y1-y2)**2
       x_out=x2
       y_out=y2
       iedge_out=-1
       DO iedge=1,grid%NEDGES
        ipoint=grid%POINTS(1,iedge)
        xe1=grid%X(1,ipoint)*grid%UNIT2SI
        ye1=grid%X(2,ipoint)*grid%UNIT2SI
        ipoint=grid%POINTS(2,iedge)
        xe2=grid%X(1,ipoint)*grid%UNIT2SI
        ye2=grid%X(2,ipoint)*grid%UNIT2SI
        IF(gridman_intersect2d(x1,y1,x2,y2,xe1,ye1,xe2,ye2,xi,yi)) THEN
         d=(x1-xi)**2+(y1-yi)**2       
         IF(d.LT.dmin) THEN
          dmin=d
          iedge_out=iedge   
          x_out=xi
          y_out=yi
         END IF
        END IF
       END DO

      END SUBROUTINE find_closest_edge

C --------------------------------------------------------------------
C     FIND BOUNDARY EDGE INTERSECTED BY A PARTICLE
C     STARTED FROM THE POINT (X1,Y1) TOWARDS (X2,Y2)
C     CLOSEST TO THE POINT (X1,Y1)
C -------------------------------------------------------------------- 
      SUBROUTINE find_closest_boundary(IEDGE_OUT,X_OUT,Y_OUT,
     s                                 iedge0,x1,y1,x2,y2)
       USE gridman_lib,ONLY:gridman_cross2d
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: IEDGE_OUT !INDEX OF INTERESECTED EDGE,
                                        
       REAL(GRIDMAN_DP),INTENT(OUT) :: X_OUT,Y_OUT  !COORDINATES OF THE INTERSECTION POINT
       INTEGER(GRIDMAN_SP),INTENT(IN) :: IEDGE0
       REAL(GRIDMAN_DP),INTENT(IN) :: X1,Y1  !COORDINATES OF THE STARTING POINT
       REAL(GRIDMAN_DP),INTENT(IN) :: X2,Y2  !COORDINATES OF THE TARGET POINT

       REAL(GRIDMAN_DP) :: XE1,YE1,XE2,YE2,XI,YI,D,DMIN,VX,VY
       INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT,ICELL1,ICELL2

       ierr=0
        
       vx=x2-x1
       vy=y2-y1
       dmin=(x1-x2)**2+(y1-y2)**2
       x_out=x2
       y_out=y2
       iedge_out=-1
       DO iedge=1,grid%NEDGES
        IF(iedge.EQ.iedge0) cycle
        icell1=grid%CELLS(1,iedge)
        icell2=grid%CELLS(2,iedge)
C       SKIP EDGES WHICH ARE NOT ON THE BOUNDARY
        IF(icell1.GT.0.AND.icell2.GT.0) cycle
        ipoint=grid%POINTS(1,iedge)
        xe1=grid%X(1,ipoint)*grid%UNIT2SI
        ye1=grid%X(2,ipoint)*grid%UNIT2SI
        ipoint=grid%POINTS(2,iedge)
        xe2=grid%X(1,ipoint)*grid%UNIT2SI
        ye2=grid%X(2,ipoint)*grid%UNIT2SI
        IF(gridman_cross2d(x1,y1,vx,vy,xe1,ye1,xe2,ye2,xi,yi)) THEN
         d=(x1-xi)**2+(y1-yi)**2
         IF(d.LT.dmin) THEN
          dmin=d
          iedge_out=iedge   
          x_out=xi
          y_out=yi
         END IF
        END IF
       END DO

      END SUBROUTINE find_closest_boundary

C --------------------------------------------------------------------
C     RETURN TRUE IF POINT (X0,Y0) IS IN ICELL,      
C     AND FALSE OTHERWISE     
C --------------------------------------------------------------------
      FUNCTION point_in_cell(X0,Y0,ICELL)
      USE gridman_lib,ONLY:gridman_point_in_polygon
      REAL(GRIDMAN_DP),INTENT(IN) :: X0,Y0
      INTEGER(GRIDMAN_SP),INTENT(IN) :: ICELL 
      INTEGER(GRIDMAN_SP) :: IP,IPOINT,N
      LOGICAL :: POINT_IN_CELL 

      n=0
      DO ip=chains%IFIRST(icell),chains%ILAST(icell)        
       ipoint=chains%IND(ip)  
       n=n+1
       xwork(n)=grid%X(1,ipoint)*grid%UNIT2SI
       ywork(n)=grid%X(2,ipoint)*grid%UNIT2SI
      END DO
      IF( gridman_point_in_polygon(xwork,ywork,n,x0,y0) ) THEN
       point_in_cell=.true.
      ELSE
       point_in_cell=.false.    
      END IF     

      END FUNCTION point_in_cell

C --------------------------------------------------------------------
C     FIND NEXT POINT OF THE TRAJECTORY
C --------------------------------------------------------------------
      SUBROUTINE find_next_point(IEDGE1,ICELL1,X1,Y1,
     f                           iedge0,icell0,x0,y0,xn,yn)   
       USE gridman_lib,ONLY:gridman_cross2d
C      INDEX OF THE NEXT EDGE AND CELL
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: IEDGE1,ICELL1
C      COORDINATES OF THE NEXT POINT
       REAL(GRIDMAN_DP),INTENT(OUT) :: X1,Y1
C      INDEX OF THE INITIAL EDGE AND CELL
       INTEGER(GRIDMAN_SP),INTENT(IN) :: IEDGE0,ICELL0
C      COORDINATES OF THE INITIAL POINT    
       REAL(GRIDMAN_DP),INTENT(IN) :: X0,Y0
C      CCORDINATES OF THE NEXT VERTEX OF THE POLYGON
       REAL(GRIDMAN_DP),INTENT(IN) :: XN,YN
       INTEGER(GRIDMAN_SP) :: IE,IEDGE,IPOINT,IMIN
       REAL(GRIDMAN_DP) :: VX,VY,D,DMIN,XE1,YE1,XE2,YE2,XI,YI,
     r                     xmin,ymin
       LOGICAL :: LERR

C      DIRECTION VECTOR
       vx=xn-x0
       vy=yn-y0        
       imin=-1         
       dmin=(xn-x0)**2+(yn-y0)**2 !DISTANCE TO THE NEXT VERTEX OF THE POLYGON          
       lerr=.true.
C      CYCLE OVER ALL EDGES OF ICELL0
       DO ie=cells%IFIRST(icell0),cells%ILAST(icell0)
        iedge=cells%IND(ie)
        IF(iedge.EQ.iedge0) cycle
        ipoint=grid%POINTS(1,iedge)
        xe1=grid%X(1,ipoint)*grid%UNIT2SI
        ye1=grid%X(2,ipoint)*grid%UNIT2SI
        ipoint=grid%POINTS(2,iedge)
        xe2=grid%X(1,ipoint)*grid%UNIT2SI
        ye2=grid%X(2,ipoint)*grid%UNIT2SI
        IF(gridman_cross2d(x0,y0,vx,vy,xe1,ye1,xe2,ye2,xi,yi)) THEN
         d=(x0-xi)**2+(y0-yi)**2
         lerr=.false.
         IF(d.LT.dmin) THEN 
C         LOOK FOR INTERSECTION WITH EDGES OF ICELL0 CLOSEST TO X0,Y0
          imin=iedge
          xmin=xi
          ymin=yi  
          dmin=d 
         END IF
        END IF
       END DO
       IF(lerr) THEN  
C       NO INTERSECTION FOUND - ERROR
        iedge1=-1
        icell1=-1
        RETURN
       END IF        
       IF(imin.GT.0) THEN
C       TAKE INTERSECTION POINT AS THE NEXT POINT
C       UPDATE IEDGE1, ICELL1
        iedge1=imin
        icell1=grid%CELLS(1,iedge1)
        IF(icell1.EQ.icell0) icell1=grid%CELLS(2,iedge1)
        x1=xmin
        y1=ymin
       ELSE
C       THEN STAY IN ICELL0, TAKE VERTEX (XN,YN) AS THE NEXT POINT
        iedge1=-1
        icell1=icell0
        x1=xn
        y1=yn
       END IF

      END SUBROUTINE find_next_point

C --------------------------------------------------------------------
C     ALLOCATE OUTPUT ARRAYS
C --------------------------------------------------------------------
      SUBROUTINE allocate_arrays(NPOL,IERR)        
       INTEGER(GRIDMAN_SP),INTENT(IN) :: NPOL
       INTEGER,INTENT(OUT) :: IERR
       INTEGER :: ST

       ierr=0

       ALLOCATE(icell_pol(npol),iedge_pol(npol),indseg(npol),
     w          x_pol(npol),y_pol(npol),stat=st)
       IF(st.NE.0) THEN
        WRITE(gridman_unit,*) 
     w       "ERROR in GRIDMAN_GRID2D_POLYGON_TRAJECTORY: ",
     w       "cannot allocate temporary arrays"
        WRITE(gridman_unit,*) " NPOL ",npol
        ierr=200
        RETURN
       END IF   

      END SUBROUTINE allocate_arrays

C --------------------------------------------------------------------
C     ALLOCATE BUFFER ARRAYS XWORK, YWORK
C --------------------------------------------------------------------
      SUBROUTINE allocate_work(CHAINS,IERR)
       TYPE(gridman_indlist) :: CHAINS
       INTEGER,INTENT(OUT) :: IERR
       INTEGER(GRIDMAN_SP) :: N,IC,ST
C      NWORK = MAXIMIM NUMBER OF POINTS IN ONE CELL   
       nwork=0
       DO ic=1,chains%N
        n=chains%ILAST(ic)-chains%IFIRST(ic)+1
        IF(n.GT.nwork)  nwork=n
       END DO
       ALLOCATE(xwork(nwork),ywork(nwork),stat=st)
       IF(st.NE.0) THEN
        WRITE(gridman_unit,*) 
     w       "ERROR in GRIDMAN_GRID2D_POLYGON_TRAJECTORY: ",
     w       "cannot allocate XWORK, YWORK"
        WRITE(gridman_unit,*) " NWORK ", nwork
        ierr=200
        RETURN
       END IF   
      END SUBROUTINE allocate_work           

C --------------------------------------------------------------------
C     DEALLOCATE TEMPORARY ARRAYS
C --------------------------------------------------------------------
      SUBROUTINE deallocate_local(IERR)
       INTEGER,INTENT(OUT) :: IERR
       INTEGER :: ST(4),IERR0

       st=0
       IF(ALLOCATED(xwork)) DEALLOCATE(xwork,stat=st(1))
       IF(ALLOCATED(xwork)) DEALLOCATE(ywork,stat=st(2))
       CALL gridman_indlist_deallocate(chains,ierr0)
       IF(ierr0.NE.0) st(3)=1
       CALL gridman_indlist_deallocate(cells,ierr0)
       IF(ierr0.NE.0) st(4)=1
       IF(any(st.NE.0)) THEN
         ierr=-200
         WRITE(gridman_unit,*) 
     w        "ERROR in GRIDMAN_GRID2D_POLYGON_TRAJECTORY: ",
     w        "problems with deallocation of temporary arrays"
       END IF
      END  SUBROUTINE deallocate_local

      END SUBROUTINE gridman_grid2d_polygon_trajectory