module dynVarTimeInterpMod

#include "shr_assert.h" 

  !---------------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Contains a derived type and associated methods that extend the base class,
  ! dyn_var_type. The type defined here is for variables that SHOULD be interpolated in
  ! time. For variables of this type, the data have the value given on year Y of the file
  ! at midnight on Jan. 1 at the start of year Y. The value then gets linearly
  ! interpolated over the year, so that by Dec. 31 of year Y, the value is close to the
  ! file's value at year Y+1. Before the start of the time series, the data are fixed at
  ! their value from the first year in the file; after the end of the time series, the
  ! data are fixed at their value from the last year in the file. If the last year on the
  ! file is X, then the data are fixed at this last value (and thus do not vary)
  ! throughout year X.
  !
  ! !USES:
  use shr_kind_mod   , only : r8 => shr_kind_r8
  use dynVarMod      , only : dyn_var_type

  implicit none
  private
  !
  ! !PUBLIC TYPES:
  public :: dyn_var_time_interp_type

  type, extends(dyn_var_type) :: dyn_var_time_interp_type
     private
     ! Note that data are stored as 1-d, then converted to the appropriate dimensionality
     ! as needed
     real(r8), allocatable :: data_at_tlower(:)  ! data at time time_index_lower
     real(r8), allocatable :: data_at_tupper(:)  ! data at time time_index_upper

     integer :: time_index_lower  ! current time_index_lower curresponding to data_at_tlower
     integer :: time_index_upper  ! current time_index_upper curresponding to data_at_tupper
   contains
     generic :: get_current_data => &      ! Get the current value of the data
          get_current_data_1d, get_current_data_2d
     procedure :: get_current_data_1d  ! Get the current value of the data, for a 1-d variable
     procedure :: get_current_data_2d  ! Get the current value of the data, for a 2-d variable
     procedure :: read_data_if_needed  ! Read the next time slice of data, if necessary
  end type dyn_var_time_interp_type

  interface dyn_var_time_interp_type
     module procedure constructor  ! initialize a new dyn_var_time_interp object
  end interface dyn_var_time_interp_type

