MODULE module_NCgeneric
! Module with generic netcdf functions

! create_NCfile: Subroutine to create a netCDF file
! handle_err: Subroutine to provide the error message when something with netCDF went wrong
! handle_errf: Subroutine to provide the error message when something with netCDF went wrong (including fname)
! isin_file: Function to tell if a given variable is inside a file
! isin_ncunit: Function to tell if a given variable is inside a netcdf file unit
! get_var2dims_file: Function to get the dimensions of a given 2D variable inside a file
! get_var3dims_file: Function to get the dimensions of a given 3D variable inside a file
! get_var4dims_file: Function to get the dimensions of a given 4D variable inside a file
! get_var1dims_ncunit: Function to get the dimensions of a given 1D variable inside a unit of a netCDF file
! get_var2dims_ncunit: Function to get the dimensions of a given 2D variable inside a unit of a netCDF file
! get_var3dims_ncunit: Function to get the dimensions of a given 3D variable inside a unit of a netCDF file
! get_varNdims_file: Function to get the number of dimensions of a given variable inside a file
! get_varNdims_ncunit: Function to get the number of dimensions of a given variable inside a unit of a netCDF file
! get_varI1D_ncunit: Subroutine to get a 1D integer variable from a netCDF file unit
! get_varI2D_ncunit: Subroutine to get a 2D integer variable from a netCDF file unit
! get_varRK0D_ncunit: Subroutine to get an scalar r_k float variable from a netCDF file unit
! get_varRK1D_ncunit: Subroutine to get a 1D r_k float variable from a netCDF file unit
! get_varRK2D_ncunit: Subroutine to get a 2D r_k float variable from a netCDF file unit
! get_varRK3D_ncunit: Subroutine to get a 3D r_k float variable from a netCDF file unit
! get_varRK2D_ncunit: Subroutine to get a 4D r_k float variable from a netCDF file unit
! put_var1D: Subroutine to write on a netCDF file a 1D float variable
! put_var2D: Subroutine to write on a netCDF file a 2D float variable
! put_var3D: Subroutine to write on a netCDF file a 3D float variable
! put_var1Dt: Subroutine to write on a netCDF file a 1D float variable at a given time-step
! put_var2Dt: Subroutine to write on a netCDF file a 2D float variable at a given time-step
! put_var3Dt: Subroutine to write on a netCDF file a 3D float variable at a given time-step

  USE module_definitions
  USE module_basic
  USE module_generic

  CONTAINS

! From UNIDATA: https://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html
  SUBROUTINE handle_err(st)
! Subroutine to provide the error message when something with netCDF went wrong

  USE netcdf

  INTEGER, INTENT(in)                                    :: st

!!!!!!! Variables
! fn: function name from which it is used

  IF (st /= nf90_noerr) THEN
    PRINT *, TRIM(emsg)
    PRINT *, '  ' // TRIM(nf90_strerror(st))
    STOP "Stopped"
  END IF

  END SUBROUTINE handle_err

! From UNIDATA: https://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html
  SUBROUTINE handle_errf(st, fn)
! Subroutine to provide the error message when something with netCDF went wrong (including fname)

  USE netcdf

  INTEGER, INTENT(in)                                    :: st
  CHARACTER(len=*), INTENT(in)                           :: fn

!!!!!!! Variables
! st: netCDF status number
! fn: function name from which it is used

  IF (st /= nf90_noerr) THEN
    PRINT *, TRIM(emsg)
    PRINT *, '  ' // TRIM(fn) // ': ' // TRIM(nf90_strerror(st))
    STOP "Stopped"
  END IF

  END SUBROUTINE handle_errf

  SUBROUTINE create_NCfile(filename, dimsfile, namelistfile, varsfile, ncid)
! Subroutine to create a netCDF file

    USE netcdf

    IMPLICIT NONE

    INCLUDE 'netcdf.inc'

    CHARACTER(LEN=*), INTENT(IN)                         :: filename, dimsfile, namelistfile, varsfile
    INTEGER, INTENT(OUT)                                 :: ncid

