! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program 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
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.

module writeOutput_module

! NetCDF types
USE netcdf
USE netcdf_util_module,only:netcdf_err                    ! netcdf error handling function

! top-level data types
USE nrtype

! missing values
USE globalData,only: integerMissing, realMissing

! provide access to global data
USE globalData,only:gru_struc                             ! gru->hru mapping structure

USE globalData,only:outputStructure

USE data_types,only:var_i

! provide access to the derived types to define the data structures
USE data_types,only:&
                    ! final data vectors
                    dlength,             & ! var%dat
                    ilength,             & ! var%dat
                    time_dlength,        & ! var(:)%tim(:)%dat (dp)
                    ! no spatial dimension
                    var_i,               & ! x%var(:)            (i4b)
                    var_i8,              & ! x%var(:)            integer(8)
                    var_d,               & ! x%var(:)            (dp)
                    var_ilength,         & ! x%var(:)%dat        (i4b)
                    var_dlength,         & ! x%var(:)%dat        (dp)
                    ! no variable dimension
                    hru_i,               & ! x%hru(:)            (i4b)
                    hru_d,               & ! x%hru(:)            (dp)
                    time_i,              &
                    ! gru dimension
                    gru_int,             & ! x%gru(:)%var(:)     (i4b)
                    gru_double,          & ! x%gru(:)%var(:)     (dp)
                    gru_intVec,          & ! x%gru(:)%var(:)%dat (i4b)
                    gru_doubleVec,       & ! x%gru(:)%var(:)%dat (dp)
                    ! gru+hru dimension
                    gru_hru_int,         & ! x%gru(:)%hru(:)%var(:)     (i4b)
                    gru_hru_int8,        & ! x%gru(:)%hru(:)%var(:)     integer(8)
                    gru_hru_double,      & ! x%gru(:)%hru(:)%var(:)     (dp)
                    gru_hru_intVec,      & ! x%gru(:)%hru(:)%var(:)%dat (i4b)
                    gru_hru_doubleVec,   & ! x%gru(:)%hru(:)%var(:)%dat (dp)
                    gru_hru_time_double, &
                    gru_hru_time_doubleVec,&
                    gru_hru_time_intVec

! vector lengths
USE var_lookup, only: maxvarFreq ! number of output frequencies
USE var_lookup, only: maxvarStat ! number of statistics

implicit none
private
public::writeParm
public::writeData
public::writeBasin
public::writeTime
public::writeRestart
! define dimension lengths
integer(i4b),parameter      :: maxSpectral=2              ! maximum number of spectral bands
contains

    ! **********************************************************************************************************
    ! public subroutine writeParm: write model parameters
    ! **********************************************************************************************************
subroutine writeParm(ncid,ispatial,struct,meta,err,message)
  USE data_types,only:var_info                    ! metadata info
  USE var_lookup,only:iLookStat                   ! index in statistics vector
  USE var_lookup,only:iLookFreq                   ! index in vector of model output frequencies
  implicit none

  ! declare input variables
  type(var_i)   ,intent(in)   :: ncid             ! file ids
  integer(i4b)  ,intent(in)   :: iSpatial         ! hydrologic response unit
  class(*)      ,intent(in)   :: struct           ! data structure
  type(var_info),intent(in)   :: meta(:)          ! metadata structure
  integer(i4b)  ,intent(out)  :: err              ! error code
  character(*)  ,intent(out)  :: message          ! error message
  ! local variables
  integer(i4b)                :: iVar             ! loop through variables

  ! initialize error control
  err=0;message="writeParm/"
  ! loop through local column model parameters
  do iVar = 1,size(meta)

    ! check that the variable is desired
    if (meta(iVar)%statIndex(iLookFREQ%timestep)==integerMissing) cycle

    ! initialize message
    message=trim(message)//trim(meta(iVar)%varName)//'/'

    ! HRU data
    if (iSpatial/=integerMissing) then
      select type (struct)
        class is (var_i)
          err = nf90_put_var(ncid%var(iLookFreq%timestep),meta(iVar)%ncVarID(iLookFreq%timestep),(/struct%var(iVar)/),start=(/iSpatial/),count=(/1/))
        class is (var_i8)
          err = nf90_put_var(ncid%var(iLookFreq%timestep),meta(iVar)%ncVarID(iLookFreq%timestep),(/struct%var(iVar)/),start=(/iSpatial/),count=(/1/))
        class is (var_d)
          err = nf90_put_var(ncid%var(iLookFreq%timestep),meta(iVar)%ncVarID(iLookFreq%timestep),(/struct%var(iVar)/),start=(/iSpatial/),count=(/1/))
        class is (var_dlength)
          print*, "Param size", size(struct%var(iVar)%dat)
          err = nf90_put_var(ncid%var(iLookFreq%timestep),meta(iVar)%ncVarID(iLookFreq%timestep),(/struct%var(iVar)%dat/),start=(/iSpatial,1/),count=(/1,size(struct%var(iVar)%dat)/))
        class default; err=20; message=trim(message)//'unknown variable type (with HRU)'; return
      end select
      call netcdf_err(err,message); if (err/=0) return

      ! GRU data
    else
      select type (struct)
        class is (var_d)
          err = nf90_put_var(ncid%var(iLookFreq%timestep),meta(iVar)%ncVarID(iLookFreq%timestep),(/struct%var(iVar)/),start=(/1/),count=(/1/))
        class is (var_i8)
          err = nf90_put_var(ncid%var(iLookFreq%timestep),meta(iVar)%ncVarID(iLookFreq%timestep),(/struct%var(iVar)/),start=(/1/),count=(/1/))
        class default; err=20; message=trim(message)//'unknown variable type (no HRU)'; return
      end select
    end if
    call netcdf_err(err,message); if (err/=0) return

    ! re-initialize message
    message="writeParm/"
  end do  ! looping through local column model parameters

