!!  Copyright (C) Stichting Deltares, 2005-2023.
!!
!!  This file is part of iMOD.
!!
!!  This program 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.
!!
!!  This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!  Contact: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.
!!
MODULE MOD_PKS

USE WINTERACTER
USE RESOURCE
USE IMODVAR
USE MOD_PMANAGER_PAR, ONLY : PBMAN !PRJIDF,PRJNLAY,PBMAN
USE MOD_UTL
USE MOD_IDF


TYPE TPKS
         !LOGICAL :: ACTIVE = .FALSE. 
         !CHARACTER(LEN=MAXLEN) :: TEXT = 'PARALLEL KRYLOV SOLVER PACKAGE'
         !INTEGER :: MXITER = 100     ! MAXIMUM NUMBER OF OUTER ITERATIONS
         !INTEGER :: INNERIT = 20     ! MAXIMUM NUMBER OF INNER ITERATIONS
         !INTEGER :: ISOLVER = 1      ! 1: PCG; 2: BICGSTAB; 3: GMRES
         !INTEGER :: NPC = 2          ! 2: INCOMPLETE LU
         !REAL    :: HCLOSEPKS = 1E-3 ! HEAD CHANGE CRITERION
         !REAL    :: RCLOSEPKS = 100. ! RESIDUAL CRITERON
         !REAL    :: RELAXPKS  = 0.98 ! RELAXATION PARAMETER
         !INTEGER :: IPRPKS = 1       ! PRINTOUT FLAG
         !INTEGER :: MUTPKS = 0       ! PRINT FLAG
         !REAL    :: DAMP  = 1.       ! DAMPING FACTOR
         !! PKSMPI PART
 INTEGER :: NRPROC = 0       ! 0 = SERIAL; > 0 = PARALLEL
         !INTEGER :: PARTOPT = 0      ! 0 = UNIFORM; 1 = RCB
 INTEGER, DIMENSION(:,:), ALLOCATABLE :: PARTMINMAX ! ! IC1, IC2, IR1, IR2
!         CHARACTER(LEN=MAXLEN)                :: LOADFILE ! ! LOADPTR
!         INTEGER :: PRESSAKEY = 0 ! FLAG REQUIRED FOR DEBUGGING WITH VISUAL STUDIO
 END TYPE TPKS
 TYPE(TPKS), PUBLIC, SAVE :: PKS
 INTEGER,PRIVATE ::  GNCOL, GNROW !, GNLAY, GNODES 
 INTEGER,PARAMETER,PRIVATE :: NOVLAPMAX=1 !novlapimpsol  = 1
 INTEGER,PARAMETER,PRIVATE :: STENIMPSOL    = 2
 INTEGER, PARAMETER :: INOVL = 1
 REAL(KIND=DP_KIND), PARAMETER :: TINY = 1.0D-20      
 
