 MODULE biolodyncellsed

#if defined key_substance && defined key_biolo && defined key_sedim_MUSTANG
   !!======================================================================
   !!                   ***  MODULE  biolodyn  ***
   !! Bilogic dynamics in sediment colomn
   !!======================================================================
   !! * Modules used
   USE parameters,   ONLY : rsh,rlg,epsilon
   USE comMUSTANG,     ONLY : poro,ksma,dzs
   USE comsubstance, ONLY : nv_state,dtbio,irkm_var_assoc,cmes_3dmgl
   USE bioloinit

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC biolo_dyn_cellsed ! routine called by biolo_sksc_sed.F90

   !! * Interface
   INTERFACE biolo_dyn_cellsed
!     MODULE PROCEDURE incellsed_dinodiatzoo
     MODULE PROCEDURE incellsed_nut4phy3zoo2
   END INTERFACE

   !! * Shared module variables

   !! * Private variables

 CONTAINS
   !!======================================================================

  SUBROUTINE incellsed_nut4phy3zoo2(km,im,jm,txfiltbenth,epn,cw,cs,dc,dcw_filtbent)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellsed_nut4phy3zoo2  ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes biologiques
   !&E              ayant lieu a l interieur d une maille dans le sediment
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  (M. Sourisseau) Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvarp2d, ONLY  : temp
   USE comvars2d, ONLY  : h0fond,dt,t,tdeb


   !! * Declaration subroutine

   !! * Arguments
   INTEGER, INTENT( in )                     :: im,jm,km    ! cell indexes
   REAL(KIND=rsh),DIMENSION(1:nv_state),INTENT(in)    :: cw
   REAL(KIND=rsh),DIMENSION(-1:nv_state),INTENT(in)    :: cs
   REAL(KIND=rsh),DIMENSION(:),INTENT(inout) :: dc,dcw_filtbent
   REAL(KIND=rsh),INTENT(in)                 :: txfiltbenth,epn

   !! * Local declarations
   REAL(KIND=rlg)      :: dtbiojour
   INTEGER               :: iv          ! loop index
                     ! variables relatives a la temperature
   REAL(KIND=rsh)        :: temper,effetchaleur
                     ! variable diatomees
   REAL(KIND=rsh)        :: diatmortsed
                     ! variables sels nutritifs
   REAL(KIND=rsh)        :: rappaz,rapsiaz,rappdiatiksmi
                     ! processus mineralisation, nitrification, dissolution
   REAL(KIND=rsh)        :: facto2remin,facto2nitrif,reminazdetsed,     &
                            xnitrifsed,reminpdetsed,dissolsilicesed
                     ! processus adsorption/desorption
   REAL(KIND=rsh)        :: varads,satursedp,adsorpsed,cmaxadsorp,      &
                            desorpsed,cmaxdesorp
                     ! cycle de l oxygene
   REAL(KIND=rsh)        :: oreminsed,onitrifsed
                     ! porosite locale
   REAL(KIND=rsh)        :: porosite
                     ! filtrage broutage benthique
   REAL(KIND=rsh)        :: filtr_benth,txfiltbenthij,factlim,txfiltbenthij_wat


   !!----------------------------------------------------------------------
   !! * Executable part

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! saisir le code apres cette limite
   ! =================================
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   dtbiojour=dtbio/86400.0_rlg
   porosite=poro(km,im,jm) 
   temper=cs(-1)
   IF (temper.gt.30.) THEN
     write(*,*)'temp sed>30:',km,im,jm,temper
     temper=15.
   END IF
   effetchaleur=exp(p_T_effect*temper)

   ! mortalite des diatomees
   ! -----------------------
   diatmortsed=p_diat_mort_sed*effetchaleur
   !diatmortsed=0.3_rsh*p_diat_mortmax*effetchaleur

   ! remineralisation dans le sediment
   ! ---------------------------------
   !facto2remin=c(iv_nutr_NH4)/(c(iv_nutr_NH4)+(p_micz_mort/2))
   !facto2nitrif=c(iv_nutr_NH4)/(c(iv_nutr_NH4)+p_micz_mort)
#if defined key_oxygen
   IF(cs(iv_oxygen)> 0.2_rsh) THEN
     facto2remin=cs(iv_oxygen)/(cs(iv_oxygen)+p_O2_kremin)
   ELSE
     facto2remin=0.0_rsh
   ENDIF
#else
   facto2remin=1.0_rsh
