subroutine mf_init(config_file,nnsub,stsOverRule)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_init
 use mf_module
 use m_main_info
 use m_vcl, only: targ
 use imod_utl, only: imod_utl_closeunits, 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 ! nam-file, components-file or run file
 integer, intent(out) :: nnsub
 logical, intent(in) :: stsOverRule
    
! functions
 integer   osd_open2,cfn_length
 integer   cfn_idx_get_i
 double precision cfn_mjd_nodata

! local variables

 character  del*1,cont*1,comm*1

 integer, dimension(1) :: m
 character (len=1)    :: cdum
 character (len=32)   :: compcls
 character (len=1024) :: compfile

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

 character tekens(1)*1

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

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

 call imod_utl_printtext('====================================================================',0)
 call imod_utl_printtext('You may use this compiled version of the iMOD-software if you are '  ,0)
 call imod_utl_printtext('entitled to this use under a iMOD software license agreement for the',0)
 call imod_utl_printtext('iMOD software executables with Deltares or with a party entitled by' ,0)
 call imod_utl_printtext('Deltares to provide sublicenses for the iMOD-software executables.'  ,0)
 call imod_utl_printtext('Otherwise use of this compiled version of the iMOD-software is'      ,0)
 call imod_utl_printtext('prohibited and illegal. If you are not allowed under a Deltares iMOD',0)
 call imod_utl_printtext('license agreement to use the iMOD-software executables, you may find',0)
 call imod_utl_printtext('a solution in compiling the open source version of the iMOD-software',0)
 call imod_utl_printtext('into an executable yourself (see oss.deltares.nl), or apply for a'   ,0)
 call imod_utl_printtext('Deltares iMOD license agreement by sending an email to'              ,0)
 call imod_utl_printtext(''                                                                    ,0)
 call imod_utl_printtext('Version 3.00.01, 18/09/14'                                           ,0)
 call imod_utl_printtext('====================================================================',0)

! ... init
 exitcode    = 0            ! when 0: ok
 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.
! call cfn_vcl_narg(ivcl,narg)
! if (narg.eq.0) then
!  call imod_utl_printtext(' ',0)
!  call imod_utl_printtext('Driver program valid arguments:',0)
!  call imod_utl_printtext(' 1: <iMOD run-file> <optional subdomain number>',0)
!  call imod_utl_printtext(' 2: <MODFLOW-2005 nam-file>',0)
!  call imod_utl_printtext(' 3: -components <components steering file>',0)
!  call exit(0)
! end if
! call cfn_vcl_arg(ivcl,1,infile,n)
 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)
    if (lrunfile) then
       write(modrecord,'(a,1x,a)') '-runfile',trim(infile)
       if (narg.eq.3) then ! runfile options
          call cfn_vcl_arg(ivcl,2,rfopt,n)
          write(modrecord,'(a,1x,a,1x,a)') trim(modrecord),'-rfopt',trim(rfopt)
          call cfn_vcl_arg(ivcl,3,rfopt,n)
          write(modrecord,'(a,1x,a)') trim(modrecord),trim(rfopt)
       end if
    end if
    wd = ''
    call rf2mf_prg(lrunfile,lipest,modrecord,usemetaswap,submstr,nsub,nsubmax,wd)
    modwd1 = trim(wd)
    call osd_s_filename(modwd1)
    if (usemetaswap) then
      simwd1 = trim(wd)
    end if
 end if

 usests        = .false.
 usestsmodflow = .false.
 usetransol    = .false.
 usemozart     = .false.
 useribasim    = .false.
 ! get explicit arguments
 call cfn_vcl_fndc(ivcl,iarg,'-components',.true.,compfile,1)
 if (iarg.gt.0) then

    ! open file and initialise components
    ios=osd_open2(lunc,0,compfile,'readonly')
    if (ios.eq.0) then

       ! read all component information and initialise all components
       call cfn_getrec(lunc,record,comm,cont)
       do while (cfn_length(record).gt.0)
          ! extract the first argument of variable record,
          ! leave the rest of the record for initialisation of the component
          tekens(1) = ' '
          call cfn_par_ext(compcls,record,tekens,1,del)
          call cfn_token(compcls,'tu')
          lc=cfn_length(compcls)

          select case( compcls(1:lc) )

          case( 'MODFLOW' )
                modrecord = record
                usemodflow=.true.
                ! set virtual command line
                call cfn_vcl_set(modrecord,jvcl)
                ! set working directory
                call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
                if (jarg.gt.0) modwd1 = trim(modwd1)//wd
                call osd_chdir(modwd1)
                ! check if state save is enabled
                call cfn_vcl_fnd(jvcl,jarg,'-sts',.true.)
                if (jarg.gt.0) then
                   usestsmodflow = .true.
!                   call sts2init(usestsmodflow,lunsts)
                end if
                ! convert iMOD run-file to MODFLOW
                call rf2mf_prg(lrf,lipest,modrecord,usemetaswap,submstr,nsub,nsubmax,modwd1)
                if (lrf) then
                   modwd1 = trim(modwd1)//'\mf2005_tmp'
                   call osd_s_filename(modwd1)
                end if
#ifdef INCLUDE_METASWAP
             case( 'METASWAP' )
                usemetaswap=.true.
                ! set working directory
                call cfn_vcl_set(record,jvcl)
                call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
                if (jarg.gt.0) simwd1 = trim(simwd1)//wd
             case( 'TRANSOL' )
                usetransol=.true.
#endif
#ifdef INCLUDE_MOZART
             case( 'MOZART' )
                mozrecord=record
                usemozart=.true.
                ! set working directory
                call cfn_vcl_set(mozrecord,jvcl)
                call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
                if (jarg.gt.0) mozwd = trim(mozwd)//wd
                ! get the river subsystems for coupling
                call cfn_vcl_fndi(jvcl,jarg,'-nhriv*sys',.true.,m,1)
                if (m(1).gt.0) then
                   nhrivsys = m(1)
                   allocate(hrivsys(nhrivsys))
                else
                   call driverChk(.false.,'Error: -nhriv*sys not found.')
                end if
                call cfn_vcl_fndi(jvcl,jarg,'-hriv*sys',.true.,hrivsys,nhrivsys)
                call cfn_vcl_fndi(jvcl,jarg,'-nwriv*sys',.true.,m,1)
                if (m(1).gt.0) then
                   nwrivsys = m(1)
                   allocate(wrivsys(nwrivsys))
                else
                   call driverChk(.false.,'-nwriv*sys not found.')
                end if
                call cfn_vcl_fndi(jvcl,jarg,'-wriv*sys',.true.,wrivsys,nwrivsys)
             case( 'STS' )
                usests = .true.
