program modflow_dll_test_1
    
    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\test1\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
!Chosen values correspond exactly to the Modflow MNW1 test example in tests/test1, and therefore allows a check whether putting these
!values via this driver program yields the same results as when they are provided using the standard Modflow input files.
    nribasimtimestep = 12
    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.le.2) then
            if(i.eq.1)then
             ribasim_end_of_timestep(i) = currenttime + 10.0*25000.0
            else
             ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) + 10.0*25000.0 
            endif
            RibaQdemand = 0.0
        else if (i.le.4) then !new Modflow stress period from i = 20         
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) + 10.0*25000.0
            RibaQdemand(1) = -.2000E+05
            RibaQdemand(2) = -.0000
            RibaQdemand(3) = -.6685E+05
            RibaQdemand(4) = -.6685E+05
            RibaQdemand(5) = -.0000E+05
            RibaQdemand(6) = -.1000E+06
            RibaQdemand(7) = -.6685E+05
            RibaQdemand(8) = -.6685E+05
            RibaQdemand(9) = -.6685E+05
            RibaQdemand(10) = -.6685E+05
            RibaQdemand(11) = -.6685E+05
            RibaQdemand(12) = -.6685E+05
            RibaQdemand(13) = -.6685E+05
            RibaQdemand(14) = -.6685E+05
            RibaQdemand(15) = -.1003E+06
            RibaQdemand(16) = -.6685E+05
            RibaQdemand(17) = -.1003E+06
        elseif (i.le.5) then !new Modflow stress period from i = 40
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) + 20.0*3.0
            RibaQdemand(1) = -.2000E+05
            RibaQdemand(2) = -.0000
            RibaQdemand(3) = -.6685E+05
            RibaQdemand(4) = -.6685E+05
            RibaQdemand(5) = -.0000E+05
            RibaQdemand(6) = -.1300E+06
            RibaQdemand(7) = -.6685E+05
            RibaQdemand(8) = -.6685E+05
            RibaQdemand(9) = -.6685E+05
            RibaQdemand(10) = -.6685E+05
            RibaQdemand(11) = -.6685E+05
            RibaQdemand(12) = -.6685E+05
            RibaQdemand(13) = -.6685E+05
            RibaQdemand(14) = -.6685E+05
            RibaQdemand(15) = -.1003E+06
            RibaQdemand(16) = -.6685E+05
            RibaQdemand(17) = -.1003E+06
        elseif (i.le.7) then !new Modflow stress period from i = 60
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) + 10.0*9.0
        elseif (i.le.12) then !new Modflow stress period from i = 80
            ribasim_end_of_timestep(i) = ribasim_end_of_timestep(i-1) + 10.0*15.0
        endif
        if (i.gt.5) then ! Qdemand values remain unchanged after Modflow stress period 3
            RibaQdemand(1) = -.2000E+05
            RibaQdemand(2) = -.000
            RibaQdemand(3) = -.6685E+05
            RibaQdemand(4) = -.6685E+05
            RibaQdemand(5) = -.0000E+05
            RibaQdemand(6) = -.1300E+06
            RibaQdemand(7) = -.6685E+05
            RibaQdemand(8) = -.6685E+05
            RibaQdemand(9) = -.6685E+05
            RibaQdemand(10) = -.6685E+05
            RibaQdemand(11) = -.6685E+05
            RibaQdemand(12) = -.6685E+05
            RibaQdemand(13) = -.6685E+05
            RibaQdemand(14) = -.6685E+05
            RibaQdemand(15) = -.1003E+06
            RibaQdemand(16) = -.6685E+05
            RibaQdemand(17) = -.1003E+06
        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_1