init_index_mpi.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE init_index_mpi (YSC, &
                           hprogram,halg,pio_frac,oshadows)
!
! 04-2014 : Modifs Matthieu Lafaysse for shadows :
!        * OSHADOWS logical to pass to get_sizes_parallel
!        * disactivate the treatments relative to PIO_FRAC in case of shadows
!        * MPI_BCAST of NIX and NIY
!
USE modd_surfex_n, ONLY : surfex_t
!
USE modd_surfex_mpi, ONLY : nrank, npio, nproc, ncomm, nindex, nsize_task, nsize, wlog_mpi
USE modd_surfex_omp, ONLY : nindx2sfx, nwork, nwork2, xwork, xwork2, xwork3, &
                            nwork_full, nwork2_full, xwork_full, xwork2_full
!
USE modd_mask, ONLY : nmask_full
!
USE modd_slope_effect, ONLY : nix,niy
!
USE mode_read_surf_cov, ONLY : read_surf_cov
!
USE modi_get_luout
USE modi_set_surfex_filein
USE modi_init_io_surf_n
USE modi_read_surf
USE modi_read_gridtype
USE modi_end_io_surf_n
USE modi_abor1_sfx
USE modi_get_sizes_parallel
USE modi_get_adjacent_meshes
USE modi_get_size_full_n
USE modi_ini_data_cover
USE modi_read_arrange_cover
USE modi_read_cover_garden
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!
TYPE(surfex_t), INTENT(INOUT) :: ysc
!
 CHARACTER(LEN=6), INTENT(IN) :: hprogram
 CHARACTER(LEN=4), INTENT(IN) :: halg
REAL, INTENT(IN) :: pio_frac
LOGICAL, INTENT(IN) :: oshadows
!
INTEGER, DIMENSION(:), ALLOCATABLE :: ileft   ! index of left   grid mesh 
INTEGER, DIMENSION(:), ALLOCATABLE :: iright  ! index of right  grid mesh 
INTEGER, DIMENSION(:), ALLOCATABLE :: itop    ! index of top    grid mesh 
INTEGER, DIMENSION(:), ALLOCATABLE :: ibottom ! index of bottom grid mesh
!
INTEGER, DIMENSION(0:NPROC-1) :: inbpts
INTEGER :: ireste, inreste, irank, iproc
INTEGER :: jj, cpt, iluout, inbmin, ip0
INTEGER :: iresp, infompi
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI',0,zhook_handle)
!
IF ( nrank==npio ) THEN
  !
  CALL get_luout(hprogram,iluout)
  !
  !*    1. Parameters of the grid
  !
  CALL set_surfex_filein(hprogram,'PGD ')
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
  !
  CALL read_surf(&
                 hprogram,'DIM_FULL  ',ysc%U%NDIM_FULL,iresp,hdir='A')
  nindx2sfx = ysc%U%NDIM_FULL
  !
  CALL read_arrange_cover(&
                          hprogram,ysc%U%LWATER_TO_NATURE,ysc%U%LTOWN_TO_ROCK,'A')
  CALL read_cover_garden(&
                         hprogram,ysc%U%LGARDEN,'A')
  !
  CALL end_io_surf_n(hprogram)
  DEALLOCATE(nmask_full)
  !
  CALL set_surfex_filein(hprogram,'PGD ') ! change input file name to pgd name
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
  !
  CALL read_surf(&
                 hprogram,'GRID_TYPE',ysc%UG%CGRID,iresp,hdir='A')
  !
  CALL read_gridtype(&
                     hprogram,ysc%UG%CGRID,ysc%UG%NGRID_PAR,ysc%U%NDIM_FULL,.false.,hdir='A')
  !
  ALLOCATE(ysc%UG%XGRID_FULL_PAR(ysc%UG%NGRID_PAR))
  !
  CALL read_gridtype(&
                     hprogram,ysc%UG%CGRID,ysc%UG%NGRID_PAR,ysc%U%NDIM_FULL,.true.,&
                     ysc%UG%XGRID_FULL_PAR,iresp,hdir='A')
  !
  !
  CALL end_io_surf_n(hprogram)
  !
  ALLOCATE(nwork_full(ysc%U%NDIM_FULL))
  ALLOCATE(xwork_full(ysc%U%NDIM_FULL))
  ALLOCATE(nwork2_full(ysc%U%NDIM_FULL,10))
  ALLOCATE(xwork2_full(ysc%U%NDIM_FULL,10))
  !
