module subcol_utils
   !---------------------------------------------------------------------------
   ! Purpose:
   !
   ! Provides utilities to support subcolumns
   !
   !---------------------------------------------------------------------------

   use shr_kind_mod,    only: r8=>shr_kind_r8, r4=>shr_kind_r4, i4=>shr_kind_i4
   use infnan,          only: nan, assignment(=)
   use physics_types,   only: physics_state, physics_ptend, physics_tend, physics_tend_alloc, physics_state_alloc
   use ppgrid,          only: pcols, psubcols, pver
   use constituents,    only: pcnst
   use cam_abortutils,  only: endrun
   use pio,             only: var_desc_t

   implicit none

   private
   save

   !! Private variable to provide default packing and unpacking of fields
   !! for use in restart functionality. Allocated as (pcols, begchunk:endchunk)
   integer, target, allocatable :: nsubcol2d(:,:)
   integer, target, allocatable :: indcol2d(:,:)
   integer, target, allocatable :: filter2d(:,:)
   real(r8),target, allocatable :: weight2d(:,:)
   logical :: weight_set, filter_set

   !! The active subcolumn scheme
   character(len=16) :: subcol_scheme = 'off'

   !! Public interface functions which do not depend on the subcolumn scheme

   public :: subcol_utils_init       ! Initialize module data (e.g., nsubcol2d)
   public :: subcol_get_nsubcol      ! Copy chunk from nsubcol2d
   public :: subcol_set_nsubcol      ! Copy chunk to nsubcol2d
   public :: subcol_get_indcol       ! Copy chunk from indcol2d
   public :: subcol_get_filter       ! return the filter values
   public :: subcol_set_filter       ! set the filter values
   public :: subcol_get_weight       ! return the weight values
   public :: subcol_set_weight       ! set the weight values

   public :: subcol_field_copy       ! copy a physics buffer field into one with subcolumn dimensions
   public :: subcol_ptend_copy       ! copy a physics_ptend object into one with subcolumn dimensions
   public :: subcol_set_subcols      ! set nsubcols and copy state & tend objects into one with subcolumn dimensions

   public :: subcol_field_avg_shr    ! Average subcol fields back into GBA fields
   public :: subcol_ptend_avg_shr    ! average subcolumn ptend to grid ptend

   public :: subcol_field_get_firstsubcol  ! Retrieve the first subcolumn and assign to grid
   public :: subcol_ptend_get_firstsubcol  ! retrieve the first subcolumn from the ptend fields and assign to grid ptend

   public :: subcol_unpack           ! Unpack a subcolumn field
   public :: subcol_pack             ! Pack a subcolumn field

   public :: subcol_utils_init_restart     ! Initialize restart with subcolumn specific fields
   public :: subcol_utils_read_restart     ! Read subcolumn specific fields from restart
   public :: subcol_utils_write_restart    ! Write subcolumn specific fields for restart
   public :: is_filter_set           ! True if filters for averaging have been set
   public :: is_weight_set           ! True if weights for averaging have been set
   public :: is_subcol_on            ! true is any subcol_scheme other than "off" is set
   public :: subcol_get_scheme       ! Return the active subcolumn scheme name
   public :: subcol_utils_readnl     ! Set the active scheme based on namelist

   interface subcol_field_avg_shr
      ! TYPE int,double,real
      ! DIMS 1,2
      module procedure subcol_field_avg_shr_{DIMS}d{TYPE}
   end interface

   interface subcol_avg_inter
      ! TYPE int,double,real
      module procedure subcol_avg_inter_{TYPE}
   end interface

   interface subcol_avg
      ! TYPE int,double,real
      module procedure subcol_avg_{TYPE}
   end interface

   interface subcol_field_get_firstsubcol
      ! TYPE int,double,real
      ! DIMS 1,2
      module procedure subcol_field_get_firstsubcol_{DIMS}d{TYPE}
   end interface

   interface subcol_state_field_copy
      ! TYPE double
      ! DIMS 1,2,3
      module procedure subcol_state_field_copy_{DIMS}d{TYPE}
   end interface

   interface subcol_field_copy
      ! TYPE int,double,real
      ! DIMS 1,2,3,4,5
      module procedure subcol_field_copy_{DIMS}d{TYPE}
   end interface

   interface subcol_pack
      ! TYPE int,double,real
      ! DIMS 1,2,3,4,5,6
      module procedure subcol_pack_{DIMS}d_{TYPE}
   end interface

   interface subcol_unpack
      ! TYPE int,double,real
      ! DIMS 1,2,3,4,5,6
      module procedure subcol_unpack_{DIMS}d_{TYPE}
   end interface

   type(var_desc_t) :: nsubcol_desc, weight2d_desc, filter2d_desc
   integer          :: subcol_dimid = -1 ! subcol dimension for restart

   integer :: ret_nan_int
   real(r8):: ret_nan_double
   real(r4):: ret_nan_real

   integer :: fillval_int
   real(r8):: fillval_double
   real(r4):: fillval_real

