!#define DEBUG 1
module physics_buffer

  !----------------------------------------------------------------------- 
  ! 
  ! Purpose: 
  !    Buffer managment for persistant variables 
  !
  ! Author: J. Edwards
  !
  ! This file is used with genf90.pl to generate buffer.F90  
  ! 
  !-----------------------------------------------------------------------

  use shr_kind_mod,   only: r8 => shr_kind_r8, r4=> shr_kind_r4, i4=> shr_kind_i4
  use ppgrid,         only: pcols, begchunk, endchunk, psubcols
  use cam_logfile,    only: iulog
  use pio,            only: var_desc_t
  use dyn_grid,       only: ptimelevels
  use cam_abortutils, only: endrun
  use buffer,         only: buffer_field_allocate, buffer_field_deallocate, buffer_get_field_ptr, buffer_set_field, &
                            dtype_i4, dtype_r4, dtype_r8, buffer_field_default_type, buffer_field_is_alloc

  implicit none
  private

  ! ngrid_types is a private parameter denoting how many types of a field
  ! (e.g., grid, subcol) a pbuf can hold
  integer, parameter :: ngrid_types = 2 

  ! --col_type parameters -- see pbuf_add_field and pbuf_get_field
  ! -- These indices should start at 1 and increment by 1 so that they
  ! -- may be used in a loop from 1 to ngrid_types

  integer, parameter, public :: col_type_grid    = 1
  integer, parameter, public :: col_type_subcol  = 2

  !
  ! PBUF field suffix strings for different grid types for restart files
  !      NB: There should be ngrid_types entries
  !
  character(len=5), parameter :: field_grid_suff(ngrid_types) = (/ "     ", "_SCOL" /)

  integer :: old_time_idx = 1

  ! The API has strings 'GLOBAL' and 'PHYSPKG' which correspond to these
  ! integer constants if global_allocate_all is false fields allocated with
  ! PHYSPKG persistence are deallocated at the end of each physics time step
  ! and reallocated at the beginning of the next.
  ! If global_allocate_all is true then all fields are allocated at model
  ! start and persist until model completion.

  integer, parameter :: persistence_global  = 1
  integer, parameter :: persistence_physpkg = 2
  logical            :: global_allocate_all = .true.

  ! 
  ! physics_buffer_hdr carries the description info for each buffer, only one header 
  ! is allocated for each field and each chunk of the field points to it. 
  ! It is carried as a linkedlist for initialization only.
  !

  type physics_buffer_hdr
     character(len=16) :: name                    = ''
     integer           :: dtype                   = -1
     integer           :: persistence             = -1
     integer           :: dimsizes(6,ngrid_types) = 0
     type(var_desc_t)  :: vardesc(ngrid_types)  ! var id for restart files
     logical           :: is_copy(ngrid_types)    = .false.
     logical           :: added                   = .false.
     type(physics_buffer_hdr), pointer :: nexthdr => null()

  end type physics_buffer_hdr

  !
  ! The default type for a buffer field is buffer_field_double (r8) since that is the 
  ! type most often used in the model.   The F90 transfer function is used to move
  ! data of other types into and out of the pbuf2d
  ! 

  type physics_buffer_desc
     private
     integer                           :: lchnk
     type(physics_buffer_hdr), pointer :: hdr => null()
     type(buffer_field_default_type)   :: bfg(ngrid_types)
  end type physics_buffer_desc

  interface pbuf_get_field
     ! TYPE int,double,real 
     ! DIMS 1,2,3,4,5
     module procedure get_pbuf1d_field_by_index_{DIMS}d_{TYPE}
     ! TYPE int,double,real 
     ! DIMS 1,2,3,4,5
     module procedure get_pbuf2d_field_by_index_{DIMS}d_{TYPE}
  end interface

  interface pbuf_get_field_restart
     ! TYPE int,double,real 
     module procedure get_pbuf2d_field_restart_{TYPE}
  end interface

  interface pbuf_set_field
     ! TYPE int,double,real 
     ! DIMS 1,2,3,4,5
     module procedure  set_pbuf2d_field_by_index_{DIMS}d_{TYPE}
     ! TYPE int,double,real 
     ! DIMS 1,2,3,4,5
     module procedure  set_pbuf1d_field_by_index_{DIMS}d_{TYPE}
     ! TYPE int,double,real 
     module procedure  set_pbuf1d_field_const_by_index_{TYPE}
     ! TYPE int,double,real 
     module procedure  set_pbuf2d_field_const_by_index_{TYPE}
  end interface

  interface pbuf_add_field
     ! TYPE int,double,real 
     module procedure pbuf_add_field_{TYPE}
  end interface

  public :: pbuf_initialize,      &
            pbuf_defaultopts,     &! set namelist defaults
            pbuf_setopts,         &! set user specified namelist values
            pbuf_init_time,       &! Initialize dyn_time_lvls
            pbuf_old_tim_idx,     &! return the index for the oldest time
            pbuf_update_tim_idx,  &! update the index for the oldest time
            pbuf_col_type_index,  &
            pbuf_get_field_name,  &
            pbuf_get_field,       &
            pbuf_add_field,       &
            pbuf_register_subcol, &
            physics_buffer_desc,  &
            pbuf_get_index,       &
            pbuf_get_chunk,       &
            pbuf_allocate,        &
            pbuf_deallocate,      &
            pbuf_set_field,       &
            pbuf_init_restart,    &
            pbuf_write_restart,   &
            pbuf_read_restart,    &
            dtype_r8, dtype_r4, dtype_i4

   integer, public :: dyn_time_lvls  ! number of time levels in physics buffer (dycore dependent)


  !
  !  Currentpbufflds is incremented in calls to pbuf_add_field and determines the size
  !  of the allocated pbuf2d
  !

  integer :: currentpbufflds=0
  type(physics_buffer_hdr), pointer :: hdrbuffertop => null()

  ! 
  ! Insures that we do not attempt to allocate physics_buffer more than once
  !

  logical :: buffer_initialized =.false.

  !
  ! private pio descriptor for time
  !

  type(var_desc_t) :: timeidx_desc

