function [f,df]=dfrac(p,p0,satmax)

fac=12000;
%saturation ratio
f=(1-satmax)*exp((p-p0)/fac)+satmax;
%saturation ratio gradient wrt pressure
df=(f-satmax)/fac;


