MODULE module_generic
! Module with generic functions

!!!!!!! Subroutines/Functions
! continguos_homogene_zones: Subroutine to look for contiguous zones by looking by continuous grid points
! freeunit: provides the number of a free unit in which open a file
! fill_matrix2DRK_winmat2D_list1D: Subroutine to fill a 2D RK matrix using a list of 1D indices from 
!   another given 2D matrix
! from_coordlist_2DRKmatrix: Subroutine to construct a 2D RK matrix from a list of values accompaigned 
!   by a list of coordinates to find i,j grid-point coordinates by minimum distance
! from_ptlist_2DRKmatrix: Subroutine to construct a 2D RK matrix from a list of values accompaigned 
!   by a list of grid-point coordinates
! from_ptlist_2DRKNmatrix: Subroutine to construct N 2D RK matrix from a list of values accompaigned 
!   by a list of grid-point coordinates
! GetInNamelist: Subroutine to get a paramter from a namelistfile
! GDATE: Subroutine to compute the gregorian calendar date (year,month,day) given the julian date (JD)
! get_xyconlimits: Subroutine for getting the limits of contiguous values from a given point in a 2D matrix
! index_list_coordsI: Function to provide the index of a given coordinate within a list of integer coordinates
! Index1DArrayI: Function to provide the first index of a given value inside a 1D integer array
! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
! Index1DArrayR_K: Function to provide the first index of a given value inside a 1D real(r_k) array
! Index1DArrayL: Function to provide the first index of a given value inside a 1D boolean array
! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array
! Index2DArrayR_K: Function to provide the first index of a given value inside a 2D real(r_k) array
! index_samevals1D_RK: Subroutine to search for the indices of the same values between 2 1D RK series 
!   of values allowing repetitions
! index_samevals2D_RK: Subroutine to search for the indices of the same values between 2 2D RK series 
!   of values allowing repetitions
! JD: Fucntion to compute the julian date (JD) given a gregorian calendar
! Nvalues_2DArrayI: Number of different values of a 2D integer array
! mat2DPosition: Function to provide the i, j indices of a given value inside a 2D matrix
! multi_Index1DArrayL: Subroutine to provide the indices of a given value inside a 1D boolean array
! numberTimes: Function to provide the number of times that a given set of characters happen within a string
! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
! RangeR: Function to provide a range of d1 values from 'iniv' to 'endv', of real values in a vector
! RangeR_K: Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector
! rm_values_vecRK: Subroutine to remove a given value (e.g. fill_Value) from a RK vector
! stoprun: Subroutine to stop running and print a message
! vectorI_S: Function to transform a vector of integers to a string of characters
! vectorR_S: Function to transform a vector of reals to a string of characters
! xzones_homogenization: Subroutine to homogenize 2D contiguous zones along x-axis. Zones might be 
!   contiguous, but with different number assigned !
! yzones_homogenization: Subroutine to homogenize 2D contiguous zones along y-axis. Zones might be 
!   contiguous, but with different number assigned !

  USE module_definitions
  USE module_basic

  CONTAINS

  SUBROUTINE Nvalues_2DArrayI(dx, dy, dxy, mat2DI, Nvals, vals)
! Subroutine to give the number of different values of a 2D integer array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dx, dy, dxy
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: mat2DI
    INTEGER, INTENT(out)                                 :: Nvals
    INTEGER, DIMENSION(dxy), INTENT(out)                 :: vals

! Local
    INTEGER                                              :: i, j, ij

!!!!!!! Variables
! dx, dy: size of the 2D space
! mat2DI: 2D integer matrix
! Nvals: number of different values
! vals: vector with the different values

  fname = 'Nvalues_2DArrayI'

  vals = 0

  Nvals = 1
  vals(1) = mat2DI(1,1)  
  DO i=1,dx
    DO j=1,dy
      IF (Index1DArrayI(vals, Nvals, mat2DI(i,j)) == -1) THEN
        Nvals = Nvals + 1
        vals(Nvals) = mat2DI(i,j)
      END IF
    END DO
  END DO

  RETURN

  END SUBROUTINE Nvalues_2DArrayI

  INTEGER FUNCTION index_list_coordsI(Ncoords, coords, icoord)
  ! Function to provide the index of a given coordinate within a list of integer coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Ncoords
    INTEGER, DIMENSION(Ncoords,2), INTENT(in)            :: coords
    INTEGER, DIMENSION(2), INTENT(in)                    :: icoord

! Local
    INTEGER, DIMENSION(Ncoords)                          :: dist
    INTEGER                                              :: i,mindist
    INTEGER, DIMENSION(1)                                :: iloc

!!!!!!! Variables
! Ncoords: number of coordinates in the list
! coords: list of coordinates
! icoord: coordinate to find

  fname = 'index_list_coordsI'

  dist = (coords(:,1)-icoord(1))**2+(coords(:,2)-icoord(2))**2

  IF (ANY(dist == 0)) THEN
    iloc = MINLOC(dist)
    index_list_coordsI = iloc(1)
  ELSE
    index_list_coordsI = -1
  END IF

  END FUNCTION index_list_coordsI

  INTEGER FUNCTION Index1DArrayI(array1D, d1, val)
! Function to provide the first index of a given value inside a 1D integer array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    INTEGER, INTENT(in)                                  :: val
    INTEGER, DIMENSION(d1), INTENT(in)                   :: array1D

! Local
    INTEGER                                              :: i

    fname = 'Index1DArrayI'

    Index1DArrayI = -1

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

  END FUNCTION Index1DArrayI

  INTEGER FUNCTION Index1DArrayR(array1D, d1, val)
