module restUtilMod

  !-----------------------------------------------------------------------
  ! provies generic routines and types for use with restart files
  !
#include "shr_assert.h"
  use shr_kind_mod, only: r8=>shr_kind_r8, r4 => shr_kind_r4, i4=>shr_kind_i4
  use shr_sys_mod,  only: shr_sys_abort
  use spmdMod,      only: masterproc
  use clm_varctl,   only: iulog, nsrest, nsrContinue, nsrBranch
  use clm_varcon,   only: spval, ispval
  use decompMod,    only: bounds_type
  use abortutils,   only: endrun
  use ncdio_pio
  use pio
  use ncdio_utils,  only: find_var_on_file
  use shr_string_mod, only: shr_string_listGetName
  !
  implicit none
  save
  private
  ! save
  !
  !-----------------------------------------------------------------------

  interface restartvar
     !DIMS 0,1,2
     !TYPE text,int,double
     module procedure restartvar_{DIMS}d_{TYPE}
     !TYPE int,double
     module procedure restartvar_2d_{TYPE}_bounds
  end interface restartvar

  ! iflag_interp => interpolate variable
  ! iflag_copy => copy variable
  ! iflag_skip => skip variable: maintain at cold start value on output file
  ! iflag_area => area-related variable: skip if we're taking areas from the surface
  !   dataset, interp if we're taking areas from the input finidat file. (The latter only
  !   works if we're running at the same resolution with a similar configuration.) This
  !   can apply to variables that are not actually areas, but for which we want this same
  !   conditional behavior.
  integer,parameter, public :: iflag_interp = 1
  integer,parameter, public :: iflag_copy   = 2
  integer,parameter, public :: iflag_skip   = 3
  integer,parameter, public :: iflag_area   = 4

  integer,parameter, public :: iflag_noswitchdim = 0
  integer,parameter, public :: iflag_switchdim   = 1

  public :: restartvar

  ! Set values of a missing restart field from a template field, with some constant
  ! multiplier
  public :: set_missing_from_template

  ! Set a gridcell-level field from a column-level field
  public :: set_grc_field_from_col_field

  ! Set any NaN or spval entries in a field to some constant value
  interface set_missing_vals_to_constant
     !DIMS 1,2
     module procedure set_missing_vals_to_constant_{DIMS}d
  end interface set_missing_vals_to_constant
  public :: set_missing_vals_to_constant

  private :: missing_field_possibly_abort
  private :: write_interpinic_flag

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

