module coupler

  !-  modules:
  use cmod_mem, only  : cmr1, cmr2, cmi1

  !-  directives:
  implicit none
  private

  !-  public coupler methods:
  public :: coupler_usage, coupler_init, coupler_dynamics, coupler_finish
  logical,    public,  parameter :: build_with_coupler = .true.

  !-  coupler index variables:
  integer(4), private, save :: nlow = 1, nup = 1

  !-  coupler time variables:
  integer(4), private, save :: out_nt
  real(8),    private, save :: dtbio

  ! These must only be modified by MASTER thread, but can be used by all:
  integer(4), private, save :: itime = 1, otime = 1

  !-  global ergom variables:
  type(cmr2), allocatable, private, save :: benthos(:)

  !-  local ergom variables:
  type(cmr1), allocatable, private, save :: npr_l(:)
  type(cmr2), allocatable, private, save :: benthos_l(:)
  type(cmr1), allocatable, private, save :: light_l(:)
  type(cmr1), allocatable, private, save :: chl_l(:)
  type(cmr1), allocatable, private, save :: kpar_l(:)
  type(cmr1), allocatable, private, save :: secchi_l(:)
  integer(4), allocatable, private, save :: nrivsegments(:)
  type(cmi1), allocatable, private, save :: rivmap(:) 
  type(cmr1), allocatable, private, save :: srfammflx_l(:)
  type(cmr1), allocatable, private, save :: srfnitflx_l(:)
  type(cmr1), allocatable, private, save :: srfphosflx_l(:)
  type(cmr1), allocatable, private, save :: srfamm_with_totatmos_Nflx_l(:)
  type(cmr1), allocatable, private, save :: srfnit_with_totatmos_Nflx_l(:)
  type(cmr1), allocatable, private, save :: srfamm_with_agri_Nflx_l(:)
  type(cmr1), allocatable, private, save :: srfnit_with_agri_Nflx_l(:)
  
  !-  private methods:
  private :: bio_limit

  !NOTES:
  ! * The implementation of the variable NP ratio is a quick-and-dirty solution.
  !   It should only allocate a 2d array if it is used, and we should try to 
  !   simply the IO, scatter and allocatation calls
  !
  ! * Consider allocating all arrays in cmod_init
contains

!===============================================================================

  subroutine coupler_usage(iu06,nt,nt2d)
    !- modules
    use bioparam,  only : ReadBioparam, ntbio, ntben, do_diag_chl, do_diag_pp, &
                          do_diag_secchi
    use constants, only : ldebug

    implicit none

    ! args
    integer(4), intent(in)  :: iu06
    integer(4), intent(out) :: nt    ! number of traces to be inc. in cmp array
    integer(4), intent(out) :: nt2d  ! number of benthos tracers
    
    ! locals
    character(len=11) :: cfmt, cfmt2
    integer(4)        :: ntdiag
    
    !- get the ergom parameters:
    call ReadBioparam(iu06,ldebug)

    cfmt  = '(a36,'//'i3,'//'a8)'
    cfmt2 = '(a4,'//'i2,'//'a23)'
    
    ntdiag = 1 !light
    if (do_diag_chl)    ntdiag = ntdiag + 1
    if (do_diag_pp)     ntdiag = ntdiag + 1
    if (do_diag_secchi) ntdiag = ntdiag + 1
    
    write(iu06,trim(cfmt))                                                     &
    'Running with dmi_ergom_coupler using', ntbio, ' tracers'
    write(iu06,trim(cfmt2)) ' and', ntben, ' sediment variable(s)  '
    write(iu06,trim(cfmt2)) ' and', ntdiag,' diagnostic variable(s)'

    nt2d = ntben
    nt   = ntbio
  
  end subroutine coupler_usage

