 MODULE biolosinksource

#if defined key_substance && defined key_biolo

   !!======================================================================
   !!                   ***  MODULE  biolosinksource ***
   !! Biologic dynamics:  - estimate sink and source terms for each variables
   !!                     - update a few variable
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters,   ONLY : rsh,limin,limax,ljmin,ljmax,ljminp2,ljmaxm1, &
                            liminm1,limaxp2,ljminm1,ljmaxp2,limaxp1,ljmaxp1, &
                            liminm2,ljminm2

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
#if defined key_sedim_MUSTANG  && defined key_biolo
   PUBLIC biolo_sksc_sed            ! routine called by stepu.F90 stepv.F90
#endif
   PUBLIC biolo_sksc_wat            ! routine called by stepu.F90 stepv.F90

   !! * Shared module variables
   REAL(kind=rsh) :: txfiltbenth

   !! * Private variables

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

#if defined key_sedim_MUSTANG  && defined key_biolo
  SUBROUTINE biolo_sksc_sed(xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE biolo_sksc_sed  ***
   !&E
   !&E ** Purpose : Estimate transformations inside sediments
   !&E              Estimate new concentrations in sediments (cv_sed)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          (A. Menesguen) Original code
   !&E       !  2006-11 (V. Garnier) module - mise en forme
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,   ONLY : jmin,jmax,kmax,nb_var,valmanq,rlg
   USE comvars2d,    ONLY : ig,id,h0,t,hm,h0fond
   USE comvars3d,    ONLY : dsigu
   USE comsubstance, ONLY : nv_state,nvpc,dtbio,cv_wat,nv_adv,unit_modif_mudbio_N2dw,irk_fil,dcdt
   USE biolodyncellsed,ONLY : biolo_dyn_cellsed
   USE bioloinit,     ONLY : l_filtbenthsinusmes,p_txfiltbenthmax
   USE comMUSTANG,    ONLY : cv_sed,ksmi,ksma,dzs,poro,c_sedtot,ros,dzsmin
   USE toolgeom,    ONLY : f_dzu
#if ! defined key_biolo_opt2
   USE comvars2d,    ONLY : jjulien
#endif
   USE_MPI toolmpi,      ONLY : ex_i_rsh,ex_j_rsh

   !! * Declaration Subroutine
    REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),     INTENT(in)    :: xe

   !! * Arguments

   !! * Local declarations
   INTEGER                :: i,j,k,iv
   REAL(KIND=rsh),DIMENSION(nv_state)  :: cw,dc,dcw_filtbent
   REAL(KIND=rsh),DIMENSION(-1:nv_state)  :: cs
   REAL(kind=rsh) :: epn_bottom,cvolp
   REAL(kind=rsh) :: foncsinus,d


   !!----------------------------------------------------------------------
   !! * Executable part
   
#if ! defined key_biolo_opt2
   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))
   ELSE
     txfiltbenth=p_txfiltbenthmax
   ENDIF
#else
     txfiltbenth=0.0_rsh
#endif

