module biodiagnostic
!-------------------------------------------------------------------------------
! Module for storing and calculating the diagnostic ecoparameters
!-------------------------------------------------------------------------------

  use constants, only : zero, one

  implicit none

  !- public routines
  public  :: calc_chlorophyll, calc_secchi_depth
 
 ! - private routines 
  

 ! FIXME NEUMANND: create real template

contains


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

  subroutine calc_chlorophyll(iw2, mcol, kh, idx, eco, light, chl)

    !---------------------------------------------------------------------------
    ! Routine from BSH used for calculating the diagnostic variables. Based on
    ! Doron et al 2013._8
    !---------------------------------------------------------------------------

    use bioparam, only : idx_dia, idx_flag, idx_cyano, Nnorm
    use dmi_omp,  only : domp_get_domain 

    implicit none

    !- args
    real(8),    intent(in)    :: light(0:)  !(0:iw3)
    integer(4), intent(in)    :: iw2
    integer(4), intent(in)    :: mcol(0:)
    integer(4), intent(in)    :: kh(0:)
    integer(4), intent(inout) :: idx(:,:)
    real(8),    intent(in)    :: eco(0:,1:) ! tracer array (0:iw3,1:ntbio)
    real(8),    intent(inout) :: chl(0:)  !(0:iw3)
   
    !- params
    real(8), parameter  :: min_chl_to_n = 2.00_8  !mg Chl/mmol N
    real(8), parameter  :: max_chl_to_n = 3.00_8 !mg Chl/mmol N 
    real(8), parameter  :: light_max = 30.0_8 !W/m2

    !- vars
    integer(4) :: nsrf, mi0, mi, kb, k, n2dl, n2du
    real(8)    :: r_chl_to_n, total_phyto, frac
   
    real(8)  :: lightk, chlk   
    

    
    
    !- constant in the loop above so pre-calculate it
    frac = max_chl_to_n / (min_chl_to_n*light_max)
    
    call domp_get_domain(kh, 1, iw2, n2dl, n2du, idx)

    do nsrf=n2dl,n2du
       total_phyto = zero
       lightk = light(nsrf) 
      !- dia + flag + cyano
      total_phyto = (eco(nsrf,idx_dia)+eco(nsrf,idx_flag)+                     &
                     eco(nsrf,idx_cyano) )*Nnorm
      !- eq. for chl to N ratio
      if(lightk < zero)then
        chlk = zero
      else
        r_chl_to_n = max( min_chl_to_n, max_chl_to_n*(one-frac*lightk) )
        !- calc chl
        chlk  = r_chl_to_n*total_phyto
        
      endif
      chl(nsrf) = chlk
      
      kb = kh(nsrf)
      if (kb > 1) mi0 = mcol(nsrf)
      do k = 2,kb
        mi = mi0+k-2
        total_phyto = zero
        lightk = light(mi)
        !- dia + flag + cyano
        total_phyto=(eco(mi,idx_dia)+eco(mi,idx_flag)+                         &
                     eco(mi,idx_cyano))*Nnorm
        !- eq. for chl to N ratio
        if(lightk < zero)then
          chlk = zero
        else
          r_chl_to_n = max( min_chl_to_n, max_chl_to_n*(one-frac*lightk) )

          !- calc chl          
          chlk    = r_chl_to_n*total_phyto
        endif
        chl(mi) = chlk
      enddo
    enddo
   
  end subroutine calc_chlorophyll
!==============================================================================
  subroutine calc_secchi_depth (iw2,kh,idx,kpar,secchi)
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! a subroutine for diagnostic calculation of Secchi depth
    ! formulation as first attempt after Neumann et al 2015
    ! proven to work for the Baltic Sea
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    use bioparam, only : secchi_exp, secchi_coeff
    use dmi_omp,  only : domp_get_domain 

    implicit none

    integer(4), intent(in)    :: iw2
    integer(4), intent(in)    :: kh(0:)
    integer(4), intent(inout) :: idx(:,:)
    real(8),    intent(in)    :: kpar(0:)    !(0:iw3)
    real(8),    intent(inout) :: secchi(0:)  !(0:iw2)
    
    integer(4) :: n2dl, n2du, nsrf

    call domp_get_domain(kh, 1, iw2, n2dl, n2du, idx)

    !FIXME: Do we need to set the land value every time - one sensible init
    !       should be enough
!$OMP MASTER
    secchi(0) = zero
!$OMP END MASTER

    do nsrf=n2dl,n2du
      if(kpar(nsrf)<zero)then
        secchi(nsrf)=zero
      else
        secchi(nsrf) = secchi_coeff*kpar(nsrf)**secchi_exp
      endif
    enddo 
    
  end subroutine calc_secchi_depth
end module biodiagnostic