! Function to provide the first index of a given value inside a 1D real array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL, INTENT(in)                                     :: val
    REAL, DIMENSION(d1), INTENT(in)                      :: array1D

! Local
    INTEGER                                              :: i

    fname = 'Index1DArrayR'

    Index1DArrayR = -1

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

  END FUNCTION Index1DArrayR

  INTEGER FUNCTION Index1DArrayR_K(array1D, d1, val)
! Function to provide the first index of a given value inside a 1D real(r_k) array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL(r_k), INTENT(in)                                :: val
    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: array1D

! Local
    INTEGER                                              :: i

    fname = 'Index1DArrayR_K'

    Index1DArrayR_K = -1

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

  END FUNCTION Index1DArrayR_K

  INTEGER FUNCTION Index1DArrayL(array1D, d1, val)
! Function to provide the first index of a given value inside a 1D boolean array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    LOGICAL, INTENT(in)                                  :: val
    LOGICAL, DIMENSION(d1), INTENT(in)                   :: array1D

! Local
    INTEGER                                              :: i
    CHARACTER(LEN=50)                                    :: fname

    fname = 'Index1DArrayL'

    Index1DArrayL = -1

    DO i=1,d1
      IF (array1d(i) .EQV. val) THEN
        Index1DArrayL = i
        EXIT
      END IF
    END DO

  END FUNCTION Index1DArrayL

  SUBROUTINE multi_Index1DArrayL(array1D, d1, val, Ntimes, pos)
! Subroutine to provide the indices of a given value inside a 1D boolean array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    LOGICAL, INTENT(in)                                  :: val
    LOGICAL, DIMENSION(d1), INTENT(in)                   :: array1D
    INTEGER, INTENT(out)                                 :: Ntimes
    INTEGER, DIMENSION(d1), INTENT(out)                  :: pos

! Local
    INTEGER                                              :: i
    CHARACTER(LEN=50)                                    :: fname

!!!!!!! Variables
! array1D: 1D Array of values
! d1: length of array
! val: Value to look for
! Ntimes: Number of times val is found within array1D
! pos: positions of the values of val

    fname = 'multi_Index1DArrayL'

    Ntimes = 0
    pos = -1

    DO i=1,d1
      IF (array1d(i) .EQV. val) THEN
        Ntimes = Ntimes + 1
        pos(Ntimes) = i
      END IF
    END DO

  END SUBROUTINE multi_Index1DArrayL

  FUNCTION Index2DArrayR(array2D, d1, d2, val)
! Function to provide the first index of a given value inside a 2D real array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2
    REAL, INTENT(in)                                     :: val
    REAL, DIMENSION(d1,d2), INTENT(in)                   :: array2D
    INTEGER, DIMENSION(2)                                :: Index2DArrayR

! Local
    INTEGER                                              :: i, j

    fname = 'Index2DArrayR'

    Index2DArrayR = -1

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

  END FUNCTION Index2DArrayR

  FUNCTION Index2DArrayR_K(array2D, d1, d2, val)
! Function to provide the first index of a given value inside a 2D real array

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2
    REAL(r_k), INTENT(in)                                :: val
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: array2D
    INTEGER, DIMENSION(2)                                :: Index2DArrayR_K

! Local
    INTEGER                                              :: i, j

    fname = 'Index2DArrayR_K'

    Index2DArrayR_K = -1

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

  END FUNCTION Index2DArrayR_K

  INTEGER FUNCTION numberTimes(String, charv)
! Function to provide the number of times that a given set of characters happen within a string

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(IN)                         :: String, charv

! Local
    INTEGER                                              :: i, Lstring, Lcharv

    numberTimes = 0

    Lstring = LEN_TRIM(String)
    Lcharv = LEN_TRIM(charv)

    DO i=1,Lstring - Lcharv
      IF (String(i:i+Lcharv-1) == TRIM(charv)) numberTimes = numberTimes + 1
    END DO

  END FUNCTION numberTimes


  FUNCTION RangeI(d1, iniv, endv)
! Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, iniv, endv
    INTEGER, DIMENSION(d1)                               :: RangeI

! Local
    INTEGER                                              :: i, intv

    fname = 'RangeI'

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

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

  END FUNCTION RangeI 

  FUNCTION RangeR(d1, iniv, endv)
! Function to provide a range of d1 from 'iniv' to 'endv', of real values in a vector

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL, INTENT(in)                                     :: iniv, endv
    REAL, DIMENSION(d1)                                  :: RangeR

! Local
    INTEGER                                              :: i
    REAL                                                 :: intv

    fname = 'RangeR'

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

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

  END FUNCTION RangeR

  FUNCTION RangeR_K(d1, iniv, endv)
! Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL(r_k), INTENT(in)                                :: iniv, endv
    REAL(r_k), DIMENSION(d1)                             :: RangeR_K

! Local
    INTEGER                                              :: i
    REAL(r_k)                                            :: intv

    fname = 'RangeR_K'

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

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

  END FUNCTION RangeR_K

  INTEGER FUNCTION freeunit()
! provides the number of a free unit in which open a file

    IMPLICIT NONE

    LOGICAL                                              :: is_used

    is_used = .true.
    DO freeunit=10,100
      INQUIRE(unit=freeunit, opened=is_used)
      IF (.not. is_used) EXIT
    END DO

    RETURN

  END FUNCTION freeunit

  SUBROUTINE GetInNamelist(namelistfile, param, kindparam, Ival, Rval, Lval, Sval)
