 MODULE contadyncellwat

#if defined key_substance && defined key_contaminant
   !!======================================================================
   !!                   ***  MODULE  contadyncellwat  ***
   !! Chemistry - Contaminants dynamics: Calcul des variations dues 
   !!                    aux phenomenes physico-, bio- chimiques
   !!                    ayant lieu a l interieur d une maille dans l eau
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters,   ONLY : rsh,rlg,kmax,epsilon
   USE comvars2d,    ONLY : h0
   USE cometemp,     ONLY   : ventvent
   USE comsubstance, ONLY : irk_mod,nv_state,ws3,nv_tot
   USE contareaddat, ONLY : dtconta_dc

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC conta_dyn_cellwat   ! routine called by conta_sksc_wat.F90
   !! * Interface
   INTERFACE conta_dyn_cellwat
     MODULE PROCEDURE conta_dyn_c
     MODULE PROCEDURE conta_dyn_dc
   END INTERFACE

   !! * Shared module variables

   !! * Private variables

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

  SUBROUTINE conta_dyn_dc(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dc_ire)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  conta_dyn_dc ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes PHY-BIO-CHI
   !&E              ayant lieu a l interieur d une maille dans l eau
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : 
   !&E
   !&E ** modified variables : 
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE contareaddat, ONLY : ityp_reac 
   
   !! * Arguments
   INTEGER, INTENT( in ) ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in ) ::   ire           ! num of reaction
   REAL(KIND=rsh),             INTENT(in) :: xeta     !ssh
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(in) :: c
   REAL(KIND=rsh),DIMENSION(nv_state),INTENT(inout) :: dc_ire
   REAL(KIND=rsh),INTENT(in)              :: sal_ijk,temp_ijk

   !! * Local declarations
   INTEGER                                :: itypr

   ! indices de boucles
   ! ------------------
   INTEGER                             :: iv

   ! profondeur colonne d eau
   ! ------------------------
   REAL(KIND=rsh)                       :: d


 
   !!----------------------------------------------------------------------
   !! * Executable part
 
   ! source and sink in the cell
   !  depend on the type of the reaction
   itypr=ityp_reac(ire)
   SELECT CASE (itypr)
    CASE (1)
      CALL incellwat_REAC1(ire,km,im,jm,temp_ijk,c,dc_ire)
    CASE (2)
      CALL incellwat_GEN(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dc_ire)
    CASE(3)
      CALL incellwat_REVERS(ire,km,im,jm,temp_ijk,xeta,c,dc_ire)
    CASE(4)
      CALL incellwat_VOLAT(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dc_ire)
   END SELECT

  END SUBROUTINE conta_dyn_dc

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

  SUBROUTINE conta_dyn_c(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  conta_dyn_cellwat ***
   !&E
   !&E ** Purpose : Calcul des nouvelles concentrations apres variations 
   !&E              dues aux phenomenes PHY-BIO-CHI rapides ou instantanees (supposes toujours a l equilibre)
   !&E              ayant lieu a l interieur d une maille dans l eau
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE contareaddat, ONLY : ityp_reac
#ifdef key_species
   USE contareaddat, ONLY : frac_cwat,ifrac_var_reac,icon_reac,irkm_var_re,nb_var_re
#endif

   !! * Arguments
   INTEGER, INTENT( in ) ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in ) ::   ire           ! num of reaction
   REAL(KIND=rsh),             INTENT(in) :: xeta     !ssh
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(inout) :: c
   REAL(KIND=rsh),INTENT(in)                :: sal_ijk,temp_ijk

   !! * Local declarations
   REAL(KIND=rsh),DIMENSION(nv_state)         :: dcc,c_anc
   INTEGER                                    :: itypr

   ! indices de boucles
   ! ------------------
   INTEGER                             :: iv,ivr

   ! profondeur colonne d eau
   ! ------------------------
   REAL(KIND=rsh)                      :: d


 
   !!----------------------------------------------------------------------
   !! * Executable part
 
   ! source and sink in the cell
   ! depend on the type of the reaction
   itypr=ityp_reac(ire)
   c_anc=c
   SELECT CASE (itypr)
    CASE (1)
      dcc(:)=0.
      CALL incellwat_REAC1(ire,km,im,jm,temp_ijk,c,dcc)
      c(1:nv_state)=c(1:nv_state)+dcc(1:nv_state)*dtconta_dc
#ifdef key_species
      IF (icon_reac(ire) /= 0) THEN
        iv=irkm_var_re(ire,1)
        frac_cwat(iv,ifrac_var_reac(ire,1),km,im,jm)=(frac_cwat(iv,ifrac_var_reac(ire,1),km,im,jm)+dcc(iv)*dtconta_dc)/(c(iv)+1.e-10)
      END IF
#endif
    CASE (2)
      dcc(:)=0.
      CALL incellwat_GEN(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dcc)
      c(1:nv_state)=c(1:nv_state)+dcc(1:nv_state)*dtconta_dc
#ifdef key_species
      IF (icon_reac(ire) /= 0) THEN
        DO ivr=1,nb_var_re(ire)
          iv=irkm_var_re(ire,ivr)
          frac_cwat(iv,ifrac_var_reac(ire,ivr),km,im,jm)=(frac_cwat(iv,ifrac_var_reac(ire,ivr),km,im,jm)+dcc(iv)*dtconta_dc)/(c(iv)+1.e-10)
        END DO
      END IF
#endif
    CASE(3)
       dcc(:)=0.
       CALL incellwat_REVERS(ire,km,im,jm,temp_ijk,xeta,c,dcc)
       c(1:nv_state)=c(1:nv_state)+dcc(1:nv_state)*dtconta_dc
#ifdef key_species
      IF (icon_reac(ire) /= 0) THEN
        DO ivr=1,nb_var_re(ire)
          iv=irkm_var_re(ire,ivr)
          frac_cwat(iv,ifrac_var_reac(ire,ivr),km,im,jm)=(frac_cwat(iv,ifrac_var_reac(ire,ivr),km,im,jm)+dcc(iv)*dtconta_dc)/(c(iv)+1.e-10)
        END DO
      END IF