end subroutine writeParm

    ! **************************************************************************************
    ! public subroutine writeData: write model time-dependent data
    ! **************************************************************************************
subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,nSteps, &
            minGRU, maxGRU, numGRU, & 
            meta,stat,dat,structName,map,indx,err,message)
  USE data_types,only:var_info                       ! metadata type
  USE var_lookup,only:maxVarStat                     ! index into stats structure
  USE var_lookup,only:iLookVarType                   ! index into type structure
  USE var_lookup,only:iLookIndex                     ! index into index structure
  USE var_lookup,only:iLookStat                      ! index into stat structure
  USE globalData,only:outFreq                        ! output file information
  USE globalData,only:outputStructure
  USE globalData,only:failedHRUs
  USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
  USE get_ixName_module,only:get_statName            ! to access type strings for error messages

  implicit none
  ! declare dummy variables
  type(var_i)   ,intent(in)        :: ncid              ! file ids
  integer(i4b)  ,intent(inout)     :: outputTimestep(:) ! output time step
  integer(i4b)  ,intent(inout)     :: outputTimestepUpdate(:) ! number of HRUs in the run domain
  integer(i4b)  ,intent(in)        :: maxLayers         ! maximum number of layers
  integer(i4b)  ,intent(in)        :: nSteps            ! number of timeSteps
  integer(i4b)  ,intent(in)        :: minGRU            ! minGRU index to write
  integer(i4b)  ,intent(in)        :: maxGRU            ! maxGRU index to write - probably not needed
  integer(i4b)  ,intent(in)        :: numGRU            ! number of GRUs to write 
  type(var_info),intent(in)        :: meta(:)           ! meta data
  class(*)      ,intent(in)        :: stat              ! stats data
  class(*)      ,intent(in)        :: dat               ! timestep data
  character(*)  ,intent(in)        :: structName
  integer(i4b)  ,intent(in)        :: map(:)            ! map into stats child struct
  type(gru_hru_time_intVec) ,intent(in) :: indx         ! index data
  integer(i4b)  ,intent(out)       :: err               ! error code
  character(*)  ,intent(out)       :: message           ! error message
  ! local variables
  integer(i4b)                     :: iVar              ! variable index
  integer(i4b)                     :: iStat             ! statistics index
  integer(i4b)                     :: iFreq             ! frequency index
  integer(i4b)                     :: ncVarID           ! used only for time
  integer(i4b)                     :: nSnow             ! number of snow layers
  integer(i4b)                     :: nSoil             ! number of soil layers
  integer(i4b)                     :: nLayers           ! total number of layers
  ! output arrays
  integer(i4b)                     :: datLength         ! length of each data vector
  integer(i4b)                     :: maxLength         ! maximum length of each data vector
  real(rkind)                      :: timeVec(nSteps)   ! timeVal to copy
  real(rkind)                      :: realVec(numGRU, nSteps)   ! real vector for all HRUs in the run domain
  real(rkind)                      :: realArray(nSteps,maxLayers+1)  ! real array for all HRUs in the run domain
  integer(i4b)                     :: intArray(nSteps,maxLayers+1)   ! integer array for all HRUs in the run domain
  integer(i4b)                     :: dataType          ! type of data
  integer(i4b),parameter           :: ixInteger=1001    ! named variable for integer
  integer(i4b),parameter           :: ixReal=1002       ! named variable for real
  integer(i4b)                     :: stepCounter       ! counter to know how much data we have to write, needed because we do not always write nSteps
  integer(i4b)                     :: gruCounter
  integer(i4b)                     :: iStep
  integer(i4b)                     :: iGRU
  integer(i4b)                     :: verifiedGRUIndex    ! index of HRU verified to not have failed
  ! initialize error control
  err=0;message="writeData/"
  ! loop through output frequencies
  do iFreq=1,maxvarFreq
    ! skip frequencies that are not needed
    if(.not.outFreq(iFreq)) cycle

    ! loop through model variables
    do iVar = 1,size(meta)
      stepCounter = 0

      if (meta(iVar)%varName=='time' .and. structName == 'forc')then
        ! get variable index
        err = nf90_inq_varid(ncid%var(iFreq),trim(meta(iVar)%varName),ncVarID)
        call netcdf_err(err,message); if (err/=0) return
        
        ! make sure the HRU we are using has not failed
        if (minGRU == maxGRU)then
          verifiedGRUIndex = minGRU
        else 
          do iGRU = minGRU, maxGRU
            if(.not.failedHRUs(iGRU))then
              verifiedGRUIndex = iGRU
              exit
            endif
          end do  
        endif

        do iStep = 1, nSteps
          ! check if we want this timestep
          if(.not.outputStructure(1)%finalizeStats(1)%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
          stepCounter = stepCounter+1
          timeVec(stepCounter) = outputStructure(1)%forcStruct(1)%gru(verifiedGRUIndex)%hru(1)%var(iVar)%tim(iStep)
        end do ! iStep

        err = nf90_put_var(ncid%var(iFreq),ncVarID,timeVec(1:stepCounter),start=(/outputTimestep(iFreq)/),count=(/stepCounter/))
        call netcdf_err(err,message); if (err/=0)then; print*, "err"; return; endif
        ! save the value of the number of steps to update outputTimestep at the end of the function
        outputTimeStepUpdate(iFreq) = stepCounter
        cycle
      end if  ! id time

      ! define the statistics index
      iStat = meta(iVar)%statIndex(iFreq)
      ! check that the variable is desired
      if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle

        ! do iHRU=1,gru_struc(iGRU)%hruCount
          ! stats output: only scalar variable type
          if(meta(iVar)%varType==iLookVarType%scalarv) then
            select type(stat)
              class is (gru_hru_time_doubleVec)
                if (minGRU == maxGRU)then
                  gruCounter = 1
                  do iStep = 1, nSteps
                    if(.not.outputStructure(1)%finalizeStats(1)%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
                    stepCounter = stepCounter + 1
                    realVec(gruCounter, stepCounter) = stat%gru(verifiedGRUIndex)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
                  end do ! iStep
                else 
                  gruCounter = 0
                  do iGRU = minGRU, maxGRU
                    stepCounter = 0
                    gruCounter = gruCounter + 1
                    do iStep = 1, nSteps
                      if(.not.outputStructure(1)%finalizeStats(1)%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
                      stepCounter = stepCounter + 1
                      realVec(gruCounter, stepCounter) = stat%gru(iGRU)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
                    end do ! iStep
                  end do ! iGRU  
                endif
                
                err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec(1:gruCounter, 1:stepCounter),start=(/minGRU,outputTimestep(iFreq)/),count=(/numGRU,stepCounter/))
                if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
                  print*, "ERROR Missmatch in Steps - stat doubleVec"
                  print*, "iFreq = ", iFreq
                  print*, "outputTimeStepUpdate(iFreq) = ", outputTimeStepUpdate(iFreq)
                  print*, "stepCounter = ", stepCounter
                  return
                endif
              class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return
            end select  ! stat
          
          ! non-scalar variables: regular data structures
          else

            ! initialize the data vectors
            select type (dat)
              class is (gru_hru_time_doubleVec); realArray(:,:) = realMissing;    dataType=ixReal
              class is (gru_hru_time_intVec);     intArray(:,:) = integerMissing; dataType=ixInteger
              class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
            end select


            ! get the model layers
            nSoil   = indx%gru(iGRU)%hru(1)%var(iLookIndex%nSoil)%tim(iStep)%dat(1)
            nSnow   = indx%gru(iGRU)%hru(1)%var(iLookIndex%nSnow)%tim(iStep)%dat(1)
            nLayers = indx%gru(iGRU)%hru(1)%var(iLookIndex%nLayers)%tim(iStep)%dat(1)

            ! get the length of each data vector
            select case (meta(iVar)%varType)
                case(iLookVarType%wLength); datLength = maxSpectral
                case(iLookVarType%midToto); datLength = nLayers
                case(iLookVarType%midSnow); datLength = nSnow
                case(iLookVarType%midSoil); datLength = nSoil
                case(iLookVarType%ifcToto); datLength = nLayers+1
                case(iLookVarType%ifcSnow); datLength = nSnow+1
                case(iLookVarType%ifcSoil); datLength = nSoil+1
                case default; cycle
            end select ! vartype

            ! get the data vectors
            select type (dat)
                class is (gru_hru_time_doubleVec)
                  do iStep = 1, nSteps
                    if(.not.outputStructure(1)%finalizeStats(1)%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
                    stepCounter = stepCounter + 1
                    realArray(stepCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(:)
                  end do

                class is (gru_hru_time_intVec)
                  do iStep = 1, nSteps
                    if(.not.outputStructure(1)%finalizeStats(1)%gru(verifiedGRUIndex)%hru(1)%tim(iStep)%dat(iFreq)) cycle
                    stepCounter = stepCounter + 1
                    intArray(stepCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(:)
                  end do
                class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
            end select

            ! get the maximum length of each data vector
            select case (meta(iVar)%varType)
              case(iLookVarType%wLength); maxLength = maxSpectral
              case(iLookVarType%midToto); maxLength = maxLayers
              case(iLookVarType%midSnow); maxLength = maxLayers-nSoil
              case(iLookVarType%midSoil); maxLength = nSoil
              case(iLookVarType%ifcToto); maxLength = maxLayers+1
              case(iLookVarType%ifcSnow); maxLength = (maxLayers-nSoil)+1
              case(iLookVarType%ifcSoil); maxLength = nSoil+1
              case default; cycle
            end select ! vartype

            ! write the data vectors
            select case(dataType)
              case(ixReal)

                err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:stepCounter,:),start=(/iGRU,1,outputTimestep(iFreq)/),count=(/1,maxLength,stepCounter/))
                if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
                  print*, "ERROR Missmatch in Steps - ixReal"
                  return
                endif
              case(ixInteger)

                err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),intArray(1:stepCounter,:),start=(/iGRU,1,outputTimestep(iFreq)/),count=(/1,maxLength,stepCounter/))
                if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
                  print*, "ERROR Missmatch in Steps - ixInteger"
                  return
                endif
              case default; err=20; message=trim(message)//'data must be of type integer or real'; return
            end select ! data type

          end if ! not scalarv
        ! end do ! HRU Loop

      ! process error code
      if (err/=0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
      call netcdf_err(err,message); if (err/=0) return

    end do ! iVar
    ! outputTimeStep(iFreq) = outputTimeStep(iFreq) + outputTimeStepUpdateVal(iFreq) 
  end do ! iFreq

end subroutine writeData

! **************************************************************************************
! public subroutine writeBasin: write basin-average variables
! **************************************************************************************
subroutine writeBasin(ncid,iGRU,outputTimestep,iStep,meta,stat,dat,map,err,message)
  USE data_types,only:var_info                       ! metadata type
  USE var_lookup,only:maxVarStat                     ! index into stats structure
  USE var_lookup,only:iLookVarType                   ! index into type structure
  USE globalData,only:outFreq                        ! output file information
  USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
  USE get_ixName_module,only:get_statName            ! to access type strings for error messages
  implicit none

  ! declare dummy variables
  type(var_i)   ,intent(in)     :: ncid              ! file ids
  integer(i4b)  ,intent(in)     :: iGRU              ! GRU index
  integer(i4b)  ,intent(inout)  :: outputTimestep(:) ! output time step
  integer(i4b)  ,intent(in)     :: iStep            ! number of steps in forcing file
  type(var_info),intent(in)     :: meta(:)           ! meta data
  type(time_dlength),intent(in) :: stat(:)           ! stats data
  type(time_dlength),intent(in) :: dat(:)            ! timestep data
  integer(i4b)  ,intent(in)     :: map(:)            ! map into stats child struct
  integer(i4b)  ,intent(out)    :: err               ! error code
  character(*)  ,intent(out)    :: message           ! error message
  ! local variables
  integer(i4b)                  :: iVar              ! variable index
  integer(i4b)                  :: iStat             ! statistics index
  integer(i4b)                  :: iFreq             ! frequency index
  ! initialize error control
  err=0;message="f-writeBasin/"

  ! loop through output frequencies
  do iFreq=1,maxvarFreq

    ! skip frequencies that are not needed
    if(.not.outFreq(iFreq)) cycle

    ! check that we have finalized statistics for a given frequency
    if(.not.outputStructure(1)%finalizeStats(1)%gru(1)%hru(1)%tim(iStep)%dat(iFreq)) cycle

    ! loop through model variables
    do iVar = 1,size(meta)

      ! define the statistics index
      iStat = meta(iVar)%statIndex(iFreq)

      ! check that the variable is desired
      if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle

      ! stats/data output - select data type
      select case (meta(iVar)%varType)

        case (iLookVarType%scalarv)
          err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(/stat(map(iVar))%tim(iStep)%dat(iFreq)/),start=(/iGRU,outputTimestep(iFreq)/),count=(/1,1/))

        case (iLookVarType%routing)
          if (iFreq==1 .and. outputTimestep(iFreq)==1) then
            err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),(/dat(iVar)%tim(iStep)%dat/),start=(/1/),count=(/1000/))
          end if
        case default
          err=40; message=trim(message)//"unknownVariableType[name='"//trim(meta(iVar)%varName)//"';type='"//trim(get_varTypeName(meta(iVar)%varType))//    "']"; return
      end select ! variable type

      ! process error code
      if (err.ne.0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
      call netcdf_err(err,message); if (err/=0) return

    end do ! iVar
  end do ! iFreq

end subroutine writeBasin

    ! **************************************************************************************
    ! public subroutine writeTime: write current time to all files
    ! **************************************************************************************
subroutine writeTime(ncid,outputTimestep,iStep,meta,dat,err,message)
USE data_types,only:var_info                       ! metadata type
USE var_lookup,only:iLookStat                      ! index into stat structure
implicit none

! declare dummy variables
type(var_i)   ,intent(in)     :: ncid              ! file ids
integer(i4b)  ,intent(inout)  :: outputTimestep(:) ! output time step
integer(i4b)  ,intent(in)     :: iStep
type(var_info),intent(in)     :: meta(:)           ! meta data
type(time_i)  ,intent(in)     :: dat(:)            ! timestep data
integer(i4b)  ,intent(out)    :: err               ! error code
character(*)  ,intent(out)    :: message           ! error message
! local variables
integer(i4b)                  :: iVar              ! variable index
integer(i4b)                  :: iFreq             ! frequency index
integer(i4b)                  :: ncVarID           ! used only for time
! initialize error control
err=0;message="f-writeTime/"
! loop through output frequencies
do iFreq=1,maxvarFreq

  ! check that we have finalized statistics for a given frequency
  if(.not.outputStructure(1)%finalizeStats(1)%gru(1)%hru(1)%tim(iStep)%dat(iFreq)) cycle

  ! loop through model variables
  do iVar = 1,size(meta)

    ! check instantaneous
  if (meta(iVar)%statIndex(iFreq)/=iLookStat%inst) cycle
    print*, "Time Data", dat(iVar)%tim(iStep)
      ! get variable id in file
    err = nf90_inq_varid(ncid%var(iFreq),trim(meta(iVar)%varName),ncVarID)
    if (err/=0) message=trim(message)//trim(meta(iVar)%varName); call netcdf_err(err,message)
    if (err/=0) then; err=20; return; end if

    ! add to file
    err = nf90_put_var(ncid%var(iFreq),ncVarID,(/dat(iVar)%tim(iStep)/),start=(/outputTimestep(iFreq)/),count=(/1/))
    if (err/=0) message=trim(message)//trim(meta(iVar)%varName);call netcdf_err(err,message)
    if (err/=0) then; err=20; return; end if

  end do ! iVar
end do ! iFreq


end subroutine writeTime

    ! *********************************************************************************************************
    ! public subroutine printRestartFile: print a re-start file
    ! *********************************************************************************************************
    subroutine writeRestart(filename,         & ! intent(in): name of restart file
                            nGRU,             & ! intent(in): number of GRUs
                            nHRU,             & ! intent(in): number of HRUs
                            prog_meta,        & ! intent(in): prognostics metadata
                            prog_data,        & ! intent(in): prognostics data
                            bvar_meta,        & ! intent(in): basin (gru) variable metadata
                            bvar_data,        & ! intent(in): basin (gru) variable data
                            maxLayers,        & ! intent(in): maximum number of layers
                            maxSnowLayers,    & ! intent(in): maximum number of snow layers
                            indx_meta,        & ! intent(in): index metadata
                            indx_data,        & ! intent(in): index data
                            err,message)        ! intent(out): error control
    ! --------------------------------------------------------------------------------------------------------
    ! --------------------------------------------------------------------------------------------------------
    ! access the derived types to define the data structures
    USE data_types,only:var_info               ! metadata
    ! access named variables defining elements in the data structures
    USE var_lookup,only:iLookINDEX             ! named variables for structure elements
    USE var_lookup,only:iLookVarType           ! named variables for structure elements
    USE var_lookup,only:iLookBVAR              ! named variables for structure elements
    ! constants
    USE globalData,only:gru_struc              ! gru-hru mapping structures
    ! external routines
    USE netcdf_util_module,only:nc_file_close  ! close netcdf file
    USE netcdf_util_module,only:nc_file_open   ! open netcdf file
    USE globalData,only:nTimeDelay             ! number of timesteps in the time delay histogram
    
    implicit none
    ! --------------------------------------------------------------------------------------------------------
    ! input
    character(len=256),intent(in)      :: filename      ! name of the restart file
    integer(i4b),intent(in)            :: nGRU          ! number of GRUs
    integer(i4b),intent(in)            :: nHRU          ! number of HRUs
    type(var_info),intent(in)          :: prog_meta(:)  ! prognostic variable metadata
    type(gru_hru_doubleVec),intent(in) :: prog_data     ! prognostic vars
    type(var_info),intent(in)          :: bvar_meta(:)  ! basin variable metadata
    type(gru_doubleVec),intent(in)     :: bvar_data     ! basin variables
    type(var_info),intent(in)          :: indx_meta(:)  ! metadata
    type(gru_hru_intVec),intent(in)    :: indx_data     ! indexing vars
    ! output: error control
    integer(i4b),intent(out)           :: err           ! error code
    character(*),intent(out)           :: message       ! error message
    ! --------------------------------------------------------------------------------------------------------
    ! dummy variables
    integer(i4b), intent(in)           :: maxLayers     ! maximum number of total layers
    integer(i4b), intent(in)           :: maxSnowLayers ! maximum number of snow layers

    ! local variables
    integer(i4b)                       :: ncid          ! netcdf file id
    integer(i4b),allocatable           :: ncVarID(:)    ! netcdf variable id
    integer(i4b)                       :: ncSnowID      ! index variable id
    integer(i4b)                       :: ncSoilID      ! index variable id

    integer(i4b)                       :: nSoil         ! number of soil layers
    integer(i4b)                       :: nSnow         ! number of snow layers
    integer(i4b)                       :: maxSnow       ! maximum number of snow layers
    integer(i4b)                       :: maxSoil       ! maximum number of soil layers
    integer(i4b)                       :: nLayers       ! number of total layers
    integer(i4b),parameter             :: nSpectral=2   ! number of spectal bands
    integer(i4b),parameter             :: nScalar=1     ! size of a scalar
    integer(i4b)                       :: nProgVars     ! number of prognostic variables written to state file

    integer(i4b)                       :: hruDimID      ! variable dimension ID
    integer(i4b)                       :: gruDimID      ! variable dimension ID
    integer(i4b)                       :: tdhDimID      ! variable dimension ID
    integer(i4b)                       :: scalDimID     ! variable dimension ID
    integer(i4b)                       :: specDimID     ! variable dimension ID
    integer(i4b)                       :: midSnowDimID  ! variable dimension ID
    integer(i4b)                       :: midSoilDimID  ! variable dimension ID
    integer(i4b)                       :: midTotoDimID  ! variable dimension ID
    integer(i4b)                       :: ifcSnowDimID  ! variable dimension ID
    integer(i4b)                       :: ifcSoilDimID  ! variable dimension ID
    integer(i4b)                       :: ifcTotoDimID  ! variable dimension ID

    character(len=32),parameter        :: hruDimName    ='hru'      ! dimension name for HRUs
    character(len=32),parameter        :: gruDimName    ='gru'      ! dimension name for GRUs
    character(len=32),parameter        :: tdhDimName    ='tdh'      ! dimension name for time-delay basin variables
    character(len=32),parameter        :: scalDimName   ='scalarv'  ! dimension name for scalar data
    character(len=32),parameter        :: specDimName   ='spectral' ! dimension name for spectral bands
    character(len=32),parameter        :: midSnowDimName='midSnow'  ! dimension name for snow-only layers
    character(len=32),parameter        :: midSoilDimName='midSoil'  ! dimension name for soil-only layers
    character(len=32),parameter        :: midTotoDimName='midToto'  ! dimension name for layered varaiables
    character(len=32),parameter        :: ifcSnowDimName='ifcSnow'  ! dimension name for snow-only layers
    character(len=32),parameter        :: ifcSoilDimName='ifcSoil'  ! dimension name for soil-only layers
    character(len=32),parameter        :: ifcTotoDimName='ifcToto'  ! dimension name for layered variables

    integer(i4b)                       :: cHRU          ! count of HRUs
    integer(i4b)                       :: iHRU          ! index of HRUs
    integer(i4b)                       :: iGRU          ! index of GRUs
    integer(i4b)                       :: iVar          ! variable index
    logical(lgt)                       :: okLength      ! flag to check if the vector length is OK
    character(len=256)                 :: cmessage      ! downstream error message
    ! --------------------------------------------------------------------------------------------------------

    ! initialize error control
    err=0; message='writeRestart/'

    ! size of prognostic variable vector
    nProgVars = size(prog_meta)
    allocate(ncVarID(nProgVars+1))     ! include 1 additional basin variable in ID array (possibly more later)

    ! maximum number of soil layers
    maxSoil = gru_struc(1)%hruInfo(1)%nSoil

    ! maximum number of snow layers
    maxSnow = maxSnowLayers
    
    ! create file
    err = nf90_create(trim(filename),nf90_classic_model,ncid)
    message='iCreate[create]'; call netcdf_err(err,message); if(err/=0)return

    ! define dimensions
                err = nf90_def_dim(ncid,trim(hruDimName)    ,nHRU       ,    hruDimID); message='iCreate[hru]'     ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(gruDimName)    ,nGRU       ,    gruDimID); message='iCreate[gru]'     ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(tdhDimName)    ,nTimeDelay ,    tdhDimID); message='iCreate[tdh]'     ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(scalDimName)   ,nScalar    ,   scalDimID); message='iCreate[scalar]'  ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(specDimName)   ,nSpectral  ,   specDimID); message='iCreate[spectral]'; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(midSoilDimName),maxSoil    ,midSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(midTotoDimName),maxLayers  ,midTotoDimID); message='iCreate[midToto]' ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(ifcSoilDimName),maxSoil+1  ,ifcSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return
                err = nf90_def_dim(ncid,trim(ifcTotoDimName),maxLayers+1,ifcTotoDimID); message='iCreate[ifcToto]' ; call netcdf_err(err,message); if(err/=0)return
    if (maxSnow>0) err = nf90_def_dim(ncid,trim(midSnowDimName),maxSnow    ,midSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return
    if (maxSnow>0) err = nf90_def_dim(ncid,trim(ifcSnowDimName),maxSnow+1  ,ifcSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return
    ! re-initialize error control
    err=0; message='writeRestart/'

    ! define prognostic variables
    do iVar = 1,nProgVars
    if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle

    ! define variable
    select case(prog_meta(iVar)%varType)
    case(iLookvarType%scalarv);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,  scalDimID /),ncVarID(iVar))
    case(iLookvarType%wLength);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,  specDimID /),ncVarID(iVar))
    case(iLookvarType%midSoil);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midSoilDimID/),ncVarID(iVar))
    case(iLookvarType%midToto);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midTotoDimID/),ncVarID(iVar))
    case(iLookvarType%ifcSoil);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcSoilDimID/),ncVarID(iVar))
    case(iLookvarType%ifcToto);                err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcTotoDimID/),ncVarID(iVar))
    case(iLookvarType%midSnow); if (maxSnow>0) err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,midSnowDimID/),ncVarID(iVar))
    case(iLookvarType%ifcSnow); if (maxSnow>0) err = nf90_def_var(ncid,trim(prog_meta(iVar)%varname),nf90_double,(/hruDimID,ifcSnowDimID/),ncVarID(iVar))
    end select

    ! check errors
    if(err/=0)then
    message=trim(message)//trim(cmessage)//' [variable '//trim(prog_meta(iVar)%varName)//']'
    return
    end if

    ! add parameter description
    err = nf90_put_att(ncid,ncVarID(iVar),'long_name',trim(prog_meta(iVar)%vardesc))
    call netcdf_err(err,message)

    ! add parameter units
    err = nf90_put_att(ncid,ncVarID(iVar),'units',trim(prog_meta(iVar)%varunit))
    call netcdf_err(err,message)

    end do ! iVar
    
    ! define selected basin variables (derived) -- e.g., hillslope routing
    err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+1))
    err = nf90_put_att(ncid,ncVarID(nProgVars+1),'long_name',trim(bvar_meta(iLookBVAR%routingRunoffFuture)%vardesc));   call netcdf_err(err,message)
    err = nf90_put_att(ncid,ncVarID(nProgVars+1),'units'    ,trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varunit));   call netcdf_err(err,message)

    ! define index variables - snow
    err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),nf90_int,(/hruDimID/),ncSnowID); call netcdf_err(err,message)
    err = nf90_put_att(ncid,ncSnowID,'long_name',trim(indx_meta(iLookIndex%nSnow)%vardesc));           call netcdf_err(err,message)
    err = nf90_put_att(ncid,ncSnowID,'units'    ,trim(indx_meta(iLookIndex%nSnow)%varunit));           call netcdf_err(err,message)

    ! define index variables - soil
    err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSoil)%varName),nf90_int,(/hruDimID/),ncSoilID); call netcdf_err(err,message)
    err = nf90_put_att(ncid,ncSoilID,'long_name',trim(indx_meta(iLookIndex%nSoil)%vardesc));           call netcdf_err(err,message)
    err = nf90_put_att(ncid,ncSoilID,'units'    ,trim(indx_meta(iLookIndex%nSoil)%varunit));           call netcdf_err(err,message)

    ! end definition phase
    err = nf90_enddef(ncid); call netcdf_err(err,message); if (err/=0) return

    ! write variables
    do iGRU = 1,nGRU
    do iHRU = 1,gru_struc(iGRU)%hruCount
    cHRU = gru_struc(iGRU)%hruInfo(iHRU)%hru_ix
    do iVar = 1,size(prog_meta)

    ! excape if this variable is not used
    if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle

    ! actual number of layers
    nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow
    nSoil = gru_struc(iGRU)%hruInfo(iHRU)%nSoil
    nLayers = nSoil + nSnow

    ! check size
    ! NOTE: this may take time that we do not wish to use
    okLength=.true.
    select case (prog_meta(iVar)%varType)
        case(iLookVarType%scalarv);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nScalar  )
        case(iLookVarType%wlength);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSpectral)
        case(iLookVarType%midSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSoil    )
        case(iLookVarType%midToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nLayers  )
        case(iLookVarType%ifcSoil);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSoil+1  )
        case(iLookVarType%ifcToto);              okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nLayers+1)
        case(iLookVarType%midSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSnow    )
        case(iLookVarType%ifcSnow); if (nSnow>0) okLength = (size(prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat) == nSnow+1  )
        case default; err=20; message=trim(message)//'unknown var type'; return
    end select

    ! error check
    if(.not.okLength)then
        message=trim(message)//'bad vector length for variable '//trim(prog_meta(iVar)%varname)
        err=20; return
    endif

    ! write data
    select case (prog_meta(iVar)%varType)
        case(iLookVarType%scalarv);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nScalar  /))
        case(iLookVarType%wlength);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSpectral/))
        case(iLookVarType%midSoil);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSoil    /))
        case(iLookVarType%midToto);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nLayers  /))
        case(iLookVarType%ifcSoil);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSoil+1  /))
        case(iLookVarType%ifcToto);              err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nLayers+1/))
        case(iLookVarType%midSnow); if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSnow    /))
        case(iLookVarType%ifcSnow); if (nSnow>0) err=nf90_put_var(ncid,ncVarID(iVar),(/prog_data%gru(iGRU)%hru(iHRU)%var(iVar)%dat/),start=(/cHRU,1/),count=(/1,nSnow+1  /))
        case default; err=20; message=trim(message)//'unknown var type'; return
    end select

    ! error check
    if (err.ne.0) message=trim(message)//'writing variable:'//trim(prog_meta(iVar)%varName)
    call netcdf_err(err,message); if (err/=0) return
    err=0; message='writeRestart/'

    end do ! iVar loop

    ! write index variables
    err=nf90_put_var(ncid,ncSnowID,(/indx_data%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSnow)%dat/),start=(/cHRU/),count=(/1/))
    err=nf90_put_var(ncid,ncSoilID,(/indx_data%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSoil)%dat/),start=(/cHRU/),count=(/1/))

    end do ! iHRU loop
    
    ! write selected basin variables
    err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingRunoffFuture)%dat/),  start=(/iGRU/),count=(/1,nTimeDelay/))
    
    end do  ! iGRU loop

    ! close file
    call nc_file_close(ncid,err,cmessage)
    if(err/=0)then;message=trim(message)//trim(cmessage);return;end if

    ! cleanup
    deallocate(ncVarID)

    end subroutine writeRestart

end module writeOutput_module