! Subroutine to get a paramter from a namelistfile

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(IN)                         :: namelistfile, param
    CHARACTER(LEN=1), INTENT(IN)                         :: kindparam
    INTEGER, OPTIONAL, INTENT(OUT)                       :: Ival
    REAL, OPTIONAL, INTENT(OUT)                          :: Rval
    LOGICAL, OPTIONAL, INTENT(OUT)                       :: Lval
    CHARACTER(LEN=200), OPTIONAL, INTENT(OUT)            :: Sval

! Local 
    INTEGER                                              :: i, funit, ios
    INTEGER                                              :: Lparam, posparam
    LOGICAL                                              :: is_used
    CHARACTER(LEN=1000)                                  :: line, message
    CHARACTER(LEN=200), DIMENSION(2)                     :: lvals
    CHARACTER(LEN=200)                                   :: pval

!!!!!!! Variables
! namelistfile: name of the namelist file
! param: parameter to get
! paramkind: kind of the parameter (I: Integer, L: boolean, R: Real, S: String)

    fname = 'GetInNamelist'

! Reading dimensions file and defining dimensions
    is_used = .true.
    DO funit=10,100
      INQUIRE(unit=funit, opened=is_used)
      IF (.not. is_used) EXIT
    END DO

    OPEN(funit, FILE=TRIM(namelistfile), STATUS='old', FORM='formatted', IOSTAT=ios)
    IF ( ios /= 0 ) CALL stoprun(message, fname)

    Lparam = LEN_TRIM(param)
  
    DO i=1,10000
      READ(funit,"(A200)",END=100)line
      posparam = INDEX(TRIM(line), TRIM(param))
      IF (posparam /= 0) EXIT

    END DO
 100 CONTINUE

    IF (posparam == 0) THEN
      message = "namelist '" // TRIM(namelistfile) // "' does not have parameter '" // TRIM(param) // &
        "' !!"
      CALL stoprun(message, fname)
    END IF

    CLOSE(UNIT=funit)

    CALL split(line, '=', 2, lvals)
    IF (kindparam /= 'S') THEN
      CALL RemoveNonNum(lvals(2), pval)
    END IF

! L. Fita, LMD. October 2015
!   Up to now, only getting scalar values
    kparam: SELECT CASE (kindparam)
      CASE ('I')
        Ival = StoI(pval)
!        PRINT *,TRIM(param),'= ', Ival
      CASE ('L')
        Lval = StoL(pval)
!        PRINT *,TRIM(param),'= ', Lval
      CASE ('R')
        Rval = StoR(pval)
!        PRINT *,TRIM(param),'= ', Rval
      CASE ('S') 
        Sval = lvals(2)

      CASE DEFAULT
        message = "type of parameter '" // kindparam // "' not ready !!"
        CALL stoprun(message, fname)

    END SELECT kparam

  END SUBROUTINE GetInNamelist

  SUBROUTINE stoprun(msg, fname)
! Subroutine to stop running and print a message

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(IN)                           :: fname
    CHARACTER(LEN=*), INTENT(IN)                           :: msg

! local
    CHARACTER(LEN=50)                                      :: errmsg, warnmsg

    errmsg = 'ERROR -- error -- ERROR -- error'

    PRINT *, TRIM(errmsg)
    PRINT *, '  ' // TRIM(fname) // ': ' // TRIM(msg)
    STOP

  END SUBROUTINE stoprun

  SUBROUTINE xzones_homogenization(dx, dy, inzones, outzones)
! Subroutine to homogenize 2D contiguous zones along x-axis, zones might be contiguous, but with 
!   different number assigned !
!   Here we have a 2D matrix of integers, with contiguous integer filled zones, zero outside any zone
!   It might be, that within the same zone, might be lines which do not share the same integer
!     0 0 0 0 0          0 0 0 0 0 
!     0 1 1 0 0          0 1 1 0 0
!     0 2 0 0 1    == >  0 1 0 0 2
!     0 1 1 0 0          0 1 1 0 0 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dx, dy
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: inzones
    INTEGER, DIMENSION(dx,dy), INTENT(out)               :: outzones

! Local
    INTEGER                                              :: i,j,k
    INTEGER                                              :: Nmaxzones, TOTzones
    LOGICAL                                              :: assigned
    INTEGER, DIMENSION(dx)                               :: prevline
    INTEGER, DIMENSION(dy)                               :: Nxzones
    INTEGER, DIMENSION(:,:,:), ALLOCATABLE               :: zones

!!!!!!! Variables
! dx, dy: Shape of the 2D space
! inzones: zones to homogenize
! outzones: zones homogenized

    fname = 'xzones_homogenization'

    ! Maximum possible number of zones
    Nmaxzones = INT((dx/2)*(dy/2))

    ! Matrix with [i,j,Nzone,izone/ezone]
    IF (ALLOCATED(zones)) DEALLOCATE(zones)
    ALLOCATE(zones(dy,Nmaxzones,3))

    zones = 0
    Nxzones = 0
    ! Getting beginning/end of y-bands
    DO j=1, dy
      k = 0
      i = 1
      IF (inzones(i,j) /= 0) THEN
        k = k + 1
        zones(j,k,1) = i
        zones(j,k,3) = k
      END IF
      DO i=2, dx
        IF ( (inzones(i,j) /= 0) .AND. (inzones(i-1,j) == 0)) THEN
          k = k+1
          zones(j,k,1) = i
          zones(j,k,3) = k
        ELSE IF ( (inzones(i-1,j) /= 0) .AND. (inzones(i,j) == 0)) THEN
          zones(j,k,2) = i-1
          zones(j,k,3) = k
        END IF
      END DO
      IF (k > 0) THEN
        IF (zones(j,k,2) == 0) zones(j,k,2) = dx
      END IF
      Nxzones(j) = k
    END DO

    ! Homogenizing contigous zones
    outzones = 0
    TOTzones = 0
    j = 1
    DO k = 1, Nxzones(j)
      TOTzones = TOTzones + 1
      DO i=zones(j,k,1), zones(j,k,2)
        outzones(i,j) = TOTzones
      END DO
    END DO

    DO j=2, dy
      prevline = outzones(:,j-1)
      DO k = 1, Nxzones(j)
        assigned = .FALSE.
        DO i=zones(j,k,1), zones(j,k,2)
          IF (prevline(i) /= 0) THEN
            outzones(zones(j,k,1):zones(j,k,2),j) = prevline(i)
            assigned = .TRUE.
            EXIT
          END IF
        END DO
        IF (.NOT.assigned) THEN
          TOTzones = TOTzones + 1
          DO i=zones(j,k,1), zones(j,k,2)
            outzones(i,j) = TOTzones
          END DO
        END IF
      END DO
    END DO

    IF (ALLOCATED(zones)) DEALLOCATE(zones)

  END SUBROUTINE xzones_homogenization

  SUBROUTINE yzones_homogenization(dx, dy, inzones, outzones)