CONTAINS
 
 !###======================================================================
 SUBROUTINE PKS_INIT(IU,PRJIDF) !,NLAY)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IU !,NLAY
 TYPE(IDFOBJ),INTENT(IN) :: PRJIDF
 INTEGER :: ICOL, IROW, IC, IR,I,J,IC1,IC2,IR1,IR2
 TYPE(IDFOBJ) :: LOADPTR
 REAL(KIND=DP_KIND) :: CHDADD,X1,X2,Y1,Y2
      
 CALL PKS7MPISETGNODES(PRJIDF%NCOL, PRJIDF%NROW) !, NLAY )

 ALLOCATE(PKS%PARTMINMAX(PBMAN%NRPROC,4)) 
            
 !## uniform
 IF (PBMAN%PARTOPT.EQ.0) THEN
     CALL PKS7MPIPARTDEF(PBMAN%NRPROC, PKS%PARTMINMAX(:,1), PKS%PARTMINMAX(:,2),& 
                                       PKS%PARTMINMAX(:,3), PKS%PARTMINMAX(:,4) )
 ELSEIF(PBMAN%PARTOPT.EQ.1.OR.PBMAN%PARTOPT.EQ.2) THEN ! RCB
   CALL IDFCOPY(PRJIDF,LOADPTR)
   IF (.NOT.IDFREADSCALE(PBMAN%LOADFILE,LOADPTR,7,0,0.0D0,0))THEN
    WRITE(*,'(/1X,A/)') 'ERROR READING '//TRIM(PBMAN%LOADFILE); RETURN
   ENDIF
   !## make sure nodata is threated as zero
   DO IROW = 1, PRJIDF%NROW; DO ICOL = 1, PRJIDF%NCOL
    IF(LOADPTR%X(ICOL,IROW).EQ.LOADPTR%NODATA)LOADPTR%X(ICOL,IROW)=0.0D0
   ENDDO; ENDDO
   LOADPTR%NODATA=0.0D0
   ! ASSUME IBOUND IN CASE OF PARTOPT = 3
   IF (PBMAN%PARTOPT.EQ.2) THEN
      ! STEP 1: SET LOAD = 1 FOR ACTIVE CELLS
      DO IROW = 1, PRJIDF%NROW
         DO ICOL = 1, PRJIDF%NCOL    
            LOADPTR%X(ICOL,IROW) = MIN(LOADPTR%X(ICOL,IROW), 1.0D0) 
         END DO
      END DO
      ! STEP 2: MODIFY FOR CONSTANT HEAD CELLS              
      CHDADD = MINVAL(LOADPTR%X)    
      IF (CHDADD.LT.0.0D0) THEN ! CONSTANT HEAD CELLS DETECTED
         CHDADD = CHDADD - 1.0D0 ! SET VALUE THAT SHOULD BE ADDED (BORDER)
         DO IROW = 1, PRJIDF%NROW
            DO ICOL = 1, PRJIDF%NCOL                   
               IF (LOADPTR%X(ICOL,IROW).GE.1.0D0) THEN
                  IF (ICOL.LT.PRJIDF%NCOL) THEN ! E
                     IC = ICOL + 1
                     IR = IROW
                     IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                        LOADPTR%X(IC,IR) = CHDADD    
                     END IF
                     IF (IROW.GT.1) THEN ! NE
                        IR = IROW-1
                        IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                           LOADPTR%X(IC,IR) = CHDADD    
                        END IF
                     END IF       
                     IF (IROW.LT.PRJIDF%NROW) THEN ! SE
                        IR = IROW+1
                        IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                           LOADPTR%X(IC,IR) = CHDADD    
                        END IF
                     END IF       
                  END IF    
                        
                  IF (ICOL.GT.1) THEN ! W
                     IC = ICOL - 1
                     IR = IROW
                     IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                        LOADPTR%X(IC,IR) = CHDADD    
                     END IF
                     IF (IROW.GT.1) THEN ! NW
                        IR = IROW-1
                        IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                           LOADPTR%X(IC,IR) = CHDADD    
                        END IF
                     END IF       
                     IF (IROW.LT.PRJIDF%NROW) THEN ! SW
                        IR = IROW+1
                        IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                           LOADPTR%X(IC,IR) = CHDADD    
                        END IF
                     END IF       
                  END IF    
                        
                  IF (IROW.GT.1) THEN ! N
                     IC = ICOL
                     IR = IROW-1
                     IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                         LOADPTR%X(IC,IR) = CHDADD    
                     END IF
                  END IF
                       
                  IF (IROW.LT.PRJIDF%NROW) THEN ! S
                     IC = ICOL
                     IR = IROW+1
                     IF (LOADPTR%X(IC,IR).LT.0.0D0) THEN
                         LOADPTR%X(IC,IR) = CHDADD    
                     END IF
                  END IF
               END IF
             END DO
         END DO
         ! REMOVE CONSTANT HEAD CELLS   
         DO IROW = 1, PRJIDF%NROW
            DO ICOL = 1, PRJIDF%NCOL                   
               IF (LOADPTR%X(ICOL,IROW).LT.0.0D0 .AND. LOADPTR%X(ICOL,IROW).GT.CHDADD) THEN
                   LOADPTR%X(ICOL,IROW) = 0.0D0
               END IF 
            END DO      
         END DO            
      END IF
   END IF
   CALL PKS7MPIPARTRCB(PBMAN%NRPROC, PKS%PARTMINMAX(:,1), PKS%PARTMINMAX(:,2),& 
                                     PKS%PARTMINMAX(:,3), PKS%PARTMINMAX(:,4),& 
                                     PRJIDF%NCOL, PRJIDF%NROW, LOADPTR%X, 0.0D0)
   CALL IDFDEALLOCATEX(LOADPTR)
 END IF    

 IF(IU.GT.0)THEN
  DO I=1,PBMAN%NRPROC
   WRITE(IU,'(4(I10,1X))') (PKS%PARTMINMAX(I,J),J=1,4)
  ENDDO
 ELSE
  DO I=1,PBMAN%NRPROC
   WRITE(ABS(IU),*) I
   IC1=PKS%PARTMINMAX(I,1)
   IC2=PKS%PARTMINMAX(I,2)
   IR1=PKS%PARTMINMAX(I,3)
   IR2=PKS%PARTMINMAX(I,4)
   
   CALL IDFGETEDGE(PRJIDF,IR1,IC1,X1,Y1,X2,Y2)
   WRITE(ABS(IU),'(2(F15.7,1X))') X1,Y2

   CALL IDFGETEDGE(PRJIDF,IR1,IC2,X1,Y1,X2,Y2)