! Local
    INTEGER                                              :: i, j, k, idimnew
    INTEGER                                              :: rcode, funit, funit2, ios
    INTEGER                                              :: Nvals, dimsize, dimid, iddimnew, Ntotdims
    INTEGER                                              :: idvarnew, vartype
    CHARACTER(LEN=200)                                   :: message, vd, vs, vdd, val
    CHARACTER(LEN=200)                                   :: vname, Lvname, vunits, coornames
    CHARACTER(LEN=200), DIMENSION(:), ALLOCATABLE        :: valsline, dimsizes
    CHARACTER(LEN=1000)                                  :: line, dimsline
    INTEGER, DIMENSION(:), ALLOCATABLE                   :: dimsvar
    INTEGER                                              :: Ldimsize, Ldimsvar, dvarL
    CHARACTER(LEN=1)                                     :: dvarn

!!!!!!! Variables
! filename: name of the file to create
! dimsfile: ASCII file with the name of the dimensions to create with ('#' for comentaries)
!   [dim name]| [dim orig in WRF]| [dim orig in namelist]| ['unlimited' also, 'namelist' (from namelist parameter)]
! namelistfile: name of the Namelist file
! varsfile: ASCII file with the name of the variables to create with ('#' for comentaries)
!   [WRFvarname]| [var name]| [long var name]| [var units]| [var dimensions]
! ncid: number assigned to the file

    fname = 'create_NCfile'

! Opening creation status
    rcode = nf90_create(TRIM(filename), NF90_NETCDF4, ncid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

! Reading dimensions file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(dimsfile), STATUS='old', FORM='formatted', IOSTAT=ios)
    message = "Problems to open varibales file '" // TRIM(varsfile) // "' IOSTAT = "//                &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(message, fname)

    Nvals = 4
    IF (ALLOCATED(valsline)) DEALLOCATE(valsline)
    ALLOCATE (valsline(Nvals))

! Creation of dimensions
    idimnew = 3
    dimsline = ''
    Ntotdims = 0
    DO i=1,1000
      READ(funit, '(A1000)', END=100)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        CALL removeChar(valsline(4),' ')
        IF (TRIM(valsline(4)) == 'unlimited') THEN
          idimnew = idimnew + 1
          dimsize = NF90_UNLIMITED
          dimid = idimnew
        ELSE IF (TRIM(valsline(4)) == 'namelist') THEN
          CALL GetInNamelist(namelistfile, valsline(3), 'I', Ival=dimsize)
          SELECT CASE (TRIM(valsline(2)))
            CASE ('i')
              dimid = 1
            CASE ('j')
              dimid = 2
            CASE ('k')
              dimid = 3
            CASE ('t')
              dimid = 4
              dimsize = NF90_UNLIMITED
            CASE DEFAULT
              idimnew = idimnew + 1
              dimid = idimnew
          END SELECT
        END IF
        rcode = nf90_def_dim(ncid, TRIM(valsline(1)), dimsize, dimid)
        IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
        vs = valsline(2)
        CALL removeChar(vs, ' ')
        CALL attachString(dimsline, TRIM(vs) // ':' // TRIM(ItoS(dimid)) // ';')
        Ntotdims = Ntotdims + 1
      END IF
    END DO

 100 CONTINUE
    CLOSE(funit)

! Sort of python dictionary for [dimn]:[dimsize]...
    IF (ALLOCATED(dimsizes)) DEALLOCATE(dimsizes)
    ALLOCATE(dimsizes(Ntotdims))
    CALL split(dimsline,';',Ntotdims,dimsizes)

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(varsfile), STATUS='old', FORM='formatted', IOSTAT=ios)
    
    message = "Problems to open varibales file '" // TRIM(varsfile) // "' IOSTAT = "//                &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(message, fname)

    Nvals = 6
    IF (ALLOCATED(valsline)) DEALLOCATE(valsline)
    ALLOCATE (valsline(Nvals))