!===============================================================================

  subroutine coupler_init(ia, narea, iu06, mainarea, ncoff, n2d, n3d,          &
                          mmx, nmx, kmx, nq, nz, kh, components,               &
                          bndq, bndz, bndzk, on_IO_task, iw2l, halo2size, iw3l,&
                          me1, mcol, me_l, mcol_l, kh_l, irc, t, s )
    ! Container for reading in initial values into the arrays and reading
    ! in boundaries (river, open). 
    ! Also, allocate/associate arrays for the DMI-ERGOM coupler so that we
    ! avoid messing with the cmp array.

    use bioinput,       only : bioread, biobndcl, bioread_npr, bioatmdep
    use biooutput,      only : biofile_new, val_benthos
    use bioparam,       only : bioInd,ntbio,ntben,biooutdt,NDTBIOMOD,          &
                               do_diag_chl, do_diag_secchi   
   !use io_control,     only : io_server, ioc_open_coupler, ioc_readrestart,   &
   !                           io_readrestart
    use bioriver,       only : bioloadings, read_nrivseg, assign_rivermap
    use cmod_miscvars,  only : dtmain, schedule, lvalidate !, tima, time
    use dmi_mpi,        only : mpi_size, dmpi_scatter, mpi_io_serial
    use dmi_mpi_global, only : mpi_io_rank
   
    implicit none

    integer(4), intent(in)    :: ia, narea, iu06, mmx, nmx, kmx, nq, nz, ncoff
    logical,    intent(in)    :: mainarea
    integer(4), intent(in)    :: n2d, n3d, kh(0:)
    integer(4), intent(in)    :: iw2l, halo2size, iw3l
    integer(4), intent(in)    :: me1(0:,0:), mcol(0:)
    integer(4), intent(in)    :: me_l(0:,0:), mcol_l(0:), kh_l(0:)
    real(8),    intent(in)    :: t(0:), s(0:) 
    real(8),    intent(inout) :: bndq(:,:,0:), bndz(:,:,0:), bndzk(:,:,0:)
    logical,    intent(in)    :: on_IO_task
    integer(4), intent(out)   :: irc

    real(8),    intent(inout) :: components(0:,1:)

    logical,    save :: firsttime = .true.
    integer(4)       :: n3s, lastts
    real(8)          :: npr(0:n2d)
    real(8)          :: srfammflx(0:n2d)
    real(8)          :: srfnitflx(0:n2d)
    real(8)          :: srfphosflx(0:n2d)
    real(8)          :: srfamm_with_totatmos_Nflx(0:n2d)
    real(8)          :: srfnit_with_totatmos_Nflx(0:n2d)
    real(8)          :: srfamm_with_agri_Nflx(0:n2d)
    real(8)          :: srfnit_with_agri_Nflx(0:n2d)
    real(8), parameter :: onehour = 3600.0_8

    irc = 0
    if (firsttime) then
      !- some private vars:
      dtbio  = NDTBIOMOD*dtmain    ! FIXME: still no stability check here
      out_nt = nint(biooutdt/dtmain)

      !- try to capture obvious time-setting bugs .-----------------------------
      lastts = nint(onehour*real(schedule,8)/dtmain,4)
      if (out_nt*(lastts/out_nt) /= lastts) then
        write(iu06,*)'Error: Inconsistent combination of biooutdt and schedule.'
        write(iu06,*)'       Modify biooutdt and/or schedule, and re-run.'
        irc = 1
      endif

      !- some index stuff:
      nlow = ncoff + 1
      nup  = ncoff + ntbio

      !- allocate the benthos variable:
      allocate( benthos(narea), benthos_l(narea), npr_l(narea), light_l(narea),&
                chl_l(narea),                                                  &
                srfammflx_l(narea),                                            &
                srfnitflx_l(narea),                                            &
                srfphosflx_l(narea),                                            &
                srfamm_with_totatmos_Nflx_l(narea),                                            &
                srfnit_with_totatmos_Nflx_l(narea),                                            &
                srfamm_with_agri_Nflx_l(narea),                                            &
                srfnit_with_agri_Nflx_l(narea),                                            &
                kpar_l(narea), secchi_l(narea))
      
      !- tracer index
      call bioInd()

      !- open output file on files entry:
      !if (io_server) then
        ! open outout files with io-server
      !  call ioc_open_coupler(tima, time)
      !else
        if (on_IO_task) call biofile_new ()
      !endif


      allocate(nrivsegments(narea), rivmap(narea))
      nrivsegments(1:narea) = 0
      
      !read number of deactivted river segments in area
      call read_nrivseg(narea,nrivsegments) 
     
      firsttime = .false.
    endif

    !- allocate/associate the biological variables:
    ! FIXME: if IO server is active, these arrays are only necesarry for 
    !        validation
    allocate( benthos(ia)%p(0:n2d,1:ntben))
    
    n3s = iw3l + halo2size
    allocate( benthos_l(ia)%p(0:iw2l,1:ntben), npr_l(ia)%p(0:iw2l),            &
              rivmap(ia)%p(0:iw2l), light_l(ia)%p(0:n3s),                      &
              srfammflx_l(ia)%p(0:iw2l),                                     &
              srfnitflx_l(ia)%p(0:iw2l),                                     &
              srfphosflx_l(ia)%p(0:iw2l),                                     &
              srfamm_with_totatmos_Nflx_l(ia)%p(0:iw2l),                                     &
              srfnit_with_totatmos_Nflx_l(ia)%p(0:iw2l),                                     &
              srfamm_with_agri_Nflx_l(ia)%p(0:iw2l),                                     &
              srfnit_with_agri_Nflx_l(ia)%p(0:iw2l),                                     &
              kpar_l(ia)%p(0:n3s))
    if (do_diag_chl) then
      allocate(chl_l(ia)%p(0:n3s))
    else
      allocate(chl_l(ia)%p(0:0))
    endif

    if (do_diag_secchi) then
      allocate(secchi_l(ia)%p(0:n3s))
    else
      allocate(secchi_l(ia)%p(0:0))
    endif

    ! Initialise rivmap
    ! FIXME: This should be done at a different place in a different way
    rivmap(ia)%p(0)  = 1
    rivmap(ia)%p(1:) = 1
     
    !- read restart without IO server
   !if (.not. io_server .or. (.not. io_readrestart)) then
      if (on_IO_task) then
        !- read and permute the initial fields:
        call bioread( ia, n3d, n2d, mmx, nmx, kmx, nz, me1, mcol, kh, narea,   &
                      iu06, benthos(ia)%p, ntben, components, t, s, bndz)

        !- print initial stats of the benthos variable:
        if (lvalidate .and. mpi_io_serial) then
          call val_benthos(ia, iu06, n2d, benthos(ia)%p, ntben, (ia==1))
        endif

        ! read spatially varying NP ratio if required
        call bioread_npr(ia, narea, iu06, npr, mmx, nmx, n2d, me1 )

        ! read spatially varying atmosph. depos. if required
        call bioatmdep(iu06, ia, nmx, mmx,    &
                       srfammflx,             &
                       srfnitflx,             &
                       srfphosflx,             &
                       srfamm_with_totatmos_Nflx,             &
                       srfnit_with_totatmos_Nflx,             &
                       srfamm_with_agri_Nflx,             &
                       srfnit_with_agri_Nflx,             &
                       me1)

        if  (mainarea) then
          !- boundary values:
! FIXME: +2 is an intermediate soulution...a better one is needed
          call biobndcl(kmx, kh, me1, mcol, nz, nlow+2, nup+2, bndzk, iu06,    &
                        components)
        endif
      endif
  
      !- FIXME: task_global2local instead of scatter would be nicer and better
      !-        fit to main HBM code
      !- scatter global2local
      if (mpi_size == 1) then
        ! local == global in serial mode
        benthos_l(ia)%p(1:n2d,1:ntben)       = benthos(ia)%p(1:n2d,1:ntben)
        npr_l(ia)%p(1:n2d)                   = npr(1:n2d)
        srfammflx_l(ia)%p(1:n2d)             = srfammflx(1:n2d)
        srfnitflx_l(ia)%p(1:n2d)             = srfnitflx(1:n2d)
        srfphosflx_l(ia)%p(1:n2d)             = srfphosflx(1:n2d)
        srfamm_with_totatmos_Nflx_l(ia)%p(1:n2d)             = srfamm_with_totatmos_Nflx(1:n2d)
        srfnit_with_totatmos_Nflx_l(ia)%p(1:n2d)             = srfnit_with_totatmos_Nflx(1:n2d)
        srfamm_with_agri_Nflx_l(ia)%p(1:n2d)             = srfamm_with_agri_Nflx(1:n2d)
        srfnit_with_agri_Nflx_l(ia)%p(1:n2d)             = srfnit_with_agri_Nflx(1:n2d)
      else
        call dmpi_scatter(ia,n2d,(1),me_l,mcol_l,benthos_l(ia)%p,kh_l,me1,mcol,          &
                          benthos(ia)%p,ntben,mpi_io_rank,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,npr_l(ia)%p,me1,mcol,                   &
                          npr,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfammflx_l(ia)%p,me1,mcol,             &
                          srfammflx,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfnitflx_l(ia)%p,me1,mcol,             &
                          srfnitflx,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfphosflx_l(ia)%p,me1,mcol,             &
                          srfphosflx,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfamm_with_totatmos_Nflx_l(ia)%p,me1,mcol,             &
                          srfamm_with_totatmos_Nflx,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfnit_with_totatmos_Nflx_l(ia)%p,me1,mcol,             &
                          srfnit_with_totatmos_Nflx,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfamm_with_agri_Nflx_l(ia)%p,me1,mcol,             &
                          srfamm_with_agri_Nflx,mpi_io_rank,kh_l,(.false.))
        call dmpi_scatter(ia,(1),n2d,me_l,mcol_l,srfnit_with_agri_Nflx_l(ia)%p,me1,mcol,             &
                          srfnit_with_agri_Nflx,mpi_io_rank,kh_l,(.false.))
      endif

    !elseif (io_server .and. ia==narea) then
      ! should only be called once...it loops through narea in the routine
    !  call ioc_readrestart(tima, .true., benthos_l)
   !endif

    !- river loadings:
! FIXME: +2 is an intermediate soulution...a better one is needed
    call bioloadings ( ia, nlow+2, nup+2, nq, bndq, kmx, me1, components )

    !get riversegment coordinates for specific area
    if (nrivsegments(ia)>0) then
      call assign_rivermap(nrivsegments(ia),ia,narea,rivmap(ia)%p,me_l)
    endif

    !- that's all, folks!

  end subroutine coupler_init

!===============================================================================

  subroutine coupler_dynamics(ia, iu06, mmx, nmx, kmx, iw2l, msrf, mcol,       &
                              kh, nz, h, t, s, components, dispv, u, v,        &
                              casus, humid, cloud, wu, wv,                     &
                              update_counters, idx, irc, n2d, n3d, msrfg,      &
                              mcolg, khg, lomp, bndz)
    ! Container for the only call in the timeloop. If the coupler needs to do
    ! different things at different timesteps, e.g. dynamics and output, then
    ! it's up to the coupler_dynamics implementation to do this differentiation.

    use dmi_mpi_global,    only : mpi_io_rank
    use dmi_mpi,           only : dmpi_gather, mpi_io_serial, mpi_size
    use dmi_omp,           only : domp_get_thread_no
    use cmod_local_arrays, only : ind_l
    use cmod_arrays,       only : xkoor, ykoor, ptc, ndtp, wrk
    use cmod_miscvars,     only : izeite, i2, lvalidate, dtday, schfactor,     &
                                  time !, ctim
    use cmod_timer,        only : timer_print
    use constants,         only : spinup, archive, izeits, tim
   !use io_control,        only : io_server, ioc_write_coupler,                &
   !                              ioc_writerestart_coupler
    use biooutput,         only : biofile_next, biofile_write,                 &
                                  write_bio_restart, val_benthos, diag_out,    &
                                  bioncdf
    use bioparam,          only : NDTBIOMOD, ntbio, ntben, do_diag_chl,        &
                                  didx_chloro, didx_secchi,   do_diag_secchi
    use biomod,            only : bio_dynamics
    use biodiagnostic,     only : calc_chlorophyll, calc_secchi_depth

    implicit none

    integer(4), intent(in)    :: iw2l
    integer(4), intent(in)    :: msrf(0:,0:), mcol(0:)
    integer(4), intent(in)    :: ia, iu06, mmx, nmx, kmx, nz, kh(0:)
    real(8),    intent(in)    :: h(0:), u(0:), v(0:)
    real(8),    intent(in)    :: wu(0:), wv(0:), humid(0:), cloud(0:)
    logical,    intent(in)    :: casus(0:), update_counters, lomp
    integer(4), intent(inout) :: idx(1:,1:)
    integer(4), intent(out)   :: irc
    real(8),    intent(in)    :: dispv(0:)
    real(8),    intent(in)    :: t(0:), s(0:)
    real(8),    intent(in)    :: bndz(:,:,0:)

    real(8),    intent(inout) :: components(0:,1:)

    !- Global fields - only need for output - maybe we should do output in a
    !- separate routine
    integer(4), intent(in)    :: n2d,n3d
    integer(4), intent(in)    :: msrfg(0:,0:), mcolg(0:), khg(0:)

    integer(4)       :: tnum

    ! these must only be used on MASTER:
    integer(4), save :: ntime = 0
    logical          :: did_gather

    irc = 0

    ! do the ergom dynamics ----------------------------------------------------
    if (mod(itime,NDTBIOMOD) == 0 .and. lomp) then
      call domp_get_thread_no (tnum)
      ndtp(tnum) = 0

      !- check eco-var limits:
      call bio_limit(iw2l, mcol, kh, components,idx)

      !- do the bio dynamics:
      ! fixme: we should not send full npr_l array if it is not used
      call bio_dynamics(iw2l, kmx, msrf, mcol, ind_l(ia)%p, kh, h,u, v, wu, wv,&
                        t, s, humid, cloud, components, benthos_l(ia)%p,       &
                        xkoor(ia)%p, ykoor(ia)%p, dispv, casus, npr_l(ia)%p,   &
                        dtbio, idx, ndtp(tnum), rivmap(ia)%p, light_l(ia)%p,   &
                        kpar_l(ia)%p                                           &
                        , srfammflx_l(ia)%p                                    &
                        , srfnitflx_l(ia)%p                                    &
                        , srfphosflx_l(ia)%p                                    &
                        , srfamm_with_totatmos_Nflx_l(ia)%p                                    &
                        , srfnit_with_totatmos_Nflx_l(ia)%p                                    &
                        , srfamm_with_agri_Nflx_l(ia)%p                                    &
                        , srfnit_with_agri_Nflx_l(ia)%p                                    &
                        )
                         
      !call only when needed
      if ((do_diag_chl).and.(mod(otime,out_nt) == 0))then
        call calc_chlorophyll(iw2l, mcol, kh, idx, components, light_l(ia)%p,  &
                              chl_l(ia)%p)
      endif
      !FIXME: Empty call so far
      if ((do_diag_secchi) .and.(mod(otime,out_nt) == 0))then
        call calc_secchi_depth (iw2l, kh, idx, kpar_l(ia)%p, secchi_l(ia)%p)
      endif

