module biooutput
  !-----------------------------------------------------------------------------
  ! stuff for output of ecomodel vars
  !-----------------------------------------------------------------------------
  ! except for the restart file, all the output is written in netcdf

implicit none 
  private

  integer(4),              private, save :: ibncdf
  integer(4), allocatable, private, save :: ncid(:) 
  integer(4), allocatable, private, save :: varid(:,:)
  integer(4), allocatable, private, save :: timvid(:)
  integer(4), allocatable, private, save :: sedvarid(:,:)
  integer(4), allocatable, private, save :: sedtimvid(:)
  integer(4), allocatable, private, save :: chlvarid(:)
  integer(4), allocatable, private, save :: chltimvid(:)
  integer(4), allocatable, private, save :: sedncid(:)
  integer(4), allocatable, private, save :: chlncid(:)
  integer(4), allocatable, private, save :: lightvarid(:)
  integer(4), allocatable, private, save :: lightncid(:)
  integer(4), allocatable, private, save :: lighttimvid(:)
  integer(4), allocatable, private, save :: secchivarid(:)
  integer(4), allocatable, private, save :: secchincid(:)
  integer(4), allocatable, private, save :: secchitimvid(:)

  logical, parameter :: bioncdf = .true. 

  public  :: biofile_new, biofile_next
  public  :: write_bio_restart, val_benthos
  public  :: biofile_close, biofile_write 
  public  :: diag_out, bioncdf
  private :: get_standard_name, get_standard_name_sed
  private :: get_var_units, get_var_units_sed, biofile_write_nc_diag
  private :: biofile_write_nc_diag_2d

