cutgrid.f

Go to the documentation of this file.

C> @file apps/cutgrid.f
C> Cut part of the grid inside or outside a contour
C>
C> Based on GRIDMAN library
C 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/>.

      PROGRAM cutgrid

      USE gridman,ONLY:gridman_grid,gridman_dp,gridman_sp,
     u                 gridman_unit,gridman_tol
      USE gridman_lib,ONLY:gridman_grid_read,gridman_grid_write,
     u                     gridman_grid2d_cut,gridman_grid_deallocate
      IMPLICIT NONE
      INTRINSIC len_trim,get_command_argument,adjustl,trim    

      INTERFACE
C
       SUBROUTINE cutgrid_read_contour(CONTOUR_IN,UNIT2METER,XP,YP,NP)
        USE gridman,ONLY:gridman_dp,gridman_sp
        CHARACTER(LEN=256),INTENT(IN) :: contour_in
        REAL(GRIDMAN_DP),INTENT(IN) :: unit2meter
        INTEGER(GRIDMAN_SP),INTENT(OUT) :: np
        REAL(GRIDMAN_DP),ALLOCATABLE :: xp(:),yp(:)   
       END SUBROUTINE cutgrid_read_contour
C
       SUBROUTINE cutrgrid_mergeindex(GRID,IEIND)
        USE gridman,ONLY:gridman_grid
        TYPE(gridman_grid) :: grid
        INTEGER,INTENT(IN) :: ieind
       END  SUBROUTINE cutrgrid_mergeindex

      END INTERFACE 

C     INPUT PARAMETERS
      CHARACTER(LEN=256) :: description,
     c                      grid_in,contour_in,grid_out
      LOGICAL :: lexclude,dbgmod,lcheck
      REAL(GRIDMAN_DP) :: unit2meter,tol
C     LOCAL
      TYPE(gridman_grid) :: grid0,grid1
      INTEGER(GRIDMAN_SP) :: np
      INTEGER ::ieind 
      INTEGER :: lt,ierr
      REAL(GRIDMAN_DP),ALLOCATABLE :: xp(:),yp(:) 
      CHARACTER(LEN=8) :: ftmp

      CALL get_command_argument(1,ftmp)
      IF(len_trim(ftmp).GT.0) THEN
       IF(adjustl(trim(ftmp)).EQ.'help') THEN
        CALL print_description 
        stop
       END IF
      ELSE
C     1. READ INPUT FROM NAMELIST
       CALL read_input
      END IF

C    2. READ INPUT GRID
      lt=len_trim(grid_in)
      IF(lt.GT.0) THEN
        WRITE(gridman_unit,*) "CUTGRID: reading file ",grid_in(1:lt)
        CALL gridman_grid_read(grid0,grid_in,ierr)
        IF(ierr.GT.0) CALL cutgrid_error("")
      ELSE 
       CALL cutgrid_error("GRID_IN is not defined")
      END IF 

C    3. READ CONTOUR
      CALL cutgrid_read_contour(contour_in,unit2meter,xp,yp,np)

C    4. PERFORM CUT
      WRITE(gridman_unit,*) "CUTGRID: cut the grid"
      CALL gridman_grid2d_cut(grid1,grid0,xp,yp,np,lexclude,ierr)
      IF(ierr.GT.0) CALL cutgrid_error("")
      lt=len_trim(description)
      IF(lt.GT.0) grid1%DESCRIPTION=description

C    5. MERGE INDEXES
      IF(ieind.GT.0.AND.
     f   grid1%NEDGEINDEX.GT.grid0%NEDGEINDEX) THEN
       WRITE(gridman_unit,*) "CUTGRID: merging edge indexes"
       CALL cutrgrid_mergeindex(grid1,ieind)
      END IF

C    6. WRITE GRID   
      lt=len_trim(grid_out)
      IF(lt.GT.0) THEN
       WRITE(gridman_unit,*) "CUTGRID: writing file ",grid_out(1:lt)
       CALL gridman_grid_write(grid1,grid_out,ierr)
       IF(ierr.NE.0) CALL cutgrid_error("")
      END IF

      CALL gridman_grid_deallocate(grid0,ierr)
      CALL gridman_grid_deallocate(grid1,ierr)

      WRITE(gridman_unit,*) "CUTGRID COMPLETED"

      CONTAINS

