! WJS (1-30-12): The following (turning optimization off) is needed as a workaround for an
! xlf compiler bug, at least in IBM XL Fortran for AIX, V12.1 on bluefire
#ifdef CPRIBM
@PROCESS OPT(0)
#endif

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                                                             
!   ncdf_template.F90.in - part of the Community Ice Sheet Model (CISM)  
!                                                              
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!   Copyright (C) 2005-2018
!   CISM contributors - see AUTHORS file for list of contributors
!
!   This file is part of CISM.
!
!   CISM is free software: you can redistribute it and/or modify it
!   under the terms of the Lesser GNU General Public License as published
!   by the Free Software Foundation, either version 3 of the License, or
!   (at your option) any later version.
!
!   CISM is distributed in the hope that it will be useful,
!   but WITHOUT ANY WARRANTY; without even the implied warranty of
!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!   Lesser GNU General Public License for more details.
!
!   You should have received a copy of the Lesser GNU General Public License
!   along with CISM. If not, see <http://www.gnu.org/licenses/>.
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#define NCO outfile%nc
#define NCI infile%nc

!GENVAR_HAVE_AVG!

module NAME_io
  ! template for creating subsystem specific I/O routines
  ! written by Magnus Hagdorn, 2004

  !WHL, Sept. 2023
  ! Moved some 'use' statements to the top to avoid redundant statements that slow the build,
  ! particularly on the Intel compiler

  use DATAMOD
  use glimmer_ncdf
  use glimmer_paramets
  use glimmer_physcon
  use glimmer_scales

  implicit none

  private :: get_xtype, is_enabled, is_enabled_0dint, is_enabled_1dint, &
       is_enabled_2dint, is_enabled_0dreal, is_enabled_1dreal, is_enabled_2dreal, is_enabled_3dreal, &
       initialize_restart_variable_list

  ! List of variables needed for a restart; there is an array of restart variable lists,
  ! one for each ice sheet instance.
  character(glimmer_nc_vars_len), allocatable, save :: restart_variable_list(:)

  interface is_enabled  ! MJH 10/21/13: Interface needed for determining if arrays have been enabled.  See notes below in NAME_io_create.
    module procedure is_enabled_0dint
    module procedure is_enabled_1dint
    module procedure is_enabled_2dint
    module procedure is_enabled_0dreal
    module procedure is_enabled_1dreal
    module procedure is_enabled_2dreal
    module procedure is_enabled_3dreal
  end interface is_enabled

contains

  !*****************************************************************************
  ! netCDF output
  !*****************************************************************************
  subroutine NAME_io_createall(model,data,outfiles)
    ! open all netCDF files for output
    use glimmer_ncio
    implicit none
    type(glide_global_type) :: model
    type(DATATYPE) :: data ! MJH 10/21/13: Making 'data' mandatory.  See notes below in NAME_io_create
    type(glimmer_nc_output),optional,pointer :: outfiles
    
    ! local variables
    type(glimmer_nc_output), pointer :: oc

    if (present(outfiles)) then
       oc => outfiles
    else
       oc=>model%funits%out_first
    end if

    do while(associated(oc))
       call NAME_io_create(oc,model,data)
       oc=>oc%next
    end do
  end subroutine NAME_io_createall

  subroutine NAME_io_writeall(data,model,atend,outfiles,time)
    ! if necessary write to netCDF files
    use glimmer_ncio
    implicit none
    type(DATATYPE) :: data
    type(glide_global_type) :: model
    logical, optional :: atend
    type(glimmer_nc_output),optional,pointer :: outfiles
    real(dp),optional :: time

    ! local variables
    type(glimmer_nc_output), pointer :: oc
    logical :: forcewrite=.false.

    if (present(outfiles)) then
       oc => outfiles
    else
       oc=>model%funits%out_first
    end if

    if (present(atend)) then
       forcewrite = atend
    end if

    do while(associated(oc))