#endif
             case default
                ! ERROR, component unknown
                write(*,*) ' Warning: component ',compcls(1:cfn_length(compcls)),' unknown.'

          end select

          ! get next component
          call cfn_getrec(lunc,record,comm,cont)
       enddo

       ! close components file
       close(lunc)
    else
       call driverChk(.false.,'Could not open component description file '//compfile(1:cfn_length(compfile)))
       call exit(1)
    endif
 endif

 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)   
 
 if (stsOverRule) then
    usests = .true.
    usestsmodflow = .true.
 end if    
    
end subroutine mf_init

subroutine mf_get_current_time(timeValue,retVal)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_get_current_time
    ! arguments
    double precision, intent(out)   :: timeValue      !> start time of current time step
    integer         , intent(out)   :: retVal         !> return code   0: OK
    call mf2005_getCurrentTime(timeValue,retVal)
end subroutine mf_get_current_time

subroutine mf_init_iter(stssaveOverRule)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_init_iter
 use mf_module

 implicit none
 
 logical, intent(in) :: stssaveOverRule
 
 integer :: idum1, idum2 !!!!!!!!!!!! WORKAROUND !!!!!!!!!!!!!!
    
 tsc=tsc+1
 timestep = timestep + 1

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

!what about steady-state solutions, currenttime.eq.0
 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 ########################################
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
    ok = mf2005_PutHeads(mf_igrid,XchModSimModCells,XchModSimModNID,XchModSimModHeads,mv); call DriverChk(ok,'mf2005_PutHeads') ! Put groundwater heads
    ok = metaswap_GetHeads(XchModSimModHeads,XchModSimModNID,XchMod2SimIdx,XchMod2SimOff,mv); call DriverChk(ok,'metaswap_GetHeads') ! Get groundwater heads
    XchModSimModHeads = mv
 end if
!####### END EXCHANGE: BeginOfTimestep ########################################

 ! perform state save, phase 1 (before reading and writing data)
 if (usests) then
    if (stssaveOverRule) stsave= .true.
    if (mozstsave) stsave = .true.
    call osd_chdir(modwd2)
    if (usestsmodflow) call sts2saverestore(currentTime,stsave,strestore,1)
    if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
       call osd_chdir(simwd2)
       call MetaSWAP_saveRestore(stsave,strestore)
       !call TRANSOL_saveRestore(stsave,strestore)
    end if
 end if

!  ... init timestep
 call mf2005_prepareTimestep(currentTime,stsave,retValMF2005,idum1,idum2)
 call osd_chdir(modwd2)
 call mf2005_initTimeStep(currentTime,retValMF2005)
 if (retValMF2005.ne.0) exitcode = -15

 ! get end of current timestep
 !!!!!!!!!!!!!!!! Only modflow information is used !!!!!!!!!!!!!!!!!!!!!!!
 call mf2005_getEndOfCurrentTimeStep(endOfCurrentTimeStepMF2005,retValMF2005)
 if (retValMF2005.ne.0) exitcode = -17
 endOfCurrentTimeStep = endOfCurrentTimeStepMF2005

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

 ! read timestep data
 if (exitcode.eq.0) then
    if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
       dtMetaSwap = endOfCurrentTimeStep - currentTime
       call osd_chdir(simwd2)
       call MetaSwap_initTimestep(dtMetaSwap)
       iterMetaSwap = 0
    end if
 end if

end subroutine mf_init_iter         
              
subroutine mf_iter(converged)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_iter
 use mf_module             
           
 implicit none
 
 logical, intent(inout) :: converged     
    
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
    call MetaSwap_prepareIter(iterMetaSwap)  ! only action: iter=iter+1
    call MetaSwap_performIter(iterMetaSwap)
!##### BEGIN EXCHANGE: AfterSolveMetaSwap #####################################
    ok = metaswap_PutModSimUnsaturatedZoneFlux(XchModSimSimUnsZFlux,mv); call DriverChk(ok,'metaswap_PutUnsaturatedZoneFlux') ! Put unsaturated zone flux
    ok = mf2005_GetUnsaturatedZoneFlux(mf_igrid,&
            XchModSimModNID,XchModSimSimUnsZFlux,XchSim2ModIdx,XchSim2ModOff,mv); call DriverChk(ok,'mf2005_GetUnsaturatedZoneFlux') ! Get unsaturated zone flux
    XchModSimSimUnsZFlux = mv
    if (.not.mf_steadystate) then
       ok = metaswap_PutStorageFactor(XchModSimSimStrFct,mv); call DriverChk(ok,'metaswap_PutStorageFactor') ! Put storage factor
       if (.not.llpf) then ! BCF
          ok = mf2005_GetStorageFactor(mf_igrid,XchModSimSimStrFct,&
                                      XchModSimModNID,XchSim2ModIdx,XchSim2ModOff,mv); call DriverChk(ok,'mf2005_GetStorageFactor') ! Get storage factor
       else ! LPF
          ok = mf2005_GetStorageFactorLPF(mf_igrid,XchModSimSimStrFct,&
                                      XchModSimModNID,XchSim2ModIdx,XchSim2ModOff,mv); call DriverChk(ok,'mf2005_GetStorageFactorLPF') ! Get storage factor
       end if
       XchModSimSimStrFct = mv
    end if
 end if
!####### END EXCHANGE: AfterSolveMetaSwap #####################################
 call mf2005_performIter(retValMF2005)
 if (retValMF2005.ne.0) exitcode = -26

!##### BEGIN EXCHANGE: AfterSolve #############################################
 if (exitcode.eq.0) then
    if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
       ok = mf2005_PutHeads(mf_igrid,XchModSimModCells,XchModSimModNID,XchModSimModHeads,mv); call DriverChk(ok,'mf2005_PutHeads') ! MODFLOW put groundwater heads
       ok = metaswap_GetHeads(XchModSimModHeads,XchModSimModNID,XchMod2SimIdx,XchMod2SimOff,mv); call DriverChk(ok,'metaswap_GetHeads') ! MetaSWAP gets the groundwater heads
       XchModSimModHeads = mv
    end if
 end if
!####### END EXCHANGE: AfterSolve #############################################

 ! finish iteration
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz)&
 call MetaSwap_finishIter(iterMetaSwap,convergedMetaSwap)
 call mf2005_finishIter(convergedMF2005,retValMF2005)
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz)&
    write(*,*) ' convergence MODFLOW - MetaSWAP: ',convergedMF2005, convergedMetaSwap

 ! get next iteration information
 converged = .true.
 converged = converged .and. convergedMF2005
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz)&
    converged = converged .and. convergedMetaSwap

 if (retValMF2005.ne.0    ) exitcode = -261
 
 if (exitcode.ne.0) converged = .true.
            