! Defining variables
    idvarnew = 1
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        vtype: SELECT CASE (TRIM(valsline(6)))
          CASE ('B')
            vartype = NF_BYTE
          CASE ('C')
            vartype = NF_CHAR
          CASE ('I')
            vartype = NF_SHORT
          CASE ('I16')
            vartype = NF_INT
          CASE ('R')
            vartype = NF_FLOAT
          CASE ('R16')
            vartype = NF_DOUBLE 
        END SELECT vtype

        vd = valsline(5)
        CALL removeChar(vd, ' ')
        Ldimsvar = LEN_TRIM(vd)
        IF (ALLOCATED(dimsvar)) DEALLOCATE(dimsvar)
        ALLOCATE(dimsvar(Ldimsvar))

! Variable's dimensions
        coornames = ''
        DO j=1, Ldimsvar
          DO k=1, Ntotdims
            IF (dimsizes(k)(1:1) == vd(j:j)) THEN
              Ldimsize = LEN_TRIM(dimsizes(k))
              vdd = dimsizes(k)(3:Ldimsize)
              dimsvar(j) = StoI(vdd)
! Too complicated to assign dimvarname... (or too lazy)
!              coornames = coornames // 
              CYCLE
            END IF
          END DO
        END DO
        vname = valsline(2)
        CALL removeChar(vname, ' ')
        vartype = 5
       
        rcode = nf90_def_var(ncid, vname, vartype, dimsvar, idvarnew)
        IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

! Adding attributes
        rcode = nf90_put_att(ncid, idvarnew, 'standard_name', TRIM(valsline(2)))
        IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
        rcode = nf90_put_att(ncid, idvarnew, 'long_name', TRIM(valsline(3)))
        IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
        rcode = nf90_put_att(ncid, idvarnew, 'units', TRIM(valsline(4)))
        IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

        idvarnew = idvarnew + 1 
      END IF
    END DO

 150 CONTINUE
    CLOSE(funit)

    rcode = NF90_ENDDEF(ncid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    DEALLOCATE(valsline)
    DEALLOCATE(dimsizes)

  END SUBROUTINE create_NCfile

  FUNCTION get_var2dims_file(filename, varname)
! Function to get the dimensions of a given 2D variable inside a file

    USE netcdf

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(in)                         :: filename, varname
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(2)                                :: get_var2dims_file

! Local
    INTEGER                                              :: nid, vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(2)                                :: dimsid

!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_var2dims_file'
    !PRINT *,TRIM(fname)

! Opening creation status
    rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    IF (Ndims /= 2) THEN
      msg = "variable '" // TRIM(varname) // "' has not 2 dimensions!!"
      CALL stoprun(msg, fname)
    END IF

    rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = get_var2dims_file(1))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = get_var2dims_file(2))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = NF90_CLOSE(nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_var2dims_file

  FUNCTION get_var3dims_file(filename, varname)
! Function to get the dimensions of a given 3D variable inside a file

    USE netcdf

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(in)                         :: filename, varname
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(3)                                :: get_var3dims_file

! Local
    INTEGER                                              :: nid, vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(3)                                :: dimsid


!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_var3dims_file'
    !PRINT *,TRIM(fname)

! Opening creation status
    rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

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

    rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = get_var3dims_file(1))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = get_var3dims_file(2))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = get_var3dims_file(3))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = NF90_CLOSE(nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_var3dims_file

  FUNCTION get_var1dims_ncunit(nid, varname)
! Function to get the dimensions of a given 1D variable inside a unit of a netCDF file

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: nid
    CHARACTER(LEN=*), INTENT(in)                         :: varname
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(1)                                :: get_var1dims_ncunit

! Local
    INTEGER                                              :: vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(1)                                :: dimsid


!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_var1dims_ncunit'
    !PRINT *,TRIM(fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    IF (Ndims /= 1) THEN
      msg = "variable '" // TRIM(varname) // "' has not 1 dimensions!!"
      CALL stoprun(msg, fname)
    END IF

    rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(1), name=msg)

    rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = get_var1dims_ncunit(1))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_var1dims_ncunit

  FUNCTION get_var2dims_ncunit(nid, varname)