#ifdef HAVE_AVG
       if (oc%do_averages) then
          call NAME_avg_accumulate(oc,data,model)
       end if
#endif
       call glimmer_nc_checkwrite(oc,model,forcewrite,time)
       if (oc%nc%just_processed) then
          ! write standard variables
          call NAME_io_write(oc,data)
#ifdef HAVE_AVG
          if (oc%do_averages) then
             call NAME_avg_reset(oc,data)
          end if
#endif
       end if
       oc=>oc%next
    end do
  end subroutine NAME_io_writeall
  
  subroutine NAME_io_create(outfile,model,data)

    use cism_parallel, only: parallel_type, &
         parallel_def_dim, parallel_inq_dimid, parallel_def_var, parallel_inq_varid, parallel_put_att
    use glimmer_ncio
    use glimmer_map_types
    use glimmer_log
    implicit none
    type(glimmer_nc_output), pointer :: outfile
    type(glide_global_type) :: model
    type(DATATYPE) :: data    ! MJH 10/21/13: Making 'data' mandatory.  See note below

    integer status,varid,pos,model_id

    ! MJH 10/21/13: Local variables needed for checking if a variable is enabled.
    real(dp) :: tavgf
    integer :: up

    !GENVAR_DIMS!

    ! Expanding restart variables: if 'restart' or 'hot' is present, we remove that
    ! word from the variable list, and flip the restartfile flag.  
    ! In CISM 2.0, 'restart' is the preferred name to represent restart variables, 
    ! but 'hot' is supported for backward compatibility.  Thus, we check for both.
    NCO%vars = ' '//trim(adjustl(NCO%vars))//' '  ! Need to maintain a space at beginning and end of list
    ! expanding restart variables
    pos = index(NCO%vars,' restart ') 
    if (pos.ne.0) then
       NCO%vars = NCO%vars(:pos)//NCO%vars(pos+8:)
       NCO%restartfile = .true.
    end if
    pos = index(NCO%vars,' hot ') 
    if (pos.ne.0) then
       NCO%vars = NCO%vars(:pos)//NCO%vars(pos+4:)
       NCO%restartfile = .true.
    end if
    ! Now apply necessary changes if the file is a restart file.
    if (NCO%restartfile) then
       model_id = model%model_id
       if (model_id > size(restart_variable_list)) then
          call write_log('model_id too large in call to NAME_io_create', GM_FATAL)
       end if
       if ((len_trim(NCO%vars) + len_trim(restart_variable_list(model_id)) + 2) >= len(NCO%vars) ) then
          call write_log('Adding restart variables has made the list of output variables too long for file ' // NCO%filename, &
               GM_FATAL)
       else
          ! Expand the restart variable list 
          ! Need to maintain a space at beginning and end of list
          NCO%vars = trim(NCO%vars) // ' ' // trim(restart_variable_list(model_id)) // ' ' ! (a module variable)
          ! Set the xtype to be double (required for an exact restart)
          outfile%default_xtype = NF90_DOUBLE   
       endif
    end if

    ! Convert temp and flwa to versions on stag grid, if needed
    ! Note: this check must occur after restart variables are expanded which happens in glimmer_nc_readparams
    call check_for_tempstag(model%options%whichdycore,NCO)

    ! checking if we need to handle time averages
    pos = index(NCO%vars,AVG_SUFF)
    if (pos.ne.0) then
       outfile%do_averages = .True.
    end if    

    ! Now that the output variable list is finalized, make sure we aren't truncating what the user intends to be output.
    ! Note: this only checks that the text in the variable list does not extend to within one character of the end of the variable.
    ! It does not handle the case where the user exactly fills the allowable length with variables or has a too-long list with more than one space between variable names.
    if ((len_trim(NCO%vars) + 1 ) >= len(NCO%vars)) then 
       call write_log('The list of output variables is too long for file ' // NCO%filename, GM_FATAL)
    endif


    ! MJH, 10/21/13: In the auto-generated code below, the creation of each output variable is wrapped by a check if the data for that 
    !   variable has a size greater than 0.  This is because of recently added checks in glide_types.F90 that don't fully allocate
    !   some variables if certain model options are disabled.  This is to lower memory requirements while running the model.
    !   The reason they have to be allocated with size zero rather than left unallocated is because the data for
    !   some netCDF output variables is defined with math, which causes an error if the operands are unallocated.
    !   Note that if a variable is not created, then it will not be subsequently written to.
    !   Also note that this change requires that data be a mandatory argument to this subroutine.

    ! Some output variables will need tavgf.  The value does not matter, but it must exist.  
    ! Nonetheless, for completeness give it the proper value that it has in NAME_io_write.
    tavgf = outfile%total_time
    if (tavgf.ne.0.d0) then
       tavgf = 1.d0/tavgf
    end if
    ! Similarly, some output variables use the variable up.  Give it value of 0 here.
    up = 0

    !GENVAR_VARDEF!
  end subroutine NAME_io_create

  subroutine NAME_io_write(outfile,data)

    use cism_parallel, only: parallel_type, parallel_inq_varid, distributed_put_var, parallel_put_var
    implicit none
    type(glimmer_nc_output), pointer :: outfile
    ! structure containg output netCDF descriptor
    type(DATATYPE) :: data
    ! the model instance

    ! local variables
    real(dp) :: tavgf
    integer status, varid
    integer up
    type(parallel_type) :: parallel

    parallel = DATAPATH%parallel
     
    tavgf = outfile%total_time
    if (tavgf.ne.0.d0) then
       tavgf = 1.d0/tavgf
    end if

    ! write variables
    !GENVAR_WRITE!

  end subroutine NAME_io_write


  subroutine NAME_add_to_restart_variable_list(vars_to_add, model_id)
    ! This subroutine adds variables to the list of variables needed for a restart for the
    ! given ice sheet instance.
    !
    ! This is a public subroutine that allows other parts of the model to modify the list,
    ! which is a module level variable.   MJH 1/17/2013

    use glimmer_log
    implicit none

    !------------------------------------------------------------------------------------
    ! Subroutine arguments
    !------------------------------------------------------------------------------------
    character(len=*), intent (in) :: vars_to_add  ! list of variable(s) to be added to the list of restart variables 
    integer, intent(in) :: model_id  ! identifier of this ice sheet instance (1 - N, where N is the total number of ice sheet models in this run)

    !------------------------------------------------------------------------------------
    ! Internal variables
    !------------------------------------------------------------------------------------

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

    ! The first time this subroutine is called, allocate and initialize the restart_variable_list
    if (.not. allocated(restart_variable_list)) then
       call initialize_restart_variable_list()
    end if

    if (model_id > size(restart_variable_list)) then
       call write_log('model_id too large in call to NAME_add_to_restart_variable_list', GM_FATAL)
    end if

    ! Add the variables to the list so long as they don't make the list too long.
    if ( (len_trim(restart_variable_list(model_id)) + 1 + len_trim(vars_to_add)) > len(restart_variable_list(model_id))) then
       call write_log('Adding restart variables has made the restart variable list too long.',GM_FATAL)
    else
       restart_variable_list(model_id) = trim(adjustl(restart_variable_list(model_id))) // ' ' // trim(vars_to_add)
       !call write_log('Adding to NAME restart variable list: ' // trim(vars_to_add) )
    endif

  end subroutine NAME_add_to_restart_variable_list


  ! Functions for the interface 'is_enabled'.  These are needed by the auto-generated code in NAME_io_create
  !   to determine if a variable is 'turned on', and should be written.

  function is_enabled_0dint(var)
    integer, intent(in) :: var
    logical :: is_enabled_0dint
    is_enabled_0dint = .true.  ! scalars are always enabled
    return
  end function is_enabled_0dint

  function is_enabled_1dint(var)
    integer, dimension(:), pointer, intent(in) :: var
    logical :: is_enabled_1dint
    if (associated(var)) then
      is_enabled_1dint = .true.
    else
      is_enabled_1dint = .false.
    endif
    return
  end function is_enabled_1dint

  function is_enabled_2dint(var)
    integer, dimension(:,:), pointer, intent(in) :: var
    logical :: is_enabled_2dint
    if (associated(var)) then
      is_enabled_2dint = .true.
    else
      is_enabled_2dint = .false.
    endif
    return
  end function is_enabled_2dint

  function is_enabled_0dreal(var)
    real(dp), intent(in) :: var
    logical :: is_enabled_0dreal
    is_enabled_0dreal = .true.  ! scalars are always enabled
    return
  end function is_enabled_0dreal

  function is_enabled_1dreal(var)
    real(dp), dimension(:), pointer, intent(in) :: var
    logical :: is_enabled_1dreal
    if (associated(var)) then
      is_enabled_1dreal = .true.
    else
      is_enabled_1dreal = .false.
    endif
    return
  end function is_enabled_1dreal

  function is_enabled_2dreal(var)
    real(dp), dimension(:,:), pointer, intent(in) :: var
    logical :: is_enabled_2dreal
    if (associated(var)) then
      is_enabled_2dreal = .true.
    else
      is_enabled_2dreal = .false.
    endif
    return
  end function is_enabled_2dreal

  function is_enabled_3dreal(var)
    real(dp), dimension(:,:,:), pointer, intent(in) :: var
    logical :: is_enabled_3dreal
    if (associated(var)) then
      is_enabled_3dreal = .true.
    else
      is_enabled_3dreal = .false.
    endif
    return
  end function is_enabled_3dreal


  !*****************************************************************************
  ! netCDF input
  !*****************************************************************************  
  subroutine NAME_io_readall(data, model, filetype)
    ! read from netCDF file
    use glimmer_ncio
    implicit none
    type(DATATYPE) :: data
    type(glide_global_type) :: model
    integer, intent(in), optional :: filetype  ! 0 for input, 1 for forcing; defaults to input

    ! local variables
    type(glimmer_nc_input), pointer :: ic
    integer :: filetype_local

    if (present(filetype)) then
      filetype_local = filetype
    else
      filetype_local = 0 ! default to input type
    end if

    if (filetype_local == 0) then
      ic=>model%funits%in_first
    else
      ic=>model%funits%frc_first
    endif
    do while(associated(ic))
       call glimmer_nc_checkread(ic,model)
       if (ic%nc%just_processed) then
          call NAME_io_read(ic,data)
       end if
       ic=>ic%next
    end do

  end subroutine NAME_io_readall


  subroutine NAME_read_forcing(data, model)

    ! Read data from forcing files
    use glimmer_log
    use cism_parallel, only: main_task

    implicit none
    type(DATATYPE) :: data
    type(glide_global_type), intent(inout) :: model

    ! Locals
    type(glimmer_nc_input), pointer :: ic
    integer :: t, t_prev
    real(dp) :: current_forcing_time   ! current time with reference to the forcing file
    real(dp) :: eps                    ! a tolerance to use for stepwise constant forcing
    logical, parameter :: verbose_read_forcing = .false.

    ! Make eps a fraction of the time step.
    eps = model%numerics%tinc * 1.0d-3

    ! read forcing files
    ic=>model%funits%frc_first
    do while(associated(ic))

!       if (main_task .and. verbose_read_forcing) write(6,*) 'possible forcing times', ic%times

       if (ic%read_once) then  ! read once at initialization; do not re-read at runtime

          ic%nc%just_processed = .true.  ! prevent the file from being read
          if (main_task .and. verbose_read_forcing) then
             write(6,*) ' '
             write(6,*) 'In NAME_read_forcing; will not re-read the read_once file ', trim(ic%nc%filename)
          endif

       else  ! not a read_once file

          ic%nc%just_processed = .true. ! until we find an acceptable time, set this to true which will prevent the file from being read.

          ! Compute the current forcing time.
          ! This is the current model time, plus any offset to be consistent with the time in the forcing file,
          !  plus a small number to allow for roundoff error.
          current_forcing_time = model%numerics%time + ic%time_offset + eps

          ! If cycling repeatedly through a subset of the forcing data, make a further correction:
          ! compute the current time relative to time_start_cycle.
          if (ic%nyear_cycle > 0 .and. current_forcing_time > ic%time_start_cycle) then
             current_forcing_time = ic%time_start_cycle &
                  + mod(current_forcing_time - ic%time_start_cycle, real(ic%nyear_cycle,dp))
          endif

          if (main_task .and. verbose_read_forcing) then
             write(6,*) ' '
             write(6,*) 'In NAME_read_forcing, model time + eps =', model%numerics%time + eps
             write(6,*) 'Forcing file nt, time_offset =', ic%nt, ic%time_offset
             write(6,*) 'time_start_cycle, nyear_cycle:', ic%time_start_cycle, ic%nyear_cycle
             write(6,*) 'current forcing time =', current_forcing_time
          endif

          ! Find the time index associated with the previous model time step
          t_prev = 0
          do t = ic%nt, 1, -1  ! look through the time array backwards
             if (ic%times(t) <= current_forcing_time - model%numerics%tinc) then
                t_prev = t
                if (main_task .and. verbose_read_forcing) write(6,*) 'Previous time index =', t_prev
                exit
             end if
          enddo

          ! Find the current time in the file
          do t = ic%nt, 1, -1  ! look through the time array backwards
             if ( ic%times(t) <= current_forcing_time) then
                ! use the largest time that is smaller or equal to the current time (stepwise forcing)
                if (main_task .and. verbose_read_forcing) &
                     write(6,*) 'Largest time less than current forcing time: t, times(t):', t, ic%times(t)

                ! If this time index (t) is larger than the previous index (t_prev), then read a new time slice.
                ! Otherwise, we already have the current slice, and there is nothing new to read.
                if (t > t_prev) then
                   ! Set the desired time to be read
                   ic%current_time = t
                   ic%nc%just_processed = .false.  ! set this to false so file will be read.
                   if (main_task .and. verbose_read_forcing) write(6,*) 'Read new forcing slice: t, times(t) =', t, ic%times(t)
                endif ! t > t_prev

                exit  ! once we find the time, exit the loop
             end if   ! ic%times(t) <= model%numerics%time + eps

          end do  ! if we get to end of loop without exiting, then this file will not be read at this time

       endif  ! read_once file

       ! move on to the next forcing file
       ic=>ic%next

    end do   ! while(associated)

    ! Now that we've updated metadata for each forcing file, actually perform the read.
    ! This call will only read forcing files where just_processed=.false.
    call NAME_io_readall(data, model, filetype=1)

  end subroutine NAME_read_forcing


  subroutine NAME_read_forcing_once(data, model)

    ! Read data from forcing files with read_once = .true.
    ! Read all time slices in a single call and write to arrays with a time index.

    use glimmer_global, only: msg_length
    use glimmer_log
    use cism_parallel, only: main_task, parallel_reduce_sum

    implicit none
    type(DATATYPE) :: data
    type(glide_global_type), intent(inout) :: model

    ! Locals
    type(glimmer_nc_input), pointer :: ic
    integer :: t                         ! time index
    integer :: nx, ny, nt                ! dimension sizes
    real(dp) :: eps                      ! a tolerance to use for stepwise constant forcing
    real(dp) :: global_sum               ! global sum of an input field
    character(len=msg_length) :: message
    logical, parameter :: verbose_read_forcing = .true.

    ! Make eps a fraction of the time step.
    eps = model%numerics%tinc * 1.0d-3

    ! read forcing files
    ic=>model%funits%frc_first
    do while(associated(ic))

       if (ic%read_once) then

          if (main_task .and. verbose_read_forcing) then
             write(6,*) ' '
             write(6,*) 'In NAME_read_forcing_once'
             write(6,*) 'Filename =', trim(ic%nc%filename)
             write(6,*) 'Number of slices =', ic%nt
          endif

          write(message,*) 'Reading', ic%nt, 'slices of file ', trim(ic%nc%filename), ' just once at initialization'
          call write_log(message)

          nt = ic%nt
          ic%nc%vars = ''

          ! Allocate 3D arrays that contain all time slices for each 2D field
          ! Note: Variables with the 'read_once' attribute must be 2D

          !GENVAR_READ_ONCE_ALLOCATE!
          ! Loop over all time slices in the file
          do t = 1, ic%nt

             if (main_task .and. verbose_read_forcing) then
                write(6,*) 'Read new forcing slice: t index, times(t) =', t, ic%times(t)
             endif

             ! Set the desired time to be read
             ic%current_time = t

             ! Read one time slice into the data derived type
             call NAME_io_read(ic,data)

             ! Copy data from this time slice into the 3D array.
             ! Once the fields have been copied, zero them out.
             ! Also increment the string ic%nc%vars.
             ! This string contains a list of field names with a space before and after each name.

             !GENVAR_READ_ONCE_COPY!
          enddo   ! ic%nt

       endif  ! read_once

       if (main_task .and. verbose_read_forcing) then
          write(6,*) 'Final ic%nc%vars = ', trim(ic%nc%vars)
       endif

       ic=>ic%next

    enddo   ! while(associated)

  end subroutine NAME_read_forcing_once


  subroutine NAME_retrieve_forcing(data, model)

    ! Retrieve a single time slice of forcing from arrays that contain all the forcing.
    ! Called repeatedly at runtime, after calling the read_forcing_once subroutine at initialization.

    use glimmer_global, only: msg_length
    use glimmer_log
    use cism_parallel, only: main_task
    implicit none
    type(DATATYPE) :: data
    type(glide_global_type), intent(inout) :: model

    ! Locals
    type(glimmer_nc_input), pointer :: ic
    integer :: t, t_prev
    real(dp) :: current_forcing_time   ! current time with reference to the forcing file
    real(dp) :: eps                    ! a tolerance to use for stepwise constant forcing
    logical :: retrieve_new_slice      ! if true, then retrieve data for this forcing time slice
    integer :: forcing_year            ! year of data from the forcing file
    integer :: this_year               ! current simulation year relative to tstart; starts at 0
    integer :: year1, year2            ! years read from the shuffle file
    real(dp) :: decimal_year           ! decimal part of the current year
    character(len=msg_length) :: message
    logical, parameter :: verbose_read_forcing = .false.

    ! Make eps a fraction of the time step
    eps = model%numerics%tinc * 1.0d-3

    ! read forcing files

    ic=>model%funits%frc_first
    do while(associated(ic))

       if (ic%read_once) then

          retrieve_new_slice = .false.   ! default is to do nothing

          ! Compute the current forcing time.
          ! This is the current model time, plus any offset to be consistent with the time in the forcing file,
          !  plus a small number to allow for roundoff error.
          ! Code adapted from the read_forcing subroutine above

          current_forcing_time = model%numerics%time + ic%time_offset + eps

          ! If cycling repeatedly through a subset of the forcing data, make a further correction:
          ! compute the current time relative to time_start_cycle.
          if (ic%nyear_cycle > 0 .and. current_forcing_time > ic%time_start_cycle) then
             current_forcing_time = ic%time_start_cycle &
                  + mod(current_forcing_time - ic%time_start_cycle, real(ic%nyear_cycle,dp))
          endif

          if (main_task .and. verbose_read_forcing) then
             write(6,*) ' '
             write(6,*) 'In NAME_retrieve_forcing, model time + eps =', model%numerics%time + eps
             write(6,*) 'Filename = ', trim(ic%nc%filename)
             write(6,*) 'Forcing file nt, time_offset =', ic%nt, ic%time_offset
             write(6,*) 'time_start_cycle, nyear_cycle:', ic%time_start_cycle, ic%nyear_cycle
             write(6,*) 'current forcing time =', current_forcing_time
             write(6,*) 'variable list:', trim(ic%nc%vars)
          endif

          ! Optionally, associate the current forcing time with a different date in the forcing file.
          ! This is done by reading a file that associates each simulation year (relative to tstart)
          !  with a year that is read from a 'shuffle file'. The shuffle file typically consists of
          !  consecutive integers (in column 1), followed by years chosen at random from all the years
          !  in the forcing file (in column 2).

          if (trim(ic%shuffle_file) /= '') then  ! shuffle_file exists
             open(unit=11, file=trim(ic%shuffle_file), status='old')
             this_year = int(current_forcing_time - model%numerics%tstart)
             if (main_task .and. verbose_read_forcing) then
                write(6,*) 'shuffle_file = ', trim(ic%shuffle_file)
                write(6,*) 'tstart, this_year =', model%numerics%tstart, this_year
             endif
             forcing_year = 0
             do while (forcing_year == 0)
                read(11,'(i6,i8)') year1, year2
                if (this_year == year1) then
                   forcing_year = year2
                   exit
                endif
             enddo
             close(11)
             decimal_year = current_forcing_time - floor(current_forcing_time)
             current_forcing_time = real(forcing_year,dp) + decimal_year
             if (main_task .and. verbose_read_forcing) then
                write(6,*) 'forcing_year, decimal =', forcing_year, decimal_year
                write(6,*) 'shuffled forcing_time =', current_forcing_time
             endif
          else
             if (main_task .and. verbose_read_forcing) write(6,*) 'no shuffle_file'
          endif  ! shuffle_file exists

          ! Find the time index associated with the previous model time step
          t_prev = 0
          do t = ic%nt, 1, -1  ! look through the time array backwards
             if (ic%times(t) <= current_forcing_time - model%numerics%tinc) then
                t_prev = t
                if (main_task .and. verbose_read_forcing) write(6,*) 'Previous time index =', t_prev
                exit
             end if
          enddo

          ! Find the current time in the file
          do t = ic%nt, 1, -1  ! look through the time array backwards
             if ( ic%times(t) <= current_forcing_time) then
                ! use the largest time that is smaller or equal to the current time (stepwise forcing)
                if (main_task .and. verbose_read_forcing) &
                     write(6,*) 'Largest time less than current forcing time: t, times(t):', t, ic%times(t)
                ! If this time index (t) is larger than the previous index (t_prev), then retrieve a new time slice.
                ! Otherwise, we already have the current slice, and there is nothing new to read.
                if (t > t_prev) then
                   ! Set the desired time to be read
                   ic%current_time = t
                   retrieve_new_slice = .true.
                   if (main_task .and. verbose_read_forcing) write(6,*) 'Retrieve new forcing slice'
                   write(message,*) &
                        'Retrieve slice', t, 'at forcing time', ic%times(t), 'from file ', trim(ic%nc%filename)
                   call write_log(message)
                endif ! t > t_prev

                exit  ! once we find the time, exit the loop
             end if   ! ic%times(t) <= model%numerics%time + eps

          end do  ! if we get to end of loop without exiting, then there is nothing to retrieve at this time

          ! Check whether each potential read_once field is part of this forcing file.
          ! If so, then copy the data for this time slice from the 3D array to the 2D array.

          if (retrieve_new_slice) then

             !GENVAR_READ_ONCE_RETRIEVE!
          endif   ! retrieve_new_slice

       endif   ! read_once

       ! move on to the next forcing file
       ic=>ic%next

    enddo   ! while(associated)

  end subroutine NAME_retrieve_forcing


  subroutine NAME_io_read(infile,data)

    ! read variables from a netCDF file
    use cism_parallel, only: parallel_type, &
         parallel_inq_varid, parallel_get_att, distributed_get_var, parallel_get_var
    use glimmer_log
    implicit none
    type(glimmer_nc_input), pointer :: infile
    ! structure containg output netCDF descriptor
    type(DATATYPE) :: data
    ! the model instance

    ! local variables
    integer status,varid
    integer up
    real(dp) :: scaling_factor
    type(parallel_type) :: parallel

    parallel = DATAPATH%parallel

    ! read variables
    !GENVAR_READ!

  end subroutine NAME_io_read

  subroutine NAME_io_checkdim(infile,model,data)

    ! check if dimension sizes in file match dims of model
    use cism_parallel, only: parallel_type, parallel_inq_dimid, parallel_inquire_dimension
    use glimmer_log
    implicit none
    type(glimmer_nc_input), pointer :: infile
    ! structure containg output netCDF descriptor
    type(glide_global_type) :: model
    type(DATATYPE), optional :: data

    integer status,dimid,dimsize
    character(len=150) message

    ! check dimensions
    !GENVAR_CHECKDIM!
  end subroutine NAME_io_checkdim

  !*****************************************************************************
  ! calculating time averages
  !*****************************************************************************  
#ifdef HAVE_AVG
  subroutine NAME_avg_accumulate(outfile,data,model)

    ! Accumulate time averages for 'tavg' variables
    ! NOTE: This subroutine works for tavg variables that are written to exactly one output file.
    !       If we try to write tavg variables to multiple output files, the accumulated values
    !        will be too large, because this subroutine will be called more than once per time step.
    ! TODO: Write code to check for doubly listed tavg variables and throw a fatal error.

    use cism_parallel, only: parallel_inq_varid
    implicit none
    type(glimmer_nc_output), pointer :: outfile
    ! structure containg output netCDF descriptor
    type(glide_global_type) :: model
    type(DATATYPE) :: data

    ! local variables
    real(dp) :: factor
    integer status, varid

    ! increase total time
    outfile%total_time = outfile%total_time + model%numerics%tinc
    factor = model%numerics%tinc

    !GENVAR_CALCAVG!
  end subroutine NAME_avg_accumulate

  subroutine NAME_avg_reset(outfile,data)

    use cism_parallel, only: parallel_inq_varid
    implicit none
    type(glimmer_nc_output), pointer :: outfile
    ! structure containg output netCDF descriptor
    type(DATATYPE) :: data

    ! local variables
    integer status, varid

    ! reset total time
    outfile%total_time = 0.d0

    !GENVAR_RESETAVG!
  end subroutine NAME_avg_reset
#endif

  !*********************************************************************
  ! some private procedures
  !*********************************************************************

  !> apply default type to be used in netCDF file
  integer function get_xtype(outfile,xtype)
    implicit none
    type(glimmer_nc_output), pointer :: outfile !< derived type holding information about output file
    integer, intent(in) :: xtype                !< the external netCDF type

    get_xtype = xtype
    
    if (xtype.eq.NF90_REAL .and. outfile%default_xtype.eq.NF90_DOUBLE) then
       get_xtype = NF90_DOUBLE
    end if
    if (xtype.eq.NF90_DOUBLE .and. outfile%default_xtype.eq.NF90_REAL) then
       get_xtype = NF90_REAL
    end if
  end function get_xtype

  !> allocate and initialize the restart_variable_list
  subroutine initialize_restart_variable_list()
    use glide_model_registry, only : get_max_models

    allocate(restart_variable_list(get_max_models()))
    restart_variable_list(:) = ''
  end subroutine initialize_restart_variable_list

  !*********************************************************************
  ! lots of accessor subroutines follow
  !*********************************************************************
  !GENVAR_ACCESSORS!

end module NAME_io