contains

  !-----------------------------------------------------------------------
  !DIMS 0
  !TYPE text,int,double
  subroutine restartvar_{DIMS}d_{TYPE}(&
       ncid, flag, varname, xtype, &
       long_name, units, interpinic_flag, data, readvar, &
       comment, flag_meanings, missing_value, fill_value, &
       imissing_value, ifill_value, flag_values, nvalid_range )

    ! Note that varname can be a colon-delimited list of possible variable names (with no
    ! spaces around the colons). In this case, when flag = 'read', the input file is
    ! searched for each possible variable name in order, starting with the first, until
    ! one is found. This mechanism supports backwards compatibility with old restart
    ! files, in case variables have been renamed. For example, if variable 'foo' was
    ! recently renamed to 'bar', then varname should be 'bar:foo'. For flag = 'write',
    ! the first name in the list is used.

    !----------------------------------------------------
    ! Arguments
    type(file_desc_t) , intent(inout)        :: ncid             ! netcdf file id
    character(len=*)  , intent(in)           :: flag             ! 'read' or 'write'
    character(len=*)  , intent(in)           :: varname          ! variable name (or colon-delimited list: see above)
    integer           , intent(in)           :: xtype            ! netcdf data type
    character(len=*)  , intent(in)           :: long_name        ! long name for variable
    character(len=*)  , intent(in)           :: interpinic_flag  ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module)
    {VTYPE}           , intent(inout)        :: data{DIMSTR}
    logical           , intent(out)          :: readvar          ! was var read?
    character(len=*)  , intent(in), optional :: units            ! long name for variable
    character(len=*)  , intent(in), optional :: comment          ! attribute
    character(len=*)  , intent(in), optional :: flag_meanings(:) ! attribute
    real(r8)          , intent(in), optional :: missing_value    ! attribute for real
    real(r8)          , intent(in), optional :: fill_value       ! attribute for real
    integer           , intent(in), optional :: imissing_value   ! attribute for int
    integer           , intent(in), optional :: ifill_value      ! attribute for int
    integer           , intent(in), optional :: flag_values(:)   ! attribute for int
    integer           , intent(in), optional :: nvalid_range(2)  ! attribute for int
    ! 
    ! Local variables
    character(len=len(varname)) :: primary_varname ! first name in the varname list
    character(len=len(varname)) :: my_varname      ! actual varname to read/write
    integer          :: ivalue
    type(var_desc_t) :: vardesc  ! local vardesc
    integer          :: status   ! return error code 
    integer          :: varid
    integer          :: lxtype   ! local external type (in case logical variable)
    !----------------------------------------------------

    call shr_string_listGetName(varname, 1, primary_varname)
    if (flag == 'read') then
       call find_var_on_file(ncid, varname, is_dim=.false., varname_on_file=my_varname)
       if ((my_varname /= primary_varname) .and. masterproc) then
          write(iulog,*) 'Restart file backwards compatibility: Translating: ', &
               trim(my_varname), ' => ', trim(primary_varname)
       end if
    else
       my_varname = primary_varname
    end if

    if (flag == 'define') then

       if ( xtype == ncd_log )then
          lxtype = ncd_int
       else
          lxtype = xtype
       end if

       call ncd_defvar(ncid=ncid, varname=trim(my_varname), xtype=lxtype, &
            long_name=trim(long_name), units=units)

       status = PIO_inq_varid(ncid, trim(my_varname), vardesc)
       varid = vardesc%varid

       call write_interpinic_flag(ncid=ncid, varid=varid, varname=varname, &
            interpinic_flag=interpinic_flag)

       ! This attribute is written in order to communicate this metadata to initInterp
       call ncd_putatt(ncid, varid, 'varnames_on_old_files', trim(varname))

       if (present(comment)) then
          call ncd_putatt(ncid, varid, 'comment', trim(comment))
       end if
       if (present(units)) then
          call ncd_putatt(ncid, varid, 'units', trim(units))
       end if
       if (present(fill_value)) then
          call ncd_putatt(ncid, varid, '_FillValue', fill_value, lxtype)
       else if (lxtype == ncd_double) then
          call ncd_putatt(ncid, varid, '_FillValue', spval, lxtype)
       end if
       if (present(missing_value)) then
          call ncd_putatt(ncid, varid, 'missing_value', missing_value, lxtype)
       else if (lxtype == ncd_double) then
          call ncd_putatt(ncid, varid, 'missing_value', spval, lxtype)
       end if
       if (present(ifill_value)) then
          call ncd_putatt(ncid, varid, '_FillValue', ifill_value, lxtype)
       else if (lxtype == ncd_int) then
          call ncd_putatt(ncid, varid, '_FillValue', ispval, lxtype)
       end if
       if (present(imissing_value)) then
          call ncd_putatt(ncid, varid, 'missing_value', imissing_value, lxtype)
       else if (lxtype == ncd_int) then
          call ncd_putatt(ncid, varid, 'missing_value', ispval, lxtype)
       end if
       if ( xtype == ncd_log )then
          status = PIO_put_att(ncid,varid,'flag_values',     (/0, 1/) )
          status = PIO_put_att(ncid,varid,'flag_meanings',  "FALSE TRUE" )
          status = PIO_put_att(ncid,varid,'valid_range',    (/0, 1/) )
       end if

    else if (flag == 'read' .or. flag == 'write') then

#if ({ITYPE}!=TYPETEXT) 
       call ncd_io(varname=trim(my_varname), data=data, &
            ncid=ncid, flag=flag, readvar=readvar)