! Function to get the dimensions of a given 2D variable inside a unit of a netCDF file

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: nid
    CHARACTER(LEN=*), INTENT(in)                         :: varname
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(2)                                :: get_var2dims_ncunit

! Local
    INTEGER                                              :: vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(2)                                :: dimsid


!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_var2dims_ncunit'
    !PRINT *,TRIM(fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    IF (Ndims /= 2) THEN
      msg = "variable '" // TRIM(varname) // "' has not 2 dimensions!!"
      CALL stoprun(msg, fname)
    END IF

    rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = get_var2dims_ncunit(1))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = get_var2dims_ncunit(2))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_var2dims_ncunit

  FUNCTION get_var3dims_ncunit(nid, varname)
! Function to get the dimensions of a given 3D variable inside a unit of a netCDF file

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: nid
    CHARACTER(LEN=*), INTENT(in)                         :: varname
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(3)                                :: get_var3dims_ncunit

! Local
    INTEGER                                              :: vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(3)                                :: dimsid


!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_var3dims_ncunit'
    !PRINT *,TRIM(fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

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

    rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = get_var3dims_ncunit(1))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = get_var3dims_ncunit(2))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = get_var3dims_ncunit(3))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_var3dims_ncunit

  FUNCTION get_var4dims_file(filename, varname)
! Function to get the dimensions of a given 4D variable inside a file

    USE netcdf

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(in)                         :: filename, varname
! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
    INTEGER, DIMENSION(4)                                :: get_var4dims_file

! Local
    INTEGER                                              :: nid, vid, Ndims
    INTEGER                                              :: rcode
    INTEGER, DIMENSION(4)                                :: dimsid


!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_var4dims_file'
    !PRINT *,TRIM(fname)

