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 (DTCO, U, UG, GCP, HPROGRAM,HINIT,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_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_grid_conf_proj_n, ONLY : grid_conf_proj_t 
!
USE modd_xios, ONLY : lxios
!
USE modd_sfx_oasis, ONLY : loasis
!
USE modd_surfex_mpi, ONLY : nrank, npio, nproc, ncomm, nindex, nsize_task, nsize, wlog_mpi, &
                            nreq, nnum, ndim_full_init
!
USE modn_io_offline, ONLY : lland_use
!
USE modd_surf_conf,       ONLY : cprogname
USE modd_mask, ONLY : nmask_full
!
USE modd_slope_effect, ONLY : nix,niy
!
USE mode_read_surf_cov, ONLY : read_surf_cov
!
USE modi_ini_csts
USE modi_pgd_grid
USE modi_get_luout
USE modi_set_surfex_filein
USE modi_init_io_surf_n
USE modi_read_surf
USE modi_make_choice_array
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_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(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
 CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM
 CHARACTER(LEN=3), INTENT(IN) :: HINIT 
 CHARACTER(LEN=4), INTENT(INOUT) :: HALG
REAL, INTENT(IN) :: PIO_FRAC
LOGICAL, INTENT(IN), OPTIONAL :: 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) :: INB
LOGICAL :: GSHADOWS, GSPEC
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)
!
gspec = (hinit=='PGD' .OR. (hinit=='PRE' .AND. loasis) .OR. &
        (hinit=='OFF' .AND. (lland_use .OR. lxios .OR. loasis)))
!
ALLOCATE(nreq((nproc-1)*2))
!
IF (PRESENT(oshadows)) THEN
  gshadows = oshadows
ELSE
  gshadows = .false.