#endif
    CASE(4)
       dcc(:)=0.0_rsh
       CALL incellwat_VOLAT(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dcc)
       c(1:nv_state)=c(1:nv_state)+dcc(1:nv_state)*dtconta_dc
#ifdef key_species
       IF (icon_reac(ire) /= 0) THEN
          iv=irkm_var_re(ire,1)
          frac_cwat(iv,ifrac_var_reac(ire,1),km,im,jm)=(frac_cwat(iv,ifrac_var_reac(ire,1),km,im,jm)+dcc(iv)*dtconta_dc)/(c(iv)+1.e-10)
       END IF
#endif
    CASE(5)
#ifdef key_species
      IF (icon_reac(ire) /= 0 ) THEN
        CALL incellwat_PART_conta(ire,km,im,jm,sal_ijk,temp_ijk,c)
      END IF
#else
#ifdef key_menor
        CALL incellwat_PART_medit(ire,km,im,jm,sal_ijk,temp_ijk,c)
     !!   CALL incellwat_PARTCIN(ire,km,im,jm,sal_ijk,temp_ijk,c)
#else
        CALL incellwat_PART(ire,km,im,jm,sal_ijk,temp_ijk,c)
#endif
#endif
   END SELECT       

  END SUBROUTINE conta_dyn_c

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

  SUBROUTINE incellwat_REAC1(ire,km,im,jm,temp_ijk,c,dc_ire)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_REAC1  ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type REAC_1 ayant lieu a l interieur d une maille dans l eau
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE contareaddat, ONLY : cin_max_re,l_arh,cf_stochio_re, tempref_arh,   &
                            irkm_var_re,cf_teta_arh,ityp_arh,cf_xetsig_arh
#ifdef key_species
   USE contareaddat, ONLY : frac_cwat,ifrac_var_reac,icon_reac
#endif

   !! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
   REAL(KIND=rsh),             INTENT(in) :: temp_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(in) :: c
   REAL(KIND=rsh),DIMENSION(nv_state),INTENT(inout) :: dc_ire

   !! * Local declarations
   REAL(KIND=rsh)             :: cin_reac,xdtemp,fspc
   INTEGER                    :: iv      

! correction de la cinetique avec la temperature
   cin_reac=cin_max_re(ire)
   IF (l_arh(ire)) then
     xdtemp=(temp_ijk-tempref_arh(ire))/cf_xetsig_arh(ire)
     cin_reac=cin_max_re(ire)*((1-ityp_arh(ire))*cf_teta_arh(ire)**xdtemp &
                              +ityp_arh(ire)*cf_teta_arh(ire)*exp(-xdtemp**2))
   END IF

! calcul du terme soupui
    iv=irkm_var_re(ire,1)
    fspc=1.
#ifdef key_species
    IF (icon_reac(ire) /= 0) THEN
      fspc= frac_cwat(iv,ifrac_var_reac(ire,1),km,im,jm)
    END IF  
#endif
    dc_ire(iv)=cf_stochio_re(ire,1)*cin_reac*fspc*c(irkm_var_re(ire,1)) 
    
  END SUBROUTINE incellwat_REAC1

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

  SUBROUTINE incellwat_GEN(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dc_ire)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_GEN  ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type GENERIC ayant lieu a l interieur d une maille dans l eau
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE userconta, ONLY : userconta_coef_re2,userconta_cin_re2
   USE comsubstance, ONLY : alumplafond
   USE contareaddat, ONLY : l_cincte_re,cin_max_re,tempref_arh, &
                            cf_xetsig_arh, ityp_arh,cf_teta_arh,   &
                            nb_flim1_re,irkm_flim1_re,cf_flim1_re,   &
                            nb_flim2_re,irkm_flim2_re,cf_flim2_re,   &
                            cf_expovar_re,cf_stochio_re,irkm_var_re,   &
                            l_coefcte_re,nb_var_re,  &
                            cf_flim3_re,nb_flim3_re,l_arh

#ifdef key_species
   USE contareaddat, ONLY : frac_cwat,ifrac_var_reac,icon_reac
#endif

   !! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
   REAL(KIND=rsh),             INTENT(in) :: xeta     !ssh
   REAL(KIND=rsh),INTENT(in)              :: sal_ijk,temp_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(in) :: c
   REAL(KIND=rsh),DIMENSION(nv_state),INTENT(inout) :: dc_ire

   !! * Local declarations
   REAL(KIND=rsh)                     :: cin_reac,xdtemp,xflim1,xflim2,xflim3
   REAL(KIND=rsh)                     :: xcv,coext
   INTEGER                            :: ivr,iv,iflim,ivfl ,ire2,k    
   REAL(KIND=rsh),DIMENSION(nv_tot)   ::   coefstochio,fspc
   
   ! profondeur colonne d eau
   ! ------------------------
   REAL(KIND=rsh) :: d
   
! calcul de la cinetique si elle varie en fonction d autres elements (sauf temperature)
   cin_reac=cin_max_re(ire)
   IF (.not.l_cincte_re(ire)) then
       CALL userconta_cin_re2(ire,km,im,jm,cin_reac)
   END IF
   
! correction de la cinetique avec la temperature
   IF (l_arh(ire)) then
    xdtemp=(temp_ijk-tempref_arh(ire))/cf_xetsig_arh(ire)
    cin_reac=cin_reac*((1-ityp_arh(ire))*cf_teta_arh(ire)**xdtemp &
                     +ityp_arh(ire)*cf_teta_arh(ire)*exp(-xdtemp**2))
   END IF

! fonctions limitantes
   xflim1=1.0
   DO iflim=1,nb_flim1_re(ire)
     ivfl=irkm_flim1_re(ire,iflim)
     IF (ivfl.eq.-1) THEN
       xflim1=xflim1*abs(sal_ijk/(sal_ijk+cf_flim1_re(ire,iflim)))
     ELSE
       xflim1=xflim1*abs(c(ivfl)/(c(ivfl)+cf_flim1_re(ire,iflim)))
     END IF
   END DO
   xflim2=1.0
   DO iflim=1,nb_flim2_re(ire)
     ivfl=irkm_flim2_re(ire,iflim)
     IF (ivfl.eq.-1) THEN
       xflim2=xflim2*abs(cf_flim2_re(ire,iflim)/(sal_ijk+cf_flim2_re(ire,iflim)))
     ELSE
       xflim2=xflim2*abs(cf_flim2_re(ire,iflim)/(c(ivfl)+cf_flim2_re(ire,iflim)))
     END IF
   END DO
   xflim3=1.0
   IF (nb_flim3_re(ire) > 0) then
     DO iflim=1,nb_flim3_re(ire)