contains

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

  subroutine biofile_new ( )

    !- Open files to save the pelagic state variables, the benthic state 
    !- variables and diagnostic chlorophyll
    use netcdf_interface, only : ncdf_create_new,ncdf_set_dimensions,         &
                                 ncdf_set_dim_vars,                           &
                                 ncdf_set_3D_variables, ncdf_set_4D_variables,&
                                 ncdf_set_dim_attr, ncdf_end_definitions,     &
                                 ncdf_set_bio_var_attrib,ncdf_save_dim_vars,  &
                                 ncdf_set_sed_var_attrib,                     &
                                 ncdf_set_chl_var_attrib,                     &
                                 ncdf_set_secchi_var_attrib,                  &
                                 ncdf_set_light_var_attrib
    use exits,            only : exitme
    use io_subs,          only : io_new_unit
    use cmod_arrays,      only : dz, xkoor, ykoor
    use writenetcdf_BSH,  only : arc2deg, calculate_depth
    use constants,        only : archive
    use cmod_params,      only : iw2, narea,mmx,nmx,kmx, cfgfile
    use bioparam,         only : ntbio,ntben, do_diag_chl, write_light,       &
                                 do_diag_secchi
   
    implicit none

     
    integer(4)    :: lonid, latid, dpthid, timid
    integer(4)    :: lonid_sed, latid_sed, dpthid_sed, timid_sed
    integer(4)    :: lonid_chl, latid_chl, dpthid_chl, timid_chl
    integer(4)    :: lonid_secchi, latid_secchi, dpthid_secchi, timid_secchi
    integer(4)    :: lonid_light, latid_light, dpthid_light, timid_light
    integer(4)    :: nlon,nlat,ndpth, ndepth_sed, ndepth_secchi
    integer(4)    :: lonvid, latvid, dpthvid
    integer(4)    :: lonvid_sed, latvid_sed, dpthvid_sed
    integer(4)    :: lonvid_chl, latvid_chl, dpthvid_chl
    integer(4)    :: lonvid_light, latvid_light, dpthvid_light
    integer(4)    :: lonvid_secchi, latvid_secchi, dpthvid_secchi
    integer(4)    :: tr,ios,iak,iben
    character(2)  :: ctmp
    character(60) :: filename, filename_sed, filename_chl, filename_light
    character(60) :: filename_secchi
    character(3)  :: c3  
    character(12) :: fnam_nml
    character(120):: std_name(0:ntbio)
    character(120):: standard_name_sed(0:ntben)
    character(120):: standard_name_chl, standard_name_light
    character(120):: standard_name_secchi
    character(17) :: var_unit(0:ntbio), var_unit_sed(0:ntben)
    logical, save :: first = .true., FirstEntry = .true.

    real(8), allocatable, save :: rlon(:), rlat(:), zlev(:),zlev_sed(:)
    integer(4), parameter      :: kmxb = 1  
    integer(4)                 :: ia
    character(2)               :: cncdf
  
    !- Since biofile_new is called from the coupler only once, the files 
    !- must be set for both areas at the same time 
    do ia =1,narea 
    
      if (first) then
        first = .false.
        allocate( ncid(narea))
        allocate( varid(ntbio,narea))
        allocate( timvid(narea))
        allocate( sedncid(narea))
        allocate( sedvarid(ntben,narea))
        allocate( sedtimvid(narea))
        if (do_diag_chl)then
          allocate( chlncid(narea))
          allocate( chlvarid(narea))
          allocate( chltimvid(narea))
          chlvarid = 0
          chltimvid = 0
          chlncid = (0)
        endif
        
        if (do_diag_secchi)then
          allocate( secchincid(narea))
          allocate( secchivarid(narea))
          allocate( secchitimvid(narea))
          secchivarid = 0
          secchitimvid = 0
          secchincid = (0)
        endif
        
        if(write_light)then
          allocate( lightncid(narea))
          allocate( lightvarid(narea))
          allocate( lighttimvid(narea))
          lightvarid = 0
          lightncid = 0
          lighttimvid = 0
        endif
        
        ncid   = 0
        varid  = 0
        timvid = 0
        sedncid = 0
        sedvarid = 0
        sedtimvid = 0
      endif

      if (FirstEntry) then
        if (ia == narea) FirstEntry = .false.
        ibncdf = 0
        if (archive) then
          filename        = 'biodat_'//trim(cfgfile(ia)(6:))//'_archive.nc'
          filename_sed    = 'bendat_'//trim(cfgfile(ia)(6:))//'_archive.nc'
          if (do_diag_chl)then
            filename_chl  = 'chldat_'//trim(cfgfile(ia)(6:))//'_archive.nc'
          endif
          if (do_diag_secchi)then
            filename_secchi  =                                                &
                       'secchidat_'//trim(cfgfile(ia)(6:))//'_archive.nc'
          endif
          if(write_light)then
            filename_light= 'light_'//trim(cfgfile(ia)(6:))//'_archive.nc'
          endif
        else
          write(cncdf,'(i2.2)') ibncdf
          filename        = 'biodat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
          filename_sed    = 'bendat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
          if (do_diag_chl)then
            filename_chl  = 'chldat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
          endif
          if (do_diag_secchi)then
            filename_secchi  =                                                &
                       'secchidat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
          endif
          if(write_light)then
            filename_light= 'light_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
          endif
        endif
      else
        if (ia == 1) ibncdf = ibncdf + 1
        write(cncdf,'(i2.2)') ibncdf
        filename        = 'biodat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
        filename_sed    = 'bendat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
        if (do_diag_chl)then
          filename_chl  = 'chldat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
        endif
        if (do_diag_secchi)then
          filename_secchi  =                                                  &
                     'secchidat_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
        endif
        if(write_light)then
          filename_light= 'light_'//trim(cfgfile(ia)(6:))//'_'//cncdf//'.nc'
        endif
      endif
     
      allocate( rlon(1:nmx(ia)),rlat(1:mmx(ia)), zlev(1:kmx(ia)))
      allocate( zlev_sed(1:kmxb))
      
      nlon = nmx(ia)
      nlat = mmx(ia)
      ndpth= kmx(ia)
      ndepth_sed = kmxb
      ndepth_secchi = kmxb

      !-------------------------------------------------------------------------
      !- Biovars
      !-------------------------------------------------------------------------
      call ncdf_create_new(filename,ncid(ia)) 
      call ncdf_set_dimensions(ncid(ia),nlon,lonid,nlat,latid,ndpth,dpthid,    &
                               timid)
      call ncdf_set_dim_vars(ncid(ia),lonid,lonvid,latid,latvid,dpthid,dpthvid,&
                           timid,timvid(ia))
      
      !- Get variable names and units
      call get_standard_name(std_name)
      call get_var_units(var_unit)
      
      do tr=1,ntbio 
       
        call ncdf_set_4D_variables(ncid(ia),trim(std_name(tr)),lonvid,         &
                                   latvid, dpthvid,timvid(ia),varid(tr,ia))  
      enddo
      
      
      call ncdf_set_dim_attr(ncid(ia),lonvid,latvid,dpthvid,timvid(ia))
      call ncdf_set_bio_var_attrib(ncid(ia), varid(:,ia), ntbio,std_name,      &
                                   var_unit)
     
      call ncdf_end_definitions(ncid(ia))
     
      !- Get lon,lat,depth for grid
      call arc2deg(mmx(ia),nmx(ia),xkoor(ia)%p,ykoor(ia)%p,rlon,         &
                   rlat)
      
      call calculate_depth( kmx(ia), dz(ia)%p, zlev)
    
      call ncdf_save_dim_vars(ncid(ia), lonvid, rlon, latvid, rlat,&
                              dpthvid, zlev)
     
      !-------------------------------------------------------------------------
      !- Sediment   
      !-------------------------------------------------------------------------
      call ncdf_create_new(filename_sed,sedncid(ia)) 
      call ncdf_set_dimensions(sedncid(ia),nlon,lonid_sed,nlat,latid_sed,      &
                             ndepth_sed,dpthid_sed,timid_sed)
      call ncdf_set_dim_vars(sedncid(ia),lonid_sed,lonvid_sed,latid_sed,       &
                           latvid_sed,dpthid_sed,dpthvid_sed,timid_sed,        &
                           sedtimvid(ia))
    
      call get_standard_name_sed(standard_name_sed)
      call get_var_units_sed(var_unit_sed)

      do tr=1,ntben
        call ncdf_set_3D_variables(sedncid(ia),trim(standard_name_sed(tr)),    &
                                 lonvid_sed,latvid_sed,sedtimvid(ia),          &
                                 sedvarid(tr,ia))  
      enddo

      call ncdf_set_dim_attr(sedncid(ia),lonvid_sed,latvid_sed,dpthvid_sed,    &
                           sedtimvid(ia))  
      !IL: 30.09.13 FIXME:These arguments and the nc routine are to be checked, 
      !                   maybe they can be replaced by a commom set_var_attr
      call ncdf_set_sed_var_attrib(sedncid(ia),sedvarid(:,ia),ntben,           &
                                   standard_name_sed,var_unit_sed)
      call ncdf_end_definitions(sedncid(ia))
    
      call calculate_depth( kmxb, dz(ia)%p, zlev_sed)
  
      call ncdf_save_dim_vars( sedncid(ia), lonvid_sed, rlon, latvid_sed,&
                             rlat,dpthvid_sed, zlev_sed )
    
     
      !-------------------------------------------------------------------------
      !- Chlorophyll   
      !-------------------------------------------------------------------------
      if (do_diag_chl)then
        call ncdf_create_new(filename_chl,chlncid(ia)) 
        call ncdf_set_dimensions(chlncid(ia),nlon,lonid_chl,nlat,latid_chl,    &
                                 ndpth,dpthid_chl,timid_chl)
        call ncdf_set_dim_vars(chlncid(ia),lonid_chl,lonvid_chl,latid_chl,     &
                            latvid_chl,dpthid_chl,dpthvid_chl,timid_chl,       &
                            chltimvid(ia))
    
        standard_name_chl='mass_concentration_of_chlorophyll_in_seawater'     
        call ncdf_set_4D_variables(chlncid(ia),trim(standard_name_chl),        &
                                   lonvid_chl,latvid_chl,dpthvid_chl,          &
                                   chltimvid(ia),chlvarid(ia))  
 
        call ncdf_set_dim_attr(chlncid(ia),lonvid_chl,latvid_chl,dpthvid_chl,  &
                               chltimvid(ia))
        call ncdf_set_chl_var_attrib(chlncid(ia),chlvarid(ia))
        call ncdf_end_definitions(chlncid(ia))
    
        call ncdf_save_dim_vars(chlncid(ia), lonvid_chl, rlon, latvid_chl, &
                            rlat,dpthvid_chl, zlev)
      endif 

      !-------------------------------------------------------------------------
      !- Secchi depth
      !-------------------------------------------------------------------------
      if (do_diag_secchi)then
        call ncdf_create_new(filename_secchi,secchincid(ia)) 
        call ncdf_set_dimensions(secchincid(ia),nlon,lonid_secchi,nlat,        &
                                 latid_secchi,ndepth_secchi,dpthid_secchi,     &
                                 timid_secchi)
        call ncdf_set_dim_vars(secchincid(ia),lonid_secchi,lonvid_secchi,      &
                               latid_secchi,latvid_secchi,dpthid_secchi,       &
                               dpthvid_secchi,timid_secchi,secchitimvid(ia))
    
        standard_name_secchi='secchi_depth_of_sea_water'     
        call ncdf_set_3D_variables(secchincid(ia),trim(standard_name_secchi),  &
                                   lonvid_secchi,latvid_secchi,                &
                                   secchitimvid(ia),secchivarid(ia))  
 
        call ncdf_set_dim_attr(secchincid(ia),lonvid_secchi,latvid_secchi,     &
                               dpthvid_secchi,secchitimvid(ia))
        call ncdf_set_secchi_var_attrib(secchincid(ia),secchivarid(ia))
        call ncdf_end_definitions(secchincid(ia))

       call calculate_depth( kmxb, dz(ia)%p, zlev_sed)

        call ncdf_save_dim_vars(secchincid(ia), lonvid_secchi, rlon,           &
                                latvid_secchi, rlat,dpthvid_secchi, zlev_sed)
      endif 
        
      !-------------------------------------------------------------------------
      !- Light
      !-------------------------------------------------------------------------
      if(write_light)then
        call ncdf_create_new(filename_light,lightncid(ia))
      
        call ncdf_set_dimensions(lightncid(ia),nlon,lonid_light,nlat,          &
                                 latid_light,ndpth,dpthid_light,timid_light)
                               
        call ncdf_set_dim_vars(lightncid(ia),lonid_light,lonvid_light,         &
                               latid_light, latvid_light, dpthid_light,        &
                               dpthvid_light, timid_light, lighttimvid(ia))
                             
        standard_name_light = 'photosynthetically_active_radiation'
      
        call ncdf_set_4D_variables(lightncid(ia),trim(standard_name_light),    &
                                 lonvid_light,latvid_light,dpthvid_light,      &
                                 lighttimvid(ia),lightvarid(ia))
                                 
        call ncdf_set_dim_attr(lightncid(ia),lonvid_light,latvid_light,        &
                               dpthvid_light, lighttimvid(ia))
      
        call ncdf_set_light_var_attrib(lightncid(ia),lightvarid(ia))
    
        call ncdf_end_definitions(lightncid(ia))
      
        call ncdf_save_dim_vars(lightncid(ia),lonvid_light, rlon, latvid_light,&
                                rlat,dpthvid_light, zlev)
      endif
    
      deallocate( rlon, rlat, zlev, zlev_sed )
    enddo
  
  end subroutine biofile_new

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

  subroutine biofile_write(ia, iu06,mmx, nmx, kmx, n2d, n3d, msrf, mcol, kh,   &
                           nlow, ptc, benthos)
  !Write all possible variables at the same time  

  use netcdf_interface, only : ncdf_save_3D_vars, ncdf_save_4D_vars,           &
                               ncdf_save_time
  use bioparam,         only : ntbio,Onorm,Nnorm,   &
                               idx_oxy            , &
                               ntben
  use writenetcdf_BSH,  only : round
  use constants,        only : tim
  use cmod_params,      only : narea
  
  implicit none 
  
  
  integer(4), intent(in) :: ia, iu06,n2d,n3d,nlow
  integer(4), intent(in) :: msrf(0:,0:), mcol(0:), kh(0:)
  real(8),    intent(in) :: ptc(0:,1:)    ! (0:iw3,1:ntbio)
  integer(4), intent(in) :: mmx,nmx,kmx
  real(8)   , intent(in) :: benthos(0:,1:)! (0:iw2,1:ntben)
    
  integer(4) :: dummy_4D(nmx,mmx,kmx)
  integer(4) :: dummy_3D(nmx,mmx)
  real(8)    :: roundtim
  integer(4) :: tr,m,n,k,mm,msf
  real(8)    :: scalebio
  integer(4), save   :: counter = 1
  real(8), parameter :: scalebioN = 1.e3_8*Nnorm
  real(8), parameter :: scalebioO = 1.e3_8*Onorm
 
  
  roundtim = tim
  
  call round(5,roundtim)
  ! BIOVARS
  call ncdf_save_time(ncid(ia),timvid(ia),roundtim,counter)
  ! BENTHOS
  call ncdf_save_time(sedncid(ia),sedtimvid(ia),roundtim,counter)
     
  
  ! BIOVARS  
  do tr=1,ntbio
    scalebio=scalebioN
    
    if (tr==idx_oxy) then
      scalebio=scalebioO
    endif
   
    dummy_4D(:,:,:) = -999
    !rearange grid
    do m=1,mmx
      do n=1,nmx
        msf = msrf(m,n)
        if (msf <= 0) cycle
        dummy_4D(n,m,1)= int(ptc(msf,tr)*scalebio,4)
        do k = 2,kh(msf)
          mm = mcol(msf)-2+k    
          dummy_4D(n,m,k)= int(ptc(mm,tr)*scalebio,4) 
        enddo
      enddo
    enddo
    call ncdf_save_4D_vars (ncid(ia), varid(tr,ia), counter, dummy_4D)
  enddo
  
   !BENTHOS
   !- rearange grid 
  do tr =1,ntben
    do m=1,mmx
      do n=1,nmx
        mm=msrf(m,n)
        if(mm==0)then
          dummy_3D(n,m) = -999
        else
          dummy_3D(n,m) = int(benthos(mm,tr)*scalebioN,4)
        endif
      enddo
    enddo
    call ncdf_save_3D_vars (sedncid(ia), sedvarid(tr,ia), counter, dummy_3D)
  enddo
  
  if(ia == narea) counter = counter + 1
 
  end subroutine biofile_write  
  
  !=============================================================================

  subroutine biofile_write_nc_diag(counter,mmx,nmx,kmx,msrf,mcol,kh,diag,      &
                                   ncid,timvid,varid)
  
  use constants,        only : tim
  use netcdf_interface, only : ncdf_save_time, ncdf_save_4D_vars
  use writenetcdf_BSH,  only : round
  
  implicit none
  
  
  integer(4), intent(in) :: counter,mmx,nmx,kmx
  integer(4), intent(in) :: msrf(0:,0:),mcol(0:),kh(0:)
  integer(4), intent(in) :: ncid,timvid,varid
  real(8),    intent(in) :: diag(0:)

  integer(4) :: mm, dummy_4D(nmx,mmx,kmx),m,n,k,msf
  real(8)    :: roundtim

  real(8), parameter :: scalebioN = 1.e3_8
  
  roundtim = tim
  call round(5,roundtim)
  call ncdf_save_time(ncid,timvid,roundtim,counter)  
  
  dummy_4D(:,:,:) = -999
  !rearrange grid 
  do m=1,mmx
    do n=1,nmx
      msf = msrf(m,n)
      if (msf <= 0) cycle
      dummy_4D(n,m,1)= int(diag(msf)*scalebioN,4)
      do k = 2,kh(msf) 
        mm = mcol(msf)-2+k    
        dummy_4D(n,m,k)= int(diag(mm)*scalebioN,4)
      enddo
    enddo
  enddo
  call ncdf_save_4D_vars (ncid, varid,counter,dummy_4D)
  
  end subroutine biofile_write_nc_diag
  

  !-------------------------------------------------------------
  subroutine biofile_write_nc_diag_2d(counter,mmx,nmx,msrf,diag,ncid,timvid,   &
                                      varid)
  
  use constants,        only : tim
  use netcdf_interface, only : ncdf_save_time, ncdf_save_3D_vars
  use writenetcdf_BSH,  only : round
  
  implicit none
  
  
  integer(4), intent(in) :: counter,mmx,nmx
  integer(4), intent(in) :: msrf(0:,0:)
  integer(4), intent(in) :: ncid,timvid,varid
  real(8),    intent(in) :: diag(0:)

  integer(4) :: mm, dummy_3D(nmx,mmx),m,n,k,msf
  real(8)    :: roundtim

  real(8), parameter :: scalebioN = 1.e3_8
  
  roundtim = tim
  call round(5,roundtim)
  call ncdf_save_time(ncid,timvid,roundtim,counter)  
  

   !- rearange grid 
  do m=1,mmx
    do n=1,nmx
      mm=msrf(m,n)
      if(mm==0)then
        dummy_3D(n,m) = -999
      else
        dummy_3D(n,m) = int(diag(mm)*scalebioN,4)
      endif
    enddo
  enddo
  call ncdf_save_3D_vars (ncid, varid, counter, dummy_3D)

  
  end subroutine biofile_write_nc_diag_2d

  !-------------------------------------------------------------
  

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

  subroutine biofile_close()
  
  use netcdf_interface, only: ncdf_close_file 
  use cmod_params,      only: narea
  use bioparam,         only: do_diag_chl, write_light, do_diag_secchi 
  
  integer(4) :: ia
  
  do ia = 1,narea
     call ncdf_close_file(ncid(ia))
     call ncdf_close_file(sedncid(ia))
     if (do_diag_chl)then
       call ncdf_close_file(chlncid(ia))
     endif
     if (do_diag_secchi)then
       call ncdf_close_file(secchincid(ia))
     endif
     if(write_light)then
       call ncdf_close_file(lightncid(ia))
     endif
     
     
  enddo
  end subroutine biofile_close 

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

  subroutine write_bio_restart(ia, iu06, doclose, mmx, nmx, kmx, iw2, iw3, nz, &
                               msrf, mcol, kh, nl, ptc, benthos, bndz)

    use bioparam,      only : Nnorm, Onorm, 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, ntben, ntbio_notag, ntben_notag,          &
                              write_binary_restart_notag
    use biopermute,    only : bio_depermute, eco_o, ben_o
    use cmod_params,   only : mainarea
    use io_subs,       only : io_new_unit
    use exits,         only : exitme
    use constants,     only : aufdat
    use cmod_miscvars, only : ihour,ctim
 
    implicit none

    integer(4), intent(in) :: ia, iu06, iw2, iw3, mmx, nmx, kmx, nz
    integer(4), intent(in) :: msrf(0:,0:), mcol(0:), kh(0:)
    integer(4), intent(in) :: nl
    logical,    intent(in) :: doclose
    real(8),    intent(in) :: ptc(0:,1:)  
    real(8),    intent(in) :: benthos(0:,1:)
    real(8),    intent(in) :: bndz(:,:,0:)

    integer(4), save :: restartunit = 0
    integer(4), save :: nz_main, kmx_main
    integer(4)       :: ios
    character(11)    :: aufdat_bio

    !  on first entry, open the restart data file:
    if (ia == 1) then
      restartunit = io_new_unit()
      open(restartunit,file='bio_restart',form='unformatted',action='write',   &
           asynchronous='yes',access='stream', iostat=ios)
      if (ios /= 0) call exitme (1, 'Could not open bio_restart' )
      aufdat_bio(1:11) = aufdat(1:11)
      aufdat_bio(2:5)  = ctim(1:4)
      aufdat_bio(6:10) = ctim(6:10)
      write(aufdat_bio(11:11),'(I1.1)') ihour/3
      write(iu06,'(a37,a11)') 'writing restart file bio_restart at: ',aufdat_bio
      write(restartunit) aufdat_bio
    endif

    !  de-permute eco and benthos before output:
    call bio_depermute( mmx,nmx,kmx,iw2,iw3,ntbio,ntben,nl,msrf,mcol,kh,ptc,   &
                        benthos )

    !  scale data before output:
    eco_o(0:iw3,1:ntbio) = eco_o(0:iw3,1:ntbio)*Nnorm
    eco_o(0:iw3,idx_oxy)     = eco_o(0:iw3,idx_oxy)/NOR
    ben_o(0:iw2,1:ntben)   = ben_o(0:iw2,1:ntben)*Nnorm

    !  write data to restart file:
    if (write_binary_restart_notag) then
      write(restartunit,iostat=ios) eco_o(0:iw3,1:ntbio_notag), ben_o(0:iw2,1:ntben_notag)
    else
      write(restartunit,iostat=ios) eco_o(0:iw3,1:ntbio), ben_o(0:iw2,1:ntben)
    endif
    if (ios /= 0) then
      write(iu06,*) 'I/O problem - iostat is:',ios
      write(iu06,*) '                   area:',ia
      call exitme (2, 'Could not write bio_restart' )
    endif

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

    if (doclose) then
      if (nz_main >= 1) then
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        write(restartunit,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 write bio sponges' )
        endif
        if (.not. write_binary_restart_notag) then
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
          write(restartunit,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 write bio sponges' )
          endif
        endif
      endif

      !  on last entry, close restart data file:
      close(restartunit)
    endif

  end subroutine write_bio_restart

  !=============================================================================
  subroutine val_benthos(ia, iu06, n2d, benthos, ntben, header)

    use constants, only : zero, one
    use io_subs,   only : flush_unit

    implicit none

    integer(4), intent(in) :: ia, iu06, n2d, ntben
    real(8),    intent(in) :: benthos(0:,1:)
    logical,    intent(in) :: header

    integer(4) :: ib, iben
    real(8)    :: bmean, fac1, fac2

    if (header) write(iu06,'(a26)') 'Benthos validation prints:'
    write(iu06,'(a72)')                                                        &
      '------------------------------------------------------------------------'
    write(iu06,'(a33,i5)') 'Statistics for domain:           ', ia
    
    do iben=1,ntben
      write(iu06,'(a33,i5)') '  Benthos variable:              ', iben
      bmean = zero
      do ib=0,n2d-1
        fac2  = real(ib,8)
        fac1  = one/(fac2 + one)
        bmean = (bmean*fac2 + benthos(ib+1,iben))*fac1
      enddo
      write(iu06,'(a33,f16.10)')           '    Average value:               ',&
                                        bmean
      write(iu06,'(a33,f16.10,5x,f16.10)') '    Min and max:                 ',&
                  minval(benthos(1:,iben)), maxval(benthos(1:,iben))
    enddo
    write(iu06,'(a72)')                                                        &
    '------------------------------------------------------------------------'
    call flush_unit(iu06)

  end subroutine val_benthos
  !=============================================================================
   
  subroutine biofile_next ( iu06 )
    implicit none

    integer(4), intent(in) :: iu06

    call biofile_close ()

    write(iu06,'(a22,a)') 'bio output stored on nc file '

    call biofile_new ()

  end subroutine biofile_next

  !=============================================================================
  subroutine get_standard_name(std_name)
  
    use bioparam, only :                                 &
                        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,            &
                        nc_varname_nml_out_amm            ,            &
                        nc_varname_nml_out_nit            ,            &
                        nc_varname_nml_out_phos           ,            &
                        nc_varname_nml_out_sil            ,            &
                        nc_varname_nml_out_dia            ,            &
                        nc_varname_nml_out_flag           ,            &
                        nc_varname_nml_out_cyano          ,            &
                        nc_varname_nml_out_mez            ,            &
                        nc_varname_nml_out_miz            ,            &
                        nc_varname_nml_out_det            ,            &
                        nc_varname_nml_out_dets           ,            &
                        nc_varname_nml_out_ldon           ,            &
                        nc_varname_nml_out_oxy            ,            &
                        nc_varname_nml_out_amm_with_totatmos_N,            &
                        nc_varname_nml_out_nit_with_totatmos_N,            &
                        nc_varname_nml_out_dia_with_totatmos_N,            &
                        nc_varname_nml_out_flag_with_totatmos_N,            &
                        nc_varname_nml_out_cyano_with_totatmos_N,            &
                        nc_varname_nml_out_mez_with_totatmos_N,            &
                        nc_varname_nml_out_miz_with_totatmos_N,            &
                        nc_varname_nml_out_det_with_totatmos_N,            &
                        nc_varname_nml_out_ldon_with_totatmos_N,            &
                        nc_varname_nml_out_amm_with_agri_N,            &
                        nc_varname_nml_out_nit_with_agri_N,            &
                        nc_varname_nml_out_dia_with_agri_N,            &
                        nc_varname_nml_out_flag_with_agri_N,            &
                        nc_varname_nml_out_cyano_with_agri_N,            &
                        nc_varname_nml_out_mez_with_agri_N,            &
                        nc_varname_nml_out_miz_with_agri_N,            &
                        nc_varname_nml_out_det_with_agri_N,            &
                        nc_varname_nml_out_ldon_with_agri_N,            &
                        ntbio, use_nc_varname_nml_out
    
    implicit none
    
    character(len=120), intent(out) :: std_name(0:)

    if (use_nc_varname_nml_out) then
      std_name(idx_amm            )  = nc_varname_nml_out_amm            
      std_name(idx_nit            )  = nc_varname_nml_out_nit            
      std_name(idx_phos           )  = nc_varname_nml_out_phos           
      std_name(idx_sil            )  = nc_varname_nml_out_sil            
      std_name(idx_dia            )  = nc_varname_nml_out_dia            
      std_name(idx_flag           )  = nc_varname_nml_out_flag           
      std_name(idx_cyano          )  = nc_varname_nml_out_cyano          
      std_name(idx_mez            )  = nc_varname_nml_out_mez            
      std_name(idx_miz            )  = nc_varname_nml_out_miz            
      std_name(idx_det            )  = nc_varname_nml_out_det            
      std_name(idx_dets           )  = nc_varname_nml_out_dets           
      std_name(idx_ldon           )  = nc_varname_nml_out_ldon           
      std_name(idx_oxy            )  = nc_varname_nml_out_oxy            
      std_name(idx_amm_with_totatmos_N)  = nc_varname_nml_out_amm_with_totatmos_N
      std_name(idx_nit_with_totatmos_N)  = nc_varname_nml_out_nit_with_totatmos_N
      std_name(idx_dia_with_totatmos_N)  = nc_varname_nml_out_dia_with_totatmos_N
      std_name(idx_flag_with_totatmos_N)  = nc_varname_nml_out_flag_with_totatmos_N
      std_name(idx_cyano_with_totatmos_N)  = nc_varname_nml_out_cyano_with_totatmos_N
      std_name(idx_mez_with_totatmos_N)  = nc_varname_nml_out_mez_with_totatmos_N
      std_name(idx_miz_with_totatmos_N)  = nc_varname_nml_out_miz_with_totatmos_N
      std_name(idx_det_with_totatmos_N)  = nc_varname_nml_out_det_with_totatmos_N
      std_name(idx_ldon_with_totatmos_N)  = nc_varname_nml_out_ldon_with_totatmos_N
      std_name(idx_amm_with_agri_N)  = nc_varname_nml_out_amm_with_agri_N
      std_name(idx_nit_with_agri_N)  = nc_varname_nml_out_nit_with_agri_N
      std_name(idx_dia_with_agri_N)  = nc_varname_nml_out_dia_with_agri_N
      std_name(idx_flag_with_agri_N)  = nc_varname_nml_out_flag_with_agri_N
      std_name(idx_cyano_with_agri_N)  = nc_varname_nml_out_cyano_with_agri_N
      std_name(idx_mez_with_agri_N)  = nc_varname_nml_out_mez_with_agri_N
      std_name(idx_miz_with_agri_N)  = nc_varname_nml_out_miz_with_agri_N
      std_name(idx_det_with_agri_N)  = nc_varname_nml_out_det_with_agri_N
      std_name(idx_ldon_with_agri_N)  = nc_varname_nml_out_ldon_with_agri_N
    else
      std_name(idx_amm            )  = 'mole_concentration_of_ammonium_in_seawater'
      std_name(idx_nit            )  = 'mole_concentration_of_nitrate_in_seawater'
      std_name(idx_phos           )  = 'mole_concentration_of_phosphate_in_seawater'
      std_name(idx_sil            )  = 'mole_concentration_of_silicate_in_seawater'
      std_name(idx_dia            )  = 'mole_concentration_of_diatoms_expressed_as_nitrogen_in_seawater'
      std_name(idx_flag           )  = 'mole_concentration_of_flagellates_expressed_as_nitrogen_in_seawater'
      std_name(idx_cyano          )  = 'mole_concentration_of_cyanobacteria_expressed_as_nitrogen_in_seawater'
      std_name(idx_mez            )  = 'mole_concentration_of_zooplankton_expressed_as_nitrogen_in_seawater'
      std_name(idx_miz            )  = 'mole_concentration_of_protozooplankton_expressed_as_nitrogen_in_seawater'
      std_name(idx_det            )  = 'mole_concentration_of_detritus_expressed_as_nitrogen_in_seawater'
      std_name(idx_dets           )  = 'mole_concentration_of_detritus_expressed_as_silicate_in_seawater'
      std_name(idx_ldon           )  = 'mole_concentration_of_labile_dissolved_organic_nitrogen_expressed_as_nitrogen_in_seawater'
      std_name(idx_oxy            )  = 'mole_concentration_of_oxygen_in_seawater'
      std_name(idx_amm_with_totatmos_N)  = 'mole_concentration_of_ammonium_from_atmospheric_deposition_in_seawater'
      std_name(idx_nit_with_totatmos_N)  = 'mole_concentration_of_nitrate_from_atmospheric_deposition_in_seawater'
      std_name(idx_dia_with_totatmos_N)  = 'mole_concentration_of_diatoms_expressed_as_nitrogen_from_atmospheric_deposition_in_seawater'
      std_name(idx_flag_with_totatmos_N)  = 'mole_concentration_of_flagellates_expressed_as_nitrogen_from_atmospheric_deposition_in_seawater'
      std_name(idx_cyano_with_totatmos_N)  = 'mole_concentration_of_cyanobacteria_expressed_as_nitrogen_from_atmospheric_deposition_in_seawater'
      std_name(idx_mez_with_totatmos_N)  = 'mole_concentration_of_zooplankton_expressed_as_nitrogen_from_atmospheric_deposition_in_seawater'
      std_name(idx_miz_with_totatmos_N)  = 'mole_concentration_of_protozooplankton_expressed_as_nitrogen_from_atmospheric_deposition_in_seawater'
      std_name(idx_det_with_totatmos_N)  = 'mole_concentration_of_detritus_expressed_as_nitrogen_from_atmospheric_deposition_in_seawater'
      std_name(idx_ldon_with_totatmos_N)  = 'mole_concentration_of_labile_dissolved_organic_nitrogen_from_atmospheric_deposition_expressed_as_nitrogen_in_seawater'
      std_name(idx_amm_with_agri_N)  = 'mole_concentration_of_ammonium_from_agricultural_activities_in_seawater'
      std_name(idx_nit_with_agri_N)  = 'mole_concentration_of_nitrate_from_agricultural_activities_in_seawater'
      std_name(idx_dia_with_agri_N)  = 'mole_concentration_of_diatoms_expressed_as_nitrogen_from_agricultural_activities_in_seawater'
      std_name(idx_flag_with_agri_N)  = 'mole_concentration_of_flagellates_expressed_as_nitrogen_from_agricultural_activities_in_seawater'
      std_name(idx_cyano_with_agri_N)  = 'mole_concentration_of_cyanobacteria_expressed_as_nitrogen_from_agricultural_activities_in_seawater'
      std_name(idx_mez_with_agri_N)  = 'mole_concentration_of_zooplankton_expressed_as_nitrogen_from_agricultural_activities_in_seawater'
      std_name(idx_miz_with_agri_N)  = 'mole_concentration_of_protozooplankton_expressed_as_nitrogen_from_agricultural_activities_in_seawater'
      std_name(idx_det_with_agri_N)  = 'mole_concentration_of_detritus_expressed_as_nitrogen_from_agricultural_activities_in_seawater'
      std_name(idx_ldon_with_agri_N)  = 'mole_concentration_of_labile_dissolved_organic_nitrogen_from_agricultural_activities_expressed_as_nitrogen_in_seawater'
    endif
  
  end subroutine get_standard_name
 
 !==============================================================================
  subroutine get_standard_name_sed(std_name_sed)
    
    use bioparam, only :                                  &
                         bidx_nitr           ,            &
                         bidx_sili           ,            &
                         bidx_nitr_with_totatmos_N,            &
                         bidx_nitr_with_agri_N,            &
                         nc_varname_nml_out_ben_nitr           ,            &
                         nc_varname_nml_out_ben_sili           ,            &
                         nc_varname_nml_out_ben_nitr_with_totatmos_N,            &
                         nc_varname_nml_out_ben_nitr_with_agri_N,            &
                         ntben, use_nc_varname_nml_out_ben
    
    implicit none
    
    character(len=120), intent(out) :: std_name_sed(0:)
    
    if (use_nc_varname_nml_out_ben) then
      std_name_sed(bidx_nitr           )  = nc_varname_nml_out_ben_nitr           
      std_name_sed(bidx_sili           )  = nc_varname_nml_out_ben_sili           
      std_name_sed(bidx_nitr_with_totatmos_N)  = nc_varname_nml_out_ben_nitr_with_totatmos_N
      std_name_sed(bidx_nitr_with_agri_N)  = nc_varname_nml_out_ben_nitr_with_agri_N
    else
      std_name_sed(bidx_nitr           )  = 'mole_concentration_of_nitrogen_in_the_sediment'
      std_name_sed(bidx_sili           )  = 'mole_concentration_of_silicate_in_the_sediment'
      std_name_sed(bidx_nitr_with_totatmos_N)  = 'mole_concentration_of_nitrogen_from_atmospheric_deposition_in_the_sediment'
      std_name_sed(bidx_nitr_with_agri_N)  = 'mole_concentration_of_nitrogen_from_agricultural_activities_in_the_sediment'
    endif
    
    
  end subroutine get_standard_name_sed
 
 !==============================================================================

 subroutine get_var_units(var_unit)
  
  use bioparam, only :                                 &
                       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  
  
  implicit none
  
  character(len=17), intent(out) :: var_unit(0:)

  var_unit(idx_amm            )  = 'micromol N/m3'
  var_unit(idx_nit            )  = 'micromol N/m3'
  var_unit(idx_phos           )  = 'micromol P/m3'
  var_unit(idx_sil            )  = 'micromol Si/m3'
  var_unit(idx_dia            )  = 'micromol N/m3'
  var_unit(idx_flag           )  = 'micromol N/m3'
  var_unit(idx_cyano          )  = 'micromol N/m3'
  var_unit(idx_mez            )  = 'micromol N/m3'
  var_unit(idx_miz            )  = 'micromol N/m3'
  var_unit(idx_det            )  = 'micromol N/m3'
  var_unit(idx_dets           )  = 'micromol N/m3'
  var_unit(idx_ldon           )  = 'micromol N/m3'
  var_unit(idx_oxy            )  = 'micromol O/m3'
  var_unit(idx_amm_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_nit_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_dia_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_flag_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_cyano_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_mez_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_miz_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_det_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_ldon_with_totatmos_N)  = 'micromol N/m3'
  var_unit(idx_amm_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_nit_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_dia_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_flag_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_cyano_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_mez_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_miz_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_det_with_agri_N)  = 'micromol N/m3'
  var_unit(idx_ldon_with_agri_N)  = 'micromol N/m3'
 
 
 end subroutine get_var_units

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

 subroutine get_var_units_sed(var_unit_sed)
  
  use bioparam, only :                                  &
                       bidx_nitr           ,            &
                       bidx_sili           ,            &
                       bidx_nitr_with_totatmos_N,            &
                       bidx_nitr_with_agri_N,            &
                       ntben
  
  implicit none
  
  character(len=17), intent(out) :: var_unit_sed(0:)

  var_unit_sed(bidx_nitr           )  = 'micromol N/m2'
  var_unit_sed(bidx_sili           )  = 'micromol Si/m2'
  var_unit_sed(bidx_nitr_with_totatmos_N)  = 'micromol N/m2'
  var_unit_sed(bidx_nitr_with_agri_N)  = 'micromol N/m2'
 
 
 end subroutine get_var_units_sed
 
 !==============================================================================

 subroutine diag_out(ia,idiag,diag,msrf,mcol,kh,mmx,nmx,kmx)
  !subroutine to coordinate output of diagnostic variables
  
  use bioparam,    only : didx_chloro, do_diag_chl, write_light, didx_secchi,  &
                          do_diag_secchi
  use cmod_params, only : narea
  
  implicit none
  
  integer(4), intent(in) :: ia
  integer(4), intent(in) :: idiag
  integer(4), intent(in) :: msrf(0:,0:),mcol(0:),kh(0:)
  integer(4), intent(in) :: mmx,nmx,kmx
  real(8),    intent(in) :: diag(0:) !(1:iw3)
  
  integer(4), save :: counterl = 1, counterc = 1, counters = 1
  
  if(write_light .and. (idiag == 1)) then
    call biofile_write_nc_diag(counterl,mmx,nmx,kmx,msrf,mcol,kh,diag,         &
                               lightncid(ia),lighttimvid(ia),lightvarid(ia))
    if (ia == narea) counterl = counterl + 1
  endif
  if (do_diag_chl .and. (didx_chloro == idiag)) then
    call biofile_write_nc_diag(counterc,mmx,nmx,kmx,msrf,mcol,kh,diag,         &
                               chlncid(ia),chltimvid(ia),chlvarid(ia))
    if (ia == narea) counterc = counterc + 1
  endif
  if (do_diag_secchi .and. (didx_secchi == idiag)) then
    call biofile_write_nc_diag_2d(counters,mmx,nmx,msrf,diag,secchincid(ia),   &
                                  secchitimvid(ia),secchivarid(ia))
    if (ia == narea) counters = counters + 1
  endif

 end subroutine diag_out 

end module biooutput