#endif
   reminazdetsed=facto2remin*p_Nbent_remin*effetchaleur
   !reminazdetsed=facto2remin*p_reminbenth*p_N_remin*effetchaleur

   reminpdetsed=facto2remin*p_Pbent_remin*effetchaleur
   !reminpdetsed=facto2remin*p_reminbenth*p_P_remin*effetchaleur

   dissolsilicesed=facto2remin*p_Sibent_dissol*effetchaleur
   !dissolsilicesed=p_reminbenth*p_Si_diss*effetchaleur

   ! nitrification dans le sediment
   ! --------------------------------
   facto2nitrif=1.0_rsh
#if defined key_oxygen
   IF(cs(iv_oxygen)< 0.2_rsh)facto2nitrif=0.0_rsh
#endif   
   xnitrifsed=facto2nitrif*p_nitrifsed*effetchaleur

   rapsiaz=p_phyto_SiNratio
   rappaz=1.0_rsh/p_phyto_NPratio

   ! Adsorption et desorption du phosphore dans le sediment
   ! ------------------------------------------------------

   ! adsorption
   ! QUESTION : memes parametres que dans leau ??
   ! ----------
   IF (cs(irkm_var_assoc(iv_nutr_Pads)).ne.0.0) THEN
     varads=min(1.0,cs(iv_nutr_Pads)/(p_P_adsormaxspim*cs(irkm_var_assoc(iv_nutr_Pads))))
   ELSE
     varads=0.0
   END IF
   satursedp=max(0.,(p_P_adsormaxspim*cs(irkm_var_assoc(iv_nutr_Pads))-cs(iv_nutr_Pads)))
   adsorpsed=p_P_adsor*satursedp               
   cmaxadsorp=adsorpsed*cs(iv_nutr_PO4)*dtbiojour
   IF (cmaxadsorp.ge.cs(iv_nutr_PO4))  adsorpsed=0.9/dtbiojour                   
   ! desorption
   ! ----------
   desorpsed=p_P_desor*varads
   cmaxdesorp=desorpsed*cs(iv_nutr_Pads)*dtbiojour
   IF (cmaxdesorp.ge.cs(iv_nutr_Pads)) desorpsed=0.9/dtbiojour

   txfiltbenthij=0.0_rsh
   filtr_benth=0.0_rsh
   
#if ! defined key_biolo_opt2
   ! pression de broutage (filtrage / broutage )
   IF(km == ksma(im,jm) .AND. epn > h0fond ) THEN

     !IF(l_filtbenthsinusmes) THEN
     ! Introduction d une pression de Broutage, d apres Marie SAVINA
     ! -------------------------------------------------------------  
     ! txfiltbenth = taux sinusoidal de filtration du benthos en m3/j/m2
     !foncsinus=(1.0_rsh+SIN(2.0_rsh*3.14159_rsh/365.0_rsh*(jjulien-125)))/2.0_rsh
     !txfiltbenth=(p_txfiltbenthmax*(0.3_rsh+0.7_rsh*foncsinus))

     ! and txfiltbenth fonction des MES (filtration baisse avec la turbidité)
       txfiltbenthij=MAX(-0.0333_rsh*cmes_3dmgl(1,im,jm)+1,0.0_rsh)*txfiltbenth  !de 0 mg/L a 30 mg/L
     ! conversion par metre3 eau
       txfiltbenthij_wat=txfiltbenthij/epn
    ! ENDIF
       ! resolution explicite pour cw ==> limitation pas de temps
       factlim=MIN(1.0_rsh,1.0_rsh/dtbiojour/(txfiltbenthij_wat+epsilon))
       IF(factlim .NE. 1.0_rsh) write(*,*)'filt',(t-tdeb)/3600._rsh,im,jm,km,epn,factlim,txfiltbenthij_wat,dtbiojour
       txfiltbenthij_wat=txfiltbenthij_wat*factlim
       txfiltbenthij=txfiltbenthij*factlim/dzs(km,im,jm)
       filtr_benth= txfiltbenthij*(cw(iv_phyto_diat_N)+cw(iv_phyto_dino_N)+cw(iv_phyto_nano_N))
    ! filtrage sur les concentrations dans l eau
       dcw_filtbent(iv_phyto_diat_N)=-txfiltbenthij_wat*cw(iv_phyto_diat_N)
       dcw_filtbent(iv_phyto_dino_N)=-txfiltbenthij_wat*cw(iv_phyto_dino_N)
       dcw_filtbent(iv_phyto_nano_N)=-txfiltbenthij_wat*cw(iv_phyto_nano_N)