! attention _ a revoir
       xflim3=max(1.e-10,xflim3*alumplafond(km,im,jm)**cf_flim3_re(ire,iflim))
     END DO
   END IF
   
! calcul du coefficient stochio si il varie en fonction d autres elements
!*********************************
   coefstochio(1:size(cf_stochio_re,2))=cf_stochio_re(ire,:)
   IF (.not.l_coefcte_re(ire)) then
     CALL userconta_coef_re2(ire,coefstochio)
   END IF
   
   xcv=1.0

#ifdef key_species
    IF (icon_reac(ire) /= 0) THEN
      ! si la reaction agit sur un contaminant et les especes sont groupees
     DO ivr=1,nb_var_re(ire)
         iv=irkm_var_re(ire,ivr)
         fspc(ivr)=frac_cwat(iv,ifrac_var_reac(ire,ivr),km,im,jm)
     END DO
    ELSE
      fspc(:)=1.
    END IF
#else
    fspc(:)=1.
#endif

! concentrations actives dans le flux d echange (si exposant different de 0)
    DO ivr=1,nb_var_re(ire)
      iv=irkm_var_re(ire,ivr)
      xcv=xcv*(fspc(ivr)*c(iv))**cf_expovar_re(ire,ivr)
    END DO
! calcul du terme soupui quand le terme xcv=C^n.C^m est calcule
    DO ivr=1,nb_var_re(ire)
      iv=irkm_var_re(ire,ivr)
      dc_ire(iv)=coefstochio(ivr)*cin_reac*xflim1*xflim2*xflim3*xcv
    END DO
    
  END SUBROUTINE incellwat_GEN
!
   !!======================================================================

  SUBROUTINE incellwat_REVERS(ire,km,im,jm,temp_ijk,xeta,c,dc_ire)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_REVERS  ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type REVERSIBLE ayant lieu a l interieur d une maille dans l eau
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : dt
   USE comsubstance, ONLY : irkm_var_assoc,itypv
   USE contareaddat, ONLY : cf_rate_re3,tempref_arh, ire_typ, &
                            cf_xetsig_arh, ityp_arh,cf_teta_arh, &
                            irkm_var_re,cv_watequ_re3,l_arh

   !! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
   REAL(KIND=rsh),             INTENT(in) :: xeta     
   REAL(KIND=rsh),INTENT(in)              :: temp_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(in) :: c
   REAL(KIND=rsh),DIMENSION(nv_state),INTENT(inout) :: dc_ire

   !! * Local declarations
   REAL(KIND=rsh)                     :: rate,xdtemp,echg,cequi
   INTEGER                            :: iv ,ire3    
   
   ire3=ire_typ(ire)
   rate=cf_rate_re3(ire3)
! correction de la constante de vitesse avec la temperature
   IF (l_arh(ire)) then
    xdtemp=(temp_ijk-tempref_arh(ire))/cf_xetsig_arh(ire)
    rate=rate*((1-ityp_arh(ire))*cf_teta_arh(ire)**xdtemp &
                     +ityp_arh(ire)*cf_teta_arh(ire)*exp(-xdtemp**2))
   END IF
   
   iv=irkm_var_re(ire,1)
   IF(itypv(iv) < 6) THEN
      cequi=MAX(0._rsh,cv_watequ_re3(ire3)*c(irkm_var_assoc(iv)))
   ELSE
      cequi=cv_watequ_re3(ire3)
   END IF
   echg=rate*(c(iv)-cequi)
! Pb du traitement des valeurs negatives
   if(echg*dt < -c(irkm_var_re(ire,2))) THEN
     echg=0.9*c(irkm_var_re(ire,2))/dt
   endif
   IF((c(iv)<=0.0_rsh .AND. echg > 0.0_rsh) .OR. &
    (c(irkm_var_re(ire,2))<=0.0_rsh .AND. echg < 0.0_rsh)  .OR. &
             echg <= epsilon) echg = 0.0_rsh
   dc_ire(iv)=-echg
   dc_ire(irkm_var_re(ire,2))=+echg
 
   
  END SUBROUTINE incellwat_REVERS
!
   !!======================================================================

  SUBROUTINE incellwat_VOLAT(ire,km,im,jm,sal_ijk,temp_ijk,xeta,c,dc_ire)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_VOLAT  ***
   !&E
   !&E ** Purpose : Calcul des variations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type VOLATIL ayant lieu a l interieur d une maille dans l eau
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-09  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : h0
   USE comvarp3d,    ONLY : uz,vz
   USE cometemp,     ONLY : ventvent
   USE comsubstance, ONLY : cini_air
   USE contareaddat, ONLY : irkm_var_re,ire_typ,cf_henry_adim,cf_henry,cst_gas_univ,  &
                            l_volkg_cor,cf_difair,massmol,massmol_h2o,cf_expvol_gas,  &
                            cf_difair_h2o,l_volkl_dif,l_volkl_sal,cf_difwat_co2,cf_difwat_o2,   &
                            ich_volkl_wind,l_volkl_cor,nschmidt,massmol_co2,  &
                            cf_expvol_liq,cf_temp_henry,cf_salt_henry,cf_difwat,thick_volat
                            
   USE toolgeom,        ONLY : f_dzu 
!! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
   REAL(KIND=rsh),             INTENT(in) :: xeta     
   REAL(KIND=rsh),INTENT(in)              :: temp_ijk,sal_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(in) :: c
   REAL(KIND=rsh),DIMENSION(nv_state),INTENT(inout) :: dc_ire

   !! * Local declarations
   REAL(KIND=rsh)                     :: temperaKi,d,uci,cf_henry_adimi,xkgh2o,sc_conta
   REAL(KIND=rsh)                     :: xkg,xkldo2,xkldco2,xklvco2,xkl,sc_co2,xkol,epn,dzu
   REAL(KIND=rsh)                     :: ztransf
   INTEGER                            :: iv,ire4,k   

