module read_input_file

!-----------------------------------------------------------------------
!
!     contains routines for reading a file to interpolate
!
!-----------------------------------------------------------------------

      use kinds_mod     ! defines common data types
      use constants     ! defines useful constants
      use grids         ! includes all grid information
      use netcdf_mod    ! module with netcdf vars and utilities
      use remap_vars    ! module for all required remapping variables

      implicit none

      real (kind=dbl_kind), dimension(:,:), allocatable ::
     &    grid1_array

      integer (kind=int_kind) ::
     &    ntime, nlat, nlon, natts, nvaratts, nglobatts

      integer (kind=int_kind) ::
     &    nc_infile_id, nc_time_id, nc_invar_id, dimtime

      real (kind=dbl_kind), dimension(:), allocatable ::
     &    time

      logical :: invertlat


      contains

!***********************************************************************

      subroutine read_input(infile, varname)

!-----------------------------------------------------------------------
!
!     input variables
!
!-----------------------------------------------------------------------

      character(char_len), intent(in) ::
     &    infile, varname       ! input filename an varname

!-----------------------------------------------------------------------
    
      character(char_len) ::
     &    timename, name2    

      real (kind=dbl_kind), dimension(:,:,:), allocatable ::
     &    grid1_array_tmp

      integer, dimension(:,:,:), allocatable ::
     &    grid1_array_tmp1 
 
      real, dimension(:,:,:), allocatable ::
     &    grid1_array_tmp2

      double precision, dimension(:,:,:), allocatable ::
     &    grid1_array_tmp3

      integer, dimension(:), allocatable ::
     &    time1

      real, dimension(:), allocatable ::
     &    time2

      double precision, dimension(:), allocatable ::
     &    time3


      integer (kind=int_kind) :: 
     &    ncstat, nc_var_type, ndims, ji

      integer (kind=int_kind), dimension(:), allocatable :: ! netCDF ids
     &    nc_dims_ids

!-----------------------------------------------------------------------
!
!     read dimension information
!
!-----------------------------------------------------------------------

      ncstat = nf_open(infile, NF_NOWRITE, nc_infile_id)
      call netcdf_error_handler(ncstat)

      ncstat = nf_inq_varid(nc_infile_id, varname, nc_invar_id)
      call netcdf_error_handler(ncstat)
   
      ncstat = nf_inq_varndims(nc_infile_id, nc_invar_id, ndims)
      call netcdf_error_handler(ncstat)

      allocate(nc_dims_ids(ndims))
      ncstat = nf_inq_vardimid(nc_infile_id, nc_invar_id, nc_dims_ids)
      call netcdf_error_handler(ncstat)
      if ( ndims < 1 .or. ndims > 3) then
        stop "Input files should have (lon) or (lon,lat) or (lon,time) 
     &        or (lon,lat,time) dimensions"
      endif

      ncstat = nf_inq_dimlen(nc_infile_id, nc_dims_ids(1), nlon)
      call netcdf_error_handler(ncstat)
      if (nlon .ne. grid1_dims(1)) then
        print*," Input file : ",nlon," longitudes"
        print*," Weight file : ",grid1_dims(1) ," longitudes"
        stop " Inconsistent number of longitudes. We stop "
      endif

      name2='none'
      nlat=1
      if ( ndims >= 2 ) then
        ncstat = nf_inq_dimname(nc_infile_id, nc_dims_ids(2), name2)
        
        if ( name2 /= 'time' .and. name2 /= 'time_counter' .and. 
     &     name2 /= 't' ) then
 
          ncstat = nf_inq_dimlen(nc_infile_id, nc_dims_ids(2), nlat)
          call netcdf_error_handler(ncstat)
          if (nlat .ne. grid1_dims(2)) then
            print*," Input file : ",nlat," latitudes"
            print*," Weight file : ",grid1_dims(2)," latitudes"
            stop " Inconsistent number of latitudes. We stop "
          endif
          dimtime=3
        else
          dimtime=2 
        endif
      endif
