MODULE interpweight
!=================================================================================================================================
! MODULE       : interpweight
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF       ! Specific module for the interpolation in ORCHIDEE
!!
!!\n DESCRIPTION: L. Fita, October 2016
!!   This module is used in the re-organization done in the section of the code
!! responsible of the interpolation of the morphology datat (PFT, albedo, LAI,
!! ...) from different files done by L. Fita in the mark of the HEAT project.
!! Previous version of the interpolation of these files, proceed by three 
!! basic steps common for all the files:
!!   1.: Modification of values from file
!!   2.: Creation of land/sea mask using values from file
!!   3.: Interpolate values
!!   This module servs to generalize all these steps and construct a common 
!! interface independently of the file. Each step is generalize and the 
!! different processes almost specific for each file in the previous form
!! are introduced as different methodologies. Thus now it provides:
!!   1. Modification of initial values: Removing that wrong values from the file
!|     `interpweight_modifying_input[1/2/3/4]D': subroutines to modify 1D, 2D, 3D
!!     or 4D input data from file, using different methods
!!   2. Masking input: Compute ’on fly’ the values of the land/sea mask taking
!!    the values in the file
!!     `interpweight_masking_input[1/2/3/4]D': subroutines to compute the mask
!!     using data from input file using different methods: nomask, mbelow, 
!!     mabove, msumrange, var
!!   3. Area weights: Get coincident areas from input file to the target 
!!    projection using aggregate_p
!!      No changes on it
!!   4. Interpolate: Use obtained areas and interpolate values at each grid
!!    point at the same format as it will be provided by XIOS
!!      `interpweight_provide_fractions[1/2/3/4]D': Perform interpolation for
!!      fraction/category data
!!      `interpweight_provide_interpolation[2/4]D': Perform interpolation for
!!      continuous data
!!      `variableusetypes': Variable to provide the values along the additional
!!      dimension to perform the interpolation. (e.g.: in case the 13 pfts
!!      are: 1,3,6,18,23,34,35,39,48,...)
!!   5. A new variable is added which explains the ‘availability’ of data to 
!!    perform the interpolation (0, 1). When it is negative it means that there
!!    was no data for that grid point
!!            availability = SUM(area_source)/Area_target_grid_point
!!   This module will disappear once interpolation will be done throuhguot XIOS
!! 
!! More details at:
!!   https://forge.ipsl.jussieu.fr/orchidee/wiki/DevelopmentActivities/ORCHIDEE-DYNAMICO
!!
!! RECENT CHANGE(S): None
!!
!! REFERENCE(S): None
!!
!! SVN         :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/trunk/ORCHIDEE/src_global/interpweight.f90 $
!! $Date: 2023-07-07 11:17:43 +0200 (ven. 07 juil. 2023) $
!! $Revision: 8086 $
!! \n
!_================================================================================================================================

  USE ioipsl_para
  USE interpol_help
  USE grid

  IMPLICIT NONE

  PRIVATE
  PUBLIC interpweight_1D, interpweight_2D, interpweight_3D, interpweight_4D, interpweight_2Dcont,     &
    interpweight_4Dcont, interpweight_get_varNdims_file, interpweight_RangeR,                         &
    interpweight_get_var2dims_file, interpweight_get_var3dims_file, interpweight_get_var4dims_file,   &
    interpweight_Index2DLonLat, interpweight_ValVecR

  CONTAINS

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_1D
!!
!>\BRIEF        ! Interpolate any input file to the grid of the model for a 1D-variable  
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): outvar1D, aoutvar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_1D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours,       &
    contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype,               &
    maskvalues, maskvarname, dim1, dim2, initime, typefrac,                                           &
    maxresollon, maxresollat, outvar1D, aoutvar)

    USE ioipsl
    USE constantes_var
    USE mod_orchidee_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !! 0.1 Input variables
    !
    INTEGER(i_std), INTENT(in)                 :: nbpt                    !! Number of points for which the data 
                                                                          !!   needs to be interpolated
    INTEGER(i_std), INTENT(in)                 :: Nvariabletypes          !! Number of types of the variable
    REAL(r_std), DIMENSION(Nvariabletypes),    &                          !! Vector of values of the types
      INTENT(in)    :: variabletypes                                      !!   (-1, not used)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo                    !! Vector of latitude and longitudes 
                                                                          !! (beware of the order = 1 : latitude, 
                                                                          !!   2 : longitude)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution              !! The size in km of each grid-box in X and Y
    !
    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours    !! Vector of neighbours for each grid point 
                                                                          !!   1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) 
                            
    REAL(r_std), DIMENSION(nbpt),  INTENT(in)  :: contfrac                !! Fraction of land in each grid box.
    CHARACTER(LEN=80), INTENT(in)              :: filename                !! name of the input map to read

    CHARACTER(LEN=80), INTENT(in)              :: varname                 !! name of the variable to interpolate
    CHARACTER(LEN=80), INTENT(in)              :: inlonname, inlatname    !! names of the longitude and latitude
    REAL(r_std), INTENT(inout)                 :: varmin, varmax          !! min/max values to use for the
                                                                          !!   renormalization
! Transformations to apply to the original read values from the file
    LOGICAL, INTENT(in)                        :: noneg                   !! no negative values
    CHARACTER(LEN=50), INTENT(in)              :: masktype                !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below maskvalues(1)
                                                                          !!   'mabove': take values above maskvalues(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
                                                                          !!      maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
                                                                          !!        (normalized by maskvalues(3))
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)      :: maskvalues              !! values to use to mask (according to 
                                                                          !!   `masktype') 
    CHARACTER(LEN=250), INTENT(in)             :: maskvarname             !! name of the variable to use to mask 
    INTEGER(i_std), INTENT(in)                 :: dim1, dim2              !! 4D complementary dimension sizes
    CHARACTER(LEN=50), INTENT(in)              :: typefrac                !! Type of fraction retrieval:
                                                                          !!   'default': standard

    INTEGER(i_std), INTENT(in)                 :: initime                 !! Initial time to take (for input data 
                                                                          !!   with time values, -1 not used)
    REAL(r_std), INTENT(in)                    :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
                                                                          !!   (-un, not used)
    !! 0.2 Modified variables
    !
    !! 0.3 output variables
    !
    REAL(r_std), DIMENSION(nbpt), INTENT(out) :: outvar1D                 !! 1D output variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: aoutvar                 !! availability of input data to
                                                                          !!   interpolate output variable 
                                                                          !!   (on the nbpt space)
    !
    !! 0.4 Local variables
    !
    INTEGER(i_std)                             :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
    INTEGER                                    :: inNdims                 !! Number of dimensions of the input data
    INTEGER                                    :: inNLldims               !! Number of dimensions of the lon/lat input data
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: lat_in, lon_in          !! latitude and
                                                                          !! longitude, extract from input map
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: sub_area                !! the area of the fine grid in the 
                                                                          !!   model grid ???
                                                                          !! cf src_global/interpol_help.f90, 
                                                                          !!   line 377, called "areaoverlap"
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)  :: sub_index             !! the indexes from the grid boxes from the 
                                                                          !!   data that go into the 
                                                                          !!   model's boxes  
                                                                          !! cf src_global/interpol_help.f90,
                                                                          !!   line 300, called "ip"

    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: invar1D               !! 1D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: invar2D               !! 2D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: invar3D               !! 3D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D               !! 4D input values

    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: resol_in
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:)    :: mask
    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: maskvar
    INTEGER(i_std)                             :: nix, njx
    INTEGER(i_std)                             :: idi, nbvmax             !! nbvmax : number of maximum vegetation map
                                                                          !! points in the GCM grid ; idi : its counter
    CHARACTER(LEN=30)                          :: callsign                !! Allows to specify which variable is beeing
                                                                          !! treated
    LOGICAL                                    :: ok_interpol             !! optionnal return of aggregate_2d
    INTEGER                                    :: ALLOC_ERR 
    CHARACTER(LEN=50)                          :: S1
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: variableusetypes
    CHARACTER(LEN=256)                         :: msg
    INTEGER, DIMENSION(:), ALLOCATABLE         :: indims

! Debug vars
    INTEGER                                    :: nb_coord, nb_var, nb_gat


!_ ================================================================================================================================
! Allocating input data

    IF (printlev >= 3) WRITE(numout,*)"In interpweight_1D filename: ",TRIM(filename),          &
      " varname:'" // TRIM(varname) // "'"
    
    IF (is_root_prc) THEN
       IF (printlev >=5) THEN
          WRITE(numout,*) "Entering 'interpweight_1D'. Debug mode."
          WRITE(numout,*)'(/," --> fliodmpf")'
          CALL fliodmpf (TRIM(filename))
          WRITE(numout,*)'(/," --> flioopfd")'
       ENDIF
       CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
       IF (printlev >=5) THEN
          WRITE (numout,'(" Number of coordinate        in the file : ",I2)') nb_coord
          WRITE (numout,'(" Number of variables         in the file : ",I2)') nb_var
          WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
       ENDIF
    ENDIF

! Getting shape of input variable from input file
    inNdims = interpweight_get_varNdims_file(filename, varname)
    IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)

    InputDims: SELECT CASE (inNdims)
      CASE (1)
        CALL bcast(iml)

        ALLOCATE(invar1D(iml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, 0, 0, 0, 1, 1, invar1D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_1D input data dimensions _______'
          WRITE(numout,*)'    iml:',iml
        END IF

      CASE DEFAULT
        WRITE (S1,'(I3)') inNdims
        CALL ipslerr_p(3,'interpweight_1D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
          '','')

    END SELECT InputDims
    CALL bcast(invar1D)

    IF (printlev >= 5) WRITE(numout,*)'  interpweight_1D getting variable ended!'

! Getting longitudes/latitudes from input file
    inNLldims = interpweight_get_varNdims_file(filename, inlonname)
    
    IF (inNLldims == 1) THEN
      ! Allocate memory for latitudes
      ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable lon_in','','')

      ALLOCATE(lat_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable lat_in','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
        CALL flinget(fid, TRIM(inlatname), iml, 0, 0, 0, 1, 1, lat_in)
      ENDIF

    END IF
    CALL bcast(lon_in)
    CALL bcast(lat_in)

    IF (noneg) THEN
      IF (is_root_prc) THEN
        SELECT CASE (inNdims)
          CASE(1)
            CALL interpweight_modifying_input1D(iml, 0, 0, 0, 1, 1, zero, invar1D)
        END SELECT       
      END IF
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_1D modifying input ended!'

    IF (is_root_prc) CALL flinclo(fid)

    !
    ! Consider all points a priori
    !
    ALLOCATE(mask(iml), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable mask','','')
    mask = zero

    IF (TRIM(masktype) == 'var') THEN
      ALLOCATE(maskvar(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable maskvar','','')
      maskvar = zero
! Reads the mask from the file
      IF (is_root_prc) THEN
! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
        CALL flinget(fid, TRIM(maskvarname), iml, 0, 0, 0, 1, 1, maskvar)
        WHERE (maskvar > zero)
          mask = un 
        END WHERE
      END IF
      CALL bcast(mask)
    ELSE
      SELECT CASE (inNdims)
        CASE(1)
          CALL interpweight_masking_input1D(iml, 0, 0, 0, 0, 0, masktype, maskvalues, invar1D,        &
            mask)
      END SELECT
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_1D masking input ended!'

    ! Computation of nbvmax
    ! In 'condveg_soilalb, slowproc_update' nbvmax=200
    !
    ! The number of maximum vegetation map points in the GCM grid is estimated.
    ! Some lmargin is taken.
    !
    IF (is_root_prc) THEN
      IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
!       Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum  &
!         lon,lat resolution given as input paramater
         nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
         njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
      ELSE
!       Assuming regular lon/lat projection for the input data
        CALL interpweight_calc_resolution_in(lon_in, lat_in, iml, jml, mincos, R_Earth, pi, resol_in)
        nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
        njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
      END IF
      nbvmax = MAX(nix*njx, 200)
    END IF  
    CALL bcast(nbvmax)
    CALL bcast(nix)
    CALL bcast(njx)

    !
    callsign = TRIM(varname) // ' map'
    !
    ok_interpol = .FALSE.
    DO WHILE ( .NOT. ok_interpol )
       IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
       IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx


       ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable sub_index','','')

       sub_index(:,:)=zero

       ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable sub_area','','')

       sub_area(:,:)=zero
       IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon_in), MAXVAL(lon_in)
       IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat_in), MAXVAL(lat_in)

       !
       CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac,  &
            &                iml, lon_in, lat_in, maxresollon, maxresollat, callsign, &
            &                nbvmax, sub_index, sub_area, ok_interpol)
       
       !
       IF ( .NOT. ok_interpol ) THEN
          DEALLOCATE(sub_area)
          DEALLOCATE(sub_index)
          nbvmax = nbvmax * 2
       ENDIF
    ENDDO
    IF (printlev>=4) WRITE (numout,*) '  interpweight_1D: read/allocate OK'

! Getting variables thresholds
    ALLOCATE(variableusetypes(Nvariabletypes), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variableusetypes','','')

    IF (variabletypes(1) == -un) THEN
      variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
      IF (printlev>=4) WRITE(numout,*) 'interpweight_1D: Note: Using default ',&
           'equivalence between dimension in the file and dimension in the variable, for file ', filename
      varmin = 1._r_std
      varmax = Nvariabletypes*1._r_std
    ELSE
      variableusetypes = variabletypes    
    END IF

    outvar1D = zero

! Providing the fractions of each type for each interpolated grid point
    CALL interpweight_provide_fractions1D(typefrac, nbpt, Nvariabletypes, variableusetypes,           &
      iml, jml, lml, tml, nbvmax, zero, invar1D, sub_area, sub_index, varmax, varmin,                 &
      deux, lalo, outvar1D, aoutvar)
    IF (printlev>=5) WRITE (numout,*)'  interpweight_1D end of interp_provide_fractions1D'

    IF (printlev>=3) THEN
       IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
          WRITE(numout,*) '  interpweight_1D total number of points: ', nbpt
          WRITE(numout,*) '  interpweight_1D interpolated: ', COUNT(aoutvar /= -1),                 &
               ' non-interpolated: ', COUNT(aoutvar == -1)
          msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) //       & 
               "' have been succesfully interpolated !!"
          WRITE(numout,*) '  interpweight_1D: ' // TRIM(msg)
       ELSE
          msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
          WRITE(numout,*) "  interpweight_1D: '" // TRIM(msg)
          WRITE(numout,*) '  interpweight_1D total number of points: ', nbpt
          WRITE(numout,*) '  interpweight_1D interpolated: ', COUNT(aoutvar /= -1),               &
               ' non-interpolated: ', COUNT(aoutvar == -1)
       END IF
    END IF

! Scaling to 1 the quality variable
    DO ip=1, nbpt
       IF (contfrac(ip) == zero) THEN
          ! HACK: This should not happen but it does sometimes.
          ! We can't calculate aoutvar, therefor assign a value -2
          aoutvar(ip) = -2
       ELSE
          IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
       END IF
    END DO

    IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
    IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
    IF (ALLOCATED(invar1D)) DEALLOCATE(invar1D)
    DEALLOCATE(mask)
    IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
    DEALLOCATE(sub_area)
    DEALLOCATE(sub_index)

    RETURN

  END SUBROUTINE interpweight_1D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_2D
!!
!>\BRIEF        ! Interpolate any input file to the grid of the model for a 2D-variable  
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): outvar2D, aoutvar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_2D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours,       &
    contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype,               &
    maskvalues, maskvarname, dim1, dim2, initime, typefrac,                                           & 
    maxresollon, maxresollat, outvar2D, aoutvar)

    USE ioipsl
    USE constantes_var
    USE mod_orchidee_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !! 0.1 Input variables
    !
    INTEGER(i_std), INTENT(in)                 :: nbpt                    !! Number of points for which the data 
                                                                          !!   needs to be interpolated
    INTEGER(i_std), INTENT(in)                 :: Nvariabletypes          !! Number of types of the variable
    REAL(r_std), DIMENSION(Nvariabletypes),    &                          !! Vector of values of the types
      INTENT(in)    :: variabletypes                                      !!   (-1, not used)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo                    !! Vector of latitude and longitudes 
                                                                          !! (beware of the order = 1 : latitude, 
                                                                          !!   2 : longitude)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution              !! The size in km of each grid-box in X and Y
    !
    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours    !! Vector of neighbours for each grid point 
                                                                          !!   1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) 
                            
    REAL(r_std), DIMENSION(nbpt),  INTENT(in)  :: contfrac                !! Fraction of land in each grid box.
    CHARACTER(LEN=80), INTENT(in)              :: filename                !! name of the input map to read

    CHARACTER(LEN=80), INTENT(in)              :: varname                 !! name of the variable to interpolate
    CHARACTER(LEN=80), INTENT(in)              :: inlonname, inlatname    !! names of the longitude and latitude
    REAL(r_std), INTENT(inout)                 :: varmin, varmax          !! min/max values to use for the
                                                                          !!   renormalization