!  profondeur a partir de la surface
   d=h0(im,jm)+xeta
   epn=0.
   DO k=kmax,km,-1
     dzu=f_dzu(h0(im,jm),xeta,k,im,jm)
     epn=epn+dzu
   END DO
   ztransf=min(d,thick_volat)
! volatilisation sur une epaisseur constante en surface
   IF (epn < thick_volat)THEN 
     temperaKi=temp_ijk+273.15_rsh
     uci=sqrt(0.0625*(uz(km,im,jm)+uz(km,im,jm+1)+uz(km,im-1,jm+1) &
              +uz(km,im-1,jm))**2+0.0625*(vz(km,im,jm)+vz(km,im+1,jm) &
              +vz(km,im,jm-1)+vz(km,im+1,jm-1))**2)

!    variable a volatiliser
     iv=irkm_var_re(ire,1)
     ire4=ire_typ(ire)
   
     IF (cf_henry_adim(iv).eq.0) THEN
   ! Henry constant  reference to temperature = 20degC
       cf_henry_adimi=cf_henry(iv)/cst_gas_univ/(20_rsh+273.15_rsh)
     ELSE
       cf_henry_adimi=cf_henry_adim(iv)
     END IF
   
! phase gazeuse : compose de reference = H2O
!=============================================
! formule de Schwartzenbach page 897
     xkgh2o=(0.2_rsh*ventvent(im,jm)+0.3_rsh)/100._rsh
       
!  si Yg=masse molaire      ==>> l_volkg_cor(iv) true
!  si Yl=diffusivite        ==>> l_volkg_cor(iv) false
     IF (l_volkg_cor(iv) .or. cf_difair(iv)==0.) THEN              
         xkg=xkgh2o*(massmol(iv)/massmol_h2o)**cf_expvol_gas(iv)
     ELSE
         xkg=xkgh2o*(cf_difair(iv)/cf_difair_h2o)**cf_expvol_gas(iv)
     END IF
            
! phase liquide : compose de reference = CO2
!=============================================
!  terme de diffusion turbulente (kld) et engendree par le vent (klv)
!  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     IF (l_volkl_dif(iv)) THEN
!! formule de Churchill 
         xkldo2=5.8e-5*uci**0.969_rsh/d**0.673_rsh
     ELSE 
!! formule de O Connor 
         xkldo2=4.24e-5*uci**0.5_rsh/d**0.5_rsh
     END IF
! formule de Wen et al, 1984 si effet du sel
     IF (l_volkl_sal(iv)) THEN
            xkldo2=xkldo2*exp(0.0127_rsh*sal_ijk)
     END IF
! effet de la temperature - formule de Elmore et West, in Wilcock 1982
     xkldo2=xkldo2*1.024_rsh**(temp_ijk-20.0_rsh)
     xkldco2=xkldo2*(cf_difwat_co2/cf_difwat_o2)**0.5_rsh

 
     IF (ich_volkl_wind(iv) == 1) THEN
! formule de  Schwartzenbach page 905
          xklvco2=4.e-6+4.e-7_rsh*ventvent(im,jm)**2 
     ELSE IF (ich_volkl_wind(iv) == 2) THEN
! formule de Wanninkhof 
          xklvco2=1.25e-6*ventvent(im,jm)**1.64_rsh
     ELSE IF (ich_volkl_wind(iv) == 3) THEN
! formule de Nightingale 
          xklvco2=7.0619e-07*ventvent(im,jm)+2.777e-06_rsh*ventvent(im,jm)**2 
     END IF
         
!  si Yl=masse molaire      ==>> l_volkl_cor : true
!  si Yl=nombre de schmidt  ==>> l_volkl_cor : false      
     IF (l_volkl_cor(iv) .or.  &
         (nschmidt(iv) == 0 .and. cf_difwat(iv) ==0.)) THEN              
         xkl=(xkldco2+xklvco2)*(massmol(iv)/massmol_co2)**cf_expvol_liq(iv)
     ELSE
!  calcul du nombre de Schmidt pour CO2 en fonction de la temperature
         sc_co2=exp(-0.052_rsh*temperaKi+21.71_rsh)
!  calcul du nombre de Schmidt pour le compose en fonction du coef de diffusion moleculaire
         sc_conta=nschmidt(iv)
!         IF (nschmidt(iv) == 0.0_rsh) THEN
!           sc_conta= viscosity/cf_difwat(iv)/density 
!         END IF
         xkl=(xkldco2+xklvco2)*(sc_conta/sc_co2)**cf_expvol_liq(iv)
     END IF
       
! calcul de la constante de Henry en fonction de la temperature 
!   formule de Staudinger et Roberts (2001)
     IF (cf_temp_henry(iv).ne.0) THEN      
          cf_henry_adimi=cf_henry_adimi*10**(-cf_temp_henry(iv)*(1./temperaKi-1./293._rsh))
     END IF
       
! correction de la constante de Henry en fonction de la salinite 
     cf_henry_adimi=cf_henry_adimi*10**(cf_salt_henry(iv)*sal_ijk/68.35_rsh)
       
! concentration dans l air supposee ne pas evoluer 
! et toujours egale a la concentration initiale
     xkol=xkl*xkg*cf_henry_adimi/(xkg*cf_henry_adimi+xkl)
     dc_ire(iv)=-xkol/ztransf*(c(iv)-cini_air(iv)/cf_henry_adimi)

   END IF
    
  END SUBROUTINE incellwat_VOLAT