contains

   subroutine subcol_allocate_internal()
      use ppgrid,         only: begchunk, endchunk
   !-----------------------------------------------------------------------
   ! Allocate nsubcol2d and indcol2d
   !-----------------------------------------------------------------------
      if (allocated(nsubcol2d)) then
         deallocate(nsubcol2d)
      end if
      allocate(nsubcol2d(pcols, begchunk:endchunk))
      nsubcol2d = 0

      if (allocated(indcol2d)) then
         deallocate(indcol2d)
      end if
      allocate(indcol2d(pcols*psubcols, begchunk:endchunk))
      indcol2d = 0

      if (allocated(filter2d)) then
         deallocate(filter2d)
      end if
      allocate(filter2d(pcols*psubcols, begchunk:endchunk))
      filter2d = 0

      if (allocated(weight2d)) then
         deallocate(weight2d)
      end if
      allocate(weight2d(pcols*psubcols, begchunk:endchunk))
      weight2d = 0._r8

    end subroutine subcol_allocate_internal

   subroutine subcol_utils_init(subcol_scheme_init)
   !-----------------------------------------------------------------------
   ! Initialize the nsubcol module variable
   !-----------------------------------------------------------------------
 
      character(len=*), optional, intent(in) :: subcol_scheme_init ! Name of subcolumn generator

      call subcol_allocate_internal()

      ret_nan_int    = 0
      ret_nan_double = nan
      ret_nan_real   = nan

      fillval_int    = 0
      fillval_double = 0._r8
      fillval_real   = 0._r4
 
      weight_set = .false.
      filter_set = .false.

   end subroutine subcol_utils_init

   subroutine subcol_get_nsubcol(lchnk, nsubcol)
   !-----------------------------------------------------------------------
   ! Retrieve a chunk from the nsubcol module variable
   !-----------------------------------------------------------------------

     integer,    intent(in)    :: lchnk
     integer,    intent(out)   :: nsubcol(:)

     if (.not. allocated(nsubcol2d)) then
       call endrun('subcol_get_nsubcol: nsubcol2d not allocated')
     end if
     nsubcol(:) = nsubcol2d(:,lchnk)
   end subroutine subcol_get_nsubcol

   subroutine subcol_get_indcol(lchnk, indcol)
   !-----------------------------------------------------------------------
   ! Retrieve a chunk from the nsubcol module variable
   !-----------------------------------------------------------------------

     integer,    intent(in)    :: lchnk
     integer,    intent(out)   :: indcol(:)

     if (.not. allocated(indcol2d)) then
       call endrun('subcol_get_indcol: indcol2d not allocated')
     end if
     indcol(:) = indcol2d(:,lchnk)
  end subroutine subcol_get_indcol

   subroutine subcol_get_filter(lchnk, filter)
   !-----------------------------------------------------------------------
   ! Retrieve the filter module variable
   !-----------------------------------------------------------------------

     integer,    intent(in)    :: lchnk
     integer,    intent(out)   :: filter(:)

     filter(:) = filter2d(:,lchnk)
   end subroutine subcol_get_filter

   subroutine subcol_get_weight(lchnk, weight)
   !-----------------------------------------------------------------------
   ! Retrieve the weight module variable
   !-----------------------------------------------------------------------

     integer,    intent(in)    :: lchnk
     real(r8),   intent(out)   :: weight(:)

     weight(:) = weight2d(:,lchnk)
   end subroutine subcol_get_weight

   integer function subcol_get_ncol(lchnk) result(ncol)
   !-----------------------------------------------------------------------
   ! Compute the number of (sub)columns for a chunk
   ! NB: This is considered an internal function so it can use nsubcol2d
   !-----------------------------------------------------------------------
     integer,    intent(in)    :: lchnk

     ncol = sum(nsubcol2d(:,lchnk))

   end function subcol_get_ncol

   subroutine subcol_set_nsubcol(lchnk, ngrdcol, nsubcol)
   !-----------------------------------------------------------------------
   ! Set a chunk of the nsubcol module variable
   ! Also, recompute indcol for lchnk
   !-----------------------------------------------------------------------

      integer,    intent(in)    :: lchnk
      integer,    intent(in)    :: ngrdcol
      integer,    intent(in)    :: nsubcol(:)

      integer                   :: i, j, indx

      if (any(nsubcol(:) > psubcols)) then
         call endrun('subcol_set_nsubcol: psubcols not set large enough to hold the number of subcolumns requested')
      end if
      if (any(nsubcol(:) < 0)) then
         call endrun('subcol_set_nsubcol: nsubcols must be non-negative')
      end if
      if (ngrdcol < pcols) then
         if (any(nsubcol(ngrdcol+1:) > 0)) then
            call endrun('subcol_set_nsubcol: Cannot set subcolumns for columns past ngrdcol')
         end if
      end if
      nsubcol2d(:, lchnk) = nsubcol(:)
      ! Recalculate indcol for the chunk
      indx = 1
      do i = 1, pcols
         do j = 1, nsubcol2d(i, lchnk)
            indcol2d(indx, lchnk) = i
            indx = indx + 1
         end do
      end do
      ! Fill with zeros
      if (indx <= pcols * psubcols) then
         indcol2d(indx:pcols*psubcols, lchnk) = 0
      end if
   end subroutine subcol_set_nsubcol

   subroutine subcol_set_filter(lchnk, filter)
   !-----------------------------------------------------------------------
   ! Set the filter module variable
   !-----------------------------------------------------------------------

     integer,    intent(in)    :: lchnk
     integer,    intent(in)    :: filter(:)

     filter2d(:,lchnk) = filter(:)
     filter_set        = .true.
   end subroutine subcol_set_filter

   subroutine subcol_set_weight(lchnk, weight)
   !-----------------------------------------------------------------------
   ! Set the weight module variable
   !-----------------------------------------------------------------------

     integer,    intent(in)    :: lchnk
     real(r8),   intent(in)    :: weight(:)

     weight2d(:,lchnk) = weight(:)
     weight_set        = .true.
   end subroutine subcol_set_weight

   logical function is_weight_set()
      is_weight_set=weight_set 
   end function is_weight_set

   logical function is_filter_set()
      is_filter_set=filter_set 
   end function is_filter_set

   logical function is_subcol_on()
      is_subcol_on = (trim(subcol_scheme) /= 'off')
   end function is_subcol_on

   character(len=16) function subcol_get_scheme()
      subcol_get_scheme = trim(subcol_scheme)
   end function subcol_get_scheme

   subroutine subcol_utils_readnl(nlfile)
      use namelist_utils,  only: find_group_name
      use units,           only: getunit, freeunit
      use spmd_utils,      only: masterproc, mpi_character, masterprocid, mpicom
      
      character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input

      ! Local variables
      integer :: unitn, ierr

      namelist /subcol_nl/ subcol_scheme

      if (masterproc) then
         unitn = getunit()
         open( unitn, file=trim(nlfile), status='old' )
         call find_group_name(unitn, 'subcol_nl', status=ierr)
         if (ierr == 0) then
            read(unitn, subcol_nl, iostat=ierr)
            if (ierr /= 0) then
               call endrun('subcol_readnl: ERROR reading namelist')
            end if
         end if
         close(unitn)
         call freeunit(unitn)
      end if

#ifdef SPMD
      ! Broadcast namelist variables
      call mpi_bcast(subcol_scheme, len(subcol_scheme), mpi_character, masterprocid, mpicom, ierr)
#endif

   end subroutine subcol_utils_readnl

   ! TYPE int,double,real
   ! DIMS 1,2,3,4,5
   subroutine subcol_field_copy_{DIMS}d{TYPE} (field, lchnk, field_sc)
   !-----------------------------------------------------------------------
   ! Copy a pbuf field dimensioned pcols to one with pcols*psubcols and fill appropriately
   !-----------------------------------------------------------------------

      {VTYPE},          intent(in)      :: field{DIMSTR}
      integer,          intent(in)      :: lchnk
      {VTYPE},          intent(out)     :: field_sc{DIMSTR}   ! intent out

      ! Local Variables
      integer                           :: ncol
      integer                           :: indcol(pcols*psubcols)

      if (size(field,dim=1) .ne. pcols) then
         call endrun('subcol_field_copy error: only fields with first dimension pcols may use this routine')
      end if

      field_sc = fillval_{TYPE}
      ncol = subcol_get_ncol(lchnk)
      call subcol_get_indcol(lchnk, indcol)

#if ({DIMS} == 1)
      field_sc(:ncol) = field(indcol(:ncol))
#endif

#if ({DIMS} == 2)
      field_sc(:ncol,:) = field(indcol(:ncol),:)
#endif

#if ({DIMS} == 3)
      field_sc(:ncol,:,:) = field(indcol(:ncol),:,:)
#endif