end subroutine mf_iter         
          
subroutine mf_finish_timestep(converged,endOfSimulation,doTimeLoop)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_finish_timestep
 use mf_module                   
       
 implicit none
 
 logical, intent(in) :: converged
 logical, intent(inout) :: endOfSimulation,doTimeLoop 
 
 if (.not.converged) write(*,*) ' Model did not converge!',exitcode

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

!##### BEGIN EXCHANGE: AfterFinishTimeStepMODFLOW #############################
    if (rt.eq.rtmodsimtranmoz) then

       ok = mf2005_PutSeepageFlux(mf_igrid,XchModTranModSeepageFlux,XchModTranModCells,XchModTranModNID,mv,.true.); call DriverChk(ok,'mf2005_PutSeepageFlux') ! MODFLOW puts the seepage flux for TRANSOL
       ok = mf2005_PutRiverFlux(mf_igrid,XchModTranModRiverFlux,XchModTranModCells,XchModTranModNID,mv,&
                                nhrivsys,hrivsys,nwrivsys,wrivsys,.true.,.false.); call DriverChk(ok,'mf2005_PutRiverFlux') ! MODFLOW puts the river flux for TRANSOL: skip H and W [m]
       ok = mf2005_PutRiverFlux(mf_igrid,XchModMozModRiverFlux,XchModMozModCells,XchModMozModNID,mv,&
                                nhrivsys,hrivsys,nwrivsys,wrivsys,.false.,.false.); call DriverChk(ok,'mf2005_PutRiverFlux') ! MODFLOW puts the river flux for MOZART: skip H and W [m3]
       ok = mf2005_PutRiverFlux(mf_igrid,XchModMozModRiverFluxWells,XchModMozModCells,XchModMozModNID,mv,&
                                nhrivsys,hrivsys,nwrivsys,wrivsys,.false.,.true.); call DriverChk(ok,'mf2005_PutRiverFlux') ! MODFLOW puts the river flux for MOZART ("wellen"): skip H; only W; [m3]
       ok = mf2005_PutDrainFlux(mf_igrid,XchModTranModDrainFlux,XchModTranModCells,XchModTranModNID,mv,.true.); call DriverChk(ok,'mf2005_PutDrainFlux') ! MODFLOW puts the drain flux for TRANSOL: don't skip; all layers; [m]
       ok = mf2005_PutDrainFlux(mf_igrid,XchModMozModDrainFlux,XchModMozModCells,XchModMozModNID,mv,.false.); call DriverChk(ok,'mf2005_PutDrainFlux') ! MODFLOW puts the drain flux for MOZART; don't skip; all layers; [m3]
       ok = mf2005_PutSaltFlux(mf_igrid,XchModMozModSalt,XchModMozModCells,XchModMozModNID,mv,nwrivsys,wrivsys); call DriverChk(ok,'mf2005_PutSaltFlux') ! MODFLOW puts the salt flux for MOZART
       ok = TRANSOL_GetSeepageRiverDrainFlux(XchModTranModSeepageFlux,XchModTranModNID,XchMod2TranIdx,XchMod2TranOff,mv,'see'); call DriverChk(ok,'TRANSOL_GetSeepageRiverDrainFlux') ! TRANSOL gets the seepage flux
       XchModTranModSeepageFlux = mv
       ok = TRANSOL_GetSeepageRiverDrainFlux(XchModTranModRiverFlux,XchModTranModNID,XchMod2TranIdx,XchMod2TranOff,mv,'riv'); call DriverChk(ok,'TRANSOL_GetSeepageRiverDrainFlux') ! TRANSOL gets the river flux
       ok = TRANSOL_GetSeepageRiverDrainFlux(XchModTranModDrainFlux,XchModTranModNID,XchMod2TranIdx,XchMod2TranOff,mv,'drn'); call DriverChk(ok,'TRANSOL_GetSeepageRiverDrainFlux') ! TRANSOL gets the drain flux
       XchModTranModDrainFlux = mv
       ok = MOZART_GetRiverDrainFlux(XchModMozModRiverFlux,XchMozNID,XchMod2MozIdx,XchMod2MozOff,mv,1); call DriverChk(ok,'MOZART_GetRiverDrainFlux') ! MOZART gets the river flux
       XchModMozModRiverFlux = mv
       ok = MOZART_GetRiverDrainFlux(XchModMozModRiverFluxWells,XchMozNID,XchMod2MozIdx,XchMod2MozOff,mv,2); call DriverChk(ok,'MOZART_GetRiverDrainFlux') ! MOZART gets the river flux (wells)
       XchModMozModRiverFluxWells = mv
       ok = MOZART_GetRiverDrainFlux(XchModMozModDrainFlux,XchMozNID,XchMod2MozIdx,XchMod2MozOff,mv,1); call DriverChk(ok,'MOZART_GetRiverDrainFlux') ! MOZART gets the drain flux
       XchModMozModDrainFlux = mv
       ok = MOZART_GetSalt(XchModMozModSalt,XchMozNID,XchMod2MozIdx,XchMod2MozOff,mv,2); call DriverChk(ok,'MOZART_GetSalt') ! MOZART get the salt flux
       XchModMozModSalt = mv
    end if
!##### BEGIN EXCHANGE: AfterFinishTimeStepMODFLOW #############################
    if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
       call osd_chdir(simwd2)
       call MetaSwap_finishTimestep(currentTime)
    end if

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

 endif

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

! ... next timestep
 if (exitcode.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) exitcode = -31
    if (exitcode.ne.0)     endOfSimulation=.true.
    ! 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 (exitcode.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) endOfSimulation=.true.
!       if (currentTime.eq.nodataTime) endOfSimulation=.true.
    endif
 endif

 if (rt.eq.rtmodsimtranmoz) then
    doTimeLoop = currenttime.lt.endOfCurrentTimeStepMozart .and. .not.endOfSimulation
 else
    doTimeLoop = .not.endOfSimulation
 end if

 if (exitcode.ne.0) doTimeLoop = .false.
   
end subroutine mf_finish_timestep

subroutine mf_init_simulation(endOfSimulation,convergedMozart,start_time)
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_init_simulation
 use mf_module
 use mod_pest, only: pest1_meteo_metaswap, pest1alpha_metaswap, pest1appendlogfile
 implicit none
 
 logical, intent(inout) :: endOfSimulation,convergedMozart
 double precision, intent(in) :: start_time
 
