program modflow_dll_test_3

    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

    double precision  :: currenttime, dummy_start_time
    double precision, dimension(:,:), allocatable :: head1
    real, dimension(:,:), allocatable :: recharge
    
    logical :: convergedPest, endOfSimulation, convergedMozart, doTimeLoop, convergedIter, enableStateSave
    integer           :: retVal = 0
    integer           :: nlay,nrow,ncol
    integer           :: isub, nnsub

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

    enableStateSave = .true.
    call mf_init(config_file_name,nnsub,enableStateSave)
    do isub = 1, nnsub ! subdomain loop
       convergedPest = .False.
       do while (.not.convergedPest) ! pest loop
          endOfSimulation   = .false.
          dummy_start_time  = -1.0d0
          call mf_init_simulation(endOfSimulation,convergedMozart, dummy_start_time)

          !!!!!! get grid dimension and allocate head/recharge !!!!!
          call mf_get_grid_sizes(nlay,nrow,ncol)
          if (allocated(head1)) deallocate(head1)
          if (allocated(recharge)) deallocate(recharge)
          allocate(head1(ncol,nrow),recharge(ncol,nrow))
          recharge = 0.
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          
          do while (.not.endOfSimulation) ! simulation loop
             call mf_init_mozart(); convergedMozart = .false.
             do while (.not.endOfSimulation .and. .not.convergedMozart) ! mozart loop
                call mf_init_timestep()
                doTimeLoop = .true.
                do while(doTimeLoop) ! time loop 
                    call mf_init_iter(enableStateSave)
                    
                    !!!!!! set recharge !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    recharge = recharge + 0.001
                    call mf_set_recharge(recharge,ncol,nrow)  
                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    
                    convergedIter = .false.
                    do while(.not.convergedIter) ! Picard iter loop
                       call  mf_iter(convergedIter)    
                    end do
                    call mf_finish_timestep(convergedIter,endOfSimulation,doTimeLoop)
                    
                    !!!!! get the head for layer 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!
                    call mf_get_heads(head1,ncol,nrow)
                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                
                end do   
             end do
             call mf_finish_simulation(endOfSimulation,convergedMozart)
          end do
          call mf_finish_pest(convergedPest)
       end do
    end do

    end program modflow_dll_test_3