#if defined key_psnz
       dcw_filtbent(iv_phyto_psnz_N)=-txfiltbenthij_wat*cw(iv_phyto_psnz_N)
       dcw_filtbent(iv_phyto_psnz_Si)=-txfiltbenthij_wat*cw(iv_phyto_psnz_Si)
#endif
#if defined key_karenia
       dcw_filtbent(iv_phyto_karenia_N)=-txfiltbenthij_wat*cw(iv_phyto_karenia_N)
       dcw_filtbent(iv_phyto_karenia_C)=-txfiltbenthij_wat*cw(iv_phyto_karenia_C)
#endif
#ifdef key_phaeocystis
       dcw_filtbent(iv_phyto_phaeocystis_cell_N)=-txfiltbenthij_wat*cw(iv_phyto_phaeocystis_cell_N)
#endif
#ifdef key_ulvas
       dcw_filtbent(iv_ulv_benth_N)=-ulvebenthmortsed*cs(iv_ulv_benth_N)/0.014_rsh
       dcw_filtbent(iv_ulv_benth_P)=-ulvebenthmortsed*cs(iv_ulv_benth_P)/0.014_rsh
#endif

   ENDIF
#endif

#if defined key_oxygen
   ! production et consommation d oxygene dans le sediment 
   ! ------------------------------------------------------

   ! remineralisation
   !oreminsed=reminazdetsed*cs(iv_detr_N)*(0.212+0.064) 
   ! QUESTIONS !!!!!!!!!
   ! 0.212+0.064 ???
   oreminsed=reminazdetsed*0.212_rsh
   ! pourquoi le N du meso qui se nitrifi ???? bizarre non ?     
   !onitrifsed=xnitrifsed*c(iv_zoo_meso_N)*0.064
   onitrifsed=xnitrifsed*0.064_rsh

   ! denitrification  ??
   ! IF((cs(iv_oxygen)/o2sat).le.0.05_rsh) THEN
   !
   !   dc(iv_nutr_NO3)=dc(iv_nutr_NO3)-((-0.2_rsh/o2sat)*cs(iv_oxygen)+0.2_rsh)*cs(iv_nutr_NO3)
   !  END IF
#endif


   ! Evolution de l ammonium
   ! -------------------------------
   dc(iv_nutr_NH4)=reminazdetsed*cs(iv_detr_N)/porosite            &
#if defined key_zostera
                  +(reminpdetsed/10._rsh)*cs(iv_detr_N)/porosite   &
#endif
                  -xnitrifsed*cs(iv_nutr_NH4)


   ! Evolution du nitrate
   ! -------------------------------
   dc(iv_nutr_NO3)=xnitrifsed*cs(iv_nutr_NH4)

   ! Evolution de la silice dissoute
   ! -------------------------------
   dc(iv_nutr_SiOH)=dissolsilicesed*cs(iv_detr_Si)/porosite

   ! Evolution du phosphate dissous
   ! -------------------------------
   dc(iv_nutr_PO4)=-adsorpsed*cs(iv_nutr_PO4)                        &
                    +desorpsed*cs(iv_nutr_Pads)/porosite             &
#if defined key_zostera
                    +(reminpdetsed/10._rsh)*cs(iv_detr_P)/porosite   &
#endif
                   +reminpdetsed*cs(iv_detr_P)/porosite

   ! Evolution du phosphate adsorbe
   ! -------------------------------
   dc(iv_nutr_Pads)=adsorpsed*cs(iv_nutr_PO4)*porosite-desorpsed*cs(iv_nutr_Pads)

   ! Evolution de l azote des diatomees
   ! -------------------------------
   dc(iv_phyto_diat_N)=-diatmortsed*cs(iv_phyto_diat_N)

   ! Evolution de l azote detritique 
   ! -------------------------------
   dc(iv_detr_N)=-reminazdetsed*cs(iv_detr_N)                          &
                 +diatmortsed*cs(iv_phyto_diat_N)                      &
#if defined key_psnz
                 +txfiltbenthij*cw(iv_phyto_psnz_N)                     &
#endif
#if defined key_karenia
                 +txfiltbenthij*cw(iv_phyto_karenia_N)                    &
#endif
#ifdef key_phaeocystis
                 +txfiltbenthij*cw(iv_phyto_phaeocystis_cell_N)         &
#endif
#ifdef key_ulvas
                 +ulvebenthmortsed*cs(iv_ulv_benth_N)/0.014_rsh         &