if (lrunfile) then
    if (nsub.gt.0) then
       write(modwd2,'(4a)') trim(modwd1),'\',trim(submstr(isub)),'\modflow\mf2005_tmp'
       if (usemetaswap) then
          write(simwd2,'(4a)') trim(simwd1),'\',trim(submstr(isub)),'\metaswap'
       end if
       write(*,'(50(''+''),/,1x,a,1x,a,1x,a,i3.3,a,i3.3,a,/,50(''+''))') 'Computing for submodel:',trim(submstr(isub)),'(',isub,'/',nsub,')'
    else
       write(modwd2,'(2a)') trim(modwd1),'\modflow\mf2005_tmp'
       if (usemetaswap) then
          write(simwd2,'(2a)') trim(simwd1),'\metaswap'
       end if
    end if
 else
    modwd2 = modwd1
    simwd2 = simwd1
 end if

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

 ! init return values
 retValMF2005      = 0

 ! time values
 currentTime       = nodataTime
 currentTimeMF2005 = 0.d0

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

 ! 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
    if (lipest) then
       ok = mf2005_GetPestFlag(lipest); call driverChk(ok,'mf2005_GetPestFlag')
    end if
 end if
 if (usemetaswap) then
    call osd_chdir(simwd2)
    if (lipest) then
       mf_nrow = 0; mf_ncol = 0; lpwt = .false.
       ok = mf2005_PutGridDimensions(mf_igrid, mf_nrow, mf_ncol, mf_nlay); call driverChk(ok,'mf2005_PutGridDimensions')
       ok = mf2005_PutPWTActive(mf_igrid, lpwt); call driverChk(ok,'mf2005_PutPWTActive')
       call pest1_meteo_metaswap()
       call pest1alpha_metaswap(mf_nrow,mf_ncol,lpwt)
    end if
    call MetaSwap_initComponent()
 end if
 if (usetransol) call TRANSOL_initComponent()
 if (usemozart) then
    call osd_chdir(mozwd)
    call mozart_initComponent(mozrecord)
 end if

 ! append the PEST log-file
 if (lipest) then
    call pest1appendlogfile()
    call pest1log()
 end if

 ! check if MODFLOW is activated
 if (retValMF2005.ne.0) call driverChk(.false.,'MODFLOW initialization failed.')
 
 write(*,*) '--> start_time=',start_time
 if (start_time > -0.01d0) 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
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
    call osd_chdir(simwd2)
    call metaswap_initSimulation(currentTime)
 end if
 call osd_chdir(modwd2)
 call mf2005_initSimulation(currentTime,retValMF2005)
 if (retValMF2005.ne.0) then
    exitcode = -12
    return
 endif
 
! #### 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)

 if (rt.eq.rtmodsimtranmoz) then
    call osd_chdir(mozwd)
    call mozart_initSimulation()
    convergedMozart = .false.
    endOfSimulation = .false.
    call mozart_prepareTimeStep(endOfSimulation,convergedMozart,currentTime) ! read signal file
 end if

! #### BEGIN EXCHANGE: BeforeTimestep #########################################
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
    ok = mf2005_PutSimulationType(mf_igrid, mf_steadystate); call driverChk(ok,'mf2005_PutSimulationType') ! MODFLOW puts the simulation type (transient/steady-state)
    if (mf_steadystate) call driverChk(ok,'MODFLOW may not be steady-state')
 end if
 ! Coupling MODFLOW-MetaSWAP
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
    ok = mf2005_PutModSimNumberOfIDs(mf_igrid, XchModSimModNID); call driverChk(ok,'mf2005_PutModSimNumberOfIDs') ! MODFLOW puts the number of MOD-SIM IDs
    ok = metaswap_PutModSimNumberOfIDs(XchModSimSimNID); call driverChk(ok,'metaswap_PutModSimNumberOfIDs')       ! MetaSWAP puts the number of MOD-SIM IDs
    call driverXchInitModSim() ! allocate and init

    ok = mf2005_PutModSimIDs(mf_igrid,XchModSimModIds); call driverChk(ok,'mf2005_PutModSimIDs')       ! MODFLOW puts the MOD-SIM exchange IDs
    ok = mf2005_PutModSimCells(mf_igrid,XchModSimModCells); call driverChk(ok,'mf2005_PutModSimCells') ! MODFLOW puts the MOD-SIM cells
    ok = metaswap_PutModSimIDs(XchModSimSimIds); call driverChk(ok,'metaswap_PutModSimIDs')            ! get exchange id's
 end if
 ! Coupling MODFLOW-MOZART/LSW, MODFLOW-MOZART/PV, MetaSWAP-MOZART/LSW
 if (rt.eq.rtmodsimtranmoz) then
    ! Main program puts the river systems to skip in the coupling
    ok = mf2005_PutModMozRiversToSkip(mf_igrid,nhrivsys,hrivsys); call driverChk(ok,'mf2005_PutModMozRiversToSkip')     ! Main program puts the river systems to skip in the coupling
    ok = mozart_PutModMozNumberOfIDs(XchMozNID); call driverChk(ok,'mozart_PutModMozNumberOfIDs')                       ! MOZART *puts* the number of LSW IDs
    ok = mozart_PutModMozPVNumberOfIDs(XchMozPVNID); call driverChk(ok,'mozart_PutModMozPVNumberOfIDs')                 ! MOZART *puts* the number of PV IDs
    ok = mf2005_PutModMozNumberOfIDs(mf_igrid, XchModMozModNID); call driverChk(ok,'mf2005_PutModMozNumberOfIDs')       ! MODFLOW *puts* the number of LSW IDs
    ok = mf2005_PutModMozPVNumberOfIDs(mf_igrid, XchModMozPVModNID); call driverChk(ok,'mf2005_PutModMozPVNumberOfIDs') ! MODFLOW *puts* the number of PV IDs
    ok = metaswap_PutModMozNumberOfIDs(XchSimMozSimNID); call driverChk(ok,'metaswap_PutModMozNumberOfIDs')             ! MetaSWAP *puts* the number of LSW IDs

    call driverXchInitModSimTranMoz() ! allocate and init

    ok = mozart_PutModMozIDs(XchMozIds,XchMozNID); call driverChk(ok,'mozart_PutModMozIDs')             ! MOZART *puts* the LSW IDs
    ok = mozart_PutModMozPVIDs(XchMozPVIds,XchMozPVNID); call driverChk(ok,'mozart_PutModMozPVIDs')     ! MOZART *puts* the PV IDs
    ok = mf2005_PutModMozIDs(mf_igrid, XchModMozModIds); call driverChk(ok,'mf2005_PutModMozIDs')       ! MODFLOW *puts* the LSW IDs
    ok = mf2005_PutModMozPVIDs(mf_igrid, XchModMozPVModIds); call driverChk(ok,'mf2005_PutModMozPVIDs') ! MODFLOW *puts* the PV IDs
    ok = mf2005_PutModMozCells(mf_igrid,XchModMozModCells); call driverChk(ok,'mf2005_PutModMozCells')  ! MODFLOW puts the MOD-MOZ cells
    ok = metaswap_PutModMozIDs(XchSimMozSimIds); call driverChk(ok,'metaswap_PutModMozIDs')             ! MetaSWAP *puts* the LSW IDs
 end if

 ! Create mappings
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) call driverXchIniMapModSim()
 if (rt.eq.rtmodsimtranmoz) then
    call driverXchIniMapModMoz()
    call driverXchIniMapModMozPV()
    call driverXchIniMapSimMoz()
    call driverXchIniMapModTran()
    call driverXchIniMapTranMoz()
 end if