!
   !!======================================================================

  SUBROUTINE incellwat_PART(ire,km,im,jm,sal_ijk,temp_ijk,c)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_PART  ***
   !&E
   !&E ** Purpose : Calcul des nouvelles concentrations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type PART_EQ ayant lieu a l interieur d une maille dans l eau
   !&E **
   !&E **           PAS DE PARTAGE NON LINEAIRE ICI
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : t
   USE comsubstance, ONLY : nvpc,imud1
   USE contareaddat, ONLY : ire_typ,irkm_var_re,nb_var_re,irkm_sub_re5,l_arh, &
                            minpt_mes_re5,tempref_arh,cf_xetsig_arh,  &
                            ityp_arh,cf_teta_arh,l_pce_re5,cf_pce_re5,l_salt_re5,&
                            cf_salt_re5,cf_distb_re5,minpt_conta,l_kds_re5,  &
                            ndg_kds_re5,cf_kd0_re5,cf_xkds_re5
  
   !! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
!   REAL(KIND=rsh),             INTENT(in) :: xeta     !ssh
   REAL(KIND=rsh),INTENT(in)            :: sal_ijk,temp_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(inout) :: c

   !! * Local declarations
   REAL(KIND=rsh)              :: contot,cmestot,ctot
   REAL(KIND=rsh)              :: fT,xdtemp,fMES,fSAL,fcorKd,alsal
   REAL(KIND=rsh)              :: som,xsalt
   INTEGER                     :: ire5,iv,ivr,isub,ivdl,ndgpoly,idg,IERR_MPI
   REAL(KIND=rsh),DIMENSION(nv_tot)   ::  cv_anc,csorb
      
   ire5=ire_typ(ire)
 ! calcul de la concentration totale de la variable a partager et des MES associees
   iv=irkm_var_re(ire,1)
   contot= c(iv)
   ctot=0.
   cv_anc(:)=c(:)
   cmestot=0.0 
   DO iv=imud1,nvpc
     cmestot= cmestot+c(iv)       
   END DO
   DO ivr=2,nb_var_re(ire)
      iv=irkm_var_re(ire,ivr)
      contot= contot+c(iv)
   END DO
    
 ! calcul des FACTEURS DE variations du coefficients de partage
 ! with Temperature
     fT=1.
     IF (l_arh(ire)) then
       xdtemp=(temp_ijk-tempref_arh(ire))/cf_xetsig_arh(ire)
       fT=(1-ityp_arh(ire))*cf_teta_arh(ire)**xdtemp &
          +ityp_arh(ire)*cf_teta_arh(ire)*exp(-xdtemp**2)
     END IF
  ! with MES (PCE effect)
     fMES=1.
     IF (l_pce_re5(ire5)) then
       fMES=cf_pce_re5(ire5,1)*cmestot**cf_pce_re5(ire5,2)
     END IF
  ! with salinity (salting out effect)
     fSAL=1.
     IF (l_salt_re5(ire5)) then
       IF (sal_ijk.gt.cf_salt_re5(ire5,1)) THEN
        xsalt=cf_salt_re5(ire5,2)*cf_salt_re5(ire5,1)
        fSAL=exp(xsalt)
       ELSE IF (sal_ijk.gt.0.0_rsh) THEN
        xsalt=cf_salt_re5(ire5,2)*sal_ijk
        fSAL=exp(xsalt)
       END IF
     END IF
     fcorKd=fT*fSAL*fMES
     
     IF (l_kds_re5(ire5)) THEN
! polynomial law function of ln(Sal+1) from MOCO or TURNER
        alsal=log(max(0.,sal_ijk)+1)
        cf_distb_re5(ire5,:)=cf_kd0_re5(ire5)
        ndgpoly=ndg_kds_re5(ire5)
        DO idg=1,ndgpoly
           cf_distb_re5(ire5,:)=cf_distb_re5(ire5,:)+   &
           cf_xkds_re5(ire5,idg)*alsal**(ndgpoly+1-idg)
        END DO
        cf_distb_re5(ire5,:)=exp(cf_distb_re5(ire5,:))
     END IF
    
   som=1
   DO ivr=2,nb_var_re(ire)
     iv=irkm_var_re(ire,ivr)
     isub=irkm_sub_re5(ire5,ivr)
     IF(isub == -1) THEN
       csorb(ivr)=max(0._rsh,fcorKd*cf_distb_re5(ire5,ivr)*cmestot) 
     ELSE
       csorb(ivr)=max(0._rsh,fcorKd*cf_distb_re5(ire5,ivr)*c(isub))
     END IF
     som=som+csorb(ivr)
   END DO

!  concentration de l espece dissoute libre
   ivdl=irkm_var_re(ire,1)
   IF (som > epsilon) THEN
     !if(im==188 .and. jm==2 .and. km==1)then
     !  write(*,*)'som,csorb',som,csorb(2:4)
     !  write(*,*)'avant c',c
     !endif
     c(ivdl)=contot/som
     !if(im==188 .and. jm==2 .and. km==1)then
     !  write(*,*)'c(ivdl)',c(ivdl)
     !endif
   ELSE
     c(ivdl)=contot
     csorb(:)=0._rsh
   ENDIF
   ctot=c(ivdl)
!  concentration des autres especes 
   DO ivr=2,nb_var_re(ire)
      iv=irkm_var_re(ire,ivr)
!      c(iv)=c(ivdl)*f_sorpt_lin(ispc)
      c(iv)=c(ivdl)*csorb(ivr)
      ctot=ctot+c(iv)
   END DO
     !if(im==188 .and. jm==2 .and. km==1)then
     !  write(*,*)'ctot,contot',ctot,contot
     !  write(*,*)'apres c',c
     !endif