!$OMP PARALLEL DEFAULT(none) shared(ig,id,h0,ksmi,ksma,dzs,cv_sed     &
!$OMP& ,txfiltbenth,cv_wat,unit_modif_mudbio_N2dw,irk_fil,xe,dsigu    &
!$OMP& ,dzsmin,c_sedtot,nvpc,ros,poro,t                               &
!$OMP& ,nv_state,nv_adv,dtbio,limin,limax,ljminp2,ljmaxm1,dcdt,hm,h0fond)    &
!$OMP& private(i,j,k,cs,cw,dc,dcw_filtbent,iv,epn_bottom,d,cvolp )

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       d=h0(i,j)+xe(i,j)
       IF (h0(i,j).GT.-valmanq)  THEN
         IF (d < hm .AND. d> h0fond) THEN
           epn_bottom=d
           DO iv=1,nv_state
              cw(iv)=SUM(cv_wat(iv,:,i,j)*dsigu(:))
              cw(iv)=0.5_rsh*(cw(iv)+ABS(cw(iv)))
           ENDDO
         ELSE IF (d > hm) THEN
           epn_bottom=f_dzu(h0(i,j),xe(i,j),1,i,j)
           cw(1:nv_state)=0.5_rsh*( cv_wat(1:nv_state,1,i,j)+   &
                              ABS(cv_wat(1:nv_state,1,i,j)) ) 
         ELSE
           cw(1:nv_state)=0.0_rsh
         ENDIF
         cw(1:nv_state)= cw(1:nv_state)/unit_modif_mudbio_N2dw(irk_fil(1:nv_state))                            
         
         DO k=ksmi(i,j),ksma(i,j)
           ! no bio transformation if very small layer
           IF (dzs(k,i,j).GT.dzsmin) THEN
             ! initialize to 0
             dc(1:nv_state)=0.0_rsh
             dcw_filtbent(1:nv_state)=0.0_rsh
             ! transfer 3D var cv_sed into 1D var c
             cs(-1:nv_state)=0.5_rsh*( cv_sed(-1:nv_state,k,i,j)+   &
                                  ABS(cv_sed(-1:nv_state,k,i,j)) )
             cs(1:nv_state)= cs(1:nv_state)/unit_modif_mudbio_N2dw(irk_fil(1:nv_state))                            

             CALL biolo_dyn_cellsed(k,i,j,txfiltbenth,epn_bottom,cw,cs,dc,dcw_filtbent)
             ! update cv_sed and transfer from 1D to 3D 
             cv_sed(1:nv_adv,k,i,j)=cv_sed(1:nv_adv,k,i,j)+  &
               dtbio/86400.0_rsh*dc(1:nv_adv)*unit_modif_mudbio_N2dw(irk_fil(1:nv_adv))
             ! if variable bio = var.part.constitutiv ==> update poro
             IF(SUM(unit_modif_mudbio_N2dw(irk_fil(1:nv_adv))) .NE. nv_adv) THEN
                c_sedtot(k,i,j)=0.0_rsh
                cvolp=0.0_rsh
                DO iv=1,nvpc
                  c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                  cvolp=cvolp+cv_sed(iv,k,i,j)/ros(iv)
                ENDDO
                poro(k,i,j)=1.0_rsh-cvolp
                IF(poro(k,i,j) .LE. 0.0_rsh) THEN
                  write(*,*)'biolosed poro<0 !! ',t,i,j,k,poro(k,i,j),cv_sed(1:nvpc,k,i,j)
                ENDIF
             ENDIF
             IF(k==ksma(i,j) .AND. d > h0fond) THEN
               ! filtrage benthique : phytos dans la couche du fond broutes
               !dcdt(1:nv_state,1,i,j)=dcdt(1:nv_state,1,i,j)+dcw_filtbent(1:nv_state)/86400.0_rsh
               IF (d < hm) THEN
                  cv_wat(1:nv_adv,1,i,j)=cv_wat(1:nv_adv,1,i,j)+  &
                      dcw_filtbent(1:nv_adv)*dtbio/86400.0_rsh/d*unit_modif_mudbio_N2dw(irk_fil(1:nv_adv))  
                  DO iv=1,nv_adv
                     cv_wat(iv,:,i,j)=cv_wat(iv,1,i,j)
                  ENDDO
               ELSE
                  cv_wat(1:nv_adv,1,i,j)=cv_wat(1:nv_adv,1,i,j)+  &
                      dcw_filtbent(1:nv_adv)*dtbio/86400.0_rsh*unit_modif_mudbio_N2dw(irk_fil(1:nv_adv))  
               ENDIF  
             ENDIF
           ENDIF
         END DO
       END IF
     END DO
   END DO
!$OMP END DO

!$OMP END PARALLEL

#if ! defined key_biolo_opt2
  IF(txfiltbenth .NE. 0.0_rsh) THEN
     CALL_MPI ex_i_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_j_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
  ENDIF
#endif
   
  END SUBROUTINE biolo_sksc_sed
#endif

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

  SUBROUTINE biolo_sksc_wat(cv_wat,sal,temp,xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE biolo_sksc_wat  ***
   !&E
   !&E ** Purpose : Estimate variations of concentrations in the water.
   !&E              Estimate new concentrations of fixed variables.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          (A. Menesguen) Original code
   !&E       !  2006-11 (V. Garnier) module - mise en forme
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,   ONLY : kmax,rsh,nb_var
   USE comvars2d,    ONLY : h0,ig,id,h0fond,hm
   USE comvars3d,    ONLY : dsigu
#if ! defined key_biolo_opt2
   USE comvars2d,    ONLY : jjulien
#endif
   USE comsubstance, ONLY : nv_adv,nv_state,dcdt,dtbio,unit_modif_mudbio_N2dw,irk_fil
   USE bioloinit,    ONLY : l_filtbenthsinusmes,p_txfiltbenthmax
   USE biolodyncellwat,ONLY : biolo_dyn_cellwat
   USE_MPI toolmpi,      ONLY : ex_i_rsh,ex_j_rsh

   !! * Declaration Subroutine

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: sal,temp
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_wat

   !! * Local declarations
   INTEGER                :: i,j,k,iv
   REAL(KIND=rsh),DIMENSION(nv_state)  :: c,dc
   REAL(kind=rsh) :: foncsinus

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

#if defined key_sedim_MUSTANG && ! defined key_sedim_partial
   txfiltbenth=0.0_rsh
#elif ! defined key_biolo_opt2
   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))
   ELSE
     txfiltbenth=p_txfiltbenthmax
   ENDIF
