anchor.diff:  !
   !          DO itrc=1,NBGCPW
   !            DO k=1,Nbed-1
  -!              cff1=porod(k)*usmd(i,k)
  -!              cff2=dens*(1.0d0-porod(k))*usmd(i,k)*Kads(i,k,itrc)
  +!              cff1=porod(i,k)*usmd(i,k)
  +!              cff2=dens*(1.0d0-porod(i,k))*usmd(i,k)*Kads(i,k,itrc)
   !              pwflux(i,k,itrc)=(cff1+cff2)*(pw(i,k,itrc)-pw(i,k+1,itrc))/dzd(i,k)
   !            END DO
   !            pwflux(i,Nbed,itrc)=pwflux(i,Nbed-1,itrc)
  
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
              
             !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)
              Rpw(iwH2S)=Rpw(iwH2S)-R(13)
              Rsm(iFeS_)=Rsm(iFeS_)+R(13)
  
negative.code:             cff1=bw(i,iwO2_)*32.0d0/1000.0d0
              cff2=11.5d0*0.748d0**cff1
              cff3=cff2*1.02d0**(tsm(i,k)-20.0d0)
              !Kads(i,k,iwPO4)=poro(i,k)/(1.0d0-poro(i,k))/cff3
            END DO
  !
            DO itrc=1,NBGCPW
  