CONTAINS

  !
  ! Purpose: Return default runtime options
  !

  subroutine pbuf_defaultopts(pbuf_global_allocate_out)
     logical, intent(out), optional :: pbuf_global_allocate_out
     if ( present(pbuf_global_allocate_out) ) then
        pbuf_global_allocate_out = global_allocate_all
     endif
  end subroutine pbuf_defaultopts


  !
  ! Purpose: Set runtime options
  !

  subroutine pbuf_setopts(pbuf_global_allocate_in)
     logical, intent(in), optional :: pbuf_global_allocate_in
     if ( present(pbuf_global_allocate_in) ) then
        global_allocate_all = pbuf_global_allocate_in
     endif
  end subroutine pbuf_setopts

  !
  ! Initialize dyn_time_lvls
  !
  subroutine pbuf_init_time()
     dyn_time_lvls = ptimelevels - 1
  end subroutine pbuf_init_time

  !
  ! Return index of oldest time sample in the physics buffer.
  !

  function pbuf_old_tim_idx()
     integer :: pbuf_old_tim_idx
     pbuf_old_tim_idx = old_time_idx
  end function pbuf_old_tim_idx

  !
  ! Update index of old time sample in the physics buffer.
  !

  subroutine pbuf_update_tim_idx()
     old_time_idx = mod(old_time_idx, dyn_time_lvls) + 1
  end subroutine pbuf_update_tim_idx

  !
  ! pbuf_col_type_index returns an index for use with pbuf calls
  !
  ! * col_type:         is set to col_type_grid (if use_subcol=.false.) or
  !                     col_type_subcol (if (use_subcol=.true.). 
  !

  subroutine pbuf_col_type_index(use_subcol, col_type)

    logical,                  intent(in)  :: use_subcol
    integer,                  intent(out) :: col_type

    if (use_subcol) then
      col_type = col_type_subcol
    else
      col_type = col_type_grid
    end if

  end subroutine pbuf_col_type_index

  ! 
  !  Return a pointer to the current chunks physics_buffer.
  !

  function pbuf_get_field_name(index)
    integer, intent(in) :: index
    character(len=16) :: pbuf_get_field_name
    integer :: i
    type(physics_buffer_hdr), pointer :: hdrbuffer

    hdrbuffer => hdrbuffertop
    do i=2,index
       hdrbuffer=>hdrbuffer%nexthdr
    end do
    pbuf_get_field_name = hdrbuffer%name

  end function pbuf_get_field_name


  function pbuf_get_chunk(pbuf2d, lchnk)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    integer, intent(in) :: lchnk

    type(physics_buffer_desc), pointer :: pbuf_get_chunk(:)

    pbuf_get_chunk => pbuf2d(:,lchnk)


  end function pbuf_get_chunk

  ! 
  !  Return .true. iff pbuf has an allocated grid field
  !

  logical function pbuf_field_has_gridcols(pbuf, index) result(rval)

    type(physics_buffer_desc), pointer                :: pbuf(:)
    integer,                            intent(in)    :: index

    ! If the field is a copy, return .false. even if it is allocated
    if (pbuf(index)%hdr%is_copy(col_type_grid)) then
      rval = .false.
    else
      rval = buffer_field_is_alloc(pbuf(index)%bfg(col_type_grid))
    end if

  end function pbuf_field_has_gridcols

  ! 
  !  Return .true. iff pbuf has an allocated subcolumn field
  !

  logical function pbuf_field_has_subcols(pbuf, index) result(rval)

    type(physics_buffer_desc), pointer                :: pbuf(:)
    integer,                            intent(in)    :: index

    ! If the field is a copy, return .false. even if it is allocated
    if (pbuf(index)%hdr%is_copy(col_type_subcol)) then
      rval = .false.
    else
      rval = buffer_field_is_alloc(pbuf(index)%bfg(col_type_subcol))
    end if

  end function pbuf_field_has_subcols

  ! 
  !  Return .true. iff pbuf has an allocated col_type column field
  !

  logical function pbuf_field_has_col_type(pbuf, index, col_type) result(rval)

    type(physics_buffer_desc), pointer                :: pbuf(:)
    integer,                            intent(in)    :: index
    integer,                            intent(in)    :: col_type

    if (col_type == col_type_grid) then
      rval = pbuf_field_has_gridcols(pbuf, index)
    else if(col_type == col_type_subcol) then
      rval = pbuf_field_has_subcols(pbuf, index)
    else
      call endrun('pbuf_field_has_col_type: Invalid col_type')
    end if

  end function pbuf_field_has_col_type

  !
  !  Initialize the buffer, should be called after all pbuf_add_field calls 
  !  have been completed and should only be called once in a run
  !
  subroutine pbuf_initialize(pbuf2d)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    
    integer :: i, c
    type(physics_buffer_hdr), pointer :: hdrbuffer

    !
    ! Allocate memory
    !
    if(buffer_initialized) return
    ! Allocate at least 1 to avoid unallocated error in ideal physics
    allocate(pbuf2d(max(1,currentpbufflds),begchunk:endchunk)) 
    if(currentpbufflds<1) return

    do c = begchunk, endchunk
      hdrbuffer => hdrbuffertop
      do i = 1, currentpbufflds
        if(.not. hdrbuffer%added) then
          call endrun('pbuf_initialize: pbuf, '//trim(hdrbuffer%name)//', not added')
        end if
        pbuf2d(i,c)%lchnk = c
        pbuf2d(i,c)%hdr => hdrbuffer
        hdrbuffer => hdrbuffer%nexthdr
      end do
    end do
    
    buffer_initialized=.true.
    call pbuf_allocate(pbuf2d, 'global')
#ifdef DEBUG
    call pbuf2d_print(pbuf2d)
#endif
    return
  end subroutine pbuf_initialize


  subroutine pbuf_allocate(pbuf2d, persistence)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    character(len=*), intent(in) :: persistence
    integer :: i
    logical :: allocate_all

    ! allocate_all is used to force allocation of all fields at same time as allocation
    ! for global scope.
    allocate_all = .false.
    if ( global_allocate_all ) then
       if ( persistence == 'global' ) then
          allocate_all = .true.
       else
          return
       end if
    end if

    if(allocate_all) then
       do i=1,currentpbufflds
          select case(pbuf2d(i,begchunk)%hdr%dtype)
          case(TYPEDOUBLE)
             call pbuf_allocate_field_double(pbuf2d, i, dtype_r8)
          case(TYPEREAL)
             call pbuf_allocate_field_real(pbuf2d, i, dtype_r4)
!          case(i8)
!             call pbuf_allocate_field_long(pbuf2d, i)
          case(TYPEINT)
             call pbuf_allocate_field_int(pbuf2d, i, dtype_i4)
          end select
       end do
    else if(persistence .eq. 'physpkg') then
       do i=1,currentpbufflds
          if(pbuf2d(i,begchunk)%hdr%persistence==persistence_physpkg) then
             select case(pbuf2d(i,begchunk)%hdr%dtype)
             case(TYPEDOUBLE)
                call pbuf_allocate_field_double(pbuf2d, i, dtype_r8)
             case(TYPEREAL)
                call pbuf_allocate_field_real(pbuf2d, i, dtype_r4)
!             case(i8)
!                call pbuf_allocate_field_long(pbuf2d, i)
             case(TYPEINT)
                call pbuf_allocate_field_int(pbuf2d, i, dtype_i4)
             end select
          end if
       end do
    else if(persistence .eq. 'global') then
       do i=1,currentpbufflds
          if(pbuf2d(i,begchunk)%hdr%persistence==persistence_global) then

             select case(pbuf2d(i,begchunk)%hdr%dtype)
             case(TYPEDOUBLE)
                call pbuf_allocate_field_double(pbuf2d, i, dtype_r8)
             case(TYPEREAL)
                call pbuf_allocate_field_real(pbuf2d, i, dtype_r4)
!             case(i8)
!                call pbuf_allocate_field_long(pbuf2d, i)
             case(TYPEINT)
                call pbuf_allocate_field_int(pbuf2d, i, dtype_i4)
             end select
          end if
       end do
    end if
  end subroutine pbuf_allocate

  ! TYPE int,double,real 
  subroutine pbuf_allocate_field_{TYPE}(pbuf2d, index, dtype)
    {VTYPE},                            intent(in)    :: dtype
    type(physics_buffer_desc), pointer                :: pbuf2d(:,:)
    integer,                            intent(in)    :: index

    integer,                   pointer                :: dimsizes(:)
    integer                                           :: c, i
    logical                                           :: is_copy

    do i = 1, ngrid_types
       ! Note - dimsizes(:)=0 is special case to indicate "do not allocate" and is not fatal
       ! Since this is called by a single thread, setting the dimsizes pointer to first chunk is okay
       dimsizes => pbuf2d(index,begchunk)%hdr%dimsizes(:,i)
       ! Note - We may have dimsizes /= 0 but still don't allocate if copy
       is_copy = pbuf2d(index,begchunk)%hdr%is_copy(i)

       if (any(dimsizes(:) < 0)) &
          call endrun('pbuf_allocate_field: dimsizes must be greater than 0 for pbuf field '&
                       //trim(pbuf2d(index,begchunk)%hdr%name))

       if (all(dimsizes(:) /= 0) .and. (.not. is_copy)) then
          do c = begchunk, endchunk
             call buffer_field_allocate(pbuf2d(index,c)%bfg(i), dimsizes, dtype)
          end do
       end if
    end do
    
  end subroutine pbuf_allocate_field_{TYPE}

  pure logical function pbuf_do_deallocate(hdr, persistence, col_type) result(rval)
    type(physics_buffer_hdr), pointer            :: hdr
    character(len=*),                 intent(in) :: persistence
    integer,                          intent(in) :: col_type

    if (persistence .eq. 'physpkg') then
      if (hdr%is_copy(col_type)) then
        ! If this is a copied field, it is always deallocated
        rval = .true.
      else if (global_allocate_all) then
        ! This is an all-global run so no deallocation
        rval = .false.
      else if(hdr%persistence == persistence_physpkg) then
        ! This pbuf is a physpkg type, deallocate
        rval = .true.
      else
        ! This pbuf is global, do not deallocate
        rval = .false.
      end if
    else
      ! We are doing a global deallocate, everything must go!
      rval = .true.
    end if
  end function pbuf_do_deallocate

  subroutine pbuf_deallocate(pbuf2d, persistence)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    character(len=*) :: persistence

    integer :: i, j, c

    do i = 1, currentpbufflds
      do j = 1, ngrid_types
        if (pbuf_do_deallocate(pbuf2d(i,begchunk)%hdr, persistence, j)) then
          if (buffer_field_is_alloc(pbuf2d(i,begchunk)%bfg(j))) then
             select case(pbuf2d(i,begchunk)%hdr%dtype)
             case(TYPEDOUBLE)
                do c = begchunk,endchunk
                  call buffer_field_deallocate(pbuf2d(i,c)%bfg(j), dtype_r8)
                end do
             case(TYPEREAL)
                do c = begchunk,endchunk
                  call buffer_field_deallocate(pbuf2d(i,c)%bfg(j), dtype_r4)
                end do
             case(TYPEINT)
                do c = begchunk,endchunk
                  call buffer_field_deallocate(pbuf2d(i,c)%bfg(j), dtype_i4)
                end do
             end select
           end if
        end if
      end do    ! ngrid_types
    end do      ! currentpbufflds
  end subroutine pbuf_deallocate

  ! Find a pbuf header pointer based on the name input
  ! Automatically tacks on an extra header if <name> is not found.
  subroutine find_pbuf_header(name, index, bufptr)
    character(len=*),                 intent(in)  :: name
    type(physics_buffer_hdr), pointer             :: bufptr
    integer,                          intent(out) :: index

    ! Local Variables
    logical                                       :: buf_found

    if(.not. associated(hdrbuffertop)) then
      ! This is the very first pbuf, allocate and set
      allocate(hdrbuffertop)
      currentpbufflds = 1
      hdrbuffertop%name = name
    end if

    bufptr=>hdrbuffertop
    buf_found = .false.
    index = 1
    do while(.not. buf_found)
      if(trim(bufptr%name) == trim(name)) then
        buf_found = .true.
      else
        if (associated(bufptr%nexthdr)) then
          bufptr=>bufptr%nexthdr
          index = index + 1
        else
          ! We ran off the end of the buffers, make a new one for <name>
          ! Sanity check, we should have checked exactly currentpbufflds
          if (index /= currentpbufflds) then
            call endrun("find_pbuf_header: currentpbufflds indexing off")
          end if
          currentpbufflds = currentpbufflds + 1
          index = currentpbufflds
          allocate(bufptr%nexthdr)
          bufptr=>bufptr%nexthdr
          bufptr%name = trim(name)
          buf_found = .true.
        end if
      end if
    end do
  end subroutine find_pbuf_header
  !
  !  Register a field in the pbuf
  !  This should be called from physics register routines.
  !  persistence must be 'global' or 'physpkg' 
  !  dtype can be any of r8, r4, i4 as defined in shr_kinds_mod.F90
  !  col_type is either col_type_grid or col_type_subcol.
  !  If no col_type, then grid field is defined (i.e., dimsizes set)

  subroutine pbuf_register_field_int(name, pname, index, persistence,         &
       dtype, dimsizes, col_type, pbuf_add)

    ! Dummy Arguments
    character(len=*),           intent(in)  :: name
    character(len=*),           intent(in)  :: pname ! name of calling parameterization
    integer,                    intent(out) :: index
    character(len=*), optional, intent(in)  :: persistence
    integer,          optional, intent(in)  :: dtype ! used to differentiate specific calls
    integer,          optional, intent(in)  :: dimsizes(:) ! dimension sizes of grid field
    integer,          optional, intent(in)  :: col_type
    logical,          optional, intent(in)  :: pbuf_add


    ! Local Variables
    type(physics_buffer_hdr), pointer :: bufptr
    integer                           :: dimcnt, col_type_use
    character(len=128)                :: errmsg

    if(buffer_initialized) then
      call endrun('Attempt to register pbuf field after buffer initialized')
    end if

    call find_pbuf_header(name, index, bufptr)

    if (present(persistence)) then
      if(persistence .eq. "global") then
        bufptr%persistence = persistence_global
      else
        bufptr%persistence = persistence_physpkg
      end if
    end if

    if (present(dtype)) then
      bufptr%dtype = dtype
    end if

    if (present(dimsizes)) then
      ! Normally, we only set buffer dimsizes if dimsizes is passed
      dimcnt = size(dimsizes)

      if (present(col_type)) then
        col_type_use = col_type
      else
        col_type_use = col_type_grid
      end if

      ! Only allow up to 5 dimensions. #6 is reserved for subcolumn umpacking
      ! and #7 is for the physics chunk index.
      if (dimcnt > 5) then
        call endrun('pbuf_register_field: Attempt to exceed maximum of 5 dimensions for '//trim(name))
      end if

      ! Assign dimensions dependent on col_types input
      ! Note that dimensions are initialized to zero and set if being used
      if (col_type_use == col_type_grid) then
        ! grid is requested
        bufptr%dimsizes(:,col_type_grid) = 1
        bufptr%dimsizes(1:dimcnt,col_type_grid) = dimsizes
        ! If someone previously registered the subcol, reset those dims
        if (bufptr%dimsizes(1,col_type_subcol) == pcols*psubcols) then
          bufptr%dimsizes(2:dimcnt,col_type_subcol) = dimsizes(2:dimcnt)
        end if
      else if (col_type_use == col_type_subcol) then
        bufptr%dimsizes(:,col_type_subcol) = 1
        bufptr%dimsizes(1,col_type_subcol) = pcols*psubcols
        ! This case should only be used for a pbuf_add_field call adding a subcolumn-only field
        if (dimcnt > 1) then
          bufptr%dimsizes(2:dimcnt,col_type_subcol) = dimsizes(2:dimcnt)
        else if (bufptr%dimsizes(1,col_type_grid) > 0) then
          ! If grid field previously registered or added, update dims
          bufptr%dimsizes(2:,col_type_subcol) = bufptr%dimsizes(2:,col_type_grid)
        end if
      end if
    end if

    if (present(pbuf_add)) then
      if (pbuf_add) then
        ! An add request has been made but we have to make sure we have all info
        if (all(bufptr%dimsizes(:,:) == 0)) then
          write(errmsg, *) 'pbuf_add_field: trying to add field with no dimensions'
          call endrun(errmsg)
        end if
        if (bufptr%dtype < 0) then
          write(errmsg, *) 'pbuf_add_field: trying to add field with no type'
          call endrun(errmsg)
        end if
        if (bufptr%persistence < 0) then
          write(errmsg, *) 'pbuf_add_field: trying to add field with bad persistence'
          call endrun(errmsg)
        end if
        bufptr%added = .true.
      end if
    end if

  end subroutine pbuf_register_field_int

  !
  !  Add a field to the pbuf, this should be called from physics register routines
  !  name is required to be unique, 
  !  persistence must be 'global' or 'physpkg' 
  !  dtype can be any of r8, r4, i4 as defined in shr_kinds_mod.F90
  !  If present, col_type must be either col_type_grid or col_type_subcol

  ! TYPE int,double,real 
  subroutine pbuf_add_field_{TYPE}(name, persistence, dtype, dimsizes, index, col_type)

    character(len=*), intent(in)  :: name, persistence
    {VTYPE}, intent(in)           :: dtype ! used only to differentiate specific calls
    integer, intent(in)           :: dimsizes(:) ! dimension sizes of grid field
    integer, intent(in), optional :: col_type
    integer, intent(out)          :: index

    ! Local Variables
    integer                           :: col_type_use

    if(buffer_initialized) then
      call endrun('Attempt to add pbuf field after buffer initialized')
    end if

    if (present(col_type)) then
      col_type_use = col_type
    else
      col_type_use = col_type_grid
    end if

    call pbuf_register_field_int(trim(name), '', index,                       &
         persistence=persistence, dtype={ITYPE},                              &
         dimsizes=dimsizes, col_type=col_type_use, pbuf_add=.true.)

  end subroutine pbuf_add_field_{TYPE}

  subroutine pbuf_register_subcol(name, pname, index)
    use subcol_utils, only: is_subcol_on

    ! Dummy Arguments
    character(len=*),           intent(in)  :: name
    character(len=*),           intent(in)  :: pname ! name of calling parameterization
    integer,                    intent(out) :: index

    ! Local variables
    integer                                 :: dimsizes(1)

    ! You really should not call this routine if subcolumns are not on
    if (.not. is_subcol_on()) then
      call endrun('pbuf_register_subcol: subcolumns are not active')
    end if
    ! Create and pass dimsizes so that subcolums are registered
    dimsizes(1) = pcols * psubcols
    call pbuf_register_field_int(trim(name), trim(pname), index,              &
         dimsizes=dimsizes, col_type=col_type_subcol)

  end subroutine pbuf_register_subcol

  subroutine pbuf2d_print(pbuf2d)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    call pbuf1d_print(pbuf_get_chunk(pbuf2d,begchunk))


  end subroutine pbuf2d_print

  subroutine pbuf1d_print(pbuf)
    type(physics_buffer_desc), pointer :: pbuf(:)
    integer :: i
    type(physics_buffer_desc), pointer :: pbufPtr

    print *,__FILE__,__LINE__,currentpbufflds,size(pbuf)
    do i=1,currentpbufflds
       pbufPtr => pbuf(i)
       print *,__FILE__,__LINE__,i,trim(pbufPtr%hdr%name),pbufPtr%hdr%dtype,pbufPtr%hdr%persistence,pbufPtr%hdr%dimsizes
    end do

  end subroutine pbuf1d_print
  !
  ! Given a pbuf field name return an integer index to the field.  
  ! This index can be used to retrieve the field and is much faster
  ! than using the name in most cases  
  ! 

  function pbuf_get_index(name, errcode) result(index)
    character(len=*), intent(in)     :: name
    integer, intent(inout), optional :: errcode
    integer :: index
    integer :: i
    type(physics_buffer_hdr), pointer :: bufptr


    bufptr=>hdrbuffertop
    index = -1
    do i=1,currentpbufflds
       if(bufptr%name .eq. name) then
          index=i
          exit
       end if
       bufptr=>bufptr%nexthdr
    end do

    if (present(errcode)) then
        errcode = index 
    else if(index<0) then
       call endrun('Attempt to find undefined name in pbuf '//trim(name))
    end if


  end function pbuf_get_index

  !=========================================================================================

  ! 
  ! Given a pbuf2d chunk and an index return a pointer to a field chunk 
  ! 
  !


  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  subroutine get_pbuf1d_field_by_index_{DIMS}d_{TYPE}(pbuf, index, field, start,kount, col_type, copy_if_needed, errcode)

    ! Get the data based on the col_type which is specified.  If no col_type, then grid field is returned
    use subcol_utils,    only: subcol_field_copy

    type(physics_buffer_desc), pointer:: pbuf(:)
    integer,             intent(in)               :: index
    {VTYPE},    pointer                           :: field{DIMSTR}
    integer,             intent(in),  optional    :: start(:),kount(:)
    integer,             intent(in),  optional    :: col_type
    logical,             intent(in),  optional    :: copy_if_needed
    integer,             intent(out), optional    :: errcode

    integer,    pointer               :: dimsizes(:)
    {VTYPE},    pointer               :: gfield{DIMSTR}
    {VTYPE},    pointer               :: sfield{DIMSTR}
    integer                           :: col_type_use
    integer, allocatable              :: kount_grid(:) ! For use in copy_if_needed
    logical                           :: subset

    ! Copy the generic type to one compatable with the request
    if(index<1 .or. index>size(pbuf)) then
      print *,__FILE__,__LINE__,index
      call endrun('index out of range')
    end if

    ! Default col_type is grid columns
    if (present(col_type)) then
      col_type_use = col_type
    else
      col_type_use = col_type_grid
    end if

    ! If there is an errcode, start with an OK status (zero)
    if (present(errcode)) then
      errcode = 0
    end if

    ! Check whether subset of data requested (default is false)
    subset = .false.
    if (present(start) .and. present(kount)) subset = .true.

    ! Check for ill-formed request
    if ( (present(start) .and. .not. present(kount)) .or.        &
         (.not. present(start) .and. present(kount)) ) then
      call endrun('pbuf_get_field: Both start and kount must be present for '//trim(pbuf(index)%hdr%name))
    end if

    ! See if we need to copy the grid field to subcolumns
    if (present(copy_if_needed)) then
      if (copy_if_needed) then
        if (col_type_use == col_type_subcol) then
          ! If a subcolumn field buffer does not exist, allocate one.
          ! Even if start and kount are being passed, allocate and copy the
          ! full grid field buffer so that a future access will succeed.
          if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_use))) then
            dimsizes => pbuf(index)%hdr%dimsizes(:,col_type_use)
            dimsizes(2:) = pbuf(index)%hdr%dimsizes(2:,col_type_grid)
            dimsizes(1) = pcols * psubcols
            select case(pbuf(index)%hdr%dtype)
            case(TYPEDOUBLE)
              call buffer_field_allocate(pbuf(index)%bfg(col_type_subcol),       &
                   dimsizes, dtype_r8)
            case(TYPEREAL)
              call buffer_field_allocate(pbuf(index)%bfg(col_type_subcol),       &
                   dimsizes, dtype_r4)
            case(TYPEINT)
              call buffer_field_allocate(pbuf(index)%bfg(col_type_subcol),       &
                   dimsizes, dtype_i4)
            end select
            ! Set copy only if we did the allocation after init time
            pbuf(index)%hdr%is_copy(col_type_subcol) = .true.
          end if
        else
          call endrun('pbuf_get_field: copy_if_needed only supported for subcolumns')
        end if
        if (pbuf(index)%hdr%is_copy(col_type_subcol)) then
          ! Only do the copy if we did the alloc (i.e., set the is_copy flag)
          ! Only copy the portion we are going to hand back (i.e., start, kount)
          ! Chances are that kount(1) = pcols*psubcols because we are looking
          ! for a subcolumn field (or we wouldn't be here). Now, 
          ! subcol_field_copy requires kount(1) = pcols for the input and
          ! therefore, kount(1) = pcols*psubcols for the output. Check and
          ! make it work

          if (subset) then

            ! Create kount array for grid field
            ! Use input subcol kount array, replacing the first dimension with pcols
            if (size(kount) > size(pbuf(index)%hdr%dimsizes(:,col_type_subcol))) then
              call endrun('pbuf_get_field: kount input has too many dimensions')
            end if
            if (kount(1) /= pcols * psubcols) then
              call endrun('pbuf_get_field: kount(1) must be pcols*psubcols when using copy_if_needed=.true.')
            endif

            allocate(kount_grid(size(kount)))
            kount_grid(2:) = kount(2:)
            kount_grid(1) = pcols

            ! Don't need to create start array for grid field as start array for subcolumn field is identical
            if (size(start) > size(pbuf(index)%hdr%dimsizes(:,col_type_subcol))) then
              call endrun('pbuf_get_field: start input has too many dimensions')
            end if
            if (start(1) /= 1) then
              call endrun('pbuf_get_field: start(1) must be 1 when using copy_if_needed=.true.')
            end if

            ! Get the grid field
            call buffer_get_field_ptr(pbuf(index)%bfg(col_type_grid),           &
                 gfield, start, kount_grid)

            deallocate(kount_grid)

          else
            ! Get the grid field
            call buffer_get_field_ptr(pbuf(index)%bfg(col_type_grid),           &
                 gfield)
          end if

          ! Get the subcol field pointer (note optional start/kount retain their status in this call)
          call buffer_get_field_ptr(pbuf(index)%bfg(col_type_subcol),            &
               sfield, start, kount)

          call subcol_field_copy(gfield, pbuf(index)%lchnk, sfield)

        end if
      end if
    end if

    ! Copy or not, retrieve the requested field pointer
    if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_use))) then
       if (present(errcode)) then
         errcode = -1
       else
         if (col_type_use == col_type_grid) then
           call endrun('pbuf_get_field: probably missing a pbuf_add_field call.  field not allocated for '&
                        //trim(pbuf(index)%hdr%name))
         else if (col_type_use == col_type_subcol) then
           call endrun('pbuf_get_field: probably missing a pbuf_register_subcol.  field not allocated for '&
                        //trim(pbuf(index)%hdr%name))
         else
           call endrun('pbuf_get_field:  field not allocated for '//trim(pbuf(index)%hdr%name))
         end if
       end if
    else
        ! Get the field pointer (note optional start/kount retain their status in this call)
        call buffer_get_field_ptr(pbuf(index)%bfg(col_type_use),field,start,kount )
    end if

  end subroutine get_pbuf1d_field_by_index_{DIMS}d_{TYPE}

  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  subroutine get_pbuf2d_field_by_index_{DIMS}d_{TYPE}(pbuf2d, lchnk, index, field, start, kount, col_type, errcode)

    ! Get the data based on the col_type which is specified.  If no col_type, then grid field is returned
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    integer, intent(in) :: lchnk
    integer, intent(in) :: index
    integer, intent(in),  optional :: start(:),kount(:)
    integer, intent(in),  optional :: col_type
    integer, intent(out), optional :: errcode

    {VTYPE}, pointer :: field{DIMSTR}

    ! Check for ill-formed request
    if ( (present(start) .and. .not. present(kount)) .or.&
         (.not. present(start) .and. present(kount)) ) then
      call endrun('pbuf_get_field: Both start and kount must be present for '//trim(pbuf2d(index,begchunk)%hdr%name))
    end if

    call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), index, field,           &
         start=start, kount=kount, col_type=col_type, errcode=errcode)
  end subroutine get_pbuf2d_field_by_index_{DIMS}d_{TYPE}

  ! TYPE int,double,real 
  subroutine get_pbuf2d_field_restart_{TYPE}(pbuf2d, lchnk, index, field, mdimsize, col_type)
    use subcol_utils,    only: subcol_unpack

    ! NB: This routine is really only useful for write_restart_field
    ! Get the data based on the col_type which is specified.
    ! If no col_type, then grid field is assumed
    ! For grid field, reference buffer and copy into chunk (field)
    ! If col_type is a subcol, unpack subcolumn data
    ! Field is then reshaped into (/pcols, psubcols*mdims/).

    ! Dummy variables
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    integer, intent(in)                :: lchnk
    integer, intent(in)                :: index
    integer, intent(in)                :: mdimsize
    integer, intent(in)                :: col_type
    {VTYPE}, intent(inout)             :: field(:,:)

    ! Local variables
    {VTYPE},      allocatable          :: exp_fld(:,:,:,:,:,:)
    {VTYPE},                   pointer :: buf(:,:,:,:,:)
    {VTYPE}                            :: fillvalue
    character(len=128)                 :: errmsg
    integer,                   pointer :: dimsizes(:)

#if ({ITYPE}==TYPEREAL)
    fillvalue = 0._r4
#elif ({ITYPE}==TYPEDOUBLE)
    fillvalue = 0._r8
#else
    fillvalue = 0
#endif

    if (col_type == col_type_grid) then
      call pbuf_get_field(pbuf2d, lchnk, index, buf, col_type=col_type)
      field(:,:) = reshape(buf, (/pcols, mdimsize/))
    else if (col_type == col_type_subcol) then
      ! Don't initialize field, unpack will fill in unused slots
      call pbuf_get_field(pbuf2d, lchnk, index, buf, col_type=col_type)
      dimsizes => pbuf2d(index, lchnk)%hdr%dimsizes(:,col_type_subcol)
      allocate(exp_fld(pcols, psubcols, dimsizes(2), dimsizes(3), dimsizes(4), dimsizes(5)))
      ! unpack the subcolumns into their own dimension
      call subcol_unpack(lchnk, buf, exp_fld, fillvalue)
      ! Reshape back to pcols for outputting
      field(:,:) = reshape(exp_fld, (/pcols, mdimsize/))
      deallocate(exp_fld)
    else
      write(errmsg, *) "get_pbuf2d_field_restart_{TYPE}: Bad col_type:",col_type
      call endrun(errmsg)
    end if
  end subroutine get_pbuf2d_field_restart_{TYPE}

  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  subroutine set_pbuf2d_field_const_by_index_{TYPE}(pbuf2d,index,const, col_type)

    ! Set the field specified by the col_type

    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    integer, intent(in)                :: index
    {VTYPE},intent(in)                 :: const
    integer,intent(in)  ,optional      :: col_type

    integer :: c
    
    do c=begchunk,endchunk
       if(present(col_type)) then
          call set_pbuf1d_field_const_by_index_{TYPE}(pbuf_get_chunk(pbuf2d,c),index,const, col_type=col_type)
       else
          call set_pbuf1d_field_const_by_index_{TYPE}(pbuf_get_chunk(pbuf2d,c),index,const)
       end if
    end do

  end subroutine set_pbuf2d_field_const_by_index_{TYPE}

  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  subroutine set_pbuf1d_field_const_by_index_{TYPE}(pbuf,index,const,start,kount, col_type)

    ! Set the field(s) specified by the col_type
    type(physics_buffer_desc), pointer :: pbuf(:)
    integer, intent(in) :: index
    {VTYPE},intent(in) :: const
    integer, intent(in), optional :: start(:),kount(:)
    integer, intent(in), optional :: col_type
   
    integer :: col_type_use
    logical :: subset

    ! Default col_type is grid
    if (present(col_type)) then
      col_type_use = col_type
    else
      col_type_use = col_type_grid
    end if

    ! Check whether subset of data requested (default is false)
    subset = .false.
    if (present(start) .and. present(kount)) subset = .true.

    ! Check for ill-formed request
    if ( (present(start) .and. .not. present(kount)) .or.&
         (.not. present(start) .and. present(kount)) ) then
       call endrun('pbuf_set_field: Both start and kount must be present for '//trim(pbuf(index)%hdr%name))
    end if
    
    ! Set the appropriate grid or sub-column field.  Check that the fields have been allocated.
    if(subset) then 

       if (col_type_use == col_type_subcol) then
         ! Set sub-column field
         if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_subcol))) &
              call endrun('pbuf_set_field: sub-column field not allocated for '//trim(pbuf(index)%hdr%name))
         call buffer_set_field(pbuf(index)%bfg(col_type_subcol),const,start,kount)
       
       else if (col_type_use == col_type_grid) then
         ! Set grid field
         if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_grid))) &
              call endrun('pbuf_set_field: grid field not allocated for '//trim(pbuf(index)%hdr%name))
         call buffer_set_field(pbuf(index)%bfg(col_type_grid),const,start,kount)
       else
         call endrun('pbuf_set_field: Trying to set '//trim(pbuf(index)%hdr%name)//&
              ' but col_type is neither col_type_grid nor col_type_subcol')
       end if

     else 

       if (col_type_use == col_type_subcol) then
         ! Set sub-column field
         if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_subcol))) &
              call endrun('pbuf_set_field: sub-column field not allocated for '//trim(pbuf(index)%hdr%name))
         call buffer_set_field(pbuf(index)%bfg(col_type_subcol),const)

       else if (col_type_use == col_type_grid) then
         ! Set grid field
         if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_grid))) &
              call endrun('pbuf_set_field: grid field not allocated for '//trim(pbuf(index)%hdr%name))
         call buffer_set_field(pbuf(index)%bfg(col_type_grid),const)
       else
         call endrun('pbuf_set_field: Trying to set '//trim(pbuf(index)%hdr%name)//&
              ' but col_type is neither col_type_grid nor col_type_subcol')
       end if

     end if

   end subroutine set_pbuf1d_field_const_by_index_{TYPE}

  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  subroutine set_pbuf2d_field_by_index_{DIMS}d_{TYPE}(pbuf2d,index,field, start, kount, col_type)

    ! Set the field(s) specified by the col_type
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    integer, intent(in) :: index
    integer,intent(in),optional :: start(:), kount(:)
    integer,intent(in),optional :: col_type

    logical :: subset

    integer :: c
    {VTYPE}, pointer :: fld{DIMSTR}


#if ({DIMS}==1)
    {VTYPE},pointer :: field(:,:)
#elif ({DIMS}==2)
    {VTYPE},pointer :: field(:,:,:)
#elif ({DIMS}==3)
    {VTYPE},pointer :: field(:,:,:,:)
#elif ({DIMS}==4)
    {VTYPE},pointer :: field(:,:,:,:,:)
#elif ({DIMS}==5)
    {VTYPE},pointer :: field(:,:,:,:,:,:)
#endif

    ! Check whether subset of data requested (default is false)
    subset = .false.
    if (present(start) .and. present(kount)) subset = .true.

    ! Check for ill-formed request
    if ( (present(start) .and. .not. present(kount)) .or.&
         (.not. present(start) .and. present(kount)) ) then
      call endrun('pbuf_set_field: Both start and kount must be present for '//trim(pbuf2d(index,begchunk)%hdr%name))
    end if

    do c=begchunk,endchunk
      fld => get_field_chunk_{DIMS}d_{TYPE}(field,c)
      if(subset .and. present(col_type)) then
        call pbuf_set_field(pbuf_get_chunk(pbuf2d,c),index,fld,start,kount, col_type)
      else if(subset) then
        call pbuf_set_field(pbuf_get_chunk(pbuf2d,c),index,fld,start,kount)
      else if(present(col_type)) then
        call pbuf_set_field(pbuf_get_chunk(pbuf2d,c),index,fld,col_type=col_type)
      else
        call pbuf_set_field(pbuf_get_chunk(pbuf2d,c),index,fld)
      end if
    end do
  end subroutine set_pbuf2d_field_by_index_{DIMS}d_{TYPE}

  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  function get_field_chunk_{DIMS}d_{TYPE}(field, c) result(fld) 
    ! module private helper function
    {VTYPE}, pointer :: fld{DIMSTR}
    integer, intent(in) :: c

#if ({DIMS}==1)
    {VTYPE},pointer :: field(:,:)
    fld => field(:,c)
#elif ({DIMS}==2)
    {VTYPE},pointer :: field(:,:,:)
    fld => field(:,:,c)
#elif ({DIMS}==3)
    {VTYPE},pointer :: field(:,:,:,:)
    fld => field(:,:,:,c)
#elif ({DIMS}==4)
    {VTYPE},pointer :: field(:,:,:,:,:)
    fld => field(:,:,:,:,c)
#elif ({DIMS}==5)
    {VTYPE},pointer :: field(:,:,:,:,:,:)
    fld => field(:,:,:,:,:,c)
#endif    

  end function get_field_chunk_{DIMS}d_{TYPE}



  ! TYPE int,double,real 
  ! DIMS 1,2,3,4,5
  subroutine set_pbuf1d_field_by_index_{DIMS}d_{TYPE}(pbuf,index,field, start, kount, col_type)

    type(physics_buffer_desc), pointer :: pbuf(:)
    integer, intent(in) :: index
    {VTYPE}, intent(in) :: field{DIMSTR}
    integer,intent(in),optional :: start(:), kount(:)
    integer,intent(in),optional :: col_type

    integer :: col_type_use
    logical :: subset

    ! Default col_type is grid only
    if (present(col_type)) then
      col_type_use = col_type
    else
      col_type_use = col_type_grid
    end if

    ! Check whether subset of data requested (default is false)
    subset = .false.
    if (present(start) .and. present(kount)) subset = .true.

    ! Check for ill-formed request
    if ( (present(start) .and. .not. present(kount)) .or.&
         (.not. present(start) .and. present(kount)) ) then
       call endrun('pbuf_set_field: Both start and kount must be present for '//trim(pbuf(index)%hdr%name))
    end if

    if(subset) then

       ! Set sub-column field
       if (col_type_use == col_type_subcol) then
          if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_subcol))) &
              call endrun('pbuf_set_field: sub-column field not allocated for '//trim(pbuf(index)%hdr%name))
          call buffer_set_field(pbuf(index)%bfg(col_type_subcol),field,start,kount)      

       ! Set grid field
       else if (col_type_use == col_type_grid) then
          if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_grid))) &
                 call endrun('pbuf_set_field: grid field not allocated for '//trim(pbuf(index)%hdr%name))
          call buffer_set_field(pbuf(index)%bfg(col_type_grid),field,start,kount)
       else
         call endrun('pbuf_set_field: Trying to set '//trim(pbuf(index)%hdr%name)//&
                    ' but col_type is neigher col_type_grid nor col_type_subcol ')
       end if
    else 

       ! Set sub-column field
       if (col_type_use == col_type_subcol) then
          if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_subcol))) &
              call endrun('pbuf_set_field: sub-column field not allocated for '//trim(pbuf(index)%hdr%name))
          call buffer_set_field(pbuf(index)%bfg(col_type_subcol),field)

       ! Set grid field
       else if (col_type_use == col_type_grid) then
          if (.not. buffer_field_is_alloc(pbuf(index)%bfg(col_type_grid))) &
                 call endrun('pbuf_set_field: grid field not allocated for '//trim(pbuf(index)%hdr%name))
          call buffer_set_field(pbuf(index)%bfg(col_type_grid),field)
       else
         call endrun('pbuf_set_field: Trying to set '//trim(pbuf(index)%hdr%name)//&
                    ' but col_type is neither col_type_grid nor col_type_subcol ')
       end if
    endif

  end subroutine set_pbuf1d_field_by_index_{DIMS}d_{TYPE}


  function pbuftype2piotype(pbuftype) result(piotype)
    use pio, only : pio_double, pio_int, pio_real

    integer, intent(in) :: pbuftype
    integer :: piotype

    select case(pbuftype)
    case (TYPEDOUBLE)
       piotype = pio_double
    case(TYPEINT)
       piotype = pio_int
    case(TYPEREAL)
       piotype = pio_real
!    case(TYPELONG)
!       piotype = pio_int
    case default
       write(iulog, *) 'Dtype = ', pbuftype
       call endrun('No restart support for dtype')
    end select
  end function pbuftype2piotype


!
! Initialize a restart file to write - all additional dims in a field are
! bundled into a single dimension for output and a dim pbuf_xxxxx is declared
! in the file if it does not already exist.   
!
  subroutine pbuf_init_restart(File, pbuf2d)
    use pio,              only: file_desc_t, pio_int
    use cam_pio_utils,    only: cam_pio_def_dim, cam_pio_def_var
    use phys_grid,        only: phys_decomp
    use cam_grid_support, only: cam_grid_get_file_dimids, cam_grid_dimensions


    ! Dummy Variables
    type(file_desc_t),                 intent(inout) :: file
    type(physics_buffer_desc), pointer               :: pbuf2d(:,:)

    ! Local Variables
    type(physics_buffer_desc), pointer               :: pbuf
    integer                                          :: i, grid_select
    integer                                          :: adimid(3) ! PIO IDs
    integer                                          :: hdimcnt ! # grid dims
    integer                                          :: dimcnt  ! # array dims
    integer                                          :: mdimsize, piodtype

    character(len=10)                                :: dimname
    character(len=24)                                :: varname

    ! Use adimid as a temp to find number of horizontal dims
    call cam_grid_dimensions(phys_decomp, adimid(1:2), hdimcnt)
    call cam_grid_get_file_dimids(phys_decomp, File, adimid)

    do i = 1, currentpbufflds
      pbuf => pbuf2d(i,begchunk)

      ! Only save global pbufs for restart
      if(pbuf%hdr%persistence /= persistence_global) cycle

      piodtype = pbuftype2piotype(pbuf%hdr%dtype)

      do grid_select = 1, ngrid_types
        ! For subcol fields, mdimsize includes psubcols in size
        mdimsize  = product(pbuf%hdr%dimsizes(:,grid_select))/pcols
        if(mdimsize > 1) then
          write(dimname,'(a,i5.5)') 'pbuf_',mdimsize
          call cam_pio_def_dim(File, dimname, mdimsize, adimid(hdimcnt+1),    &
               existOK=.true.)
          dimcnt = hdimcnt + 1
        else
          dimcnt = hdimcnt
        end if
        if (mdimsize > 0) then
          varname = trim(pbuf%hdr%name)//trim(field_grid_suff(grid_select))
          call cam_pio_def_var(File, varname, piodtype, adimid(1:dimcnt),     &
               pbuf%hdr%vardesc(grid_select), existOK=.false.)
        end if
      end do

    end do

    call cam_pio_def_var(File, 'pbuf_time_idx', pio_int, timeidx_desc,        &
         existOK=.false.)

  end subroutine pbuf_init_restart


  subroutine pbuf_write_restart(File, pbuf2d)
    use pio,              only: file_desc_t, pio_put_var
    use cam_grid_support, only: cam_grid_dimensions
    use phys_grid,        only: phys_decomp

    ! Dummy Variables
    type(file_desc_t),                 intent(inout) :: file
    type(physics_buffer_desc), pointer               :: pbuf2d(:,:)

    ! Local Variables
    type(physics_buffer_desc), pointer                :: pbufhdr(:)

    integer                                           :: index, dtype, ierr
    integer                                           :: gdimlens(2)
    integer                                           :: grank

    pbufhdr => pbuf_get_chunk(pbuf2d, begchunk)
    gdimlens = 1
    call cam_grid_dimensions(phys_decomp, gdimlens(1:2), grank)

    do index = 1, currentpbufflds
      if(pbufhdr(index)%hdr%persistence == persistence_global) then
        dtype = pbufhdr(index)%hdr%dtype
        select case(dtype)
        case (TYPEDOUBLE)
          call write_restart_field_double(File, pbuf2d, index, gdimlens, grank)
        case (TYPEREAL)
          call write_restart_field_real(File, pbuf2d, index, gdimlens, grank)
        case (TYPEINT)
          call write_restart_field_int(File, pbuf2d, index, gdimlens, grank)
        end select
      end if
    end do
    ierr = pio_put_var(File, timeidx_desc, (/old_time_idx/))

  end subroutine pbuf_write_restart

  subroutine pbuf_restart_dimsizes(pbuf, gdimlens, grank, gridnum, adimlens,  &
       fdimlens, frank)

    ! Dummy arguments
    type(physics_buffer_desc), pointer    :: pbuf
    integer,                intent(in)    :: gdimlens(2) ! Grid horiz. dim sizes
    integer,                intent(in)    :: grank       ! Array rank on file
    integer,                intent(in)    :: gridnum     ! pbuf grid selector
    integer,                intent(out)   :: adimlens(3) ! Array dims
    integer,                intent(out)   :: fdimlens(3) ! Array dims on file
    integer,                intent(out)   :: frank       ! array rank on file

    ! Local variable
    integer                             :: mdimsize

    mdimsize = PRODUCT(pbuf%hdr%dimsizes(:,gridnum)) / pcols
    fdimlens(1:2) = gdimlens
    if (grank == 1) then
      ! Unstructured grid
      fdimlens(2) = mdimsize
      fdimlens(3) = 1
      frank       = 2
    else
      fdimlens(3) = mdimsize
      frank       = 3
    end if
    ! Restart does not write a dimension if mdimsize == 1
    if (mdimsize == 1) then
      frank = frank - 1
    end if

    adimlens(1) = pcols
    adimlens(2) = mdimsize
    adimlens(3) = endchunk - begchunk + 1

  end subroutine pbuf_restart_dimsizes

  ! TYPE int,double,real 
  subroutine write_restart_field_{TYPE}(File, pbuf2d, index, gdimlens, grank)
    use cam_grid_support, only: cam_grid_write_dist_array
    use phys_grid,        only: phys_decomp
    use pio,              only: file_desc_t

    ! Dummy Variables
    type(file_desc_t),                 intent(inout) :: File
    type(physics_buffer_desc), pointer               :: pbuf2d(:,:)
    integer,                           intent(in)    :: index
    integer,                           intent(in)    :: gdimlens(2)
    integer,                           intent(in)    :: grank

    ! Local Variables
    type(physics_buffer_desc), pointer               :: pbuf
    {VTYPE},       allocatable                       :: field(:,:,:)

    integer                                          :: grid_select ! 1=grid, 2=subcol
    integer                                          :: c           ! Chunk index
    integer                                          :: fdimlens(3) ! Array dimensions on NetCDF file
    integer                                          :: adimlens(3) ! Array dimensions
    integer                                          :: frank

    pbuf => pbuf2d(index, begchunk)

    do grid_select = 1, ngrid_types
      call pbuf_restart_dimsizes(pbuf, gdimlens, grank, grid_select,          &
           adimlens, fdimlens, frank)
      if ((PRODUCT(adimlens) > 0) .and. (.not. pbuf%hdr%is_copy(grid_select))) then
        allocate(field(adimlens(1), adimlens(2), begchunk:endchunk))
        do c = begchunk, endchunk
          call pbuf_get_field_restart(pbuf2d, c, index, field(:,:,c),         &
               adimlens(2), grid_select)
        end do
        if (size(field,2) == 1) then
          ! Special case for 2-D pbuf fields
          adimlens(2) = adimlens(3)
          call cam_grid_write_dist_array(File, phys_decomp, adimlens(1:2),    &
               fdimlens(1:frank), field(:,1,:), pbuf%hdr%vardesc(grid_select))
        else
          call cam_grid_write_dist_array(File, phys_decomp, adimlens,         &
               fdimlens(1:frank), field, pbuf%hdr%vardesc(grid_select))
        end if
        deallocate(field)
      end if
    end do

  end subroutine write_restart_field_{TYPE}


  subroutine pbuf_read_restart(File, pbuf2d)
    use cam_grid_support, only: cam_grid_dimensions
    use phys_grid,        only: phys_decomp
    use pio,              only: file_desc_t, pio_inq_varid, pio_get_var

    ! Dummy Variables
    type(File_desc_t),                 intent(inout) :: File
    type(physics_buffer_desc), pointer               :: pbuf2d(:,:)

    ! Local Variables
    type(physics_buffer_desc), pointer               :: pbufhdr(:)

    integer                                          :: index, dtype, ierr
    integer                                          :: gdimlens(2) ! Horizontal grid dimensions
    integer                                          :: grank

    call pbuf_initialize(pbuf2d)
    
    pbufhdr => pbuf_get_chunk(pbuf2d, begchunk)

    gdimlens = 1
    call cam_grid_dimensions(phys_decomp, gdimlens(1:2), grank)
    ierr = pio_inq_varid(File, 'pbuf_time_idx', timeidx_desc)
    ierr = pio_get_var(File, timeidx_desc, old_time_idx)

    do index = 1, currentpbufflds
      if(pbufhdr(index)%hdr%persistence == persistence_global) then
        dtype = pbufhdr(index)%hdr%dtype
        select case(dtype)
        case (TYPEDOUBLE)
          call read_restart_field_double(File, pbuf2d, index, gdimlens, grank)
        case (TYPEREAL)
          call read_restart_field_real(File, pbuf2d, index, gdimlens, grank)
        case (TYPEINT)
          call read_restart_field_int(File, pbuf2d, index, gdimlens, grank)
        end select
      end if
    end do

  end subroutine pbuf_read_restart

  ! TYPE int,double,real 
  subroutine read_restart_field_{TYPE} (File, pbuf2d, index, gdimlens, grank)
    use pio,              only: file_desc_t
    use pio,              only: pio_inq_varid
    use cam_pio_utils,    only: cam_pio_handle_error
    use cam_grid_support, only: cam_grid_read_dist_array
    use subcol_utils,     only: subcol_pack
    use phys_grid,        only: phys_decomp

    ! Dummy Variables
    type(file_desc_t),                 intent(inout) :: File
    type(physics_buffer_desc), pointer               :: pbuf2d(:,:)
    integer,                           intent(in)    :: index
    integer,                           intent(in)    :: gdimlens(2)
    integer,                           intent(in)    :: grank

    ! Local Variables
    type(physics_buffer_desc), pointer :: pbuf
    {VTYPE},                   pointer :: fld6(:,:,:,:,:,:)
    {VTYPE},       allocatable         :: fld5_pack(:,:,:,:,:)
    {VTYPE},       allocatable         :: field(:,:,:)

    integer                            :: grid_select ! (1=grid, 2=subcol)
    integer                            :: ierr, c
    integer                            :: start(6)
    integer                            :: dimsizes(6)
    integer                            :: fdimlens(3) ! Array dimensions on NetCDF file
    integer                            :: frank       ! Array rank on file
    integer                            :: adimlens(3) ! Array dimensions

    character(len=24)                  :: varname
    character(len=*), parameter        :: subname = 'read_restart_field_{TYPE}'

    pbuf => pbuf2d(index, begchunk)
    start = 1

    do grid_select = 1, ngrid_types

      dimsizes(:) = pbuf%hdr%dimsizes(:, grid_select)
      if(all(dimsizes(:) == 0)) then
        ! None of this grid type for this variable
        cycle
      end if

      call pbuf_restart_dimsizes(pbuf, gdimlens, grank, grid_select,          &
           adimlens, fdimlens, frank)
      ! Fix up dimensions for subcolumn field
      if (grid_select == col_type_subcol) then
        ! Field stored as  (pcols, psubcols*restDims, chunks)
        ! Read then pack to (pcols*psubcols,restDims, chunks)
        allocate(fld5_pack(dimsizes(1),dimsizes(2),dimsizes(3),dimsizes(4),dimsizes(5)))
        do c = 5, 2, -1
          dimsizes(c + 1) = dimsizes(c)
        end do
        dimsizes(1) = pcols
        dimsizes(2) = psubcols
      end if

      varname = trim(pbuf%hdr%name)//trim(field_grid_suff(grid_select))
      ierr = pio_inq_varid(File, varname, pbuf%hdr%vardesc(grid_select))
      call cam_pio_handle_error(ierr, trim(subname)//': '//trim(varname)//' not found')

      allocate(field(adimlens(1), adimlens(2), begchunk:endchunk))
      if (size(field,2) == 1) then
        ! Special case for 2-D pbuf fields
        adimlens(2) = adimlens(3)
        call cam_grid_read_dist_array(File, phys_decomp, adimlens(1:2),       &
             fdimlens(1:frank), field(:,1,:), pbuf%hdr%vardesc(grid_select))
      else
        call cam_grid_read_dist_array(File, phys_decomp, adimlens,            &
             fdimlens(1:frank), field, pbuf%hdr%vardesc(grid_select))
      end if

      do c = begchunk, endchunk
        if (grid_select == col_type_grid) then
          call buffer_get_field_ptr(pbuf2d(index,c)%bfg(grid_select), fld6,   &
               start, dimsizes)
          fld6 = reshape(field(:,:,c), dimsizes)
        else if (grid_select == col_type_subcol) then
          call subcol_pack(c, reshape(field(:,:,c), dimsizes), fld5_pack)
          call buffer_set_field(pbuf2d(index,c)%bfg(grid_select), fld5_pack)
        else
          call endrun('read_restart_field_{TYPE}: invalid grid selector - must be either 1 or 2')
        end if
      end do
      nullify(fld6)
      deallocate(field)
      if (allocated(fld5_pack)) then
        deallocate(fld5_pack)
      end if
    end do

  end subroutine read_restart_field_{TYPE}

end module physics_buffer
