subroutine modflow_dll_init(config_file, start_time)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_init
! modules
use modflow_dll_module
use m_main_info
use m_vcl, only: targ
use imod_utl, only: imod_utl_has_ext, imod_utl_printtext
use mod_pest, only: pest1_meteo_metaswap, pest1alpha_metaswap, pest1appendlogfile, pestnext

implicit none

! arguments
character(len=*), intent(in) :: config_file
double precision, intent(in) :: start_time

! general
 logical :: mf_steadystate

! PEST with MetaSWAP
 integer :: mf_nrow, mf_ncol, mf_nlay, kper
 
! parameters
 integer, parameter :: nsubmax = 1000

! functions
 integer   osd_open2,cfn_length
 integer   cfn_idx_get_i
 double precision cfn_mjd_nodata

! local variables

 logical :: ok
 character  del*1,cont*1,comm*1

 integer, allocatable, dimension(:) :: hrivsys, wrivsys
 integer :: nhrivsys, nwrivsys
 integer, dimension(1) :: m
 character (len=1)    :: cdum
 character (len=1024) :: record, modrecord, sobrecord, mozrecord
 character (len=32)   :: compcls
 character (len=1024) :: compfile

 double precision :: dtMetaSwap            ! variable for metaswap
 double precision :: dtMozart
 integer          :: iterMetaSwap          ! variable for metaswap
 double precision :: stsTestTime

 integer                         :: stsFtest,stsBtest,stsItest,i,j,n
 double precision, allocatable   :: stsTtest(:)

 character tekens(1)*1

 logical :: lrunfile, lnamfile, llpf, lipest, lpwt, lrf

 integer :: isub, nsub, nnsub
 character(len=50), dimension(nsubmax) :: submstr

! debug
 integer lswid, js, je, k, ilay, irow, icol, lun, cfn_getlun

! program section
! ------------------------------------------------------------------------------

! ... init
 lunc        = 9            ! unit number for components file
 lunsts      = 8            ! <=0: sts2init() will assign a free unit number
                            !  >0: unit number will be used
 comm        = '!'          ! comment sign for components input file
 cont        = char(92)     ! continuation mark for components input file records
 nodataTime  = cfn_mjd_nodata()

 nsub = 0

! ... process initialisation
! ... init components
!     ios: 0: OK
!         -1: commandline argument not found
!        <-1: ERROR
 ! define virtual commandline
 call cfn_vcl_set(' ',ivcl)

 ! get the root name
 call osd_getcwd(root)
 modwd1  = root
 simwd1  = root
 mozwd   = root

 ! run-file
 lrunfile = .false.
 lnamfile = .false.
 usemodflow    = .false.
 usemetaswap   = .false.
 lipest = .false.
 infile = config_file
 if (imod_utl_has_ext(infile,'nam')) lnamfile= .true.
 if (imod_utl_has_ext(infile,'run')) lrunfile= .true.
 if (lrunfile .or. lnamfile) then ! run-file or nam-file
    usemodflow=.true.
    modrecord=''
    if (lnamfile) write(modrecord,'(a,1x,a)') '-namfile',trim(infile)
    wd = ''
    call rf2mf_prg(lrunfile,lipest,modrecord,usemetaswap,submstr,nsub,nsubmax,wd)
    modwd1 = trim(wd)
    call osd_s_filename(modwd1)
 end if

 usests        = .true.
 usestsmodflow = .true.
 usetransol    = .false.
 usemozart     = .false.
 useribasim    = .true.
 ! get explicit arguments

 call cfn_vcl_fndc(ivcl,iarg,'-components',.true.,compfile,1)
 
 rt = driverGetRunType(usemodflow,usemetaswap,usetransol,usemozart,usests,useribasim)
 write(*,'(50(''*''),/,2(1x,a),/,50(''*''))') 'Running mode:', trim(rtdesc(rt))

! check state-save options
 if (usests.and..not.usestsmodflow) call driverChk(.false.,'MODFLOW state-save not activated.')
 if (.not.usests.and.usestsmodflow) call driverChk(.false.,'MODFLOW state-save activated.')

 nnsub = max(1,nsub)
 
 modwd2 = modwd1
 simwd2 = simwd1
 

 timestep    = 1
 stsave      = .false.
 strestore   = .false.

 ! init return values
 retValMF2005      = 0

 ! time values
 currentTime       = nodataTime
 currentTimeMF2005 = 0.d0

 ! init all components TODO!!!!!!!
 if (usemodflow) then
    call osd_chdir(modwd2)
    if (usestsmodflow) call sts2init(usestsmodflow,lunsts)
    call mf2005_initComponent(modrecord,retValMF2005)
    ! get number of grids
    ok = mf2005_PutNumberOfGrids(mf_ngrid); call driverChk(ok,'mf2005_PutNumberOfGrids')
    if (mf_ngrid.ne.1) call driverChk(.false.,'More than one MODFLOW grids is not supported.'); mf_igrid = 1    
 end if

 ! check if MODFLOW is activated
 if (retValMF2005.ne.0) call driverChk(.false.,'MODFLOW initialization failed.')

 if (start_time > -1) then
    ! outside world provides starting date/time
    call mf2005_setCurrentTime(start_time,retValMF2005)
    currentTime = start_time
 else
    ! get starting date/time from modflow
    call mf2005_getCurrentTime(currentTimeMF2005,retValMF2005)
    currentTime = currentTimeMF2005
 endif