! ###### END EXCHANGE: BeforeTimestep #########################################

! ... timestep
 timestep = 0

 tsc=0 ! debug: TimeStepCycle

 stsave    = .false.
 mozstsave = .false.
 strestore = .false.
 
end subroutine mf_init_simulation

subroutine mf_init_mozart
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_init_mozart
 use mf_module
 implicit none
 iterMozart = 0.
end subroutine mf_init_mozart

subroutine mf_init_timestep
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_init_timestep
 use mf_module
 
 implicit none
 
 if (rt.eq.rtmodsimtranmoz) then
    iterMozart = iterMozart + 1
    if (iterMozart.eq.1) mozstsave = .true.
    call mozart_getCurrentTime(currentTimeMozart) ! get the current time
    currentTime = currentTimeMozart
    call mozart_getEndOfCurrentTimeStep(endOfCurrentTimeStepMozart) ! get the end of the currentMozart time step
    call osd_chdir(mozwd)
    call mozart_initTimeStep(iterMozart)  ! initialize time step

!BEGIN EXCHANGE: BeginOfMozartTimestep ##################################

    ok = mozart_PutLSWLevels(XchModMozMozLevels,mv); call DriverChk(ok,'mozart_PutLSWLevels') ! MOZART puts the LSW levels
    ok = mozart_PutPVLevels(XchModMozPVMozPVLevels,mv); call DriverChk(ok,'mozart_PutPVLevels') ! MOZART puts the PV levels
    ok = mozart_PutLSWFractions(XchSimMozMozFractions,mv); call DriverChk(ok,'mozart_PutLSWFractions') ! MOZART puts the LSW fractions
    ok = mf2005_GetLSWLevels(mf_igrid,XchModMozMozLevels,&
                             XchModMozModNID,XchMoz2ModIdx,XchMoz2ModOff,mv); call DriverChk(ok,'mf2005_GetLSWLevels') ! MODFLOW gets LSW levels
    XchModMozMozLevels = mv
    ok = mf2005_GetPVLevels(mf_igrid,XchModMozPVMozPVLevels,&
                            XchModMozPVModNID,XchMozPV2ModIdx,XchMozPV2ModOff,mv); call DriverChk(ok,'mf2005_GetPVLevels') ! MODFLOW gets PV levels
    XchModMozPVMozPVLevels = mv
    ok = metaswap_GetFractions(XchSimMozMozFractions,XchSimMozSimNID,XchMoz2SimIdx,XchMoz2ModOff,mv); call DriverChk(ok,'metaswap_GetFractions') ! MetaSWAP gets the LSW fractions
    XchSimMozMozFractions = mv
!# END EXCHANGE: BeginOfMozartTimestep ##################################
 end if
  
end subroutine mf_init_timestep

subroutine mf_finish_simulation(endOfSimulation,convergedMozart)
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_finish_simulation
 use mf_module
 
 implicit none

 logical, intent(out) :: endOfSimulation,convergedMozart

!##### BEGIN EXCHANGE: BeforeTRANSOLIteration #################################
 if (rt.eq.rtmodsimtranmoz .and. exitcode.eq.0) then
    ok = mozart_PutLSWSalt(XchTranMozMozSalt,mv); call DriverChk(ok,'Error mozart_PutLSWSalt') ! MOZART puts the LSW salt concentrations
    ok = TRANSOL_GetSalt(XchTranMozMozSalt,&
                         XchSimMozSimNID,XchMoz2SimIdx,XchMoz2ModOff,mv); call DriverChk(ok,'Error TRANSOL_GetSalt') ! TRANSOL gets the LSW salt concentrations
    XchTranMozMozSalt = mv
 end if

!####### END EXCHANGE: BeforeTRANSOLIteration #################################

 if (rt.eq.rtmodsimtranmoz .and. exitcode.eq.0) then
    call TRANSOL_performIter()
 end if
!##### BEGIN EXCHANGE: AfterMozartIteration ###################################
 if (rt.eq.rtmodsimtranmoz .and. exitcode.eq.0) then
    ok = MetaSWAP_PutCumSWSprinklingDemandFluxes(XchSimMozSimCuSWSprinklingFlux,mv); call DriverChk(ok,'MetaSWAP_PutCumSWSprinklingDemandFluxes') ! MetaSWAP puts the cumulative surface water sprinkling demand fluxes
    ok = MetaSWAP_PutCumRunonFluxes(XchSimMozSimCuRunonFlux,mv); call DriverChk(ok,'MetaSWAP_PutCumRunonFluxes') ! MetaSWAP puts the cumulative runon fluxes
    ok = TRANSOL_PutSalt(XchTranMozTranSalt,mv); call DriverChk(ok,'TRANSOL_PutSalt') ! TRANSOL puts the concentrations
    ok = MOZART_GetCumSWSprinklingDemandFluxes(XchSimMozSimCuSWSprinklingFlux,&
                                            XchMozNID,XchSim2MozIdx,XchSim2MozOff,mv); call DriverChk(ok,'MOZART_GetCumSWSprinklingDemandFluxes') ! MOZART gets the surface water sprinkling demand fluxes
    XchSimMozSimCuSWSprinklingFlux = mv
    ok = MOZART_GetCumRunonFluxes(XchSimMozSimCuRunonFlux,&
                                  XchMozNID,XchSim2MozIdx,XchSim2MozOff,mv); call DriverChk(ok,'MOZART_GetCumRunonFluxes') ! MOZART gets the cumulative runon fluxes
    ok = MOZART_GetSalt(XchTranMozTranSalt,&
                        XchMozNID,XchTran2MozIdx,XchTran2MozOff,mv,1); call DriverChk(ok,'MOZART_GetSalt') ! MOZART gets the concentrations
    XchTranMozTranSalt = mv
 end if