!  conservativity of the partition must be exact to within 0.1% of ctot
!  otherwise, there is a problem and we stop the simulation
   IF(abs(ctot-contot) > contot*0.001_RSH .and. ctot > epsilon*0.00001_rsh) THEN
      write(*,*)'ERROR : partition non conservative at t=',t
      write(*,*)'Ctot before =',contot,' Ctot after =',ctot
      write(*,*)'at km,im,jm = ',km,im,jm,'  som=',som
      write(*,*)'salinity = ',sal_ijk,'temperature = ',temp_ijk
      write(*,*)'csorb(2:nb_var_re(ire))',csorb(2:nb_var_re(ire))
      write(*,*)'cv_anc(1:nv_tot)',cv_anc
      write(*,*)'c(1:nv_tot)',c
      CALL_MPI MPI_FINALIZE(IERR_MPI)
      STOP
   END IF
  
  END SUBROUTINE incellwat_PART
   !!======================================================================
  SUBROUTINE incellwat_PART_conta(ire,km,im,jm,sal_ijk,temp_ijk,c)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_PART_conta ***
   !&E
   !&E ** Purpose : Calcul des nouvelles concentrations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type PART_EQ ayant lieu a l interieur d une maille dans l eau
   !&E **           POUR CONTAMINANT AVEC PARTAGE DE DIFFERENTES ESPECES
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comsubstance, ONLY : nvpc,irk_mod,imud1
   USE contareaddat, ONLY : ire_typ,icon_reac,nb_spcmax,ivar_spc,ispc_var,  &
                            nb_spc_diss,nb_spc_part_eq,isit_spc,nb_spc_tot,      &
                            nb_spc_part_eq,nb_spc_tot,l_arh,icon_var, &
                            tempref_arh,cf_xetsig_arh,ityp_arh,cf_teta_arh, &
                            l_pce_re5,cf_pce_re5,l_salt_re5,cf_salt_re5,l_metal,&
                            l_kds_re5,cf_distb_re5,cf_kd0_re5,ndg_kds_re5,      &
                            cf_xkds_re5,minpt_mes_re5,minpt_conta,irkm_var_re,   &
                            l_sorpt_nlin,l_langmuir,cf_kfkl_nlin,cf_expqmax_nlin,irkm_sub_re5
                           
#ifdef key_species
   USE contareaddat, ONLY : frac_cwat,ifrac_var_spc,nb_frac_var
#endif

   !! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
   REAL(KIND=rsh),INTENT(in)            :: sal_ijk,temp_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(inout) :: c

   !! * Local declarations
   REAL(KIND=rsh)              :: contot,cmestot,csubtot,ctot
   REAL(KIND=rsh)              :: fT,xdtemp,fMES,fSAL,fcorKd,alsal
   REAL(KIND=rsh)              :: somd,somp,som,xsalt,ff
   INTEGER                     :: ire5,iv,ivr,isub,ndgpoly,idg,ispc,ico,ifracv  
   REAL(KIND=rsh),DIMENSION(nb_spcmax)   ::  cspc_anc,cspc,csorb
   
#ifdef key_species
   ire5=ire_typ(ire)
   ico=icon_reac(ire)
   contot=0.
   csorb(:)=0.
   cspc(:)=0.
   cspc_anc(:)=0.
   cmestot=0.0 
   DO iv=imud1,nvpc
     cmestot= cmestot+c(iv)       
   END DO
   csubtot=0.0 
 ! calcul de la concentration totale de la variable a partager et des MES associees
 ! c est a dire du contaminant echangeable = somme des especes dissoutes et particulaires a l equilibre
   ispc=1
   iv=irk_mod(ivar_spc(ico,ispc))
   cspc_anc(ispc)=frac_cwat(iv,ifrac_var_spc(ico,ispc),km,im,jm)*c(iv)
   contot= contot+cspc_anc(ispc)
   DO ispc=2,nb_spc_diss(ico)+nb_spc_part_eq(ico)
     iv=irk_mod(ivar_spc(ico,ispc))
     cspc_anc(ispc)=frac_cwat(iv,ifrac_var_spc(ico,ispc),km,im,jm)*c(iv)
     contot= contot+cspc_anc(ispc)
     isub=irk_mod(isit_spc(ico,ispc))
     csubtot= csubtot+c(isub)
   END DO
   DO ispc=nb_spc_diss(ico)+nb_spc_part_eq(ico)+1,nb_spc_tot(ico)
     iv=irk_mod(ivar_spc(ico,ispc))
     cspc_anc(ispc)=frac_cwat(iv,ifrac_var_spc(ico,ispc),km,im,jm)*c(iv)
     cspc(ispc)=cspc_anc(ispc)
   END DO
    
 ! calcul des FACTEURS DE variations du coefficients de partage
 ! with Temperature
     fT=1.
     IF (l_arh(ire)) then
       xdtemp=(temp_ijk-tempref_arh(ire))/cf_xetsig_arh(ire)
       fT=(1-ityp_arh(ire))*cf_teta_arh(ire)**xdtemp &
          +ityp_arh(ire)*cf_teta_arh(ire)*exp(-xdtemp**2)
     END IF
  ! with MES (PCE effect)
     fMES=1.
     IF (l_pce_re5(ire5)) then
       fMES=cf_pce_re5(ire5,1)*cmestot**cf_pce_re5(ire5,2)
     END IF
  ! with salinity (salting out effect)
     fSAL=1.
     IF (l_salt_re5(ire5)) then
       IF (sal_ijk.gt.cf_salt_re5(ire5,1)) THEN
        xsalt=cf_salt_re5(ire5,2)*cf_salt_re5(ire5,1)
        fSAL=exp(xsalt)
       ELSE IF (sal_ijk.gt.0.) THEN
        xsalt=cf_salt_re5(ire5,2)*sal_ijk
        fSAL=exp(xsalt)
       END IF
     END IF
     fcorKd=fT*fSAL*fMES
     
     IF (l_metal(ico) .and.l_kds_re5(ire5)) THEN
! polynomial law function of ln(Sal+1) from MOCO or TURNER
        alsal=log(max(0.,sal_ijk)+1)
        cf_distb_re5(ire5,:)=cf_kd0_re5(ire5)
        ndgpoly=ndg_kds_re5(ire5)
        DO idg=1,ndgpoly
           cf_distb_re5(ire5,:)=cf_distb_re5(ire5,:)+   &
           cf_xkds_re5(ire5,idg)*alsal**(ndgpoly+1-idg)
        END DO
        cf_distb_re5(ire5,:)=exp(cf_distb_re5(ire5,:))
     END IF
      
   somd=1
   somp=0
!    somnl=0
   som=1
! ispc est ici = ivr numero de variable active
   IF (cmestot.gt.minpt_mes_re5.and.contot.gt.minpt_conta) THEN
     DO ispc=2,nb_spc_diss(ico)+nb_spc_part_eq(ico)
       iv=irkm_var_re(ire,ispc)
       isub=irk_mod(isit_spc(ico,ispc))
