anchor.diff:            DO itrc=1,NBGCSM
               Thomas_a(1)=0.0d0
               DO k=2,Nbed
  -              fac=dens*(1.0d0-poro(k))
  +              fac=dens*(1.0d0-poro(i,k))
                 cff=dtBgc/fac/dz(i,k)
  -              cff3=dens*(1.0d0-porod(k-1))*DBsd(i,k-1)
  +              cff3=dens*(1.0d0-porod(i,k-1))*DBsd(i,k-1)
                 Thomas_a(k)=cff*cff3/dzd(i,k-1)
               END DO
   !
               Thomas_c(Nbed)=0.0d0
               DO k=1,Nbed-1
  -              fac=dens*(1.0d0-poro(k))
  +              fac=dens*(1.0d0-poro(i,k))
                 cff=dtBgc/fac/dz(i,k)
  -              cff4=dens*(1.0d0-porod(k))*DBsd(i,k)
  +              cff4=dens*(1.0d0-porod(i,k))*DBsd(i,k)
                 Thomas_c(k)=cff*cff4/dzd(i,k)
               END DO
   !
               DO k=1,Nbed
  -              fac=dens*(1.0d0-poro(k))
  +              fac=dens*(1.0d0-poro(i,k))
                 Thomas_b(k)=1.0d0+Thomas_a(k)+Thomas_c(k)
                 Thomas_d(k)=sm(i,k,itrc)+                                 &
        &                    (smflux(i,k-1,itrc)-smflux(i,k,itrc))*        &
  
positive.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(k)/(1.0d0-poro(k))/cff3
            END DO
  !
            DO itrc=1,NBGCPW
  
negative.code: !
  !          DO itrc=1,NBGCPW
  !            DO k=1,Nbed-1
  !              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)
  