!####### END EXCHANGE: AfterMozartIteration ###################################

 if (rt.eq.rtmodsimtranmoz .and. exitcode.eq.0) then
    call TRANSOL_finishTimeStep(currentTime)
    dtMozart = endOfCurrentTimeStepMozart - currentTimeMozart
    call osd_chdir(mozwd)
    call mozart_finishTimeStep(iterMozart,dtMozart,usetransol)
    call mozart_prepareTimeStep(endOfSimulation,convergedMozart,currentTime) ! read signal file
 end if

 ! make sure MODFLOW-MetaSWAP terminates normally
 if (rt.eq.rtmodsim) then
    convergedMozart = .true.
 end if
 
 if (exitcode.ne.0) endOfSimulation = .True.
 
end subroutine mf_finish_simulation

subroutine mf_finish_pest(convergedPest)
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_finish_pest
 use mf_module
 use imod_utl, only: imod_utl_closeunits
 use mod_pest, only: pestnext
 
 implicit none
 
 logical, intent(inout) :: convergedPest
 
! ... end
 call osd_chdir(modwd2)
 if (exitcode.eq.0)  call mf2005_finishSimulation(retValMF2005)
 if (retValMF2005.ne.0 ) exitcode = -61
 if (rt.eq.rtmodsim .or. rt.eq.rtmodsimtranmoz) then
    call osd_chdir(simwd2)
    call MetaSWAP_finishSimulation()
 end if
!#### 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

! next pest iteration
 if (.not.lipest) then
    convergedPest=.true.
 else
    convergedPest=pestnext()
 end if
 call imod_utl_closeunits() 
 
end subroutine mf_finish_pest

subroutine mf_set_timestep_done(done_with_current_timestep, Reset_time)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_set_timestep_done
 use mf_module
 implicit none
 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
 else
     ! time step done
 endif
end subroutine mf_set_timestep_done
 
!!!
!!! DLL-methods for coupling with WFLOW
!!!

function mf_init_no_sub_pest_moz(config_file_name, endOfSimulation, startTime,stsOverRule) result(ret_val)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_init_no_sub_pest_moz
 use mf_module
 implicit none
 integer              :: ret_val
 character(len=*)     :: config_file_name
 logical, intent(out) :: endOfSimulation
 logical, intent(in) :: stsOverRule
 double precision, intent(in) :: startTime
 logical              :: dummyConvergedMozart
 integer :: dummynnsub
 
 call mf_init(config_file_name,dummynnsub,stsOverRule)
 if (exitcode == 0) then
     endOfSimulation   = .false.
     call mf_init_simulation(endOfSimulation,dummyConvergedMozart,startTime)
 endif
 ret_val = exitcode

end function  mf_init_no_sub_pest_moz

subroutine mf_update_no_sub_pest_moz(endOfSimulation, doTimeLoop, stssaveOverRule)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_update_no_sub_pest_moz
 implicit none
 logical, intent(inout) :: endOfSimulation
 logical, intent(inout) :: doTimeLoop
 logical, intent(in) :: stssaveOverRule
 
 logical :: convergedIter
 call mf_init_iter(stssaveOverRule)
 convergedIter = .false.
 do while(.not.convergedIter) ! Picard iter loop
    call mf_iter(convergedIter)    
 end do
 call mf_finish_timestep(convergedIter,endOfSimulation,doTimeLoop)
end subroutine mf_update_no_sub_pest_moz

subroutine mf_finish_no_sub_pest_moz()
!DEC$ ATTRIBUTES DLLEXPORT :: mf_finish_no_sub_pest_moz
 implicit none
 logical :: dummyConvergedPest
 call mf_finish_pest(dummyConvergedPest)
end subroutine mf_finish_no_sub_pest_moz

subroutine mf_get_grid_sizes(nlay, nrow, ncol)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_get_grid_sizes
use mf_module
 implicit none
 integer, intent(out) :: nlay, nrow, ncol
 integer :: igrid = 1
 ok = mf2005_PutGridDimensions(igrid,nrow,ncol,nlay)  
end subroutine mf_get_grid_sizes
    
subroutine mf_set_recharge(recharge, ncol, nrow)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_set_recharge
use mf_module
  implicit none
  integer, intent(in)                        :: ncol, nrow
  real   , intent(in), dimension(ncol, nrow) :: recharge
  integer :: igrid = 1
  ok = mf2005_GetRecharge(recharge,ncol,nrow,igrid)
end subroutine mf_set_recharge
    
subroutine mf_get_heads(head,ncol,nrow)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_get_heads
  use mf_module 
  implicit none
  integer, intent(in)                                 :: ncol,nrow
  double precision, intent(out), dimension(ncol,nrow) :: head

  integer :: ilay = 1
  integer :: igrid = 1
  
  ok = mf2005_PutHeadsForLayer(head,ncol,nrow,ilay,igrid)
  
end subroutine mf_get_heads

!!!
!!! DLL-methods for coupling with RIBASIM
!!!

function mf_end_of_ribasim_timestep_reached(ribasim_end_of_timestep) result(end_of_ribasim_timestep_reached)
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_end_of_ribasim_timestep_reached
     use mf_module
     implicit none
 
     logical end_of_ribasim_timestep_reached
     double precision ribasim_end_of_timestep
 
     if(currenttime.ge.ribasim_end_of_timestep)then
         end_of_ribasim_timestep_reached = .true.
     else
         end_of_ribasim_timestep_reached = .false.
     endif
 end function

subroutine set_well_values(RibaQdemand, numWells)
!DEC$ ATTRIBUTES DLLEXPORT :: set_well_values
    use gwfmnw1module, only: demandbuff
    implicit none
    real, dimension(numWells), intent(in) :: RibaQdemand
    integer, intent(in) :: numWells

    integer, parameter :: igrid=1
    integer :: i
    
    call sgwf2mnw1pnt(igrid)
    
    if(.not.associated(demandbuff)) then
        allocate(demandbuff(numWells))
    end if    
    
    do i = 1, numWells
       demandbuff(i) = RibaQdemand(i)   
    end do
        
    call sgwf2mnw1psv(igrid)    
    
!ok = mf2005_GetWelldemand(RibaQdemand); call DriverChk(ok,'mf2005_GetQdemand') ! Modflow gets Qdemand 

end subroutine set_well_values

subroutine get_well_values(RibaQrealized,numWells)
!DEC$ ATTRIBUTES DLLEXPORT :: get_well_values
    use gwfmnw1module, only: well2
    implicit none
    real, dimension(numWells), intent(out) :: RibaQrealized
    integer, intent(in) :: numWells

    integer, parameter :: igrid=1 
    integer :: i
    
    call sgwf2mnw1pnt(igrid)
    
    do i = 1, numWells
        ribaqrealized(i) = well2(3,i)
    end do        
    
    !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