#endif
    end if

    if (flag == 'read') then
       if (.not. readvar) then
          call missing_field_possibly_abort(my_varname)
       end if
    end if

  end subroutine restartvar_{DIMS}d_{TYPE}

  !-----------------------------------------------------------------------
  !DIMS 1,2
  !TYPE text,int,double
  subroutine restartvar_{DIMS}d_{TYPE}(&
       ncid, flag, varname, xtype, dim1name, dim2name, &
       long_name, units, interpinic_flag, data, readvar, &
       comment, flag_meanings, missing_value, fill_value, &
       imissing_value, ifill_value, flag_values, nvalid_range, is_spatial )

    ! Note that varname can be a colon-delimited list of possible variable names (with no
    ! spaces around the colons). In this case, when flag = 'read', the input file is
    ! searched for each possible variable name in order, starting with the first, until
    ! one is found. This mechanism supports backwards compatibility with old restart
    ! files, in case variables have been renamed. For example, if variable 'foo' was
    ! recently renamed to 'bar', then varname should be 'bar:foo'. For flag = 'write',
    ! the first name in the list is used.

    !----------------------------------------------------
    ! Arguments
    type(file_desc_t) , intent(inout)        :: ncid             ! netcdf file id
    character(len=*)  , intent(in)           :: flag             ! 'read' or 'write'
    character(len=*)  , intent(in)           :: varname          ! variable name (or colon-delimited list: see above)
    integer           , intent(in)           :: xtype            ! netcdf data type
    character(len=*)  , intent(in)           :: long_name        ! long name for variable
    character(len=*)  , intent(in)           :: interpinic_flag  ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module)
    {VTYPE}           , pointer              :: data{DIMSTR}
    logical           , intent(inout)        :: readvar          ! was var read?
    character(len=*)  , intent(in), optional :: dim1name         ! dimension name
    character(len=*)  , intent(in), optional :: dim2name         ! dimension name
    character(len=*)  , intent(in), optional :: units            ! long name for variable
    character(len=*)  , intent(in), optional :: comment          ! attribute
    character(len=*)  , intent(in), optional :: flag_meanings(:) ! attribute
    real(r8)          , intent(in), optional :: missing_value    ! attribute for real
    real(r8)          , intent(in), optional :: fill_value       ! attribute for real
    integer           , intent(in), optional :: imissing_value   ! attribute for int
    integer           , intent(in), optional :: ifill_value      ! attribute for int
    integer           , intent(in), optional :: flag_values(:)   ! attribute for int
    integer           , intent(in), optional :: nvalid_range(2)  ! attribute for int
    logical           , intent(in), optional :: is_spatial       ! is this a spatial variable (with a dimension of gridcell, column, etc.)? If not present, assumed true.
    ! 
    ! Local variables
    character(len=len(varname)) :: primary_varname ! first name in the varname list
    character(len=len(varname)) :: my_varname      ! actual varname to read/write
    logical          :: l_is_spatial               ! local version of is_spatial
    integer          :: ivalue
    type(var_desc_t) :: vardesc  ! local vardesc
    integer          :: status   ! return error code 
    integer          :: varid
    integer          :: lxtype   ! local external type (in case logical variable)
    !----------------------------------------------------

    if (present(is_spatial)) then
       l_is_spatial = is_spatial
    else
       l_is_spatial = .true.
    end if

    call shr_string_listGetName(varname, 1, primary_varname)
    if (flag == 'read') then
       call find_var_on_file(ncid, varname, is_dim=.false., varname_on_file=my_varname)
       if ((my_varname /= primary_varname) .and. masterproc) then
          write(iulog,*) 'Restart file backwards compatibility: Translating: ', &
               trim(my_varname), ' => ', trim(primary_varname)
       end if
    else
       my_varname = primary_varname
    end if

    if (flag == 'define') then

       if ( xtype == ncd_log )then
          lxtype = ncd_int
       else
          lxtype = xtype
       end if

       if (.not. present(dim1name)) then
          call ncd_defvar(ncid=ncid, varname=trim(my_varname), xtype=lxtype, &
               long_name=trim(long_name), units=units)
       else if (.not. present(dim2name)) then 
          call ncd_defvar(ncid=ncid, varname=trim(my_varname), xtype=lxtype, &
               dim1name=trim(dim1name), &
               long_name=trim(long_name), units=units)
       else if (present(dim2name)) then
          call ncd_defvar(ncid=ncid, varname=trim(my_varname), xtype=lxtype, &
               dim1name=trim(dim1name), dim2name=trim(dim2name), &
               long_name=trim(long_name), units=units)
       end if

       status = PIO_inq_varid(ncid, trim(my_varname), vardesc)
       varid = vardesc%varid

       call write_interpinic_flag(ncid=ncid, varid=varid, varname=varname, &
            interpinic_flag=interpinic_flag)

       ! This attribute is written in order to communicate this metadata to initInterp
       call ncd_putatt(ncid, varid, 'varnames_on_old_files', trim(varname))

       if (present(comment)) then
          call ncd_putatt(ncid, varid, 'comment', trim(comment))
       end if
       if (present(units)) then
          call ncd_putatt(ncid, varid, 'units', trim(units))
       end if

       if (present(fill_value)) then
          call ncd_putatt(ncid, varid, '_FillValue', fill_value, lxtype)
       else if (lxtype == ncd_double) then
          call ncd_putatt(ncid, varid, '_FillValue', spval, lxtype)
       end if
       if (present(missing_value)) then
          call ncd_putatt(ncid, varid, 'missing_value', missing_value, lxtype)
       else if (lxtype == ncd_double) then
          call ncd_putatt(ncid, varid, 'missing_value', spval, lxtype)
       end if
       if (present(ifill_value)) then
          call ncd_putatt(ncid, varid, '_FillValue', ifill_value, lxtype)
       else if (lxtype == ncd_int) then
          call ncd_putatt(ncid, varid, '_FillValue', ispval, lxtype)
       end if
       if (present(imissing_value)) then
          call ncd_putatt(ncid, varid, 'missing_value', imissing_value, lxtype)
       else if (lxtype == ncd_int) then
          call ncd_putatt(ncid, varid, 'missing_value', ispval, lxtype)
       end if
       if (present(nvalid_range)) then
          status = PIO_put_att(ncid,varid,'valid_range', nvalid_range )
       end if
       if ( xtype == ncd_log )then
          status = PIO_put_att(ncid,varid,'flag_values',     (/0, 1/) )
          status = PIO_put_att(ncid,varid,'flag_meanings',  "FALSE TRUE" )
          status = PIO_put_att(ncid,varid,'valid_range',    (/0, 1/) )
       end if

    else if (flag == 'read' .or. flag == 'write') then