! Transformations to apply to the original read values from the file
    LOGICAL, INTENT(in)                        :: noneg                   !! no negative values
    CHARACTER(LEN=50), INTENT(in)              :: masktype                !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below maskvalues(1)
                                                                          !!   'mabove': take values above maskvalues(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
                                                                          !!      maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
                                                                          !!        (normalized by maskvalues(3))
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)      :: maskvalues              !! values to use to mask (according to 
                                                                          !!   `masktype') 
    CHARACTER(LEN=250), INTENT(in)             :: maskvarname             !! name of the variable to use to mask 
    CHARACTER(LEN=50), INTENT(in)              :: typefrac                !! Type of fraction retrieval:
                                                                          !!   'default': standard
    INTEGER(i_std), INTENT(in)                 :: dim1, dim2              !! 4D complementary dimension sizes
    INTEGER(i_std), INTENT(in)                 :: initime                 !! Initial time to take (for input data 
                                                                          !!   with time values, -1 not used)
    REAL(r_std), INTENT(in)                    :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
    !! 0.2 Modified variables
    !
    !  0.3 OUTPUT
    !
!    REAL(r_std), INTENT(out)                   ::  laimap(nbpt,nvm,12)    !! lai read variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt,Nvariabletypes), INTENT(out) :: outvar2D  !! 2D output variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: aoutvar                 !! availability of input data to
                                                                          !!   interpolate output variable 
                                                                          !!   (on the nbpt space)
    !
    !  0.4 LOCAL
    !
    INTEGER(i_std)                             :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
    INTEGER                                    :: inNdims                 !! Number of dimensions of the input data
    INTEGER                                    :: inNLldims               !! Number of dimensions of the lon/lat input data
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: lat_in, lon_in          !! latitude and
                                                                          !! longitude, extract from input map
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: lat, lon                !! en 2D ???
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: sub_area                !! the area of the fine grid in the 
                                                                          !!   model grid ???
                                                                          !! cf src_global/interpol_help.f90, 
                                                                          !!   line 377, called "areaoverlap"
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index           !! the indexes from the grid boxes from the 
                                                                          !!   data that go into the 
                                                                          !!   model's boxes  
                                                                          !! cf src_global/interpol_help.f90,
                                                                          !!   line 300, called "ip"

!    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: invar1D               !! 1D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: invar2D               !! 2D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: invar3D               !! 3D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D               !! 4D input values

    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: resol_in
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)  :: mask
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: maskvar
    INTEGER(i_std)                             :: nix, njx
    INTEGER(i_std)                             :: idi, nbvmax             !! nbvmax : number of maximum vegetation map
                                                                          !! points in the GCM grid ; idi : its counter
    CHARACTER(LEN=30)                          :: callsign                !! Allows to specify which variable is beeing
                                                                          !! treated
    LOGICAL                                    :: ok_interpol             !! optionnal return of aggregate_2d
    !
    INTEGER                                    :: ALLOC_ERR 
    CHARACTER(LEN=50)                          :: S1
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: variableusetypes
    CHARACTER(LEN=256)                         :: msg
    INTEGER, DIMENSION(:), ALLOCATABLE         :: indims

! Debug vars
    INTEGER                                    :: nb_coord, nb_var, nb_gat


!_ ================================================================================================================================
! Allocating input data

    IF (printlev >= 3) WRITE(numout,*)"In interpweight_2D filename: ",TRIM(filename),          &
      " varname:'" // TRIM(varname) // "'"
    !
    IF (is_root_prc) THEN
       IF (printlev >=5) THEN
          WRITE(numout,*) "Entering 'interpweight_2D'. Debug mode."
          WRITE(numout,*)'(/," --> fliodmpf")'
          CALL fliodmpf (TRIM(filename))
          WRITE(numout,*)'(/," --> flioopfd")'
       ENDIF
       CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
       IF (printlev >=5) THEN
          WRITE (numout,'(" Number of coordinate        in the file : ",I2)') nb_coord
          WRITE (numout,'(" Number of variables         in the file : ",I2)') nb_var
          WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
       ENDIF
    ENDIF

! Getting shape of input variable from input file
    inNdims = interpweight_get_varNdims_file(filename, varname)
    IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
    
    ! Read variable (varname) from input file
    InputDims: SELECT CASE (inNdims)
      CASE (2)
        CALL bcast(iml)
        CALL bcast(jml)

        ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_2D input data dimensions: '
          WRITE(numout,*)'    iml, jml:',iml, jml
        END IF

      CASE (3)
        ALLOCATE(indims(3), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of indims','','')
        indims = interpweight_get_var3dims_file(filename, varname)

        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)

        ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_2D input data dimensions:'
          WRITE(numout,*)'    iml, jml, lml:',iml, jml, lml
        END IF

      CASE (4)
        ALLOCATE(indims(4), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of indims','','')
        indims = interpweight_get_var4dims_file(filename, varname)
        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)
!       Taken from `slowproc_update'
        IF (initime /= -1) THEN
          inNdims = 3
          IF (printlev >= 3)  WRITE(numout,*) &
               'interpweight_2D: taking into account the initial reading time:',initime
          IF (initime <= 0) tml = 0

          CALL bcast(tml)
          ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF ( initime > 0 ) THEN
            IF (is_root_prc) THEN
              IF (initime <= tml) THEN
                 CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
              ELSE
                 CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
              ENDIF
            ENDIF
          ELSE
            IF (is_root_prc) THEN
              CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
            ENDIF
          ENDIF
        ELSE
          CALL bcast(tml)
          ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
        END IF

        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_2D input data dimensions'
          WRITE(numout,*)'    iml, jml, lml, tml :',iml, jml, lml, tml
        END IF

      CASE DEFAULT
        WRITE (S1,'(I3)') inNdims
        CALL ipslerr_p(3,'interpweight_2D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
          '','')

    END SELECT InputDims
    CALL bcast(invar2D)

    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2D getting variable ended!'

! Getting longitudes/latitudes from input file
    inNLldims = interpweight_get_varNdims_file(filename, inlonname)
    
    IF (inNLldims == 1) THEN
      ! Allocate memory for latitudes
      ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lon_in','','')

      ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lat_in','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
        CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
      ENDIF
      CALL bcast(lon_in)
      CALL bcast(lat_in)

      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lon','','')
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lat','','')

      DO ip=1,iml
         lat(ip,:) = lat_in(:)
      ENDDO
      DO jp=1,jml
         lon(:,jp) = lon_in(:)
      ENDDO
    ELSE
      ! Allocate memory for longitude
      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Pb in allocation for lon','','')
      ! Allocate memory for latitudes
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Pb in allocation for lat','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
        CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
      ENDIF
    END IF
    CALL bcast(lon)
    CALL bcast(lat)

    IF (is_root_prc) THEN
      ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr(3,'interpweight_2D','Problem in allocation of variable resol_in','','')
    END IF

    IF (noneg) THEN
      IF (is_root_prc) THEN
        SELECT CASE (inNdims)
          CASE(2)
            CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
          CASE(3)
            CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
          CASE(4)
            CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
        END SELECT       
      END IF
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2D modifying input ended!'

    !
    ! Consider all points a priori
    !
    ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable mask','','')
    mask = zero

    IF (TRIM(masktype) == 'var') THEN
      ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable maskvar','','')
      maskvar(:,:) = 0
! Reads the mask from the file
      IF (is_root_prc) THEN
! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
        CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
        WHERE (maskvar > zero)
          mask = un 
        END WHERE
      END IF
      CALL bcast(mask)
    ELSE
      ! Create land-sea mask based on values from the variable according to maskvalues and masktype
      SELECT CASE (inNdims)
        CASE(2)
          CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D,      &
            mask)
        CASE(3)
          CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D,    &
            mask)
        CASE(4)
          CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D,  &
            mask)
      END SELECT
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2D masking input ended!'

    IF (is_root_prc) CALL flinclo(fid)

! Assuming regular lon/lat projection for the input data
    IF (is_root_prc) THEN
      CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
    END IF

    ! Computation of nbvmax
    ! In 'condveg_soilalb, slowproc_update' nbvmax=200
    !
    ! The number of maximum vegetation map points in the GCM grid is estimated.
    ! Some lmargin is taken.
    !
    IF (is_root_prc) THEN
      IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
!       Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum  &
!         lon,lat resolution given as input paramater
         nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
         njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
      ELSE
!       Assuming regular lon/lat projection for the input data
        CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
        nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
        njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
      END IF
      nbvmax = MAX(nix*njx, 200)
    END IF  
    CALL bcast(nbvmax)
    CALL bcast(nix)
    CALL bcast(njx)

    !
    callsign = TRIM(varname) // ' map'
    !
    ok_interpol = .FALSE.
    DO WHILE ( .NOT. ok_interpol )
       IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
       IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx


       ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable sub_index','','')

       sub_index(:,:,:)=zero

       ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable sub_area','','')

       sub_area(:,:)=zero
       IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon), MAXVAL(lon)
       IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat), MAXVAL(lat)

       !
       CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
            &                iml, jml, lon, lat, mask, callsign, &
            &                nbvmax, sub_index, sub_area, ok_interpol)
       
       !
       IF ( .NOT. ok_interpol ) THEN
          DEALLOCATE(sub_area)
          DEALLOCATE(sub_index)
          nbvmax = nbvmax * 2
       ENDIF
    ENDDO
    IF (printlev>=4) WRITE (numout,*) '  interpweight_2D: read/allocate OK'

! Getting variables thresholds
    ALLOCATE(variableusetypes(Nvariabletypes), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variableusetypes','','')
    

    IF (variabletypes(1) == -un) THEN
      variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
      IF (printlev>=4) THEN
         WRITE(numout,*)'  interpweight_2D: Note: Using default equivalence between dimension ',&
              'in the file and dimension in the variable, for file ', filename
      END IF
      varmin = 1._r_std
      varmax = Nvariabletypes*1._r_std
    ELSE
      variableusetypes = variabletypes    
    END IF

    outvar2D = zero

! Providing the fractions of each type for each interpolated grid point
    CALL interpweight_provide_fractions2D(typefrac, nbpt, Nvariabletypes, variableusetypes,           &
      iml, jml, lml, tml, nbvmax, zero, invar2D, sub_area, sub_index, varmax, varmin,                 &
      deux, lalo, outvar2D, aoutvar)
    IF (printlev>=5) WRITE (numout,*)'  interpweight_2D end of provide_fractions2D'

    IF (printlev>=3) THEN
       IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
          WRITE(numout,*) '  interpweight_2D total number of points: ', nbpt
          WRITE(numout,*) '  interpweight_2D interpolated: ', COUNT(aoutvar /= -1),                   &
               ' non-interpolated: ', COUNT(aoutvar == -1)
          msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) //       & 
               "' have been succesfully interpolated !!"
          WRITE(numout,*) '  interpweight_2D: ' // TRIM(msg)
       ELSE
          msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
          WRITE(numout,*) "interpweight_2D: '" // TRIM(msg)
          WRITE(numout,*) '  Total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
       END IF
    END IF

! Scaling to 1 the quality variable
    DO ip=1, nbpt
       IF (contfrac(ip) == zero) THEN
          ! HACK: This should not happen but it does sometimes.
          ! We can't calculate aoutvar, therefor assign a value -2
          aoutvar(ip) = -2
       ELSE
          IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
       END IF
    END DO

    IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
    IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
    DEALLOCATE(lon)
    DEALLOCATE(lat)
    IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
    IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
    IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
    DEALLOCATE(mask)
    IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
    DEALLOCATE(sub_area)
    DEALLOCATE(sub_index)
    IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)

    RETURN

  END SUBROUTINE interpweight_2D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_3D
!!
!>\BRIEF        ! Interpolate any input file to the grid of the model for a 3D-variable  
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): outvar3D, aoutvar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_3D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours,       &
    contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype,               &
    maskvalues, maskvarname, dim1, dim2, initime, typefrac,                                           &
    maxresollon, maxresollat, outvar3D, aoutvar)

    USE ioipsl
    USE constantes_var
    USE mod_orchidee_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !  0.1 INPUT
    !
    INTEGER(i_std), INTENT(in)                 :: nbpt                    !! Number of points for which the data 
                                                                          !!   needs to be interpolated
    INTEGER(i_std), INTENT(in)                 :: Nvariabletypes          !! Number of types of the variable
    REAL(r_std), DIMENSION(Nvariabletypes),    &                          !! Vector of values of the types
      INTENT(in)    :: variabletypes                                      !!   (-1, not used)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo                    !! Vector of latitude and longitudes 
                                                                          !! (beware of the order = 1 : latitude, 
                                                                          !!   2 : longitude)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution              !! The size in km of each grid-box in X and Y
    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours    !! Vector of neighbours for each grid point 
                                                                          !!   1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) 
                            
    REAL(r_std), DIMENSION(nbpt), INTENT(in)   :: contfrac                !! Fraction of land in each grid box.
    CHARACTER(LEN=80), INTENT(in)              :: filename                !! name of the input map to read

    CHARACTER(LEN=80), INTENT(in)              :: varname                 !! name of the variable to interpolate
    CHARACTER(LEN=80), INTENT(in)              :: inlonname, inlatname    !! names of the longitude and latitude
                                                                          !!   in the input file
! Transformations to apply to the original read values from the file
    LOGICAL, INTENT(in)                        :: noneg                   !! no negative values
    CHARACTER(LEN=50), INTENT(in)              :: masktype                !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below maskvalues(1)
                                                                          !!   'mabove': take values above maskvalues(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
                                                                          !!      maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
                                                                          !!        (normalized by maskvalues(3))
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)      :: maskvalues              !! values to use to mask (according to 
                                                                          !!   `masktype') 
    CHARACTER(LEN=250), INTENT(in)             :: maskvarname             !! name of the variable to use to mask 
    INTEGER, INTENT(in)                        :: dim1, dim2              !! 3/4D dimensions of the output variable
    REAL(r_std), DIMENSION(dim1), INTENT(inout) :: varmin, varmax          !! min/max values to use for the renormalization
    CHARACTER(LEN=50), INTENT(in)              :: typefrac                !! Type of fraction retrieval:
                                                                          !!   'default': standard
    INTEGER(i_std), INTENT(in)                 :: initime                 !! Initial time to take (for input data 
                                                                          !!   with time values, -1 not used)
    REAL(r_std), INTENT(in)                    :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
                                                                          !!   (-un, not used)
    !! 0.2 Modified variables
    !
    !  0.3 OUTPUT
    !
    REAL(r_std), DIMENSION(nbpt,dim1), INTENT(out) :: outvar3D            !! 3D output variable and re-dimensioned
! aovar: availability of input data to interpolate output variable (on the nbpt space)
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: aoutvar                 !! availability of input data to
                                                                          !!   interpolate output variable 
                                                                          !!   (on the nbpt space)
    !
    !  0.4 Local variables
    !
    INTEGER(i_std)                             :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
    INTEGER                                    :: inNdims                 !! Number of dimensions of the input data
    INTEGER                                    :: inNLldims               !! Number of dimensions of the lon/lat input data
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: lat_in, lon_in          !! latitude and
                                                                          !! longitude, extract from input map
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: lat, lon                !! en 2D ???
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: sub_area                !! the area of the fine grid in the 
                                                                          !!   model grid ???
                                                                          !! cf src_global/interpol_help.f90, 
                                                                          !!   line 377, called "areaoverlap"
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index           !! the indexes from the grid boxes from the 
                                                                          !!   data that go into the 
                                                                          !!   model's boxes  
                                                                          !! cf src_global/interpol_help.f90,
                                                                          !!   line 300, called "ip"

!    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: invar1D               !! 1D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: invar2D               !! 2D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: invar3D               !! 3D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D               !! 4D input values

    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: resol_in
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)  :: mask
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: maskvar
    INTEGER(i_std)                             :: nix, njx
    INTEGER(i_std)                             :: idi, nbvmax             !! nbvmax : number of maximum vegetation map
                                                                          !! points in the GCM grid ; idi : its counter
    CHARACTER(LEN=30)                          :: callsign                !! Allows to specify which variable is beeing
                                                                          !! treated
    LOGICAL                                    :: ok_interpol             !! optionnal return of aggregate_2d
    INTEGER                                    :: ALLOC_ERR 

    CHARACTER(LEN=50)                          :: S1
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: variableusetypes
    CHARACTER(LEN=256)                         :: msg

! Debug vars
    INTEGER                                    :: nb_coord, nb_var, nb_gat
    CHARACTER(LEN=5)                           :: iS
    CHARACTER(LEN=50)                          :: suminS

    INTEGER, ALLOCATABLE, DIMENSION(:)         :: indims


!_ ================================================================================================================================
    IF (printlev >= 3) WRITE(numout, *)"In interpweight_3D filename: ",TRIM(filename),         &
      " varname:'" // TRIM(varname) // "'"
! Allocating input data
    IF (is_root_prc) THEN
       IF (printlev >=5) THEN
          WRITE(numout,*) "Entering 'interpweight_3D'. Debug mode."
          WRITE(numout,*)'(/," --> fliodmpf")'
          CALL fliodmpf (TRIM(filename))
          WRITE(numout,*)'(/," --> flioopfd")'
       ENDIF
       CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
       IF (printlev >=5) THEN
          WRITE (numout,'(" Number of coordinate        in the file : ",I2)') nb_coord
          WRITE (numout,'(" Number of variables         in the file : ",I2)') nb_var
          WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
       ENDIF
    ENDIF

! Getting shape of input variable from input file
    inNdims = interpweight_get_varNdims_file(filename, varname)
    IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)

    InputDims: SELECT CASE (inNdims)
      CASE (2)
        CALL bcast(iml)
        CALL bcast(jml)

        ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_3D input data dimensions'
          WRITE(numout,*)'    iml, jml:',iml, jml
        END IF

      CASE (3)
        ALLOCATE(indims(3), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of indims','','')
        indims = interpweight_get_var3dims_file(filename, varname)

        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)

        ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','') 
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 2, 1, invar3D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_3D input data dimensions'
          WRITE(numout,*)'    iml, jml, lml :',iml, jml, lml
        END IF

      CASE (4)
        ALLOCATE(indims(4), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of indims','','')
        
        indims = interpweight_get_var4dims_file(filename, varname)
        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)
!       Taken from `slowproc_update'
        IF (initime /= -1) THEN
          inNdims = 3
          IF (printlev >= 3) WRITE (numout,*) &
               'interpweight_3D: taking into account the initial reading time:',initime
          IF (initime <= 0) tml = 0

          CALL bcast(tml)
          ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable invar3D','','')
          IF ( initime > 0 ) THEN
            IF (is_root_prc) THEN
              IF (initime <= tml) THEN
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
              ELSE
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
              ENDIF
            ENDIF
          ELSE
            IF (is_root_prc) THEN
              CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
            ENDIF
          ENDIF
        ELSE
          CALL bcast(tml)
          ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable invar4D','','')
          IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
        END IF

        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_3D input data dimensions'
          WRITE(numout,*)'    iml, jml, lml, tml:',iml, jml, lml, tml
        END IF

      CASE DEFAULT
        WRITE (S1,'(I3)') inNdims
        CALL ipslerr_p(3,'interpweight_3D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
          '','')

    END SELECT InputDims
    CALL bcast(invar3D)
  
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_3D getting variable ended!'
! Test lecture
    DO ip=1,lml
      IF ( SUM(invar3D(:,:,ip)) < 0.00001 ) THEN
        WRITE(iS, '(i5)')ip 
        WRITE(suminS, '(f20.10)')SUM(invar3d(:,:,ip)) 
        msg = 'All values are zero or smaller than 0.00001 for data for type ' // TRIM(iS)
        CALL ipslerr(1,'interpweight_3D',TRIM(msg),'','')
      END IF
    END DO

! Getting longitudes/latitudes from input file
    inNLldims = interpweight_get_varNdims_file(filename, inlonname)
    IF (inNLldims == 1) THEN
      ! Allocate memory for latitudes
      ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lon_in','','')

      ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lat_in','','')

      IF (is_root_prc) THEN
        CALL flininspect(fid)
        CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
        CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
      ENDIF
      CALL bcast(lon_in)
      CALL bcast(lat_in)

      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lon','','')
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lat','','')

      DO ip=1,iml
         lat(ip,:) = lat_in(:)
      ENDDO
      DO jp=1,jml
         lon(:,jp) = lon_in(:)
      ENDDO
    ELSE
      ! Allocate memory for longitude
      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Pb in allocation for lon','','')
      ! Allocate memory for latitudes
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Pb in allocation for lat','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
        CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
      ENDIF
    END IF
    CALL bcast(lon)
    CALL bcast(lat)

    ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable resol_in','','')

    IF (noneg) THEN
      IF (is_root_prc) THEN
        SELECT CASE (inNdims)
          CASE(2)
            CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
          CASE(3)
            CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
          CASE(4)
            CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
        END SELECT       
      END IF
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_3D modifying input ended!'

    ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable mask','','')
    mask(:,:) = 0

    IF (TRIM(masktype) == 'var') THEN
      ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable maskvar','','')
      maskvar(:,:) = 0.
! Reads the mask from the file
      IF (is_root_prc) THEN
! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
        CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
        WHERE (maskvar > zero)
          mask = un 
        END WHERE
      END IF
      CALL bcast(mask)
    ELSE
      SELECT CASE (inNdims)
        CASE(2)
          CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D,      &
            mask)
        CASE(3)
          CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D,    &
            mask)
        CASE(4)
          CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D,  &
            mask)
      END SELECT
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_3D masking input ended!'

    IF (is_root_prc) CALL flinclo(fid)

    ! Computation of nbvmax
    ! In 'condveg_soilalb, slowproc_update' nbvmax=200
    !
    ! The number of maximum vegetation map points in the GCM grid is estimated.
    ! Some lmargin is taken.
    !
    IF (is_root_prc) THEN
      IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