!     csorb(ispc)=f_sorpt_lin(ispc)*fcorKd*cf_distb_re5(ire5,ivr)*c(isub)
       csorb(ispc)=fcorKd*cf_distb_re5(ire5,ispc)*c(isub)
       som=som+csorb(ispc)
       IF (itypv(iv).eq.6) THEN
         somd=somd+csorb(ispc)
       ELSE 
         somp= somp+csorb(ispc)
!  pas de domaine non lineaire
!       IF (l_sorpt_nlin(icon_var(iv))) THEN
!c MES associe
!           ipartesp=ipart_sub(isub_esp(ispc))
!c il existe un ou des domaine non lineaires                 
!c modele DMM : un seul domaine non lineaire possible 
!c  et une seule espece avec domaine non lineaire
!           IF (l_langmuir) THEN
!            somnl=cf_kfkl_nlin(ispc,1)*f_sorpt_nlin(ispc,1)*  &
!                cf_expqmax_nlin(ispc,1)*cpart_v(ipartesp,i)
!           END IF
!       END IF
        END IF
      END DO

!  concentration de l espece dissoute libre
     cspc(1)=contot/som
     ctot=cspc(1)
!  concentration des autres especes 
     DO ivr=2,nb_var_re(ire)
        isub=irk_mod(isit_spc(ico,ivr))
!           cspc(ivr)=cspc(1)*f_sorpt_lin(ispc)
         cspc(ivr)=cspc(1)*csorb(ivr)
         ctot=ctot+cspc(ivr)
     END DO
            
   ELSE IF (contot<0.0_rsh) THEN
! no partition if contot<0.

   ELSE IF (contot.le.minpt_conta) THEN
!  concentration dans l espece dissoute libre
      cspc(1)=contot
      cspc(2:nb_spc_diss(ico)+nb_spc_part_eq(ico))=0.
         
   ELSE 
   ! si MES faible < seuil : tout est dans le dissous
      somd=1
      somp=0
      som=1
      DO ivr=2,nb_spc_diss(ico)
       isub=irk_mod(isit_spc(ico,ivr))
       csorb(ivr)=fcorKd*cf_distb_re5(ire5,ivr)*c(isub)
       som=som+csorb(ivr)
       somd=somd+csorb(ivr)
      END DO

!  concentration de l espece dissoute libre
      cspc(1)=contot/som
      ctot=cspc(1)
!  concentration des autres especes dissoutes
      DO ivr=2,nb_spc_diss(ico)
        isub=irkm_sub_re5(ire5,ivr)
        cspc(ivr)=cspc(1)*cf_distb_re5(ire5,ivr)*fcorKd*c(isub)
        ctot=ctot+cspc(ivr)
      END DO
      cspc(nb_spc_diss(ico)+1:nb_spc_diss(ico)+nb_spc_part_eq(ico))=0.
   END IF
   
! repartition dans les differentes variables et fractions
   DO iv=1,nv_state
      IF (icon_var(iv).eq.ico) THEN
        c(iv)=0.
        DO ifracv=1,nb_frac_var(iv)
          c(iv)=c(iv)+cspc(ispc_var(iv,ifracv))
        END DO
        ff=0.
        DO ifracv=1,nb_frac_var(iv)-1
          frac_cwat(iv,ifrac_var_spc(ico,ifracv),km,im,jm)=cspc(ispc_var(iv,ifracv))/(c(iv)+1.e-10)
          ff=ff+frac_cwat(iv,ifrac_var_spc(ico,ifracv),km,im,jm)
        END DO
        frac_cwat(iv,nb_frac_var(iv),km,im,jm)=1.-ff
      END IF
    END DO
  
#endif
  END SUBROUTINE incellwat_PART_conta
   !!======================================================================

  SUBROUTINE incellwat_PART_medit(ire,km,im,jm,sal_ijk,temp_ijk,c)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE incellwat_PART_medit  ***
   !&E
   !&E ** Purpose : Calcul des nouvelles concentrations dues aux phenomenes PHY-BIO-CHI
   !&E **           du type PART_EQ ayant lieu a l interieur d une maille dans l eau
   !&E **
   !&E **  formulation differente pour faibles concentrations en substrats
   !&E **
   !&E **          PAS DE PARTAGE NON LINEAIRE ICI
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          2008-06  B. thouvenin
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : t
   USE comsubstance, ONLY : nvpc,imud1
   USE contareaddat, ONLY : ire_typ,irkm_var_re,nb_var_re,irkm_sub_re5,l_arh, &
                            minpt_mes_re5,tempref_arh,cf_xetsig_arh,  &
                            ityp_arh,cf_teta_arh,l_pce_re5,cf_pce_re5,l_salt_re5,&
                            cf_salt_re5,cf_distb_re5,minpt_conta,l_kds_re5,  &
                            ndg_kds_re5,cf_kd0_re5,cf_xkds_re5
  
   !! * Arguments
   INTEGER, INTENT( in )             ::   km,im,jm      ! ocean spatial-step index
   INTEGER, INTENT( in )             ::   ire      