!   CALL IDFGETLOC(PRJIDF,IR1,IC2,X,Y)
   WRITE(ABS(IU),'(2(F15.7,1X))') X2,Y2

   CALL IDFGETEDGE(PRJIDF,IR2,IC2,X1,Y1,X2,Y2)
!   CALL IDFGETLOC(PRJIDF,IR2,IC2,X,Y)
   WRITE(ABS(IU),'(2(F15.7,1X))') X2,Y1

   CALL IDFGETEDGE(PRJIDF,IR2,IC1,X1,Y1,X2,Y2)
!   CALL IDFGETLOC(PRJIDF,IR2,IC1,X,Y)
   WRITE(ABS(IU),'(2(F15.7,1X))') X1,Y1

   CALL IDFGETEDGE(PRJIDF,IR1,IC1,X1,Y1,X2,Y2)
!   CALL IDFGETLOC(PRJIDF,IR1,IC1,X,Y)
   WRITE(ABS(IU),'(2(F15.7,1X))') X1,Y2   
   WRITE(ABS(IU),'(A)') 'END'
  ENDDO
  WRITE(ABS(IU),'(A)') 'END'
 ENDIF
 
 DEALLOCATE(PKS%PARTMINMAX)
 
 END SUBROUTINE PKS_INIT
       
 !###======================================================================
 SUBROUTINE PKS7MPISETGNODES( NCOL, NROW) !, NLAY)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NCOL  !< (I) NUMBER OF GLOBAL COLUMNS
 INTEGER, INTENT(IN) :: NROW  !< (I) NUMBER OF GLOBAL ROWS
! INTEGER, INTENT(IN) :: NLAY  !< (I) NUMBER OF GLOBAL LAYS
 
 GNCOL  = NCOL
 GNROW  = NROW
! GNLAY  = NLAY
! GNODES = NCOL*NROW*NLAY
      
 END SUBROUTINE PKS7MPISETGNODES

 !###======================================================================
 SUBROUTINE PKS7MPIPARTDEF(NRPROC, PROC_ICOLMIN, PROC_ICOLMAX,PROC_IROWMIN, PROC_IROWMAX )
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NRPROC
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_ICOLMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_ICOLMAX
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_IROWMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_IROWMAX     
 INTEGER :: IPROC, I, J, N, NC, NRC, NRR, BLKSIZE, NPROC
 INTEGER :: WNCOL, WNROW
 INTEGER, DIMENSION(:), ALLOCATABLE :: DIVIDE
 INTEGER                    :: NRPROC_NCOL !< NUMBER OF COLUMN PROCESSES
 INTEGER                    :: NRPROC_NROW !< NUMBER OF ROW PROCESSES
 INTEGER, DIMENSION(:,:), ALLOCATABLE :: PROC_NCOL   !< NUMBER OF COLUMNS FOR EACH PROCESS
 INTEGER, DIMENSION(:,:), ALLOCATABLE :: PROC_NROW   !< NUMBER OF ROWS FOR EACH PROCESS