!       Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum  &
!         lon,lat resolution given as input paramater
         nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
         njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
      ELSE
!       Assuming regular lon/lat projection for the input data
        CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
        nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
        njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
      END IF
      nbvmax = MAX(nix*njx, 200)
    END IF  
    CALL bcast(nbvmax)
    CALL bcast(nix)
    CALL bcast(njx)

    !
    callsign = TRIM(varname) // ' map'
    !
    ok_interpol = .FALSE.
    DO WHILE ( .NOT. ok_interpol )
       IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
       IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx

       ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable sub_index','','')

       sub_index(:,:,:)=zero

       ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable sub_area','','')

       sub_area(:,:)=zero

       IF (printlev >= 3) THEN
         WRITE(numout,*) 'interpweight_3D:'
         WRITE(numout,*) '  Carteveg range LON:', MINVAL(lon), MAXVAL(lon)
         WRITE(numout,*) '  Carteveg range LAT:', MINVAL(lat), MAXVAL(lat)
       END IF
       
       CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
            &                iml, jml, lon, lat, mask, callsign, &
            &                nbvmax, sub_index, sub_area, ok_interpol)
       
       IF ( .NOT. ok_interpol ) THEN
          DEALLOCATE(sub_area)
          DEALLOCATE(sub_index)
          nbvmax = nbvmax * 2
       ENDIF
    ENDDO

! Getting variables thresholds
    ALLOCATE(variableusetypes(Nvariabletypes),STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable variableusetypes','','')

    IF (variabletypes(1) == -un) THEN
      variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
      IF (printlev>=4) WRITE(numout,*)'interpweight_3D: Note: Using default equivalence between ',&
           'dimension in the file and dimension in the variable, for file ', filename
      varmin = 1._r_std
      varmax = Nvariabletypes*1._r_std
    ELSE
      variableusetypes = variabletypes
    END IF
    
    outvar3D = zero

! Providing the fractions of each type for each interpolated grid point
    CALL interpweight_provide_fractions3D(typefrac, nbpt, Nvariabletypes, variableusetypes,           &
      iml, jml, lml, tml, nbvmax, zero, invar3D, sub_area, sub_index, varmax, varmin,                 &
      deux, lalo, outvar3D, aoutvar)
    IF (printlev >= 5) WRITE(numout,*)'interpweight_3D end of interpweight_provide_fractions3D'

    IF (printlev>=2) THEN
       IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
          WRITE(numout,*) 'interpweight_3D total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
          msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) //       & 
               "' have been succesfully interpolated on the current process !!"
          WRITE(numout,*)'   ' // TRIM(msg)
       ELSE
          msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done"
          WRITE(numout,*) "interpweight_3D: '" // TRIM(msg)
          WRITE(numout,*) '  Total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
       END IF
    END IF

! Scaling to 1 the quality variable
    DO ip=1, nbpt
       IF (contfrac(ip) == zero) THEN
          ! HACK: This should not happen but it does sometimes.
          ! We can't calculate aoutvar, therefor assign a value -2
          aoutvar(ip) = -2
       ELSE
          IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
       END IF
    END DO

    IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
    IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
    DEALLOCATE(lon)
    DEALLOCATE(lat)
    IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
    IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
    IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
    DEALLOCATE(mask)
    IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
    DEALLOCATE(sub_area)
    DEALLOCATE(sub_index)
    IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)

    RETURN

  END SUBROUTINE interpweight_3D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_4D
!!
!>\BRIEF        ! Interpolate any input file to the grid of the model for a 4D-variable  
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): outvar4D, aoutvar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_4D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours,       &
    contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype,               &
    maskvalues, maskvarname, dim1, dim2, initime, typefrac,                                           &
    maxresollon, maxresollat, outvar4D, aoutvar)

    USE ioipsl
    USE constantes_var
    USE mod_orchidee_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !! 0.1 Input variables
    !
    INTEGER(i_std), INTENT(in)                 :: nbpt                    !! Number of points for which the data 
                                                                          !!   needs to be interpolated
    INTEGER(i_std), INTENT(in)                 :: Nvariabletypes          !! Number of types of the variable
    REAL(r_std), DIMENSION(Nvariabletypes),    &                          !! Vector of values of the types
      INTENT(in)    :: variabletypes                                      !!   (-1, not used)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo                    !! Vector of latitude and longitudes 
                                                                          !! (beware of the order = 1 : latitude, 
                                                                          !!   2 : longitude)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution              !! The size in km of each grid-box in X and Y
    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours    !! Vector of neighbours for each grid point 
                                                                          !!   1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) 
                            
    REAL(r_std), DIMENSION(nbpt), INTENT(in)   :: contfrac                !! Fraction of land in each grid box.
    CHARACTER(LEN=80), INTENT(in)              :: filename                !! name of the input map to read

    CHARACTER(LEN=80), INTENT(in)              :: varname                 !! name of the variable to interpolate
    CHARACTER(LEN=80), INTENT(in)              :: inlonname, inlatname    !! names of the longitude and latitude
! Transformations to apply to the original read values from the file
    LOGICAL, INTENT(in)                        :: noneg                   !! no negative values
    CHARACTER(LEN=50), INTENT(in)              :: masktype                !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below maskvalues(1)
                                                                          !!   'mabove': take values above maskvalues(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
                                                                          !!      maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
                                                                          !!        (normalized by maskedvalues(3))
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)      :: maskvalues              !! values to use to mask (according to 
                                                                          !!   `masktype') 
    CHARACTER(LEN=250), INTENT(in)             :: maskvarname             !! name of the variable to use to mask 
    INTEGER, INTENT(in)                        :: dim1, dim2              !! 3/4D dimensions of the output variable
                                                                          !! dim1 = 1 fpr outvar2D
                                                                          !!   in the input file
    REAL(r_std), DIMENSION(dim1), INTENT(inout) :: varmin, varmax         !! min/max values to use for the renormalization
    CHARACTER(LEN=50), INTENT(in)              :: typefrac                !! Type of fractions retrieval:
                                                                          !!   'default': standard

    INTEGER(i_std), INTENT(in)                 :: initime                 !! Initial time to take (for input data 
                                                                          !!   with time values, -1 not used)
    REAL(r_std), INTENT(in)                    :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
                                                                          !!   (-un, not used)
    !! 0.2 Modified variables
    !
    !! 0.3 Output variables
    !
!    REAL(r_std), INTENT(out)                   ::  laimap(nbpt,nvm,12)    !! lai read variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt,dim1,dim2), INTENT(out) :: outvar4D       !! 4D output variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: aoutvar                 !! availability of input data to
                                                                          !!   interpolate output variable 
                                                                          !!   (on the nbpt space)
    !
    !! 0.4 Local variables
    !
    INTEGER(i_std)                             :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
    INTEGER                                    :: inNdims                 !! Number of dimensions of the input data
    INTEGER                                    :: inNLldims               !! Number of dimensions of the lon/lat input data
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: lat_in, lon_in          !! latitude and
                                                                          !! longitude, extract from input map
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: lat, lon                !! en 2D ???
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: sub_area                !! the area of the fine grid in the 
                                                                          !!   model grid ???
                                                                          !! cf src_global/interpol_help.f90, 
                                                                          !!   line 377, called "areaoverlap"
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index           !! the indexes from the grid boxes from the 
                                                                          !!   data that go into the 
                                                                          !!   model's boxes  
                                                                          !! cf src_global/interpol_help.f90,
                                                                          !!   line 300, called "ip"

    INTEGER                                      :: ALLOC_ERR
!    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: invar1D               !! 1D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: invar2D               !! 2D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: invar3D               !! 3D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D               !! 4D input values

    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: resol_in
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)  :: mask
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: maskvar
    INTEGER(i_std)                             :: nix, njx
    INTEGER(i_std)                             :: idi, nbvmax             !! nbvmax : number of maximum vegetation map
                                                                          !! points in the GCM grid ; idi : its counter
    CHARACTER(LEN=30)                          :: callsign                !! Allows to specify which variable is beeing
                                                                          !! treated
    LOGICAL                                    :: ok_interpol             !! optionnal return of aggregate_2d
    CHARACTER(LEN=50)                          :: S1
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: variableusetypes
    CHARACTER(LEN=256)                         :: msg

! Debug vars
    INTEGER                                    :: nb_coord, nb_var, nb_gat

    INTEGER, ALLOCATABLE, DIMENSION(:)         :: indims


!_ ================================================================================================================================
    IF (printlev >= 3) WRITE(numout, *)"In interpweight_4D filename: ",TRIM(filename),         &
      " varname:'" // TRIM(varname) // "'"

! Allocating input data
    IF (is_root_prc) THEN
       IF (printlev >=5) THEN
          WRITE(numout,*) "Entering 'interpweight_4D'. Debug mode."
          WRITE(numout,*)'(/," --> fliodmpf")'
          CALL fliodmpf (TRIM(filename))
          WRITE(numout,*)'(/," --> flioopfd")'
       ENDIF
       CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
       IF (printlev >=5) THEN
          WRITE (numout,'(" Number of coordinate        in the file : ",I2)') nb_coord
          WRITE (numout,'(" Number of variables         in the file : ",I2)') nb_var
          WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
       ENDIF
    ENDIF

! Getting shape of input variable from input file
    inNdims = interpweight_get_varNdims_file(filename, varname)
    IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)

    InputDims: SELECT CASE (inNdims)
      CASE (2)
        CALL bcast(iml)
        CALL bcast(jml)

        ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_4D input data dimensions _______'
          WRITE(numout,*)'    iml, jml:',iml, jml
        END IF

      CASE (3)
        ALLOCATE(indims(3), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of indims','','')
        indims = interpweight_get_var3dims_file(filename, varname)
        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)

        ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_4D input data dimensions:'
          WRITE(numout,*)'    iml, jml, lml:',iml, jml, lml
        END IF

      CASE (4)
        ALLOCATE(indims(4), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of indims','','')
        indims = interpweight_get_var4dims_file(filename, varname)
        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)
!       Taken from `slowproc_update'
        IF (initime /= -1) THEN
          inNdims = 3
          IF (printlev >= 3) WRITE(numout,*) &
               'interpweight_4D: taking into account the initial reading time:',initime
          IF (initime <= 0) tml = 0

          CALL bcast(tml)
          ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF ( initime > 0 ) THEN
            IF (is_root_prc) THEN
              IF (initime <= tml) THEN
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
              ELSE
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
              ENDIF
            ENDIF
          ELSE
            IF (is_root_prc) THEN
              CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
            ENDIF
          ENDIF
        ELSE
          CALL bcast(tml)
          ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
        END IF

        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_4D input data dimensions'
          WRITE(numout,*)'    iml, jml, lml, tml:',iml, jml, lml, tml
        END IF

      CASE DEFAULT
        WRITE (S1,'(I3)') inNdims
        CALL ipslerr_p(3,'interpweight_4D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
          '','')

    END SELECT InputDims
    CALL bcast(invar4D)
  
! Getting longitudes/latitudes from input file
    inNLldims = interpweight_get_varNdims_file(filename, inlonname)
    IF (inNLldims == 1) THEN
      ! Allocate memory for latitudes
      ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lon_in','','')

      ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lat_in','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
        CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
      ENDIF
      CALL bcast(lon_in)
      CALL bcast(lat_in)

      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lon','','')
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lat','','')

      DO ip=1,iml
         lat(ip,:) = lat_in(:)
      ENDDO
      DO jp=1,jml
         lon(:,jp) = lon_in(:)
      ENDDO
    ELSE
      ! Allocate memory for longitude
      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Pb in allocation for lon','','')
      ! Allocate memory for latitudes
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Pb in allocation for lat','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
        CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
      ENDIF
    END IF
    CALL bcast(lon)
    CALL bcast(lat)

    IF (is_root_prc) THEN
      ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr(3,'interpweight_4D','Problem in allocation of variable resol_in','','')
    END IF

    IF (noneg) THEN
      IF (is_root_prc) THEN
        SELECT CASE (inNdims)
          CASE(2)
            CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
          CASE(3)
            CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
          CASE(4)
            CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
        END SELECT       
      END IF
    END IF

    IF (printlev >= 5) WRITE(numout,*)'  interpweight_4D modifying input ended!'

    ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable mask','','')

    IF (TRIM(masktype) == 'var') THEN
      ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable maskvar','','')
      maskvar(:,:) = 0
! Reads the mask from the file
      IF (is_root_prc) THEN
! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
        CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
        WHERE (maskvar > zero)
          mask = un 
        END WHERE
      END IF
      CALL bcast(mask)
    ELSE
      SELECT CASE (inNdims)
        CASE(2)
          CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D,      &
            mask)
        CASE(3)
          CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D,    &
            mask)
        CASE(4)
          CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D,  &
            mask)
      END SELECT
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_4D masking input ended!'

    IF (is_root_prc) CALL flinclo(fid)

    ! Computation of nbvmax
    ! In 'condveg_soilalb, slowproc_update' nbvmax=200
    !
    ! The number of maximum vegetation map points in the GCM grid is estimated.
    ! Some lmargin is taken.
    !
    IF (is_root_prc) THEN
      IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
!       Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum  &
!         lon,lat resolution given as input paramater
         nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
         njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
      ELSE
!       Assuming regular lon/lat projection for the input data
        CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
        nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
        njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
      END IF
      nbvmax = MAX(nix*njx, 200)
    END IF  
    CALL bcast(nbvmax)
    CALL bcast(nix)
    CALL bcast(njx)

    callsign = TRIM(varname) // ' map'

    ok_interpol = .FALSE.
    DO WHILE ( .NOT. ok_interpol )
       IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
       IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx

       ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable sub_index','','')

       sub_index(:,:,:)=zero

       ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable sub_area','','')

       sub_area(:,:)=zero
       CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
            &                iml, jml, lon, lat, mask, callsign, &
            &                nbvmax, sub_index, sub_area, ok_interpol)
       
       IF ( .NOT. ok_interpol ) THEN
          DEALLOCATE(sub_area)
          DEALLOCATE(sub_index)
          nbvmax = nbvmax * 2
       ENDIF
    ENDDO

! Getting variables thresholds
    ALLOCATE(variableusetypes(Nvariabletypes), STAT=ALLOC_ERR)
            IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable variableusetypesr','','')

    IF (variabletypes(1) == -un) THEN
      variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
      IF (printlev>=4) WRITE(numout,*)'  interpweight_4D: Note: Using default ',&
           'equivalence between dimension in the file and dimension in the variable, for file ', filename
      varmin = 1._r_std
      varmax = Nvariabletypes*1._r_std
    ELSE
      variableusetypes = variabletypes
    END IF
    
    outvar4D = zero

! Doing the interpolation as function of the type fo the values
    CALL interpweight_provide_fractions4D(typefrac, nbpt, Nvariabletypes, variableusetypes,           &
      iml, jml, lml, tml, nbvmax, zero, invar4D, sub_area, sub_index, varmax, varmin,                 &
      deux, lalo, outvar4D, aoutvar)
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_4D end of interpweight_provide_fractions4D'

    IF (printlev>=2) THEN
       IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
          WRITE(numout,*) 'interpweight_4D total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
          msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) //       & 
               "' have been succesfully interpolated !!"
          WRITE(numout,*)'  ' // TRIM(msg)
       ELSE
          msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
          WRITE(numout,*) "interpweight_4D: '" // TRIM(msg)
          WRITE(numout,*) '  Total number of points: ', nbpt
          WRITE(numout,*) '  Interpweight_4D interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
       END IF
    END IF

! Scaling to 1 the quality variable
    DO ip=1, nbpt
       IF (contfrac(ip) == zero) THEN
          ! HACK: This should not happen but it does sometimes.
          ! We can't calculate aoutvar, therefor assign a value -2
          aoutvar(ip) = -2
       ELSE
          IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
       END IF
    END DO

    IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
    IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
    DEALLOCATE(lon)
    DEALLOCATE(lat)
    IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
    IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
    IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
    DEALLOCATE(mask)
    IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
    DEALLOCATE(sub_area)
    DEALLOCATE(sub_index)
    IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)

    RETURN

  END SUBROUTINE interpweight_4D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_2Dcont
!!
!>\BRIEF        ! Interpolate any input file to the grid of the model for a continuos 2D-variable  
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): outvar2D, aoutvar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_2Dcont(nbpt, dim1, dim2, lalo, resolution, neighbours,                      &
    contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype,               &
    maskvalues, maskvarname, initime, typefrac, defaultvalue, defaultNOvalue,                         &
    outvar2D, aoutvar)

    USE ioipsl
    USE constantes_var
    USE mod_orchidee_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !! 0.1 Input variables
    !
    INTEGER(i_std), INTENT(in)                 :: nbpt                    !! Number of points for which the data 
                                                                          !!   needs to be interpolated
    INTEGER(i_std), INTENT(in)                 :: dim1, dim2              !! 4D complementary dimension sizes
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo                    !! Vector of latitude and longitudes 
                                                                          !! (beware of the order = 1 : latitude, 
                                                                          !!   2 : longitude)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution              !! The size in km of each grid-box in X and Y
    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours    !! Vector of neighbours for each grid point 
                                                                          !!   1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) 
                            
    REAL(r_std), DIMENSION(nbpt), INTENT(in)   :: contfrac                !! Fraction of land in each grid box.
    CHARACTER(LEN=80), INTENT(in)              :: filename                !! name of the input map to read

    CHARACTER(LEN=80), INTENT(in)              :: varname                 !! name of the variable to interpolate
    CHARACTER(LEN=80), INTENT(in)              :: inlonname, inlatname    !! names of the longitude and latitude
    REAL(r_std), INTENT(in)                    :: varmin, varmax          !! min/max values to use for the
                                                                          !!   renormalization
    REAL(r_std), INTENT(in)                    :: defaultvalue            !! default interpolated value 
    REAL(r_std), INTENT(in)                    :: defaultNOvalue          !! default interpolated NO-value 

