Old code:      

      denom1 = drp(i)*drm(i)*dr2
      denom2 = rho(i)*cp(i)*r2(i)*dr2
C
C   Set things up in tridiagonal form
      a(i) = u(i)*rip2/denom1 + 2.*rim2*appakm/drm(i)/denom2
      b(i) = - u(i)*(rip2-rim2)/denom1 - (rim2*appakm/drm(i)
     2   + rip2*appakp/drp(i))*2./denom2
      c(i) = - u(i)*rim2/denom1 + 2.*rip2*appakp/drp(i)/denom2
      d = q(i)/(rho(i)*cp(i)) - u(i)*g(i)/cp(i)
      rhs(i) = a(i)*T(i-1) + b(i)*T(i) + c(i)*T(i+1) + d

Corrected code:

      denom1 = drp(i)*drm(i)*dr2
      denom2 = rho(i)*cp(i)*r2(i)*dr2
C
C   Set things up in tridiagonal form
      a(i) = u(i)*drp2/denom1 + 2.*rim2*appakm/drm(i)/denom2
      b(i) = - u(i)*(drp2-drm2)/denom1 - (rim2*appakm/drm(i)
     2   + rip2*appakp/drp(i))*2./denom2
      c(i) = - u(i)*drm2/denom1 + 2.*rip2*appakp/drp(i)/denom2
      d = q(i)/(rho(i)*cp(i)) - u(i)*g(i)/cp(i)
      rhs(i) = a(i)*T(i-1) + b(i)*T(i) + c(i)*T(i+1) + d

More corrections: Jeans escape

Old code:

      ujeans(i) = rootpi/2. * exp(-alamp(i))*sqrt(2.*u02(i))*
     2  (1.+alamp(i))

Corrected code:

      ujeans(i) = 1./(2.*rootpi) * exp(-alamp(i))*sqrt(2.*u02(i))*
     2  (1.+alamp(i))