ELSE
  ALLOCATE(nwork_full(0))
  ALLOCATE(xwork_full(0))
  ALLOCATE(nwork2_full(0,0))
  ALLOCATE(xwork2_full(0,0))
ENDIF
!

IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_bcast(ysc%U%NDIM_FULL,kind(ysc%U%NDIM_FULL)/4,mpi_integer,npio,ncomm,infompi)
#endif
 nindx2sfx = ysc%U%NDIM_FULL
ENDIF
!
!*     3. Points by task
!
IF (.NOT.ALLOCATED(nindex)) ALLOCATE(nindex(ysc%U%NDIM_FULL))
nindex(:) = -1
!
IF (nrank==npio) THEN
  !
  IF (halg=='LIN ') THEN
    !
    CALL set_nb_points_lin(nproc,nproc-1,ysc%U%NDIM_FULL,nindex,oshadows)
    !
  ELSEIF (halg=='ADJ ' .OR. halg=='TILA' .OR. halg=='TILL') THEN
    !
    IF (halg=='ADJ ' .OR. halg=='TILA') THEN
      ALLOCATE(ileft(ysc%U%NDIM_FULL))
      ALLOCATE(iright(ysc%U%NDIM_FULL))
      ALLOCATE(ibottom(ysc%U%NDIM_FULL))
      ALLOCATE(itop(ysc%U%NDIM_FULL))
      CALL get_adjacent_meshes(ysc%UG%CGRID,ysc%UG%NGRID_PAR,ysc%U%NDIM_FULL,&
                                ysc%UG%XGRID_FULL_PAR,ileft,iright,itop,ibottom)
    ELSE
      ALLOCATE(ileft(0))
      ALLOCATE(iright(0))
      ALLOCATE(ibottom(0))
      ALLOCATE(itop(0))
    ENDIF
    !
    IF (halg=='ADJ ') THEN
      CALL set_nb_points_adj(nproc,ysc%U%NDIM_FULL,ysc%U%NDIM_FULL,ileft,iright,itop,ibottom,nindex,oshadows)
    ELSEIF (halg=='TILA' .OR. halg=='TILL') THEN
      CALL ini_data_cover(ysc%DTCO, ysc%U)
      CALL set_nb_points_til(hprogram,halg,nproc,ysc%U%NDIM_FULL,ileft,iright,itop,ibottom,nindex,oshadows)
    ENDIF
    !
    DEALLOCATE(ileft  )
    DEALLOCATE(iright )
    DEALLOCATE(itop   )
    DEALLOCATE(ibottom)
    !
  ELSE
    !
    CALL abor1_sfx("INIT_INDEX_MPI: ALG="//halg//" not defined for the moment")
    !
  ENDIF
  !
  IF (nproc>1) THEN
    !
    inbpts(:) = 0
    DO jj=1,SIZE(nindex)
      inbpts(nindex(jj)) = inbpts(nindex(jj)) + 1
    ENDDO
    !   
    inbmin = minval(inbpts)
    ip0 = maxval(minloc(inbpts)) - 1
    !   
    IF (.NOT. oshadows) THEN
    ! Matthieu Lafaysse :
    ! With shadows we don't want the repartition of points to be modified by the following instructions
    
      DO WHILE( inbpts(npio) > nint(pio_frac*inbmin) )
        !
        DO jj=1,SIZE(nindex)
          IF (nindex(jj)==npio) THEN
            nindex(jj) = ip0
            inbpts(npio) = inbpts(npio) - 1
            inbpts(ip0) = inbpts(ip0) + 1
            EXIT
          ENDIF
        ENDDO
      !
        inbmin = maxval(inbpts)
        ip0 = 0
        DO jj=0,nproc-1
          IF (jj/=npio .AND. inbpts(jj)<inbmin) THEN
            inbmin = inbpts(jj)
            ip0 = jj
          ENDIF
        ENDDO
      !
      END DO
    
    END IF
    !
  ENDIF
  !
ENDIF
!

IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_bcast(nindex,SIZE(nindex)*kind(nindex)/4,mpi_integer,npio,ncomm,infompi) 
  IF (oshadows) THEN
  ! Matthieu Lafaysse :
  ! Each thread need to know NIX and NIY
    CALL mpi_bcast(nix,kind(nix)/4,mpi_integer,npio,ncomm,infompi)
    CALL mpi_bcast(niy,kind(niy)/4,mpi_integer,npio,ncomm,infompi)
  ENDIF
#endif
ENDIF
!
ALLOCATE(nsize_task(0:nproc-1))
nsize_task(:) = 0
!
 CALL get_size_full_n(ysc%U, &
                      hprogram,ysc%U%NDIM_FULL,nsize_task(nrank))
nindx2sfx = nsize_task(nrank)
!
nsize = 0
IF (nproc>1) THEN
  DO jj=0,nproc-1
#ifdef SFX_MPI
    CALL mpi_bcast(nsize_task(jj),kind(nsize_task)/4,mpi_integer,jj,ncomm,infompi)
#endif
    IF ( nsize_task(jj)>nsize ) nsize = nsize_task(jj)
    CALL wlog_mpi('SIZE_TASK ',klog=jj,klog2=nsize_task(jj))
  ENDDO
ELSE
  nsize = nsize_task(0)
ENDIF
!
ysc%UG%XGRID_FULL_PAR=>null()
!
ALLOCATE(nwork(nsize))
ALLOCATE(xwork(nsize))
ALLOCATE(nwork2(nsize,10))
ALLOCATE(xwork2(nsize,10))
ALLOCATE(xwork3(nsize,10,10))
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI',1,zhook_handle)
!
 CONTAINS
!
!***************************************************************
!
SUBROUTINE set_nb_points_lin(KPROC,KPROCMIN,KSIZE,KINDEX,OSHADOWS)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kproc
INTEGER, INTENT(IN) :: kprocmin
INTEGER, INTENT(IN) :: ksize
INTEGER, DIMENSION(:), INTENT(INOUT) :: kindex
LOGICAL, INTENT(IN) :: oshadows
!
INTEGER, DIMENSION(0:KPROC-1) :: isize_task
INTEGER :: ji, jj, cpt, iproc1
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_LIN',0,zhook_handle)
!
 CALL get_sizes_parallel(ysc%DTCO, ysc%DGU, ysc%UG, ysc%U, &
                         kproc,ksize,kprocmin,isize_task,oshadows)
!
iproc1 = kprocmin
!
 cpt   = 0
!
DO jj=1,SIZE(kindex) ! boucle sur les points du domaine
  !
  IF (kindex(jj)==-1) THEN ! si le point doit être affecté à cette itération
    cpt = cpt + 1 ! on augmente le nombre de points affectés de 1
    DO WHILE( cpt>isize_task(iproc1) ) ! si on est hors les bornes permises par le proc en cours
      IF ( iproc1.GE.kprocmin ) THEN ! d'abord, on va de IPROCMIN à IPROC-1
        iproc1 = iproc1+1
        IF ( iproc1==kproc ) iproc1 = kprocmin-1 ! une fois qu'on est à IPROC-1,
      ELSE                                       ! on redescend sous IPROC1
        iproc1 = iproc1-1
      ENDIF
      cpt = 1
    ENDDO
    kindex(jj) = iproc1
  ENDIF
  !
ENDDO
!

IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_LIN',1,zhook_handle)
!
END SUBROUTINE set_nb_points_lin
!
!**************************************************************************
!
SUBROUTINE set_nb_points_adj(KPROC,KSIZE,KSIZE_FULL,KLEFT,KRIGHT,KTOP,KBOTTOM,KINDEX,OSHADOWS)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kproc
INTEGER, INTENT(IN) :: ksize
INTEGER, INTENT(IN) :: ksize_full
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: kleft
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: kright
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: ktop 
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: kbottom
INTEGER, DIMENSION(KSIZE_FULL), INTENT(INOUT) :: kindex
LOGICAL, INTENT(IN) :: oshadows
!
INTEGER, DIMENSION(KSIZE_FULL) :: ipoint
INTEGER, DIMENSION(0:KPROC-1) :: isize_task
INTEGER, DIMENSION(4) :: inear
INTEGER :: cpt_tot, cpt_loc, cpt_inter, ji, jj, jk
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_ADJ',0,zhook_handle)
!
 CALL get_sizes_parallel(ysc%DTCO, ysc%DGU, ysc%UG, ysc%U, &
                         kproc,ksize,npio,isize_task,oshadows)
!
DO ji=1,kproc-1
  !
  cpt_inter = 0
  cpt_loc   = 0
  cpt_tot   = 0
  !
  DO WHILE ( cpt_tot < isize_task(ji) )
    !
    IF (cpt_loc < 1 ) THEN ! if no free point has been found in neighbours
      !
      DO jj=1,ksize_full
        IF ( kindex(jj)==-1 ) THEN
          cpt_tot = cpt_tot + 1
          cpt_loc = 1
          ipoint(1) = jj
          kindex(jj) = ji
          EXIT
        ENDIF
      ENDDO
      !
    ENDIF
    !
    ipoint(1:cpt_loc-cpt_inter) = ipoint(cpt_inter+1:cpt_loc)
    cpt_loc   = cpt_loc - cpt_inter
    cpt_inter = cpt_loc
    !
    b1 : DO jj=1,cpt_inter
      !
      inear(1) = kbottom(ipoint(jj))
      inear(2) = kleft(ipoint(jj))
      inear(3) = kright(ipoint(jj)) 
      inear(4) = ktop(ipoint(jj))    
      !
      DO jk=1,4
        !
        IF ( inear(jk).NE.0 ) THEN
          IF( kindex(inear(jk)).EQ.-1 ) THEN
            cpt_tot = cpt_tot + 1
            cpt_loc = cpt_loc + 1
            IF (cpt_tot.GT.isize_task(ji)) EXIT b1
            kindex(inear(jk)) = ji
            ipoint(cpt_loc) = inear(jk)
          ENDIF
        ENDIF
        !
      ENDDO
      !
    ENDDO b1
    !
  ENDDO
  !
ENDDO
!
WHERE (kindex(:)==-1) kindex(:) = 0
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_ADJ',1,zhook_handle)
!
END SUBROUTINE set_nb_points_adj
!
!*****************************************************************
!
SUBROUTINE set_nb_points_til(HPROGRAM,HALG,KPROC,KSIZE_FULL,KLEFT,KRIGHT,KTOP,KBOTTOM,KINDEX,OSHADOWS)
!
! modif 
! 05/14 B. Decharme Partion done according to patch and not vegtype (vegtype can be > to patch grid)
!
USE modd_data_cover_par, ONLY : nvegtype, jpcover
USE modd_surf_par,       ONLY : xundef
!
USE modi_av_pgd
USE modi_surf_patch
!
IMPLICIT NONE
!
 CHARACTER(LEN=6), INTENT(IN) :: hprogram
 CHARACTER(LEN=4), INTENT(IN) :: halg
INTEGER, INTENT(IN) :: kproc
INTEGER, INTENT(IN) :: ksize_full
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: kleft
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: kright
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: ktop 
INTEGER, DIMENSION(KSIZE_FULL), INTENT(IN) :: kbottom
INTEGER, DIMENSION(KSIZE_FULL), INTENT(INOUT) :: kindex
LOGICAL, INTENT(IN) :: oshadows
!
REAL, DIMENSION(KSIZE_FULL) :: zsea
REAL, DIMENSION(KSIZE_FULL) :: zwater
REAL, DIMENSION(KSIZE_FULL) :: znature
REAL, DIMENSION(KSIZE_FULL) :: ztown
REAL, DIMENSION(KSIZE_FULL, JPCOVER)  :: zcover
REAL, DIMENSION(KSIZE_FULL, NVEGTYPE) :: zvegtype
INTEGER, DIMENSION(KSIZE_FULL,2) :: itype 
INTEGER,DIMENSION(4,17) :: itype0
INTEGER, DIMENSION(KSIZE_FULL) :: iindex
INTEGER, DIMENSION(0:KPROC-1) :: isize_task
INTEGER, DIMENSION(0:KPROC-1) :: inbpts
INTEGER, DIMENSION(1) :: iprocmin
INTEGER :: ifull, iresp, jj, ji, jk, jvegtype, cpt, in1, in2
 CHARACTER(LEN=6) :: ynature
LOGICAL, DIMENSION(JPCOVER) :: gcover
LOGICAL :: gdata_vegtype
!
 CHARACTER(LEN=12) :: yrecfm         ! Name of the article to be read
INTEGER           :: iversion, ipatch
REAL, DIMENSION(:,:), ALLOCATABLE :: zpatch
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_TIL',0,zhook_handle)
!
gdata_vegtype = .false.
!
! Full read
!
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
!
 CALL read_surf(&
                 hprogram,'FRAC_SEA   ',zsea,   iresp, hdir='A')
 CALL read_surf(&
                 hprogram,'FRAC_NATURE',znature,iresp, hdir='A')
 CALL read_surf(&
                 hprogram,'FRAC_WATER ',zwater, iresp, hdir='A')
 CALL read_surf(&
                 hprogram,'FRAC_TOWN  ',ztown,  iresp, hdir='A')
 CALL read_surf(&
                 hprogram,'NATURE     ',ynature,iresp, hdir='A')
 CALL read_surf(&
                 hprogram,'VERSION    ',iversion,iresp,hdir='A')
!
 CALL end_io_surf_n(hprogram)
!
! Nature read
!
IF (ynature/='NONE' .AND. sum(znature)>0.) THEN
!
  !CALL SET_SURFEX_FILEIN(HPROGRAM,'PGD')
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'NATURE','SURF  ','READ ')
!
  yrecfm='PATCH_NUMBER'
  CALL read_surf(&
                 hprogram,yrecfm,ipatch,iresp,hdir='A')
  ALLOCATE(zpatch(ksize_full,ipatch))
  zpatch(:,:)=0.0