ENDIF
!
IF ( nrank==npio ) THEN
  !
  CALL get_luout(hprogram,iluout)
  !
  !*    1. Parameters of the grid
  !
  IF (hinit=='PGD') THEN
    !
    cprogname=hprogram
    !
    CALL ini_csts
    !
    CALL pgd_grid (ug, u, gcp, hprogram,'                            ','      ',&
                   .false.,'A')
    !
  ELSE
    !
    CALL set_surfex_filein(hprogram,'PGD ')
    CALL init_io_surf_n(dtco, u, hprogram,'FULL  ','SURF  ','READ ')
    !
    CALL read_surf(hprogram, 'DIM_FULL  ', u%NDIM_FULL, iresp, hdir='A')
    !
    CALL read_arrange_cover(hprogram, u%LWATER_TO_NATURE, u%LTOWN_TO_ROCK,'A')
    CALL read_cover_garden(hprogram, 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(dtco, u, hprogram,'FULL  ','SURF  ','READ ')
    !
    CALL read_surf(hprogram, 'GRID_TYPE', ug%G%CGRID, iresp, hdir='A')
    !
    CALL read_gridtype(hprogram, ug%G%CGRID, ug%NGRID_FULL_PAR, u%NDIM_FULL, .false., hdir='A')
    !
    ALLOCATE(ug%XGRID_FULL_PAR(ug%NGRID_FULL_PAR))
    !
    CALL read_gridtype(hprogram, ug%G%CGRID, ug%NGRID_FULL_PAR, u%NDIM_FULL, .true.,&
                       ug%XGRID_FULL_PAR, iresp ,hdir='A')
    !
    !
    CALL end_io_surf_n(hprogram)
    !
  ENDIF
  ! 
  nsize = u%NDIM_FULL
  !
ELSE
  u%NDIM_FULL = 0
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI
  IF (gspec) THEN
    CALL mpi_bcast(ug%NGRID_FULL_PAR,kind(ug%NGRID_FULL_PAR)/4,mpi_integer,npio,ncomm,infompi)
    IF (nrank/=npio) ALLOCATE(ug%XGRID_FULL_PAR(ug%NGRID_FULL_PAR))
    CALL mpi_bcast(ug%XGRID_FULL_PAR,SIZE(ug%XGRID_FULL_PAR)*kind(ug%XGRID_FULL_PAR)/4,mpi_real,npio,ncomm,infompi)
  ENDIF
  CALL mpi_bcast(u%NDIM_FULL,kind(u%NDIM_FULL)/4,mpi_integer,npio,ncomm,infompi)
#endif
ENDIF
!
ndim_full_init = u%NDIM_FULL
!
!*     3. Points by task
!
IF (nrank==npio .OR. gspec) THEN
  IF (.NOT.ALLOCATED(nindex)) ALLOCATE(nindex(u%NDIM_FULL))
  nindex(:) = -1
ELSEIF (.NOT.ALLOCATED(nindex)) THEN
  ALLOCATE(nindex(0))
ENDIF
!
IF (hinit=='PGD' .AND. (halg=='TILA' .OR. halg=='TILL')) THEN
  halg = 'LIN '
  IF (nrank==npio) THEN
    WRITE(*,*) 'INIT_INDEX_MPI: for PGD, TILA and TILL are forbidden, forced to LIN'
    WRITE(iluout,*) 'INIT_INDEX_MPI: for PGD, TILA and TILL are forbidden, forced to LIN'
  ENDIF
ENDIF
!
IF (halg=='LIN ') THEN
  !
  IF (nrank==npio) CALL set_nb_points_lin(nproc,nproc-1,u%NDIM_FULL,nindex,oshadows)
  !
ELSEIF (halg=='ADJ ' .OR. halg=='TILA' .OR. halg=='TILL') THEN
  !
  IF (nrank==npio .AND. (halg=='ADJ ' .OR. halg=='TILA')) THEN
    ALLOCATE(ileft(u%NDIM_FULL))
    ALLOCATE(iright(u%NDIM_FULL))
    ALLOCATE(ibottom(u%NDIM_FULL))
    ALLOCATE(itop(u%NDIM_FULL))
    CALL get_adjacent_meshes(ug%G%CGRID,ug%NGRID_FULL_PAR,u%NDIM_FULL,&
                             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
    IF (nrank==npio) CALL set_nb_points_adj(nproc,u%NDIM_FULL,u%NDIM_FULL,&
                                            ileft,iright,itop,ibottom,nindex,oshadows)
  ELSEIF (halg=='TILA' .OR. halg=='TILL') THEN
    CALL set_nb_points_til(hprogram,halg,nproc,u%NDIM_FULL,ileft,iright,itop,ibottom,nindex,oshadows)
  ENDIF
  !
ELSE
  !
  CALL abor1_sfx("INIT_INDEX_MPI: ALG="//halg//" not defined for the moment")
  !
ENDIF
!
IF (gspec) THEN
  ALLOCATE(nnum(u%NDIM_FULL))
  nnum(:) = 0
ENDIF  
!
IF (nrank==npio) THEN
  !
  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. gshadows) 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
  !
  IF (gspec) THEN
    inb(:) = 0
    DO jj=1,u%NDIM_FULL
      inb(nindex(jj)) = inb(nindex(jj))+1
      nnum(jj) = inb(nindex(jj))
    ENDDO
  ENDIF
  !  
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI
  IF (gspec) THEN
    CALL mpi_bcast(nindex,SIZE(nindex)*kind(nindex)/4,mpi_integer,npio,ncomm,infompi) 
    CALL mpi_bcast(nnum,SIZE(nnum)*kind(nnum)/4,mpi_integer,npio,ncomm,infompi)
  ENDIF
  IF (gshadows) 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
!
IF (nrank==npio) THEN
  DO jj=1,u%NDIM_FULL
    nsize_task(nindex(jj)) = nsize_task(nindex(jj)) + 1
  ENDDO
ENDIF
!
nsize = 0
IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_bcast(nsize_task,SIZE(nsize_task)*kind(nsize_task)/4,mpi_integer,npio,ncomm,infompi)
#endif
  nsize = maxval(nsize_task)
  DO jj=0,nproc-1
    CALL wlog_mpi('SIZE_TASK ',klog=jj,klog2=nsize_task(jj))
  ENDDO
ELSE
  nsize = nsize_task(0)
ENDIF
!
IF (.NOT.gspec) ug%XGRID_FULL_PAR=>null()
!
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(dtco, ug, 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,KDIM_FULL,KLEFT,KRIGHT,KTOP,KBOTTOM,KINDEX,OSHADOWS)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KPROC
INTEGER, INTENT(IN) :: KSIZE
INTEGER, INTENT(IN) :: KDIM_FULL
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KLEFT
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KRIGHT
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KTOP 
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KBOTTOM
INTEGER, DIMENSION(KDIM_FULL), INTENT(INOUT) :: KINDEX
LOGICAL, INTENT(IN) :: OSHADOWS
!
INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK
INTEGER, DIMENSION(KDIM_FULL) :: IPOINT
INTEGER, DIMENSION(0:KPROC-1) :: ISIZE_TASK
INTEGER, DIMENSION(4) :: INEAR
INTEGER :: ICOUNT,ICPT,IDEB
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(dtco, ug, u, &
                         kproc,ksize,npio,isize_task,oshadows)
!
icount = count(kindex(:)==-1)
ALLOCATE(imask(icount))
icpt = 0
DO jj=1,kdim_full
  IF (kindex(jj)==-1) THEN
    icpt = icpt + 1
    imask(icpt) = jj
  ENDIF
ENDDO
ideb = 1
!
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=ideb,icount
        IF ( kindex(imask(jj))==-1 ) THEN
          cpt_tot = cpt_tot + 1
          cpt_loc = 1
          ipoint(1) = imask(jj)
          kindex(imask(jj)) = ji
          EXIT
        ENDIF
      ENDDO
      ideb = jj+1
      !
    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
!
DEALLOCATE(imask)
!
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,KDIM_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 modn_io_offline, ONLY : lwr_vegtype
!
USE modd_data_cover_par, ONLY : nvegtype, jpcover
!
USE modd_surf_par,       ONLY : xundef
!
USE modi_av_pgd
!
IMPLICIT NONE
!
 CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM
 CHARACTER(LEN=4), INTENT(IN) :: HALG
INTEGER, INTENT(IN) :: KPROC
INTEGER, INTENT(IN) :: KDIM_FULL
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KLEFT
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KRIGHT
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KTOP 
INTEGER, DIMENSION(KDIM_FULL), INTENT(IN) :: KBOTTOM
INTEGER, DIMENSION(KDIM_FULL), INTENT(INOUT) :: KINDEX
LOGICAL, INTENT(IN) :: OSHADOWS
!
REAL, DIMENSION(KDIM_FULL) :: ZSEA
REAL, DIMENSION(KDIM_FULL) :: ZWATER
REAL, DIMENSION(KDIM_FULL) :: ZNATURE
REAL, DIMENSION(KDIM_FULL) :: ZTOWN
REAL, DIMENSION(:,:), ALLOCATABLE  :: ZCOVER
REAL, DIMENSION(:,:), ALLOCATABLE :: ZVEGTYPE
INTEGER, DIMENSION(KDIM_FULL,2) :: ITYPE 
INTEGER, DIMENSION(KDIM_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
INTEGER :: IDIM_NATURE
 CHARACTER(LEN=6) :: YNATURE
LOGICAL, DIMENSION(:), ALLOCATABLE :: GCOVER
LOGICAL :: GDATA_VEGTYPE, GDIM, GDIM2
!
 CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
INTEGER           :: IVERSION, IBUGFIX
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!
IF (lhook) CALL dr_hook('INIT_INDEX_MPI:SET_NB_POINTS_TIL',0,zhook_handle)
!
gdata_vegtype = .false.
!
! Full read
!
 CALL set_surfex_filein(hprogram,'PGD ')

 CALL init_io_surf_n(dtco, 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='-')
 CALL read_surf(hprogram,'VERSION    ',iversion,iresp, hdir='-')
 CALL read_surf(hprogram,'BUG        ',ibugfix,iresp, hdir='-')
 CALL read_surf(hprogram,'DIM_NATURE ',idim_nature,iresp,hdir='-')
!
gdim = (iversion>8 .OR. (iversion==8 .AND. ibugfix>=1))
!
IF (gdim) THEN
  CALL read_surf(hprogram,'ECOSG',u%LECOSG,iresp,hdir='-')
ELSE
  u%LECOSG = .false.
ENDIF
!
 CALL end_io_surf_n(hprogram)
!
 CALL ini_data_cover(dtco, u)
!
! Algo
!
IF (nrank==npio) THEN
  !
  itype(:,:) = 0
  !
  DO jj=1,kdim_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
  !
ENDIF  
!
! Nature read
!
IF (ynature/='NONE' .AND. idim_nature>0) THEN
  !
  ALLOCATE(zvegtype(kdim_full,nvegtype))
  !
  CALL init_io_surf_n(dtco, u, hprogram,'NATURE','SURF  ','READ ')
  CALL read_surf(hprogram,'L_VEGTYPE',gdata_vegtype,iresp,hdir='-')
  CALL end_io_surf_n(hprogram)
  !
  CALL init_io_surf_n(dtco, u, hprogram,'FULL  ','SURF  ','READ ')
  !  
  IF (gdata_vegtype) THEN
    IF (gdim) THEN
      yrecfm='D_VEGTY_'
    ELSE
      yrecfm='D_VEGTYPE'
    ENDIF        
    gdim2 = gdim 
    IF (gdim) CALL read_surf(hprogram,'SPLIT_PATCH',gdim2,iresp) 
    CALL make_choice_array(hprogram, nvegtype, gdim2, yrecfm, zvegtype(:,:),hdir='A')
  ELSE
    IF (lwr_vegtype) THEN
      CALL read_surf(hprogram,'VEGTYPE',zvegtype(:,:),iresp,hdir='E')
    ELSE
      ALLOCATE(gcover(jpcover))
      CALL read_surf(hprogram,'COVER_LIST',gcover(:),iresp,hdir='-')
      IF (nrank==npio) THEN
        ALLOCATE(zcover(kdim_full,count(gcover)))
      ELSE
        ALLOCATE(zcover(0,0))
      ENDIF
      CALL read_surf_cov(hprogram,'COVER',zcover(:,:),gcover,iresp,hdir='E')

      IF (nrank==npio) THEN
        DO jvegtype=1,nvegtype
          CALL av_pgd(dtco,zvegtype(:,jvegtype),zcover,dtco%XDATA_VEGTYPE(:,jvegtype),'NAT','ARI',gcover)
        END DO
      ENDIF
      DEALLOCATE(zcover)
      DEALLOCATE(gcover)

    ENDIF

  ENDIF
!
  CALL end_io_surf_n(hprogram)
!
  IF (nrank==npio) THEN
    !we give numbers to the 40 types of vegtypes
    DO jj=1,kdim_full
      in1 = count(zvegtype(jj,1:3)/=0.)
      in2 = count(zvegtype(jj,4:nvegtype)/=0.)
      itype(jj,2) = in2*4 + in1 + 1
    ENDDO
  ENDIF
  !
ELSE
  !
  IF (nrank==npio) itype(:,2) = 1
  !
ENDIF
!
!
!we give numbers to the 40 types of vegtypes
!
IF (nrank==npio) THEN
  !
  kindex(:) = -1
  !
  DO ji = 1,15
    !
    DO jj = 1,maxval(maxval(itype,2))
      !
      !count of the points in the couple of types
      iindex(:) = -2
      DO jk=1,kdim_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,kdim_full,kleft,kright,ktop,kbottom,iindex(:),oshadows)
        ENDIF
        !
        DO jk=1,kdim_full
          IF ( iindex(jk)> -1 ) kindex(jk) = iindex(jk)
        ENDDO
        !
      ENDIF
      !
    ENDDO
    !
  ENDDO
  !
ENDIF
!
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