! Opening creation status
    rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

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

    rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = get_var4dims_file(1))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = get_var4dims_file(2))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = get_var4dims_file(3))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_dimension(nid, dimsid(4), LEN = get_var4dims_file(4))
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = NF90_CLOSE(nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_var4dims_file

  INTEGER FUNCTION get_varNdims_file(filename, varname)
! Function to get the number of dimensions of a given variable inside a file

    USE netcdf

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(in)                         :: filename, varname

! Local
    INTEGER                                              :: nid, vid
    INTEGER                                              :: rcode

!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_varNdims_file'
    !PRINT *,TRIM(fname)

! Opening creation status
    rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = get_varNdims_file)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = NF90_CLOSE(nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_varNdims_file

  INTEGER FUNCTION get_varNdims_ncunit(nid, varname)
! Function to get the number of dimensions of a given variable inside a unit of a netCDF file

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: nid
    CHARACTER(LEN=*), INTENT(in)                         :: varname

! Local
    INTEGER                                              :: vid
    INTEGER                                              :: rcode

!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'get_varNdims_ncunit'
    !PRINT *,TRIM(fname)

    rcode = nf90_inq_varid(nid, varname, vid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire_variable(nid, vid, NDIMS = get_varNdims_ncunit)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END FUNCTION get_varNdims_ncunit

LOGICAL FUNCTION isin_file(filename, varname)
! Function to tell if a given variable is inside a file

    USE netcdf

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(in)                         :: filename, varname

! Local
    INTEGER                                              :: nid, vid, Ndims, Nvars
    INTEGER                                              :: iv, rcode
    CHARACTER(LEN=1000)                                  :: varinfile

!!!!!!! Variables
! filename: name of the file to open
! varname: name of the variable 

    fname = 'isin_file'

! Opening creation status
    rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_inquire(nid, Ndims, Nvars)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    DO iv=1, Nvars
      rcode = nf90_inquire_variable(nid, iv, name=varinfile)
      IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
      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 handle_errf(rcode, fname)
    
  END FUNCTION isin_file

LOGICAL FUNCTION isin_ncunit(nid, varname)
! Function to tell if a given variable is inside a netcdf file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: nid
    CHARACTER(LEN=*), INTENT(in)                         :: varname

! Local
    INTEGER                                              :: vid, Ndims, Nvars
    INTEGER                                              :: iv, rcode
    CHARACTER(LEN=1000)                                  :: varinfile

!!!!!!! Variables
! nid: number of the opened netCDF
! varname: name of the variable 

    fname = 'isin_ncunit'

    rcode = nf90_inquire(nid, Ndims, Nvars)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    DO iv=1, Nvars
      rcode = nf90_inquire_variable(nid, iv, name=varinfile)
      IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
      IF (TRIM(varinfile) == TRIM(varname)) THEN
        isin_ncunit = .TRUE.
        EXIT
      ELSE
        isin_ncunit = .FALSE.
      END IF
    END DO
    
  END FUNCTION isin_ncunit

  SUBROUTINE put_var1D(ncid, d1, vals, vname, filevarn)
! Subroutine to write on a netCDF file a 1D float variable

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(IN)                                  :: ncid, d1
    REAL, DIMENSION(d1), INTENT(IN)                      :: vals
    CHARACTER(LEN=*), INTENT(IN)                         :: vname, filevarn
  
! Local
    INTEGER                                              :: funit, i, idvarnew, ios
    INTEGER                                              :: Nvals, rcode, varid
    CHARACTER(LEN=50)                                    :: ncvarname
    CHARACTER(LEN=1000)                                  :: line
    CHARACTER(LEN=200), DIMENSION(6)                     :: valsline
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to include
! vname: name of the variable in the model to be included
! filevarn: name of the ASCII file with the information about the variables
    fname = 'put_var1D'

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(filevarn), STATUS='old', FORM='formatted', IOSTAT=ios)
    msg = "Problems to open variables file '" // TRIM(filevarn) // "' IOSTAT = "//                    &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(msg, fname)

    Nvals = 6

    idvarnew = 1
    vfound = .FALSE.
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        IF (TRIM(vname) == TRIM(valsline(1))) THEN
          ncvarname = TRIM(valsline(2))
          CALL removeChar(ncvarname, ' ')
          rcode = nf90_inq_varid(ncid, ncvarname, varid)
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

          rcode = nf90_put_var(ncid, varid, vals)          
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
          vfound = .TRUE.
          CYCLE
        END IF
      END IF
    END DO

 150 CONTINUE

    CLOSE(funit)
    IF (.NOT.vfound) THEN
      msg = "variables file '" // TRIM(filevarn) // "' does not have varible '" // TRIM(vname) //     &
        "' !!"
      CALL stoprun(msg, fname)
    END IF

  END SUBROUTINE put_var1D

  SUBROUTINE put_var2D(ncid, d1, d2, vals, vname, filevarn)
! Subroutine to write on a netCDF file a 2D float variable

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(IN)                                  :: ncid, d1, d2
    REAL, DIMENSION(d1,d2), INTENT(IN)                   :: vals
    CHARACTER(LEN=*), INTENT(IN)                         :: vname, filevarn
  
! Local
    INTEGER                                              :: funit, i, idvarnew, ios
    INTEGER                                              :: Nvals, rcode, varid
    CHARACTER(LEN=50)                                    :: ncvarname
    CHARACTER(LEN=1000)                                  :: line
    CHARACTER(LEN=200), DIMENSION(6)                     :: valsline
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1,d2: shape of the matrix
! vals: values to include
! vname: name of the variable in the model to be included
! filevarn: name of the ASCII file with the information about the variables
    fname = 'put_var2D'

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(filevarn), STATUS='old', FORM='formatted', IOSTAT=ios)
    msg = "Problems to open variables file '" // TRIM(filevarn) // "' IOSTAT = "//                    &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(msg, fname)

    Nvals = 6

    idvarnew = 1
    vfound = .FALSE.
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        IF (TRIM(vname) == TRIM(valsline(1))) THEN
          ncvarname = TRIM(valsline(2))
          CALL removeChar(ncvarname, ' ')
          rcode = nf90_inq_varid(ncid, ncvarname, varid)
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

          rcode = nf90_put_var(ncid, varid, vals)          
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
          vfound = .TRUE.
          CYCLE
        END IF
      END IF
    END DO

 150 CONTINUE

    CLOSE(funit)
    IF (.NOT.vfound) THEN
      msg = "variables file '" // TRIM(filevarn) // "' does not have varible '" // TRIM(vname) //     &
        "' !!"
      CALL stoprun(msg, fname)
    END IF

  END SUBROUTINE put_var2D

  SUBROUTINE put_var3D(ncid, d1, d2, d3, vals, vname, filevarn)
