program modflow_dll_test_2
    
    implicit none

    integer, external           :: nargs
    integer, parameter          :: max_path_len = 256

    character(len=max_path_len) :: config_file_name    ! config (*.run / *.nam / components file)
    integer                     :: num_args            ! #program arguments
    
    integer           :: num_wells
    integer           :: nribasimtimestep
    real, allocatable :: RibaQdemand(:)
    real, allocatable :: RibaQrealized(:)
    double precision, allocatable :: ribasim_end_of_timestep(:)
    double precision  :: starttime, currenttime 
    
    integer, external :: modflow_dll_get_number_of_wells, modflow_dll_timestep
    logical, external :: mf_end_of_simulation_reached
    logical, external :: mf_external_timestep_reached

    integer           :: i
    integer           :: retVal = 0
    
    ! run the test with tests\test2\Coupled as working directory
    starttime = 57136.0e+00 ! 2015-04-24
    call modflow_dll_init('mnw1.nam', starttime)
    num_wells = modflow_dll_get_number_of_wells()

    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
        
        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. mf_end_of_simulation_reached() .and. retVal == 0)  
            !modflow-ribasim          
            call set_well_values(RibaQdemand) !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.
            retVal = modflow_dll_timestep()
            call get_well_values(RibaQrealized)
        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. mf_end_of_simulation_reached() .and. retVal == 0)  
                call set_well_values(RibaQdemand)
                retVal = modflow_dll_timestep()
                call get_well_values(RibaQrealized)
            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 modflow_dll_finish()

    end program modflow_dll_test_2