function mf_dll_get_number_of_mnwells() result(number_of_wells)
!DEC$ ATTRIBUTES DLLEXPORT :: mf_dll_get_number_of_mnwells
use gwfmnw1module, only:mxwel2
    implicit none
    
    integer, parameter :: igrid=1
    integer :: number_of_wells
    
    call sgwf2mnw1pnt(igrid)
    number_of_wells = mxwel2
end function mf_dll_get_number_of_mnwells

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

function mf_external_timestep_reached(external_end_of_timestep) result(external_timestep_reached)
 !DEC$ ATTRIBUTES DLLEXPORT :: mf_external_timestep_reached
     use mf_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 mf_external_timestep_reached

subroutine modflow_dll_set_timestep_done(done_with_current_timestep, Reset_time)
!DEC$ ATTRIBUTES DLLEXPORT :: modflow_dll_set_timestep_done
    use mf_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
 
 
 !======================
 ! BMI
 !======================
 
 module mf_bmi_module
    use iso_c_binding

    ! NOTE: these 'MAX'-sizes must be the same as those in
    !       https://github.com/openearth/bmi-python/blob/master/bmi/wrapper.py
    integer(c_int)            :: MAXDIMS = 6
    integer(c_int), parameter :: MAXSTRINGLEN = 1024
    
    ! various variables for calling modflow dll methods
    double precision            :: mf_start_time = -0.1d+0 ! <-0.01 = take starttime from model 
    logical                     :: endOfSimulation = .false.
    integer                     :: mf_ncol, mf_nrow, mf_nper ! TODO: use mf-dll-method
    double precision, &
          dimension(:,:), pointer :: mf_layer1_head
    real, dimension(:,:), pointer :: mf_recharge
    double precision :: mf_end_time
    real, dimension(:), pointer :: mf_delt

 contains
    
    pure function char_array_to_string(char_array, length)
        integer(c_int), intent(in) :: length
        character(c_char),intent(in) :: char_array(length)
        character(len=length) :: char_array_to_string
        integer :: i
        do i = 1, length
           char_array_to_string(i:i) = char_array(i)
        enddo
    end function char_array_to_string
    
    integer(c_int) pure function strlen(char_array)
        character(c_char), intent(in) :: char_array(MAXSTRINGLEN)
        integer :: inull, i
        strlen = 0
        do i = 1, size(char_array)
           if (char_array(i) .eq. C_NULL_CHAR) then
              strlen = i-1
              exit
           end if
        end do
    end function strlen

    pure function string_to_char_array(string) result(char_array)
      ! pass only trimmed strings to this one
      character(len=*), intent(in) :: string
      character(kind=c_char,len=1) :: char_array(MAXSTRINGLEN)
      integer :: i
      do i = 1, len(string)
         char_array(i) = string(i:i)
      enddo
      char_array(len(string)+1) = C_NULL_CHAR
    end function string_to_char_array

end module mf_bmi_module
 
   !> The initialize() function accepts a string argument that
  !! gives the name (and path) of its "main input file", called
  !! a configuration file. This function should perform all tasks
  !! that are to take place before entering the model's time loop.
  integer(c_int) function initialize(c_config_file) result(c_iresult) bind(C, name="initialize")
    !DEC$ ATTRIBUTES DLLEXPORT :: initialize
    
    use mf_bmi_module
    use global
    use rf2mf_module, only: dis, starttime
    
    ! return value (c_iresult == 0: everything ok)

    ! arguments
    character(kind=c_char),intent(in)    :: c_config_file(MAXSTRINGLEN)
    character(len=strlen(c_config_file)) :: config_file

    ! parameters
    integer, parameter :: igrid = 1
    
    ! locals
    integer                      :: ret_val
    character(len=MAXSTRINGLEN)  :: config_file_name
    logical                      :: enableStateSave
    integer                      :: retValMF2005
    double precision             :: tstart, factor
    integer                      :: iper, year, month, day, date
    integer                      :: perln, nsp
    real                         :: tsmlt
 
    ! get the config file name name
    config_file = char_array_to_string(c_config_file, strlen(c_config_file))

    ! Now we can initialize with the config_file
    
    endOfSimulation = .false.
    enableStateSave = .true.
    
    c_iresult = mf_init_no_sub_pest_moz (config_file, endOfSimulation, mf_start_time, enableStateSave)
    if (c_iresult == 0) then
        call mf2005_getCurrentTime(mf_start_time, ret_val)
        c_iresult = ret_val
    endif
    
    if (c_iresult == 0) then
        mf_ncol = ncol
        mf_nrow = nrow
        allocate(mf_layer1_head(ncol, nrow))
        allocate(mf_recharge(ncol, nrow))
    endif
    
    ! determine starttime/endtime
    read(starttime(1:4),*) year
    read(starttime(5:6),*) month
    read(starttime(7:8),*) day
    date = year*10000 + month*100 + day
    call cfn_datehms2mjd(date,0,0,0,tstart)
    factor = 1.d0
    if (itmuni.eq.1) factor = 1/86400.d0 ! seconds
    if (itmuni.eq.2) factor = 1/1440.d0  ! minutes
    if (itmuni.eq.3) factor = 1/24.d0    ! hours    
    
    allocate(mf_delt(dis%nper))
    one=1.
    mf_end_time = tstart  
    mf_nper = dis%nper
    do kper = 1, dis%nper 
       perln  = dis%sp(kper)%perlen
       nsp    = dis%sp(kper)%nstp
       tsmlt  = dis%sp(kper)%tsmult  
       mf_delt(kper) = perln/float(nsp)
       if(tsmlt.ne.one)then
          mf_delt(kper)=perln*(one-tsmlt)/(one-tsmlt)**nsp
       end if    
       mf_end_time = mf_end_time+mf_delt(kper)*factor ! [d]    
    end do     
    
  end function initialize

  subroutine get_start_time(t) bind(C, name="get_start_time")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_start_time

    use mf_bmi_module

    real(c_double), intent(out) :: t
    t = mf_start_time
  end subroutine get_start_time

  subroutine get_end_time(t) bind(C, name="get_end_time")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_end_time

    use mf_bmi_module

    real(c_double), intent(out) :: t
    t = mf_end_time
  end subroutine get_end_time

  subroutine get_time_step(dt) bind(C, name="get_time_step")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_time_step

    use mf_bmi_module

    real(c_double), intent(out) :: dt

    dt = -1d+0  ! TODO: implement
    
  end subroutine get_time_step

  subroutine get_current_time(t) bind(C, name="get_current_time")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_current_time

    use mf_bmi_module

    ! arguments
    real(c_double), intent(out) :: t

    ! locals
    double precision :: time
    integer :: ret_val

    call mf2005_getCurrentTime(time, ret_val)
    t = time
    c_iresult = ret_val

  end subroutine get_current_time

  subroutine get_time_units(c_unit)  bind(C, name="get_time_units")

    use mf_bmi_module

    ! arguments
    character(kind=c_char), intent(out) :: c_unit(MAXSTRINGLEN)
    
    ! locals
    character(len=MAXSTRINGLEN) :: unit

    unit = 'Modified Julian Date'
    
    c_unit = string_to_char_array(trim(unit))
  
  end subroutine get_time_units

  integer(c_int) function update(dt) result(c_iresult) bind(C, name="update")
    !DEC$ ATTRIBUTES DLLEXPORT :: update

    use mf_bmi_module

    real(c_double),intent(in) :: dt

    logical :: stssave
    logical :: doTimeLoop
    integer :: ret_val
    double precision :: tt,endTime, currentTime, nextTime
    
    stssave=.true.
    
    call mf2005_getCurrentTime(currentTime, ret_val)
    endTime = currentTime + dt
    nextTime = currentTime
    do while (nextTime < endTime)
        call mf_update_no_sub_pest_moz(endOfSimulation, doTimeLoop, stssave) 
        call mf2005_getBeginOfNextTimeStep(nextTime,ret_val) 
    end do
        
    c_iresult = 0
  end function update
  
  subroutine update_until(c_time) bind(C, name="update_until")

    use mf_bmi_module

    real(c_double),intent(in) :: c_time

    ! locals
    logical :: stssave
    logical :: doTimeLoop
    logical, external :: mf_external_timestep_reached
    double precision :: time

    time = c_time

    do while (.not. mf_external_timestep_reached(time) .and. .not. endOfSimulation)  
        stssave=.true.
        call mf_update_no_sub_pest_moz(endOfSimulation, doTimeLoop, stssave) 
    enddo

  end subroutine update_until

  subroutine get_var_type(c_var_name, c_type)  bind(C, name="get_var_type")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_var_type

    use mf_bmi_module

    character(kind=c_char), intent(in) :: c_var_name(*)
    character(kind=c_char), intent(out) :: c_type(MAXSTRINGLEN)
    character(len=MAXSTRINGLEN) :: type_name, var_name

    ! Use one of the following types
    ! BMI datatype        C datatype        NumPy datatype
    ! BMI_STRING          char*             S<
    ! BMI_INT             int               int16
    ! BMI_LONG            long int          int32
    ! BMI_FLOAT           float             float32
    ! BMI_DOUBLE          double            float64
    var_name = char_array_to_string(c_var_name, strlen(c_var_name))

    select case(var_name)
    case("head")
        type_name = "double"
    case("recharge","delt")
        type_name = "float"
    end select

    c_type = string_to_char_array(trim(type_name))

  end subroutine get_var_type