C********************************************************************
C Read namelist from standard input
C********************************************************************
       SUBROUTINE read_input
       namelist /cutgrid/ description,grid_in,grid_out,contour_in,
     n                    lexclude,ieind,dbgmod,lcheck,tol,unit2meter
      INTEGER :: io      

      description=''
      grid_in=''
      contour_in=''
      grid_out=''
      lexclude=.false.
      ieind=-1
      dbgmod=.false.
      lcheck=.false.
      tol=-1.
      unit2meter=-1.

      rewind(5)
      READ(5,nml=cutgrid,iostat=io)
      IF(io.NE.0) CALL cutgrid_error("can not read namelist CUTGRID")

      IF(tol.GT.0) gridman_tol=tol      

      END SUBROUTINE read_input

C
      SUBROUTINE print_description
       USE gridman,ONLY:gridman_dbg
       INTRINSIC index,len,trim
       INTEGER,PARAMETER :: file_length=1024
       CHARACTER(LEN=FILE_LENGTH) :: path
       CHARACTER(LEN=128) :: str
       INTEGER :: i,io
       CALL get_command_argument(0,path)
       i=index(path,'/',.true.) 
       IF(i.GT.len(path)) THEN 
        path=''
       ELSE 
        path=path(1:i)
       END IF

       IF(gridman_dbg)  WRITE(gridman_unit,*) " PATH:",trim(path)

       OPEN(unit=3,file=trim(path)//'cutgrid.parameters.description',
     o      status='OLD',iostat=io) 
       IF(io.NE.0) GOTO 100

       DO
        READ(3,'(A)',iostat=io,end=200) str
        IF(io.NE.0) GOTO 200
        WRITE(*,*) trim(str)
       END DO
 200   RETURN

 100   WRITE(gridman_unit,*)
     w      "Could not find cutgrid.parameters.description",
     w      trim(path)
       WRITE(gridman_unit,*)
     w  "Use 'which cutgrid' to invoke via the full path"

       END SUBROUTINE print_description

      END PROGRAM cutgrid

C***************************************************************************************
C Read contour
C***************************************************************************************
      SUBROUTINE cutgrid_read_contour(CONTOUR_IN,UNIT2METER,XP,YP,NP)
       USE gridman,ONLY:gridman_dp,gridman_sp,gridman_unit,gridman_tol
       USE gridman_lib,ONLY:gridman_template_read
       IMPLICIT NONE
       INTRINSIC len_trim,abs,tiny

       CHARACTER(LEN=256),INTENT(IN) :: contour_in
       REAL(GRIDMAN_DP),INTENT(IN) :: unit2meter
       INTEGER(GRIDMAN_SP),INTENT(OUT) :: np
       REAL(GRIDMAN_DP),ALLOCATABLE :: xp(:),yp(:) 

       INTEGER(GRIDMAN_SP),ALLOCATABLE :: n(:)
       REAL(GRIDMAN_DP),ALLOCATABLE :: x0(:),y0(:)
       REAL(GRIDMAN_DP) :: dx,dy
       INTEGER(GRIDMAN_SP) :: m,l,n0
       INTEGER :: lt,st,ierr

       lt=len_trim(contour_in)
       IF(lt.GT.0) THEN
        WRITE(gridman_unit,*) "CUTGRID: reading file ",contour_in(1:lt)
        CALL gridman_template_read(contour_in,m,n,l,x0,y0,ierr)
        IF(ierr.GT.0) CALL cutgrid_error("")
        IF(m.GT.1) THEN
         WRITE(gridman_unit,*) "WARNING from CUTGRID:"
         WRITE(gridman_unit,*) 
     w        " more then one contour is defined in the file ",
     w         contour_in(1:lt)
         WRITE(gridman_unit,*) " First contour will be taken" 
        END IF        
        n0=n(1) 
        dx=gridman_tol*(abs(x0(1))+abs(x0(n0)))+10.*tiny(dx) 
        dy=gridman_tol*(abs(y0(1))+abs(y0(n0)))+10.*tiny(dy) 
        IF( abs(x0(1)-x0(n0)).GT.dx.OR.
     f      abs(y0(1)-y0(n0)).GT.dy) THEN 