#else
  txfiltbenth=0.0_rsh
#endif

!$OMP PARALLEL DEFAULT(none) shared(ig,id,h0,xe,h0fond,cv_wat   &
!$OMP& ,txfiltbenth,unit_modif_mudbio_N2dw,irk_fil,hm,dsigu     &
!$OMP& ,nv_state,nv_adv,dcdt,dtbio,limin,limax,ljminp2,ljmaxm1) &
!$OMP& private(i,j,k,c,dc )

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       IF (h0(i,j)+xe(i,j).GT.h0fond)  THEN
         IF(h0(i,j)+xe(i,j) < hm) THEN
             ! initialize to 0
             dc(1:nv_state)=0.0_rsh
             ! 3D to 1D for variable cv_wat(nb_var,kmax,limin:limax,ljmin:ljmax)
             DO iv=1,nv_adv
              c(iv)=SUM(cv_wat(iv,:,i,j)*dsigu(:))
              c(iv)=0.5_rsh*(c(iv)+ABS(c(iv)))
             ENDDO
             c(nv_adv+1:nv_state)=0.5_rsh*( cv_wat(nv_adv+1:nv_state,1,i,j)+   &
                                    ABS(cv_wat(nv_adv+1:nv_state,1,i,j)) )
             c(1:nv_state)=c(1:nv_state)/unit_modif_mudbio_N2dw(irk_fil(1:nv_state))

             CALL biolo_dyn_cellwat(1,i,j,xe(i,j),txfiltbenth,c,dc)
             ! 1D to 3D or variable dcdt(nv_state,kmax,limin:limax,ljmin:ljmax)
             dcdt(1:nv_state,1,i,j)=dc(1:nv_state)/86400.0_rsh  &
                                 *unit_modif_mudbio_N2dw(irk_fil(1:nv_state))
             dcdt(1:nv_state,2:kmax,i,j)=0.0_rsh
             ! evolution des variables fixees :
             cv_wat(nv_adv+1:nv_state,1,i,j)=cv_wat(nv_adv+1:nv_state,1,i,j) &
                                      +dtbio*dcdt(nv_adv+1:nv_state,1,i,j)
 
          ELSE
           DO k=1,kmax
             ! initialize to 0
             dc(1:nv_state)=0.0_rsh
             ! 3D to 1D for variable cv_wat(nb_var,kmax,limin:limax,ljmin:ljmax)
             c(1:nv_state)=0.5_rsh*( cv_wat(1:nv_state,k,i,j)+   &
                                    ABS(cv_wat(1:nv_state,k,i,j)) )
             c(1:nv_state)=c(1:nv_state)/unit_modif_mudbio_N2dw(irk_fil(1:nv_state))

             CALL biolo_dyn_cellwat(k,i,j,xe(i,j),txfiltbenth,c,dc)
             ! 1D to 3D or variable dcdt(nv_state,kmax,limin:limax,ljmin:ljmax)
             dcdt(1:nv_state,k,i,j)=dc(1:nv_state)/86400.0_rsh  &
                                 *unit_modif_mudbio_N2dw(irk_fil(1:nv_state))
             ! evolution des variables fixees :
             cv_wat(nv_adv+1:nv_state,k,i,j)=cv_wat(nv_adv+1:nv_state,k,i,j) &
                                      +dtbio*dcdt(nv_adv+1:nv_state,k,i,j)
           END DO
         ENDIF
       END IF
     END DO
   END DO
!$OMP END DO

!$OMP END PARALLEL
   IF(nv_state > nv_adv) THEN
     CALL_MPI ex_i_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_j_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
   ENDIF

  END SUBROUTINE biolo_sksc_wat
#endif

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

 END MODULE