! start sts
 if (usestsmodflow) call sts2start(currentTime)

! read data for entire simulation
 call osd_chdir(modwd2)
 call mf2005_initSimulation(currentTime,retValMF2005) !Allocate and read MF data constant for entire simulation
 if (retValMF2005.ne.0) exitcode = -12
! #### BEGIN EXCHANGE: BeforeInitSimulationMozart #############################

 ok = mf2005_PutLPFActive(mf_igrid, llpf); call driverChk(ok,'mf2005_PutLPFActive') ! get flag if BCF is used or LPF
 ok = mf2005_PutSimulationType(mf_igrid, lss); call driverChk(ok,'mf2005_PutSimulationType')  ! get flag if simulation is steady-state or transient
 ok = mf2005_PutHeadNoFlo(mf_igrid, hnoflo);  call driverChk(ok,'mf2005_PutHeadNoFlo')  ! get hnoflo

 !#### TIMESERIES ####
 call tserie1init1(lipest,lss,hnoflo)
 ok = mf2005_TimeserieInit(mf_igrid)
 call driverChk(ok,'mf2005_TimeserieInit')
 call tserie1init2(lipest,lss,hnoflo)

! #### BEGIN EXCHANGE: BeforeTimestep #########################################
 ! Coupling MODFLOW-Ribasim
 if (rt.eq.rtmodriba) then
    ok = mf2005_PutModRibaNumberOfIDs(mf_igrid, XchModRibaModNID); call driverChk(ok,'mf2005_PutModRibaNumberOfIDs') ! MODFLOW puts the number of MOD-RIBA IDs
!    ok = ribasim_PutModRibaNumberOfIDs(XchModRibaRibaNID); call driverChk(ok,'ribasim_PutModRibaNumberOfIDs')       ! Ribasim puts the number of MOD-RIBA IDs
     XchModRibaRibaNID = 17 !temp
    call driverXchInitModRiba() ! allocate and init

    ok = mf2005_PutModRibaIDs(mf_igrid,XchModRibaModIds); call driverChk(ok,'mf2005_PutModRibaIDs')       ! MODFLOW puts the MOD-Riba exchange IDs
    ok = mf2005_PutModRibaCells(mf_igrid,XchModRibaModCells); call driverChk(ok,'mf2005_PutModRibaCells') ! MODFLOW puts the MOD-Riba cells (per ID the lay,row,col)
!    ok = ribasim_PutModRibaIDs(XchModRibaRibaIds); call driverChk(ok,'ribasim_PutModRibaIDs')             ! get exchange id's   
    do i=1,XchModRibaRibaNID; XchModRibaRibaIds(i)=i; enddo !temp
 end if

  if (rt.eq.rtmodriba) call driverXchIniMapModRiba()
! ###### END EXCHANGE: BeforeTimestep #########################################

! ... timestep
 timestep = 0

 tsc=0 ! debug: TimeStepCycle

 end subroutine modflow_dll_init

subroutine modflow_dll_get_grid_sizes(mmax, nmax, kmax)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_get_grid_sizes
  integer, intent(out) :: mmax, nmax, kmax
  mmax = -1
  nmax = -1
  kmax = -1
  ! TODO
end subroutine modflow_dll_get_grid_sizes
    
subroutine modflow_dll_set_recharge(recharge, mmax, nmax)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_set_recharge
  integer, intent(in)                        :: mmax, nmax
  real   , intent(in), dimension(mmax, nmax) :: recharge
  ! TODO
end subroutine modflow_dll_set_recharge
    
subroutine modflow_dll_get_heads(head, mmax, nmax)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_get_heads
  integer, intent(in)                         :: mmax, nmax
  real   , intent(out), dimension(mmax, nmax) :: head
  head = -888.0
  ! TODO
