module bioinput

  use constants,  only : zero, half, one

  implicit none
  private

  public  :: bioread, biobndcl, bioread_npr, bioatmdep

  real(8), parameter, private :: hpm = 720.0_8
  real(8), parameter, private :: onemill = 1000000.0_8


contains

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

  subroutine bioread( iact, n3d, n2d, mmx, nmx, kmx, nz, msrf, mcol, kh, narea,&
                      iu06, benthos, ntben, ptc, t, s, bndz )

    use io_subs,         only : io_new_unit
    use exits,           only : exitme
    use cmod_params,     only : mainarea
    use constants,       only : aufdat
    use bioparam,        only : Nnorm_inv, Onorm_inv, NOR,    &
                                idx_amm            ,          &
                                idx_nit            ,          &
                                idx_phos           ,          &
                                idx_sil            ,          &
                                idx_dia            ,          &
                                idx_flag           ,          &
                                idx_cyano          ,          &
                                idx_mez            ,          &
                                idx_miz            ,          &
                                idx_det            ,          &
                                idx_dets           ,          &
                                idx_ldon           ,          &
                                idx_oxy            ,          &
                                idx_amm_with_totatmos_N,          &
                                idx_nit_with_totatmos_N,          &
                                idx_dia_with_totatmos_N,          &
                                idx_flag_with_totatmos_N,          &
                                idx_cyano_with_totatmos_N,          &
                                idx_mez_with_totatmos_N,          &
                                idx_miz_with_totatmos_N,          &
                                idx_det_with_totatmos_N,          &
                                idx_ldon_with_totatmos_N,          &
                                idx_amm_with_agri_N,          &
                                idx_nit_with_agri_N,          &
                                idx_dia_with_agri_N,          &
                                idx_flag_with_agri_N,          &
                                idx_cyano_with_agri_N,          &
                                idx_mez_with_agri_N,          &
                                idx_miz_with_agri_N,          &
                                idx_det_with_agri_N,          &
                                idx_ldon_with_agri_N,          &
                                ntbio_notag, ntben_notag,     &
                                ntbio, ibiofields, biocoldst, &
                                read_binary_restart_notag
    use biogas,          only : oxygensat
    use biopermute,      only : bio_permute
    use bionetcdf_input, only : netcdf_init

    implicit none

    integer(4), intent(in)    :: iact, narea, iu06, n3d, n2d, ntben,           &
                                 mmx, nmx, kmx, nz
    integer(4), intent(in)    :: msrf(0:,0:), mcol(0:), kh(0:)
    real(8),    intent(in)    :: t(0:), s(0:)
    real(8),    intent(out)   :: benthos(0:,1:) 
    real(8),    intent(inout) :: ptc(0:,1:), bndz(:,:,0:)

    integer(4)       :: ios, i
    integer(4), save :: lun, kmx_main, nz_main
    character(11)    :: starta
    character(60)    :: fnam
    character(3)     :: c3


    if (.not.biocoldst) then
      ! read from a bio_start file (i.e. like a hotstarted simulation)

      if (iact == 1) then
        lun = io_new_unit()
        fnam='bio_initial.bin'
        write(iu06,*) 'Reading bio restart binary file: '//trim(fnam)
        open(lun,file=fnam,form='unformatted',iostat=ios,         &
             action='read',asynchronous='yes',access='stream',status='old')
        if (ios /= 0) call exitme (1, 'Could not open bio_start' )
        read(lun,iostat=ios) starta
        if (ios /= 0) then
          write(iu06,*) 'I/O problem - iostat is:',ios
          call exitme (2, 'Could not read starta from bio_start' )
        endif
        if (starta(2:11) /= aufdat(2:11)) then
          call exitme(3,'Time stamp in bio_start does not match starting time')
        endif
      endif
      
      ! NEUMANND: TODO, REVERT?
      ! read(lun,iostat=ios) ptc(:,1:ntbio), benthos(:,:)
      if (read_binary_restart_notag .and. (ntbio > ntbio_notag .or. ntben > ntben_notag)) then
        read(lun,iostat=ios) ptc(:,1:ntbio_notag), benthos(:,1:ntben_notag)
        if (ntbio > ntbio_notag) then
          ptc(:,ntbio_notag+1:ntbio) = zero
        endif
        if (ntben > ntben_notag) then
          benthos(:,ntben_notag+1:ntben) = zero
        endif
      else
        read(lun,iostat=ios) ptc(:,1:ntbio), benthos(:,1:ntben)
      endif
      ! NEUMANND: END TODO
      if (ios /= 0) then
        write(iu06,*) 'I/O problem - iostat is:',ios
        write(iu06,*) '                   area:',iact
        call exitme (4, 'Could not read bio_start' )
      endif
      ! NEUMANND: special treatment of oxygen
