anchor.diff:              
               F(i,35)=F(i,35)-dz(i,k)*R(13)         !Fe 
               F(i,36)=F(i,36)-dz(i,k)*R(13)         !H2S
  -
  +            
  +            R(14)=K14*cff1*sm(i,k,iFeS_)*sm(i,k,iS0__)
  +            Rsm(iFeS_)=Rsm(iFeS_)-R(14)
  +            Rsm(iS0__)=Rsm(iS0__)-R(14)
  +            Rsm(iFeS2)=Rsm(iFeS2)+R(14)
  +            
  +            F(i,37)=F(i,37)-dz(i,k)*R(14)         !FeS
  +            F(i,38)=F(i,38)-dz(i,k)*R(14)         !S0  
  +            
               R(15)=K15*cff1*sm(i,k,iFeS_)*pw(i,k,iwH2S)*pw(i,k,iwSO4)
               Rsm(iFeS_)=Rsm(iFeS_)-R(15)*4.0d0
               Rpw(iwH2S)=Rpw(iwH2S)-R(15)*3.0d0
  
positive.code:             
              F(i,33)=F(i,33)-dz(i,k)*R(12)         !H2S
              F(i,34)=F(i,34)-dz(i,k)*R(12)         !MnOA
  !
  !  - Re-oxydations using FeOOH (A,B,=PO4)
  !
             !R(10)=K10*cff1*pw(i,k,iwH2S)*(sm(i,k,iFeOA)+sm(i,k,iFeOB)) ! (B)
              R(10)=K10*cff1*pw(i,k,iwH2S)*sm(i,k,iFeOA) ! (A)
              Rpw(iwH2S)=Rpw(iwH2S)-R(10)
              Rsm(iFeOA)=Rsm(iFeOA)-R(10)*2.0d0
              Rpw(iwFe_)=Rpw(iwFe_)+R(10)*2.0d0
              Rsm(iS0__)=Rsm(iS0__)+R(10)
              
              F(i,27)=F(i,27)-dz(i,k)*R(10)         !H2S
              F(i,28)=F(i,28)-dz(i,k)*R(10)*2.0d0   !FeOA
  
              R(25)=K25*cff1*pw(i,k,iwH2S)*sm(i,k,iFeOP)
              Rpw(iwH2S)=Rpw(iwH2S)-R(25)
              Rsm(iFeOP)=Rsm(iFeOP)-R(25)*2.0d0
              Rpw(iwFe_)=Rpw(iwFe_)+R(25)*2.0d0
              Rsm(iS0__)=Rsm(iS0__)+R(25)
              Rpw(iwPO4)=Rpw(iwPO4)+R(25)*2.0d0
              
              F(i,29)=F(i,29)-dz(i,k)*R(25)         !H2S
              F(i,30)=F(i,30)-dz(i,k)*R(25)*2.0d0   !FeOP
  !
  !  - Re-oxydations using SO4
  !
             !R(23)=K23*poro(k)*pw(i,k,iwCH4)*pw(i,k,iwSO4)
  !
  !  - Settlings
  !
             !R(13)=K13*poro(k)*pw(i,k,iwFe_)*HS/H1/KsFeS-1.0d0 ! (W)
              R(13)=K13*poro(k)*pw(i,k,iwFe_)*pw(i,k,iwH2S) ! (F)(A)
              Rpw(iwFe_)=Rpw(iwFe_)-R(13)
  
negative.code:             
              F(i,19)=F(i,19)-dz(i,k)*R(6)          !NH4
              F(i,20)=F(i,20)-dz(i,k)*R(6)*2.0d0    !O2
              
             !R(7)=K07*cff1*(sm(i,k,iFeOA)+sm(i,k,iFeOB))*pw(i,k,iwPO4) ! (B)
              R(7)=K07*cff1*sm(i,k,iFeOA)*pw(i,k,iwPO4) ! (A)
              Rsm(iFeOA)=Rsm(iFeOA)-R(7)
              Rpw(iwPO4)=Rpw(iwPO4)-R(7)
              Rsm(iFeOP)=Rsm(iFeOP)+R(7)
              
              F(i,21)=F(i,21)-dz(i,k)*R(7)         !FeOA
              F(i,22)=F(i,22)-dz(i,k)*R(7)         !PO4
  !
  !  - Re-oxydations using MnO2
  !
             !R(8)=K08*cff1*pw(i,k,iwFe_)*(sm(i,k,iMnOA)+sm(i,k,iMnOB)) ! (B)
              R(8)=K08*cff1*pw(i,k,iwFe_)*sm(i,k,iMnOA) ! (A)
  