! Transformations to apply to the original read values from the file
    LOGICAL, INTENT(in)                        :: noneg                   !! no negative values
    CHARACTER(LEN=50), INTENT(in)              :: masktype                !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below maskvalues(1)
                                                                          !!   'mabove': take values above maskvalues(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
                                                                          !!      maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
                                                                          !!        (normalized by maskedvalues(3))
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (1/0)
    REAL(r_std), DIMENSION(3), INTENT(in)      :: maskvalues              !! values to use to mask (according to 
                                                                          !!   `masktype') 
    CHARACTER(LEN=250), INTENT(in)             :: maskvarname             !! name of the variable to use to mask 
                                                                          !!   (according to `masktype') 
    CHARACTER(LEN=50), INTENT(in)              :: typefrac                !! Type of fraction retrieval:
                                                                          !!   'default': standard

    INTEGER(i_std), INTENT(in)                 :: initime                 !! Initial time to take (for input data 
                                                                          !!   with time values, -1: not used)
    !! 0.2 Modified variables
    !
    !! 0.3 Output variables
    !
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: outvar2D                !! 2D output variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: aoutvar                 !! availability of input data to
                                                                          !!   interpolate output variable 
                                                                          !!   (on the nbpt space)
    !
    !! 0.4 Local variables
    !
    INTEGER(i_std)                             :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
    INTEGER                                    :: inNdims                 !! Number of dimensions of the input data
    INTEGER                                    :: inNLldims               !! Number of dimensions of the lon/lat input data
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: lat_in, lon_in          !! latitude and
                                                                          !! longitude, extract from input map
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: lat, lon                !! en 2D ???
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: sub_area                !! the area of the fine grid in the 
                                                                          !!   model grid ???
                                                                          !! cf src_global/interpol_help.f90, 
                                                                          !!   line 377, called "areaoverlap"
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index           !! the indexes from the grid boxes from the 
                                                                          !!   data that go into the 
                                                                          !!   model's boxes  
                                                                          !! cf src_global/interpol_help.f90,
                                                                          !!   line 300, called "ip"

!    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: invar1D               !! 1D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: invar2D               !! 2D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: invar3D               !! 3D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D               !! 4D input values

    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: resol_in
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)  :: mask
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: maskvar
    INTEGER(i_std)                             :: nix, njx
    INTEGER(i_std)                             :: idi, nbvmax             !! nbvmax : number of maximum vegetation map
                                                                          !! points in the GCM grid ; idi : its counter
    CHARACTER(LEN=30)                          :: callsign                !! Allows to specify which variable is beeing
                                                                          !! treated
    LOGICAL                                    :: ok_interpol             !! optionnal return of aggregate_2d
    INTEGER                                    :: ALLOC_ERR 

    CHARACTER(LEN=50)                          :: S1
    CHARACTER(LEN=250)                         :: msg


! Debug vars
    INTEGER                                    :: nb_coord, nb_var, nb_gat

!_ ================================================================================================================================
    IF (printlev >= 3) WRITE(numout,*)"In interpweight_2Dcont filename: ",TRIM(filename),          &
      " varname:'" // TRIM(varname) // "'"
! Allocating input data
    IF (is_root_prc) THEN
       IF (printlev >=5) THEN
          WRITE(numout,*) "Entering 'interpweight_2Dcont'. Debug mode."
          WRITE(numout,*)'(/," --> fliodmpf")'
          CALL fliodmpf (TRIM(filename))
          WRITE(numout,*)'(/," --> flioopfd")'
       ENDIF
       CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
       IF (printlev >=5) THEN
          WRITE (numout,'(" Number of coordinate        in the file : ",I2)') nb_coord
          WRITE (numout,'(" Number of variables         in the file : ",I2)') nb_var
          WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
       ENDIF
    ENDIF

! Getting shape of input variable from input file
    inNdims = interpweight_get_varNdims_file(filename, varname)
    IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)

    InputDims: SELECT CASE (inNdims)
      CASE (2)
        CALL bcast(iml)
        CALL bcast(jml)

        ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_2Dcont input data dimensions _______'
          WRITE(numout,*)'    iml, jml:',iml, jml
        END IF

      CASE (3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)

        ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_2Dcont input data dimensions:'
          WRITE(numout,*)'    iml, jml, lml:',iml, jml, lml
        END IF

      CASE (4)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)
!       Taken from `slowproc_update'
        IF (initime /= -1) THEN
          inNdims = 3
          IF (printlev >= 3) WRITE(numout,*) &
               'interpweight_2Dcont: taking into account the initial reading time:',initime
          IF (initime <= 0) tml = 0
          
          CALL bcast(tml)
          ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF ( initime > 0 ) THEN
            IF (is_root_prc) THEN
              IF (initime <= tml) THEN
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
              ELSE
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
              ENDIF
            ENDIF
          ELSE
            IF (is_root_prc) THEN
              CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
            ENDIF
          ENDIF
        ELSE
          CALL bcast(tml)
          ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
        END IF

        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_2Dcont input data dimensions:'
          WRITE(numout,*)'    iml, jml, lml, tml :',iml, jml, lml, tml
        END IF

      CASE DEFAULT
        WRITE (S1,'(I3)') inNdims
        CALL ipslerr_p(3,'interpweight_2Dcont',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
          '','')

    END SELECT InputDims
    CALL bcast(invar2D)

    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2Dcont getting variable ended!'

! Getting longitudes/latitudes from input file
    inNLldims = interpweight_get_varNdims_file(filename, inlonname)

    IF (inNLldims == 1) THEN
      ! Allocate memory for latitudes
      ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lon_in','','')

      ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lat_in','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
        CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
      ENDIF
      CALL bcast(lon_in)
      CALL bcast(lat_in)

      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lon','','')
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lat','','')

      DO ip=1,iml
         lat(ip,:) = lat_in(:)
      ENDDO
      DO jp=1,jml
         lon(:,jp) = lon_in(:)
      ENDDO
    ELSE
      ! Allocate memory for longitude
      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Pb in allocation for lon','','')
      ! Allocate memory for latitudes
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Pb in allocation for lat','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
        CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
      ENDIF
    END IF
    CALL bcast(lon)
    CALL bcast(lat)

    IF (is_root_prc) THEN
      ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr(3,'interpweight_2Dcont','Problem in allocation of variable resol_in','','')
    END IF

    IF (noneg) THEN
      IF (is_root_prc) THEN
        SELECT CASE (inNdims)
          CASE(2)
            CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
          CASE(3)
            CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
          CASE(4)
            CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
        END SELECT       
      END IF
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2Dcont modifying input ended!'

    !
    ! Consider all points a priori
    !
    ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable mask','','')
    mask(:,:) = 0

    IF (TRIM(masktype) == 'var') THEN
      ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable maskvar','','')
      maskvar(:,:) = 0
! Reads the mask from the file
      IF (is_root_prc) THEN
! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
        CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
        WHERE (maskvar > zero)
          mask = un 
        END WHERE
      END IF
      CALL bcast(mask)
    ELSE
      SELECT CASE (inNdims)
        CASE(2)
          CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D,      &
            mask)
        CASE(3)
          CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D,    &
            mask)
        CASE(4)
          CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D,  &
            mask)
      END SELECT
    END IF
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2Dcont masking input ended!'
    IF (is_root_prc) CALL flinclo(fid)

    IF (printlev>=5) WRITE(numout,*)"  interpweight_2Dcont '" // TRIM(typefrac) // "' interpolation"

! Assuming regular lon/lat projection for the input data
    IF (is_root_prc) THEN
      CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
    END IF

!   Computation of nbvmax
!   In 'condveg_soilalb, slowproc_update' nbvmax=200

!
!   The number of maximum vegetation map points in the GCM grid is estimated.
!   Some lmargin is taken.
!
    IF (is_root_prc) THEN
       nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
       njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
       nbvmax = MAX(nix*njx, 200)
    ENDIF
    CALL bcast(nbvmax)
    CALL bcast(nix)
    CALL bcast(njx)

    callsign = TRIM(varname) // ' map'

    ok_interpol = .FALSE.
    DO WHILE ( .NOT. ok_interpol )
       IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
       IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx

       ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable sub_index','','')

       sub_index(:,:,:)=0

       ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable sub_area','','')

       sub_area(:,:)=zero
       IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon), MAXVAL(lon)
       IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat), MAXVAL(lat)

       CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
            &                iml, jml, lon, lat, mask, callsign, &
            &                nbvmax, sub_index, sub_area, ok_interpol)
       
       IF ( .NOT. ok_interpol ) THEN
          DEALLOCATE(sub_area)
          DEALLOCATE(sub_index)
          nbvmax = nbvmax * 2
       ENDIF
    ENDDO
    IF (printlev>=4) WRITE (numout,*) '  interpweight_2Dcont: read/allocate OK'

    outvar2D = zero
! Providing the interpolated grid point

    CALL interpweight_provide_interpolation2D(typefrac, nbpt, 0, 0, iml, jml, lml, tml,               &
      nbvmax, zero, invar2D, sub_area, sub_index, varmax, varmin, un, deux, lalo,                     &
      defaultvalue, defaultNOvalue, outvar2D, aoutvar)
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_2Dcont end of interpweight_provide_interpolation2D'

    IF (printlev>=2) THEN
       IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
          WRITE(numout,*) 'interpweight_2Dcont total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' non-interpolated: ', COUNT(aoutvar == -1)
          msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) //       &
               "' have been succesfully interpolated !!"
          WRITE(numout,*) '  ' // TRIM(msg)
       ELSE
          msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
          WRITE(numout,*) "interpweight_2Dcont: '" // TRIM(msg)
          WRITE(numout,*) '  Total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
       END IF
    END IF

    ! Scaling to 1 the quality variable
    DO ip=1, nbpt
      
       !+++HACK+++
       IF (aoutvar(ip) /= -1) THEN
          IF (contfrac(ip) == zero) THEN
             ! HACK: This should not happen but it does sometimes.
             ! We can't calculate aoutvar, therefor assign a value -2
             aoutvar(ip) = -2
          ELSE
             ! The original code assumes that for every pixels that needs
             ! to be considered by ORCHIDEE the land-sea mask has a value
             ! This is a fair assumption. This assumption should be checked
             ! earlier on in the code. This is the original code.
             aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
          END IF
       END IF
       !++++++++++++
    END DO

    IF (printlev>=5) WRITE(numout,*) '  interpweight_2Dcont: Interpolation Done'

    IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
    IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
    DEALLOCATE(lon)
    DEALLOCATE(lat)
    IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
    IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
    IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
    DEALLOCATE(mask)
    IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
    IF (ALLOCATED(invar4D)) DEALLOCATE(sub_area)
    IF (ALLOCATED(invar4D)) DEALLOCATE(sub_index)
    IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)

    RETURN

  END SUBROUTINE interpweight_2Dcont

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_4Dcont
!!
!>\BRIEF        ! Interpolate any input file to the grid of the model for a continuos 4D-variable  
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): outvar4D, aoutvar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_4Dcont(nbpt, dim1, dim2, lalo, resolution, neighbours,                      &
    contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype,               &
    maskvalues, maskvarname, initime, typefrac, defaultvalue, defaultNOvalue,                         &
    outvar4D, aoutvar)

    USE ioipsl
    USE constantes_var
    USE mod_orchidee_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !  0.1 INPUT
    !
    INTEGER(i_std), INTENT(in)                 :: nbpt                    !! Number of points for which the data 
                                                                          !!   needs to be interpolated
    INTEGER(i_std), INTENT(in)                 :: dim1, dim2              !! 4D complementary dimension sizes
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo                    !! Vector of latitude and longitudes 
                                                                          !! (beware of the order = 1 : latitude, 
                                                                          !!   2 : longitude)
    REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution              !! The size in km of each grid-box in X and Y
    INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours    !! Vector of neighbours for each grid point 
                                                                          !!   1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) 
                            
    REAL(r_std), DIMENSION(nbpt), INTENT(in)   :: contfrac                !! Fraction of land in each grid box.
    CHARACTER(LEN=80), INTENT(in)              :: filename                !! name of the input map to read

    CHARACTER(LEN=80), INTENT(in)              :: varname                 !! name of the variable to interpolate
    CHARACTER(LEN=80), INTENT(in)              :: inlonname, inlatname    !! names of the longitude and latitude
    REAL(r_std), INTENT(in)                    :: varmin, varmax          !! min/max values to use for the
                                                                          !!   renormalization
    REAL(r_std), DIMENSION(dim1,dim2), INTENT(in) :: defaultvalue         !! default interpolated value 
    REAL(r_std), INTENT(in)                    :: defaultNOvalue          !! default interpolated NO-value 

! Transformations to apply to the original read values from the file
    LOGICAL, INTENT(in)                        :: noneg                   !! no negative values
    CHARACTER(LEN=50), INTENT(in)              :: masktype                !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below maskvalues(1)
                                                                          !!   'mabove': take values above maskvalues(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
                                                                          !!      maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
                                                                          !!        (normalized by maskedvalues(3))
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (1/0)
    REAL(r_std), DIMENSION(3), INTENT(in)      :: maskvalues              !! values to use to mask (according to 
                                                                          !!   `masktype') 
    CHARACTER(LEN=250), INTENT(in)             :: maskvarname             !! name of the variable to use to mask 
                                                                          !!   (according to `masktype') 
    CHARACTER(LEN=50), INTENT(in)              :: typefrac                !! Type of fraction retrieval:
                                                                          !!   'default': standard

    INTEGER(i_std), INTENT(in)                 :: initime                 !! Initial time to take (for input data 
                                                                          !!   with time values, -1: not used)
    !! 0.2 Modified variables
    !
    !! 0.3 Output variables
    !
    REAL(r_std), DIMENSION(nbpt,dim1,dim2), INTENT(out) :: outvar4D       !! 2D output variable and re-dimensioned
    REAL(r_std), DIMENSION(nbpt), INTENT(out)  :: aoutvar                 !! availability of input data to
                                                                          !!   interpolate output variable 
                                                                          !!   (on the nbpt space)
    !
    !! 0.4 Local variables
    !
    INTEGER(i_std)                             :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
    INTEGER                                    :: inNdims                 !! Number of dimensions of the input data
    INTEGER                                    :: inNLldims               !! Number of dimensions of the lon/lat input data
    REAL(r_std), ALLOCATABLE, DIMENSION(:)     :: lat_in, lon_in          !! latitude and
                                                                          !! longitude, extract from input map
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: lat, lon                !! en 2D ???
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)   :: sub_area                !! the area of the fine grid in the 
                                                                          !!   model grid ???
                                                                          !! cf src_global/interpol_help.f90, 
                                                                          !!   line 377, called "areaoverlap"
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index           !! the indexes from the grid boxes from the 
                                                                          !!   data that go into the 
                                                                          !!   model's boxes  
                                                                          !! cf src_global/interpol_help.f90,
                                                                          !!   line 300, called "ip"

!    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: invar1D               !! 1D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: invar2D               !! 2D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: invar3D               !! 3D input values
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D               !! 4D input values

    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)   :: resol_in
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: maskvar
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)  :: mask
    INTEGER(i_std)                             :: nix, njx
    INTEGER(i_std)                             :: idi, nbvmax             !! nbvmax : number of maximum vegetation map
                                                                          !! points in the GCM grid ; idi : its counter
    CHARACTER(LEN=30)                          :: callsign                !! Allows to specify which variable is beeing
                                                                          !! treated
    LOGICAL                                    :: ok_interpol             !! optionnal return of aggregate_2d
    INTEGER                                    :: ALLOC_ERR 

    CHARACTER(LEN=50)                          :: S1
    INTEGER, ALLOCATABLE, DIMENSION(:)         :: indims
    CHARACTER(LEN=250)                         :: msg

! Debug vars
    INTEGER                                    :: nb_coord, nb_var, nb_gat
    INTEGER                                    :: nblL
    INTEGER, DIMENSION(2)                      :: lLpt

!_ ================================================================================================================================
    IF (printlev >= 3) WRITE(numout, *)"In interpweight_4Dcont filename: '",TRIM(filename),         &
      " varname:'" // TRIM(varname) // "'"
! Allocating input data
    IF (is_root_prc) THEN
       IF (printlev >=5) THEN
          WRITE(numout,*) "Entering 'interpweight_4Dcont'. Debug mode."
          WRITE(numout,*)'(/," --> fliodmpf")'
          CALL fliodmpf (TRIM(filename))
          WRITE(numout,*)'(/," --> flioopfd")'
       ENDIF
       CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
       IF (printlev >=5) THEN
          WRITE (numout,'(" Number of coordinate        in the file : ",I2)') nb_coord
          WRITE (numout,'(" Number of variables         in the file : ",I2)') nb_var
          WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
       ENDIF
    ENDIF

! Getting shape of input variable from input file
    inNdims = interpweight_get_varNdims_file(filename, varname)
    IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)

    InputDims: SELECT CASE (inNdims)
      CASE (2)
        CALL bcast(iml)
        CALL bcast(jml)

        ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_4Dcont input data dimensions:'
          WRITE(numout,*)'    iml, jml:',iml, jml
        END IF

      CASE (3)
        ALLOCATE(indims(3),STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable indims','','')
        indims = interpweight_get_var3dims_file(filename, varname)
        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)

        ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" //     &
          TRIM(varname) ,'','')
        IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
        IF (printlev >= 5) THEN
          WRITE(numout,*)'  interpweight_4Dcont input data dimensions:'
          WRITE(numout,*)'    iml, jml, lml:',iml, jml, lml
        END IF

      CASE (4)
        ALLOCATE(indims(4), STAT=ALLOC_ERR)
        IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of indims','','')
        indims = interpweight_get_var4dims_file(filename, varname)
        lml = indims(3)
        CALL bcast(iml)
        CALL bcast(jml)
        CALL bcast(lml)
!       Taken from `slowproc_update'
        IF (initime /= -1) THEN
          inNdims = 3
          IF (printlev >= 3) WRITE(numout,*) &
               'interpweight_4Dcont: taking into account the initial reading time:',initime
          IF (initime <= 0) tml = 0

          CALL bcast(tml)
          ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF ( initime > 0 ) THEN
            IF (is_root_prc) THEN
              IF (initime <= tml) THEN
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
              ELSE
                CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
              ENDIF
            ENDIF
          ELSE
            IF (is_root_prc) THEN
              CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
            ENDIF
          ENDIF
        ELSE
          CALL bcast(tml)
          ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" //   &
            TRIM(varname) ,'','')
          IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
        END IF

        IF (printlev >= 5) THEN
          WRITE(numout,*)'interpweight_4Dcont input data dimensions:'
          WRITE(numout,*)'    iml, jml, lml, tml:',iml, jml, lml, tml
        END IF

      CASE DEFAULT
        WRITE (S1,'(I3)') inNdims
        CALL ipslerr_p(3,'interpweight_4Dcont',"Input number of dimensions: '" // S1 // "' not ready !!",'','')

    END SELECT InputDims
    CALL bcast(invar4D)

    IF (printlev >= 5) WRITE(numout,*)'  interpweight_4Dcont getting variable ended!'

! Getting longitudes/latitudes from input file
    inNLldims = interpweight_get_varNdims_file(filename, inlonname)

    IF (inNLldims == 1) THEN
      ! Allocate memory for latitudes
      ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lon_in','','')

      ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lat_in','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
        CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
      ENDIF
      CALL bcast(lon_in)
      CALL bcast(lat_in)

      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lon','','')
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lat','','')

      DO ip=1,iml
         lat(ip,:) = lat_in(:)
      ENDDO
      DO jp=1,jml
         lon(:,jp) = lon_in(:)
      ENDDO
    ELSE
      ! Allocate memory for longitude
      ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Pb in allocation for lon','','')
      ! Allocate memory for latitudes
      ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Pb in allocation for lat','','')

      IF (is_root_prc) THEN
        CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
        CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
      ENDIF
    END IF
    CALL bcast(lon)
    CALL bcast(lat)

    IF (is_root_prc) THEN
      ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
      IF (ALLOC_ERR /= 0) CALL ipslerr(3,'interpweight_4Dcont','Problem in allocation of variable resol_in','','')
    END IF

    IF (noneg) THEN
      IF (is_root_prc) THEN
        SELECT CASE (inNdims)
          CASE(2)
            CALL interpweight_modifying_input2D(iml, jml, 0, 0, 0, 0, zero, invar2D)
          CASE(3)
            CALL interpweight_modifying_input3D(iml, jml, lml, 0, 0, 0, zero, invar3D)
          CASE(4)
            CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
        END SELECT       
      END IF
    END IF
    CALL bcast(invar4D)
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_4Dcont modifying input ended!'

    !
    ! Consider all points a priori
    !
    ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable mask','','')
    mask(:,:) = zero

    IF (TRIM(masktype) == 'var') THEN
      ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
      maskvar(:,:) = zero
      IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable maskvar','','')
      ! Reads the mask from the file
      IF (is_root_prc) THEN
         ! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
        CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
        WHERE (maskvar > zero)
          mask = un 
        ENDWHERE
      END IF
      CALL bcast(mask)
    ELSE
      SELECT CASE (inNdims)
        CASE(2)
          CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D,      &
            mask)
        CASE(3)
          CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D,    &
            mask)
        CASE(4)
          CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D,  &
            mask)
      END SELECT
    END IF

    IF (is_root_prc) CALL flinclo(fid)