!-----------------------------------------------------------------------
!
!     read time dimension 
!
!-----------------------------------------------------------------------

      nc_time_id=-1
      if ( ndims == 3 .or. dimtime == 2 ) then
        ncstat = nf_inq_dim(nc_infile_id, nc_dims_ids(dimtime), 
     &           timename, ntime)
        call netcdf_error_handler(ncstat)
        
        allocate(time(ntime))

        ncstat = nf_inq_varid(nc_infile_id, timename, nc_time_id)
        if (ncstat == NF_NOERR) then

          ncstat = nf_inq_vartype(nc_infile_id, nc_time_id, nc_var_type)
          call netcdf_error_handler(ncstat)
        
          allocate(time1(ntime))
          allocate(time2(ntime))
          allocate(time3(ntime))
          if ( nc_var_type == NF_INT ) then
            ncstat = nf_get_var_int(nc_infile_id, nc_time_id, time1)
            if ( ncstat == NF_NOERR ) then
              time = real(time1,kind=dbl_kind)
            else
              stop "Problem reading input data"
            endif
          elseif (nc_var_type == NF_REAL) then
            ncstat = nf_get_var_real(nc_infile_id, nc_time_id, time2)
            if ( ncstat == NF_NOERR ) then
              time = real(time2,kind=dbl_kind)
            else
              stop "Problem reading input data"
            endif
          elseif ( nc_var_type == NF_DOUBLE) then
            ncstat = nf_get_var_double(nc_infile_id,nc_time_id, time3)
            if ( ncstat == NF_NOERR ) then
              time = real(time3,kind=dbl_kind)
            else
              stop"Problem reading input data"
            endif
          else
            stop"Problem with input data type"
          endif 
          deallocate(time1, time2, time3)

          ncstat = nf_inq_varnatts(nc_infile_id, nc_time_id, natts)
          call netcdf_error_handler(ncstat)
        else
          do ji=1,ntime
            time(ji)=ji
          enddo
          natts=0
        endif

      else
        ntime=1
        allocate(time(1))
        time=1.
        natts=0
      endif

!-----------------------------------------------------------------------
!
!     allocate arrays
!
!-----------------------------------------------------------------------

      allocate( grid1_array_tmp (nlon,nlat,ntime))
      allocate( grid1_array_tmp1 (nlon,nlat,ntime))
      allocate( grid1_array_tmp2 (nlon,nlat,ntime))
      allocate( grid1_array_tmp3 (nlon,nlat,ntime))

!-----------------------------------------------------------------------
!
!     read variable 
!
!-----------------------------------------------------------------------
 
      ncstat = nf_inq_vartype(nc_infile_id, nc_invar_id, nc_var_type) 
      call netcdf_error_handler(ncstat)

      if ( nc_var_type == NF_INT ) then
        ncstat = nf_get_var_int(nc_infile_id, nc_invar_id, 
     &                          grid1_array_tmp1)
        if ( ncstat == NF_NOERR ) then
          grid1_array_tmp = real(grid1_array_tmp1,kind=dbl_kind)
        else
          stop "Problem reading input data"
        endif
      elseif (nc_var_type == NF_REAL) then
        ncstat = nf_get_var_real(nc_infile_id, nc_invar_id, 
     &                           grid1_array_tmp2)
        if ( ncstat == NF_NOERR ) then
          grid1_array_tmp = real(grid1_array_tmp2,kind=dbl_kind)
        else
          stop "Problem reading input data"
        endif
      elseif ( nc_var_type == NF_DOUBLE) then
        ncstat = nf_get_var_double(nc_infile_id,nc_invar_id, 
     &                             grid1_array_tmp3)
        if ( ncstat == NF_NOERR ) then
          grid1_array_tmp = real(grid1_array_tmp3,kind=dbl_kind)
        else
          stop"Problem reading input data"
        endif
      else
        stop"Problem with input data type"
      endif

      ncstat = nf_inq_varnatts(nc_infile_id, nc_invar_id, nvaratts)
      call netcdf_error_handler(ncstat)
 
      deallocate(grid1_array_tmp1, grid1_array_tmp2, grid1_array_tmp3)

!-----------------------------------------------------------------------
!
!     reshape input file
!
!-----------------------------------------------------------------------
     
      allocate( grid1_array (nlon*nlat,ntime)) 
      if (invertlat) then
        grid1_array = reshape ( grid1_array_tmp (:,nlat:1:-1,:), 
     &                         (/ nlat*nlon, ntime /) )
      else
        grid1_array = reshape (grid1_array_tmp, (/ nlat*nlon, ntime /) )
      endif 
      deallocate(grid1_array_tmp)

      ncstat = nf_inq_varnatts(nc_infile_id, NF_GLOBAL, nglobatts)
      call netcdf_error_handler(ncstat)
!-----------------------------------------------------------------------

      end subroutine read_input

!***********************************************************************

      end module read_input_file

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!