!  
  IF (iversion<8) CALL read_surf(&
                 hprogram,'L_VEGTYPE',gdata_vegtype,iresp,hdir='A')
!
ELSE
!
  ALLOCATE(zpatch(ksize_full,1))
  zpatch(:,:)=0.0
!
ENDIF
!
IF (iversion >= 8 .AND. ynature/='NONE' .AND. sum(znature)>0.) THEN
!
   CALL end_io_surf_n(hprogram)
!
   CALL set_surfex_filein(hprogram,'PREP')
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'NATURE','SURF  ','READ ')
   yrecfm = 'PATCH'
   CALL read_surf(&
                 hprogram,yrecfm,zpatch(:,:),iresp,hdir='A')
   WHERE(zpatch(:,:)==xundef)zpatch=0.0
   CALL end_io_surf_n(hprogram)
!
ELSE
!
  IF (gdata_vegtype) THEN
    CALL read_surf(&
                 hprogram,'D_VEGTYPE',zvegtype(:,:),iresp,hdir='A')
    CALL end_io_surf_n(hprogram)
  ELSE
    CALL end_io_surf_n(hprogram)
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')    
    CALL read_surf(&
                 hprogram,'COVER_LIST',gcover(:),iresp,hdir='A')
    CALL read_surf_cov(&
                       hprogram,'COVER',zcover(:,:),gcover,iresp,hdir='A')
    DO jvegtype=1,nvegtype
      CALL av_pgd(ysc%DTCO, &
                   zvegtype(:,jvegtype),zcover,ysc%DTCO%XDATA_VEGTYPE(:,jvegtype),'NAT','ARI',gcover)
    END DO
  ENDIF
