I k1local,k2local,
k12 = k1local*k2local
b = x2 + k1local*x + k12
db = 2.0*x + k1local
fn = k1local*x*diclocal/b +
df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -
co2star=diclocal*htotal2/(htotal2 + k1local*htotal
& + k1local*k2local)
_RL k1local, k2local
stuff = (1.0-gamm)*(1.0-gamm)*k1local*k1local
& - 4.0*k1local*k2local*(1.0-2.0*gamm)
hnew = 0.5*( (gamm-1.0)*k1local + sqrt(stuff) )
& (1.0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))