! Assuming regular lon/lat projection for the input data
    IF (is_root_prc) THEN
      CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
    END IF

    ! Computation of nbvmax
    ! In 'condveg_soilalb, slowproc_update' nbvmax=200

    !
    ! The number of maximum vegetation map points in the GCM grid is estimated.
    ! Some lmargin is taken.
    !
    IF (is_root_prc) THEN
       nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
       njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
       nbvmax = MAX(nix*njx, 200)
    ENDIF
    CALL bcast(nbvmax)
    CALL bcast(nix)
    CALL bcast(njx)

    !
    callsign = TRIM(varname) // ' map'
    !
    ok_interpol = .FALSE.
    DO WHILE ( .NOT. ok_interpol )
       IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
       IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx


       ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable sub_index','','')

       sub_index(:,:,:)=0

       ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
       IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable sub_area','','')

       sub_area(:,:)=zero
       IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon), MAXVAL(lon)
       IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat), MAXVAL(lat)

       CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
            &                iml, jml, lon, lat, mask, callsign, &
            &                nbvmax, sub_index, sub_area, ok_interpol)
       
       IF ( .NOT. ok_interpol ) THEN
          DEALLOCATE(sub_area)
          DEALLOCATE(sub_index)
          nbvmax = nbvmax * 2
       ENDIF
    ENDDO
    IF (printlev>=4) WRITE (numout,*) '  interpweight_4Dcont: read/allocate OK'

    outvar4D = zero

! Providing the interpolated grid point
    CALL interpweight_provide_interpolation4D(typefrac, nbpt, dim1, dim2, iml, jml, lml, tml,         &
      nbvmax, zero, invar4D, sub_area, sub_index, varmax, varmin, un, deux, lalo,                     &
      defaultvalue, defaultNOvalue, outvar4D, aoutvar)
    IF (printlev >= 5) WRITE(numout,*)'  interpweight_4Dcont end of interpweight_provide_interpolation4D'
    
    IF (printlev>=2) THEN
       IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
          WRITE(numout,*) 'interpweight_4Dcont total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' non-interpolated: ', COUNT(aoutvar == -1)
          msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) //       &
               "' have been succesfully interpolated !!"
          WRITE(numout,*) '  ' // TRIM(msg)
       ELSE
          msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
          WRITE(numout,*) "interpweight_4Dcont: '" // TRIM(msg)
          WRITE(numout,*) '  Total number of points: ', nbpt
          WRITE(numout,*) '  Interpolated: ', COUNT(aoutvar /= -1),                 &
               ' Non-interpolated: ', COUNT(aoutvar == -1)
       END IF
    END IF

! Scaling to 1 the quality variable
    DO ip=1, nbpt
       IF (contfrac(ip) == zero) THEN
          ! HACK: This should not happen but it does sometimes.
          ! We can't calculate aoutvar, therefor assign a value -2
          aoutvar(ip) = -2
       ELSE
          IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
       END IF
    END DO

    IF (printlev>=5) WRITE(numout,*) '  interpweight_4Dcont: Interpolation Done'

    IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
    IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
    DEALLOCATE(lon)
    DEALLOCATE(lat)
    IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
    IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
    IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
    DEALLOCATE(mask)
    IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
    DEALLOCATE(sub_area)
    DEALLOCATE(sub_index)
    IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)

    RETURN

  END SUBROUTINE interpweight_4Dcont

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_calc_resolution_in
!!
!>\BRIEF        ! Subroutine to compute the resolution of the input data
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): resolin
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_calc_resolution_in(lons, lats, dx, dy, mcos, REarth, piv, resolin)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !
    !! 0.1 Input variables
    !
    INTEGER(i_std), INTENT(in)                           :: dx, dy        !! Resolution input data
    REAL(r_std), INTENT(in)                              :: mcos          !! cosinus
    REAL(r_std), INTENT(in)                              :: REarth        !! Earth Radius
    REAL(r_std), INTENT(in)                              :: piv           !! pi
    REAL(r_std), DIMENSION(dx,dy), INTENT(in)            :: lons, lats    !! longitude and latitudes (degrees)
    !! 0.2 Modified variables
    !
    !! 0.3 Output variables
    !
    REAL(r_std), DIMENSION(dx,dy,2), INTENT(out)         :: resolin       !! distance between grid points (km)
    !
    !! 0.4 Local variables
    !
    INTEGER                                              :: ip,jp
    REAL(r_std)                                          :: coslat

    DO ip=1,dx
       DO jp=1,dy
          !
          ! Resolution in longitude
          !
          coslat = MAX( COS( lats(ip,jp) * piv/180. ), mcos )     
          IF ( ip .EQ. 1 ) THEN
             resolin(ip,jp,1) = ABS( lons(ip+1,jp) - lons(ip,jp) ) * piv/180. * REarth * coslat
          ELSEIF ( ip .EQ. dx ) THEN
             resolin(ip,jp,1) = ABS( lons(ip,jp) - lons(ip-1,jp) ) * piv/180. * REarth * coslat
          ELSE
             resolin(ip,jp,1) = ABS( lons(ip+1,jp) - lons(ip-1,jp) )/2. * piv/180. * REarth * coslat
          ENDIF
          !
          ! Resolution in latitude
          !
          IF ( jp .EQ. 1 ) THEN
             resolin(ip,jp,2) = ABS( lats(ip,jp) - lats(ip,jp+1) ) * piv/180. * REarth
          ELSEIF ( jp .EQ. dy ) THEN
             resolin(ip,jp,2) = ABS( lats(ip,jp-1) - lats(ip,jp) ) * piv/180. * REarth
          ELSE
             resolin(ip,jp,2) =  ABS( lats(ip,jp-1) - lats(ip,jp+1) )/2. * piv/180. * REarth
          ENDIF
          !
       ENDDO
    ENDDO

  END SUBROUTINE interpweight_calc_resolution_in

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_modifying_input1D
!!
!>\BRIEF        ! modify the initial 1D values from the given file
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ivar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_modifying_input1D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)

    USE constantes_var
    USE ioipsl_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dim1, dim2, dim3, dim4  !! size of the dimensions to get
    INTEGER, INTENT(in)                                  :: stime, etime  !! starting/ending time step to get
    REAL(r_std), INTENT(in)                              :: zeroval       !! zero value
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dim1), INTENT(inout)          :: ivar          !! modified input data


    WHERE (ivar(:) < zeroval )
      ivar(:) = zeroval
    END WHERE

  END SUBROUTINE interpweight_modifying_input1D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_modifying_input2D
!!
!>\BRIEF        ! modify the initial 2D values from the given file
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ivar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_modifying_input2D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)

    USE constantes_var
    USE ioipsl_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dim1, dim2, & !! size of the dimensions to get
      dim3, dim4
    INTEGER, INTENT(in)                                  :: stime, etime  !! starting/ending time step to get
    REAL(r_std), INTENT(in)                              :: zeroval       !! zero value
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dim1,dim2), INTENT(inout)     :: ivar          !! modified input data

    WHERE (ivar(:,:) < zeroval )
      ivar(:,:) = zeroval
    END WHERE

  END SUBROUTINE interpweight_modifying_input2D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_modifying_input3D
!!
!>\BRIEF        ! modify the initial 3D values from the given file
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ivar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_modifying_input3D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)

    USE constantes_var
    USE ioipsl_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dim1, dim2, & !! size of the dimensions to get
      dim3, dim4
    INTEGER, INTENT(in)                                  :: stime, etime  !! starting/ending time step to get
    REAL(r_std), INTENT(in)                              :: zeroval       !! zero value
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dim1,dim2,dim3), INTENT(inout):: ivar          !! modified input data

    WHERE (ivar(:,:,:) < zeroval )
      ivar(:,:,:) = zeroval
    END WHERE

  END SUBROUTINE interpweight_modifying_input3D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_modifying_input4D
!!
!>\BRIEF        ! modify the initial 4D values from the given file
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ivar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_modifying_input4D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)

    USE constantes_var
    USE ioipsl_para

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dim1, dim2, & !! size of the dimensions to get
      dim3, dim4
    INTEGER, INTENT(in)                                  :: stime, etime  !! starting/ending time step to get
    REAL(r_std), INTENT(in)                              :: zeroval       !! zero value
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dim1,dim2,dim3,dim4), INTENT(inout) :: ivar          !! modified input data

    WHERE (ivar(:,:,:,:) < zeroval )
      ivar(:,:,:,:) = zeroval
    END WHERE

  END SUBROUTINE interpweight_modifying_input4D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_interpweight_1D