! Subroutine to write on a netCDF file a 3D float variable

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(IN)                                  :: ncid, d1, d2, d3
    REAL, DIMENSION(d1,d2,d3), INTENT(IN)                :: vals
    CHARACTER(LEN=*), INTENT(IN)                         :: vname, filevarn
  
! Local
    INTEGER                                              :: funit, i, idvarnew, ios
    INTEGER                                              :: Nvals, rcode, varid
    CHARACTER(LEN=50)                                    :: ncvarname
    CHARACTER(LEN=1000)                                  :: line
    CHARACTER(LEN=200), DIMENSION(6)                     :: valsline
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1,d2,d3: shape of the matrix
! vals: values to include
! vname: name of the variable in the model to be included
! filevarn: name of the ASCII file with the information about the variables
    fname = 'put_var3D'

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(filevarn), STATUS='old', FORM='formatted', IOSTAT=ios)
    msg = "Problems to open variables file '" // TRIM(filevarn) // "' IOSTAT = "//                    &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(msg, fname)

    Nvals = 6

    idvarnew = 1
    vfound = .FALSE.
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        IF (TRIM(vname) == TRIM(valsline(1))) THEN
          ncvarname = TRIM(valsline(2))
          CALL removeChar(ncvarname, ' ')
          rcode = nf90_inq_varid(ncid, ncvarname, varid)
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

          rcode = nf90_put_var(ncid, varid, vals)          
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
          vfound = .TRUE.
          CYCLE
        END IF
      END IF
    END DO

 150 CONTINUE
    CLOSE(funit)

    IF (.NOT.vfound) THEN
      msg = "variables file '" // TRIM(filevarn) // "' does not have varible '" // TRIM(vname) //     &
        "' !!"
      CALL stoprun(msg, fname)
    END IF

  END SUBROUTINE put_var3D

  SUBROUTINE put_var1Dt(ncid, d1, vals, vname, filevarn, it)
! Subroutine to write on a netCDF file a 1D float variable at a given time-step

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(IN)                                  :: ncid, d1, it
    REAL, DIMENSION(d1), INTENT(IN)                      :: vals
    CHARACTER(LEN=*), INTENT(IN)                         :: vname, filevarn
  
! Local
    INTEGER                                              :: funit, i, idvarnew, ios
    INTEGER                                              :: Nvals, rcode, varid
    CHARACTER(LEN=50)                                    :: ncvarname
    CHARACTER(LEN=1000)                                  :: line
    CHARACTER(LEN=200), DIMENSION(6)                     :: valsline
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to include
! vname: name of the variable in the model to be included
! filevarn: name of the ASCII file with the information about the variables
! it: time-step to add

    fname = 'put_var1Dt'

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(filevarn), STATUS='old', FORM='formatted', IOSTAT=ios)
    msg = "Problems to open variables file '" // TRIM(filevarn) // "' IOSTAT = "//                    &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(msg, fname)

    Nvals = 6

    idvarnew = 1
    vfound = .FALSE.
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        IF (TRIM(vname) == TRIM(valsline(1))) THEN
          ncvarname = TRIM(valsline(2))
          CALL removeChar(ncvarname, ' ')
          rcode = nf90_inq_varid(ncid, ncvarname, varid)
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

          rcode = nf90_put_var(ncid, varid, vals, start=(/1,it/), count=(/d1,1/))
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
          vfound = .TRUE.
          CYCLE
        END IF
      END IF
    END DO

 150 CONTINUE

    CLOSE(funit)
    IF (.NOT.vfound) THEN
      msg = "variables file '" // TRIM(filevarn) // "' does not have varible '" // TRIM(vname) //     &
        "' !!"
      CALL stoprun(msg, fname)
    END IF

  END SUBROUTINE put_var1Dt

  SUBROUTINE put_var2Dt(ncid, d1, d2, vals, vname, filevarn, it)