!      ptc(:,1:ntbio-1) = ptc(:,1:ntbio-1)/Nnorm
!      ptc(:,ntbio)     = ptc(:,ntbio)/Onorm
      ptc(:,:                  ) = ptc(:,:)*Nnorm_inv
      ptc(:,idx_oxy            )     = ptc(:,idx_oxy            )*NOR ! instead of *Nnorm/Onorm
      benthos(:,:              ) = benthos(:,:)*Nnorm_inv

      if (iact == mainarea) then
        kmx_main = kmx
        nz_main  = nz
      endif

      if (iact == narea) then
        if (nz_main >= 1) then
          read(lun,iostat=ios) bndz(idx_amm+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_amm
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_nit+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_nit
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_phos+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_phos
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_sil+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_sil
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_dia+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_dia
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_flag+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_flag
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_cyano+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_cyano
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_mez+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_mez
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_miz+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_miz
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_det+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_det
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_dets+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_dets
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_ldon+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_ldon
            call exitme (4, 'Could not read bio sponges' )
          endif
          read(lun,iostat=ios) bndz(idx_oxy+2,1:kmx_main,0:)
          if (ios /= 0) then
            write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_oxy
            call exitme (4, 'Could not read bio sponges' )
          endif
          if (read_binary_restart_notag) then
            bndz(ntbio_notag+3:ntbio+2,1:kmx_main,0:) = zero
          else
            read(lun,iostat=ios) bndz(idx_amm_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_amm_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_nit_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_nit_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_dia_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_dia_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_flag_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_flag_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_cyano_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_cyano_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_mez_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_mez_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_miz_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_miz_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_det_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_det_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_ldon_with_totatmos_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_ldon_with_totatmos_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_amm_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_amm_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_nit_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_nit_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_dia_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_dia_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_flag_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_flag_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_cyano_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_cyano_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_mez_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_mez_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_miz_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_miz_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_det_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_det_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
            read(lun,iostat=ios) bndz(idx_ldon_with_agri_N+2,1:kmx_main,0:)
            if (ios /= 0) then
              write(iu06,*) 'I/O problem - ia, iostat and tracer index are:',mainarea,ios,idx_ldon_with_agri_N
              call exitme (4, 'Could not read bio sponges' )
            endif
          endif
        endif
      endif

      if (iact == narea) close(lun)

    else
      ! assign bio vars as constants or read from a file (like coldstarted)

!      ptc(0,1:ntbio-1)  = onemill/Nnorm
!      ptc(0,ntbio)      = onemill/Onorm
!      ptc(1:,1:ntbio-1) = one
      ! TODO, NEUMANND how to deal with colored tracers?
      ptc(0,idx_amm            )  = onemill*Nnorm_inv
      ptc(0,idx_nit            )  = onemill*Nnorm_inv
      ptc(0,idx_phos           )  = onemill*Nnorm_inv
      ptc(0,idx_sil            )  = onemill*Nnorm_inv
      ptc(0,idx_dia            )  = onemill*Nnorm_inv
      ptc(0,idx_flag           )  = onemill*Nnorm_inv
      ptc(0,idx_cyano          )  = onemill*Nnorm_inv
      ptc(0,idx_mez            )  = onemill*Nnorm_inv
      ptc(0,idx_miz            )  = onemill*Nnorm_inv
      ptc(0,idx_det            )  = onemill*Nnorm_inv
      ptc(0,idx_dets           )  = onemill*Nnorm_inv
      ptc(0,idx_ldon           )  = onemill*Nnorm_inv
      ptc(0,idx_oxy            )  = onemill*Onorm_inv
      ptc(0,idx_amm_with_totatmos_N)  = zero
      ptc(0,idx_nit_with_totatmos_N)  = zero
      ptc(0,idx_dia_with_totatmos_N)  = zero
      ptc(0,idx_flag_with_totatmos_N)  = zero
      ptc(0,idx_cyano_with_totatmos_N)  = zero
      ptc(0,idx_mez_with_totatmos_N)  = zero
      ptc(0,idx_miz_with_totatmos_N)  = zero
      ptc(0,idx_det_with_totatmos_N)  = zero
      ptc(0,idx_ldon_with_totatmos_N)  = zero
      ptc(0,idx_amm_with_agri_N)  = zero
      ptc(0,idx_nit_with_agri_N)  = zero
      ptc(0,idx_dia_with_agri_N)  = zero
      ptc(0,idx_flag_with_agri_N)  = zero
      ptc(0,idx_cyano_with_agri_N)  = zero
      ptc(0,idx_mez_with_agri_N)  = zero
      ptc(0,idx_miz_with_agri_N)  = zero
      ptc(0,idx_det_with_agri_N)  = zero
      ptc(0,idx_ldon_with_agri_N)  = zero

      ptc(1:,idx_amm            )  = one
      ptc(1:,idx_amm_with_totatmos_N)  = zero
      ptc(1:,idx_amm_with_agri_N)  = zero
      ptc(1:,idx_nit            )  = one
      ptc(1:,idx_nit_with_totatmos_N)  = zero
      ptc(1:,idx_nit_with_agri_N)  = zero
      ptc(1:,idx_phos           )  = one
      ptc(1:,idx_sil            )  = one
      ptc(1:,idx_dia            )  = one
      ptc(1:,idx_dia_with_totatmos_N)  = zero
      ptc(1:,idx_dia_with_agri_N)  = zero
      ptc(1:,idx_flag           )  = one
      ptc(1:,idx_flag_with_totatmos_N)  = zero
      ptc(1:,idx_flag_with_agri_N)  = zero
      ptc(1:,idx_cyano          )  = one
      ptc(1:,idx_cyano_with_totatmos_N)  = zero
      ptc(1:,idx_cyano_with_agri_N)  = zero
      ptc(1:,idx_mez            )  = one
      ptc(1:,idx_mez_with_totatmos_N)  = zero
      ptc(1:,idx_mez_with_agri_N)  = zero
      ptc(1:,idx_miz            )  = one
      ptc(1:,idx_miz_with_totatmos_N)  = zero
      ptc(1:,idx_miz_with_agri_N)  = zero
      ptc(1:,idx_det            )  = one
      ptc(1:,idx_det_with_totatmos_N)  = zero
      ptc(1:,idx_det_with_agri_N)  = zero
      ptc(1:,idx_dets           )  = one
      ptc(1:,idx_ldon           )  = one
      ptc(1:,idx_ldon_with_totatmos_N)  = zero
      ptc(1:,idx_ldon_with_agri_N)  = zero
      ptc(1:,idx_oxy            )  = one

      ! Initialize oxygen at 100% saturation
      do i=1,n3d
        ptc(i,idx_oxy            ) = oxygensat(t(i),s(i))*Onorm_inv
      enddo

      benthos(:,:               )  = zero

      if (ibiofields == 1) then
        !- read from a dummy restart file containing only biological tracers:
        if (iact == 1) then
          lun = io_new_unit()
          fnam='restart.bio'
          write(iu06,*) 'Reading bio coldstart binary file: '//trim(fnam)
          open(lun,file=fnam,form='unformatted',iostat=ios,           &
               action='read',status='old')
          if (ios /= 0) call exitme (5, 'Could not open restart.bio' )
        endif

        read(lun,iostat=ios) ptc(:,1:ntbio), benthos(:,:)
        if (ios /= 0) then
          write(iu06,*) 'I/O problem - iostat is:',ios
          write(iu06,*) '                   area:',iact
          call exitme (6, 'Could not read restart.bio' )
        endif
        if (iact == narea) close(lun)

      elseif (ibiofields == 2) then
        !- read biological tracers (excl. benthos) from an ascii data file:

        if (iact < 10) then
          write(c3,'(i1)') iact
        elseif (iact < 100) then
          write(c3,'(i2)') iact
        elseif (iact < 1000) then
          write(c3,'(i3)') iact
        else
          call exitme(7,'Cannot handle more than 999 areas.')
        endif

        fnam = 'bio.init'//trim(c3)
        lun = io_new_unit()
        write(iu06,*) 'Reading biological fields from ',fnam
        open(lun,file=fnam,form="formatted",status="old",iostat=ios)
        if (ios /= 0) call exitme(8,'Could not open '//trim(fnam))

        ptc(0,1:ntbio) = onemill
        do i=1,n3d
          read(lun,*) ptc(i,1:ntbio)
        enddo
        ptc(:,1:ntbio) = ptc(:,1:ntbio)*Nnorm_inv
        ptc(:,idx_oxy)     = ptc(:,idx_oxy)*NOR ! instead of *Nnorm/Onorm

        if (iact == narea) close(lun)

      endif
    endif

    !- permute the fields
    call bio_permute(mmx, nmx, kmx, n2d, n3d, msrf, mcol, kh, ptc, benthos,    &
                     1, ntbio, ntben)

    if (biocoldst .and. ibiofields == 3) then
      ! Read biological tracers from netcdf files:
      call netcdf_init (iact, mmx, nmx, kmx, msrf, mcol, kh, ntbio, ntben,     &
                        ptc, benthos, iu06)
    endif

  end subroutine bioread

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

  subroutine biobndcl(kmx, kh, msrf, mcol, nz, nlow, nup, bndzk, iu06, ptc)

    use bionetcdf_input, only : ReadNetcdfBounds
    use bioparam,        only : netcdf_bounds, read_binary_rwzk_notag, ntbio_notag
    use cmod_miscvars,   only : jstuna
    use cmod_params,     only : krz, mainarea
    use exits,           only : exitme
    use io_subs,         only : io_new_unit
    implicit none

    integer(4), intent(in)    :: kmx, kh(0:), nz, nlow, nup, iu06
    integer(4), intent(in)    :: msrf(0:,0:), mcol(0:)
    real(8),    intent(in)    :: ptc(0:,1:)
    real(8),    intent(inout) :: bndzk(:,:,0:)        ! (ntracers,kmx,0:nz)

    character(20)  :: fnam
    character(256) :: cformat
    integer(4)     :: lune, nbp, k, t, mon, monat, kb, n, ios, ibnd, jbnd, ntr
                      ! NOTE: monat and nbp are only read in but never used
    integer(4)     :: msf, mi
    integer(4)     :: khbnd(nz)
    real(8)        :: pqr, w1, w2
    real(8), dimension(nup-nlow+1,kmx,0:nz,12) :: rwzke_a

    !- initialise --------------------------------------------------------------
    ntr = nup-nlow+1
    rwzke_a(:,:,:,:) = zero
    khbnd(:) = 0
    do n=1,nz
      ibnd = krz(mainarea)%p(1,n)
      jbnd = krz(mainarea)%p(2,n)
      msf = msrf(ibnd,jbnd)
      khbnd(n) = kh(msf)
      do mon = 1,12
        rwzke_a(1:ntr,1,n,mon) = ptc(msf,1:ntr)
      enddo
      mi = msf
      do k=2,kmx
        if (k <= khbnd(n)) mi = mcol(msf) - 2 + k
        do mon = 1,12
          rwzke_a(1:ntr,k,n,mon) = ptc(mi,1:ntr)
        enddo
      enddo
    enddo

    !- read climatology ecological or passive tracers --------------------------
    if (nz > 0) then
      if (netcdf_bounds) then
        call ReadNetcdfBounds (nz, iu06, ntr, kmx, rwzke_a)
      else
        fnam = 'rwzkebnd'
        write(iu06,*) 'Reading bio boundary binary file: '//trim(fnam)
        lune = io_new_unit()
        open (lune,file=fnam,status='old',iostat=ios,pad='NO')
        if (ios /= 0) then
          write(iu06,*) 'Cannot open rwzkebnd ... continue without'
        else
          if     (kmx < 10) then
            write(cformat,'(i1)') kmx
          elseif (kmx < 100) then
            write(cformat,'(i2)') kmx
          elseif (kmx < 1000) then
            write(cformat,'(i3)') kmx
          else
            call exitme(2,'Cannot handle > 999 layers for tracer climatology.')
          endif

          cformat = '(i4,'//trim(cformat)//'f8.3)'

          if (read_binary_rwzk_notag .and. ntr > ntbio_notag) then
            do t=1,ntbio_notag
              do n=1,nz
                read (lune,*,iostat=ios)
                if (ios /= 0) call exitme (2,'Cannot read rwzkebnd')
                read (lune,'(i6)',iostat=ios) nbp
                if (ios /= 0) call exitme (3,'Cannot read rwzkebnd')
                kb = khbnd(n)
                do mon=1,12
                  read (lune,cformat,iostat=ios) monat,(rwzke_a(t,k,n,mon),k=1,kb)
                  if (ios /= 0) call exitme (4,'Cannot read rwzkebnd')
                enddo
                read (lune,*,iostat=ios)
                if (ios /= 0) call exitme (5,'Cannot read rwzkebnd')
              enddo
            enddo
            ! set remaining values to zero
            rwzke_a(ntbio_notag+1:ntr,:,:,:) = zero
          else
            do t=1,ntr
              do n=1,nz
                read (lune,*,iostat=ios)
                if (ios /= 0) call exitme (2,'Cannot read rwzkebnd')
                read (lune,'(i6)',iostat=ios) nbp
                if (ios /= 0) call exitme (3,'Cannot read rwzkebnd')
                kb = khbnd(n)
                do mon=1,12
                  read (lune,cformat,iostat=ios) monat,(rwzke_a(t,k,n,mon),k=1,kb)
                  if (ios /= 0) call exitme (4,'Cannot read rwzkebnd')
                enddo
                read (lune,*,iostat=ios)
                if (ios /= 0) call exitme (5,'Cannot read rwzkebnd')
              enddo
            enddo
          endif
          close(lune)
        endif
      endif

      !- interpolation to actual day:
      mon = (jstuna+720/2)/720       !       mon=0,1,2,...,11,12
      pqr = real(jstuna-mon*720,8)/hpm
      w1  = half - pqr
      w2  = half + pqr

      do n=1,nz
        kb = khbnd(n)
        do k=1,kb
          bndzk(nlow:nup,k,n) = w1*rwzke_a(1:ntr,k,n,mod(mon+11,12)+1)         &
                              + w2*rwzke_a(1:ntr,k,n,mod(mon+12,12)+1)
        enddo
        do k=kb+1,kmx
          bndzk(nlow:nup,k,n) = bndzk(nlow:nup,kb,n)
        enddo
      enddo
      bndzk(nlow:nup,:,0) = bndzk(nlow:nup,:,1)


    !- Dummy values if no bnd points -------------------------------------------
    else
      bndzk(:,:,:) = zero
    endif

  end subroutine biobndcl

  !=============================================================================
  
  subroutine bioread_npr( iact, narea, iu06, npr, mi, ni, n2d, msrf )
    use io_subs,   only : io_new_unit
    use exits,     only : exitme
    use bioparam,  only : rfr, npread
    use biopermute,only : npr_task_perm

    implicit none

    integer(4), intent(in)    :: iact, narea, iu06, mi, ni, n2d
    integer(4), intent(in)    :: msrf(0:,0:)
    real(8),    intent(out)   :: npr(0:)

    integer(4)       :: ios
    logical, save    :: fopened = .false.
    integer(4), save :: lun
    character(20)    :: fnam

    if (iact == 1 .and. npread) then
      fopened = .true.
      lun = io_new_unit()
      fnam='NPr.bin'
      write(iu06,*) 'Reading bio NP-ratio binary file: '//trim(fnam)
      open(lun,file=fnam,form='unformatted',iostat=ios,action='read',     &
           access='stream',status='old')
      if (ios /= 0) then
        write(iu06,*) 'WARNING: Cannot open '//trim(fnam)//', using a fixed NP-ratio'
        fopened = .false.
      endif
    endif

    if(fopened) then
      read(lun,iostat=ios) npr(0:)
      if (ios /= 0) then
        write(iu06,*) 'ERROR: cannot read '//trim(fnam)//' ... stop'
        call exitme(1)
      endif
      call npr_task_perm(mi,ni,n2d,msrf,npr)
    else
      npr(:)=rfr
    endif

    if (iact == narea .and. fopened) close(lun)

  end subroutine bioread_npr

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

  subroutine bioatmdep(iu06, ia, nmx, mmx,    &
                       srfammflx,             &
                       srfnitflx,             &
                       srfphosflx,             &
                       srfamm_with_totatmos_Nflx,             &
                       srfnit_with_totatmos_Nflx,             &
                       srfamm_with_agri_Nflx,             &
                       srfnit_with_agri_Nflx,             &
                       msrf)
   
    use bionetcdf_input, only: ReadNETCDFatm_dep 
    use bioparam,        only: netcdf_atm,               &
                               srfnredflx_const,          &
                               srfnoxflx_const,          &
                               srfpoxflx_const,          &
                               adread, tag_atmos_dep
    use io_subs,         only: io_new_unit
    use exits,           only: exitme
    implicit none
    
    integer(4), intent(in)  :: iu06, ia, nmx, mmx, msrf(0:,0:) 
    real(8),    intent(out) :: srfammflx(0:) ! atmospheric deposition flux: reduced nitrogen
    real(8),    intent(out) :: srfnitflx(0:) ! atmospheric deposition flux: oxidized nitrogen
    real(8),    intent(out) :: srfphosflx(0:) ! atmospheric deposition flux: oxidized phosphorus (phosphate)
    real(8),    intent(out) :: srfamm_with_totatmos_Nflx(0:) ! atmospheric deposition flux: reduced nitrogen; containing 
    real(8),    intent(out) :: srfnit_with_totatmos_Nflx(0:) ! atmospheric deposition flux: oxidized nitrogen; containing 
    real(8),    intent(out) :: srfamm_with_agri_Nflx(0:) ! atmospheric deposition flux: reduced nitrogen; containing 
    real(8),    intent(out) :: srfnit_with_agri_Nflx(0:) ! atmospheric deposition flux: oxidized nitrogen; containing 

    character(3)           :: c3
    character(13)          :: fnameamm
    integer(4), save       :: lunamm
    character(13)          :: fnamenit
    integer(4), save       :: lunnit
    character(13)          :: fnamephos
    integer(4), save       :: lunphos
    character(13)          :: fnameamm_with_totatmos_N
    integer(4), save       :: lunamm_with_totatmos_N
    character(13)          :: fnamenit_with_totatmos_N
    integer(4), save       :: lunnit_with_totatmos_N
    character(13)          :: fnameamm_with_agri_N
    integer(4), save       :: lunamm_with_agri_N
    character(13)          :: fnamenit_with_agri_N
    integer(4), save       :: lunnit_with_agri_N
    integer(4)             :: ios, i,j,idx
    real(8)                :: srfflx_in(mmx,nmx)
    real(8), parameter     :: missing = -5.0_8

    if (.not. adread) then
      srfammflx(0:) =  srfnredflx_const
      srfnitflx(0:) =  srfnoxflx_const
      srfphosflx(0:) =  srfpoxflx_const
      srfamm_with_totatmos_Nflx(0:) =  zero
      srfnit_with_totatmos_Nflx(0:) =  zero
      srfamm_with_agri_Nflx(0:) =  zero
      srfnit_with_agri_Nflx(0:) =  zero
      return
    endif    

    if (netcdf_atm) then
      call ReadNETCDFatm_dep(iu06,ia,                         &
                             srfnredflx_const,                &
                             srfnoxflx_const,                &
                             srfpoxflx_const,                &
                             nmx,mmx,                         &
                             srfammflx,                      &
                             srfnitflx,                      &
                             srfphosflx,                      &
                             srfamm_with_totatmos_Nflx,                      &
                             srfnit_with_totatmos_Nflx,                      &
                             srfamm_with_agri_Nflx,                      &
                             srfnit_with_agri_Nflx,                      &
                             msrf)
    else
      if (ia < 10) then
        write(c3,'(i1)') ia
        fnameamm='atm_nred_0'//trim(c3)
        fnamenit='atm_nox_0'//trim(c3)
        fnamephos='atm_pox_0'//trim(c3)
        fnameamm_with_totatmos_N='atm_nred_with_totatmos_N_0'//trim(c3)
        fnamenit_with_totatmos_N='atm_nox_with_totatmos_N_0'//trim(c3)
        fnameamm_with_agri_N='atm_nred_with_agri_N_0'//trim(c3)
        fnamenit_with_agri_N='atm_nox_with_agri_N_0'//trim(c3)
      elseif (ia < 100) then
        write(c3,'(i2)') ia
        fnameamm='atm_nred_'//trim(c3)
        fnamenit='atm_nox_'//trim(c3)
        fnamephos='atm_pox_'//trim(c3)
        fnameamm_with_totatmos_N='atm_nred_with_totatmos_N_'//trim(c3)
        fnamenit_with_totatmos_N='atm_nox_with_totatmos_N_'//trim(c3)
        fnameamm_with_agri_N='atm_nred_with_agri_N_'//trim(c3)
        fnamenit_with_agri_N='atm_nox_with_agri_N_'//trim(c3)
      else
        call exitme(7,'Cannot handle more than 99 areas.')
      endif

      ! NOTE: We need to distinguish between untagged/uncolored deposition
      !       and tagged/colored deposition, i.e. from shipping. The latter
      !       is set to 0 (zero) if no data are read in from files and the
      !       former is set to srf*flx_const if no data are read in from
      !       files.
      ! untagged/bulk atmospheric deposition (reduce and oxidized nitrogen; phosphate)
      lunamm=io_new_unit()
      write(iu06,*) 'Reading atmos dep binary file: '//trim(fnameamm)
      open(lunamm, file = trim(fnameamm),action='read',status='old',iostat=ios)
      if (ios /= 0) then
        write(iu06,*) 'No file: Running without variable nred Deposition'
        srfammflx(0:) =  srfnredflx_const
      else
        read(lunamm,'(10f8.4)',iostat=ios) ((srfflx_in(i,j), i=1,mmx), j=1,nmx)
        if (ios /= 0) then
          write(iu06,*) 'Cannot read nred Deposition'
          srfammflx(0:) =  srfnredflx_const
        else
          srfammflx(0) =  srfnredflx_const
          do j=1,nmx
            do i=1,mmx
              idx = 0
              if (msrf(i,j) > 0) then
                idx = msrf(i,j)
                if (srfflx_in(i,j) <= missing) then
                  write(iu06,'(a,a,a,3i5)') 'Missing value in ',fnameamm,' at ', &
                                            i,j
                  srfammflx(idx) =  srfnredflx_const
                else
                  srfammflx(idx) = srfflx_in(i,j)
                endif
              endif
            enddo
          enddo
        endif
        close(lunamm)
      endif
      
      lunnit=io_new_unit()
      write(iu06,*) 'Reading atmos dep binary file: '//trim(fnamenit)
      open(lunnit, file = trim(fnamenit),action='read',status='old',iostat=ios)
      if (ios /= 0) then
        write(iu06,*) 'No file: Running without variable nox Deposition'
        srfnitflx(0:) =  srfnoxflx_const
      else
        read(lunnit,'(10f8.4)',iostat=ios) ((srfflx_in(i,j), i=1,mmx), j=1,nmx)
        if (ios /= 0) then
          write(iu06,*) 'Cannot read nox Deposition'
          srfnitflx(0:) =  srfnoxflx_const
        else
          srfnitflx(0) =  srfnoxflx_const
          do j=1,nmx
            do i=1,mmx
              idx = 0
              if (msrf(i,j) > 0) then
                idx = msrf(i,j)
                if (srfflx_in(i,j) <= missing) then
                  write(iu06,'(a,a,a,3i5)') 'Missing value in ',fnamenit,' at ', &
                                            i,j
                  srfnitflx(idx) =  srfnoxflx_const
                else
                  srfnitflx(idx) = srfflx_in(i,j)
                endif
              endif
            enddo
          enddo
        endif
        close(lunnit)
      endif
      
      lunphos=io_new_unit()
      write(iu06,*) 'Reading atmos dep binary file: '//trim(fnamephos)
      open(lunphos, file = trim(fnamephos),action='read',status='old',iostat=ios)
      if (ios /= 0) then
        write(iu06,*) 'No file: Running without variable pox Deposition'
        srfphosflx(0:) =  srfpoxflx_const
      else
        read(lunphos,'(10f8.4)',iostat=ios) ((srfflx_in(i,j), i=1,mmx), j=1,nmx)
        if (ios /= 0) then
          write(iu06,*) 'Cannot read pox Deposition'
          srfphosflx(0:) =  srfpoxflx_const
        else
          srfphosflx(0) =  srfpoxflx_const
          do j=1,nmx
            do i=1,mmx
              idx = 0
              if (msrf(i,j) > 0) then
                idx = msrf(i,j)
                if (srfflx_in(i,j) <= missing) then
                  write(iu06,'(a,a,a,3i5)') 'Missing value in ',fnamephos,' at ', &
                                            i,j
                  srfphosflx(idx) =  srfpoxflx_const
                else
                  srfphosflx(idx) = srfflx_in(i,j)
                endif
              endif
            enddo
          enddo
        endif
        close(lunphos)
      endif
      

      ! atmospheric nitrogen deposition from shipping
        ! For ERGOM CGT users: One could also use the "celements" + "containingstracers"
        !                      loops instead of the "tracers" loop. However, one needed
        !                      to replace the tags. E.g. <tracerBelow> in the tracers
        !                      loop is <ctBelow> in the containingtracers loop. Since
        !                      the code written above for uncolored tracers is quite 
        !                      similar (note the different standard values assigned), it
        !                      is reasonable to choose the loop type for which the least
        !                      changes are needed.
        ! Alternatives:
        !   <  celements color=ship; atmosdep=1>
        !     <  containingtracers vertLoc=SUR>
        !
        ! START: If you see no code below, the code was generated without "N from ship"
      if ( tag_atmos_dep ) then
        ! END: If you see no code above, the code was generated without "N from ship"
        ! 
        ! Alternatives:
        !   <  /containingtracers>
        !     <  /celements>
      else
      endif
        
    endif

  end subroutine bioatmdep

end module bioinput