!!
!>\BRIEF        ! mask 1D input values
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): msk
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_masking_input1D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var,                &
    msk)

    USE constantes_var

    IMPLICIT NONE
  
    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dx, dy, d1, & !! shapes of the input variables
      d2, d3, d4
    CHARACTER(LEN=50), INTENT(in)                        :: mtype         !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below values(1)
                                                                          !!   'mabove': take values above values(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      values(2) <= SUM(vals(k)) <= values(1)
                                                                          !!      values(1) < SUM(vals(k)) <= values(3): 
                                                                          !!        renormalize
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)                :: mvalues       !! values to use to mask (according to `mtype') 
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dx), INTENT(inout)            :: var           !! variable from which provide the mask
    INTEGER(i_std), DIMENSION(dx), INTENT(inout)         :: msk           !! mask to retrieve
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,l

    SELECT CASE (mtype)
      CASE('nomask')
        msk=un
      CASE('mbelow')
        ! e.g.: 
        ! Exclude the points where there is never a LAI value. It is probably 
        ! an ocean point.
        !
        DO i=1,dx
          IF ( var(i) < mvalues(1) ) msk(i) = un
        END DO
      CASE('mabove')
        DO i=1,dx
          IF ( var(i) > mvalues(1) ) msk(i) = un 
        ENDDO
      CASE('msumrange')
        CALL ipslerr_p(3, 'interpweight_masking_input1D',"'msumrange' no sens on a 1D variable",'','')
      CASE DEFAULT
        CALL ipslerr_p(3, 'interpweight_masking_input1D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
    END SELECT

  END SUBROUTINE interpweight_masking_input1D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_masking_input2D
!!
!>\BRIEF        ! mask 2D input values
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): msk
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_masking_input2D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var,                &
    msk)

    USE constantes_var

    IMPLICIT NONE
  
    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dx, dy, d1, & !! shapes of the input variables
      d2, d3, d4
    CHARACTER(LEN=50), INTENT(in)                        :: mtype         !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below values(1)
                                                                          !!   'mabove': take values above values(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      values(2) <= SUM(vals(k)) <= values(1)
                                                                          !!      values(1) < SUM(vals(k)) <= values(3): 
                                                                          !!        renormalize
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)                :: mvalues       !! values to use to mask (according to `mtype') 
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dx,dy), INTENT(inout)         :: var           !! variable from which provide the mask
    INTEGER(i_std), DIMENSION(dx,dy), INTENT(inout)      :: msk           !! mask to retrieve
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,l
    REAL(r_std)                                          :: sumval


    SELECT CASE (mtype)
      CASE('nomask')
        msk=un
      CASE('mbelow')
        ! e.g.: 
        ! Exclude the points where there is never a LAI value. It is probably 
        ! an ocean point.
        !
        DO i=1,dx
          DO j=1,dy
            IF ( var(i,j) < mvalues(1) ) msk(i,j) = un
          END DO
        END DO
      CASE('mabove')
        DO i=1,dx
          DO j=1,dy
            IF ( var(i,j) > mvalues(1) ) msk(i,j) = un 
          ENDDO
        ENDDO
      CASE('msumrange')
        IF (printlev>=3) WRITE(numout,*) &
             'interpweight_masking_input2D: masking with mask sum range. Interval:', mvalues
         DO i=1,dx
           DO j=1,dy
              sumval=SUM(var(i,:))
              IF ( sumval .GE. mvalues(2) .AND. sumval .LE. mvalues(1) ) THEN
                 msk(i,j) = un 
              ELSEIF ( sumval .GT. mvalues(1) .AND. sumval .LE. mvalues(3) ) THEN
                 ! normalization
                 var(i,j) = var(i,j) / sumval
                 msk(i,j) = un 
              ENDIF
           ENDDO
        ENDDO
      CASE DEFAULT
        CALL ipslerr_p(3, 'interpweight_masking_input2D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
    END SELECT

  END SUBROUTINE interpweight_masking_input2D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_masking_input3D
!!
!>\BRIEF        ! mask 3D input values
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): msk
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_masking_input3D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var,                &
    msk)

    USE constantes_var

    IMPLICIT NONE
  
    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dx, dy, d1, & !! shapes of the input variables
      d2, d3, d4
    CHARACTER(LEN=50), INTENT(in)                        :: mtype         !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below values(1)
                                                                          !!   'mabove': take values above values(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      values(2) <= SUM(vals(k)) <= values(1)
                                                                          !!      values(1) < SUM(vals(k)) <= values(3): 
                                                                          !!        renormalize
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)                :: mvalues       !! values to use to mask (according to `mtype') 
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dx,dy,d3), INTENT(inout)      :: var           !! variable from which provide the mask
    INTEGER(i_std), DIMENSION(dx,dy), INTENT(inout)      :: msk           !! mask to retrieve
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,l
    REAL(r_std)                                          :: sumval


    SELECT CASE (mtype)
      CASE('nomask')
        msk=un
      CASE('mbelow')
        ! e.g.: 
        ! Exclude the points where there is never a LAI value. It is probably 
        ! an ocean point.
        !
        DO i=1,dx
          DO j=1,dy
            IF ( ANY(var(i,j,:) < mvalues(1)) ) msk(i,j) = un
          END DO
        END DO
      CASE('mabove')
        DO i=1,dx
          DO j=1,dy
            IF ( ANY(var(i,j,:) > mvalues(1)) ) msk(i,j) = un 
          ENDDO
        ENDDO
      CASE('msumrange')
         DO i=1,dx
           DO j=1,dy
             sumval=SUM(var(i,j,:))
             IF ( sumval .GE. mvalues(2) .AND. sumval .LE. mvalues(1)) THEN
                msk(i,j) = un 
             ELSEIF ( sumval .GT. mvalues(1) .AND. sumval .LE. mvalues(3)) THEN
                ! normalization
                var(i,j,:) = var(i,j,:) / sumval
                msk(i,j) = un 
             ENDIF
           ENDDO
        ENDDO
      CASE DEFAULT
        CALL ipslerr_p(3, 'interpweight_masking_input3D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
    END SELECT
  END SUBROUTINE interpweight_masking_input3D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_masking_input4D
!!
!>\BRIEF        ! mask 4D input values
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): msk
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_masking_input4D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var,                &
    msk)

    USE constantes_var

    IMPLICIT NONE
  
    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: dx, dy, d1, & !! shapes of the input variables
      d2, d3, d4
    CHARACTER(LEN=50), INTENT(in)                        :: mtype         !! Type of masking
                                                                          !!   'nomask': no-mask is applied
                                                                          !!   'mbelow': take values below values(1)
                                                                          !!   'mabove': take values above values(1)
                                                                          !!   'msumrange': take values within 2 ranges;
                                                                          !!      values(2) <= SUM(vals(k)) <= values(1)
                                                                          !!      values(1) < SUM(vals(k)) <= values(3): 
                                                                          !!        renormalize
                                                                          !!   'var': mask values are taken from a 
                                                                          !!     variable (>0)
    REAL(r_std), DIMENSION(3), INTENT(in)                :: mvalues       !! values to use to mask (according to `mtype') 
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(dx,dy,d3,d4), INTENT(inout)   :: var           !! variable from which provide the mask
    INTEGER(i_std), DIMENSION(dx,dy), INTENT(inout)      :: msk           !! mask to retrieve
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,l
    REAL(r_std)                                          :: sumval


    SELECT CASE (mtype)
      CASE('nomask')
        msk=un
      CASE('mbelow')
        ! e.g.: 
        ! Exclude the points where there is never a LAI value. It is probably 
        ! an ocean point.
        !
        DO i=1,dx
          DO j=1,dy
            IF ( ANY(var(i,j,:,:) < mvalues(1)) ) msk(i,j) = un
          END DO
        END DO
      CASE('mabove')
        DO i=1,dx
          DO j=1,dy
            IF ( ANY(var(i,j,:,:) > mvalues(1)) ) msk(i,j) = un 
          ENDDO
        ENDDO
      CASE('msumrange')
        IF (printlev>=3) WRITE(numout,*) &
             'interpweight_masking_input4D: masking with mask sum range. Interval:', mvalues
         DO i=1,dx
           DO j=1,dy
             DO l=1,d4
               sumval=SUM(var(i,j,:,l))
               IF ( sumval .GE. mvalues(2) .AND. sumval .LE. mvalues(1) .AND. msk(i,j) /= un) THEN
                  msk(i,j) = un 
               ELSEIF ( sumval .GT. mvalues(1) .AND. sumval .LE. mvalues(3)) THEN
                  ! normalization
                  var(i,j,:,l) = var(i,j,:,l) / sumval
                  msk(i,j) = un 
               ENDIF
             ENDDO
           ENDDO
        ENDDO
      CASE DEFAULT
        CALL ipslerr_p(3, 'interpweight_masking_input4D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
    END SELECT

  END SUBROUTINE interpweight_masking_input4D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_provide_fractions1D
!!
!>\BRIEF        ! provide fractions from a 1D incoming variable
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ovar, aovar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_provide_fractions1D(tint, npt, Ntypes, valstypes,                           &
    dx, dy, d3, d4, nbmax, zeroval, ivar1D, sarea, sindex, vmax, vmin,                                &
    two, latlon, ovar, aovar)

    USE constantes_var

    IMPLICIT NONE
  
    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=50), INTENT(in)                        :: tint          !! tint: type of interpolation
                                                                          !!   'XYKindTime': Input values are kinds 
                                                                          !!     of something with a temporal 
                                                                          !!     evolution on the dx*dy matrix
    INTEGER, INTENT(in)                                  :: npt           !! total number of grid points
    INTEGER, INTENT(in)                                  :: dx,dy,d3,d4   !! shape of the input values
    INTEGER, INTENT(in)                                  :: nbmax         !! maximum input grid cells to use 
                                                                          !!   for in 1 grid point
    REAL(r_std), INTENT(in)                              :: zeroval, two  !! zero and 2. real values
    INTEGER, INTENT(in)                                  :: Ntypes        !! number of types
    REAL(r_std), DIMENSION(Ntypes), INTENT(in)           :: valstypes     !! value for each type 
    REAL(r_std), DIMENSION(dx), INTENT(in)               :: ivar1D        !! values of the input data
    REAL(r_std), DIMENSION(npt,nbmax), INTENT(in)        :: sarea         !! overlaping area of a given input 
                                                                          !!   grid point
    INTEGER(i_std), DIMENSION(npt,nbmax), INTENT(in)     :: sindex        !! index of the overlaping input 
                                                                          !!   grid point
    REAL(r_std), DIMENSION(npt,2), INTENT(in)            :: latlon        !! longitude and latitude of the destiny
                                                                          !!   grid points
    !! 0.2 Modified variables
    REAL(r_std), INTENT(inout)                           :: vmax, vmin    !! thresholds of the values
    !! 0.3 Output variables
    REAL(r_std), DIMENSION(npt,Ntypes), INTENT(out)      :: ovar          !! output variable (on the nbpt space)
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: aovar         !! availability of input data to 
                                                                          !!   interpolate output variable (on 
                                                                          !!   the nbpt space)
                                                                          !!     >0.: As the total sum of the 
                                                                          !!       space fraction used from the 
                                                                          !!       source
                                                                          !!     -1: Did not find any required values 
                                                                          !!       from input data
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
    CHARACTER(LEN=250)                                   :: msg


! Initialization
!!
    ovar = zeroval

    SELECT CASE (tint)
      CASE ('default')
        IF (printlev>=4) WRITE(numout,*)'  interpweight_provide_fractions1D Fractions following default1D method'
        IF (ALL(ivar1D == zeroval)) THEN
          msg = "' wrong 1D input variable"
          CALL ipslerr_p(3,'interpweight_provide_fractions1D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF 
        DO ib=1,npt
          aovar(ib) = zeroval
          idi = COUNT(sarea(ib,:) > zeroval)
          IF ( idi > 0 ) THEN
            DO jj=1,idi
               ip = sindex(ib,jj)
               DO it = 1, Ntypes
                 IF (ivar1d(ip) == valstypes(it)*un) THEN
                   ovar(ib,it) = ovar(ib,it) + sarea(ib,jj)
                 END IF
               END DO
               aovar(ib) = aovar(ib) + sarea(ib,jj)
            ENDDO
            !
            ! Normalize
            !
            ovar(ib,:) = ovar(ib,:)/aovar(ib)
            !
            ! Quality variable
          ELSE
             ! No points found.
            aovar(ib) = -un
            ovar(ib,INT((vmax+vmin)/two)) = un
          ENDIF
          IF (SUM(ovar(ib,:)) > 1.0000001) THEN
            IF (printlev>=3) WRITE(numout,*) '  interpweight_provide_fractions1D On point ', ib, &
                 ' location: ', latlon(ib,2), latlon(ib,1),' total fractions=', SUM(ovar(ib,:))
            IF (printlev>=3) THEN
               WRITE(numout,*) '  type fraction: '
               DO i=1, Ntypes
                  WRITE(numout,*) i, ovar(ib,i)
               END DO
            END IF
            msg='total of fractions above 1.!'
            CALL ipslerr(3,'interpweight_provide_fractions1D',TRIM(msg),'','')
          END IF
        END DO

        IF (printlev >= 3 .AND. ANY(aovar == -un)) THEN
           WRITE(numout,*) 'interpweight_provide_fractions1D: '
           WRITE(numout,*) 'Some points on the model grid did not have any points on the source grid for interpolation.'
           WRITE(numout,*) 'These points are initialized by putting the average VAR ', INT((vmax+vmin)/two), ' to 1.'
           WRITE(numout,*) 'The points are (ib, lon, lat) :'
           DO ib=1, npt
              IF (aovar(ib) == -un) WRITE(numout,*) ib, latlon(ib,2), latlon(ib,1)
           END DO
        END IF
      CASE DEFAULT
        CALL ipslerr_p(3,'interpweight_provide_fractions1D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
    END SELECT

  END SUBROUTINE interpweight_provide_fractions1D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_provide_fractions2D
!!
!>\BRIEF        ! provide fractions from a 2D incoming variable
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ovar, aovar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_provide_fractions2D(tint, npt, Ntypes, valstypes,                           &
    dx, dy, d3, d4, nbmax, zeroval, ivar2D, sarea, sindex, vmax, vmin,                                &
    two, latlon, ovar, aovar)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=50), INTENT(in)                        :: tint          !! tint: type of interpolation
                                                                          !!   'XYKindTime': Input values are kinds 
                                                                          !!     of something with a temporal 
                                                                          !!     evolution on the dx*dy matrix
    INTEGER, INTENT(in)                                  :: npt           !! total number of grid points
    INTEGER, INTENT(in)                                  :: dx,dy,d3,d4   !! shape of the input values
    INTEGER, INTENT(in)                                  :: nbmax         !! maximum input grid cells to use 
                                                                          !!   for in 1 grid point
    REAL(r_std), INTENT(in)                              :: zeroval, two  !! zero and 2. real values
    INTEGER, INTENT(in)                                  :: Ntypes        !! number of types
    REAL(r_std), DIMENSION(Ntypes), INTENT(in)           :: valstypes     !! value for each type 
    REAL(r_std), DIMENSION(dx,dy), INTENT(in)            :: ivar2D        !! values of the input data
    REAL(r_std), DIMENSION(npt,nbmax), INTENT(in)        :: sarea         !! overlaping area of a given input 
                                                                          !!   grid point
    INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in)   :: sindex        !! index of the overlaping input 
                                                                          !!   grid point
    REAL(r_std), DIMENSION(npt,2), INTENT(in)            :: latlon        !! longitude and latitude of the destiny
                                                                          !!   grid points
    !! 0.2 Modified variables
    REAL(r_std), INTENT(inout)                           :: vmax, vmin    !! thresholds of the values
    !! 0.3 Output variables
    REAL(r_std), DIMENSION(npt,Ntypes), INTENT(out)      :: ovar          !! output variable (on the nbpt space)
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: aovar         !! availability of input data to 
                                                                          !!   interpolate output variable (on 
                                                                          !!   the nbpt space)
                                                                          !!     >0.: As the total sum of the 
                                                                          !!       space fraction used from the 
                                                                          !!       source
                                                                          !!     -1: Found any required values 
                                                                          !!       from input data
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
    CHARACTER(LEN=250)                                   :: msg


! Initialization
!!
    ovar = zeroval

    SELECT CASE (tint)
      CASE ('default')
        IF (printlev>=3) WRITE(numout,*)'  interpweight_provide_fractions2D Fractions following default2D method'
        IF (ALL(ivar2D == zeroval)) THEN
          msg = "' wrong 2D input variable"
          CALL ipslerr_p(3,'interpweight_provide_fractions2D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF 
        DO ib=1,npt
          aovar(ib) = zeroval
          idi = COUNT(sarea(ib,:) > zeroval)
          IF ( idi > 0 ) THEN
            DO jj=1,idi
               ip = sindex(ib,jj,1)
               jp = sindex(ib,jj,2)
               DO it = 1, Ntypes
                 IF (ivar2d(ip,jp) == valstypes(it)*un) THEN
                   ovar(ib,it) = ovar(ib,it) + sarea(ib,jj)
                 END IF
               END DO
               aovar(ib) = aovar(ib) + sarea(ib,jj)
            ENDDO
            !
            ! Normalize
            !
            ovar(ib,:) = ovar(ib,:)/aovar(ib)
            !
            ! Quality variable
          ELSE
            ! No points on the source grid were found for the current point. 
            ! We provide here a default value.
            aovar(ib) = -un
            ovar(ib,INT((vmax+vmin)/two)) = un
          ENDIF

          IF (SUM(ovar(ib,:)) > 1.0000001) THEN
            WRITE(numout,*) '  interpweight_provide_fractions2D On point ', ib, ' location: ', latlon(ib,2),     &
              latlon(ib,1),' total fractions=', SUM(ovar(ib,:))
            WRITE(numout,*) '  type fraction _______'
            DO i=1, Ntypes
              WRITE(numout,*) i, ovar(ib,i)
            END DO
            msg='total of fractions above 1.!'
            CALL ipslerr(3,'interpweight_provide_fractions2D',TRIM(msg),'','')
          END IF
        END DO

        IF (printlev >= 3 .AND. ANY(aovar == -un)) THEN
           WRITE(numout,*) 'interpweight_provide_fractions2D: '
           WRITE(numout,*) 'Some points on the model grid did not have any points on the source grid for interpolation.'
           WRITE(numout,*) 'These points are initialized by putting the average VAR ', INT((vmax+vmin)/two), ' to 1.'
           WRITE(numout,*) 'The points are (ib, lon, lat) :'
           DO ib=1, npt
              IF (aovar(ib) == -un) WRITE(numout,*) ib, latlon(ib,2), latlon(ib,1)
           END DO
        END IF
 
      CASE DEFAULT
        CALL ipslerr_p(3,'interpweight_provide_fractions2D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
    END SELECT

  END SUBROUTINE interpweight_provide_fractions2D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_provide_fractions3D
!!
!>\BRIEF        ! provide fractions from a 3D incoming variable
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ovar, aovar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_provide_fractions3D(tint, npt, Ntypes, valstypes,                           &
    dx, dy, d3, d4, nbmax, zeroval, ivar3D, sarea, sindex, vmax, vmin,                                &
    two, latlon, ovar, aovar)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=50), INTENT(in)                        :: tint          !! tint: type of interpolation
                                                                          !!   'XYKindTime': Input values are kinds 
                                                                          !!     of something with a temporal 
                                                                          !!     evolution on the dx*dy matrix
    INTEGER, INTENT(in)                                  :: npt           !! total number of grid points
    INTEGER, INTENT(in)                                  :: dx,dy,d3,d4   !! shape of the input values
    INTEGER, INTENT(in)                                  :: nbmax         !! maximum input grid cells to use 
                                                                          !!   for in 1 grid point
    REAL(r_std), INTENT(in)                              :: zeroval, two  !! zero and 2. real values
    INTEGER, INTENT(in)                                  :: Ntypes        !! number of types
    REAL(r_std), DIMENSION(Ntypes), INTENT(in)           :: valstypes     !! value for each type 
    REAL(r_std), DIMENSION(dx,dy,d3), INTENT(in)         :: ivar3D        !! values of the input data
    REAL(r_std), DIMENSION(npt,nbmax), INTENT(in)        :: sarea         !! overlaping area of a given input 
                                                                          !!   grid point
    INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in)   :: sindex        !! index of the overlaping input 
                                                                          !!   grid point
    REAL(r_std), DIMENSION(npt,2), INTENT(in)            :: latlon        !! longitude and latitude of the destiny
                                                                          !!   grid points
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(Ntypes), INTENT(inout)        :: vmax, vmin    !! thresholds of the values
    !! 0.3 Output variables
    REAL(r_std), DIMENSION(npt,Ntypes), INTENT(out)      :: ovar          !! output variable (on the nbpt space)
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: aovar         !! availability of input data to 
                                                                          !!   interpolate output variable (on 
                                                                          !!   the nbpt space)
                                                                          !!     >0.: As the total sum of the 
                                                                          !!       space fraction used from the 
                                                                          !!       source
                                                                          !!     -1: Found any required values 
                                                                          !!       from input data
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
    CHARACTER(LEN=250)                                   :: msg
    INTEGER                                              :: Nzeros
    INTEGER                                              :: ALLOC_ERR 
    REAL(r_std), ALLOCATABLE, DIMENSION(:)               :: subfrac
    REAL(r_std)                                          :: sumpt


! Initialization
!!
    ovar = zeroval
    
    IF (Ntypes /= d3) THEN
       WRITE(numout,*) 'Ntypes, ',Ntypes
       WRITE(numout,*) 'd3, ', d3
       msg = 'Different number of types than third dimension in the input data!!'
       CALL ipslerr_p(3,'interpweight_provide_fractions3D',TRIM(msg),'','')
    END IF

    Nzeros = 0
    SELECT CASE (tint)
      CASE ('default')
        IF (printlev>=3) WRITE(numout,*)'Fractions following default3D method'
        IF (ALL(ivar3D == zeroval)) THEN
          msg = "' wrong 3D input variable"
          CALL ipslerr_p(3,'interpweight_provide_fractions3D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF
        DO ib=1,npt
          aovar(ib) = zeroval
          idi = COUNT(sarea(ib,:) > zeroval)
          IF ( idi > 0 ) THEN
            sumpt = zeroval
            DO jj=1,idi
               ip = sindex(ib,jj,1)
               jp = sindex(ib,jj,2)
               DO it = 1, Ntypes
! Do not get that areas with zero input value and with missing value
                 IF (ivar3D(ip,jp,it) > zeroval .AND. ivar3D(ip,jp,it) < two) THEN
                   ovar(ib,it) = ovar(ib,it) + ivar3D(ip,jp,it)*sarea(ib,jj)
                   sumpt = sumpt + ivar3D(ip,jp,it)*sarea(ib,jj)
                 END IF
               END DO
               aovar(ib) = aovar(ib) + sarea(ib,jj)
            ENDDO
            !
            ! Normalize
            !
            ovar(ib,:) = ovar(ib,:)/aovar(ib)
            !
            ! Quality variable
          ELSE
            Nzeros = Nzeros + 1
            aovar(ib) = -un

            DO it = 1, d3
              ovar(ib,INT((vmax(it)+vmin(it))/two)) = un
            END DO
          ENDIF
! Tests
!
          IF (SUM(ovar(ib,:)) > 1.0000001) THEN
            idi = COUNT(sarea(ib,:) > zeroval)
            ALLOCATE(subfrac(idi), STAT=ALLOC_ERR)
            IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_provide_fractions3D','Problem in allocation of subfrac','','')
            IF (printlev>=3) THEN
               WRITE(numout,*) 'interpweight_provide_fractions3D On point ', ib, &
                    ' location: ', latlon(ib,2), latlon(ib,1),' total fractions=', SUM(ovar(ib,:)), &
                    ' total points from invar: ', idi, ' total invar area: ', aovar(ib)
               WRITE(numout,*) '  type ivar3D jj sarea fraction: '

               DO it = 1, Ntypes
                  DO jj=1,idi
                     ip = sindex(ib,jj,1)
                     jp = sindex(ib,jj,2)
                     IF (ivar3D(ip,jp,it) > zeroval) THEN
                        WRITE(numout,*) it, ivar3D(ip,jp,it), jj, sarea(ib,jj), ovar(ib,it)
                     END IF
                  END DO
               END DO
            END IF
            msg='total of fractions above 1.!'
            CALL ipslerr(3,'interpweight_provide_fractions3D',TRIM(msg),'','')
         END IF

        END DO

        IF (printlev >= 3 .AND. ANY(aovar == -un)) THEN
           WRITE(numout,*) 'interpweight_provide_fractions3D: '
           WRITE(numout,*) 'Some points on the model grid did not have any points on the source grid for interpolation.'
           WRITE(numout,*) 'These points are initialized by putting the average VAR ', INT((vmax+vmin)/two), ' to 1.'
           WRITE(numout,*) 'The points are (ib, lon, lat) :'
           DO ib=1, npt
              IF (aovar(ib) == -un) WRITE(numout,*) ib, latlon(ib,2), latlon(ib,1)
           END DO
        END IF

        IF (COUNT(aovar == -1) == npt) THEN
          msg='No grid points could be interpolated from the source grid for the current process'
          CALL ipslerr_p(2,'interpweight_provide_fractions3D',TRIM(msg),'','')
        END IF

        IF (printlev >=3) WRITE(numout,*) 'interpweight_provide_fractions3D: Nzeros=', &
             Nzeros,' aovar=', COUNT(aovar == -1), ' nbpt=', npt

        IF (COUNT(aovar == -1) /= Nzeros) THEN
          msg='Something went wrong COUNT(aovar == -1) /= Nzeros'
          CALL ipslerr_p(3,'interpweight_provide_fractions3D',TRIM(msg),'','')
        END IF
 
      CASE DEFAULT
        CALL ipslerr_p(3,'interpweight_provide_fractions3D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
    END SELECT

  END SUBROUTINE interpweight_provide_fractions3D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_provide_fractions4D
!!
!>\BRIEF        ! provide fractions from a 4D incoming variable
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ovar, aovar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_provide_fractions4D(tint, npt, Ntypes, valstypes,                           &
    dx, dy, d3, d4, nbmax, zeroval, ivar4D, sarea, sindex, vmax, vmin,                                &
    two, latlon, ovar, aovar)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=50), INTENT(in)                        :: tint          !! tint: type of interpolation
                                                                          !!   'XYKindTime': Input values are kinds 
                                                                          !!     of something with a temporal 
                                                                          !!     evolution on the dx*dy matrix
    INTEGER, INTENT(in)                                  :: npt           !! total number of grid points
    INTEGER, INTENT(in)                                  :: dx,dy,d3,d4   !! shape of the input values
    INTEGER, INTENT(in)                                  :: nbmax         !! maximum input grid cells to use 
                                                                          !!   for in 1 grid point
    REAL(r_std), INTENT(in)                              :: zeroval, two  !! zero and 2. real values
    INTEGER, INTENT(in)                                  :: Ntypes        !! number of types
    REAL(r_std), DIMENSION(Ntypes), INTENT(in)           :: valstypes     !! value for each type 
    REAL(r_std), DIMENSION(dx,dy,d3,d4), INTENT(in)      :: ivar4D        !! values of the input data
    REAL(r_std), DIMENSION(npt,nbmax), INTENT(in)        :: sarea         !! overlaping area of a given input 
                                                                          !!   grid point
    INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in)   :: sindex        !! index of the overlaping input 
                                                                          !!   grid point
    REAL(r_std), DIMENSION(npt,2), INTENT(in)            :: latlon        !! longitude and latitude of the destiny
                                                                          !!   grid points
    !! 0.2 Modified variables
    REAL(r_std), DIMENSION(Ntypes), INTENT(inout)        :: vmax, vmin    !! thresholds of the values
    !! 0.3 Output variables
    REAL(r_std), DIMENSION(npt,Ntypes,d4), INTENT(out)   :: ovar          !! output variable (on the nbpt space)
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: aovar         !! availability of input data to 
                                                                          !!   interpolate output variable (on 
                                                                          !!   the nbpt space)
                                                                          !!     >0.: As the total sum of the 
                                                                          !!       space fraction used from the 
                                                                          !!       source
                                                                          !!     -1: Found any required values 
                                                                          !!       from input data
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it,itt
    CHARACTER(LEN=250)                                   :: msg


! Initialization
!!
    ovar = zeroval

    IF (Ntypes /= d3) THEN
      msg = 'Different number of types than third dimension in the input data!!'
      CALL ipslerr_p(3,'interpweight_provide_fractions4D',TRIM(msg),'','')
    END IF

    SELECT CASE (tint)
      CASE ('default')
        IF (printlev>=3) WRITE(numout,*)'Fractions following default4D method'
        IF (ALL(ivar4D == zeroval)) THEN
          msg = "' wrong 4D input variable"
          CALL ipslerr_p(3,'interpweight_provide_fractions4D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF
        DO ib=1,npt
          aovar(ib) = zeroval
          idi = COUNT(sarea(ib,:) > zeroval)
          IF ( idi > 0 ) THEN
            DO jj=1,idi
               ip = sindex(ib,jj,1)
               jp = sindex(ib,jj,2)
               DO it = 1, Ntypes
                 DO itt = 1, d4
                   ! Do not get that areas with zero input value and with missing value
                   IF (ivar4D(ip,jp,it,itt) > zeroval .AND. ivar4D(ip,jp,it,itt) < 20.) THEN
                     ovar(ib,it,itt) = ovar(ib,it,itt) + ivar4D(ip,jp,it,itt)*sarea(ib,jj)
                   END IF
                 END DO
               END DO
               aovar(ib) = aovar(ib) + sarea(ib,jj)
            ENDDO
            !
            ! Normalize
            !
            ovar(ib,:,:) = ovar(ib,:,:)/aovar(ib)
            !
            ! Quality variable
          ELSE
            ! No points found. Set the average of input threashold (vmax, vmin) as value. 
            aovar(ib) = -un
            IF (printlev >= 3) WRITE(numout,*) 'interpweight_provide_fractions4D: ',&
                 'No points found for interpolating of point ib=',ib,&
                 ' with location lon, lat = ',latlon(ib,2), latlon(ib,1)
            DO it=1, Ntypes
              DO itt=1, d4
                ovar(ib,INT((vmax(it)+vmin(it))/two),itt) = un
              END DO
            END DO
          ENDIF
!          DO itt = 1, d4
!             IF (SUM(ovar(ib,:,itt)) > 1.0000001) THEN
!                IF (printlev>=3) THEN
!                   WRITE(numout,*) '  interpweight_provide_fractions4D On point ', ib, &
!                        ' location: ', latlon(ib,2), latlon(ib,1),' itt=', itt, &
!                        ' total fractions=', SUM(ovar(ib,:,itt))
!                   WRITE(numout,*) '  type fraction:'
!                   DO i=1, Ntypes
!                      WRITE(numout,*) i, ovar(ib,:,itt)
!                   END DO
!                END IF
!                msg='total of fractions above 1.!'
!                CALL ipslerr_p(3,'interpweight_provide_fractions4D',TRIM(msg),'','')
!             END IF
!          END DO
       END DO
       
      CASE DEFAULT
        CALL ipslerr_p(3,'interpweight_provide_fractions4D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
    END SELECT

  END SUBROUTINE interpweight_provide_fractions4D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_provide_interpolation2D
!!
!>\BRIEF        ! perform the interpolation for a continuos 2D field
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ovar, aovar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_provide_interpolation2D(tint, npt, d1, d2, dx, dy, d3, d4,                  &
    nbmax, zeroval, ivar2D, sarea, sindex, vmax, vmin, one, two, latlon,                              &
    defaultval, defaultNOval, ovar, aovar)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=50), INTENT(in)                        :: tint          !! tint: type of interpolation
                                                                          !!   'XYKindTime': Input values are kinds 
                                                                          !!     of something with a temporal 
                                                                          !!     evolution on the dx*dy matrix
    INTEGER, INTENT(in)                                  :: npt           !! total number of grid points
    INTEGER, INTENT(in)                                  :: d1, d2, d3, & !! shape of the input values
      d4, dx, dy
    INTEGER, INTENT(in)                                  :: nbmax         !! maximum input grid cells to use 
                                                                          !!   for in 1 grid point
    REAL(r_std), INTENT(in)                              :: zeroval,one,two !! zero 1. and 2. real values
    REAL(r_std), DIMENSION(dx,dy), INTENT(in)            :: ivar2D        !! values of the input data
    REAL(r_std), DIMENSION(npt,nbmax), INTENT(in)        :: sarea         !! overlaping area of a given input 
                                                                          !!   grid point
    INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in)   :: sindex        !! index of the overlaping input 
                                                                          !!   grid point
    REAL(r_std), DIMENSION(npt,2), INTENT(in)            :: latlon        !! longitude and latitude of the destiny
                                                                          !!   grid points
    REAL(r_std), INTENT(in)                              :: defaultval    !! default interpolated value, only usef if tint='default'
    REAL(r_std), INTENT(in)                              :: defaultNOval  !! default interpolated NO value, only used if tint='slopecalc'
    !! 0.2 Modified variables
    REAL(r_std), INTENT(in)                              :: vmax, vmin    !! thresholds of the values
    !! 0.3 Output variables
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: ovar          !! output variable (on the nbpt space)
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: aovar         !! availability of input data to 
                                                                          !!   interpolate output variable (on 
                                                                          !!   the nbpt space)
                                                                          !!     >0.: As the total sum of the 
                                                                          !!       space fraction used from the 
                                                                          !!       source
                                                                          !!     -1: Found any required values 
                                                                          !!       from input data
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
    INTEGER                                              :: idi_last
    REAL(r_std)                                          :: int1D
    CHARACTER(LEN=250)                                   :: msg


! Initialization
!!
    ovar = zeroval

    SELECT CASE (tint)
      CASE ('default')
        IF (printlev>=3) WRITE(numout,*)'Fractions following default method'
        IF (ALL(ivar2D == zeroval)) THEN
          msg = "' wrong 2D input variable"
          CALL ipslerr_p(3,'interpweight_provide_interpolation2D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF

        DO ib = 1, npt
          !-
          !- Reinfiltration coefficient due to the slope: Calculation with parameteres maxlope_ro 
          !-
          int1D = zeroval

          ! Initialize last index to the highest possible 
          idi_last=nbmax
          aovar(ib) = zero
          DO idi=1, nbmax
             ! Leave the do loop if all sub areas are treated, sub_area <= 0
             IF ( sarea(ib,idi) <= zeroval ) THEN
                ! Set last index to the last one used
                idi_last=idi-1
                ! Exit do loop
                EXIT
             END IF

             ip = sindex(ib,idi,1)
             jp = sindex(ib,idi,2)

             int1D = int1D + ivar2D(ip,jp) * sarea(ib,idi)
             aovar(ib) = aovar(ib) + sarea(ib,idi)
          ENDDO

          IF ( idi_last >= 1 ) THEN
             ovar(ib) = int1D / aovar(ib)
          ELSE
             ovar(ib) = defaultval
            !
            ! Quality variable
            aovar(ib) = -1.
          ENDIF
        ENDDO
      CASE ('slopecalc')
        IF (printlev>=3) WRITE(numout,*)'Fractions following slopecalc method'
        IF (ALL(ivar2D == zeroval)) THEN
          msg = "' requires a 2D input variable"
          CALL ipslerr_p(3,'interpweight_provide_interpolation2D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF

        DO ib = 1, npt
          !-
          !- Reinfiltration coefficient due to the slope: Calculation with parameteres maxlope_ro 
          !-
          int1D = zeroval

          ! Initialize last index to the highest possible 
          idi_last=nbmax
          aovar(ib) = zero
          DO idi=1, nbmax
             ! Leave the do loop if all sub areas are treated, sub_area <= 0
             IF ( sarea(ib,idi) <= zeroval ) THEN
                ! Set last index to the last one used
                idi_last=idi-1
                ! Exit do loop
                EXIT
             END IF

             ip = sindex(ib,idi,1)
             jp = sindex(ib,idi,2)

             int1D = int1D + MIN(ivar2D(ip,jp)/defaultNOval,un) * sarea(ib,idi)
             aovar(ib) = aovar(ib) + sarea(ib,idi)
          ENDDO

          IF ( idi_last >= 1 ) THEN
             ovar(ib) = un - int1D / aovar(ib)
          ELSE
             ovar(ib) = defaultval
            !
            ! Quality variable
            aovar(ib) = -1.
          ENDIF
        ENDDO

      CASE DEFAULT
        CALL ipslerr_p(3,'interpweight_provide_interpolation2D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
    END SELECT

  END SUBROUTINE interpweight_provide_interpolation2D

!! ================================================================================================================================
!! SUBROUTINE   : interpweight_provide_interpolation4D
!!
!>\BRIEF        ! perform the interpolation for a continuos 4D field
!!
!! DESCRIPTION  : (definitions, functional, design, flags): 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ovar, aovar
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  SUBROUTINE interpweight_provide_interpolation4D(tint, npt, d1, d2, dx, dy, d3, d4,                  &
    nbmax, zeroval, ivar4D, sarea, sindex, vmax, vmin, one, two, latlon,                              &
    defaultval, defaultNOval, ovar, aovar)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=50), INTENT(in)                        :: tint          !! tint: type of interpolation
                                                                          !!   'XYKindTime': Input values are kinds 
                                                                          !!     of something with a temporal 
                                                                          !!     evolution on the dx*dy matrix
    INTEGER, INTENT(in)                                  :: npt           !! total number of grid points
    INTEGER, INTENT(in)                                  :: d1, d2, d3, & !! shape of the input values
      d4, dx, dy
    INTEGER, INTENT(in)                                  :: nbmax         !! maximum input grid cells to use 
                                                                          !!   for in 1 grid point
    REAL(r_std), INTENT(in)                              :: zeroval,two,one  !! zero 1. and 2. real values
    REAL(r_std), DIMENSION(dx,dy,d3,d4), INTENT(in)      :: ivar4D        !! values of the input data
    REAL(r_std), DIMENSION(npt,nbmax), INTENT(in)        :: sarea         !! overlaping area of a given input 
                                                                          !!   grid point
    INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in)   :: sindex        !! index of the overlaping input 
                                                                          !!   grid point
    REAL(r_std), DIMENSION(npt,2), INTENT(in)            :: latlon        !! longitude and latitude of the destiny
                                                                          !!   grid points
    REAL(r_std), DIMENSION(d1,d2), INTENT(in)            :: defaultval    !! default interpolated value
    REAL(r_std), INTENT(in)                              :: defaultNOval  !! default interpolated NO value
    !! 0.2 Modified variables
    REAL(r_std), INTENT(in)                              :: vmax, vmin    !! thresholds of the values
    !! 0.3 Output variables
    REAL(r_std), DIMENSION(npt,d1,d2), INTENT(out)       :: ovar          !! output variable (on the nbpt space)
    REAL(r_std), DIMENSION(npt), INTENT(out)             :: aovar         !! availability of input data to 
                                                                          !!   interpolate output variable (on 
                                                                          !!   the nbpt space)
                                                                          !!     >0.: As the total sum of the 
                                                                          !!       space fraction used from the 
                                                                          !!       source
                                                                          !!     -1: Found any required values 
                                                                          !!       from input data
    !! 0.4 Local variables
    INTEGER                                              :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
    INTEGER                                              :: ilf, fopt
    REAL(r_std)                                          :: totarea
    CHARACTER(LEN=250)                                   :: msg
! Debug
    INTEGER                                              :: nbptlL
    REAL(r_std), DIMENSION(npt,d1,d2)                    :: ovartst


! Initialization
!!
    ovar = zero
    ovartst = zero

    SELECT CASE (tint)
      CASE ('default')
        IF (printlev>=3) WRITE(numout,*)'Fractions following default method'
        IF ( ALL(ivar4D == zeroval) ) THEN
          msg = "' wrong 4D input variable"
          CALL ipslerr_p(3,'provide_interpolation4D',"'" // TRIM(tint) // TRIM(msg),'','')
        END IF

        DO ib = 1, npt

          fopt = COUNT(sarea(ib,:) > zero)
          aovar(ib) = zero

          IF ( fopt > 0 ) THEN
            !-
            !- Reinfiltration coefficient
            !-
            totarea = zeroval
 
            DO ilf=1, fopt
               ! Leave the do loop if all sub areas are treated, sub_area <= 0
               ip = sindex(ib,ilf,1)
               jp = sindex(ib,ilf,2)

               DO jv=1,d1
                 DO it=1,d2
                   ovar(ib,jv,it) = ovar(ib,jv,it) + ivar4D(ip,jp,jv,it)*sarea(ib,ilf)
                 END DO
               END DO
               totarea = totarea + sarea(ib,ilf)
            ENDDO
            ! Normalize
            ovar(ib,:,:) = ovar(ib,:,:)/totarea
            aovar(ib) = totarea
          ELSE
        ! Set defalut value for points where the interpolation fail
            IF (printlev>=3) THEN
               WRITE(numout,*) 'provide_interpolation4D: no point found !!'
               WRITE(numout,*) 'On point ', ib, ' no points were found for interpolation data. ' //      &
                    'Default values are used.'
               WRITE(numout,*) 'Location : ', latlon(ib,2), latlon(ib,1)
            END IF
            ovar(ib,ivis,:) = defaultval(ivis,:)
            ovar(ib,inir,:) = defaultval(inir,:) 
            aovar(ib) = -un
          END IF
        ENDDO

      CASE DEFAULT
        CALL ipslerr_p(3,'provide_interpolation4D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
    END SELECT

  END SUBROUTINE interpweight_provide_interpolation4D


!!================================================================================================================================
!! FUNCTION     : interpweight_get_var2dims_file
!!
!>\BRIEF        ! Function to get the dimensions of a given 2D variable inside a file
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_get_var2dims_file
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_get_var2dims_file(filename, varname)

    USE netcdf

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=*), INTENT(in)                         :: filename      !! filename: name of the file to open
    CHARACTER(LEN=*), INTENT(in)                         :: varname       !! varname: name of the variable 
    !! 0.2 Modified variables
    !! 0.3 Output variables
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(2)                                :: interpweight_get_var2dims_file
    !! 0.4 Local variables
    INTEGER                                              :: nid, vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(2)                                :: dimsid
    CHARACTER(LEN=250)                                   :: msg


    IF (LEN_TRIM(varname) < 1) THEN
      msg = "  interpweight_get_var2dims_file: any variable name '" // TRIM(varname) // "' was provided !!"
      CALL ipslerr_p(3,'interpweight_get_var2dims_file',TRIM(msg), '', '')
    END IF

    IF (is_root_prc) THEN
       ! Open file for read only mode
       rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_open', &
            TRIM(nf90_strerror(rcode)),'')

       rcode = nf90_inq_varid(nid, varname, vid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_inq_varid', &
            TRIM(nf90_strerror(rcode)),'')

       rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_variable 1', &
            TRIM(nf90_strerror(rcode)),'')
       
       IF (Ndims /= 2) THEN
          msg = "variable '" // TRIM(varname) // "' has not 2 dimensions!!"
          CALL ipslerr(3,'interpweight_get_var2dims_file',TRIM(msg), '', '')
       END IF
       
       rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_variable 2', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = interpweight_get_var2dims_file(1))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_dimension 1', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = interpweight_get_var2dims_file(2))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_dimension 2', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = NF90_CLOSE(nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3, 'interpweight_get_var2dims_file', 'Error in ng90_close', &
            TRIM(nf90_strerror(rcode)),'')
    END IF

    CALL bcast(interpweight_get_var2dims_file)

  END FUNCTION interpweight_get_var2dims_file