#endif
         +filtr_benth

   ! Evolution de la silice particulaire
   ! -------------------------------
   dc(iv_detr_Si)=-dissolsilicesed*cs(iv_detr_Si)                  &
                  +diatmortsed*cs(iv_phyto_diat_N)*p_phyto_SiNratio  &
#if defined key_psnz
                  +txfiltbenthij*cw(iv_phyto_psnz_Si)                     &
#endif
                  +rapsiaz*txfiltbenthij*cw(iv_phyto_diat_N)

   ! Evolution du phosphore detritique
   ! -------------------------------
   dc(iv_detr_P)=-reminpdetsed*cs(iv_detr_P)                            &
                +diatmortsed*cs(iv_phyto_diat_N)*rappaz       &
#if defined key_psnz
                +txfiltbenthij*cw(iv_phyto_psnz_N)*rappaz                   &
#endif
#if defined key_karenia
                +txfiltbenthij*cw(iv_phyto_karenia_C)/(p_phyto_NPratio*p_phyto_CNratio)  &
#endif
#ifdef key_phaeocystis
                +txfiltbenthij*cw(iv_phyto_phaeocystis_cell_N)*rappaz    &
#endif
#ifdef key_ulvas
                +ulvebenthmortsed*cw(iv_ulv_benth_P)/0.031_rsh         &
#endif 
     +filtr_benth*rappaz



#if defined key_oxygen
   ! Evolution de l oxygene dissous 
   ! ------------------------------------------
   dc(iv_oxygen)=-oreminsed*cs(iv_detr_N)/porosite     &
                 -onitrifsed*cs(iv_nutr_NH4)
#endif

#if defined key_zostera
   ! Evolution de l azote detritique peu labile (mmol/m2/j)
   ! ------------------------------------------------------
   dc(iv_detr_zost_N)=-(reminazdetsed/10._rsh)*cs(iv_detr_zost_N)

   ! Evolution du phosphore detritique peu labile (mmol/m2/j)
   ! ------------------------------------------------------
   dc(iv_detr_zost_P)=-(reminazdetsed/10._rsh)*cs(iv_detr_zost_P)
                           
#endif

  END SUBROUTINE incellsed_nut4phy3zoo2

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

  SUBROUTINE incellsed_dinodiatzoo(km,im,jm,c,dc)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellsed_dinodiatzoo  ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes biologiques
   !&E              ayant lieu a l interieur d une maille dans le sediment
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  (A. Menesguen P. Cugier) Original code
   !&E       !  2006-11 (V. Garnier) module - mise en forme
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvarp2d, ONLY  : temp
   USE bioloinit, ONLY : iv_nutr_NO3,iv_nutr_SiOH,iv_nutr_PO4,iv_nutr_Pads, &
                         iv_phyto_diat_N,iv_detr_N,iv_detr_Si,iv_detr_P,     &
                         iv_nutr_NH4,iv_zoo_meso_N,iv_phyto_diat_pp

   !! * Declaration subroutine

   !! * Arguments
   INTEGER, INTENT( in ) :: im,jm,km    ! cell indexes
   REAL(KIND=rsh),DIMENSION(:),INTENT(in) :: c
   REAL(KIND=rsh),DIMENSION(:),INTENT(inout) :: dc

   !! * Local declarations
   INTEGER               :: iv          ! loop index
                     ! variables relative to temperature
   REAL(KIND=rsh)        :: temper,effetchaleur
                     ! variable diatomees
   REAL(KIND=rsh)        :: diatmortsed
                     ! variables sels nutritifs
   REAL(KIND=rsh)        :: rapsiaz,rappaz
                     ! processus mineralisation, nitrification, dissolution
   REAL(KIND=rsh)        :: facto2remin,facto2nitrif,reminazdetsed,     &
                            xnitrifsed,reminpdetsed,dissolsilicesed
                     ! processus adsorption/desorption
   REAL(KIND=rsh)        :: varads,satursedp,adsorpsed,cmaxadsorp,      &
                            desorpsed,cmaxdesorp
                     ! cycle de l oxygene
   REAL(KIND=rsh)        :: oreminsed,onitrifsed