#if ({DIMS} == 4)
      field_sc(:ncol,:,:,:) = field(indcol(:ncol),:,:,:)
#endif

#if ({DIMS} == 5)
      field_sc(:ncol,:,:,:,:) = field(indcol(:ncol),:,:,:,:)
#endif

   end subroutine subcol_field_copy_{DIMS}d{TYPE}

   ! TYPE double
   ! DIMS 1,2,3
   subroutine subcol_state_field_copy_{DIMS}d{TYPE} (field, state, field_sc)
   !-----------------------------------------------------------------------
   ! Copy a state field dimensioned with pcols to one with pcols*psubcols and fill appropriately
   !-----------------------------------------------------------------------

      {VTYPE}, intent(in)                     :: field{DIMSTR}
      type(physics_state), intent(in)         :: state
      {VTYPE},          allocatable           :: field_sc{DIMSTR}   ! intent out

      integer :: dim2, dim3
      integer :: indcol(pcols*psubcols)

      if (size(field,dim=1) .ne. pcols) then
         call endrun('subcol_state_field_copy error: only fields with first dimension pcols may use this routine')
      end if

      call subcol_get_indcol(state%lchnk, indcol)

#if ({DIMS} == 1)
      if (.not. allocated(field_sc)) then
         allocate(field_sc(pcols*psubcols))
      end if

      field_sc = 0._r8
      field_sc(:state%ncol) = field(indcol(:state%ncol))
#endif

#if ({DIMS} == 2)
      if (.not. allocated(field_sc)) then
         dim2 = size(field,dim=2)
         allocate(field_sc(pcols*psubcols,dim2))
      end if

      field_sc = 0._r8
      field_sc(:state%ncol,:) = field(indcol(:state%ncol),:)
#endif

#if ({DIMS} == 3)
      if (.not. allocated(field_sc)) then
         dim2 = size(field,dim=2)
         dim3 = size(field,dim=3)
         allocate(field_sc(pcols*psubcols,dim2,dim3))
      end if

      field_sc = 0._r8
      field_sc(:state%ncol,:,:) = field(indcol(:state%ncol),:,:)