! Subroutine to homogenize 2D contiguous zones along y-axis, zones might be contiguous, but with 
!   different number assigned !
!   Here we have a 2D matrix of integers, with contiguous integer filled zones, zero outside any zone
!   It might be, that within the same zone, might be lines which do not share the same integer
!     0 0 0 0 0          0 0 0 0 0 
!     0 1 1 0 0          0 1 1 0 0
!     0 2 0 0 1    == >  0 1 0 0 2
!     0 1 1 0 0          0 1 1 0 0 

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dx, dy
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: inzones
    INTEGER, DIMENSION(dx,dy), INTENT(out)               :: outzones

! Local
    INTEGER                                              :: i,j,k
    INTEGER                                              :: Nmaxzones, TOTzones
    LOGICAL                                              :: assigned
    INTEGER, DIMENSION(dy)                               :: prevline
    INTEGER, DIMENSION(dx)                               :: Nyzones
    INTEGER, DIMENSION(:,:,:), ALLOCATABLE               :: zones

!!!!!!! Variables
! dx, dy: Shape of the 2D space
! inzones: zones to homogenize
! outzones: zones homogenized

    fname = 'yzones_homogenization'

    ! Maximum possible number of zones
    Nmaxzones = INT((dx/2)*(dy/2))

    ! Matrix with [i,j,Nzone,izone/ezone]
    IF (ALLOCATED(zones)) DEALLOCATE(zones)
    ALLOCATE(zones(dx,Nmaxzones,3))

    zones = 0
    Nyzones = 0
    ! Getting beginning/end of y-bands
    DO i=1, dx
      k = 0
      j = 1
      IF (inzones(i,j) /= 0) THEN
        k = k + 1
        zones(i,k,1) = j
        zones(i,k,3) = k
      END IF
      DO j=2, dy
        IF ( (inzones(i,j) /= 0) .AND. (inzones(i,j-1) == 0)) THEN
          k = k+1
          zones(i,k,1) = j
          zones(i,k,3) = k
        ELSE IF ( (inzones(i,j-1) /= 0) .AND. (inzones(i,j) == 0)) THEN
          zones(i,k,2) = j-1
          zones(i,k,3) = k
        END IF
      END DO
      IF (k > 0) THEN
        IF (zones(i,k,2) == 0) zones(i,k,2) = dy
      END IF
      Nyzones(i) = k
    END DO

    ! Homogenizing contigous zones
    outzones = 0
    TOTzones = 0
    i = 1
    DO k = 1, Nyzones(i)
      TOTzones = TOTzones + 1
      DO j=zones(i,k,1), zones(i,k,2)
        outzones(i,j) = TOTzones
      END DO
    END DO

    DO i=2, dx
      prevline = outzones(i-1,:)
      DO k = 1, Nyzones(i)
        assigned = .FALSE.
        DO j=zones(i,k,1), zones(i,k,2)
          IF (prevline(j) /= 0) THEN
            outzones(i,zones(i,k,1):zones(i,k,2)) = prevline(j)
            assigned = .TRUE.
            EXIT
          END IF
        END DO
        IF (.NOT.assigned) THEN
          TOTzones = TOTzones + 1
          DO j=zones(i,k,1), zones(i,k,2)
            outzones(i,j) = TOTzones
          END DO
        END IF
      END DO
    END DO

    IF (ALLOCATED(zones)) DEALLOCATE(zones)

  END SUBROUTINE yzones_homogenization

  CHARACTER(len=1000) FUNCTION vectorI_S(d1, vector)
  ! Function to transform a vector of integers to a string of characters

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    INTEGER, DIMENSION(d1), INTENT(in)                   :: vector

! Local
    INTEGER                                              :: iv
    CHARACTER(len=50)                                    :: IS

!!!!!!! Variables
! d1: length of the vector
! vector: values to transform

    fname = 'vectorI_S'

    vectorI_S = ''
    DO iv=1, d1
      WRITE(IS, '(I50)')vector(iv)
      IF (iv == 1) THEN
        vectorI_S = TRIM(IS)
      ELSE
        vectorI_S = TRIM(vectorI_S) // ', ' // TRIM(IS)
      END IF
    END DO    

  END FUNCTION vectorI_S

  CHARACTER(len=1000) FUNCTION vectorR_S(d1, vector)
  ! Function to transform a vector of reals to a string of characters

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL, DIMENSION(d1), INTENT(in)                      :: vector

