 MODULE contasinksource

#if defined key_substance && defined key_contaminant

   !!======================================================================
   !!                   ***  MODULE  contasinksource ***
   !!                 Chimestry - Contaminants dynamics:  - 
   !!            estimate sink and source terms for each variables
   !!======================================================================
   !! * Modules used
   USE parameters,   ONLY : rsh,rlg,kmax,imin,imax,jmin,jmax,limin,limax,        &
                            ljmin,ljmax,liminm1,limaxp2,ljminm1,ljmaxp2,ljminp2, &
                            limaxp2,ljmaxp2,ljmaxm1,liminm2,ljminm2

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
!   PUBLIC conta_sksc_sed            ! routine called by stepu.F90 stepv.F90
   PUBLIC conta_sksc_wat                ! routine called by stepu.F90 stepv.F90
   PUBLIC conta_reac_after              ! routine called by step.F90 
#ifdef key_species
   PUBLIC conta_fspc_wat               ! routine called by stepu.F90 stepv.F90
#endif
          
   !! * Shared module variables        
  
   !! * Private variables

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

   SUBROUTINE conta_reac_after(cv_watloc,sali,tempe,xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_reac_after  ***
   !&E
   !&E ** Purpose : fast reactions or at equilibrium
   !&E              modification of concentrations after transport (splitting) 
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : cv_wat,sal,temp,xe,h0,flux_reac_wat,surf,iarea_balanc
   !&E                   if key_species : dc_spc_wat
   !&E
   !&E ** modified variables : cv_wat,flux_reac_wat
   !&E                   if key_species : dc_spc_wat
   !&E
   !&E ** Called by : stepu, stepv
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06        (B. Thouvenin) Original code 
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : h0,ig,id,h0fond,t,dt,surf,fwetp
   USE toolgeom,     ONLY : f_dzu
   USE comvars3d,    ONLY : dsigu
#ifdef key_species
   USE comsubstance, ONLY : nv_state,nv_tot                          
!   USE comsubstance, ONLY : nv_state,nv_tot,l_output_balpv                          
   USE contareaddat, ONLY : dtconta_dc,nreacmax,l_treat_re,l_conta,l_gpspecies  &
                            nb_var_re,l_sorfluxreac,irkm_var_re
                                                     
#else
!   USE comsubstance, ONLY : nv_state,nv_tot,l_output_balpv
   USE comsubstance, ONLY : nv_state,nv_tot
   USE contareaddat, ONLY : dtconta_dc,nreacmax,l_treat_re,l_conta, &
                            nb_var_re,irkm_var_re,l_sorfluxreac,ityp_reac
#endif
   USE contadyncellwat,ONLY : conta_dyn_cellwat
!   USE contareaddat, ONLY  : flux_reac_wat
   USE submassbalance,  ONLY: flx_reac_zone,nb_zone_budget,mask_bud_in

   !! * Declaration Subroutine

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_watloc
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: sali,tempe
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe

   !! * Local declarations
   INTEGER                                     :: i,j,k,ireac,ivr,iv,iz
   REAL(KIND=rsh),DIMENSION(nv_tot)            :: c
   REAL(KIND=rsh),DIMENSION(nv_state)          :: c_old
   REAL(KIND=rsh)                              :: dzu,voltot

   !!----------------------------------------------------------------------
   !! * Executable part
   !!======================================================================
    dtconta_dc=dt

#ifdef key_species      
     !   evolution at each half time step of contaminant species into advected variables
     IF(l_conta .and. l_gpspecies)  &
             CALL conta_fspc_wat(cv_watloc,xe)
#endif
     !   new concentrations after each half time step of state variables 
     !   which are submitted to rapid reactions (time splitting)
    DO j=ljmin,ljmax
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       IF (h0(i,j)+xe(i,j).gt.h0fond)  THEN
         DO k=1,kmax
         ! 3D to 1D for variable cv_watloc(nv_tot,kmax,imin:imax,jmin:jmax)
           c(1:nv_tot)=0.5_rsh*(cv_watloc(1:nv_tot,k,i,j)+   &
                                    abs(cv_watloc(1:nv_tot,k,i,j)) )
           c_old(1:nv_state)=c(1:nv_state)
           ! contribution of each reaction in c
             DO ireac=1,nreacmax   
             ! Call to generic procedure  conta_dyn_cellwat, without argument 'dc' 
              IF(.not.l_treat_re(ireac)) THEN
               ! si reaction de partage et variable particulaire negative
                IF(ityp_reac(ireac)==5) THEN
                  DO ivr=2,nb_var_re(ireac)
                    iv=irkm_var_re(ireac,ivr)
                    IF(cv_watloc(iv,k,i,j).lt.0.0_rsh)then
                       IF(-cv_watloc(iv,k,i,j).gt.c(irkm_var_re(ireac,1))) then
                         cv_watloc(iv,k,i,j)=c(irkm_var_re(ireac,1))+cv_watloc(iv,k,i,j)
                         c(irkm_var_re(ireac,1))=0.0_rsh
                       ELSE
                         c(irkm_var_re(ireac,1))=c(irkm_var_re(ireac,1))+cv_watloc(iv,k,i,j)
                         cv_watloc(iv,k,i,j)=0.0_rsh
                       ENDIF
                    ENDIF
                 ENDDO 
               ENDIF
                call conta_dyn_cellwat(ireac,k,i,j,  &
                                              sali(k,i,j),tempe(k,i,j),xe(i,j),c)
                !  1D to 3D 
                DO ivr=1,nb_var_re(ireac)
                  iv=irkm_var_re(ireac,ivr)
                  IF(cv_watloc(iv,k,i,j) >= 0.0_rsh) THEN
                     cv_watloc(iv,k,i,j)=c(iv)
                  ELSE IF(c(iv) >0.0_rsh) THEN
                     cv_watloc(iv,k,i,j)=c(iv)+cv_watloc(iv,k,i,j)
                  ENDIF
                ENDDO 
     !           If(j==84 .and. i==44 .and. k==1)then
     !              write(*,*)'c_old',k,c_old(1:nv_state)
     !              write(*,*)'c',k,c(1:nv_state)
     !              write(*,*)'cv',k,cv_watloc(1:nv_state,k,i,j)
     !           endif         
                !    calcul des flux reactifs entre variables
!     !           IF(l_sorfluxreac) THEN
!     !             DO ivr=1,nb_var_re(ireac)
!     !              iv=irkm_var_re(ireac,ivr)
!     !              dzu=f_dzu(h0(i,j),xe(i,j),k,i,j)
!     !              voltot=REAL(dzu*(h0(i,j)+xe(i,j))*surf(i,j)*fwetp(i,j),rlg)
!     !              flux_reac_wat(ireac,ivr,k,i,j)=flux_reac_wat(ireac,ivr,k,i,j)+   &
!     !                                               c(iv)-c_old(iv)
!     !             ENDDO
!     !           ENDIF
                  DO ivr=1,nb_var_re(ireac)
                   iv=irkm_var_re(ireac,ivr)
                   voltot=REAL(f_dzu(h0(i,j),xe(i,j),k,i,j)*surf(i,j)*fwetp(i,j),rlg)
                   DO iz=1,nb_zone_budget
                      flx_reac_zone(iz,ireac,iv)=flx_reac_zone(iz,ireac,iv)+   &
                                   mask_bud_in(iz,i,j)*(c(iv)-c_old(iv))*voltot
                   ENDDO
                  ENDDO
              ENDIF
             ENDDO
         END DO
       END IF
     END DO
   END DO

   END SUBROUTINE conta_reac_after
  !!======================================================================

   SUBROUTINE conta_sksc_wat(cv_watloc,sali,tempe,xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_sksc_wat  ***
   !&E
   !&E ** Purpose : Estimate variations of concentrations in the water dcdt if call before adv
   !&E              Estimate new concentrations after rapid reactions if call after adv.
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : cv_wat,sal,temp,xe,h0,flux_reac_wat,surf,iarea_balanc
   !&E                   if key_species : dc_spc_wat
   !&E
   !&E ** modified variables : cv_wat,flux_reac_wat
   !&E                   if key_species : dc_spc_wat
   !&E
   !&E ** Called by : stepu, stepv
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06        (B. Thouvenin) Original code 
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,   ONLY : lchain
   USE comvars2d,    ONLY : h0,ig,id,h0fond,t,dt,dx,dy,surf,fwetp
   USE toolgeom,     ONLY : f_dzu 
   USE comvars3d,    ONLY : dsigu
   USE comsubstance, ONLY : nv_adv,nv_state,dcdt,irk_fil,irk_mod,nv_tot
!  USE comsubstance, ONLY : l_output_balpv
   USE contareaddat, ONLY : dtconta_dc,nreacmax,ityp_reac,l_treat_re,l_conta, &
                           icon_reac,nb_var_re,irkm_var_re,l_sorfluxreac
!   USE contareaddat, ONLY : tconta_sksc0,tconta_sksc,f_dtconta_sksc
                         
#ifdef key_species
   USE contareaddat, ONLY : ifrac_var_reac,dc_spc_wat,l_gpspecies                          
#endif
   USE contadyncellwat,ONLY : conta_dyn_cellwat
!   USE contareaddat, ONLY  : flux_reac_wat
   USE submassbalance,  ONLY: flx_reac_zone,nb_zone_budget,mask_bud_in

   !! * Declaration Subroutine

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_watloc
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: sali,tempe
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe

   !! * Local declarations
   INTEGER                                     :: i,j,k,ireac,ivr,iv,iz
   REAL(KIND=rsh),DIMENSION(nv_tot)            :: c
   REAL(KIND=rsh),DIMENSION(nv_state)          :: dc
   REAL(KIND=rsh),DIMENSION(nv_state)          :: dc_ire
   real(kind=rsh)                              :: cvwatmax,voltot,dzu

   !!----------------------------------------------------------------------
   !! * Executable part
   dtconta_dc=dt*0.5_rlg

!!!$OMP PARALLEL DEFAULT(none) shared(ig,id,h0,xe,h0fond,cv_watloc,dtconta_dc &
!!!$OMP& ,flux_reac_wat,surf,fwetp,irkm_var_re, &
!!#ifdef key_species
!!$OMP& ,dc_spc_wat,icon_reac,l_gpspecies,ifrac_var_reac &
!!#endif
!!!$OMP& ,nv_tot,nreacmax,l_treat_re,sali,tempe,l_sorfluxreac   &
!!!$OMP& ,nb_var_re,irk_fil,nv_state,nv_adv,irk_mod,dcdt,ljminp2,ljmaxm1)  &
!!!$OMP& private(i,j,k,c,dc,dc_ire,ireac,voltot,ivr,c_old )

!!!$OMP DO SCHEDULE(RUNTIME)
    dcdt(:,:,:,:)=0.0_rsh
    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
         DO k=1,kmax
           dc(1:nv_state)=0.
         ! 3D to 1D for variable cv_watloc(nv_tot,kmax,imin:imax,jmin:jmax)
           c(1:nv_tot)=0.5*( cv_watloc(1:nv_tot,k,i,j)+   &
                                    abs(cv_watloc(1:nv_tot,k,i,j)) )
                   
         ! contribution of each reaction in dc
#ifdef key_species
           dc_spc_wat(:,:,k,i,j)=0.
#endif
           DO ireac=1,nreacmax 
              dc_ire(:)=0.
           ! Call to generic procedure  conta_dyn_cellwat, with argument dc 
              IF(l_treat_re(ireac)) THEN
                call conta_dyn_cellwat(ireac,k,i,j,   &
                                  sali(k,i,j),tempe(k,i,j),xe(i,j),c,dc_ire)
                dc(1:nv_state)=dc(1:nv_state)+dc_ire(1:nv_state)
#ifdef key_species
                IF(icon_reac(ireac) /= 0 .and. l_gpspecies) THEN
                  DO ivr=1,nb_var_re(ireac)
                    iv=irkm_var_re(ireac,ivr)
                    dc_spc_wat(iv,ifrac_var_reac(ireac,ivr),k,i,j)=  &
                        dc_spc_wat(iv,ifrac_var_reac(ireac,ivr),k,i,j)+dc_ire(iv)
                  ENDDO
                ENDIF
#endif
           ! calcul des flux reactifs entre variables
!                IF(l_sorfluxreac) THEN
!                  dzu=f_dzu(h0(i,j),xe(i,j),k,i,j)
!                  voltot=REAL(dzu*(h0(i,j)+xe(i,j))*surf(i,j)*fwetp(i,j),rlg)
!                  DO ivr=1,nb_var_re(ireac)
!                    iv=irkm_var_re(ireac,ivr)
!                    flux_reac_wat(ireac,ivr,k,i,j)=flux_reac_wat(ireac,ivr,k,i,j)+   &
!                                                    dc_ire(iv)*dtconta_dc*voltot
!                  ENDDO
!                ENDIF
                 DO ivr=1,nb_var_re(ireac)
                   iv=irkm_var_re(ireac,ivr)
                   voltot=REAL(f_dzu(h0(i,j),xe(i,j),k,i,j)*surf(i,j)*fwetp(i,j),rlg)
                   DO iz=1,nb_zone_budget
                      flx_reac_zone(iz,ireac,iv)=flx_reac_zone(iz,ireac,iv)+   &
                                 mask_bud_in(iz,i,j)*dc_ire(iv)*dtconta_dc*voltot
                   ENDDO
                 ENDDO
              ENDIF
            ENDDO
            ! 1D to 3D for variable dcdt(nv_state,kmax,imin:imax,jmin:jmax)
            dcdt(1:nv_state,k,i,j)=dc(1:nv_state)
          END DO
       END IF
      END DO
     END DO
    
!!!$OMP END DO

#ifdef key_species
     !  mise a jour de chaque espece de contaminant dans les variables advectees
     IF(l_conta .and.l_gpspecies)dc_spc_wat(1:nv_adv,:,:,:,:)=dc_spc_wat(1:nv_adv,:,:,:,:)*dt*0.5_rlg
#endif
!     tconta_sksc0=t
!     tconta_sksc=t+f_dtconta_sksc*dt

!!!$OMP END PARALLEL
   END SUBROUTINE conta_sksc_wat
   
#if defined key_contaminant && defined key_species
  !!======================================================================

   SUBROUTINE conta_fspc_wat(cv_watloc,xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_fspc_wat  ***
   !&E
   !&E ** Purpose : Estimate variations of species of a contaminant 
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : if key_species : dc_spc_wat,cv_wat,dcspc_tot,h0,xe,frac_cwat
   !&E
   !&E ** modified variables : if key_species : dcspc_tot,frac_cwat
   !&E
   !&E ** Called by : stepu, stepv
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06        (B. Thouvenin) Original code, 
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comsubstance, ONLY : nv_adv
   USE contareaddat, ONLY : nb_spcmax
   USE comsubstance, ONLY : dc_spc_wat,nb_frac_var,frac_cwat
   USE comvars2d,    ONLY : h0,ig,id,h0fond

   !! * Declaration Subroutine

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_watloc
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: xe

   !! * Local declarations
   INTEGER                                               :: iv,ifracv,i,j,k
   REAL(KIND=rsh),DIMENSION(nv_adv,kmax,limin:limax,ljmin:ljmax)  :: dcspc_tot
   REAL(KIND=rsh) ::  cvwatmax

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

   dcspc_tot(:,:,:,:)=0.0_rsh
   DO iv=1,nv_adv
     DO ifracv=1,nb_frac_var(iv)
           dcspc_tot(iv,:,:,:)=dcspc_tot(iv,:,:,:)+dc_spc_wat(iv,ifracv,:,:,:)
     ENDDO
   ENDDO
   DO iv=1,nv_adv
    IF(nb_frac_var(iv) > 1)then
     DO ifracv=1,nb_frac_var(iv)
      DO j=ljminp2,ljmaxm1
       DO i=ig(j)+1,id(j)-1
        IF (h0(i,j)+xe(i,j).gt.h0fond)  THEN
          DO k=1,kmax
            IF(cv_watloc(iv,k,i,j) <= 1.e-10)then
              frac_cwat(iv,:,k,i,j)=0.
            ELSE
              frac_cwat(iv,ifracv,k,i,j)=(dc_spc_wat(iv,ifracv,k,i,j) +  &
                    frac_cwat(iv,ifracv,k,i,j)*(cv_watloc(iv,k,i,j)-dcspc_tot(iv,k,i,j))) &
                    /cv_watloc(iv,k,i,j)
            ENDIF
          END DO
        END IF
       END DO
      END DO
     ENDDO
    ELSE
     frac_cwat(iv,1,:,:,:)=1.
     frac_cwat(iv,2:nb_spc_max,:,:,:)=0.
    ENDIF
   ENDDO
   
   END SUBROUTINE conta_fspc_wat
#endif

#endif
   !!======================================================================
 END MODULE contasinksource