!
  CALL end_io_surf_n(hprogram)
!
  CALL surf_patch(ipatch,zvegtype,zpatch)
!
ENDIF
!
! Algo
!
itype(:,:) = 0
!
DO jj=1,ksize_full
  !
  IF (zsea(jj)/=0. .AND. znature(jj)/=0. .AND. zwater(jj)/=0. .AND. ztown(jj)/=0.) THEN ! 1 2 3 4
    itype(jj,1) = 1 
  ELSEIF (zsea(jj)/=0. .AND. znature(jj)/=0. .AND. zwater(jj)/=0.) THEN ! 1 2 3
    itype(jj,1) = 2
  ELSEIF (zsea(jj)/=0. .AND. znature(jj)/=0. .AND. ztown(jj)/=0.) THEN ! 1 2 4
    itype(jj,1) = 3
  ELSEIF (zsea(jj)/=0. .AND. zwater(jj)/=0. .AND. ztown(jj)/=0.) THEN ! 1 3 4 
    itype(jj,1) = 4
  ELSEIF (znature(jj)/=0. .AND. zwater(jj)/=0. .AND. ztown(jj)/=0.) THEN ! 2 3 4
    itype(jj,1) = 5
  ELSEIF (zsea(jj)/=0.) THEN 
    IF (znature(jj)/=0.) THEN ! 1 2
      itype(jj,1) = 6
    ELSEIF (zwater(jj)/=0.) THEN ! 1 3
      itype(jj,1) = 7
    ELSEIF (ztown(jj)/=0.) THEN ! 1 4
      itype(jj,1) = 8
    ELSE ! 1
      itype(jj,1) = 12
    ENDIF
  ELSEIF (znature(jj)/=0.) THEN
    IF (zwater(jj)/=0.) THEN ! 2 3
      itype(jj,1) = 9
    ELSEIF (ztown(jj)/=0.) THEN ! 2 4
      itype(jj,1) = 10
    ELSE ! 2
      itype(jj,1) = 13
    ENDIF
  ELSEIF (zwater(jj)/=0.) THEN
    IF (ztown(jj)/=0.) THEN ! 3 4
      itype(jj,1) = 11
    ELSE ! 3
      itype(jj,1) = 14
    ENDIF
  ELSE ! 4
    itype(jj,1) = 15
  ENDIF
  !