! Local
    INTEGER                                              :: iv
    CHARACTER(len=50)                                    :: RS

!!!!!!! Variables
! d1: length of the vector
! vector: values to transform

    fname = 'vectorR_S'

    vectorR_S = ''
    DO iv=1, d1
      WRITE(RS, '(F50.25)')vector(iv)
      IF (iv == 1) THEN
        vectorR_S = TRIM(RS)
      ELSE
        vectorR_S = TRIM(vectorR_S) // ', ' // TRIM(RS)
      END IF
    END DO    

  END FUNCTION vectorR_S

  SUBROUTINE continguos_homogene_zones(dx, dy, matvals, Nzones, contzones)
  ! Subroutine to look for contiguous zones by looking by continuous grid points

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dx, dy
    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: matvals
    INTEGER, INTENT(out)                                 :: Nzones
    INTEGER, DIMENSION(dx,dy), INTENT(out)               :: contzones
! Local
    INTEGER                                              :: i,j, k
    INTEGER                                              :: ii, ei, ij, ej
    INTEGER                                              :: Ncont, Nassigned, Ncontmin
    LOGICAL, DIMENSION(dx,dy)                            :: assigned, notdone
    INTEGER, DIMENSION(:), ALLOCATABLE                   :: pzones, allzones
    LOGICAL, DIMENSION(:), ALLOCATABLE                   :: passigns

!!!!!!! Variables
! dx, dy: shape of the matrix
! matvals: matrix with the values
! contzones: homogeneous zones found

    fname = 'continguos_homogene_zones'

    ! Vector to keep track of all zone values
    IF (ALLOCATED(allzones)) DEALLOCATE(allzones)
    ALLOCATE(allzones(dx*dy/4))
    allzones = 0

    assigned = .FALSE.
    notdone = .TRUE.
    contzones = -1
    Nzones = 0
    DO i=1, dx
      DO j=1, dy
        ! First
        IF (matvals(i,j) /= 0 .AND. notdone(i,j)) THEN
          CALL get_xyconlimits(dx, dy, matvals, 0, i, j, ii, ei, ij, ej)

          ! Has any point of the rays already been assigned?
          ! Along x
          IF (ALLOCATED(pzones)) DEALLOCATE(pzones)
          ALLOCATE(pzones(ei-ii+1))
          IF (ALLOCATED(passigns)) DEALLOCATE(passigns)
          ALLOCATE(passigns(ei-ii+1))
          passigns = assigned(ii:ei,j)
          CALL multi_Index1DArrayL(passigns, ei-ii+1, .TRUE., Nassigned, pzones)
          IF (Nassigned /= 0) THEN
            Ncontmin = 10000000
            DO k=1, Nassigned
              IF (contzones(ii+pzones(k)-1,j) < Ncontmin) Ncontmin = contzones(ii+pzones(k)-1,j)
            END DO
            ! If there is more than one assigned value change all values to the minimum one
            DO k=1, Nassigned
              IF (contzones(ii+pzones(k)-1,j) /= Ncontmin) THEN
                WHERE (contzones == contzones(ii+pzones(k)-1,j) .AND. contzones /= Ncontmin)
                  contzones = Ncontmin
                END WHERE
                allzones(contzones(ii+pzones(k)-1,j)) = 0
              END IF
            END DO
            Ncont = Ncontmin
          END IF
          ! Along y
          IF (ALLOCATED(pzones)) DEALLOCATE(pzones)
          ALLOCATE(pzones(ej-ij+1))
          IF (ALLOCATED(passigns)) DEALLOCATE(passigns)
          ALLOCATE(passigns(ej-ij+1))
          passigns = assigned(i,ij:ej)
          CALL multi_Index1DArrayL(passigns, ej-ij+1, .TRUE., Nassigned, pzones)
          IF (Nassigned /= 0) THEN
            Ncontmin = 10000000
            DO k=1, Nassigned
              IF (contzones(i,ij+pzones(k)-1) < Ncontmin) Ncontmin = contzones(i,ij+pzones(k)-1)
            END DO
            ! If there is more than one assigned value change all values to the minimum one
            DO k=1, Nassigned
              IF (contzones(i,ij+pzones(k)-1) /= Ncontmin) THEN
                WHERE (contzones == contzones(i,ij+pzones(k)-1) .AND. contzones /= Ncontmin)
                  contzones = Ncontmin
                END WHERE
                allzones(contzones(i,ij+pzones(k)-1)) = 0
              END IF
            END DO
            Ncont = Ncontmin
          END IF

          IF (.NOT.assigned(i,j)) THEN
            Nzones = Nzones + 1
            Ncont = Nzones
            allzones(Nzones) = 1
          ELSE
            Ncont = contzones(i,j)
          END IF
          contzones(i,j) = Ncont
          contzones(ii:ei,j) = Ncont
          contzones(i,ij:ej) = Ncont
          notdone(i,j) = .FALSE.
          assigned(i,j) = .TRUE.
          assigned(ii:ei,j) = .TRUE.
          assigned(i,ij:ej) = .TRUE.
        END IF
      END DO
    END DO

    ! Using allzones to provide continuous assigned values
    Nzones = 0
    DO k=1, dx*dy/4
      IF (allzones(k) /= 0) THEN
        Nzones = Nzones + 1
        IF (allzones(k) /= Nzones) THEN
          WHERE(contzones == allzones(k))
            contzones = Nzones
          END WHERE
        END IF
      END IF
    END DO

    RETURN

  END SUBROUTINE continguos_homogene_zones

  SUBROUTINE get_xyconlimits(d1, d2, matv, NOval, i, j, ix, ex, iy, ey)
  ! Subroutine for getting the limits of contiguous values from a given point in a 2D matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, i, j, NOval
    INTEGER, DIMENSION(d1,d2), INTENT(in)                :: matv
    INTEGER, INTENT(out)                                 :: ix, ex, iy, ey

    ! Local
    INTEGER                                              :: i1, j1
    LOGICAL                                              :: found