#if ({ITYPE}!=TYPETEXT) 
       if (.not. present(dim1name) .or. .not. l_is_spatial) then
          call ncd_io(varname=trim(my_varname), data=data, &
               ncid=ncid, flag=flag, readvar=readvar)
       else 
          call ncd_io(varname=trim(my_varname), data=data, &
            dim1name=trim(dim1name), ncid=ncid, flag=flag, readvar=readvar)
       end if
#endif
    end if

    if (flag == 'read') then
       if (.not. readvar) then
          call missing_field_possibly_abort(my_varname)
       end if
    end if

  end subroutine restartvar_{DIMS}d_{TYPE}

  !-----------------------------------------------------------------------
  !TYPE int,double
  subroutine restartvar_2d_{TYPE}_bounds(ncid, flag, varname, xtype, &
       dim1name, dim2name, switchdim, lowerb2, upperb2, &
       long_name, units, interpinic_flag, data, readvar, &
       comment, flag_meanings, missing_value, fill_value, &
       imissing_value, ifill_value, flag_values, nvalid_range )

    ! Note that varname can be a colon-delimited list of possible variable names (with no
    ! spaces around the colons). In this case, when flag = 'read', the input file is
    ! searched for each possible variable name in order, starting with the first, until
    ! one is found. This mechanism supports backwards compatibility with old restart
    ! files, in case variables have been renamed. For example, if variable 'foo' was
    ! recently renamed to 'bar', then varname should be 'bar:foo'. For flag = 'write',
    ! the first name in the list is used.

    !----------------------------------------------------
    ! Arguments
    type(file_desc_t), intent(inout)        :: ncid             ! netcdf file id
    character(len=*) , intent(in)           :: flag             ! 'read' or 'write'
    character(len=*) , intent(in)           :: varname          ! variable name (or colon-delimited list: see above)
    integer          , intent(in)           :: xtype            ! netcdf data type
    character(len=*) , intent(in)           :: dim1name         ! dimension name
    character(len=*) , intent(in)           :: dim2name         ! dimension name
    logical          , intent(in)           :: switchdim
    character(len=*) , intent(in)           :: long_name        ! long name for variable
    character(len=*) , intent(in)           :: interpinic_flag  ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module)
    {VTYPE}          , pointer              :: data(:,:)        ! raw data
    logical          , intent(out)          :: readvar          ! was var read?
    integer          , intent(in), optional :: lowerb2
    integer          , intent(in), optional :: upperb2 
    character(len=*) , intent(in), optional :: units            ! long name for variable
    character(len=*) , intent(in), optional :: comment          ! attribute
    character(len=*) , intent(in), optional :: flag_meanings(:) ! attribute
    real(r8)         , intent(in), optional :: missing_value    ! attribute for real
    real(r8)         , intent(in), optional :: fill_value       ! attribute for real
    integer          , intent(in), optional :: imissing_value   ! attribute for int
    integer          , intent(in), optional :: ifill_value      ! attribute for int
    integer          , intent(in), optional :: flag_values(:)   ! attribute for int
    integer          , intent(in), optional :: nvalid_range(2)  ! attribute for int
    ! 
    ! Local variables
    character(len=len(varname)) :: primary_varname ! first name in the varname list
    character(len=len(varname)) :: my_varname      ! actual varname to read/write
    integer          :: ivalue
    type(var_desc_t) :: vardesc  ! local vardesc
    integer          :: status   ! return error code 
    integer          :: varid    ! returned var id
    integer          :: lxtype   ! local external type (in case logical variable)
    !----------------------------------------------------

    call shr_string_listGetName(varname, 1, primary_varname)
    if (flag == 'read') then
       call find_var_on_file(ncid, varname, is_dim=.false., varname_on_file=my_varname)
       if ((my_varname /= primary_varname) .and. masterproc) then
          write(iulog,*) 'Restart file backwards compatibility: Translating: ', &
               trim(my_varname), ' => ', trim(primary_varname)
       end if
    else
       my_varname = primary_varname
    end if

    if (flag == 'define') then

       if ( xtype == ncd_log )then
          lxtype = ncd_int
       else
          lxtype = xtype
       end if

       if (switchdim) then
          call ncd_defvar(ncid=ncid, varname=trim(my_varname), xtype=lxtype, &
               dim1name=trim(dim2name), dim2name=trim(dim1name), &
               long_name=trim(long_name), units=units)
       else
          call ncd_defvar(ncid=ncid, varname=trim(my_varname), xtype=lxtype, &
               dim1name=trim(dim1name), dim2name=trim(dim2name), &
               long_name=trim(long_name), units=units)
       end if
       status = PIO_inq_varid(ncid, trim(my_varname), vardesc)

       varid = vardesc%varid

       call write_interpinic_flag(ncid=ncid, varid=varid, varname=varname, &
            interpinic_flag=interpinic_flag)

       ! This attribute is written in order to communicate this metadata to initInterp
       call ncd_putatt(ncid, varid, 'varnames_on_old_files', trim(varname))

       if (switchdim) then
          status = PIO_put_att(ncid, vardesc%varid, 'switchdim_flag', 1)
       else
          status = PIO_put_att(ncid, vardesc%varid, 'switchdim_flag', 0)
       end if
       status = PIO_put_att(ncid, vardesc%varid, 'switchdim_flag_values', (/0,1/))
       status = PIO_put_att(ncid, vardesc%varid, 'switchdim_flag_is_0', &
            "1st and 2nd dims are same as model representation")
       status = PIO_put_att(ncid, vardesc%varid, 'switchdim_flag_is_1', &
            "1st and 2nd dims are switched from model representation")

       if (present(comment)) then
          call ncd_putatt(ncid, varid, 'comment', trim(comment))
       end if
       if (present(units)) then
          call ncd_putatt(ncid, varid, 'units', trim(units))
       end if
       if (present(fill_value)) then
          call ncd_putatt(ncid, varid, '_FillValue', fill_value, lxtype)
       else if (lxtype == ncd_double) then
          call ncd_putatt(ncid, varid, '_FillValue', spval, lxtype)
       end if
       if (present(missing_value)) then
          call ncd_putatt(ncid, varid, 'missing_value', missing_value, lxtype)
       else if (lxtype == ncd_double) then
          call ncd_putatt(ncid, varid, 'missing_value', spval, lxtype)
       end if
       if (present(ifill_value)) then
          call ncd_putatt(ncid, varid, '_FillValue', ifill_value, lxtype)
       else if (lxtype == ncd_int) then
          call ncd_putatt(ncid, varid, '_FillValue', ispval, lxtype)
       end if
       if (present(imissing_value)) then
          call ncd_putatt(ncid, varid, 'missing_value', imissing_value, lxtype)
       else if (lxtype == ncd_int) then
          call ncd_putatt(ncid, varid, 'missing_value', ispval, lxtype)
       end if
       if (present(nvalid_range)) then
          status = PIO_put_att(ncid,varid,'valid_range', nvalid_range )
       end if
       if ( xtype == ncd_log )then
          status = PIO_put_att(ncid,varid,'flag_values',     (/0, 1/) )
          status = PIO_put_att(ncid,varid,'flag_meanings',  "FALSE TRUE" )
          status = PIO_put_att(ncid,varid,'valid_range',    (/0, 1/) )
       end if

    else 

       if (present(lowerb2) .and. present(upperb2)) then
          call ncd_io(varname=trim(my_varname), data=data, &
               dim1name=trim(dim1name), switchdim=switchdim, &
               lowerb2=lowerb2, upperb2=upperb2, &
               ncid=ncid, flag=flag, readvar=readvar)
       else
          call ncd_io(varname=trim(my_varname), data=data, &
               dim1name=trim(dim1name), switchdim=switchdim, &
               ncid=ncid, flag=flag, readvar=readvar)
       end if

    end if

    if (flag == 'read') then
       if (.not. readvar) then
          call missing_field_possibly_abort(my_varname)
       end if
    end if

  end subroutine restartvar_2d_{TYPE}_bounds

  !-----------------------------------------------------------------------
  subroutine set_missing_from_template(my_var, template_var, multiplier)
    !
    ! !DESCRIPTION:
    ! Set values of a missing restart field from a template field, with some constant
    ! multiplier
    !
    ! !USES:
    use shr_infnan_mod, only : isnan => shr_infnan_isnan, nan => shr_infnan_nan, assignment(=)
    !
    ! !ARGUMENTS:
    real(r8), intent(out) :: my_var(:)
    real(r8), intent(in)  :: template_var(:)
    real(r8), intent(in)  :: multiplier
    !
    ! !LOCAL VARIABLES:
    integer :: i

    character(len=*), parameter :: subname = 'set_missing_from_template'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(template_var) == ubound(my_var)), sourcefile, __LINE__)

    do i = 1, size(my_var)
       if (isnan(template_var(i))) then
          my_var(i) = nan
       else if (template_var(i) == spval) then
          my_var(i) = spval
       else
          my_var(i) = template_var(i) * multiplier
       end if
    end do

  end subroutine set_missing_from_template

  !-----------------------------------------------------------------------
  subroutine set_grc_field_from_col_field(bounds, ncid, varname, data_grc, readvar)
    !
    ! !DESCRIPTION:
    ! Set a gridcell-level field from a column-level field on the restart file.
    !
    ! If the column-level field isn't found, then 'data' remains unchanged, and readvar
    ! is set to .false.
    !
    ! This can be useful for backwards compatibility: If a field that was previously on
    ! the column level has been moved to the gridcell-level, then this routine can be
    ! used to read the old column-level field and then average it to the gridcell-level.
    !
    ! We may want to make this more general, so that a single routine can set a
    ! gridcell-level field from a patch-, column- or landunit-level field - and maybe
    ! also set a column-level field from a patch-level field, etc.
    !
    ! !USES:
    use subgridAveMod, only : c2g
    !
    ! !ARGUMENTS:
    type(bounds_type) , intent(in)    :: bounds
    type(file_desc_t) , intent(inout) :: ncid        ! netcdf file id
    character(len=*)  , intent(in)    :: varname     ! variable name
    real(r8)          , intent(inout) :: data_grc( bounds%begg: ) ! gridcell-level data to set
    logical           , intent(out)   :: readvar     ! was var read?
    !
    ! !LOCAL VARIABLES:
    real(r8), pointer :: data_col(:)

    character(len=*), parameter :: subname = 'set_grc_field_from_col_field'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(data_grc) == (/bounds%endg/)), sourcefile, __LINE__)

    allocate(data_col(bounds%begc:bounds%endc))
    call ncd_io(varname=trim(varname), data=data_col, &
         dim1name='column', &
         ncid=ncid, flag='read', readvar=readvar)

    if (readvar) then
       call c2g(bounds, data_col, data_grc, &
            c2l_scale_type = 'unity', &
            l2g_scale_type = 'unity')
    end if

    deallocate(data_col)
  end subroutine set_grc_field_from_col_field

  !-----------------------------------------------------------------------
  ! DIMS 1,2
  subroutine set_missing_vals_to_constant_{DIMS}d(data, val)
    !
    ! !DESCRIPTION:
    ! Set missing values (NaN or spval) of a restart field to some constant value
    !
    ! !USES:
    use shr_infnan_mod, only : isnan => shr_infnan_isnan
    !
    ! !ARGUMENTS:
    real(r8), intent(inout) :: data{DIMSTR}  ! variable to modify
    real(r8), intent(in)    :: val           ! value to replace any NaNs or spvals
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'set_missing_vals_to_constant_{DIMS}d'
    !-----------------------------------------------------------------------

    where (isnan(data))
       data = val
    elsewhere (data == spval)
       data = val
    end where

  end subroutine set_missing_vals_to_constant_{DIMS}d

  !-----------------------------------------------------------------------
  subroutine missing_field_possibly_abort(varname)
    !
    ! !DESCRIPTION:
    ! This should be called if a field is missing from the restart file. It aborts with a
    ! helpful error message if this run is one where a missing field is treated as an
    ! error (e.g., a restart run).
    !
    ! !ARGUMENTS:
    character(len=*), intent(in) :: varname ! name of variable that triggered the error
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'missing_field_possibly_abort'
    !-----------------------------------------------------------------------

    if (nsrest == nsrContinue .or. &
         nsrest == nsrBranch) then
       if (masterproc) then
          write(iulog,*) 'ERROR: Field missing from restart file: ', trim(varname)
          write(iulog,*) 'Missing fields are not allowed in branch or continue (restart) runs.'
          write(iulog,*) ' '
          write(iulog,*) 'This can happen when the restart file is from a different'
          write(iulog,*) 'model configuration or different code base, which did not'
          write(iulog,*) 'contain all of the restart fields needed for the current'
          write(iulog,*) 'code base and configuration.'
          write(iulog,*) ' '
          write(iulog,*) 'You can get around this problem by doing a startup or hybrid run'
          write(iulog,*) 'that points to this initial condition file, instead of attempting'
          write(iulog,*) 'a branch or continue run. However, note that, in many cases,'
          write(iulog,*) 'this will result in the missing fields being filled in with'
          write(iulog,*) 'their cold start initialization values, which may or may not'
          write(iulog,*) 'be what you want.'
       end if
       call shr_sys_abort('Field missing from restart file.')
    end if
  end subroutine missing_field_possibly_abort

  !-----------------------------------------------------------------------
  subroutine write_interpinic_flag(ncid, varid, varname, interpinic_flag)
    !
    ! !DESCRIPTION:
    ! Write interpinic_flag metadata for the given variable
    !
    ! !ARGUMENTS:
    type(file_desc_t) , intent(inout) :: ncid            ! netcdf file id
    integer           , intent(in)    :: varid           ! variable id
    character(len=*)  , intent(in)    :: varname         ! variable name (just used for output on error)
    character(len=*)  , intent(in)    :: interpinic_flag ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module)
    !
    ! !LOCAL VARIABLES:
    integer :: status   ! return error code

    character(len=*), parameter :: subname = 'write_interpinic_flag'
    !-----------------------------------------------------------------------

    if (interpinic_flag == 'interp') then
       status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_interp)
    else if (interpinic_flag == 'copy') then
       status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_copy)
    else if (interpinic_flag == 'skip') then
       status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_skip)
    else if (interpinic_flag == 'area') then
       status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_area)
    else
       write(iulog,*) 'ERROR in restartvar for ', trim(varname)
       write(iulog,*) 'Unknown interpinic_flag: ', trim(interpinic_flag)
       write(iulog,*) 'Allowed values are: interp, copy, skip, area'
       call endrun(msg='Unknown interpinic_flag')
    end if
    status = PIO_put_att(ncid, varid, 'interpinic_flag_meanings', &
         "1=nearest neighbor, 2=copy directly, 3=skip, 4=area")

  end subroutine write_interpinic_flag


end module restUtilMod