!   REAL(KIND=rsh),             INTENT(in) :: xeta     !ssh
   REAL(KIND=rsh),INTENT(in)            :: sal_ijk,temp_ijk
   REAL(KIND=rsh),DIMENSION(nv_tot),INTENT(inout) :: c

   !! * Local declarations
   REAL(KIND=rsh)              :: contot,cmestot,ctot,contot_ads
   REAL(KIND=rsh)              :: fT,xdtemp,fMES,fSAL,fcorKd,alsal
   REAL(KIND=rsh)              :: som,xsalt
   INTEGER                     :: ire5,iv,ivr,isub,ivdl,ndgpoly,idg,IERR_MPI
   REAL(KIND=rsh),DIMENSION(nv_tot)   ::  cv_anc,csorb,ctemp,cdl
      
   ire5=ire_typ(ire)
 ! calcul de la concentration totale de la variable a partager et des MES associees
   iv=irkm_var_re(ire,1)
   contot= c(iv)
  ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART cdiss',iv,c(iv)
   ctot=0.
   cv_anc(:)=c(:)
   cmestot=0.0 
   DO iv=imud1,nvpc
     cmestot= cmestot+c(iv)       
   END DO
   DO ivr=2,nb_var_re(ire)
      iv=irkm_var_re(ire,ivr)
      contot= contot+c(iv)
    ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART c',iv,c(iv)
   END DO
    ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART contot',contot
 ! calcul des FACTEURS DE variations du coefficients de partage
 ! with Temperature
     fT=1.
     IF (l_arh(ire)) then
       xdtemp=(temp_ijk-tempref_arh(ire))/cf_xetsig_arh(ire)
       fT=(1-ityp_arh(ire))*cf_teta_arh(ire)**xdtemp &
          +ityp_arh(ire)*cf_teta_arh(ire)*exp(-xdtemp**2)
     END IF
  ! with MES (PCE effect)
     fMES=1.
     IF (l_pce_re5(ire5)) then
       fMES=cf_pce_re5(ire5,1)*cmestot**cf_pce_re5(ire5,2)
     END IF
  ! with salinity (salting out effect)
     fSAL=1.
     IF (l_salt_re5(ire5)) then
       IF (sal_ijk.gt.cf_salt_re5(ire5,1)) THEN
        xsalt=cf_salt_re5(ire5,2)*cf_salt_re5(ire5,1)
        fSAL=exp(xsalt)
       ELSE IF (sal_ijk.gt.0.0_rsh) THEN
        xsalt=cf_salt_re5(ire5,2)*sal_ijk
        fSAL=exp(xsalt)
       END IF
     END IF
     fcorKd=fT*fSAL*fMES
     
     IF (l_kds_re5(ire5)) THEN
! polynomial law function of ln(Sal+1) from MOCO or TURNER
        alsal=log(max(0.,sal_ijk)+1)
        cf_distb_re5(ire5,:)=cf_kd0_re5(ire5)
        ndgpoly=ndg_kds_re5(ire5)
        DO idg=1,ndgpoly
           cf_distb_re5(ire5,:)=cf_distb_re5(ire5,:)+   &
           cf_xkds_re5(ire5,idg)*alsal**(ndgpoly+1-idg)
        END DO
        cf_distb_re5(ire5,:)=exp(cf_distb_re5(ire5,:))
     END IF

! calcul pour chaque particule l equilibre avec le dissous
! equilibre entre les especes qui sont sont en adsorption
   ivdl=irkm_var_re(ire,1)
   contot_ads=c(ivdl)
   som=1
   DO ivr=2,nb_var_re(ire)
     iv=irkm_var_re(ire,ivr)
     isub=irkm_sub_re5(ire5,ivr)
     IF(isub == -1) THEN
       csorb(ivr)=max(0._rsh,fcorKd*cf_distb_re5(ire5,ivr)*cmestot) 
     ELSE
       csorb(ivr)=max(0._rsh,fcorKd*cf_distb_re5(ire5,ivr)*c(isub))
     END IF
     ctemp(iv)=csorb(ivr)*(c(iv)+c(ivdl))/(1+csorb(ivr))
     if(ctemp(iv)<c(iv))then
        ! desorption
        cdl(iv)=c(iv)-ctemp(iv)
        c(iv)=0.0_rsh
        c(isub)=0.0_rsh
     else
         ! adsorption
        cdl(iv)=0.0_rsh
        som=som+csorb(ivr)
          ! attention ne marche pas pour isub=-1
     endif
     contot_ads=contot_ads+cdl(iv)+c(iv)
    ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART csorb,cdl,ctemp,c',iv,ivr,isub,c(isub),csorb(ivr),cdl(iv),ctemp(iv),c(iv)
   END DO
    ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART contot_ads som',contot_ads,som

!  concentration de l espece dissoute libre
   IF (som > epsilon) THEN
     !if(im==188 .and. jm==2 .and. km==1)then
     !  write(*,*)'som,csorb',som,csorb(2:4)
     !  write(*,*)'avant c',c
     !endif
     c(ivdl)=contot_ads/som
     !if(im==188 .and. jm==2 .and. km==1)then
     !  write(*,*)'c(ivdl)',c(ivdl)
     !endif
   ELSE
     c(ivdl)=contot_ads
     csorb(:)=0._rsh
   ENDIF
   ctot=c(ivdl)
  ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART cnouvdiss',c(ivdl)

!  concentration des autres especes (en remettant les valeurs des especes en desorption)
   DO ivr=2,nb_var_re(ire)
      iv=irkm_var_re(ire,ivr)
      if(cdl(iv)==0.0_rsh) then
        ! adsorption
        c(iv)=c(ivdl)*csorb(ivr)
      else
        ! desorption
        c(iv)=ctemp(iv)
      endif
      ctot=ctot+c(iv)
   !  if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART cnouv',c(iv)
   END DO
     !if(im==188 .and. jm==2 .and. km==1)then
     !  write(*,*)'ctot,contot',ctot,contot
     !  write(*,*)'apres c',c
     !endif
    ! if(im==320.and.jm==163.and.km==kmax)write(*,*)'PART ctotnouv',ctot


!  conservativity of the partition must be exact to within 0.1% of ctot
!  otherwise, there is a problem and we stop the simulation
   IF(abs(ctot-contot) > contot*0.001_RSH .and. ctot > epsilon*0.00001_rsh) THEN
      write(*,*)'ERROR : partition non conservative at t=',t
      write(*,*)'Ctot before =',contot,' Ctot after =',ctot
      write(*,*)'at km,im,jm = ',km,im,jm,'  som=',som
      write(*,*)'salinity = ',sal_ijk,'temperature = ',temp_ijk
      write(*,*)'csorb(2:nb_var_re(ire))',csorb(2:nb_var_re(ire))
      write(*,*)'cv_anc(1:nv_tot)',cv_anc
      write(*,*)'c(1:nv_tot)',c
      CALL_MPI MPI_FINALIZE(IERR_MPI)
      STOP
   END IF
  
  END SUBROUTINE incellwat_PART_medit
   !!======================================================================

!
#endif
 END MODULE