!!!!!!! Variables
! d1, d2: Shape of input 2D values
! i, j: grid point from which get the contiguous values
! NOval: value for grid points without data
! matv: 2D matrx values
! ix, ex, iy, ey: limits from a grid point

    fname = 'get_xyconlimits'

    ix = -1
    ex = -1
    iy = -1
    ey = -1

    ! Before i
    IF (i > 1) THEN
      ix = i
      found = .FALSE.
      DO i1=i-1,1,-1
        IF (matv(i1,j) == NOval) THEN
          ix = i1 + 1
          found = .TRUE.
          EXIT
        END IF
      END DO
      IF (.NOT.found) ix=1
    ELSE
      ix = 1
    END IF
    ! After i
    IF (i < d1) THEN
      ex = i
      found = .FALSE.
      DO i1=i+1,d1
        IF (matv(i1,j) == NOval) THEN
          ex = i1 - 1
          found = .TRUE.
          EXIT
        END IF
      END DO
      IF (.NOT.found) ex=d1
    ELSE
      ex = d1
    END IF

    ! Before j
    IF (j > 1) THEN
      iy = j
      found = .FALSE.
      DO j1=j-1,1,-1
        IF (matv(i,j1) == NOval) THEN
          iy = j1 + 1
          found = .TRUE.
          EXIT
        END IF
      END DO
      IF (.NOT.found) iy=1
    ELSE
      iy = 1
    END IF
    ! After j
    IF (j < d2) THEN
      ey = j
      found = .FALSE.
      DO j1=j+1,d2
        IF (matv(i,j1) == NOval) THEN
          ey = j1 - 1
          found = .TRUE.
          EXIT
        END IF
      END DO
      IF (.NOT.found) ey=d2
    ELSE
      ey = d2
    END IF

    RETURN

  END SUBROUTINE get_xyconlimits

  SUBROUTINE from_ptlist_2DRKmatrix(Npts, pts, Flike, vals, dx, dy, NOval, matrix)
  ! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of grid-point 
  !   coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Npts, dx, dy
    LOGICAL, INTENT(in)                                  :: Flike
    REAL(r_k), INTENT(in)                                :: NOval
    INTEGER, DIMENSION(Npts,2), INTENT(in)               :: pts
    REAL(r_k), DIMENSION(Npts), INTENT(in)               :: vals
    REAL(r_k), DIMENSION(dx,dy), INTENT(out)             :: matrix

! Local
    INTEGER                                              :: iv, i, j

!!!!!!! Variables
! Npts: Number of values of the list
! pts: 2D matrix coordinates of the values
! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
! vals: list of values correspondant to the coordinates
! dx, dy: shape of the 2D matrix
! NOval: Value to assign when there is no coordinate
! matrix: resultant matrix

    fname = 'from_ptlist_2DRKmatrix'

    matrix = NOval

    IF (Flike) THEN
      DO iv=1, Npts
        i = pts(iv,1)
        j = pts(iv,2)
        matrix(i,j) = vals(iv)
      END DO
    ELSE
      DO iv=1, Npts
        i = pts(iv,2)+1
        j = pts(iv,1)+1
        matrix(i,j) = vals(iv)
      END DO
    END IF

    RETURN

  END SUBROUTINE from_ptlist_2DRKmatrix

  SUBROUTINE from_ptlist_2DRKNmatrix(Npts, pts, Flike, Nvals, vals, dx, dy, NOval, Nmatrix)
  ! Subroutine to construct N 2D RK matrix from a list of values accompaigned by a list of grid-point 
  !   coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Npts, dx, dy, Nvals
    LOGICAL, INTENT(in)                                  :: Flike
    REAL(r_k), INTENT(in)                                :: NOval
    INTEGER, DIMENSION(Npts,2), INTENT(in)               :: pts
    REAL(r_k), DIMENSION(Npts,Nvals), INTENT(in)         :: vals
    REAL(r_k), DIMENSION(dx,dy,Nvals), INTENT(out)       :: Nmatrix

! Local
    INTEGER                                              :: iv, i, j, ic

!!!!!!! Variables
! Npts: Number of values of the list
! pts: 2D matrix coordinates of the values
! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
! Nvals: number of different values
! vals: list of N-values correspondant to the coordinates
! dx, dy: shape of the 2D matrix
! NOval: Value to assign when there is no coordinate
! Nmatrix: resultant N-matrix

    fname = 'from_ptlist_2DRKNmatrix'

    Nmatrix = NOval

    IF (Flike) THEN
      DO iv=1, Npts
        i = pts(iv,1)
        j = pts(iv,2)
        Nmatrix(i,j,:) = vals(iv,:)
      END DO
    ELSE
      DO iv=1, Npts
        i = pts(iv,2)+1
        j = pts(iv,1)+1
        Nmatrix(i,j,:) = vals(iv,:)
      END DO
    END IF

    RETURN

  END SUBROUTINE from_ptlist_2DRKNmatrix

  SUBROUTINE from_coordlist_2DRKmatrix(Npts, coords, Flike, xcoord, ycoord, vals, dx, dy, NOval,      &
    matrix)
  ! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of coordinates
  !   to find i,j grid-point coordinates by minimum distance

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Npts, dx, dy
    LOGICAL, INTENT(in)                                  :: Flike
    REAL(r_k), INTENT(in)                                :: NOval
    REAL(r_k), DIMENSION(Npts,2), INTENT(in)             :: coords
    REAL(r_k), DIMENSION(Npts), INTENT(in)               :: vals
    REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: xcoord, ycoord
    REAL(r_k), DIMENSION(dx,dy), INTENT(out)             :: matrix