! Subroutine to write on a netCDF file a 2D float variable at a given time-step

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(IN)                                  :: ncid, d1, d2, it
    REAL, DIMENSION(d1), INTENT(IN)                      :: vals
    CHARACTER(LEN=*), INTENT(IN)                         :: vname, filevarn
  
! Local
    INTEGER                                              :: funit, i, idvarnew, ios
    INTEGER                                              :: Nvals, rcode, varid
    CHARACTER(LEN=50)                                    :: ncvarname
    CHARACTER(LEN=1000)                                  :: line
    CHARACTER(LEN=200), DIMENSION(6)                     :: valsline
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to include
! vname: name of the variable in the model to be included
! filevarn: name of the ASCII file with the information about the variables
! it: time-step to add

    fname = 'put_var2Dt'

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(filevarn), STATUS='old', FORM='formatted', IOSTAT=ios)
    msg = "Problems to open variables file '" // TRIM(filevarn) // "' IOSTAT = "//                    &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(msg, fname)

    Nvals = 6

    idvarnew = 1
    vfound = .FALSE.
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        IF (TRIM(vname) == TRIM(valsline(1))) THEN
          ncvarname = TRIM(valsline(2))
          CALL removeChar(ncvarname, ' ')
          rcode = nf90_inq_varid(ncid, ncvarname, varid)
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

          rcode = nf90_put_var(ncid, varid, vals, start=(/1,1,it/), count=(/d1,d2,1/))
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
          vfound = .TRUE.
          CYCLE
        END IF
      END IF
    END DO

 150 CONTINUE

    CLOSE(funit)
    IF (.NOT.vfound) THEN
      msg = "variables file '" // TRIM(filevarn) // "' does not have varible '" // TRIM(vname) //     &
        "' !!"
      CALL stoprun(msg, fname)
    END IF

  END SUBROUTINE put_var2Dt

  SUBROUTINE put_var3Dt(ncid, d1, d2, d3, vals, vname, filevarn, it)
! Subroutine to write on a netCDF file a 3D float variable at a given time-step

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(IN)                                  :: ncid, d1, d2, d3, it
    REAL, DIMENSION(d1), INTENT(IN)                      :: vals
    CHARACTER(LEN=*), INTENT(IN)                         :: vname, filevarn
  
! Local
    INTEGER                                              :: funit, i, idvarnew, ios
    INTEGER                                              :: Nvals, rcode, varid
    CHARACTER(LEN=50)                                    :: ncvarname
    CHARACTER(LEN=1000)                                  :: line
    CHARACTER(LEN=200), DIMENSION(6)                     :: valsline
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1,d2,d3: shape of the matrix
! vals: values to include
! vname: name of the variable in the model to be included
! filevarn: name of the ASCII file with the information about the variables
! it: time-step to add

    fname = 'put_var3Dt'

