program modflow_dll_test_5
    
    implicit none

    integer, parameter          :: max_path_len = 256

    character(len=max_path_len) :: config_file_name    ! config (*.run / *.nam / components file)
    integer                     :: num_args            ! #program arguments

    logical, external :: mf_external_timestep_reached
    integer, external :: mf_dll_get_number_of_mnwells    
    
    logical :: endOfSimulation, doTimeLoop, enableStateSave, stssave
    integer           :: num_wells
    integer           :: nribasimtimestep
    real, allocatable :: RibaQdemand(:)
    real, allocatable :: RibaQrealized(:)
    double precision, allocatable :: ribasim_end_of_timestep(:)
    double precision  :: currenttime 
    integer, external :: mf_init_no_sub_pest_moz

    integer           :: i
    integer           :: retVal = 0

    num_args = nargs()
    if (num_args == 2) then
       call getarg(1, config_file_name)
    else
       stop 'provide *.run / *.nam'
    endif

    ! initialize model
    enableStateSave = .true.
    retVal = mf_init_no_sub_pest_moz(config_file_name, endOfSimulation, -1.d0, enableStateSave)
    if (retVal /= 0) then
        stop 'modflow did not initialize correctly, see modflow log file'
    else if (endOfSimulation) then
        stop 'modflow did not initialize correctly, see modflow log file'
    endif

    ! get number of multi-node wells
    ! *** NOTE: at this stage, this is equal to MXWEL2! ***
    num_wells = mf_dll_get_number_of_mnwells()

    allocate(RibaQdemand(num_wells))
    allocate(RibaQRealized(num_wells))
    
!----------------------------------------------
!Info that should come from Ribasim, but for now is hard-coded here
    nribasimtimestep = 24
    allocate(ribasim_end_of_timestep(nribasimtimestep))
    call get_currenttime(currenttime) !this sets the starting time of Ribasim to the starting time of Modflow.
    
    do i = 1, nribasimtimestep
        if (i.eq.1)then
            ribasim_end_of_timestep(i) = currenttime +  15.
        elseif(i.eq.5.or.i.eq.7.or.i.eq.8.or.i.eq.9.or.i.eq.11.or.i.eq.12.or.i.eq.13.or.i.eq.15.or.i.eq.17.or.i.eq.18.or.i.eq.19.or.i.eq.21.or.i.eq.22.or.i.eq.23)Then
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) +  15.
        elseif(i.eq.2.or.i.eq.3.or.i.eq.6.or.i.eq.10.or.i.eq.14.or.i.eq.16.or.i.eq.20.or.i.eq.24)then
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) +  16.           
        elseif(i.eq.4)then
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) +  13.
        endif
        
        RibaQdemand = 0.0
        RibaQrealized = 0.0
        if(i.eq.1)then
            RibaQdemand(1) = 0.0
            RibaQdemand(2) = 0.0
        else
            RibaQdemand(1) = -259200.0
            RibaQdemand(2) = -1296000.0
        endif 
 
!----------------------------------------------
    
        do while (.not. mf_external_timestep_reached(ribasim_end_of_timestep(i)) .and. .not.endOfSimulation)          
            stssave=.true.
            call set_well_values(RibaQdemand,num_wells) !Note: this only has effect in the first Modflow time step of a new Modflow stress period. Qdemand values subsequently apply to that entire stressperiod.
            call mf_update_no_sub_pest_moz(endOfSimulation, doTimeLoop, stssave) 
            call get_well_values(RibaQrealized,num_wells)
        enddo !Loop over Modflow timesteps
        !
        if (i == 5) then    !Note: the Qdemand adjustment intended here can only have effect if i is chosen such that it corresponds with the first timestep of a new Modflow stress period
            call modflow_dll_set_timestep_done(.false., ribasim_end_of_timestep(i-1)) !Sets the current time back one ribasim timestep
            RibaQdemand = RibaQdemand/2.0 !Adjust Qdemand, to check if statesave/restore and subsequent adjustment of Qdemand by Ribasim will work 
            do while (.not. mf_external_timestep_reached(ribasim_end_of_timestep(i)) .and. .not. endOfSimulation)  
                stssave=.true.
                call set_well_values(RibaQdemand,num_wells)
                call mf_update_no_sub_pest_moz(endOfSimulation, doTimeLoop, stssave) 
                call get_well_values(RibaQrealized,num_wells)
            enddo
        endif
        call modflow_dll_set_timestep_done(.true., ribasim_end_of_timestep(i)) ! time step done, we arrived at the end of the current ribasim time step
    
    enddo !Loop over Ribasim Timesteps
    
    call mf_finish_no_sub_pest_moz()

    end program modflow_dll_test_5