! INTEGER, DIMENSION(:,:), ALLOCATABLE :: PROC_JCOLMIN
! INTEGER, DIMENSION(:,:), ALLOCATABLE :: PROC_JCOLMAX
! INTEGER, DIMENSION(:,:), ALLOCATABLE :: PROC_JROWMIN
! INTEGER, DIMENSION(:,:), ALLOCATABLE :: PROC_JROWMAX    

 NPROC = NRPROC 
 ALLOCATE(PROC_NCOL(NPROC,1),PROC_NROW(NPROC,1))

!C...     ALLOCATE WORK ARRAY
 ALLOCATE( DIVIDE(NPROC) )

 DO I = 1, NPROC
    J = NPROC/I
    IF (I*J.EQ.NPROC) THEN
       DIVIDE(I) = 1
    ELSE
       DIVIDE(I) = 0
    END IF
 END DO
      
 WNCOL = MAX(GNCOL-2*NOVLAPMAX,1)
 WNROW = MAX(GNROW-2*NOVLAPMAX,1)
      
 BLKSIZE = SQRT(REAL(WNCOL*WNROW/NPROC))
 N  = MIN( WNCOL, WNROW )
 NC = MAX( 1, INT(N/BLKSIZE) )
 J  = 1

 DO WHILE (DIVIDE(NC).EQ.0)
    IF(DIVIDE(NC+J).EQ.0) THEN
        IF (DIVIDE(NC-J).EQ.0) THEN
           J=J-1
        ELSE
           NC=NC+J
        ENDIF
    ELSE
      NC=NC+J
    ENDIF
 ENDDO

 IF (WNCOL.LT.WNROW) THEN
    NRPROC_NCOL = NC
    NRPROC_NROW = NPROC/NC
 ELSE
    NRPROC_NCOL = NPROC/NC
    NRPROC_NROW = NC
 END IF

 IPROC = 0
 NRR = WNROW
 DO I = NRPROC_NROW, 1, -1
     NRC = WNCOL
     DO J = NRPROC_NCOL, 1, -1
        IPROC = IPROC + 1
        PROC_NCOL(IPROC,INOVL) = INT(NRC/J)
        PROC_NROW(IPROC,INOVL) = INT(NRR/I)
        NRC = NRC - PROC_NCOL(IPROC,INOVL)
     END DO
     NRR = NRR - PROC_NROW(IPROC,INOVL)
 END DO 
      
 DEALLOCATE( DIVIDE )

!C FILL PARTITIONS
 CALL PKS7MPIMINMAX(NRPROC,NRPROC_NROW, NRPROC_NCOL, &
  PROC_ICOLMIN, PROC_ICOLMAX, PROC_IROWMIN, PROC_IROWMAX, &
  PROC_NCOL, PROC_NROW )

!C...     CORRECT FOR PARTITION OVERLAP
 DO IPROC = 1, NRPROC
      PROC_ICOLMIN(IPROC,INOVL) = PROC_ICOLMIN(IPROC,INOVL)+NOVLAPMAX
      PROC_ICOLMAX(IPROC,INOVL) = PROC_ICOLMAX(IPROC,INOVL)+NOVLAPMAX
      PROC_IROWMIN(IPROC,INOVL) = PROC_IROWMIN(IPROC,INOVL)+NOVLAPMAX
      PROC_IROWMAX(IPROC,INOVL) = PROC_IROWMAX(IPROC,INOVL)+NOVLAPMAX
      IF (PROC_ICOLMIN(IPROC,INOVL).EQ.NOVLAPMAX+1)PROC_ICOLMIN(IPROC,INOVL) = 1   
      IF (PROC_ICOLMAX(IPROC,INOVL).EQ.NOVLAPMAX+WNCOL)PROC_ICOLMAX(IPROC,INOVL) = GNCOL
      IF (PROC_IROWMIN(IPROC,INOVL).EQ.NOVLAPMAX+1)PROC_IROWMIN(IPROC,INOVL) = 1
      IF (PROC_IROWMAX(IPROC,INOVL).EQ.NOVLAPMAX+WNROW)PROC_IROWMAX(IPROC,INOVL) = GNROW
 END DO      
      
 DEALLOCATE(PROC_NCOL,PROC_NROW)

 END SUBROUTINE PKS7MPIPARTDEF

 !###======================================================================
 SUBROUTINE PKS7MPIMINMAX(NRPROC,NRPROC_NROW, NRPROC_NCOL, &
   PROC_ICOLMIN, PROC_ICOLMAX, PROC_IROWMIN, PROC_IROWMAX, &
   PROC_NCOL, PROC_NROW ) 
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NRPROC
 INTEGER, INTENT(IN) :: NRPROC_NROW
 INTEGER, INTENT(IN) :: NRPROC_NCOL
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_ICOLMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_ICOLMAX
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_IROWMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_IROWMAX
 INTEGER, DIMENSION(NRPROC,1), INTENT(IN) :: PROC_NCOL
 INTEGER, DIMENSION(NRPROC,1), INTENT(IN) :: PROC_NROW
