module clm_mod

  use netcdf

  use kind_mod

  implicit none

  integer, parameter :: namelen = 128   ! string length for paths

  integer, parameter :: varnamelen = 32 ! string length for varnames

  integer, parameter :: posnamelen =  4 ! string length for positive attribute

  integer, parameter :: varattlen = 128 ! string length for variable attributes

  character (len=*), parameter  :: xname  = 'lon'

  character (len=*), parameter  :: yname  = 'lat'

  character (len=*), parameter  :: zsoi_name  = 'levsoi'

  character (len=*), parameter  :: zlak_name  = 'levlak'

  character (len=*), parameter  :: tname  = 'time'

  character (len=*), parameter  :: hr_name  = 'hour'

  character (len=*), parameter  :: nsamples_name  = 'nsamples'

  character (len=*), parameter  :: hiname = 'hist_interval'

  character (len=*), parameter  :: tbname = 'time_bounds'

  character (len=*), parameter  :: aname  = 'area'

  character (len=*), parameter  :: lname  = 'landfrac'

  character (len=*), parameter  :: mname  = 'landmask'

  character (len=*), parameter  :: oname  = 'topo'

  character (len=*), parameter  :: long_name_name = 'long_name'

  character (len=*), parameter  :: missing_value_name = 'missing_value'

  character (len=*), parameter  :: units_name = 'units'

  character (len=*), parameter  :: case_id_name = 'case_id'

  character (len=*), parameter  :: inidat_name = 'Initial_conditions_dataset'

  character (len=*), parameter  :: surdat_name = 'Surface_dataset'

  character (len=*), parameter  :: pftdat_name = 'PFT_physiological_constants_dataset'

  character (len=*), parameter  :: rtmdat_name = 'RTM_input_datset'

  type surface

    integer  :: xsize     ! number of longitude points (lon)

    integer  :: ysize     ! number of latitutde points (lat)

    integer  :: psize     ! number of pfts

    real(r8), pointer :: x(:)    ! longitude values (lon)

    real(r8), pointer :: y(:)    ! latitude values (lat)

    real(r8), pointer :: p(:)    ! pft number

    real(r4), pointer :: pft(:,:,:) ! percent of each pft

  end type surface

  type coordinates

    integer  :: xsize     ! number of longitude points (lon)

    integer  :: ysize     ! number of latitutde points (lat)

    integer  :: zsoi_size ! number of soil z-levels (levsoi)

    integer  :: zlak_size ! number of lake z-levels (levlak)

    integer  :: tsize     ! number of time points (time)

    integer  :: hr_size   ! number of hour points (hour = 24)

    integer  :: nsamples  ! number of time samples in post-processing statistics

    real(r8), pointer :: x(:)    ! longitude values (lon)

    real(r8), pointer :: y(:)    ! latitude values (lat)

    real(r8), pointer :: zsoi(:) ! soil z-level values (levsoi)

    real(r8), pointer :: zlak(:) ! lake z-level values (levlak)

    real(r8), pointer :: hr(:)   ! hour-of-day values (hour)

    real(r8), pointer :: area(:,:)     ! gridcell areas

    real(r8), pointer :: landfrac(:,:) ! gridcell land fractions

    integer, pointer  :: landmask(:,:) ! gridcell land/ocean mask

    real(r8), pointer :: orog(:,:)     ! gridcell orography/topography

    character (len=namelen) :: case_id ! case name of run

    character (len=namelen) :: inidat  ! initial dataset filename

    character (len=namelen) :: surdat  ! surface dataset filename

    character (len=namelen) :: pftdat  ! pft dataset filename

    character (len=namelen) :: rtmdat  ! runoff dataset filename

    logical  :: hour_flag ! flag for hour-of-day statistical summary files

  end type coordinates

  type variable_data

    character (len=varnamelen) :: varname

    integer           :: ndims

    integer, pointer  :: int2d(:,:)

    real(r4), pointer :: var2d(:,:)

    real(r4), pointer :: var3d(:,:,:)

    real(r4), pointer :: var4d(:,:,:,:)

    real(r8), pointer :: time(:)

    real(r8), pointer :: time_bounds(:,:)

    logical  :: int_type ! set if type is integer

    logical  :: soil_layer_flag

    logical  :: hour_flag ! flag for hour-of-day statistical summary variable

    character (len=varattlen)  :: long_name

    real(r4)                   :: missing_value

    character (len=varattlen)  :: units

    character (len=posnamelen) :: positive

    character (len=varnamelen) :: out_varname

  end type variable_data

  type(surface), target       :: surf

  type(coordinates), target   :: coord

  type(variable_data), target :: var_data

  save

  public  clm_read_surf

  public  clm_free_surf

  public  clm_read_coord

  public  clm_copy_grid_data

  public  clm_read_data

  public  clm_convert_data

  public  clm_free_data

  private nc_err

  contains

    !

    !--------------------------------------------------------------------------

    !

    subroutine clm_read_surf(fsurdat)

    implicit none

    character (len=namelen), intent(in) :: fsurdat ! input surface dataset

    integer :: i

    integer :: ierr

    integer :: ncid

    integer :: xdimid, xvarid, xlen

    integer :: ydimid, yvarid, ylen

    integer :: pdimid, pvarid, plen

    real(r8), allocatable :: longxy(:,:)

    real(r8), allocatable :: latixy(:,:)

    real(r4), allocatable :: pctpft(:,:,:)

    ! Read coordinated from surface dataset

    print *, 'Opening ',trim(fsurdat)

    call nc_err(nf90_open(path=trim(fsurdat), mode=nf90_nowrite, ncid=ncid))

    call nc_err(nf90_inq_dimid(ncid, 'lsmlon', xdimid))

    call nc_err(nf90_inq_dimid(ncid, 'lsmlat', ydimid))

    call nc_err(nf90_inq_dimid(ncid, 'lsmpft', pdimid))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=xdimid, len=xlen))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=ydimid, len=ylen))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=pdimid, len=plen))

    call nc_err(nf90_inq_varid(ncid, 'LONGXY', xvarid))

    call nc_err(nf90_inq_varid(ncid, 'LATIXY', yvarid))

    call nc_err(nf90_inq_varid(ncid, 'PCT_PFT', pvarid))

    allocate(longxy(xlen, ylen), latixy(xlen, ylen), pctpft(xlen, ylen, plen), stat=ierr)

    if (ierr /= 0) then

       print *, 'clm_read_surf: Cannot allocate longxy, latixy, pctpft variables'

       stop

    end if

    call nc_err(nf90_get_var(ncid, xvarid, longxy))

    call nc_err(nf90_get_var(ncid, yvarid, latixy))

    call nc_err(nf90_get_var(ncid, pvarid, pctpft))

    print *, 'Closing ',trim(fsurdat)

    call nc_err(nf90_close(ncid))

    ! Allocate and fill in the surface derived data type

    surf%xsize = xlen

    surf%ysize = ylen

    surf%psize = plen

    allocate(surf%x(xlen), surf%y(ylen), surf%p(plen), surf%pft(xlen, ylen, plen), stat=ierr)

    if (ierr /= 0) then

       print *, 'clm_read_surf: Cannot allocate surface type variables'

       stop

    end if

    surf%x(:) = longxy(:,1)

    surf%y(:) = latixy(1,:)

    do i = 1, plen

       surf%p(i) = real(i)

    end do

    surf%pft  = pctpft

    deallocate(longxy, latixy, pctpft)

    end subroutine clm_read_surf

    !

    !--------------------------------------------------------------------------

    !

    subroutine clm_free_surf()

    implicit none

    if (allocated(surf%x))   deallocate(surf%x)

    if (allocated(surf%y))   deallocate(surf%y)

    if (allocated(surf%p))   deallocate(surf%p)

    if (allocated(surf%pft)) deallocate(surf%pft)

    end subroutine clm_free_surf

    !

    !--------------------------------------------------------------------------

    !

    subroutine clm_read_coord(numf, input_path, fnames, read_grid_all)

    implicit none

    integer, intent(in)                 :: numf          ! number of filenames

    character (len=namelen), intent(in) :: input_path    ! input path

    character (len=namelen), intent(in) :: fnames(numf)  ! input filenames

    logical, intent(in)                 :: read_grid_all ! flag to read area,...

    integer :: ncid

    integer :: xdimid, xvarid, xlen

    integer :: ydimid, yvarid, ylen

    integer :: zsoi_dimid, zsoi_varid, zsoi_len

    integer :: zlak_dimid, zlak_varid, zlak_len

    integer :: tdimid, tvarid, tlen

    integer :: hr_dimid, hr_varid, hr_len

    integer :: nsamples_dimid, nsamples_len

    integer :: avarid ! area (when read_grid_all true)

    integer :: lvarid ! landfrac (when read_grid_all true)

    integer :: mvarid ! landmask (when read_grid_all true)

    integer :: ovarid ! orog/topo (when read_grid_all true)

    integer(i8) :: i

    integer(i8) :: time_size

    real(r4), allocatable :: xvals(:)

    real(r4), allocatable :: yvals(:)

    real(r4), allocatable :: zsoi_vals(:)

    real(r4), allocatable :: zlak_vals(:)

    real(r4), allocatable :: hr_vals(:)

    real(r4), allocatable :: avals(:,:) ! area (when read_grid_all true)

    real(r4), allocatable :: lvals(:,:) ! landfrac (when read_grid_all true)

    integer, allocatable  :: mvals(:,:) ! landmask (when read_grid_all true)

    real(r4), allocatable :: ovals(:,:) ! orog/topo (when read_grid_all true)

    character (len=namelen) :: case_id = ''

    character (len=namelen) :: inidat = ''

    character (len=namelen) :: surdat = ''

    character (len=namelen) :: pftdat = ''

    character (len=namelen) :: rtmdat = ''

    logical :: hrflag = .false.

    integer :: ierr

    ! Read coordinates from first file

    print *, 'Opening ',trim(input_path)//'/'//trim(fnames(1))

    call nc_err(nf90_open(path=trim(input_path)//'/'//trim(fnames(1)), mode=nf90_nowrite, ncid=ncid))

    call nc_err(nf90_inq_dimid(ncid, xname, xdimid))

    call nc_err(nf90_inq_dimid(ncid, yname, ydimid))

    call nc_err(nf90_inq_dimid(ncid, zsoi_name, zsoi_dimid))

    call nc_err(nf90_inq_dimid(ncid, zlak_name, zlak_dimid))

    call nc_err(nf90_inq_dimid(ncid, tname, tdimid))

    if (nf90_inq_dimid(ncid, hr_name, hr_dimid) == nf90_noerr) then

       call nc_err(nf90_inq_dimid(ncid, nsamples_name, nsamples_dimid))

       hrflag = .true.

    end if

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=xdimid, len=xlen))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=ydimid, len=ylen))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=zsoi_dimid, len=zsoi_len))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=zlak_dimid, len=zlak_len))

    call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=tdimid, len=tlen))

    if (hrflag) then

       call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=hr_dimid, len=hr_len))

       call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=nsamples_dimid, len=nsamples_len))

    else

       hr_len = 0

       nsamples_len = 0

    end if

    print *, tlen,' values'

    allocate(xvals(xlen), yvals(ylen), zsoi_vals(zsoi_len), zlak_vals(zlak_len), hr_vals(hr_len), stat=ierr)

    if (ierr /= 0) then

       print *, 'clm_read_coord: Cannot allocate grid variables'

       stop

    end if

    if (read_grid_all) then

       allocate(avals(xlen, ylen), lvals(xlen, ylen), mvals(xlen, ylen), ovals(xlen, ylen), stat=ierr)

       if (ierr /= 0) then

          print *, 'clm_read_coord: Cannot allocate extra grid variables'

          stop

       end if

    end if

    call nc_err(nf90_inq_varid(ncid, xname, xvarid))

    call nc_err(nf90_inq_varid(ncid, yname, yvarid))

    call nc_err(nf90_inq_varid(ncid, zsoi_name, zsoi_varid))

    call nc_err(nf90_inq_varid(ncid, zlak_name, zlak_varid))

    call nc_err(nf90_inq_varid(ncid, tname, tvarid))

    if (hrflag) call nc_err(nf90_inq_varid(ncid, hr_name, hr_varid))

    if (read_grid_all) then

       call nc_err(nf90_inq_varid(ncid, aname, avarid))

       call nc_err(nf90_inq_varid(ncid, lname, lvarid))

       call nc_err(nf90_inq_varid(ncid, mname, mvarid))

       call nc_err(nf90_inq_varid(ncid, oname, ovarid))

    end if

    call nc_err(nf90_get_var(ncid, xvarid, xvals))

    call nc_err(nf90_get_var(ncid, yvarid, yvals))

    call nc_err(nf90_get_var(ncid, zsoi_varid, zsoi_vals))

    call nc_err(nf90_get_var(ncid, zlak_varid, zlak_vals))

    if (hrflag) call nc_err(nf90_get_var(ncid, hr_varid, hr_vals))

    if (read_grid_all) then

       call nc_err(nf90_get_var(ncid, avarid, avals))

       call nc_err(nf90_get_var(ncid, lvarid, lvals))

       call nc_err(nf90_get_var(ncid, mvarid, mvals))

       call nc_err(nf90_get_var(ncid, ovarid, ovals))

    end if

    ! Read global attributes

    call nc_err(nf90_get_att(ncid, nf90_global, case_id_name, case_id))

    call nc_err(nf90_get_att(ncid, nf90_global, inidat_name, inidat))

    call nc_err(nf90_get_att(ncid, nf90_global, surdat_name, surdat))

    call nc_err(nf90_get_att(ncid, nf90_global, pftdat_name, pftdat))

    call nc_err(nf90_get_att(ncid, nf90_global, rtmdat_name, rtmdat))

    call nc_err(nf90_close(ncid))

    time_size = tlen

    ! Loop over all remaining files, reading size of time dimension

    do i = 2, numf

      print *, 'Opening ',trim(input_path)//'/'//trim(fnames(i))

      call nc_err(nf90_open(path=trim(input_path)//'/'//trim(fnames(i)), &

         mode=nf90_nowrite, ncid=ncid))

      call nc_err(nf90_inq_dimid(ncid, tname, tdimid))

      call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=tdimid, len=tlen))

      call nc_err(nf90_close(ncid))

      print *, tlen,' values'

      time_size = time_size + tlen

    end do

    ! Allocate and fill the coordinates derived data type

    coord%xsize = xlen

    coord%ysize = ylen

    coord%zsoi_size = zsoi_len

    coord%zlak_size = zlak_len

    coord%tsize = time_size

    coord%hr_size = hr_len

    coord%nsamples = nsamples_len

    allocate(coord%x(xlen), coord%y(ylen), coord%zsoi(zsoi_len), coord%zlak(zlak_len), coord%hr(hr_len), stat=ierr)

    if (ierr /= 0) then

       print *, 'clm_read_coord: Cannot allocate grid variables'

       stop

    end if

    if (read_grid_all) then

       allocate(coord%area(xlen, ylen), coord%landfrac(xlen, ylen), coord%landmask(xlen, ylen), coord%orog(xlen, ylen), stat=ierr)

       if (ierr /= 0) then

          print *, 'clm_read_coord: Cannot allocate extra grid variables'

          stop

       end if

    end if

    coord%x = xvals

    coord%y = yvals

    coord%zsoi = zsoi_vals

    coord%zlak = zlak_vals

    if (hrflag) coord%hr = hr_vals

    if (read_grid_all) then

       coord%area = avals

       coord%landfrac = lvals

       coord%landmask = mvals

       coord%orog = ovals

    end if

    deallocate(xvals, yvals, zsoi_vals, zlak_vals, hr_vals)

    if (read_grid_all) then

       deallocate(avals, lvals, mvals, ovals)

    end if

    coord%case_id = trim(case_id)

    coord%inidat = trim(inidat)

    coord%surdat = trim(surdat)

    coord%pftdat = trim(pftdat)

    coord%rtmdat = trim(rtmdat)

    coord%hour_flag = hrflag

    end subroutine clm_read_coord

    !

    !--------------------------------------------------------------------------

    !

    subroutine clm_copy_grid_data(varname)

    implicit none

    character (len=varnamelen), intent(in) :: varname      ! variable name

    integer :: ierr

    var_data%varname       = varname

    var_data%ndims         = 2

    var_data%missing_value = 1.0e+36

    if (varname == 'landmask') then

       var_data%int_type = .true.

       allocate(var_data%int2d(coord%xsize, coord%ysize), stat=ierr)

       if (ierr /= 0) then

         print *, 'clm_copy_grid_data: Cannot allocate 2d integer variable'

         stop

       end if

    else

       var_data%int_type = .false.

       allocate(var_data%var2d(coord%xsize, coord%ysize), stat=ierr)

       if (ierr /= 0) then

         print *, 'clm_copy_grid_data: Cannot allocate 2d real variable'

         stop

       end if

    end if

    select case (varname)

       case ('area')

          var_data%var2d = coord%area

       case ('landfrac')

          var_data%var2d = coord%landfrac

       case ('landmask')

          var_data%int2d = coord%landmask

       case ('topo')

          var_data%var2d = coord%orog

       case default

          print *,'clm_copy_grid_data: ',varname,' is not a known grid variable'

          stop

    end select

    end subroutine clm_copy_grid_data

    !

    !--------------------------------------------------------------------------

    !

    subroutine clm_read_data(numf, input_path, fnames, tshift, varname)

    implicit none

    integer, intent(in)                    :: numf         ! number of filenames

    character (len=namelen), intent(in)    :: input_path   ! input path

    character (len=namelen), intent(in)    :: fnames(numf) ! input filenames

    integer, intent(in)                    :: tshift(numf) ! time shift by file

    character (len=varnamelen), intent(in) :: varname      ! variable name

    integer :: ncid

    integer :: tdimid, tvarid, tlen

    integer :: hidimid, hilen

    integer :: tbvarid

    integer :: varid, dimids(nf90_max_var_dims), dimlen

    integer :: zsize

    integer(i8) :: i, j, k, l

    integer(i8) :: time_size

    real(r4), allocatable :: tvals(:)    ! temporary time values

    real(r4), allocatable :: tbvals(:,:) ! temporary time bounds values

    real(r4), allocatable :: vals(:,:)   ! temporary data values

    character (len=varnamelen) :: dimname      ! dimension name

    integer :: ierr

    ! Figure out the dimensionality of the variable and variable attributes

    ! from the first file

    print *, 'Opening ',trim(input_path)//'/'//trim(fnames(1)), &

       ' to check dimensionality of ',trim(varname)

    call nc_err(nf90_open(path=trim(input_path)//'/'//trim(fnames(1)), mode=nf90_nowrite, ncid=ncid))

    call nc_err(nf90_inq_varid(ncid, trim(varname), varid))

    call nc_err(nf90_inquire_variable(ncid, varid, name=var_data%varname, ndims=var_data%ndims, dimids=dimids))

    ! Figure out if there is a Z-dimension and/or an hour-dimension

    ierr = 0

    var_data%int_type = .false.

    var_data%hour_flag = .false.

    var_data%soil_layer_flag = .false.

    do i = 1, var_data%ndims

       call nc_err(nf90_inquire_dimension(ncid, dimids(i), name=dimname, len=dimlen))

       if (i == 1 .and. dimname /= xname) ierr = ierr + 1

       if (i == 2 .and. dimname /= yname) ierr = ierr + 1

       if (i > 2 .and. i < var_data%ndims) then

          if (dimname == hr_name) then

             var_data%hour_flag = .true.

          else if (dimname == zsoi_name) then

             var_data%soil_layer_flag = .true.

          !else if (dimname == zlak_name) then

          !   var_data%zdim_name = dimname

          else

             ierr = ierr + 1

          end if

       end if

       ! This is not a good assumption (that time is always there), so skip it

       !if (i == var_data%ndims .and. dimname /= tname) ierr = ierr + 1

       print *,'Dimension ',i,' is ',dimname

    end do

    if (ierr /= 0) then

       print *,'clm_read_data: Unexpected dimension combination'

       stop

    end if

    !

    call nc_err(nf90_get_att(ncid, varid, long_name_name, var_data%long_name))

    call nc_err(nf90_get_att(ncid, varid, missing_value_name, var_data%missing_value))

    call nc_err(nf90_get_att(ncid, varid, units_name, var_data%units))

    call nc_err(nf90_close(ncid))

    ! Allocate space for the variable of the correct dimensionality

    select case (var_data%ndims)

      case (2)

        allocate(var_data%var2d(coord%xsize, coord%ysize), stat=ierr)

        if (ierr /= 0) then

          print *, 'clm_read_data: Cannot allocate 2d variable'

          stop

        end if

      case (3)

        allocate(var_data%var3d(coord%xsize, coord%ysize, coord%tsize), stat=ierr)

        if (ierr /= 0) then

          print *, 'clm_read_data: Cannot allocate 3d variable'

          stop

        end if

      case (4)

        if (var_data%hour_flag) then

          zsize = coord%hr_size

        else if (var_data%soil_layer_flag) then

          zsize = coord%zsoi_size

        !else if (var_data%zdim_name == zlak_name) then

        !   zsize = coord%zlak_size

        else

          print *,'clm_read_data: Unexpected fourth dimension while allocating'

          stop

        end if

        allocate(var_data%var4d(coord%xsize, coord%ysize, zsize, coord%tsize), stat=ierr)

        if (ierr /= 0) then

          print *, 'clm_read_data: Cannot allocate 4d variable'

          stop

        end if

      case default

        print *, 'clm_read_data: Cannot handle variables of ', var_data%ndims, ' dimensions'

        stop

    end select

    ! Allocate space for time stamps and bounds

    allocate(var_data%time(coord%tsize), stat=ierr)

    if (ierr /= 0) then

      print *, 'clm_read_data: Cannot allocate time variable'

      stop

    end if

    allocate(var_data%time_bounds(2,coord%tsize), stat=ierr)

    if (ierr /= 0) then

      print *, 'clm_read_data: Cannot allocate time bounds variable'

      stop

    end if

    ! Allocate temporary space for the variable level-slice/time-slice */

    allocate(vals(coord%xsize, coord%ysize), stat=ierr)

    if (ierr /= 0) then

      print *, 'clm_read_data: Cannot allocate variable time-slice'

      stop

    end if

    if (var_data%ndims == 2) then

      ! No time axis, so just read the first instance from the first file

      print *, 'Opening ',trim(input_path)//'/'//trim(fnames(1))

      call nc_err(nf90_open(path=trim(input_path)//'/'//trim(fnames(1)), mode=nf90_nowrite, ncid=ncid))

      call nc_err(nf90_inq_varid(ncid, trim(varname), varid))

      call nc_err(nf90_get_var(ncid, varid, vals, start = (/ 1, 1 /), count = (/ coord%xsize, coord%ysize /)))

      var_data%var2d(:,:) = vals(:,:)

      call nc_err(nf90_close(ncid))

    else

      ! Otherwise, loop over all times

      time_size = 0

      ! Loop over all files, reading data

      do i = 1, numf

        print *, 'Opening ',trim(input_path)//'/'//trim(fnames(i))

        call nc_err(nf90_open(path=trim(input_path)//'/'//trim(fnames(i)), mode=nf90_nowrite, ncid=ncid))

        call nc_err(nf90_inq_dimid(ncid, tname, tdimid))

        call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=tdimid, len=tlen))

        call nc_err(nf90_inq_dimid(ncid, hiname, hidimid))

        call nc_err(nf90_inquire_dimension(ncid=ncid, dimid=hidimid, len=hilen))

        call nc_err(nf90_inq_varid(ncid, tname, tvarid))

        call nc_err(nf90_inq_varid(ncid, tbname, tbvarid))

        call nc_err(nf90_inq_varid(ncid, trim(varname), varid))

        ! Allocate temporary space for the time stamps and bounds

        allocate(tvals(tlen), stat=ierr)

        if (ierr /= 0) then

          print *, 'clm_read_data: Cannot allocate temporary time variable'

          stop

        end if

        allocate(tbvals(hilen,tlen), stat=ierr)

        if (ierr /= 0) then

          print *, 'clm_read_data: Cannot allocate temporary time bounds variable'

          stop

        end if

        ! Get the time stamps and bounds

        call nc_err(nf90_get_var(ncid, tvarid, tvals))

        call nc_err(nf90_get_var(ncid, tbvarid, tbvals))

        do j = 1, tlen

          ! Get and store a slab of data

          select case (var_data%ndims)

            case (3)

              !print *, 'Reading slice of data'

              call nc_err(nf90_get_var(ncid, varid, vals, start = (/ 1, 1, j /), count = (/ coord%xsize, coord%ysize, 1 /)))

              var_data%var3d(:,:,time_size + j) = vals(:,:)

            case (4)

              do k = 1, zsize

                 call nc_err(nf90_get_var(ncid, varid, vals, start = (/ 1, 1, k, j /), count = (/ coord%xsize, coord%ysize, 1, 1 /)))

                 var_data%var4d(:,:,k,time_size + j) = vals(:,:)

              end do

            case default

              print *, 'clm_read_data: Cannot handle variables of ', var_data%ndims, ' dimensions'

              stop

          end select

          ! Store each time stamp and bounds, including the time offset

          var_data%time(time_size + j) = tvals(j) + tshift(i)

          var_data%time_bounds(:,time_size + j) = tbvals(:,j) + tshift(i)

        end do

        call nc_err(nf90_close(ncid))

        deallocate(tvals)

        deallocate(tbvals)

        time_size = time_size + tlen

      end do

    end if

    deallocate(vals)

    end subroutine clm_read_data

    !

    !--------------------------------------------------------------------------

    !

    subroutine clm_convert_data(casa_flux_bug)

    implicit none

    logical, intent(in)  :: casa_flux_bug

    real(r8)             :: sfactor

    real(r8)             :: offset

    integer              :: ierr

    logical, allocatable :: mask2d(:,:)

    logical, allocatable :: mask3d(:,:,:)

    logical, allocatable :: mask4d(:,:,:,:)

    ! Default offset is zero

    offset = 0.

    ! Default var_data structure values

    var_data%positive = ''

    ! Set scale factor according to which variable is to be converted

    select case (var_data%varname)

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

      ! Atmospheric forcing

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

      ! husf

      case ('QBOT')

        ! QBOT: kg kg-1 -> kg kg-1

        sfactor = 1.0

        var_data%units = 'kg kg-1'

        var_data%out_varname = 'husf'

      ! prra

      case ('RAIN')

        ! RAIN: mm s-1 -> kg m-2 s-1

        sfactor = 1.0

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'prra'

      ! prsn

      case ('SNOW')

        ! SNOW: mm s-1 -> kg m-2 s-1

        sfactor = 1.0

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'prsn'

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

      ! Biogeochemistry

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

      ! agbc

      ! aglbc

      ! agnpp

      case ('AGNPP')

        ! AGNPP: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'agnpp'

      ! ar

      case ('AR')

        ! AR: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'ar'

      ! bco

      case ('BIOGENCO')

        ! BIOGENCO: ug m-2 h-1 -> kg m-2 s-1

        sfactor = 1.0e-9 * (1.0 / 3600.)

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'bco'

      ! bgbc

      ! bgnpp

      case ('BGNPP')

        ! BGNPP: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'bgnpp'

      ! bipn

      case ('ISOPRENE')

        ! ISOPRENE: ug m-2 h-1 -> kg m-2 s-1

        sfactor = 1.0e-9 * (1.0 / 3600.)

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'bipn'

      ! bmtpn

      case ('MONOTERP')

        ! MONOTERP: ug m-2 h-1 -> kg m-2 s-1

        sfactor = 1.0e-9 * (1.0 / 3600.)

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'bmtpn'

      ! borvoc

      case ('ORVOC')

        ! ORVOC: ug m-2 h-1 -> kg m-2 s-1

        sfactor = 1.0e-9 * (1.0 / 3600.)

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'borvoc'

      ! bovoc

      case ('OVOC')

        ! OVOC: ug m-2 h-1 -> kg m-2 s-1

        sfactor = 1.0e-9 * (1.0 / 3600.)

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'bovoc'

      ! cwdc

      case ('CWDC')

        ! CWDC: g m-2 -> kg m-2

        sfactor = 1.0e-3

        var_data%units = 'kg m-2'

        var_data%out_varname = 'cwdc'

      ! cwdc_hr

      case ('CWDC_HR')

        ! CWDC_HR: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'cwdc_hr'

      ! cwdc_loss

      case ('CWDC_LOSS')

        ! CWDC_LOSS: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'cwdc_loss'

      ! froot_litterc

      ! froot_litterc_hr

      ! froot_litterc_loss

      ! frootc

      case ('FROOTC')

        ! FROOTC: g m-2 -> kg m-2

        sfactor = 1.0e-3

        var_data%units = 'kg m-2'

        var_data%out_varname = 'frootc'

      ! frootc_alloc

      case ('FROOTC_ALLOC')

        ! FROOTC_ALLOC: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'frootc_alloc'

      ! frootc_loss

      case ('FROOTC_LOSS')

        ! FROOTC_LOSS: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'frootc_loss'

      ! gpp

      case ('GPP')

        ! GPP: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'gpp'

      ! gnmin

      case ('GROSS_NMIN')

        ! GROSS_NMIN: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'gnmin'

      ! hr

      case ('HR')

        ! HR: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'hr'

      ! leaf_litterc

      ! leaf_litterc_hr

      ! leaf_litterc_loss

      ! leafc

      case ('LEAFC')

        ! LEAFC: g m-2 -> kg m-2

        sfactor = 1.0e-3

        var_data%units = 'kg m-2'

        var_data%out_varname = 'leafc'

      ! leafc_alloc

      case ('LEAFC_ALLOC')

        ! LEAFC_ALLOC: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'leafc_alloc'

      ! leafc_loss

      case ('LEAFC_LOSS')

        ! LEAFC_LOSS: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'leafc_loss'

      ! litterc

      case ('LITTERC')

        ! LITTERC: g m-2 -> kg m-2

        sfactor = 1.0e-3

        var_data%units = 'kg m-2'

        var_data%out_varname = 'litterc'

      ! litterc_hr

      case ('LITTERC_HR')

        ! LITTERC_HR: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'litterc_hr'

      ! litterc_loss

      case ('LITTERC_LOSS')

        ! LITTERC_LOSS: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'litterc_loss'

      ! nee

      case ('NEE')

        ! NEE: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'nee'

      ! nep

      case ('NEP')

        ! NEP: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'nep'

      ! nnmin

      case ('NET_NMIN')

        ! NEP: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'nnmin'

      ! npp

      case ('NPP')

        ! NPP: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'npp'

      ! soilc

      case ('SOILC')

        ! SOILC: g m-2 -> kg m-2

        sfactor = 1.0e-3

        var_data%units = 'kg m-2'

        var_data%out_varname = 'soilc'

      ! soilc_hr

      case ('SOILC_HR')

        ! SOILC_HR: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'soilc_hr'

      ! soilc_loss

      case ('SOILC_LOSS')

        ! SOILC_LOSS: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'soilc_loss'

      ! tbvoc

      case ('VOCFLXT')

        ! VOCFLXT: ug m-2 h-1 -> kg m-2 s-1

        sfactor = 1.0e-9 * (1.0 / 3600.)

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'tbvoc'

      ! woodc

      case ('WOODC')

        ! WOODC: g m-2 -> kg m-2

        sfactor = 1.0e-3

        var_data%units = 'kg m-2'

        var_data%out_varname = 'woodc'

      ! woodc_alloc

      case ('WOODC_ALLOC')

        ! WOODC_ALLOC: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'woodc_alloc'

      ! woodc_loss

      case ('WOODC_LOSS')

        ! WOODC_LOSS: g m-2 s-1 -> kg m-2 s-1

        sfactor = 1.0e-3

        if (casa_flux_bug) then

           print *,'WARNING: casa_flux_bug set to true; dividing data by 1200'

           sfactor = sfactor * (1. / 1200.)

        end if

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'woodc_loss'

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

      ! Energy

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

      ! hfg

      case ('FGR')

        ! FGR: W m-2 -> W m-2

        sfactor = 1.0

        var_data%units = 'W m-2'

        var_data%positive = 'down'

        var_data%out_varname = 'hfg'

      ! hfls

      case ('LATENT')

        ! LATENT: W m-2 -> W m-2

        sfactor = 1.0

        var_data%units = 'W m-2'

        var_data%positive = 'up'

        var_data%out_varname = 'hfls'

      ! hfss

      case ('FSH')

        ! FSH: W m-2 -> W m-2

        sfactor = 1.0

        var_data%units = 'W m-2'

        var_data%positive = 'up'

        var_data%out_varname = 'hfss'

      ! hfssg

      case ('FSH_G')

        ! FSH_G: W m-2 -> W m-2

        sfactor = 1.0

        var_data%units = 'W m-2'

        var_data%positive = 'up'

        var_data%out_varname = 'hfssg'

      ! hfssv

      case ('FSH_V')

        ! FSH_V: W m-2 -> W m-2

        sfactor = 1.0

        var_data%units = 'W m-2'

        var_data%positive = 'up'

        var_data%out_varname = 'hfssv'

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

      ! Grid

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

      ! gca

      case ('area')

        ! area: km2 -> m2

        sfactor = 1.0e+3

        var_data%units = 'm2'

        var_data%out_varname = 'gca'

      ! lbm

      case ('landmask')

        ! landmask: 1 -> 1

        var_data%units = '1'

        var_data%out_varname = 'lbm'

        ! Do not convert integer landmask, just return

        return

      ! lcf

      ! lcn

      ! orog

      case ('topo')

        ! topo: m -> m

        sfactor = 1.0

        var_data%units = 'm'

        var_data%out_varname = 'orog'

      ! sftlf

      case ('landfrac')

        ! landfrac: 1 -> 1

        sfactor = 1.0

        var_data%units = '1'

        var_data%out_varname = 'sftlf'

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

      ! Humidity

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

      ! hur2m

      case ('RH2M')

        ! RH2M: % -> %

        sfactor = 1.0

        var_data%units = '%'

        var_data%out_varname = 'hur2m'

      ! hus2m

      case ('Q2M')

        ! Q2M: kg kg-1 -> kg kg-1

        sfactor = 1.0

        var_data%units = 'kg kg-1'

        var_data%out_varname = 'hus2m'

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

      ! Hydrology

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

      ! btran

      case ('BTRAN')

        ! BTRAN: 1 -> 1

        sfactor = 1.0

        var_data%units = '1'

        var_data%out_varname = 'btran'

      ! et

      case ('ET')

        ! ET: m s-1 -> kg m-2 s-1

        sfactor = 1.0e+3

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'et'

      ! evspsblsoi

      case ('QSOIL')

        ! QSOIL: mm s-1 -> kg m-2 s-1

        sfactor = 1.0

        var_data%units = 'kg m-2 s-1'

        var_data%out_varname = 'evspsblsoi'