!!================================================================================================================================
!! FUNCTION     : interpweight_get_var3dims_file
!!
!>\BRIEF        ! Function to get the dimensions of a given 3D variable inside a file
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_get_var3dims_file
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_get_var3dims_file(filename, varname)

    USE netcdf

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=*), INTENT(in)                         :: filename      !! filename: name of the file to open
    CHARACTER(LEN=*), INTENT(in)                         :: varname       !! varname: name of the variable 
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(3)                                :: interpweight_get_var3dims_file
    !! 0.4 Local variables
    INTEGER                                              :: nid, vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(3)                                :: dimsid
    CHARACTER(LEN=250)                                   :: msg


    IF (LEN_TRIM(varname) < 1) THEN
      msg = "  interpweight_get_var3dims_file: any variable name '" // TRIM(varname) // "' was provided !!"
      CALL ipslerr_p(3,'interpweight_get_var3dims_file',TRIM(msg), '', '')
    END IF

    IF (is_root_prc) THEN
       ! Open file for read only mode
       rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_open', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inq_varid(nid, varname, vid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_inq_varid', &
            TRIM(nf90_strerror(rcode)),'')

       rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_inquire_variable 1', &
            TRIM(nf90_strerror(rcode)),'')
       
       IF (Ndims /= 3) THEN
          msg = "variable '" // TRIM(varname) // "' has not 3 dimensions!!"
          CALL ipslerr(3,'interpweight_get_var3dims_file',TRIM(msg), '', '')
       END IF
       
       rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_inquire_variable 2', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = interpweight_get_var3dims_file(1))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_dimesion 1', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = interpweight_get_var3dims_file(2))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_dimension 2', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = interpweight_get_var3dims_file(3))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_dimension 3', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = NF90_CLOSE(nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var3dims_file', 'Error in ng90_close', &
            TRIM(nf90_strerror(rcode)),'')
    END IF
    CALL bcast(interpweight_get_var3dims_file)

  END FUNCTION interpweight_get_var3dims_file

!!================================================================================================================================
!! FUNCTION     : interpweight_get_var4dims_file
!!
!>\BRIEF        ! Function to get the dimensions of a given 4D variable inside a file
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_get_var4dims_file
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_get_var4dims_file(filename, varname)

    USE netcdf

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=*), INTENT(in)                         :: filename      !! filename: name of the file to open
    CHARACTER(LEN=*), INTENT(in)                         :: varname       !! varname: name of the variable 
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(4)                                :: interpweight_get_var4dims_file
    !! 0.4 Local variables
    INTEGER                                              :: nid, vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(4)                                :: dimsid
    CHARACTER(LEN=250)                                   :: msg

    IF (LEN_TRIM(varname) < 1) THEN
      msg = "  interpweight_get_var4dims_file: any variable name '" // TRIM(varname) // "' was provided !!"
      CALL ipslerr_p(3,'interpweight_get_var4dims_file',TRIM(msg), '', '')
    END IF

    IF (is_root_prc) THEN
       ! Open file for read only mode
       rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_open', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inq_varid(nid, varname, vid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inq_varid', &
            TRIM(nf90_strerror(rcode)),'')

       rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_variable 1', &
            TRIM(nf90_strerror(rcode)),'')

       IF (Ndims /= 4) THEN
          msg = "variable '" // TRIM(varname) // "' has not 4 dimensions!!"
          CALL ipslerr(3,'interpweight_get_var4dims_file',TRIM(msg), '', '')
       END IF

       rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_variable 2', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = interpweight_get_var4dims_file(1))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 1', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = interpweight_get_var4dims_file(2))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 2', &
            TRIM(nf90_strerror(rcode)),'')

       rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = interpweight_get_var4dims_file(3))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 3', &
            TRIM(nf90_strerror(rcode)),'')
       
       rcode = nf90_inquire_dimension(nid, dimsid(4), LEN = interpweight_get_var4dims_file(4))
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 4', &
            TRIM(nf90_strerror(rcode)),'')

       rcode = NF90_CLOSE(nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_var4dims_file','Problem with nf90_close', &
            TRIM(nf90_strerror(rcode)),'')
    END iF
    CALL bcast(interpweight_get_var4dims_file)

  END FUNCTION interpweight_get_var4dims_file