! INTEGER, DIMENSION(NRPROC,2), INTENT(OUT) :: PROC_ICOLMIN
! INTEGER, DIMENSION(NRPROC,2), INTENT(OUT) :: PROC_ICOLMAX
! INTEGER, DIMENSION(NRPROC,2), INTENT(OUT) :: PROC_IROWMIN
! INTEGER, DIMENSION(NRPROC,2), INTENT(OUT) :: PROC_IROWMAX
! INTEGER, DIMENSION(NRPROC,2), INTENT(IN) :: PROC_NCOL
! INTEGER, DIMENSION(NRPROC,2), INTENT(IN) :: PROC_NROW
 INTEGER :: I, J, NC, NR, IPROC

!C...     DETERMINE START COLUMN/ROW INDICES
 DO I = 1, NRPROC_NROW
    NC = 0
    DO J = 1, NRPROC_NCOL
       IPROC = J + (I-1)*NRPROC_NCOL
       PROC_ICOLMIN(IPROC,INOVL) = NC + 1
       NC = NC + PROC_NCOL(IPROC,INOVL)
       PROC_ICOLMAX(IPROC,INOVL) = NC
    END DO
 END DO
 DO J = 1, NRPROC_NCOL
    NR = 0
    DO I = 1, NRPROC_NROW
       IPROC = J + (I-1)*NRPROC_NCOL
       PROC_IROWMIN(IPROC,INOVL) = NR + 1
       NR = NR + PROC_NROW(IPROC,INOVL)
       PROC_IROWMAX(IPROC,INOVL) = NR
    END DO
 END DO
        
 END SUBROUTINE      
 
 !###======================================================================
 SUBROUTINE PKS7MPIPARTRCB(NRPROC, PROC_ICOLMIN, PROC_ICOLMAX, &
                        PROC_IROWMIN, PROC_IROWMAX, &
                        NCOL, NROW, LOADPTR, NODATA )
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NRPROC
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_ICOLMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_ICOLMAX
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_IROWMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(OUT) :: PROC_IROWMAX
 INTEGER, INTENT(IN) :: NCOL
 INTEGER, INTENT(IN) :: NROW
 REAL(KIND=DP_KIND), DIMENSION(NCOL,NROW), INTENT(IN) :: LOADPTR
 REAL(KIND=DP_KIND), INTENT(IN) :: NODATA
 INTEGER :: IPROC, MAXLEV, ILEV
 INTEGER :: NPROC