ENDDO
!
!we give numbers to the 40 types of vegtypes
 cpt = 0
DO ji = 1,SIZE(itype0,1)
  DO jj = 1,SIZE(itype0,2)
    cpt = cpt + 1
    itype0(ji,jj) = cpt
  ENDDO
ENDDO
!
DO jj=1,ksize_full
  in1 = 0
  in2 = 0
  DO ji=1,min(3,ipatch)
    IF (zpatch(jj,ji)/=0.) in1 = in1 +1
  ENDDO
  DO ji=4,ipatch
    IF (zpatch(jj,ji)/=0.) in2 = in2 +1
  ENDDO
  itype(jj,2) = itype0(in1+1,in2+1)
ENDDO
!
kindex(:) = -1
!
DO ji = 1,15
  !
  DO jj = 1,68
    !
    !count of the points in the couple of types
    iindex(:) = -2
    DO jk=1,ksize_full
      IF (itype(jk,1)==ji .AND. itype(jk,2)==jj) iindex(jk) = -1
    ENDDO
    !
    ifull = count(iindex(:)==-1)
    !
    IF (ifull.NE.0) THEN
      !
      inbpts(:) = 0
      DO jk=1,SIZE(kindex)
        IF (kindex(jk)>-1) inbpts(kindex(jk)) = inbpts(kindex(jk)) + 1
      ENDDO
      IF (kproc>1) THEN
        iprocmin = minloc(inbpts(0:kproc-1))-1
      ELSE
        iprocmin(:) = 0
      ENDIF
      !
      IF (halg=='TILL') THEN
        CALL set_nb_points_lin(kproc,iprocmin(1),ifull,iindex(:),oshadows)
      ELSEIF (halg=='TILA') THEN
        CALL set_nb_points_adj(kproc,ifull,ksize_full,kleft,kright,ktop,kbottom,iindex(:),oshadows)
      ENDIF
      !
      DO jk=1,ksize_full
        IF ( iindex(jk)> -1 ) kindex(jk) = iindex(jk)
      ENDDO
      !
   
    ENDIF
    !
  ENDDO
  !
ENDDO
!
DEALLOCATE(zpatch)
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_TIL',1,zhook_handle)
!
END SUBROUTINE set_nb_points_til
!
!
!
END SUBROUTINE init_index_mpi