! Reading variables file
    funit = freeunit()
    OPEN(funit, FILE=TRIM(filevarn), STATUS='old', FORM='formatted', IOSTAT=ios)
    msg = "Problems to open variables file '" // TRIM(filevarn) // "' IOSTAT = "//                    &
      TRIM(ItoS(ios)) // " !!"
    IF ( ios /= 0 ) CALL stoprun(msg, fname)

    Nvals = 6

    idvarnew = 1
    vfound = .FALSE.
    DO i=1,1000
      READ(funit, '(A1000)', END=150)line
      IF (line(1:1) /= '#' .AND. LEN_TRIM(line) > 1) THEN
        CALL split(line,'|',Nvals,valsline)
        IF (TRIM(vname) == TRIM(valsline(1))) THEN
          ncvarname = TRIM(valsline(2))
          CALL removeChar(ncvarname, ' ')
          rcode = nf90_inq_varid(ncid, ncvarname, varid)
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

          rcode = nf90_put_var(ncid, varid, vals, start=(/1,1,1,it/), count=(/d1,d2,d3,1/))
          IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)
          vfound = .TRUE.
          CYCLE
        END IF
      END IF
    END DO

 150 CONTINUE

    CLOSE(funit)
    IF (.NOT.vfound) THEN
      msg = "variables file '" // TRIM(filevarn) // "' does not have varible '" // TRIM(vname) //     &
        "' !!"
      CALL stoprun(msg, fname)
    END IF

  END SUBROUTINE put_var3Dt

  SUBROUTINE get_varI1D_ncunit(ncid, d1, vname, vals)
! Subroutine to get a 1D integer variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    INTEGER, DIMENSION(d1), INTENT(out)                  :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to get
! vname: name of the variable to getºº

    fname = 'get_varI1D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varI1D_ncunit

  SUBROUTINE get_varI2D_ncunit(ncid, d1, d2, vname, vals)
! Subroutine to get a 2D integer variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1, d2
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    INTEGER, DIMENSION(d1,d2), INTENT(out)               :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varI2D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varI2D_ncunit

  SUBROUTINE get_varI3D_ncunit(ncid, d1, d2, d3, vname, vals)
! Subroutine to get a 2D integer variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1, d2, d3
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    INTEGER, DIMENSION(d1,d2,d3), INTENT(out)            :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varI3D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varI3D_ncunit

  SUBROUTINE get_varRK0D_ncunit(ncid, vname, vals)
! Subroutine to get an scalar r_k float variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    REAL, INTENT(out)                                    :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varRK0D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varRK0D_ncunit

  SUBROUTINE get_varRK1D_ncunit(ncid, d1, vname, vals)
! Subroutine to get a 1D r_k float variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    REAL, DIMENSION(d1), INTENT(out)                     :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1: shape of the matrix
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varRK1D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varRK1D_ncunit

  SUBROUTINE get_varRK2D_ncunit(ncid, d1, d2, vname, vals)
! Subroutine to get a 2D r_k float variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1, d2
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    REAL, DIMENSION(d1,d2), INTENT(out)                  :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1,d2: shape of the matrix
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varRK2D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varRK2D_ncunit

  SUBROUTINE get_varRK3D_ncunit(ncid, d1, d2, d3, vname, vals)
! Subroutine to get a 3D r_k float variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1, d2, d3
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    REAL, DIMENSION(d1,d2,d3), INTENT(out)               :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1,d2,d3: shape of the matrix
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varRK3D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varRK3D_ncunit

  SUBROUTINE get_varRK4D_ncunit(ncid, d1, d2, d3, d4, vname, vals)
! Subroutine to get a 4D r_k float variable from a netCDF file unit

    USE netcdf

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: ncid, d1, d2, d3, d4
    CHARACTER(LEN=*), INTENT(in)                         :: vname
    REAL, DIMENSION(d1,d2,d3,d4), INTENT(out)            :: vals
  
! Local
    INTEGER                                              :: rcode, varid
    LOGICAL                                              :: vfound

!!!!!!! Variables
! ncid: netCDF file identifier
! d1,d2,d3,d4: shape of the matrix
! vals: values to get
! vname: name of the variable to get

    fname = 'get_varRK4D_ncunit'

    vfound = isin_ncunit(ncid, vname)

    IF (.NOT.vfound) THEN
      msg = "Unit file does not have variable '" // TRIM(vname) // "'"
      CALL ErrMsg(msg, fname, -1)
    END IF

    rcode = nf90_inq_varid(ncid, vname, varid)
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

    rcode = nf90_get_var(ncid, varid, vals)          
    IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname)

  END SUBROUTINE get_varRK4D_ncunit

END MODULE module_NCgeneric