subroutine get_var(c_var_name, x) bind(C, name="get_var")
  !DEC$ ATTRIBUTES DLLEXPORT :: get_var
  ! Return a pointer to the variable

  use mf_bmi_module

  character(kind=c_char), intent(in) :: c_var_name(*)
  type(c_ptr), intent(inout) :: x
  integer(c_int), target, allocatable, save :: xi(:,:)

  integer :: ilay = 1
  integer :: igrid = 1
  logical :: ret_val

  ! The fortran name of the attribute name
  character(len=strlen(c_var_name)) :: var_name
  ! Store the name
  var_name = char_array_to_string(c_var_name, strlen(c_var_name))

  select case(var_name)
  case("head")
     ret_val = mf2005_PutHeadsForLayer(mf_layer1_head, mf_ncol, mf_nrow, ilay, igrid)
     x = c_loc(mf_layer1_head)
  case("delt")
     x = c_loc(mf_delt) 
  end select

end subroutine get_var

subroutine set_var(c_var_name, xptr) bind(C, name="set_var")
  !DEC$ ATTRIBUTES DLLEXPORT :: set_var
  ! Return a pointer to the variable

  use mf_bmi_module

  character(kind=c_char), intent(in) :: c_var_name(*)
  type(c_ptr), value, intent(in) :: xptr

  real(c_float), pointer :: x_2d_float_ptr(:,:)

  integer :: igrid = 1
  logical :: ret_val

  ! The fortran name of the attribute name
  character(len=strlen(c_var_name)) :: var_name
  ! Store the name
  var_name = char_array_to_string(c_var_name, strlen(c_var_name))

  select case(var_name)
  case("recharge")
     call c_f_pointer(xptr, x_2d_float_ptr, shape(mf_recharge))
     mf_recharge(:,:) = x_2d_float_ptr   
     ret_val = mf2005_GetRecharge(mf_recharge, mf_ncol, mf_nrow, igrid)
     write(*,*) 'Done setting recharge.'
  end select

end subroutine set_var

  subroutine get_var_rank(c_var_name, rank) bind(C, name="get_var_rank")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_var_rank

    use mf_bmi_module

    use iso_c_binding, only: c_int, c_char
    character(kind=c_char), intent(in) :: c_var_name(*)
    integer(c_int), intent(out) :: rank

    ! The fortran name of the attribute name
    character(len=strlen(c_var_name)) :: var_name
    ! Store the name
    var_name = char_array_to_string(c_var_name, strlen(c_var_name))

    select case(var_name)
    case("delt")
       rank = 1 
    case default    
       rank = 2
    end select
    
  end subroutine get_var_rank
  
  subroutine get_var_shape(c_var_name, shape) bind(C, name="get_var_shape")
    !DEC$ ATTRIBUTES DLLEXPORT :: get_var_shape

    use mf_bmi_module

    use iso_c_binding, only: c_int, c_char, c_loc

    character(kind=c_char), intent(in) :: c_var_name(*)
    integer(c_int), intent(inout) :: shape(MAXDIMS)
    
    character(len=strlen(c_var_name)) :: var_name

    var_name = char_array_to_string(c_var_name, strlen(c_var_name))
    shape = (/0, 0, 0, 0, 0, 0/)

    !shape(1) = mf_ncol
    !shape(2) = mf_nrow
    select case(var_name)
    case("delt")
       shape(1) = mf_nper
    case default    
       shape(1) = mf_nrow
       shape(2) = mf_ncol
    end select
    
  end subroutine get_var_shape
    
  integer(c_int) function finalize() result(c_iresult) bind(C, name="finalize")
    !DEC$ ATTRIBUTES DLLEXPORT :: finalize

    use mf_bmi_module

    call mf_finish_no_sub_pest_moz()
    c_iresult = 0
  end function finalize