end subroutine modflow_dll_get_heads
    
 function modflow_dll_timestep() result(retVal)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_timestep
 use modflow_dll_module
 use m_main_info

 implicit none

 integer :: retVal
 integer :: date, hour, minute, second

 logical :: converged,convergedMF2005,convergedMetaSwap,convergedMozart,convergedPest
 logical :: ok ! doTimeLoop
 integer :: kper, kstp

 
 ! convergence
 convergedMF2005   = .true.
 convergedMetaSwap = .true.
 convergedMozart   = .true.

 retVal  = 0            ! when 0: ok
 stsave    = .true.
 ribastsave = .true.
 strestore = .false.

 tsc=tsc+1
 timestep = timestep + 1

 call cfn_mjd2datehms(currentTime,date,hour,minute,second)

 write(*,'(5x,a,1x,i5,1x,a,1x,i8,3(a,i2.2))')&
 'Timestep     ',tsc,':',date,' ',abs(hour),':',minute,':',second

 ! one timestep for each cycle

!##### BEGIN EXCHANGE: BeginOfTimestep ########################################

!####### END EXCHANGE: BeginOfTimestep ########################################

 ! perform state save, phase 1 (before reading and writing data)
 if (usests) then
    if (ribastsave) stsave = .true. !
    call osd_chdir(modwd2)
    if (usestsmodflow) call sts2saverestore(currentTime,stsave,strestore,1) !external decision whether stsave is needed, internal decision whether strestore is needed.
 end if

 ! ... init timestep
 call mf2005_prepareTimestep(currentTime,stsave,retValMF2005,kper,kstp) !checks o.a. whether new stress period and/or new time step
 call osd_chdir(modwd2)
          
 call mf2005_initTimeStep(currentTime,retValMF2005,rt,rtmodriba) !Read and prepare for new stress period and new timestep
 if (retValMF2005.ne.0) retVal = -15
          
 ! get end of current timestep
 !!!!!!!!!!!!!!!! Only modflow information is used !!!!!!!!!!!!!!!!!!!!!!!
 call mf2005_getEndOfCurrentTimeStep(endOfCurrentTimeStepMF2005,retValMF2005)
 if (retValMF2005.ne.0) retVal = -17
 endOfCurrentTimeStep = endOfCurrentTimeStepMF2005

 ! extra check
 if (endOfCurrentTimeStep.eq.nodataTime) call driverChk(.false.,'endOfCurrentTimeStep = nodataTime')

! ... iteration
 converged =.false.
 do while (.not. converged .and. retVal.eq.0)
     call mf2005_performIter(retValMF2005)
     if (retValMF2005.ne.0) retVal = -26

!##### BEGIN EXCHANGE: AfterSolve #############################################
             !(No exchange with Ribasim within iterations)
!####### END EXCHANGE: AfterSolve #############################################

     ! finish iteration
     call mf2005_finishIter(convergedMF2005,retValMF2005)

     ! get next iteration information
     converged = .true.
     converged = converged .and. convergedMF2005

     if (retValMF2005.ne.0    ) retVal = -261
 enddo
 if (.not.converged) write(*,*) ' Model did not converge!',retVal

! ... write results
 if (converged .and. retVal.eq.0) then
     call osd_chdir(modwd2)
     call mf2005_finishTimestep(retValMF2005)
     if (retValMF2005.ne.0) retVal = -31

!##### BEGIN EXCHANGE: AfterFinishTimeStepMODFLOW #############################

!##### BEGIN EXCHANGE: AfterFinishTimeStepMODFLOW #############################

     !#### TIMESERIES #####
     ok = mf2005_TimeserieGetHead(mf_igrid); call DriverChk(ok,'mf2005_TimeserieGetHead')
    call tserie1write(0,lss,currentTime,hnoflo,usests)
     !#### TIMESERIES #####

 endif

 if (retVal.eq.0) then
     ! perform state save, phase 2 (after reading and writing data)
     if (usestsmodflow) call sts2saverestore(currentTime,stsave,strestore,2)
 endif

! ... next timestep
 if (retVal.ne.0) then
     ! ERROR, do not continue
     endOfSimulation=.true.
 else
     ! get new time information
     ! Let Modflow be leading. don't ask the other components
     call mf2005_getBeginOfNextTimeStep(BeginOfNextTimeStepMF2005,retValMF2005)
     ! find errors
     if (retValMF2005.eq.-1) then
         ! end of simulation for Modflow
         endOfSimulation=.true.
         retValMF2005   =0
     endif
     if (retValMF2005.ne.0) then
         retVal = -31
     endif
     if (retVal.ne.0) then
         endOfSimulation=.true.
     endif
     ! find start time of next timestep
     !   this may be a time in the 'future', the 'past' or the 'current' time
     ! nodata value of time-value means: this was the last timestep for this component
     
     if (retVal.eq.0) then
         currentTime=nodataTime

         ! check for each component (ok, only modflow available now)
         if (BeginOfNextTimeStepMF2005.ne.nodataTime) then
             if (currentTime.eq.nodataTime) then
                 currentTime=BeginOfNextTimeStepMF2005
             else
                 currentTime=min(currentTime,BeginOfNextTimeStepMF2005)
             endif
         endif

         if (currentTime.eq.nodataTime.and.rt.ne.rtmod) then
             endOfSimulation=.true.
         endif