!$OMP BARRIER
      irc = maxval( ndtp )
      if (irc /= 0) then
        if (ndtp(tnum) /= 0) then
          write(iu06,'(a38,i3)')                                               &
            ' Error in bio_dynamics, error flag was ',ndtp(tnum)
          write(iu06,'(a38,i3)') 'Error in domain                        ', ia
        endif
      endif
    endif

    if (lomp) return
    ! PLEASE NOTE: SERIAL BELOW HERE

    ! do the ergom outout ------------------------------------------------------
   !if (io_server) then
      ! writing with IO-server
   !  call ioc_write_coupler(izeits, ia, components, benthos_l(ia)%p)
   !else
      ! writing without IO-server
      did_gather = .false.
      if (mod(otime,out_nt) == 0) then
        
        if (mpi_size == 1) then
          ! global==local in serial mode
          benthos(ia)%p(1:n2d,1:ntben) = benthos_l(ia)%p(1:n2d,1:ntben)
          ! fixme: maybe we don't need this one:
          ptc(ia)%p(1:,1:) = components(1:,1:)
        else
          call dmpi_gather(ia,components,ptc(ia)%p,ntbio,mpi_io_rank,kmx)
          call dmpi_gather(ia,benthos_l(ia)%p,benthos(ia)%p,ntben,mpi_io_rank)
        endif

        did_gather = .true.
        if (mpi_io_serial) then
          call biofile_write(ia, iu06, mmx, nmx, kmx, n2d, n3d, msrfg, mcolg,  &
                             khg, nlow, ptc(ia)%p, benthos(ia)%p)
        endif

        !output diagnostic variables, only active with netcdf output:
        if (bioncdf) then
          !LIGHT
          if (mpi_size == 1) then
            wrk(ia)%p = light_l(ia)%p
          else
            call dmpi_gather(ia,kmx,msrf,mcol,kh,light_l(ia)%p,wrk(ia)%p,      &
                             mpi_io_rank)
          endif
          if (mpi_io_serial) then
            call diag_out( ia, 1, wrk(ia)%p, msrfg, mcolg, khg, mmx, nmx, kmx)
          endif

          !CHLOROPHYLL
          if (do_diag_chl) then
            if (mpi_size == 1) then
              wrk(ia)%p = chl_l(ia)%p
            else
              call dmpi_gather(ia, kmx, msrf, mcol, kh, chl_l(ia)%p, wrk(ia)%p,&
                               mpi_io_rank)
            endif
            if (mpi_io_serial) then
              call diag_out(ia, didx_chloro, wrk(ia)%p, msrfg, mcolg, khg, mmx,&
                            nmx, kmx)
            endif
          endif

          !SECCHI DEPTH
          if (do_diag_secchi) then
            if (mpi_size == 1) then
              wrk(ia)%p = secchi_l(ia)%p
            else
              call dmpi_gather(ia, kmx, msrf, mcol, kh, secchi_l(ia)%p,        &
                               wrk(ia)%p, mpi_io_rank)
            endif
            !FIXME: Nothing happening yet
            if (mpi_io_serial) then
              call diag_out(ia, didx_secchi, wrk(ia)%p, msrfg, mcolg, khg, mmx,&
                            nmx, kmx)
            endif
          endif
        endif

      endif
   !endif

    ! update time step counters ------------------------------------------------
    if (update_counters) then
      if (mod(itime,NDTBIOMOD) == 0) then
        itime = 1
      else
        itime = itime + 1
      endif
      if (mod(otime,out_nt) == 0) then
        otime = 1
      else
        otime = otime + 1
      endif
      ntime = ntime + 1
      if (ntime == izeite .and. (tim+dtday-dtday*schfactor <= time)) then
        ! open the next biodat file. 
        !   biodat files in numbered 6 or 12 hour chunks:  
        !       the first file is 'biodat.00', and subsequent files are 
        !       'biodat.01', 'biodat.02', ... .  .
       !if (mpi_io_serial .and. .not. io_server) call biofile_next ( iu06 )
        call biofile_next ( iu06 )
        ntime = 0
      endif
    endif

    if (izeits == izeite) then
      ! print some timings
      call timer_print('TS - biobenthic took:','biobenthic',50)
      call timer_print('TS - biosolver took:','biosolver',50)

      ! print stats of the benthos variable ------------------------------------
      if (.not.did_gather) then
        if (mpi_size == 1) then
          ! global==local in serial mode
          benthos(ia)%p(1:n2d,1:ntben) = benthos_l(ia)%p(1:n2d,1:ntben)
        else
          call dmpi_gather(ia,benthos_l(ia)%p,benthos(ia)%p,ntben,mpi_io_rank)
        endif
      endif
      if (lvalidate .and. mpi_io_serial) then
        call val_benthos(ia, iu06, n2d, benthos(ia)%p,ntben, (ia==1))
      endif

      ! writing restart --------------------------------------------------------
      if ((i2 == 1 .and. archive) .or. (spinup .and. (.not.archive)) .or.      &
          (tim+dtday-dtday*schfactor > time .and. (.not.archive))) then
   !    if (io_server) then
          !- let IO-server write the restart file
   !      call ioc_writerestart_coupler(ctim, ia, components(0:,1:),           &
   !                                    benthos_l(ia)%p(:,:))

   !    elseif (.not. io_server) then
          ! write bio-restart file with mpi_io_serial

          if (.not.did_gather) then
            if (mpi_size == 1) then
              ! global==local in serial mode
              benthos(ia)%p(1:n2d,1:ntben) = benthos_l(ia)%p(1:n2d,1:ntben)
              ! fixme: maybe we don't need this one:
              ptc(ia)%p(1:,1:) = components(1:,1:)
            else
              call dmpi_gather(ia,components,ptc(ia)%p,ntbio,mpi_io_rank,kmx)
              call dmpi_gather(ia,benthos_l(ia)%p,benthos(ia)%p,ntben,         &
                               mpi_io_rank)
            endif
          endif
          if (mpi_io_serial) then
            call write_bio_restart(ia, iu06, update_counters, mmx, nmx, kmx,   &
                                   n2d, n3d, nz, msrfg, mcolg, khg, nlow,      &
                                   ptc(ia)%p, benthos(ia)%p, bndz)
          endif
   !    endif
      endif
    endif

  end subroutine coupler_dynamics

