anchor.diff:    double precision, intent(in) :: Nwimps
     double precision, intent(out) :: EvapRate
     double precision :: Earg(nlines), muarray(niso_in),nabund(niso_in,nlines),sigma_N(niso_in)
  -  double precision :: suparg,trapz
  +  double precision :: suparg,trapz,Knud,rchi
     double precision :: mdmg,mnucg, Tc, rhoc,mn, Tw,nin,escFrac(nlines),vescc(nlines),nxIso(nlines),mfp(nlines),scatrate(nlines)
  -  double precision, parameter :: kBeV=8.617e-5
  +  double precision, parameter :: kBeV=8.617d-5, GN =6.674d-8, kB=1.3806d-16
   
     integer :: i, j, k
     mdmg = mdm*1.78d-24
  
positive.code:   nxIso = exp(mdm*vescc**2/2./Tw);
    nin = 4.d0*pi*trapz(tab_r,tab_r**2.*nxIso,nlines) !%niso norm
    nxIso = nxIso/nin
    !Fraction of the kinetic distribution above the local escape velocity
    escFrac = sqrt(2.d0/pi)*vescc*sqrt(mdm/Tw)*exp(-mdm*vescc**2/Tw/2.d0) - derf(sqrt(mdm/Tw/2.d0)*vescc) + 1.d0;
    ! print*,"EscFrac = ", escFrac,
  
negative.code: 
      end do
              Tw_out = sum(Tw_out_num)/sum(Tw_out_denom);
              dT = abs(Tw - Tw_out)/Tw;
              Tw = Tw_out;
    end do
    Tw = Tw/TcGeV;
  
    ! open(55,file = "svTw.dat")
    ! do j=1,nlines
    ! write(55,*) tab_r(j),sv(j), nxIso(j), nabund(1,j),TGeV(j),tab_vesc(j),TwK
    ! end do
    ! close(55)
  
  
  
  