! LOGICAL :: OK
! REAL(KIND=DP_KIND), DIMENSION(:,:), ALLOCATABLE :: PARTITIONS 
! INTEGER :: F
 INTEGER, DIMENSION(:), ALLOCATABLE :: NCUT

 NPROC = NRPROC

 ALLOCATE(NCUT(NPROC/2))
 CALL PKS7MPIPRIMEFACTORS(NPROC,NCUT,MAXLEV)
      
 IF (.TRUE.) THEN
  ILEV   = 0
  IPROC  = 0
  CALL PKS7MPIRCB(LOADPTR,NODATA,NCOL,NROW,1,NCOL,1,NROW, &
                  ILEV,MAXLEV,IPROC,NPROC,PROC_ICOLMIN,   &
                  PROC_ICOLMAX,PROC_IROWMIN, PROC_IROWMAX,NCUT)
 END IF

 !IF (.FALSE. .AND. MYRANK.EQ.0) THEN ! DEBUG
 !
 ! PROC_JCOLMIN = 0
 ! PROC_JCOLMAX = 0
 ! PROC_JROWMIN = 0
 ! PROC_JROWMAX = 0
 !     
 ! ILEV   = 0
 ! IPROC  = 0
 ! CALL PKS7MPIRCB(LOADPTR,NODATA,NCOL,NROW,1,NCOL,1,NROW, &
 !  ILEV,MAXLEV,IPROC, NPROC, &
 !  PROC_JCOLMIN, PROC_JCOLMAX, PROC_JROWMIN, PROC_JROWMAX,NCUT)
 !
 ! ALLOCATE(PARTITIONS(NCOL,NROW))
 ! PARTITIONS = -9999.0D0
 ! OK = .TRUE.
 ! DO IPROC = 1, NPROC
 !    IC1 = PROC_JCOLMIN(IPROC,INOVL)
 !    IC2 = PROC_JCOLMAX(IPROC,INOVL)
 !    IR1 = PROC_JROWMIN(IPROC,INOVL)
 !    IR2 = PROC_JROWMAX(IPROC,INOVL)
 !    LOAD = 0
 !    DO IROW = IR1, IR2
 !       DO ICOL = IC1, IC2
 !          IF (PARTITIONS(ICOL,IROW).GT.-9999.0D0) THEN
 !              WRITE(*,*) 'ERROR, NON OVERLAP DETECTED.' 
 !              WRITE(*,*) 'PARTITION OVERWRITTEN', PARTITIONS(ICOL,IROW) 
 !              OK = .FALSE.
 !           ELSE
 !              PARTITIONS(ICOL,IROW) = IPROC
 !              IF (LOADPTR(ICOL,IROW).GT.0.0D0) LOAD = LOAD + 1
 !          END IF 
 !       END DO
 !    END DO
 !    WRITE(*,*) 'LOAD IPROC:',IPROC,LOAD
 ! END DO
 ! WRITE(FNAME,'(A,I3.3,A)') 'PARTITIONS_NP',NPROC,'.ASC'
 !
 ! DEALLOCATE(PARTITIONS) 
 ! ENDIF
  
  DEALLOCATE(NCUT)

  END SUBROUTINE PKS7MPIPARTRCB
      
 !###======================================================================
 RECURSIVE SUBROUTINE PKS7MPIRCB(A,NODATA,NCOL,NROW, &
          IC1,IC2,IR1,IR2,ILEV,MAXLEV,IPROC, &
          NRPROC,PROC_ICOLMIN, PROC_ICOLMAX, PROC_IROWMIN, PROC_IROWMAX,NCUT)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(INOUT) :: ILEV
 INTEGER, INTENT(IN)    :: MAXLEV
 INTEGER, INTENT(INOUT) :: IPROC
 INTEGER, INTENT(IN)    :: NRPROC
 INTEGER, INTENT(IN)    :: NCOL, NROW
 REAL(KIND=DP_KIND), INTENT(IN) :: NODATA
 INTEGER, INTENT(IN)    :: IC1,IC2,IR1,IR2
 REAL(KIND=DP_KIND), DIMENSION(NCOL,NROW), INTENT(IN) :: A 
 INTEGER, DIMENSION(NRPROC,1), INTENT(INOUT) :: PROC_ICOLMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(INOUT) :: PROC_ICOLMAX
 INTEGER, DIMENSION(NRPROC,1), INTENT(INOUT) :: PROC_IROWMIN
 INTEGER, DIMENSION(NRPROC,1), INTENT(INOUT) :: PROC_IROWMAX  
 INTEGER, DIMENSION(MAXLEV), INTENT(IN) :: NCUT   
 INTEGER, DIMENSION(NRPROC) :: CUTIDX
  INTEGER :: IC1B, IC2B, IR1B, IR2B, NCOLB, NROWB, ICUT, &
             IC1C, IC2C, IR1C, IR2C, NCOLC, NROWC, JLEV     
 REAL(KIND=DP_KIND), DIMENSION(NCOL) :: ROWSUM
 REAL(KIND=DP_KIND), DIMENSION(NROW) :: COLSUM 

 ILEV = ILEV + 1 
 IF (ILEV <= MAXLEV) THEN
     ! DETERMINE BOUNDING BOX
     CALL PKS7MPIBBOX(A,NCOL,NROW,IC1,IC2,IR1,IR2, &
                 IC1B,IC2B,IR1B,IR2B,NCOLB,NROWB,NODATA)
     IF (NROWB.GT.NCOLB) THEN ! HORIZONTAL CUT(S)
        CALL PKS7MPICOLSUM(A,NCOL,NROW,IC1B,IC2B,IR1B,IR2B,NROWB,COLSUM,NODATA)
        CALL PKS7MPISPLIT(COLSUM,NROWB,CUTIDX,NCUT(ILEV))
        IR1C = IR1B
        DO ICUT = 1, NCUT(ILEV)
           IR2C = IR1B + CUTIDX(ICUT) - 1
           NROWC = IR2C - IR1C + 1
           JLEV = ILEV
           CALL PKS7MPIRCB(A,NODATA,NCOL,NROW,  &
                      IC1B,IC2B,IR1C,IR2C,      &
                      JLEV,MAXLEV,IPROC,NRPROC, &
                      PROC_ICOLMIN,PROC_ICOLMAX, &
                      PROC_IROWMIN,PROC_IROWMAX,NCUT)
           IR1C = IR2C + 1
        END DO
     ELSE ! VERTICAL CUT(S)
        CALL PKS7MPIROWSUM(A,NCOL,NROW,IC1B,IC2B,IR1B,IR2B,NCOLB,ROWSUM,NODATA)
        CALL PKS7MPISPLIT(ROWSUM, NCOLB, CUTIDX, NCUT(ILEV))
        IC1C = IC1B
        DO ICUT = 1, NCUT(ILEV)
           IC2C = IC1B + CUTIDX(ICUT) - 1
           NCOLC = IC2C - IC1C + 1
           JLEV = ILEV
           CALL PKS7MPIRCB(A,NODATA,NCOL,NROW, &
                      IC1C,IC2C,IR1B,IR2B, &
                      JLEV,MAXLEV,IPROC,NRPROC, &
                      PROC_ICOLMIN,PROC_ICOLMAX, &
                       PROC_IROWMIN,PROC_IROWMAX,NCUT)
           IC1C = IC2C + 1
        END DO    
      END IF
 ELSE
    IPROC = IPROC + 1
    PROC_ICOLMIN(IPROC,INOVL) = IC1
    PROC_ICOLMAX(IPROC,INOVL) = IC2
    PROC_IROWMIN(IPROC,INOVL) = IR1
    PROC_IROWMAX(IPROC,INOVL) = IR2
 END IF

 END SUBROUTINE PKS7MPIRCB

 !###======================================================================
  SUBROUTINE PKS7MPIBBOX(A,NCOL,NROW,IC1I,IC2I,IR1I,IR2I, &
                     IC1B,IC2B,IR1B,IR2B,NCOLO,NROWO,NODATA)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: NCOL, NROW
 REAL(KIND=DP_KIND), DIMENSION(NCOL,NROW), INTENT(IN) :: A
 REAL(KIND=DP_KIND), INTENT(IN) :: NODATA
 INTEGER, INTENT(IN)  :: IC1I, IC2I, IR1I, IR2I
 INTEGER, INTENT(OUT) :: IC1B, IC2B, IR1B, IR2B
 INTEGER, INTENT(OUT) :: NCOLO, NROWO
 INTEGER :: IROW, ICOL

 IC1B = NCOL + 1
 IC2B = 0
 IR1B = NROW + 1
 IR2B = 0

 DO IROW = IR1I, IR2I
    DO ICOL = IC1I, IC2I
       IF (ABS(A(ICOL,IROW)-NODATA).GT.TINY) THEN
          IC1B = MIN(IC1B,ICOL)  
          IC2B = MAX(IC2B,ICOL)  
          IR1B = MIN(IR1B,IROW)  
          IR2B = MAX(IR2B,IROW)  
       END IF
    END DO
 END DO

 NCOLO = IC2B-IC1B+1
 NROWO = IR2B-IR1B+1

 END SUBROUTINE PKS7MPIBBOX
      
 !###======================================================================
 SUBROUTINE PKS7MPIROWSUM(A,NCOL,NROW,IC1,IC2,IR1,IR2,NCOL2,ROWSUM,NODATA)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NCOL, NROW, IC1, IC2, IR1, IR2, NCOL2
 REAL(KIND=DP_KIND), DIMENSION(NCOL,NROW) :: A
 REAL(KIND=DP_KIND), DIMENSION(NCOL2) :: ROWSUM
 REAL(KIND=DP_KIND), INTENT(IN) :: NODATA
 INTEGER :: ICOL, IROW
 REAL(KIND=DP_KIND) :: VAL

 ROWSUM = 0
 DO IROW = IR1, IR2
    DO ICOL = IC1, IC2 
        VAL = A(ICOL,IROW)
        IF (ABS(VAL-NODATA).GT.TINY)ROWSUM(ICOL-IC1+1) = ROWSUM(ICOL-IC1+1) + ABS(VAL)
    END DO
 END DO

 END SUBROUTINE PKS7MPIROWSUM
 
 !###======================================================================
 SUBROUTINE PKS7MPICOLSUM(A,NCOL,NROW, IC1, IC2, IR1, IR2,NROW2,COLSUM,NODATA)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NCOL, NROW, IC1, IC2, IR1, IR2, NROW2
 REAL(KIND=DP_KIND), DIMENSION(NCOL,NROW),INTENT(IN) :: A
 REAL(KIND=DP_KIND), DIMENSION(NROW2),INTENT(OUT) :: COLSUM
 REAL(KIND=DP_KIND), INTENT(IN) :: NODATA
 INTEGER :: ICOL, IROW
 REAL(KIND=DP_KIND) :: VAL

 COLSUM = 0
 DO IROW = IR1, IR2
    DO ICOL = IC1, IC2 
        VAL = A(ICOL,IROW) 
        IF (ABS(VAL-NODATA).GT.TINY)COLSUM(IROW-IR1+1) = COLSUM(IROW-IR1+1) + ABS(VAL)
    END DO
 END DO  

 END SUBROUTINE PKS7MPICOLSUM

 !###======================================================================
 SUBROUTINE PKS7MPISPLIT(X, XN, Y, YN)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: XN
 REAL(KIND=DP_KIND), DIMENSION(XN), INTENT(IN) :: X
 INTEGER, INTENT(IN) :: YN
 INTEGER, DIMENSION(YN), INTENT(OUT) :: Y 
 INTEGER :: IX, IY
 REAL(KIND=DP_KIND) :: SUM, TOT 

 TOT = 0
 DO IX = 1, XN
    TOT = TOT + X(IX)   
 END DO
 TOT = TOT/YN
 
 Y = XN
 DO IY = 1, YN-1
    SUM = 0.0D0
    DO IX = 1, XN-1
       IF ((SUM + X(IX+1)) > TOT*IY) THEN
          Y(IY) = IX; EXIT       
       ELSE
          SUM = SUM + X(IX)
       END IF 
    END DO
 END DO
 
 END SUBROUTINE PKS7MPISPLIT   

 !###======================================================================
 SUBROUTINE PKS7MPIPRIMEFACTORS(NUM, FACTORS, F)
 !###======================================================================
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NUM
 INTEGER,INTENT(OUT), DIMENSION((NUM/2))::FACTORS 
 INTEGER, INTENT(INOUT) :: F
 INTEGER :: I, N

 I = 2  
 F = 1
 N = NUM 
 DO
    IF (MOD(N,I) == 0) THEN
        FACTORS(F) = I
        F = F+1
        N = N/I
    ELSE
        I = I+1
    END IF
    IF (N == 1) THEN       
        F = F-1
        EXIT
    END IF
  END DO
 
  END SUBROUTINE
      
  END MODULE MOD_PKS