#endif

   end subroutine subcol_state_field_copy_{DIMS}d{TYPE}

   subroutine subcol_tend_copy(tend, state_sc, tend_sc)
   !-----------------------------------------------------------------------
   ! Copy all of tend to subcolumns in tend_sc, allocating tend_sc if necessary
   !-----------------------------------------------------------------------
      type(physics_tend),  intent(inout) :: tend
      type(physics_state), intent(in)    :: state_sc        ! subcolumn state
      type(physics_tend),  intent(inout) :: tend_sc         ! subcolumn tend


      if (.not. allocated(tend%dtdt)) then
         call physics_tend_alloc(tend_sc, state_sc%psetcols)
      end if

      tend_sc%psetcols   = pcols*psubcols
      call subcol_state_field_copy(tend%dtdt,       state_sc, tend_sc%dtdt)
      call subcol_state_field_copy(tend%dudt,       state_sc, tend_sc%dudt)
      call subcol_state_field_copy(tend%dvdt,       state_sc, tend_sc%dvdt)
      call subcol_state_field_copy(tend%flx_net,    state_sc, tend_sc%flx_net)
      call subcol_state_field_copy(tend%te_tnd,     state_sc, tend_sc%te_tnd)
      call subcol_state_field_copy(tend%tw_tnd,     state_sc, tend_sc%tw_tnd)

   end subroutine subcol_tend_copy

   subroutine subcol_state_copy(state, state_sc)
   !-----------------------------------------------------------------------
   ! Copy all of state to subcolumns in state_sc, allocating state_sc if necessary
   !-----------------------------------------------------------------------
      type(physics_state), intent(in)    :: state
      type(physics_state), intent(inout) :: state_sc        ! subcolumn state

      !
      ! Local variables
      !
      integer :: ngrdcol

      if (.not. allocated(state_sc%lat)) then
         call endrun('subcol_state_copy: user must allocate state_sc prior to calling this routine')
      end if

      ngrdcol       = state%ngrdcol

      call subcol_state_hdrinit(state, state_sc)

      call subcol_state_field_copy(state%lat,       state_sc, state_sc%lat)
      call subcol_state_field_copy(state%lon,       state_sc, state_sc%lon)
      call subcol_state_field_copy(state%ps,        state_sc, state_sc%ps)
      call subcol_state_field_copy(state%psdry,     state_sc, state_sc%psdry)
      call subcol_state_field_copy(state%phis,      state_sc, state_sc%phis)
      call subcol_state_field_copy(state%te_ini,    state_sc, state_sc%te_ini)
      call subcol_state_field_copy(state%te_cur,    state_sc, state_sc%te_cur)
      call subcol_state_field_copy(state%tw_ini,    state_sc, state_sc%tw_ini)
      call subcol_state_field_copy(state%tw_cur,    state_sc, state_sc%tw_cur)
      call subcol_state_field_copy(state%t,         state_sc, state_sc%t)
      call subcol_state_field_copy(state%u,         state_sc, state_sc%u)
      call subcol_state_field_copy(state%v,         state_sc, state_sc%v)
      call subcol_state_field_copy(state%s,         state_sc, state_sc%s)
      call subcol_state_field_copy(state%omega,     state_sc, state_sc%omega)
      call subcol_state_field_copy(state%pmid,      state_sc, state_sc%pmid)
      call subcol_state_field_copy(state%pdel,      state_sc, state_sc%pdel)
      call subcol_state_field_copy(state%rpdel,     state_sc, state_sc%rpdel)
      call subcol_state_field_copy(state%lnpmid,    state_sc, state_sc%lnpmid)
      call subcol_state_field_copy(state%exner,     state_sc, state_sc%exner)
      call subcol_state_field_copy(state%zm,        state_sc, state_sc%zm)
      call subcol_state_field_copy(state%pint,      state_sc, state_sc%pint)
      call subcol_state_field_copy(state%lnpint,    state_sc, state_sc%lnpint)
      call subcol_state_field_copy(state%zi,        state_sc, state_sc%zi)
      call subcol_state_field_copy(state%lnpmiddry, state_sc, state_sc%lnpmiddry)
      call subcol_state_field_copy(state%pmiddry,   state_sc, state_sc%pmiddry)
      call subcol_state_field_copy(state%pdeldry,   state_sc, state_sc%pdeldry)
      call subcol_state_field_copy(state%rpdeldry,  state_sc, state_sc%rpdeldry)
      call subcol_state_field_copy(state%pintdry,   state_sc, state_sc%pintdry)
      call subcol_state_field_copy(state%lnpintdry, state_sc, state_sc%lnpintdry)
      call subcol_state_field_copy(state%q,         state_sc, state_sc%q)

   end subroutine subcol_state_copy

   subroutine subcol_set_subcols(state, tend, nsubcol, state_sc, tend_sc)
   !-----------------------------------------------------------------------
   ! Propogate nsubcol information to common areas such as state, tend,
   ! nsubcol2d, and indcol2d
   !-----------------------------------------------------------------------
      type(physics_state), intent(in)    :: state
      type(physics_tend),  intent(inout) :: tend
      integer,             intent(in)    :: nsubcol(pcols)
      type(physics_state), intent(inout) :: state_sc        ! subcolumn state
      type(physics_tend),  intent(inout) :: tend_sc         ! subcolumn tend

      call subcol_set_nsubcol(state%lchnk, state%ngrdcol, nsubcol)
      call subcol_state_copy(state, state_sc)
      call subcol_tend_copy(tend, state_sc, tend_sc)
   end subroutine subcol_set_subcols


   subroutine subcol_ptend_copy(ptend, state, ptend_cp)
   !-----------------------------------------------------------------------
   ! Copy a physics_ptend object into one which has subcolumns
   !-----------------------------------------------------------------------
      use physics_types, only: physics_ptend_init

      type(physics_ptend), intent(in)    :: ptend    ! ptend source, dimensioned with grid columns
      type(physics_state), intent(in)    :: state    ! state with subcolumns
      type(physics_ptend), intent(out)   :: ptend_cp ! copy of ptend, dimensioned with subcolumns

      ! Local Variables
      integer                            :: ncol
      integer                            :: indcol(pcols*psubcols)

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

      if (ptend%psetcols .ne. pcols) then
         call endrun('subcol_ptend_copy: ptend must be dimensioned pcols to use this routine')
      end if

      call physics_ptend_init(ptend_cp,state%psetcols, ptend%name, ls=ptend%ls, lu=ptend%lu,  &
           lv=ptend%lv, lq=ptend%lq)

      ptend_cp%top_level = ptend%top_level
      ptend_cp%bot_level = ptend%bot_level

      ncol = subcol_get_ncol(state%lchnk)
      call subcol_get_indcol(state%lchnk, indcol)

      ! Copy the grid column data into each of the subcolumns - indcol contains the grid index for each subcolumn
      if (ptend_cp%ls) then
         ptend_cp%s(:ncol,:) = ptend%s(indcol(:ncol),:)
         ptend_cp%hflux_srf(:ncol) = ptend%hflux_srf(indcol(:ncol))
         ptend_cp%hflux_top(:ncol) = ptend%hflux_top(indcol(:ncol))
      end if

      if (ptend_cp%lu) then
         ptend_cp%u(:ncol,:) = ptend%u(indcol(:ncol),:)
         ptend_cp%taux_srf(:ncol)  = ptend%taux_srf(indcol(:ncol))
         ptend_cp%taux_top(:ncol)  = ptend%taux_top(indcol(:ncol))
      end if

      if (ptend_cp%lv) then
         ptend_cp%v(:ncol,:) = ptend%v(indcol(:ncol),:)
         ptend_cp%tauy_srf(:ncol)  = ptend%tauy_srf(indcol(:ncol))
         ptend_cp%tauy_top(:ncol)  = ptend%tauy_top(indcol(:ncol))
      end if

      if (any(ptend_cp%lq(:))) then
         ptend_cp%q(:ncol,:,:) = ptend%q(indcol(:ncol),:,:)
         ptend_cp%cflx_srf(:ncol,:)  = ptend%cflx_srf(indcol(:ncol),:)
         ptend_cp%cflx_top(:ncol,:)  = ptend%cflx_top(indcol(:ncol),:)
      end if

   end subroutine subcol_ptend_copy
  
   subroutine subcol_state_hdrinit(state, state_sc)
   !-----------------------------------------------------------------------
   ! Initialize the subcolumn state header variables
   !----------------------------------------------------------------------
      type(physics_state), intent(in)    :: state
      type(physics_state), intent(inout) :: state_sc        ! subcolumn state

      integer :: isize
      integer :: nsubcol(pcols)

      call subcol_get_nsubcol(state%lchnk, nsubcol)
      isize = state%ngrdcol
      if (size(nsubcol) < isize) then
         call endrun('subcol_state_hdrinit error: input nsubcol array must be dimensioned at least as large as state%ngrdcol')
      end if

      state_sc%ngrdcol    = state%ngrdcol
      state_sc%lchnk      = state%lchnk

      state_sc%psetcols   = pcols*psubcols
      ! Set count to a too-large value. It should be correctly initialized in check_energy_timestep_init
      state_sc%count      = pcols*psubcols*2

      ! Set the number of set subcolumns to the total count
      state_sc%ncol       = subcol_get_ncol(state%lchnk)

   end subroutine subcol_state_hdrinit

   ! TYPE int,double,real
   ! DIMS 1,2
   subroutine subcol_field_avg_shr_{DIMS}d{TYPE}(field_sc, ngrdcol, lchnk, field,  usefilter, useweight)
   !-----------------------------------------------------------------------
   ! This is the high level subcol field averaging routine which 
   ! averages a field dimensioned pcols*psubcols to a grid one dimensioned pcols 
   !  
   ! Its purpose is to check filter and weight settings and to subset the
   ! appropriate subsection of the field array to pass on the averaging routine
   !-----------------------------------------------------------------------

      {VTYPE},intent(in)                      :: field_sc{DIMSTR} ! intent in
      integer, intent(in)                     :: ngrdcol          ! # grid cols
      integer, intent(in)                     :: lchnk            ! chunk index
      {VTYPE}, intent(out)                    :: field{DIMSTR}
      logical, intent(in)                     :: usefilter
      logical, intent(in)                     :: useweight

      !
      ! Local variables
      !
      integer           :: indx, nsubcol, i, j
      integer           :: nsubcols(pcols)

      field = ret_nan_{TYPE}
     
      if (usefilter .and. .not. filter_set) then
        call endrun('subcol_field_avg_shr error: Cannot specify using filters when none set')
      end if

      if (useweight .and. .not. weight_set) then
        call endrun('subcol_field_avg_shr error: Cannot specify using weights when none set')
      end if

      call subcol_get_nsubcol(lchnk, nsubcols)
      indx = 1
      do i = 1, ngrdcol
         if (nsubcols(i) .gt. 0) then
           nsubcol = nsubcols(i)
#if ({DIMS} >=2)
           do j=1,size(field_sc,dim=2)