C        THE CONTOUR IS NOT CLOSED, ONE MORE POINT WILL BE ADDED 
         np=n0+1
        ELSE
         np=n0
        END IF
        ALLOCATE(xp(np),yp(np),stat=st)  
        IF(st.NE.0) THEN
         WRITE(gridman_unit,*)  " NP ",np          
         CALL cutgrid_error("Cannot allocate temporary arrays")
        END IF
        IF(unit2meter.GT.0.) THEN
         xp(1:n0)=x0(1:n0)*unit2meter
         yp(1:n0)=y0(1:n0)*unit2meter
        ELSE
         WRITE(gridman_unit,*) "WARNING from CUTGRID: ",
     w                         " UNIT2METER is not defined"
         WRITE(gridman_unit,*) 
     w  " I assume that the coordinates in file ",contour_in(1:lt)
         WRITE(gridman_unit,*) " are defined in Millimeter"
         xp(1:n0)=x0(1:n0)*1e-3
         yp(1:n0)=y0(1:n0)*1e-3
        END IF
        IF(np.GT.n0) THEN
C        CLOSE THE CONTOUR
         xp(np)=xp(1)
         yp(np)=yp(1)
         WRITE(gridman_unit,*) 
     w        "CUTGRID: point added to close the contour"
        END IF
        DEALLOCATE(n,x0,y0)
       ELSE
        CALL cutgrid_error("CONTOUR_IN is not defined")
       END IF !IF(LT.GT.0) THEN

      END SUBROUTINE cutgrid_read_contour


C***************************************************************************************
C Merge selected edge index into indexes of the polygon segments
C***************************************************************************************
      SUBROUTINE cutrgrid_mergeindex(GRID,IEIND)
       USE gridman,ONLY:gridman_grid,gridman_index,gridman_unit
       USE gridman_lib,ONLY:gridman_index_allocate,gridman_index_merge
       IMPLICIT NONE 
       INTRINSIC trim
       TYPE(gridman_grid) :: grid
       INTEGER,INTENT(IN) :: ieind

       TYPE(gridman_index) :: ind1,indseg,ind0,
     t                        indtmp(grid%NEDGEINDEX)
       INTEGER :: ierr

C      After call of cutgrid the last index is the 
C      index of the polygon segments, otherwise
C      this subroutine is not called at all 
       indseg=grid%EDGEINDEX(grid%NEDGEINDEX)
       IF(ieind.LT.1.OR.ieind.GT.grid%NEDGEINDEX) THEN
        WRITE(gridman_unit,*) " IEIND, NEDGEINDEX ",
     w                          ieind, grid%NEDGEINDEX
        CALL cutgrid_error("IEIND is out of range")
       END IF
       ind0=grid%EDGEINDEX(ieind)
       IF(ind0%NINDEX.NE.1)  THEN
        WRITE(gridman_unit,*) " NINDEX ",ind0%NINDEX
        CALL cutgrid_error("NINDEX of IEIND must be 1")
       END IF 

       CALL gridman_index_merge(ind1,indseg,ind0,ierr)
       IF(ierr.NE.0)  CALL cutgrid_error("")

C      Move indices
       indtmp=grid%EDGEINDEX
       DEALLOCATE(grid%EDGEINDEX)
       ALLOCATE(grid%EDGEINDEX(grid%NEDGEINDEX+1))
       grid%EDGEINDEX(1)=ind1
       grid%EDGEINDEX(2:grid%NEDGEINDEX+1)=indtmp
       grid%NEDGEINDEX=grid%NEDGEINDEX+1

      END SUBROUTINE cutrgrid_mergeindex

C***************************************************************************************
C Print error message and exit
C***************************************************************************************
       SUBROUTINE cutgrid_error(STR)           
        USE gridman,ONLY:gridman_unit
        IMPLICIT NONE
        CHARACTER(*),INTENT(IN) :: str
        INTEGER :: lt
        INTRINSIC len_trim
        lt=len_trim(str)
        IF(lt.GT.0) 
     w     WRITE(gridman_unit,*) "ERROR in CUTGRID: "//str(1:lt)
        stop "ERROR in CUTGRID - see log output"
       END SUBROUTINE cutgrid_error