! Local
    INTEGER                                              :: iv, i, j
    REAL(r_k)                                            :: xi, yi
    REAL(r_k), DIMENSION(dx,dy)                          :: diff
    INTEGER, DIMENSION(2)                                :: minpt

!!!!!!! Variables
! Npts: Number of values of the list
! coords: 2D coordinates of the values
! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
! vals: list of values correspondant to the coordinates
! xcoord, ycoord: matrix of values correspondant to the each coordinate
! dx, dy: shape of the 2D matrix
! NOval: Value to assign when there is no coordinate
! matrix: resultant matrix

    fname = 'from_coordlist_2DRKmatrix'

    matrix = NOval

    IF (Flike) THEN
      DO iv=1, Npts
        xi = coords(iv,1)
        yi = coords(iv,2)
        diff = SQRT((xcoord-xi)**2+(ycoord-yi)**2)
        minpt = MINLOC(diff)
        matrix(minpt(1),minpt(2)) = vals(iv)
      END DO
    ELSE
      DO iv=1, Npts
        xi = coords(iv,2)
        yi = coords(iv,1)
        diff = SQRT((xcoord-xi)**2+(ycoord-yi)**2)
        minpt = MINLOC(diff)
        matrix(minpt(1),minpt(2)) = vals(iv)
      END DO
    END IF

    RETURN

  END SUBROUTINE from_coordlist_2DRKmatrix

  INTEGER FUNCTION JD (year,month,day)
! Fucntion to compute the julian date (JD) given a gregorian calendar
!   FROM: https://aa.usno.navy.mil/faq/docs/JD_Formula.php

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: year,month,day