#endif
#if ({DIMS} == 1)
             field(i)=subcol_avg_inter(field_sc(indx:indx+nsubcol-1), lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
#endif
#if ({DIMS} == 2)
             field(i,j)=subcol_avg_inter(field_sc(indx:indx+nsubcol-1,j), lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
#endif
#if ({DIMS} >=2)
           end do
#endif
           indx = indx + nsubcol
         end if
      end do
   end subroutine subcol_field_avg_shr_{DIMS}d{TYPE}


   ! TYPE int,double,real
   ! DIMS 1,2
   subroutine subcol_field_get_firstsubcol_{DIMS}d{TYPE}(field_sc, docheck, ngrdcol, lchnk, field)
   !-----------------------------------------------------------------------
   ! Retrieve the first subcolumn from a field dimensioned pcols*psubcols 
   ! and assign to one with pcols, performing optional checking that all other subcolumns are identical
   !-----------------------------------------------------------------------

      {VTYPE}, intent(in)                     :: field_sc{DIMSTR} ! intent in
      logical, intent(in)                     :: docheck          ! true=check, false=no check
      integer, intent(in)                     :: ngrdcol          ! # grid cols
      integer, intent(in)                     :: lchnk            ! chunk index
      {VTYPE}, intent(out)                    :: field{DIMSTR}

      !
      ! Local variables
      !
      integer           :: indx, nsubcol, i, j, l
      integer           :: nsubcols(pcols)

      call subcol_get_nsubcol(lchnk, nsubcols)
      field = 0
      indx = 1
      do i = 1, ngrdcol
         if (nsubcols(i) .gt. 0) then
            nsubcol = nsubcols(i)
#if ({DIMS}>=2)
            do j=1,size(field_sc,dim=2)
#endif

#if ({DIMS} == 1)
               field(i) = field_sc(indx)
#elif ({DIMS} == 2)
               field(i,j) = field_sc(indx,j)
#endif
               if (docheck) then
                  do l=1,nsubcol-1
#if ({DIMS} == 1)
                    if (field_sc(indx) /= field_sc(indx+l)) then
#elif ({DIMS} == 2)
                    if (field_sc(indx,j) /= field_sc(indx+l,j)) then
#endif
                      call endrun('subcol_field_get_firstsubcol error: Not all subcolumn fields are identical')
                    end if
                  end do
               end if
#if ({DIMS}>=2)
            end do
#endif
            indx = indx + nsubcol
         end if
      end do
   end subroutine subcol_field_get_firstsubcol_{DIMS}d{TYPE}

   subroutine subcol_ptend_avg_shr(ptend_sc, ngrdcol, lchnk, ptend, usefilter, useweight)
   !-----------------------------------------------------------------------
   ! Average a subcolumn ptend to a grid ptend
   !-----------------------------------------------------------------------

      type(physics_ptend),  intent(in)    :: ptend_sc      ! subcolumn ptend 
      integer,              intent(in)    :: ngrdcol       ! # grid cols
      integer,              intent(in)    :: lchnk         ! chunk index
      type(physics_ptend),  intent(inout) :: ptend         ! grid ptend
      logical,              intent(in)    :: usefilter
      logical,              intent(in)    :: useweight

      !
      ! Local variables
      !
      integer           :: indx, i, j, k, nsubcol
      integer           :: nsubcols(pcols)

      if (usefilter .and. .not. filter_set) then
         call endrun('subcol_ptend_avg_shr error: Cannot specify using filters when none set')
      end if

      if (useweight .and. .not. weight_set) then
         call endrun('subcol_ptend_avg_shr error: Cannot specify using weights when none set')
      end if

      call subcol_get_nsubcol(lchnk, nsubcols)
      ! physics_ptend_init has already been called by the master interface
      if (ptend%ls) then
         ptend%s(:,:)        = 0._r8
         ptend%hflux_srf(:)  = 0._r8
         ptend%hflux_top(:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            if (nsubcols(i) > 0) then
               nsubcol = nsubcols(i)
               do j=1,pver
                  ptend%s(i,j)=subcol_avg_inter(ptend_sc%s(indx:indx+nsubcol-1,j), &
                       lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               end do
               ptend%hflux_srf(i) = subcol_avg_inter(ptend_sc%hflux_srf(indx:indx+nsubcol-1), &
                    lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               ptend%hflux_top(i) = subcol_avg_inter(ptend_sc%hflux_top(indx:indx+nsubcol-1), &
                    lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               indx = indx + nsubcol
            end if
         end do
      end if
      if (ptend%lu) then
         ptend%u(:,:)       = 0._r8
         ptend%taux_srf(:)  = 0._r8
         ptend%taux_top(:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            if (nsubcols(i) > 0) then
               nsubcol = nsubcols(i)
               do j=1,pver
                  ptend%u(i,j)=subcol_avg_inter(ptend_sc%u(indx:indx+nsubcol-1,j), &
                       lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               end do
               ptend%taux_srf(i) = subcol_avg_inter(ptend_sc%taux_srf(indx:indx+nsubcol-1), &
                    lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               ptend%taux_top(i) = subcol_avg_inter(ptend_sc%taux_top(indx:indx+nsubcol-1), &
                    lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               indx = indx + nsubcol
            end if
         end do
      end if
      if (ptend%lv) then
         ptend%v(:,:)       = 0._r8
         ptend%tauy_srf(:)  = 0._r8
         ptend%tauy_top(:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            if (nsubcols(i) > 0) then
               nsubcol = nsubcols(i)
               do j=1,pver
                  ptend%v(i,j)=subcol_avg_inter(ptend_sc%v(indx:indx+nsubcol-1,j), &
                       lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               end do
               ptend%tauy_srf(i) = subcol_avg_inter(ptend_sc%tauy_srf(indx:indx+nsubcol-1), &
                    lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               ptend%tauy_top(i) = subcol_avg_inter(ptend_sc%tauy_top(indx:indx+nsubcol-1), &
                    lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               indx = indx + nsubcol
            end if
         end do
      end if
      if (any(ptend%lq(:))) then
         ptend%q(:,:,:)       = 0._r8
         ptend%cflx_srf(:,:)  = 0._r8
         ptend%cflx_top(:,:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            if (nsubcols(i) > 0) then
               nsubcol = nsubcols(i)
               do j=1,pver
                  do k=1, pcnst
                     ptend%q(i,j,k)=subcol_avg_inter(ptend_sc%q(indx:indx+nsubcol-1,j,k), &
                          lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
                  end do
               end do
               do k=1,pcnst
                  ptend%cflx_srf(i,k) = subcol_avg_inter(ptend_sc%cflx_srf(indx:indx+nsubcol-1,k), &
                       lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
                  ptend%cflx_top(i,k) = subcol_avg_inter(ptend_sc%cflx_top(indx:indx+nsubcol-1,k), &
                       lchnk, i, indx, indx+nsubcol-1, usefilter, useweight)
               end do
               indx = indx + nsubcol
            end if
         end do
      end if

   end subroutine subcol_ptend_avg_shr

   subroutine subcol_ptend_get_firstsubcol(ptend_sc, docheck, ngrdcol, lchnk, ptend)
   !-----------------------------------------------------------------------
   ! Retrieve the first subcolumn from a ptend field dimensioned pcols*psubcols 
   ! and assign to one with pcols, performing optional check that all other subcolumns are identical
   !-----------------------------------------------------------------------

      type(physics_ptend),  intent(in)    :: ptend_sc      ! subcolumn ptend 
      logical,              intent(in)    :: docheck       ! perform check that all subcolumns match
      integer,              intent(in)    :: ngrdcol       ! # grid cols
      integer,              intent(in)    :: lchnk         ! chunk index
      type(physics_ptend),  intent(inout) :: ptend         ! grid ptend

      !
      ! Local variables
      !
      integer           :: indx, i, j, l
      integer           :: nsubcols(pcols)

      call subcol_get_nsubcol(lchnk, nsubcols)
      ! physics_ptend_init has already been called by the master interface
      if (ptend%ls) then
         ptend%s(:,:)       = 0._r8
         ptend%hflux_srf(:) = 0._r8
         ptend%hflux_top(:) = 0._r8
         indx = 1
         do i = 1, ngrdcol
            if (docheck) then
               do l=1,nsubcols(i)-1
                 if (any(ptend_sc%s(indx,:) /= ptend_sc%s(indx+l,:))) &
                    call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%s')
                 if (ptend_sc%hflux_srf(indx) /= ptend_sc%hflux_srf(indx+l)) &
                    call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%hflux_srf')
                 if (ptend_sc%hflux_top(indx) /= ptend_sc%hflux_top(indx+l)) &
                    call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%hflux_top')
               end do
            end if
            ptend%s(i,:)          = ptend_sc%s(indx,:)
            ptend%hflux_srf(i)    = ptend_sc%hflux_srf(indx)
            ptend%hflux_top(i)    = ptend_sc%hflux_top(indx)
            indx = indx + nsubcols(i)
         end do
      end if
      if (ptend%lu) then
         ptend%u(:,:)       = 0._r8
         ptend%taux_srf(:)  = 0._r8
         ptend%taux_top(:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            do l=1,nsubcols(i)-1
              if (any(ptend_sc%u(indx,:) /= ptend_sc%u(indx+l,:))) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%u')
              if (ptend_sc%taux_srf(indx) /= ptend_sc%taux_srf(indx+l)) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%taux_srf')
              if (ptend_sc%taux_top(indx) /= ptend_sc%taux_top(indx+l)) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%taux_top')
            end do
            ptend%u(i,:)      = ptend_sc%u(indx,:)
            ptend%taux_srf(i) = ptend_sc%taux_srf(indx)
            ptend%taux_top(i) = ptend_sc%taux_top(indx)
            indx = indx + nsubcols(i)
         end do
      end if
      if (ptend%lv) then
         ptend%v(:,:)       = 0._r8
         ptend%tauy_srf(:)  = 0._r8
         ptend%tauy_top(:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            do l=1,nsubcols(i)-1
              if (any(ptend_sc%v(indx,:) /= ptend_sc%v(indx+l,:))) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%v')
              if (ptend_sc%tauy_srf(indx) /= ptend_sc%tauy_srf(indx+l)) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%tauy_srf')
              if (ptend_sc%tauy_top(indx) /= ptend_sc%tauy_top(indx+l)) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%tauy_top')
            end do
            ptend%v(i,:)         = ptend_sc%v(indx,:)
            ptend%tauy_srf(i)    = ptend_sc%tauy_srf(indx)
            ptend%tauy_top(i)    = ptend_sc%tauy_top(indx)
            indx = indx + nsubcols(i)
         end do
      end if
      if (any(ptend%lq(:))) then
         ptend%q(:,:,:)       = 0._r8
         ptend%cflx_srf(:,:)  = 0._r8
         ptend%cflx_top(:,:)  = 0._r8
         indx = 1
         do i = 1, ngrdcol
            do l=1,nsubcols(i)-1
              if (any(ptend_sc%q(indx,:,:) /= ptend_sc%q(indx+l,:,:))) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%q')
              if (any(ptend_sc%cflx_srf(indx,:) /= ptend_sc%cflx_srf(indx+l,:))) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%cflx_srf')
              if (any(ptend_sc%cflx_top(indx,:) /= ptend_sc%cflx_top(indx+l,:))) &
                   call endrun('subcol_ptend_get_firstsubcol error: Not all subcolumn fields are identical for ptend%cflx_top')
            end do
            ptend%q(i,:,:)         = ptend_sc%q(indx,:,:)
            ptend%cflx_srf(i,:)    = ptend_sc%cflx_srf(indx,:)
            ptend%cflx_top(i,:)    = ptend_sc%cflx_top(indx,:)
            indx = indx + nsubcols(i)
         end do
      end if

   end subroutine subcol_ptend_get_firstsubcol


   ! TYPE int,double,real
   ! DIMS 1,2,3,4,5,6
   subroutine subcol_pack_{DIMS}d_{TYPE}(lchnk, field, field_sc)
   !-----------------------------------------------------------------------
   ! Pack the field defined on (pcols, psubcols, *) into (pcols*psubcols, *)
   ! Packing is done accoding to the values in the proper chunk from nsubcol2d
   !-----------------------------------------------------------------------

      integer, intent(in)                     :: lchnk            ! Chunk index
#if ({DIMS} == 1)
      {VTYPE}, intent(in)                     :: field(:,:)           ! grid
#elif ({DIMS} == 2)
      {VTYPE}, intent(in)                     :: field(:,:,:)         ! grid
#elif ({DIMS} == 3)
      {VTYPE}, intent(in)                     :: field(:,:,:,:)       ! grid
#elif ({DIMS} == 4)
      {VTYPE}, intent(in)                     :: field(:,:,:,:,:)     ! grid
#elif ({DIMS} == 5)
      {VTYPE}, intent(in)                     :: field(:,:,:,:,:,:)   ! grid
#elif ({DIMS} == 6)
      {VTYPE}, intent(in)                     :: field(:,:,:,:,:,:,:) ! grid
#endif
      {VTYPE}, intent(out)                    :: field_sc{DIMSTR}     ! subcols

      !
      ! Local variables
      !
      integer           :: indx, i, j
      integer           :: nsubcol(pcols)

      call subcol_get_nsubcol(lchnk, nsubcol)
      indx = 1
      do i=1, pcols
         do j = 1, nsubcol(i)
#if ({DIMS} == 1)
            field_sc(indx) = field(i, j)
#elif ({DIMS} == 2)
            field_sc(indx, :) = field(i, j, :)
#elif ({DIMS} == 3)
            field_sc(indx, :, :) = field(i, j, :, :)
#elif ({DIMS} == 4)
            field_sc(indx, :, :, :) = field(i, j, :, :, :)
#elif ({DIMS} == 5)
            field_sc(indx, :, :, :, :) = field(i, j, :, :, :, :)
#elif ({DIMS} == 6)
            field_sc(indx, :, :, :, :, :) = field(i, j, :, :, :, :, :)
#endif
            indx = indx + 1
         end do
      end do
   end subroutine subcol_pack_{DIMS}d_{TYPE}

   ! TYPE int,double,real
   ! DIMS 1,2,3,4,5,6
   subroutine subcol_unpack_{DIMS}d_{TYPE}(lchnk, field_sc, field, fillvalue)
   !-----------------------------------------------------------------------
   ! UnPack the field defined on (pcols*psubcols, *) into (pcols, psubcols, *)
   ! Unpacking is done accoding to the values in the proper chunk from nsubcol2d
   ! If fillvalue is present, unused entries in field are set.
   ! NB: The output field is not initialized, if fillvalue is not passed, it
   !     will end up with undefined values for columns where nsubcol < psubcols
   !-----------------------------------------------------------------------

      integer, intent(in)                     :: lchnk            ! Chunk index
      {VTYPE}, intent(in)                     :: field_sc{DIMSTR}     ! subcols
#if ({DIMS} == 1)
      {VTYPE}, intent(out)                    :: field(:,:)           ! grid
#elif ({DIMS} == 2)
      {VTYPE}, intent(out)                    :: field(:,:,:)         ! grid
#elif ({DIMS} == 3)
      {VTYPE}, intent(out)                    :: field(:,:,:,:)       ! grid
#elif ({DIMS} == 4)
      {VTYPE}, intent(out)                    :: field(:,:,:,:,:)     ! grid
#elif ({DIMS} == 5)
      {VTYPE}, intent(out)                    :: field(:,:,:,:,:,:)   ! grid
#elif ({DIMS} == 6)
      {VTYPE}, intent(out)                    :: field(:,:,:,:,:,:,:) ! grid
#endif
      {VTYPE}, intent(in), optional           :: fillvalue            ! fil

      !
      ! Local variables
      !
      integer           :: indx, i, j
      integer           :: nsubcol(pcols)

      call subcol_get_nsubcol(lchnk, nsubcol)
      indx = 1
      do i=1, pcols
         do j = 1, nsubcol(i)
#if ({DIMS} == 1)
            field(i, j) = field_sc(indx)
#elif ({DIMS} == 2)
            field(i, j, :) = field_sc(indx, :)
#elif ({DIMS} == 3)
            field(i, j, :, :) = field_sc(indx, :, :)
#elif ({DIMS} == 4)
            field(i, j, :, :, :) = field_sc(indx, :, :, :)
#elif ({DIMS} == 5)
            field(i, j, :, :, :, :) = field_sc(indx, :, :, :, :)
#elif ({DIMS} == 6)
            field(i, j, :, :, :, :, :) = field_sc(indx, :, :, :, :, :)
#endif
            indx = indx + 1
         end do
         if (present(fillvalue)) then
            do j = nsubcol(i) + 1, psubcols
#if ({DIMS} == 1)
               field(i, j) = fillvalue
#elif ({DIMS} == 2)
               field(i, j, :) = fillvalue
#elif ({DIMS} == 3)
               field(i, j, :, :) = fillvalue
#elif ({DIMS} == 4)
               field(i, j, :, :, :) = fillvalue
#elif ({DIMS} == 5)
               field(i, j, :, :, :, :) = fillvalue
#elif ({DIMS} == 6)
               field(i, j, :, :, :, :, :) = fillvalue
#endif
            end do
         end if
      end do
   end subroutine subcol_unpack_{DIMS}d_{TYPE}

   ! TYPE int,double,real
   {VTYPE} function subcol_avg_inter_{TYPE}(vals, lchnk, icol, indx1, indx2, usefilter, useweight) result(avgs)
   !------------------------------------------------------------------
   ! This function handles the transformation of the usefilter and useweight logicals to passing the
   ! actual fields which subcol_avg requires based on the values of the logicals
   !------------------------------------------------------------------

      {VTYPE},intent(in)  :: vals(:)
      integer, intent(in) :: lchnk
      integer, intent(in) :: icol
      integer, intent(in) :: indx1
      integer, intent(in) :: indx2
      logical, intent(in) :: usefilter
      logical, intent(in) :: useweight

      integer :: nsubcol(pcols)

      call subcol_get_nsubcol(lchnk, nsubcol)
      if (usefilter .and. useweight) then
         avgs = subcol_avg(vals,nsubcol(icol),filter=filter2d(indx1:indx2,lchnk),weight=weight2d(indx1:indx2,lchnk))
      else if (useweight) then 
         avgs = subcol_avg(vals,nsubcol(icol),weight=weight2d(indx1:indx2,lchnk))
      else if (usefilter) then
         avgs = subcol_avg(vals, nsubcol(icol),filter=filter2d(indx1:indx2,lchnk))
      else 
         avgs = subcol_avg(vals,nsubcol(icol))
      end if

   end function subcol_avg_inter_{TYPE}

   ! TYPE int,double,real
   {VTYPE} function subcol_avg_{TYPE}(vals, nsubcol, filter, weight) result(avgs)
   !------------------------------------------------------------------
   ! This function performs the averaging of subcolumn fields, using the optional &
   ! filters and weights appropriately
   !------------------------------------------------------------------

     {VTYPE},  intent(in)           :: vals(:)
     integer,  intent(in)           :: nsubcol
     integer,  intent(in), optional :: filter(:)
     real(r8), intent(in), optional :: weight(:)

     integer :: icnt
     {VTYPE} :: fillval

     fillval = fillval_{TYPE}

     if (present(filter) .and. present(weight)) then
        if (any(filter==1).and. sum(weight,mask=(filter==1)) /=0 ) then
           avgs = sum(vals*weight, mask=(filter==1)) /  sum(weight, mask=(filter==1))
        else
           avgs = fillval
        end if
     else if (present(weight)) then 
        if (sum(weight) /=0 ) then
           avgs = sum(vals*weight) / sum(weight)
        else
           avgs = fillval
        end if
     else if (present(filter)) then
        if (any(filter==1)) then
           icnt = count(filter==1)
           avgs = sum(vals, mask=(filter==1)) / icnt
        else
           avgs = fillval
        end if
     else if (nsubcol /= 0) then
        avgs = sum(vals) / nsubcol
     else
        avgs = fillval
     end if

   end function subcol_avg_{TYPE}

   subroutine subcol_utils_init_restart(File, hdimids)

     use pio,           only: file_desc_t, pio_int, pio_double
     use cam_pio_utils, only: cam_pio_def_dim, cam_pio_def_var

     ! Dummy arguments
     type(file_desc_t), intent(inout) :: File
     integer,           intent(in)    :: hdimids(:)

     ! Local variable
     integer, allocatable          :: adimids(:)

     call cam_pio_def_var(File, 'NSUBCOL2D', pio_int, hdimids, nsubcol_desc)
     ! Storing filter and weight data even if not being filled by the current
     !     subcolumn generator. While these are 2-D arrays, they don't match
     !     the grid so we will have to recast them as 3-D

     ! We will need a dimid for the subcol dimension
     call cam_pio_def_dim(File, 'psubcols', psubcols, subcol_dimid,           &
          existOK=.true.)
     allocate(adimids(size(hdimids,1) + 1))
     adimids(1:size(hdimids)) = hdimids(:)
     adimids(size(hdimids) + 1) = subcol_dimid

     call cam_pio_def_var(File, 'FILTER2D', pio_int,    adimids, filter2d_desc)
     call cam_pio_def_var(File, 'WEIGHT2D', pio_double, adimids, weight2d_desc)

   end subroutine subcol_utils_init_restart

   subroutine subcol_utils_write_restart(File)
     use cam_grid_support, only: cam_grid_write_dist_array
     use cam_grid_support, only: cam_grid_id, cam_grid_dimensions
     use ppgrid,           only: begchunk, endchunk
     use pio,              only: file_desc_t

     ! Dummy argument
     type(file_desc_t), intent(inout) :: File

     ! Local variables
     integer                          :: c
     integer                          :: adimlens(3)
     integer                          :: fdimlens(3)
     integer                          :: frank
     integer                          :: grid_id
     integer,  allocatable            :: unpacked_i(:,:,:)
     real(r8), allocatable            :: unpacked_r(:,:,:)
     character(len=*), parameter      :: subname = 'SUBCOL_UTILS_WRITE_RESTART'

     ! File dimensions
     grid_id = cam_grid_id('physgrid')
     call cam_grid_dimensions(grid_id, fdimlens(1:2), frank)

     ! Write nsubcol2d
     adimlens(1) = size(nsubcol2d, 1)
     adimlens(2) = endchunk - begchunk + 1
     call cam_grid_write_dist_array(File, grid_id, adimlens(1:2),             &
             fdimlens(1:frank), nsubcol2d, nsubcol_desc)

     ! filter2d and weight2d are 3-D variables
     fdimlens(frank + 1) = psubcols
     frank = frank + 1

     ! Write filter2d
     adimlens(1) = pcols
     adimlens(2) = psubcols
     adimlens(3) = endchunk - begchunk + 1
     if ((pcols * psubcols) /= size(filter2d, 1)) then
       call endrun(trim(subname)//": Unsupported size for FILTER2D")
     end if
     ! Unpack filter2d for proper output
     allocate(unpacked_i(pcols, psubcols, begchunk:endchunk))
     do c = begchunk, endchunk
       call subcol_unpack(c, filter2d(:,c), unpacked_i(:,:,c), 0)
     end do
     call cam_grid_write_dist_array(File, grid_id, adimlens,                  &
             fdimlens(1:frank), unpacked_i, filter2d_desc)
     deallocate(unpacked_i)

     ! Write weight2d
     adimlens(1) = pcols
     adimlens(2) = psubcols
     adimlens(3) = endchunk - begchunk + 1
     if ((pcols * psubcols) /= size(weight2d, 1)) then
       call endrun(trim(subname)//": Unsupported size for WEIGHT2D")
     end if
     ! Unpack weight2d for proper output
     allocate(unpacked_r(pcols, psubcols, begchunk:endchunk))
     do c = begchunk, endchunk
       call subcol_unpack(c, weight2d(:,c), unpacked_r(:,:,c), 0.0_r8)
     end do
     call cam_grid_write_dist_array(File, grid_id, adimlens,                  &
             fdimlens(1:frank), unpacked_r, weight2d_desc)
     deallocate(unpacked_r)

   end subroutine subcol_utils_write_restart

   subroutine subcol_utils_read_restart(File)
     use pio,              only: file_desc_t, pio_inq_varid
     use cam_pio_utils,    only: cam_pio_handle_error
     use cam_grid_support, only: cam_grid_id, cam_grid_read_dist_array
     use cam_grid_support, only: cam_grid_dimensions
     use ppgrid,           only: begchunk, endchunk

     ! Dummy argument
     type(file_desc_t), intent(inout) :: File

     integer                          :: ierr, c
     integer                          :: adimlens(3)
     integer                          :: fdimlens(3)
     integer                          :: grid_id
     integer                          :: frank
     integer,  allocatable            :: unpacked_i(:,:,:)
     real(r8), allocatable            :: unpacked_r(:,:,:)
     character(len=*), parameter      :: subname = 'SUBCOL_UTILS_READ_RESTART'

     call subcol_allocate_internal()

     ! Array dimensions
     adimlens(1) = size(nsubcol2d, 1)
     adimlens(2) = endchunk - begchunk + 1
     ! File dimensions
     grid_id = cam_grid_id('physgrid')
     call cam_grid_dimensions(grid_id, fdimlens(1:2), frank)
     ierr = pio_inq_varid(File, 'NSUBCOL2D', nsubcol_desc)
     call cam_pio_handle_error(ierr, trim(subname)//': NSUBCOL2D not found')
     call cam_grid_read_dist_array(File, grid_id, adimlens(1:2),              &
          fdimlens(1:frank), nsubcol2d, nsubcol_desc)

     ! We need to update indcol2d so set nsubcol2d to itself
     do c = begchunk, endchunk
       call subcol_set_nsubcol(c, pcols, nsubcol2d(:, c))
     end do

     ierr = pio_inq_varid(File, 'FILTER2D', filter2d_desc)
     call cam_pio_handle_error(ierr, trim(subname)//': FILTER2D not found')
     ! Array dimensions
     adimlens(1) = pcols
     adimlens(2) = psubcols
     adimlens(3) = endchunk - begchunk + 1
     if ((pcols * psubcols) /= size(filter2d, 1)) then
       call endrun(trim(subname)//": Unsupported size for FILTER2D")
     end if
     allocate(unpacked_i(pcols, psubcols, begchunk:endchunk))
     ! File dimensions (good for both filter2d and weight2d)
     frank = frank + 1
     fdimlens(frank) = psubcols
     call cam_grid_read_dist_array(File, grid_id, adimlens,                   &
          fdimlens(1:frank), unpacked_i, filter2d_desc)
     ! Pack filter2d for proper output
     do c = begchunk, endchunk
       call subcol_pack(c, unpacked_i(:,:,c), filter2d(:,c))
     end do
     deallocate(unpacked_i)

     ierr = pio_inq_varid(File, 'WEIGHT2D', weight2d_desc)
     adimlens(1) = pcols
     adimlens(2) = psubcols
     adimlens(3) = endchunk - begchunk + 1
     if ((pcols * psubcols) /= size(weight2d, 1)) then
       call endrun(trim(subname)//": Unsupported size for WEIGHT2D")
     end if
     allocate(unpacked_r(pcols, psubcols, begchunk:endchunk))
     call cam_pio_handle_error(ierr, trim(subname)//': WEIGHT2D not found')
     call cam_grid_read_dist_array(File, grid_id, adimlens,                   &
          fdimlens(1:frank), unpacked_r, weight2d_desc)
     ! Pack weight2d for proper output
     do c = begchunk, endchunk
       call subcol_pack(c, unpacked_r(:,:,c), weight2d(:,c))
     end do
     deallocate(unpacked_r)

   end subroutine subcol_utils_read_restart

end module subcol_utils