module bioriver

  use exits,   only : exitme
  use io_subs, only : io_new_unit

  implicit none
  private
  
  public :: bioloadings, read_nrivseg, assign_rivermap

  private :: ReadLoadings
  
  integer(4), save, private ::  lun
  logical,    save, private :: fileopen

contains

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

  subroutine bioloadings ( ia, nlow, nup, nq, bndq, kmx, msrf, ptc )

    use dmi_mpi,        only : mpi_rank
    use dmi_mpi_global, only : iu06
    use cmod_params,    only : narea
    use cmod_arrays,    only : nqlist, nqfile, krq
    use constants,      only : zero

    implicit none

    integer(4), intent(in)    :: ia, nlow, nup, nq, kmx
    integer(4), intent(in)    :: msrf(0:,0:)
    real(8),    intent(in)    :: ptc(0:,1:)
    real(8),    intent(inout) :: bndq(:,:,0:)  

    integer(4),    save :: lun, ntc
    integer(4)          :: ios, iq, iqa, i, j, k
    character(1)        :: c1
    character(2)        :: c2
    character(4),  save :: ifmt
    character(11), save :: cfmt
    character(16)       :: qlogfile
    logical,       save :: firsttime = .true.


    if (firsttime) then
      firsttime = .false.
      lun = io_new_unit()
      write (qlogfile,'(a,i3.3)') 'loadings.out.',mpi_rank
      open(lun,file=qlogfile,status='unknown',iostat=ios)
      if (ios /= 0) call exitme (1,'input error - could open '//qlogfile)

      ntc = nup - nlow + 1
      if (ntc < 10) then
        ifmt = '(i1)'
        write(c1,ifmt) ntc
        cfmt = '(i3, '//c1//'f7.1)'
      elseif (ntc < 100) then
        ifmt = '(i2)'
        write(c2,ifmt) ntc
        cfmt = '(i3,'//c2//'f7.1)'
      else 
        call exitme (1,'dmi-ergom not prepared for so many tracers')
      endif
    endif

    ! FIXME: pls, come up with a better climatology, especially for 
    !        fall-back during operational runs, but until then:
    ! better safe than sorry:
    bndq(nlow:nup,:,0:) = zero
    do iq=1,nqfile(ia)
      i = krq(ia)%p(1,iq)
      j = krq(ia)%p(2,iq)
      if (i == -99 .or. j == -99) cycle
      do k = 1,kmx
        bndq(nlow:nup,k,iq)= ptc(msrf(i,j),1:ntc)
      enddo
    enddo

    ! read the data:
    if (nq > 0) then
      write(iu06,'(a32,i3)') 'Reading river loadings for area ',ia
      call ReadLoadings( ia, nqfile(ia), nlow, nup, bndq )

      ! shuffle the read source data:
      do iq=1,nqfile(ia)
        iqa = nqlist(ia)%p(iq)
        if (iqa > 0) bndq(nlow:nup,1,iqa) = bndq(nlow:nup,1,iq)
      enddo

      ! print loadings:
      write(lun,'(/,'' Loadings in area '',i3)') ia
      do iq=1,nqfile(ia)
        iqa = nqlist(ia)%p(iq)
        if (iqa > 0) write(lun,cfmt) iqa,bndq(nlow:nup,1,iqa)
      enddo
    else
      write(lun,'(/,'' No rivers in area '',i3)') ia
    endif

    ! clean up:
    if (ia == narea) close (lun)

  end subroutine bioloadings

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

  subroutine ReadLoadings( ia, nr, ncl, ncu, load )
    ! straight forward read of nutrient loadings from an ascii file

    use constants,      only : zeitstu
    use cmod_date,      only : jultimea
    use constants,      only : TIM, twenty
    use bioparam,       only : Nnorm_inv, NOR, Onorm_inv, &
!~ <tracers name=oxy; childOf=none>
                               !~ idx_BSH_ERGOM,   &
  !~ <children>
                               !~ idx_<childName>, &
  !~ </children>
!~ </tracers>
                               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_ship_N,   &
                               idx_nit_with_ship_N,   &
                               idx_dia_with_ship_N,   &
                               idx_flag_with_ship_N,   &
                               idx_cyano_with_ship_N,   &
                               idx_mez_with_ship_N,   &
                               idx_miz_with_ship_N,   &
                               idx_det_with_ship_N,   &
                               idx_ldon_with_ship_N,   &
                               idx_amm_with_river_N,   &
                               idx_nit_with_river_N,   &
                               idx_dia_with_river_N,   &
                               idx_flag_with_river_N,   &
                               idx_cyano_with_river_N,   &
                               idx_mez_with_river_N,   &
                               idx_miz_with_river_N,   &
                               idx_det_with_river_N,   &
                               idx_ldon_with_river_N,   &
                               ntbio, netcdf_rivers
    use bionetcdf_input,only : ReadNetcdfClimaloads
    use dmi_mpi_global, only : iu06

    implicit none

    integer(4), intent(in) :: ia, nr, ncl, ncu
    real(8), intent(inout) :: load(:,:,0:)    ! (1:ntracers,1:kmax,0:nr)

    character(1)  :: c1
    character(2)  :: c2
    character(3)  :: c3
    character(4)  :: c4
    character(4)  :: iformat
    character(14) :: cformat
    character(11) :: fname
    character(16) :: fnamec
    character(19) :: cdat
    integer(4)    :: lun, iq, ic, ios, nrc, rclbioloads,ii,day,d
    real(8)       :: timfl

    !- set up file name, format, etc ... ---------------------------------------
    nrc = nr*(ncu-ncl+1)
    if (nrc < 1) then
      return
    elseif (nrc < 10) then
      iformat = '(i1)'
      write(c1,iformat) nrc
      cformat = '(a19,   '//c1//'f8.1)'
    elseif (nrc < 100) then
      iformat = '(i2)'
      write(c2,iformat) nrc
      cformat = '(a19,  '//c2//'f8.1)'
    elseif (nrc < 1000) then
      iformat = '(i3)'
      write(c3,iformat) nrc
      cformat = '(a19, '//c3//'f8.1)'
    elseif (nrc < 10000) then
      iformat = '(i4)'
      write(c4,iformat) nrc
      cformat = '(a19,'//c4//'f8.1)'
    else
      ! trap
      call exitme(1,'ERROR: Number of bio sources too large.')
    endif
    if (ia < 10) then
      iformat = '(i1)'
      write(c1,iformat) ia
      fname   = 'bioloads_0'//c1
      fnamec  = 'bioclimaloads_0'//c1
    elseif (ia < 100) then
      iformat = '(i2)'
      write(c2,iformat) ia
      fname   = 'bioloads_'//c2
      fnamec  = 'bioclimaloads_'//c2
    else
      ! trap
      call exitme(1,'ERROR: Number of areas too large.')
    endif

    !- determine recordlength --------------------------------------------------
    rclbioloads = 1024
    ii = 19+nrc*8+1 
    do while (ii > rclbioloads)
      ! adjust record length of bioloads file
      rclbioloads = rclbioloads + 1024
    enddo

    !- open climatological river input and read it -----------------------------
    day = floor(zeitstu/24.0_8,4)+1
    if (netcdf_rivers) then
      call ReadNetcdfClimaloads( day, nr, iu06, ia, load )    
    else
      lun = io_new_unit()
      write(iu06,*) 'Reading clim river loads binary file: '//trim(fnamec)
      open (lun,file=fnamec,action='read',status='old',iostat=ios,             &
            recl=rclbioloads)
      if (ios /= 0) then
        write(iu06,*) 'Unable to open '//fnamec//', continuing without'  
      else
        do d = 1,day-1
          read (lun,'(1X)',iostat=ios)
          if (ios /= 0) exit
        enddo
        if (ios /= 0) then
          write(iu06,*) 'Unable to space in '//fnamec//', continuing without'
        else
          read(lun,cformat,iostat=ios) cdat,((load(ic,1,iq),ic=ncl,ncu),iq=1,nr)
          if (ios /= 0) then
            write(iu06,*) 'Unable to read '//fnamec//', continuing without'
          endif
        endif 
        ! load(ncl:ncl-2+ntbio,1,:) = load(ncl:ncl-2+ntbio,1,:) / Nnorm
        ! load(    ncl-1+ntbio,1,:) = load(    ncl-1+ntbio,1,:) / Onorm
          load(ncl-1+idx_amm            ,1,:) = load(ncl-1+idx_amm            ,1,:) * Nnorm_inv
          load(ncl-1+idx_nit            ,1,:) = load(ncl-1+idx_nit            ,1,:) * Nnorm_inv
          load(ncl-1+idx_phos           ,1,:) = load(ncl-1+idx_phos           ,1,:) * Nnorm_inv
          load(ncl-1+idx_sil            ,1,:) = load(ncl-1+idx_sil            ,1,:) * Nnorm_inv
          load(ncl-1+idx_dia            ,1,:) = load(ncl-1+idx_dia            ,1,:) * Nnorm_inv
          load(ncl-1+idx_flag           ,1,:) = load(ncl-1+idx_flag           ,1,:) * Nnorm_inv
          load(ncl-1+idx_cyano          ,1,:) = load(ncl-1+idx_cyano          ,1,:) * Nnorm_inv
          load(ncl-1+idx_mez            ,1,:) = load(ncl-1+idx_mez            ,1,:) * Nnorm_inv
          load(ncl-1+idx_miz            ,1,:) = load(ncl-1+idx_miz            ,1,:) * Nnorm_inv
          load(ncl-1+idx_det            ,1,:) = load(ncl-1+idx_det            ,1,:) * Nnorm_inv
          load(ncl-1+idx_dets           ,1,:) = load(ncl-1+idx_dets           ,1,:) * Nnorm_inv
          load(ncl-1+idx_ldon           ,1,:) = load(ncl-1+idx_ldon           ,1,:) * Nnorm_inv
          load(ncl-1+idx_amm_with_ship_N,1,:) = load(ncl-1+idx_amm_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_nit_with_ship_N,1,:) = load(ncl-1+idx_nit_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_dia_with_ship_N,1,:) = load(ncl-1+idx_dia_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_flag_with_ship_N,1,:) = load(ncl-1+idx_flag_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_cyano_with_ship_N,1,:) = load(ncl-1+idx_cyano_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_mez_with_ship_N,1,:) = load(ncl-1+idx_mez_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_miz_with_ship_N,1,:) = load(ncl-1+idx_miz_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_det_with_ship_N,1,:) = load(ncl-1+idx_det_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_ldon_with_ship_N,1,:) = load(ncl-1+idx_ldon_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_amm_with_river_N,1,:) = load(ncl-1+idx_amm_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_nit_with_river_N,1,:) = load(ncl-1+idx_nit_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_dia_with_river_N,1,:) = load(ncl-1+idx_dia_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_flag_with_river_N,1,:) = load(ncl-1+idx_flag_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_cyano_with_river_N,1,:) = load(ncl-1+idx_cyano_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_mez_with_river_N,1,:) = load(ncl-1+idx_mez_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_miz_with_river_N,1,:) = load(ncl-1+idx_miz_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_det_with_river_N,1,:) = load(ncl-1+idx_det_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_ldon_with_river_N,1,:) = load(ncl-1+idx_ldon_with_river_N,1,:) * Nnorm_inv
          load(    ncl-1+idx_oxy            ,1,:) = load(    ncl-1+idx_oxy,1,:) * Onorm_inv
          !~ load(ncl:ncl-1+ntbio              ,1,:) = load(ncl:ncl-1+ntbio,1,:) * Nnorm_inv
!~ <tracers name=oxy; childOf=none>
          !~ load(    ncl-1+idx_BSH_ERGOM,1,:) = load(    ncl-1+idx_BSH_ERGOM,1,:) * NOR
  !~ <children>
          !~ load(    ncl-1+idx_<childName>,1,:) = load(    ncl-1+idx_<trimChildName>,1,:) * NOR
  !~ </children>
!~ </tracers>
        close(lun)
      endif 
    endif

    !- open the file containing actual bioloadings and read it -----------------
    lun = io_new_unit()
    write(iu06,*) 'Reading river loads binary file: '//trim(fname)
    open (lun,file=fname,action='read',status='old',iostat=ios,recl=rclbioloads)
    if (ios /= 0) then
      write(iu06,*) 'Unable to open '//fname//', continuing without'
    else
      read (lun,'(a1)',iostat=ios) c1
      if (ios /= 0) then
        write(iu06,*) 'Unable to read '//fname//', continuing without'
      else
        read (lun,'(a19)',iostat=ios) cdat
        if (ios /= 0) then
          write(iu06,*) 'Unable to read date in '//fname//', continuing without'
        else
          call jultimea(cdat,timfl)
          if ( abs(TIM-timfl) <= twenty ) then
            backspace (lun)
            read (lun,cformat,iostat=ios) cdat,                                &
                                         ((load(ic,1,iq), ic=ncl,ncu), iq=1,nr)
            if (ios /= 0) then
              write(iu06,*) 'Unable to read loadings in '//fname//             &
                            ', continuing without'
            else
            ! load(ncl:ncl-2+ntbio,1,:) = load(ncl:ncl-2+ntbio,1,:) / Nnorm
            ! load(    ncl-1+ntbio,1,:) = load(    ncl-1+ntbio,1,:) / Onorm
          load(ncl-1+idx_amm            ,1,:) = load(ncl-1+idx_amm            ,1,:) * Nnorm_inv
          load(ncl-1+idx_nit            ,1,:) = load(ncl-1+idx_nit            ,1,:) * Nnorm_inv
          load(ncl-1+idx_phos           ,1,:) = load(ncl-1+idx_phos           ,1,:) * Nnorm_inv
          load(ncl-1+idx_sil            ,1,:) = load(ncl-1+idx_sil            ,1,:) * Nnorm_inv
          load(ncl-1+idx_dia            ,1,:) = load(ncl-1+idx_dia            ,1,:) * Nnorm_inv
          load(ncl-1+idx_flag           ,1,:) = load(ncl-1+idx_flag           ,1,:) * Nnorm_inv
          load(ncl-1+idx_cyano          ,1,:) = load(ncl-1+idx_cyano          ,1,:) * Nnorm_inv
          load(ncl-1+idx_mez            ,1,:) = load(ncl-1+idx_mez            ,1,:) * Nnorm_inv
          load(ncl-1+idx_miz            ,1,:) = load(ncl-1+idx_miz            ,1,:) * Nnorm_inv
          load(ncl-1+idx_det            ,1,:) = load(ncl-1+idx_det            ,1,:) * Nnorm_inv
          load(ncl-1+idx_dets           ,1,:) = load(ncl-1+idx_dets           ,1,:) * Nnorm_inv
          load(ncl-1+idx_ldon           ,1,:) = load(ncl-1+idx_ldon           ,1,:) * Nnorm_inv
          load(ncl-1+idx_amm_with_ship_N,1,:) = load(ncl-1+idx_amm_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_nit_with_ship_N,1,:) = load(ncl-1+idx_nit_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_dia_with_ship_N,1,:) = load(ncl-1+idx_dia_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_flag_with_ship_N,1,:) = load(ncl-1+idx_flag_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_cyano_with_ship_N,1,:) = load(ncl-1+idx_cyano_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_mez_with_ship_N,1,:) = load(ncl-1+idx_mez_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_miz_with_ship_N,1,:) = load(ncl-1+idx_miz_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_det_with_ship_N,1,:) = load(ncl-1+idx_det_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_ldon_with_ship_N,1,:) = load(ncl-1+idx_ldon_with_ship_N,1,:) * Nnorm_inv
          load(ncl-1+idx_amm_with_river_N,1,:) = load(ncl-1+idx_amm_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_nit_with_river_N,1,:) = load(ncl-1+idx_nit_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_dia_with_river_N,1,:) = load(ncl-1+idx_dia_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_flag_with_river_N,1,:) = load(ncl-1+idx_flag_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_cyano_with_river_N,1,:) = load(ncl-1+idx_cyano_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_mez_with_river_N,1,:) = load(ncl-1+idx_mez_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_miz_with_river_N,1,:) = load(ncl-1+idx_miz_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_det_with_river_N,1,:) = load(ncl-1+idx_det_with_river_N,1,:) * Nnorm_inv
          load(ncl-1+idx_ldon_with_river_N,1,:) = load(ncl-1+idx_ldon_with_river_N,1,:) * Nnorm_inv
          load(    ncl-1+idx_oxy            ,1,:) = load(    ncl-1+idx_oxy,1,:) * Onorm_inv
              !~ load(ncl:ncl-1+ntbio              ,1,:) = load(ncl:ncl-1+ntbio,1,:) * Nnorm_inv
!~ <tracers name=oxy; childOf=none>
              !~ load(    ncl-1+idx_BSH_ERGOM,1,:) = load(    ncl-1+idx_BSH_ERGOM,1,:) * NOR
  !~ <children>
              !~ load(    ncl-1+idx_<childName>,1,:) = load(    ncl-1+idx_<trimChildName>,1,:) * NOR
  !~ </children>
!~ </tracers>
            endif
          else 
            write(iu06,*) "Warning, no river loadings found for this day in ", &
                          fname
          endif
        endif
      endif
      close (lun)
    endif
    

  end subroutine ReadLoadings

!===============================================================================
!-----------------------------------------------------------------------------
 subroutine read_nrivseg(narea,nrivsegments) 
   ! read in how many special areas there are for each grid
    use exits,          only : exitme
    use io_subs,        only : io_new_unit
    use dmi_mpi_global, only : iu06
    use cmod_params,    only : maxareas 

    implicit none

    integer(4), intent(in) :: narea
    integer(4), intent(out):: nrivsegments(1:)                 
    integer(4) :: ios,i
    integer(4) :: nsegments(maxareas) 

    namelist /nsegnml /nsegments

  
    nsegments(:) = 0
    fileopen = .false.
    lun = io_new_unit()
    ios = 0 
    open (unit=lun, file='biorivermaps.nml', status='old', iostat=ios)
    
    if (ios /= 0) then
      !  apply default settings:
      write(iu06,'(a71)')                                                      &
      'Warning: File biorivermaps.nml not found. Continue with default settings'
    else
      fileopen = .true.
     
      read (unit=lun, nml=nsegnml, iostat=ios)
      
      if (ios /= 0) then
        write(iu06,*) 'Error: Cannot read basic numbers from biorivermaps.nml.'
        write(iu06,*) 'Error: The iostat failure was: ', ios
        call exitme(1)
      endif
      do i=1,narea
        if (nsegments(i) < 0) call exitme(1,'Invalid No. of segments')
        nrivsegments(i) = nsegments(i)
      enddo
    endif
    
 end subroutine read_nrivseg
 
 !-----------------------------------------------------------------------------

  subroutine assign_rivermap (nrivsegments,ia,narea,rivmap,msrf)
    
    use exits,          only : exitme
    use dmi_mpi_global, only : iu06

    implicit none

    integer(4), intent(in)      :: nrivsegments,ia,narea, msrf(0:,0:)
    integer(4), intent(inout)   :: rivmap(0:)
    integer(4) :: is, ios, ijsegment(1:4), nsurf, i, j
    

    namelist /ijsegnml / ijsegment
    
    
    if(fileopen)then
      do is=1,nrivsegments
        read (unit=lun, nml=ijsegnml, iostat=ios)
        if (ios /= 0) then
          write(iu06,*) 'Error: Cannot read ij-segment from biorivermaps.nml.'
          write(iu06,*) 'Error: The iostat failure was: ', ios
          call exitme(1)
        endif    
       
        !set index of riverarea to zero in rivermap index list
        do i = ijsegment(1),ijsegment(2)
          do j = ijsegment(3),ijsegment(4)
            nsurf= msrf(i,j)
            rivmap(nsurf) = 0
          enddo
        enddo
        
      enddo

      if (ia == narea) close(lun)
    endif
  end subroutine assign_rivermap


end module bioriver