! Local
  INTEGER                                                :: I,J,K

  I= year
  J= month
  K= day

  JD= K-32075+1461*(I+4800+(J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4

  END FUNCTION JD

  SUBROUTINE GDATE(jd,year,month,day)
  ! Subroutine to compute the gregorian calendar date (year,month,day) given the julian date (JD)

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: jd
    INTEGER, INTENT(out)                                 :: year,month,day
! Local
    INTEGER                                              ::  I,J,K,L,N

    L= JD+68569
    N= 4*L/146097
    L= L-(146097*N+3)/4
    I= 4000*(L+1)/1461001
    L= L-1461*I/4+31
    J= 80*L/2447
    K= L-2447*J/80
    L= J/11
    J= J+2-12*L
    I= 100*(N-49)+I+L

    year= I
    month= J
    day= K

    RETURN

  END SUBROUTINE GDATE

  SUBROUTINE rm_values_vecRK(d1, vector, rmvalue, rmd1, rmvector)
  ! Subroutine to remove a given value (e.g. fill_Value) from a RK vector

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1
    REAL(r_k), INTENT(in)                                :: rmvalue
    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: vector
    INTEGER, INTENT(out)                                 :: rmd1
    REAL(r_k), DIMENSION(d1), INTENT(out)                :: rmvector

! Local
    INTEGER                                              :: i

!!!!!!! Variables
! d1: length of the vector
! vector: vector with values
! rmvalue: RK value to remove
! rmd1: length of the new vector
! rmvector: new vector

    fname = 'rm_values_vecRK'

    rmd1 = 0
    rmvector = zeroRK

    DO i=1, d1
      IF (vector(i) /= rmvalue) THEN
        rmd1 = rmd1 + 1
        rmvector(rmd1) = vector(i)
      END IF
    END DO

  END SUBROUTINE rm_values_vecRK

  SUBROUTINE index_samevals1D_RK(d1r, refv, d1v, vals, ii, indices, missval, samevalues)
  ! Subroutine to search for the indices of the same values between 2 1D RK series of values allowing 
  !   repetitions

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1r, d1v
    REAL(r_k), INTENT(in)                                :: missval
    REAL(r_k), DIMENSION(d1r), INTENT(in)                :: refv
    REAL(r_k), DIMENSION(d1v), INTENT(in)                :: vals
    INTEGER, INTENT(out)                                 :: ii
    INTEGER, DIMENSION(d1v,2), INTENT(out)               :: indices
    REAL(r_k), DIMENSION(d1v), INTENT(out)               :: samevalues

! Local
    INTEGER                                              :: ir, iv, iiv

!!!!!!! Variables
! d1r: size of the reference data
! refv: reference values
! d1v: size of the values
! vals: values to look in
! ii: quantity of same values found
! indices: output
! missval: missing value
! samevalues: values where coincidence is found

    fname = 'index_samevals1D_RK'

    indices = 0
    samevalues = zeroRK
    ii = 0
    DO ir=1, d1r
      IF (refv(ir) /= missval) THEN
        DO iv=1, d1v
          IF (vals(iv) == refv(ir)) THEN
            ii = ii + 1
            indices(ii,:) = (/ ir, iv /)
            samevalues(ii) = refv(ir)
          END IF
        END DO
      END IF
    END DO

    RETURN

  END SUBROUTINE index_samevals1D_RK

  SUBROUTINE index_samevals2D_RK(d1r, d2r, refv, d1v, d2v, d12v, vals, ii, indices, missval,          &
    samevalues)
  ! Subroutine to search for the indices of the same values between 2 2D RK series of values allowing 
  !   repetitions

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1r, d2r, d1v, d2v, d12v
    REAL(r_k), INTENT(in)                                :: missval
    REAL(r_k), DIMENSION(d1r,d2r), INTENT(in)            :: refv
    REAL(r_k), DIMENSION(d1v,d2v), INTENT(in)            :: vals
    INTEGER, INTENT(out)                                 :: ii
    INTEGER, DIMENSION(d12v,2,2), INTENT(out)            :: indices
    REAL(r_k), DIMENSION(d12v), INTENT(out)              :: samevalues

! Local
    INTEGER                                              :: ir1, ir2, iv1, iv2

!!!!!!! Variables
! d1r, d2r: size of the reference data
! refv: reference values
! d1v, d2v: size of the values
! vals: values to look in
! ii: quantity of same values found
! indices: output
! missval: missing value
! samevalues: values where coincidence is found

    fname = 'index_samevals2D_RK'

    indices = 0
    samevalues = zeroRK
    ii = 0
    DO ir1=1, d1r
      DO ir2=1, d2r
        IF (refv(ir1,ir2) /= missval) THEN
          DO iv1=1, d1v
            DO iv2=1, d2v
              IF (vals(iv1,iv2) == refv(ir1,ir2)) THEN
                ii = ii + 1
                indices(ii,1,1) = ir1
                indices(ii,1,2) = ir2
                indices(ii,2,1) = iv1
                indices(ii,2,2) = iv2
                samevalues(ii) = refv(ir1,ir2)
              END IF
            END DO
          END DO
        END IF
      END DO
    END DO

    RETURN

  END SUBROUTINE index_samevals2D_RK

  SUBROUTINE fill_matrix2DRK_winmat2D_list1D(d1i, d2i, inmatrix, ind, dlist, inlist, olist, missval,  &
    od, d1o, d2o, omat)
! Subroutine to fill a 2D RK matrix using a list of 1D indices from another given 2D matrix

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1i, d2i, ind, dlist, od, d1o, d2o
    REAL(r_k), DIMENSION(d1i, d2i), INTENT(in)           :: inmatrix
    INTEGER, DIMENSION(dlist), INTENT(in)                :: inlist, olist
    REAL(r_k), INTENT(in)                                :: missval
    REAL(r_k), DIMENSION(d1o, d2o), INTENT(out)          :: omat

! Local
    INTEGER                                              :: ii, ij, oi, oj, il
    INTEGER                                              :: isame, osame, irun, orun
    INTEGER                                              :: ilx, olx, iln, oln
    CHARACTER(len=3)                                     :: isS, osS

!!!!!!! Variables
! d1i, d2i: size of the input matrix
! inmatrix: input matrix with the values to fill the output matrix
! ind: dimension of the input matrix that the list of indices refer to
! dlist: number of indices from the list
! inlist: list of indices from the input matrix
! olist: list of indices from the output matrix
! missval: missing value
! od: dimension of the output matix to which assign the indices of the list
! d1o, d2o: size of the output matrix
! omat: output matrix

    fname = 'fill_matrix2DRK_winmat2D_list1D'

    omat = missval

    IF (ind == 1) THEN
      isame = d2i
      irun = d1i
    ELSE
      isame = d1i
      irun = d2i
    END IF

    IF (od == 1) THEN
      osame = d2o
      orun = d1o
    ELSE
      osame = d1o
      orun = d2o
    END IF

    IF (isame /= osame) THEN
      WRITE(isS,'(I3)')isame
      WRITE(osS,'(I3)')osame
      msg = 'Resultant working size from input ' // isS // ' to output ' // osS // ' differ !!'
      CALL ErrMsg(msg, fname, -1)
    END IF

    iln = MINVAL(inlist)
    oln = MINVAL(olist)

    IF (iln < 1) THEN
      WRITE(isS,'(I3)')iln
      msg = 'Incorrect minimum input index: ' // isS // ' Negative or zero indices not allowed !!'
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (oln < 1) THEN
      WRITE(isS,'(I3)')oln
      msg = 'Incorrect minimum output index: ' // isS // ' Negative or zero indices not allowed !!'
      CALL ErrMsg(msg, fname, -1)
    END IF

    ilx = MAXVAL(inlist)
    olx = MAXVAL(olist)

    IF (ilx > irun) THEN
      WRITE(isS,'(I3)')ilx
      WRITE(osS,'(I3)')irun
      msg = 'Maximum value in input indices ' // isS // ' larger than assigned dimension ' // osS //  &
        ' !!'
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (olx > orun) THEN
      WRITE(isS,'(I3)')olx
      WRITE(osS,'(I3)')orun
      msg = 'Maximum value in output indices ' // isS // ' larger than assigned dimension ' // osS // &
        ' !!'
      CALL ErrMsg(msg, fname, -1)
    END IF

    IF (od == 1) THEN
      IF (ind == 1) THEN
        DO il=1, dlist
          omat(olist(il),:) = inmatrix(inlist(il),:)
        END DO
      ELSE
        DO il=1, dlist
          omat(olist(il),:) = inmatrix(:,inlist(il))
        END DO
      END IF
    ELSE
      IF (ind == 1) THEN
        DO il=1, dlist
          omat(:,olist(il)) = inmatrix(inlist(il),:)
        END DO
      ELSE
        DO il=1, dlist
          omat(:,olist(il)) = inmatrix(:,inlist(il))
        END DO
      END IF
    END IF

    RETURN

  END SUBROUTINE fill_matrix2DRK_winmat2D_list1D


END MODULE module_generic