contains

  ! ======================================================================
  ! Constructors
  ! ======================================================================

  !-----------------------------------------------------------------------
  type(dyn_var_time_interp_type) function constructor( &
       dyn_file, varname, dim1name, conversion_factor, &
       do_check_sums_equal_1, data_shape)
    !
    ! !DESCRIPTION:
    ! Creates an object of type dyn_var_time_interp_type. This also reads the first
    ! set of data.
    !
    ! Assumes that dyn_file has already been initialized.
    !
    ! !USES:
    use dynFileMod     , only : dyn_file_type
    use dynTimeInfoMod , only : time_info_type
    !
    ! !ARGUMENTS:
    type(dyn_file_type), target, intent(in) :: dyn_file ! file containing this variable
    character(len=*), intent(in) :: varname             ! variable name on file
    character(len=*), intent(in) :: dim1name            ! dim1name on file
    real(r8), intent(in) :: conversion_factor           ! data are DIVIDED by conversion_factor immediately after reading them

     ! Only relevant for 2-d variables: should we check to make sure that all sums equal 1?
    logical, intent(in) :: do_check_sums_equal_1
    
    ! Shape of data; max number of dimensions is given by dyn_var_max_dims in dynVarMod.
    ! First dimension is the spatial dimension.
    integer, intent(in) :: data_shape(:)
    !-----------------------------------------------------------------------

    ! Set metadata common to all dyn_var_type objects
    call constructor%set_metadata( &
         dyn_file=dyn_file, &
         varname=varname, &
         dim1name=dim1name, &
         conversion_factor=conversion_factor, &
         do_check_sums_equal_1=do_check_sums_equal_1, &
         data_shape=data_shape)

    ! Allocate space for data
    allocate(constructor%data_at_tlower(product(data_shape)))
    allocate(constructor%data_at_tupper(product(data_shape)))

    ! Read first set of data
    constructor%time_index_lower = dyn_file%time_info%get_time_index_lower()
    constructor%time_index_upper = dyn_file%time_info%get_time_index_upper()
    call constructor%read_variable(constructor%time_index_lower, constructor%data_at_tlower)
    call constructor%read_variable(constructor%time_index_upper, constructor%data_at_tupper)
  end function constructor

  ! ======================================================================
  ! Public methods
  ! ======================================================================

  ! The following specific procedures are NOT actually public, but they can be accessed
  ! via the generic type-bound procedure, get_current_data

  ! DIMS 1,2
  !-----------------------------------------------------------------------
  subroutine get_current_data_{DIMS}d(this, cur_data)
    !
    ! !DESCRIPTION:
    ! Get the current, interpolated value of the data, in cur_data. cur_data should have
    ! the same dimensionality as the underlying data, as given by the data_shape argument
    ! that was passed to the constructor.
    ! 
    ! If necessary, new data are read from the file.
    !
    ! Should be called once per time step, AFTER calling set_current_year on the
    ! underlying dyn_file variable
    !
    ! !USES:
    use dynFileMod     , only : dyn_file_type
    !
    ! !ARGUMENTS:
    class(dyn_var_time_interp_type) , intent(inout) :: this             ! this object
    real(r8)                        , intent(out)   :: cur_data{DIMSTR} ! current value of data
    !
    ! !LOCAL VARIABLES:
    type(dyn_file_type), pointer :: dyn_file ! the dyn_file of this object
    integer               :: ndims          ! ndims of data in 'this'
    real(r8)              :: wt1            ! weight of time1 (the left-hand time point)
    real(r8), allocatable :: cur_data_1d(:) ! 1-d version of data at the current time
    character(len=*), parameter :: subname = 'get_current_data_{DIMS}d'
    !-----------------------------------------------------------------------

    ! Do some error checking
    ndims = size(this%get_data_shape())
    SHR_ASSERT({DIMS} == ndims, subname//' ERROR: # dims of output argument must match ndims')
    SHR_ASSERT_ALL((shape(cur_data) == this%get_data_shape()), subname//' ERROR: shape of cur_data must match shape of data')

    ! Get current data, using a temporal weighting of the data at time 1 and the data at
    ! time 2
    call this%read_data_if_needed()
    allocate(cur_data_1d(size(this%data_at_tlower)))

    dyn_file => this%get_dyn_file()
    wt1 = 1.0_r8 - dyn_file%time_info%get_yearfrac()

    cur_data_1d(:) = this%data_at_tupper(:) + wt1*(this%data_at_tlower(:) - this%data_at_tupper(:))
    cur_data = reshape(cur_data_1d, shape(cur_data))
    deallocate(cur_data_1d)

  end subroutine get_current_data_{DIMS}d

  ! ======================================================================
  ! Private methods
  ! ======================================================================

  !-----------------------------------------------------------------------
  subroutine read_data_if_needed(this)
    !
    ! !DESCRIPTION:
    ! Determine if new data need to be read from the file; if so, read them.
    ! 
    ! We need to read new data (or at least copy them) if the current time on dyn_file
    ! disagrees with the time for which we currently have stored data, for either time time_index_lower
    ! or time time_index_upper
    !
    ! !USES:
    use dynFileMod     , only : dyn_file_type
    use dynTimeInfoMod , only : time_info_type
    !
    ! !ARGUMENTS:
    class(dyn_var_time_interp_type), intent(inout) :: this   ! this object
    !
    ! !LOCAL VARIABLES:
    type(dyn_file_type), pointer :: dyn_file ! the dyn_file of this object
    integer :: time_index_lower_cur                       ! current value of time_index_lower on dyn_fileb
    integer :: time_index_upper_cur                       ! current value of time_index_upper on dyn_fileb
    
    character(len=*), parameter :: subname = 'read_data_if_needed'
    !-----------------------------------------------------------------------

    dyn_file => this%get_dyn_file()
    time_index_lower_cur = dyn_file%time_info%get_time_index_lower()
    time_index_upper_cur = dyn_file%time_info%get_time_index_upper()

    ! If time_index_lower time has changed, get a new set of data for time time_index_lower
    if (time_index_lower_cur /= this%time_index_lower) then

       ! The typical case is that we have moved forward by a single time; thus we can
       ! avoid an extra read by simply setting the new data at time_index_lower equal to the old data
       ! at time_index_upper
       if (time_index_lower_cur == this%time_index_upper) then
          this%data_at_tlower(:) = this%data_at_tupper(:)

       ! Otherwise, handle the general (but atypical) case where we have not moved
       ! forward by a single time
       else
          call this%read_variable(time_index_lower_cur, this%data_at_tlower)
       end if

       this%time_index_lower = time_index_lower_cur
    end if

    ! If time_index_upper time has changed, read a new set of data for time time_index_upper
    if (time_index_upper_cur /= this%time_index_upper) then
       call this%read_variable(time_index_upper_cur, this%data_at_tupper)
       this%time_index_upper = time_index_upper_cur
    end if

  end subroutine read_data_if_needed


end module dynVarTimeInterpMod