!        if (currentTime.eq.nodataTime) endOfSimulation=.true.
     endif !retVal.eq.0
 endif !retVal.ne.0

!   doTimeLoop = .not.endOfSimulation
!   enddo ! time loop

 end function modflow_dll_timestep
 
 subroutine modflow_dll_set_timestep_done(done_with_current_timestep, Reset_time)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_set_timestep_done
    use modflow_dll_module
    implicit none
    integer nribasimtimestep
    double precision Reset_Time,LastSavedTime
    logical done_with_current_timestep
    if (.not. done_with_current_timestep) then
        ! jump back to begin of time step
        currenttime = Reset_time
        call set_last_saved_time(Reset_time)
    else
        ! time step done
    endif
    
 end subroutine modflow_dll_set_timestep_done
 
 subroutine set_last_saved_time(LastSavedTime)
 !DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_set_timestep_done
     use m_sts2
     implicit none
     double precision LastSavedTime
 
     currenttime = LastSavedTime
 !    Last_Saved_Time = currenttime
 end subroutine
 
 function mf_end_of_simulation_reached() result(end_of_simulation_reached)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_end_of_simulation_reached
     use modflow_dll_module
     implicit none

     logical end_of_simulation_reached
     end_of_simulation_reached = endOfSimulation
 end function mf_end_of_simulation_reached
 
 function mf_external_timestep_reached(external_end_of_timestep) result(external_timestep_reached)
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_external_timestep_reached
     use modflow_dll_module
     implicit none
 
     logical external_timestep_reached
     double precision external_end_of_timestep
 
     if(currenttime.ge.external_end_of_timestep)then
         external_timestep_reached = .true.
     else
         external_timestep_reached = .false.
     endif
 end function
     
 
 subroutine modflow_dll_finish
 !DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_finish
 use modflow_dll_module
 use imod_utl, only: imod_utl_closeunits

 implicit none

 ! ... end
 call osd_chdir(modwd2)
 if (exitcode.eq.0)  call mf2005_finishSimulation(retValMF2005)
 if (retValMF2005.ne.0 ) exitcode = -61

!#### TIMESERIES #####
 call tserie1write(1,lss,currentTime,hnoflo,usests)
 call tserie1close()
!#### TIMESERIES #####

! list routine timer info
 call cfn_rtt_list(0)

! Deallocate coupling arrays
 call driverXchDeallocate()

! exit
 if (exitcode.ne.0) then
    write(unit=*,fmt=*) 'ERROR, exit code: ',exitcode
    call exit(10)
 endif

 call imod_utl_closeunits()

end subroutine modflow_dll_finish

function modflow_dll_get_number_of_wells() result(number_of_wells)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_get_number_of_wells
use modflow_dll_module
    implicit none
    integer :: number_of_wells
    number_of_wells = XchModRibaModNID
end function modflow_dll_get_number_of_wells

subroutine get_currenttime(current_time)
!DEC$ ATTRIBUTES DLLEXPORT :: get_currenttime
use modflow_dll_module
    implicit none
    double precision current_time
    current_time = currenttime


end subroutine    


subroutine set_well_values(RibaQdemand)
!DEC$ ATTRIBUTES DLLEXPORT :: set_well_values
use modflow_dll_module
    implicit none
    real, dimension(17), intent(in) :: RibaQdemand
logical :: ok

ok = mf2005_GetQdemand(mf_igrid,XchModRibaModNID,RibaQdemand,XchRiba2ModIdx,XchRiba2ModOff,mv); call DriverChk(ok,'mf2005_GetQdemand') ! Modflow gets Qdemand 

end subroutine set_well_values


subroutine get_well_values(RibaQrealized)
!DEC$ ATTRIBUTES DLLEXPORT :: get_well_values
use modflow_dll_module
real, dimension(17), intent(out) :: RibaQrealized
logical :: ok

ok = mf2005_PutQrealized(mf_igrid,RibaQrealized,XchModRibaModNID,XchMod2RibaIdx,XchMod2RibaOff,mv); call DriverChk(ok,'mf2005_PutQrealized') ! MODFLOW puts the realized extraction rates flux for Ribasim
end subroutine get_well_values