#ifdef key_sedim_anc

   !!----------------------------------------------------------------------
   !! * Executable part

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! saisir le code apres cette limite
   ! =================================
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   temper=temp(km,im,jm)
   IF (temper.gt.30.) THEN
     write(*,*)'tempsed:',km,im,jm,temper
     temper=15.
   END IF
   effetchaleur=exp(p_T_effect*temper)

   ! mortalite des diatomees
   ! -----------------------
   diatmortsed=10.*p_diat_mort*effetchaleur
       
   ! remineralisation dans le sediment
   ! ---------------------------------
   facto2remin=c(iv_nutr_NH4)/(c(iv_nutr_NH4)+(p_micz_mort/2))
   facto2nitrif=c(iv_nutr_NH4)/(c(iv_nutr_NH4)+p_micz_mort)
   reminazdetsed=facto2remin*p_micz_kgraz*effetchaleur
   xnitrifsed=facto2nitrif*bioparam(66)*effetchaleur
   reminpdetsed=facto2remin*p_micz_assim*effetchaleur
   dissolsilicesed=p_micz_excret*effetchaleur
        
   rapsiaz=p_diat_iksmith
   rappaz=1/bioparam(27)

   ! Adsorption et desorption du phosphore dans le sediment
   ! ------------------------------------------------------

   ! adsorption
   ! ----------
   IF (c(iv_phyto_diat_pp).ne.0.0) THEN
     varads=min(1.0,c(iv_nutr_Pads)/(p_P_adsormaxspim*c(iv_phyto_diat_pp)))
   ELSE
     varads=0.0
   END IF
   satursedp=max(0.,(p_P_adsormaxspim*c(iv_phyto_diat_pp)-c(iv_nutr_Pads)))
   adsorpsed=p_P_adsor*satursedp               
   cmaxadsorp=adsorpsed*c(iv_nutr_PO4)*dtbio/86400.
   IF (cmaxadsorp.ge.c(iv_nutr_PO4))  adsorpsed=0.9*86400./dtbio                   
   ! desorption
   ! ----------
   desorpsed=p_P_desor*varads
   cmaxdesorp=desorpsed*c(iv_nutr_Pads)*dtbio/86400.
   IF (cmaxdesorp.ge.c(iv_nutr_Pads)) desorpsed=0.9*86400./dtbio

   ! production et consommation d oxygene dans le sediment
   ! -----------------------------------------------------

   ! remineralisation
   oreminsed=reminazdetsed*c(iv_detr_N)*(0.212+0.064)      
   onitrifsed=xnitrifsed*c(iv_zoo_meso_N)*0.064

   ! Evolution du nitrate
   ! -------------------------------
   dc(iv_nutr_NO3)=xnitrifsed*c(iv_zoo_meso_N)

   ! Evolution de la silice dissoute
   ! -------------------------------
   dc(iv_nutr_SiOH)=dissolsilicesed*c(iv_detr_Si)/poro(km,im,jm)

   ! Evolution du phosphate dissous
   ! -------------------------------
   dc(iv_nutr_PO4)=-adsorpsed*c(iv_nutr_PO4)+desorpsed*c(iv_nutr_Pads)/poro(km,im,jm)             &
             +reminpdetsed*c(iv_detr_P)/poro(km,im,jm)

   ! Evolution du phosphate adsorbe
   ! -------------------------------
   dc(iv_nutr_Pads)=adsorpsed*c(iv_nutr_PO4)*poro(km,im,jm)-desorpsed*c(iv_nutr_Pads)

   ! Evolution de l azote des diatomees
   ! -------------------------------
   dc(iv_phyto_diat_N)=-diatmortsed*c(iv_phyto_diat_N)

   ! Evolution de l azote detritique 
   ! -------------------------------
   dc(iv_detr_N)=-reminazdetsed*c(iv_detr_N)+diatmortsed*c(iv_phyto_diat_N)

   ! Evolution de la silice particulaire
   ! -------------------------------
   dc(iv_detr_Si)=-dissolsilicesed*c(iv_detr_Si)+rapsiaz*diatmortsed*c(iv_phyto_diat_N)

   ! Evolution du phosphore detritique
   ! -------------------------------
   dc(iv_detr_P)=rappaz*diatmortsed*c(iv_phyto_diat_N)-reminpdetsed*c(iv_detr_P)

   ! Evolution de l oxygene dissous
   ! -------------------------------
   dc(iv_nutr_NH4)=-oreminsed/poro(km,im,jm)-onitrifsed

   ! Evolution de l ammonium
   ! -------------------------------
   dc(iv_zoo_meso_N)=reminazdetsed*c(iv_detr_N)/poro(km,im,jm)-xnitrifsed*c(iv_zoo_meso_N)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
   ! saisir le code avant cette limite
   ! =================================
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#else
        write(*,*)  'incellsed_dinodiatzoo : You should not enter this routine'
#endif

  END SUBROUTINE incellsed_dinodiatzoo

#endif
 END MODULE