!!================================================================================================================================
!! FUNCTION     : interpweight_get_varNdims_file
!!
!>\BRIEF        ! Function to get the rank(number of dimensions) of a given variable inside a file
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_get_varNdims_file
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  INTEGER FUNCTION interpweight_get_varNdims_file(filename, varname)

    USE netcdf

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=*), INTENT(in)                         :: filename      !! filename: name of the file to open
    CHARACTER(LEN=*), INTENT(in)                         :: varname       !! varname: name of the variable 
    !! 0.2 Modified variables
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: nid, vid
    INTEGER                                              :: rcode
    CHARACTER(LEN=256)                                   :: msg


    IF (LEN_TRIM(varname) < 1) THEN
      msg = "  interpweight_get_varNdims_file: any variable name '" // TRIM(varname) // "' was provided !!"
      CALL ipslerr_p(3,'interpweight_get_varNdims_file',TRIM(msg), '', '')
    END IF

    IF (is_root_prc) THEN
       
       ! Open file for read only mode
       rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_varNdims_file','Problem with nf90_open', &
            TRIM(nf90_strerror(rcode)), '')
       
       rcode = nf90_inq_varid(nid, varname, vid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_varNdims_file','Problem with nf90_inq_varid', &
            TRIM(nf90_strerror(rcode)), '')
       
       rcode = nf90_inquire_variable(nid, vid, NDIMS = interpweight_get_varNdims_file)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_varNdims_file','Problem with nf90_inquire_variable', &
            TRIM(nf90_strerror(rcode)), '')

       rcode = NF90_CLOSE(nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'interpweight_get_varNdims_file','Problem with nf90_close', &
            TRIM(nf90_strerror(rcode)), '')
    END IF
    CALL bcast(interpweight_get_varNdims_file)

  END FUNCTION interpweight_get_varNdims_file

!!================================================================================================================================
!! FUNCTION     : isin_file
!!
!>\BRIEF        ! Function to tell if a given variable is inside a file
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : isin_file
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
LOGICAL FUNCTION isin_file(filename, varname)

    USE netcdf

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    CHARACTER(LEN=*), INTENT(in)                         :: filename      !! filename: name of the file to open
    CHARACTER(LEN=*), INTENT(in)                         :: varname       !! varname: name of the variable 
    !! 0.2 Modified variables
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: nid, vid, Ndims, Nvars
    INTEGER                                              :: iv, rcode
    CHARACTER(LEN=1000)                                  :: varinfile
    CHARACTER(LEN=250)                                   :: msg


    IF (is_root_prc) THEN
       rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'isin_file','Problem with nf90_open',TRIM(nf90_strerror(rcode)),'')

       rcode = nf90_inquire(nid, Ndims, Nvars)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'isin_file','Problem with nf90_inquire',TRIM(nf90_strerror(rcode)),'')

       DO iv=1, Nvars
          rcode = nf90_inquire_variable(nid, iv, name=varinfile)
          IF (rcode /= NF90_NOERR) CALL ipslerr(3,'isin_file','Problem with nf90_inquire_variable',TRIM(nf90_strerror(rcode)),'')
          
          IF (TRIM(varinfile) == TRIM(varname)) THEN
             isin_file = .TRUE.
             EXIT
          ELSE
             isin_file = .FALSE.
          END IF
       END DO

       rcode = NF90_CLOSE(nid)
       IF (rcode /= NF90_NOERR) CALL ipslerr(3,'isin_file','Problem with nf90_close',TRIM(nf90_strerror(rcode)),'')
    END IF

    CALL bcast(isin_file)

  END FUNCTION isin_file

!!!!!!! !!!!!! !!!!! !!!! !!! !! !
! Generic functions no netCDF

!!================================================================================================================================
!! FUNCTION     : interpweight_Index1DArrayI
!!
!>\BRIEF        ! Function to provide the first index of a given value inside a 1D intger array
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_Index1DArrayI
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  INTEGER FUNCTION interpweight_Index1DArrayI(array1D, d1, val)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! size of the array
    INTEGER, INTENT(in)                                  :: val           !! value to search
    INTEGER, DIMENSION(d1), INTENT(in)                   :: array1D       !! values of the array
    !! 0.2 Modified variables
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i


    interpweight_Index1DArrayI = -1

    DO i=1,d1
      IF (array1d(i) == val) THEN
        interpweight_Index1DArrayI = i
        EXIT
      END IF
    END DO

  END FUNCTION interpweight_Index1DArrayI

!!================================================================================================================================
!! FUNCTION     : interpweight_Index1DArrayR
!!
!>\BRIEF        ! Function to provide the first index of a given value inside a 1D real array
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_Index1DArrayR
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  INTEGER FUNCTION interpweight_Index1DArrayR(array1D, d1, val)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! size of the array 
    REAL, INTENT(in)                                     :: val           !! value to search
    REAL, DIMENSION(d1), INTENT(in)                      :: array1D       !! values of the array
    !! 0.2 Modified variables
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i

    interpweight_Index1DArrayR = -1

    DO i=1,d1
      IF (array1d(i) == val) THEN
        interpweight_Index1DArrayR = i
        EXIT
      END IF
    END DO

  END FUNCTION interpweight_Index1DArrayR

!!================================================================================================================================
!! FUNCTION     : interpweight_Index2DArrayI
!!
!>\BRIEF        ! Function to provide the first index of a given value inside a 2D integer array
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_Index2DArrayI
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_Index2DArrayI(array2D, d1, d2, val)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1, d2        !! size of the array 
    INTEGER, INTENT(in)                                  :: val           !! value to search
    INTEGER, DIMENSION(d1,d2), INTENT(in)                :: array2D       !! values of the array
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(2)                                :: interpweight_Index2DArrayI
    !! 0.4 Local variables
    INTEGER                                              :: i, j


    interpweight_Index2DArrayI = -1

    DO i=1,d1
      DO j=1,d2
        IF (array2d(i,j) == val) THEN
          interpweight_Index2DArrayI(1) = i
          interpweight_Index2DArrayI(2) = j
          EXIT
        END IF
      END DO
    END DO

  END FUNCTION interpweight_Index2DArrayI

!!================================================================================================================================
!! FUNCTION     : interpweight_Index2DArrayR
!!
!>\BRIEF        ! Function to provide the first index of a given value inside a 2D real array
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : bm_vol
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_Index2DArrayR(array2D, d1, d2, val)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1, d2        !! size of the array 
    REAL, INTENT(in)                                     :: val           !! value to search
    REAL, DIMENSION(d1,d2), INTENT(in)                   :: array2D       !! values of the array
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(2)                                :: interpweight_Index2DArrayR
    !! 0.4 Local variables
    INTEGER                                              :: i, j

    interpweight_Index2DArrayR = -1

    DO i=1,d1
      DO j=1,d2
        IF (array2d(i,j) == val) THEN
          interpweight_Index2DArrayR(1) = i
          interpweight_Index2DArrayR(2) = j
          EXIT
        END IF
      END DO
    END DO

  END FUNCTION interpweight_Index2DArrayR

!!================================================================================================================================
!! FUNCTION     : interpweight_Index1DLonLat
!!
!>\BRIEF        ! Function to provide the first index of a given pair of 
!! longitude and latitude from a set of lon, lat 1D arrays
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_Index1DLonLat
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  INTEGER FUNCTION interpweight_Index1DLonLat(lon, lat, d1, lonval, latval)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! size of the array 
    REAL(r_std), INTENT(in)                              :: lonval,latval !! longitude and latitude value
    REAL(r_std), DIMENSION(d1), INTENT(in)               :: lon, lat      !! longitude and latitudes to search in
    !! 0.2 Modified variables
    !! 0.3 Output variables
    !! 0.4 Local variables
    INTEGER                                              :: i
    REAL(r_std)                                          :: mindist
    REAL(r_std), DIMENSION(d1)                           :: dist


    interpweight_Index1DLonLat = -1

    dist = SQRT((lon - lonval)**2. + (lat - latval)**2.)
!    mindist = MINVAL(ABS(dist))
    mindist = 0.

    DO i=1,d1
      IF (dist(i) == mindist) THEN
        interpweight_Index1DLonLat = i
        EXIT
      END IF
    END DO

  END FUNCTION interpweight_Index1DLonLat

!!================================================================================================================================
!! FUNCTION     : interpweight_Index2DLonLat
!!
!>\BRIEF          Function to provide the first index of a given pair of longitude and latitude
!!                from a set of lon, lat 2D arrays
!!
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_Index2DLonLat
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_Index2DLonLat(lon, lat, d1, d2, lonval, latval)

    USE constantes_var

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1, d2        !! size of the array 
    REAL(r_std), INTENT(in)                              :: lonval,latval !! longitude and latitude value
    REAL(r_std), DIMENSION(d1,d2), INTENT(in)            :: lon, lat      !! longitude and latitudes to search in
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(2)                                :: interpweight_Index2DLonLat
    !! 0.4 Local variables
    INTEGER                                              :: i, j
    REAL(r_std)                                          :: mindist
    REAL(r_std), DIMENSION(d1,d2)                        :: dist


    interpweight_Index2DLonLat = -1

    dist = SQRT((lon - lonval)**2. + (lat - latval)**2.)
    mindist = 0.

    DO i=1,d1
      DO j=1,d2
        IF (dist(i,j) == mindist) THEN
          interpweight_Index2DLonLat(1) = i
          interpweight_Index2DLonLat(2) = j
          EXIT
        END IF
      END DO
    END DO

  END FUNCTION interpweight_Index2DLonLat

!!================================================================================================================================
!! FUNCTION     : interpweight_RangeI
!!
!>\BRIEF        ! Function to provide a range of d1 values from 'iniv' to 
!!  'endv', of integer values in a vector
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_RangeI
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_RangeI(d1, iniv, endv)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! number of values
    INTEGER, INTENT(in)                                  :: iniv          !! initial value
    INTEGER, INTENT(in)                                  :: endv          !! end value
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(d1)                               :: interpweight_RangeI
    !! 0.4 Local variables
    INTEGER                                              :: i, intv

    intv = (endv - iniv) / (d1*1 - 1)

    interpweight_RangeI(1) = iniv
    DO i=2,d1
      interpweight_RangeI(i) = interpweight_RangeI(i-1) + intv
    END DO

  END FUNCTION interpweight_RangeI 

!!================================================================================================================================
!! FUNCTION     : interpweight_RangeR
!!
!>\BRIEF        ! Function to provide a range of d1 values from 'iniv' to 
!!  'endv', of real values in a vector
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_RangeR
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_RangeR(d1, iniv, endv)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! number of values
    REAL, INTENT(in)                                     :: iniv          !! initial value
    REAL, INTENT(in)                                     :: endv          !! end value
    !! 0.2 Modified variables
    !! 0.3 Output variables
    REAL, DIMENSION(d1)                                  :: interpweight_RangeR
    !! 0.4 Local variables
    INTEGER                                              :: i
    REAL                                                 :: intv

    intv = (endv - iniv) / (d1*1. - 1.)

    interpweight_RangeR(1) = iniv
    DO i=2,d1
      interpweight_RangeR(i) = interpweight_RangeR(i-1) + intv
    END DO

  END FUNCTION interpweight_RangeR

!!================================================================================================================================
!! FUNCTION     : interpweight_ValVecI
!!
!>\BRIEF        ! Function to provide the number of times and where that a 
!!  given value 'oper' on a vector of integers
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_ValVecI
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_ValVecI(vec, d1, val, oper)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! Number of values
    INTEGER, DIMENSION(d1), INTENT(in)                   :: vec           !! vector of integer values
    INTEGER, INTENT(in)                                  :: val           !! value to search
    CHARACTER(LEN=*), INTENT(in)                         :: oper          !! operation:
                                                                          !!   'eq': equal
                                                                          !!   'ge': greater or equal
                                                                          !!   'le': less or equal
                                                                          !!   'ne': not equal
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(d1)                               :: interpweight_ValVecI !! Vector with positions 
                                                                                 !! (NtimesVal, pos1, pos2, ..., posNtimesVal, ...)
                                                                                 !! NtimesVal=number of correspondencies with 'oper'
    !! 0.4 Local variables
    INTEGER                                              :: i, NtimesVal
    CHARACTER(LEN=50)                                    :: errormsg

    errormsg = 'ERROR -- error -- ERROR -- error'

    NtimesVal = 0
    DO i=1,d1
      SELECT CASE(oper)
        CASE ('eq')
          IF (vec(i) == val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
          END IF
        CASE ('ge')
          IF (vec(i) >= val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
          END IF
        CASE ('le')
          IF (vec(i) <= val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
          END IF
        CASE ('neq')
          IF (vec(i) /= val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
          END IF
        CASE DEFAULT
          WRITE(numout,*)TRIM(errormsg)
          WRITE(numout,*)"  interpweight_ValVecI: operation '" // TRIM(oper) // "' not ready!!"
          WRITE(numout,*)"    only available: 'eq', 'ge', 'le', 'neq'"
          CALL ipslerr_p (3, 'interpweight_ValVecI', 'wrong operation', 'provide another one', '')
      END SELECT
    END DO
    interpweight_ValVecI(1) = NtimesVal


  END FUNCTION interpweight_ValVecI

!!================================================================================================================================
!! FUNCTION     : interpweight_ValVecR
!!
!>\BRIEF         Function to provide the number of times and where that a 
!!  given value 'oper' on a vector of reals
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_ValVecR
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_ValVecR(vec, d1, val, oper)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: d1            !! Number of values
    REAL, DIMENSION(:), INTENT(in)                       :: vec           !! vector of integer values
    REAL, INTENT(in)                                     :: val           !! value to search
    CHARACTER(LEN=*), INTENT(in)                         :: oper          !! operation:
                                                                          !!   'eq': equal
                                                                          !!   'ge': greater or equal
                                                                          !!   'le': less or equal
                                                                          !!   'neq': not equal
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(d1)                               :: interpweight_ValVecR !! Vector with positions 
                                                                                 !! (NtimesVal, pos1, pos2, ..., posNtimesVal, ...)
                                                                                 !! NtimesVal=number of correspondencies with 'oper'
    !! 0.4 Local variables
    INTEGER                                              :: i, NtimesVal
    CHARACTER(LEN=50)                                    :: errormsg

    errormsg = 'ERROR -- error -- ERROR -- error'

    NtimesVal = 0
    DO i=1,d1
      SELECT CASE(oper)
        CASE ('eq')
          IF (vec(i) == val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
          END IF
        CASE ('ge')
          IF (vec(i) >= val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
          END IF
        CASE ('le')
          IF (vec(i) <= val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
          END IF
        CASE ('neq')
          IF (vec(i) /= val) THEN
            NtimesVal = NtimesVal + 1
            IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
          END IF
        CASE DEFAULT
          WRITE(numout,*)TRIM(errormsg)
          WRITE(numout,*)"  interpweight_ValVecR: operation '" // TRIM(oper) // "' not ready!!"
          WRITE(numout,*)"    only available: 'eq', 'ge', 'le', 'neq'"
          CALL ipslerr_p (3, 'interpweight_ValVecR', 'wrong operation', 'provide another one', '')
      END SELECT
    END DO
    interpweight_ValVecR(1) = NtimesVal


  END FUNCTION interpweight_ValVecR

!!================================================================================================================================
!! FUNCTION     : interpweight_From2DmatTo1Dvec
!!
!>\BRIEF        ! Suroutine to provide value on a transformation from a 2D 
!!  matrix to a 1D vector
!!    e.g.: from T2(i,j) to T2M(ijklon)
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_From2DmatTo1Dvec
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  INTEGER FUNCTION interpweight_From2DmatTo1Dvec(ix,iy,dx,dy)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: ix,iy         !! location within the 2D matrix
    INTEGER, INTENT(in)                                  :: dx,dy         !! dimension of the 2D matrix

    IF (ix > dx) THEN
      WRITE(numout,*) "Error: "
      WRITE(numout,*)"  interpweight_From2DmatTo1Dvec: 'ix' too large ix > dx ", ix, '>', dx, '!!'
      CALL ipslerr_p (3, 'interpweight_From2DmatTo1Dvec', 'ix value too big', 'out of bounds', '')
    END IF

    IF (iy > dy) THEN
      WRITE(numout,*) "Error: "
      WRITE(numout,*)"  interpweight_From2DmatTo1Dvec: 'iy' too large iy > dy ", iy, '>', dy, '!!'
      CALL ipslerr_p (3, 'interpweight_From2DmatTo1Dvec', 'iy value too big', 'out of bounds', '')
    END IF

    interpweight_From2DmatTo1Dvec = (ix-1)*dy + iy

  END FUNCTION interpweight_From2DmatTo1Dvec

!!================================================================================================================================
!! FUNCTION     : interpweight_From1DvecTo2Dmat
!!
!>\BRIEF        ! Suroutine to provide value on a transformation from a 1D 
!!  vector to a 2D matrix 
!!     e.g.: from T2M(ijklon) to T2(i,j)
!!
!! DESCRIPTION : None
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : interpweight_From1DvecTo2Dmat
!!
!! REFERENCE(S) : See above, module description.
!!
!! FLOWCHART    : None
!! \n
!_================================================================================================================================
  FUNCTION interpweight_From1DvecTo2Dmat(ivec,dx,dy)

    IMPLICIT NONE

    !! 0. Variables and parameter declaration
    !! 0.1 Input variables
    INTEGER, INTENT(in)                                  :: ivec              !! position within the 1D vec
    INTEGER, INTENT(in)                                  :: dx,dy             !! dimension of the 2D matrix
    !! 0.2 Modified variables
    !! 0.3 Output variables
    INTEGER, DIMENSION(2)                                :: interpweight_From1DvecTo2Dmat


    IF (ivec > dx*dy) THEN
      WRITE(numout,*) "Error: "
      WRITE(numout,*)"  interpweight_From1DvecTo2Dmat: 'ivec' too large ivec > dx*dy ", ivec, '>', dx*dy, '!!'
      CALL ipslerr_p (3, 'interpweight_From1DvecTo2Dmat', "'ivec' too large", 'out of bounds!!', '')
    END IF

    interpweight_From1DvecTo2Dmat(1) = INT((ivec-1)/dy) + 1
    interpweight_From1DvecTo2Dmat(2) = ivec - (interpweight_From1DvecTo2Dmat(1)-1)*dy

  END FUNCTION interpweight_From1DvecTo2Dmat

END MODULE interpweight