!===============================================================================

  subroutine coupler_finish(irc)
    ! Container for the call after the timeloop. Can be used to close files and
    ! deallocate memory and to do final outputs.
    use biooutput,  only : biofile_close
   !use io_control, only : io_server, ioc_close_coupler
    implicit none
    integer(4), intent(out)  :: irc

    !- close output file:
   !if (io_server) then
   !  call ioc_close_coupler()
   !else
      call biofile_close()
   !endif

    irc = 0
  end subroutine coupler_finish

!===============================================================================

  subroutine bio_limit( n2d, mcol, kh, ptc, idx)

    use dmi_omp,   only : domp_get_domain
    use bioparam,  only : zero,                  &
                          idx_amm            ,   &
                          idx_nit            ,   &
                          idx_phos           ,   &
                          idx_sil            ,   &
                          idx_dia            ,   &
                          idx_flag           ,   &
                          idx_cyano          ,   &
                          idx_mez            ,   &
                          idx_miz            ,   &
                          idx_det            ,   &
                          idx_dets           ,   &
                          idx_ldon           ,   &
                          ntbio


    implicit none

    integer(4), intent(in)    :: n2d, kh(0:), mcol(0:)
    real(8),    intent(inout) :: ptc(0:,1:)
    integer(4), intent(inout) :: idx(1:,1:)

    integer(4)         :: n, nl, nu, k, mi, nupm1, kb, mi0, mil, miu

    call domp_get_domain( kh, 1, n2d, nl, nu, idx )

    ! avoid cutting oxygene 'coz unlike the physical oxygene, the model 
    ! state variable called oxygene can indeed have negative values

    ptc(nl:nu,idx_amm) = max( zero, ptc(nl:nu,idx_amm) )
    ptc(nl:nu,idx_nit) = max( zero, ptc(nl:nu,idx_nit) )
    ptc(nl:nu,idx_phos) = max( zero, ptc(nl:nu,idx_phos) )
    ptc(nl:nu,idx_sil) = max( zero, ptc(nl:nu,idx_sil) )
    ptc(nl:nu,idx_dia) = max( zero, ptc(nl:nu,idx_dia) )
    ptc(nl:nu,idx_flag) = max( zero, ptc(nl:nu,idx_flag) )
    ptc(nl:nu,idx_cyano) = max( zero, ptc(nl:nu,idx_cyano) )
    ptc(nl:nu,idx_mez) = max( zero, ptc(nl:nu,idx_mez) )
    ptc(nl:nu,idx_miz) = max( zero, ptc(nl:nu,idx_miz) )
    ptc(nl:nu,idx_det) = max( zero, ptc(nl:nu,idx_det) )
    ptc(nl:nu,idx_dets) = max( zero, ptc(nl:nu,idx_dets) )
    ptc(nl:nu,idx_ldon) = max( zero, ptc(nl:nu,idx_ldon) )

    do n = nl, nu
      kb = kh(n)
      if (kb <= 1) cycle
      ! mi0 = mcol(n)
      ! do k = 2,kb
      !   mi = mi0+k-2
      !   ptc(mi,idx_amm) = max( zero, ptc(mi,idx_amm) )
      !   ptc(mi,idx_nit) = max( zero, ptc(mi,idx_nit) )
      !   ptc(mi,idx_phos) = max( zero, ptc(mi,idx_phos) )
      !   ptc(mi,idx_sil) = max( zero, ptc(mi,idx_sil) )
      !   ptc(mi,idx_dia) = max( zero, ptc(mi,idx_dia) )
      !   ptc(mi,idx_flag) = max( zero, ptc(mi,idx_flag) )
      !   ptc(mi,idx_cyano) = max( zero, ptc(mi,idx_cyano) )
      !   ptc(mi,idx_mez) = max( zero, ptc(mi,idx_mez) )
      !   ptc(mi,idx_miz) = max( zero, ptc(mi,idx_miz) )
      !   ptc(mi,idx_det) = max( zero, ptc(mi,idx_det) )
      !   ptc(mi,idx_dets) = max( zero, ptc(mi,idx_dets) )
      !   ptc(mi,idx_ldon) = max( zero, ptc(mi,idx_ldon) )
      ! enddo
      mil = mcol(n)
      miu = mil+kb-2
      ptc(mil:miu,idx_amm) = max( zero, ptc(mil:miu,idx_amm) )
      ptc(mil:miu,idx_nit) = max( zero, ptc(mil:miu,idx_nit) )
      ptc(mil:miu,idx_phos) = max( zero, ptc(mil:miu,idx_phos) )
      ptc(mil:miu,idx_sil) = max( zero, ptc(mil:miu,idx_sil) )
      ptc(mil:miu,idx_dia) = max( zero, ptc(mil:miu,idx_dia) )
      ptc(mil:miu,idx_flag) = max( zero, ptc(mil:miu,idx_flag) )
      ptc(mil:miu,idx_cyano) = max( zero, ptc(mil:miu,idx_cyano) )
      ptc(mil:miu,idx_mez) = max( zero, ptc(mil:miu,idx_mez) )
      ptc(mil:miu,idx_miz) = max( zero, ptc(mil:miu,idx_miz) )
      ptc(mil:miu,idx_det) = max( zero, ptc(mil:miu,idx_det) )
      ptc(mil:miu,idx_dets) = max( zero, ptc(mil:miu,idx_dets) )
      ptc(mil:miu,idx_ldon) = max( zero, ptc(mil:miu,idx_ldon) )
    enddo

  end subroutine bio_limit

end module coupler

