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'
      ! evspsblveg
      case ('QVEGE')
        ! QVEGE: mm s-1 -> kg m-2 s-1
        sfactor = 1.0
        var_data%units = 'kg m-2 s-1'
        var_data%out_varname = 'evspsblveg'
      ! mrfsl (multi-level)
      case ('SOILICE')
        ! SOILICE: kg m-2 -> kg m-2
        sfactor = 1.0
        var_data%units = 'kg m-2'
        var_data%out_varname = 'mrfsl'
      ! mrfso
      case ('TOTAL_SOILICE')
        ! TOTAL_SOILICE: kg m-2 -> kg m-2
        sfactor = 1.0
        var_data%units = 'kg m-2'
        var_data%out_varname = 'mrfso'
      ! mrlsl (multi-level)
      case ('SOILLIQ')
        ! SOILLIQ: kg m-2 -> kg m-2
        sfactor = 1.0
        var_data%units = 'kg m-2'
        var_data%out_varname = 'mrlsl'
      ! mrlso
      case ('TOTAL_SOILLIQ')
        ! TOTAL_SOILLIQ: kg m-2 -> kg m-2
        sfactor = 1.0
        var_data%units = 'kg m-2'
        var_data%out_varname = 'mrlso'
      ! mrros
      case ('QOVER')
        ! QOVER: mm s-1 -> kg m-2 s-1
        sfactor = 1.0
        var_data%units = 'kg m-2 s-1'
        var_data%out_varname = 'mrros'
      ! mrross
      case ('QDRAI')
        ! QDRAI: mm s-1 -> kg m-2 s-1
        sfactor = 1.0
        var_data%units = 'kg m-2 s-1'
        var_data%out_varname = 'mrross'
      ! prveg
      case ('QINTR')
        ! QINTR: mm s-1 -> kg m-2 s-1
        sfactor = 1.0
        var_data%units = 'kg m-2 s-1'
        var_data%out_varname = 'prveg'
      ! snd
      case ('H2OSNO')
        ! H2OSNO: mm -> m
        sfactor = 1.0e-3
        var_data%units = 'm'
        var_data%out_varname = 'snd'
      ! tran
      case ('QVEGT')
        ! QVEGT: mm s-1 -> kg m-2 s-1
        sfactor = 1.0
        var_data%units = 'kg m-2 s-1'
        var_data%out_varname = 'tran'
      ! vfwsl (multi-level)
      case ('H2OSOI')
        ! H2OSOI: mm3 mm-3 -> m3 m-3
        sfactor = 1.0
        var_data%units = 'm3 m-3'
        var_data%out_varname = 'vfwsl'
      ! wpsl (multi-level)
      case ('SOILPSI')
        ! SOILPSI: Pa -> Pa
        sfactor = 1.0
        var_data%units = 'Pa'
        var_data%out_varname = 'wpsl'

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! Radiation
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! alsas
      case ('ALL_SKY_ALBEDO')
        ! ALL_SKY_ALBEDO: 1 -> 1
        sfactor = 1.0
        var_data%units = '1'
        var_data%out_varname = 'alsas'
      ! alsbs
      case ('BLACK_SKY_ALBEDO')
        ! BLACK_SKY_ALBEDO: 1 -> 1
        sfactor = 1.0
        var_data%units = '1'
        var_data%out_varname = 'alsbs'
      ! laisha
      case ('LAISHA')
        ! LAISHA: m2 m-2 -> m2 m-2
        sfactor = 1.0
        var_data%units = 'm2 m-2'
        var_data%out_varname = 'laisha'
      ! laisun
      case ('LAISUN')
        ! LAISUN: m2 m-2 -> m2 m-2
        sfactor = 1.0
        var_data%units = 'm2 m-2'
        var_data%out_varname = 'laisun'
      ! rlds
      case ('FLDS')
        ! FLDS: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rlds'
      ! rls
      case ('FIRA')
        ! FIRA: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rls'
      ! rlus
      case ('FIRE')
        ! FIRE: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'up'
        var_data%out_varname = 'rlus'
      ! rs
      case ('NETRAD')
        ! NETRAD: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rs'
      ! rsas
      case ('FSA')
        ! FSA: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rsas'
      ! rsasg
      case ('SABG')
        ! SABG: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rsasg'
      ! rsasv
      case ('SABV')
        ! SABV: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rsasv'
      ! rsds
      case ('FSDS')
        ! FSDS: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'down'
        var_data%out_varname = 'rsds'
      ! rsus
      case ('FSR')
        ! FSR: W m-2 -> W m-2
        sfactor = 1.0
        var_data%units = 'W m-2'
        var_data%positive = 'up'
        var_data%out_varname = 'rsus'
      ! snc
      case ('FSNO')
        ! FSNO: 1 -> 1
        sfactor = 1.0
        var_data%units = '1'
        var_data%out_varname = 'snc'

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! Temperature
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! degday
      case ('DEGDAY')
        ! DEGDAY: C -> K
        sfactor = 1.0
        offset  = 273.15
        var_data%units = 'K'
        var_data%out_varname = 'degday'
      ! ta2m
      case ('TSA')
        ! TSA: K -> K
        sfactor = 1.0
        var_data%units = 'K'
        var_data%out_varname = 'ta2m'
      ! ta2mdmax
      case ('TREFMXAV')
        ! TREFMXAV: K -> K
        sfactor = 1.0
        var_data%units = 'K'
        var_data%out_varname = 'ta2mdmax'
      ! ta2mdmin
      case ('TREFMNAV')
        ! TREFMNAV: K -> K
        sfactor = 1.0
        var_data%units = 'K'
        var_data%out_varname = 'ta2mdmin'
      ! tsl (multi-level)
      case ('TSOI')
        ! TSOI: K -> K
        sfactor = 1.0
        var_data%units = 'K'
        var_data%out_varname = 'tsl'
      ! tv
      case ('TV')
        ! TV: K -> K
        sfactor = 1.0
        var_data%units = 'K'
        var_data%out_varname = 'tv'

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! Vegetation phenology
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! laiexp
      case ('ELAI')
        ! ELAI: m2 m-2 -> m2 m-2
        sfactor = 1.0
        var_data%units = 'm2 m-2'
        var_data%out_varname = 'laiexp'
      ! laitot
      case ('TLAI')
        ! TLAI: m2 m-2 -> m2 m-2
        sfactor = 1.0
        var_data%units = 'm2 m-2'
        var_data%out_varname = 'laitot'
      ! saiexp
      case ('ESAI')
        ! ESAI: m2 m-2 -> m2 m-2
        sfactor = 1.0
        var_data%units = 'm2 m-2'
        var_data%out_varname = 'saiexp'
      ! saitot
      case ('TSAI')
        ! TSAI: m2 m-2 -> m2 m-2
        sfactor = 1.0
        var_data%units = 'm2 m-2'
        var_data%out_varname = 'saitot'

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! Vegetation physiology
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! psn
      case ('FPSN')
        ! FPSN: umol m-2 s-1 -> kg m-2 s-1
        sfactor = 12.0 * 1.0e-6 * 1.0e-3
        var_data%units = 'kg m-2 s-1'
        var_data%out_varname = 'psn'
      ! rssha
      case ('RSSHA')
        ! RSSHA: s m-1 -> s m-1
        sfactor = 1.0
        var_data%units = 's m-1'
        var_data%out_varname = 'rssha'
      ! rssun
      case ('RSSUN')
        ! RSSUN: s m-1 -> s m-1
        sfactor = 1.0
        var_data%units = 's m-1'
        var_data%out_varname = 'rssun'

      case default
        print *, 'clm_convert_data: Unknown variable name ', var_data%varname
        stop
    end select

    ! Convert to the desired units using the offset and scale factor (sfactor)
    select case (var_data%ndims)
      case (2)
        allocate(mask2d(coord%xsize, coord%ysize), stat=ierr)
        if (ierr /= 0) then
          print *, 'clm_convert_data: Cannot allocate 2d mask'
          stop
        end if
        mask2d = .false.
        if (var_data%int_type) then
          where (var_data%int2d /= var_data%missing_value)
            var_data%int2d = int(var_data%int2d * sfactor + offset)
            mask2d = .true.
          end where
          print *, var_data%out_varname, &
            ' min:', minval(var_data%int2d, mask=mask2d), &
            ' max:', maxval(var_data%int2d, mask=mask2d)
        else
          where (var_data%var2d /= var_data%missing_value)
            var_data%var2d = var_data%var2d * sfactor + offset
            mask2d = .true.
          end where
          print *, var_data%out_varname, &
            ' min:', minval(var_data%var2d, mask=mask2d), &
            ' max:', maxval(var_data%var2d, mask=mask2d)
        end if
        if (allocated(mask2d)) deallocate(mask2d)
      case (3)
        allocate(mask3d(coord%xsize, coord%ysize, coord%tsize), stat=ierr)
        if (ierr /= 0) then
          print *, 'clm_convert_data: Cannot allocate 3d mask'
          stop
        end if
        mask3d = .false.
        where (var_data%var3d /= var_data%missing_value)
           var_data%var3d = (var_data%var3d + offset) * sfactor
           mask3d = .true.
        end where
        print *, var_data%out_varname, &
          ' min:', minval(var_data%var3d, mask=mask3d), &
          ' max:', maxval(var_data%var3d, mask=mask3d)
        if (allocated(mask3d)) deallocate(mask3d)
      case (4)
        if (var_data%soil_layer_flag) then
          allocate(mask4d(coord%xsize, coord%ysize, coord%zsoi_size, coord%tsize), stat=ierr)
          if (ierr /= 0) then
            print *, 'clm_convert_data: Cannot allocate 4d mask'
            stop
          end if
          mask4d = .false.
          where (var_data%var4d /= var_data%missing_value)
             var_data%var4d = (var_data%var4d + offset) * sfactor
             mask4d = .true.
          end where
          print *, var_data%out_varname, &
            ' min:', minval(var_data%var4d, mask=mask4d), &
            ' max:', maxval(var_data%var4d, mask=mask4d)
          if (allocated(mask4d)) deallocate(mask4d)
        else
          print *,'clm_convert_data: Unknown data dimension'
          stop
        end if
      case default
        print *,'clm_convert_data: Dimension count not supported'
        stop
    end select

    end subroutine clm_convert_data
    !
    !--------------------------------------------------------------------------
    !
    subroutine clm_free_data

    implicit none

    var_data%varname = ''
    var_data%ndims = 0
    if (allocated(var_data%int2d))        deallocate(var_data%int2d)
    if (allocated(var_data%var2d))        deallocate(var_data%var2d)
    if (allocated(var_data%var3d))        deallocate(var_data%var3d)
    if (allocated(var_data%var4d))        deallocate(var_data%var4d)
    if (allocated(var_data%time))         deallocate(var_data%time)
    if (allocated(var_data%time_bounds))  deallocate(var_data%time_bounds)
    var_data%int_type        = .false.
    var_data%soil_layer_flag = .false.
    var_data%hour_flag       = .false.
    var_data%long_name       = ''
    var_data%missing_value   = 0.
    var_data%units           = ''
    var_data%positive        = ''
    var_data%out_varname     = ''

    end subroutine clm_free_data
    !
    !--------------------------------------------------------------------------
    !
    subroutine nc_err(status)

    implicit none
    integer :: status ! status code

    if (status /= nf90_noerr) then
       print *, nf90_strerror(status)
       stop 'Program stopped after encountering netCDF error.'
    end if

    end subroutine nc_err
    !
    !--------------------------------------------------------------------------
    !

end module clm_mod
