anchor.diff:              Rdomf=KDOMf*cff*pw(i,k,iwDOMf)                       !Nchange
               Rdoms=KDOMs*cff*pw(i,k,iwDOMs)
                
  -            F(i,5)=F(i,5)-dz(i,k)*Rdomf
  -            F(i,6)=F(i,6)-dz(i,k)*Rdoms
  +            F(i,5)=F(i,5)+dz(i,k)*Rdomf
  +            F(i,6)=F(i,6)+dz(i,k)*Rdoms
   !
   !  - Wijsman
   !
  
positive.code:             Rpw(iwFe_)=Rpw(iwFe_)+R(18)
              Rpw(iwSO4)=Rpw(iwSO4)+R(18)*2.0d0
              
              F(i,46)=F(i,46)-dz(i,k)*R(18)         !FeS2
              F(i,47)=F(i,47)-dz(i,k)*R(18)*3.5d0   !O2
              
              R(19)=K19*cff1*sm(i,k,iS0__)
              IF (pw(i,k,iwH2S).lt.H2Sstop) THEN
  
negative.code:             Rpw(iwH2S)=Rpw(iwH2S)-R(13)
              Rsm(iFeS_)=Rsm(iFeS_)+R(13)
              
              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
  