function dTdp=moistlapse2(p,T,p0,T0,satmax,icefracFinal,TFinal,loadopt,fracload,Tunload)

%Input arguments, MKS units unless noted:
%p: current pressure
%T: current temperature
%p0: cloud base pressure
%T0: cloud base temperature
%satmax: equilibrium supersaturation
%icefracFinal: fraction of condensate that eventually freezes. This is
%              always set to 1 in the GRL paper
%TFinal: temperature at which freezing ends. Either -4 or -40C. Enter in deg C.
%loadopt: loading treatment. See getvapor.m
%fracload: final fraction of condensate that is loaded. Used in conjunction
%          with loadopt. See getvapor.m
%Tunload: Temperature at which unloading begins

g=9.8;    %acceleration due to gravity
cpd=1005; %heat capacity of dry air
cpv=1870; %heat capacity of water vapor
cpl=4200; %heat capacity of liquid water
cpi=2100; %heat capacity of ice
Lv=2.5e6; %latent heat of vaporization
Lf=3.34e5;%latent heat of freezing
Ls=Lv+Lf; %latent heat of sublimation
R1=461.5; %water vapor gas constant
R0=287;   %dry air gas constant
ep=R0/R1; 

%saturation vapor pressure equations, liquid and ice
es=@(T) 100*10.^(-7.90298*(373.16./T-1)+5.02808*log10(373.16./T)...
    -1.3816e-7*(10.^(11.344*(1-T/373.16))-1)+...
    8.1328e-3*(10.^(-3.49149*(373.16./T-1))-1)+log10(1013.246));
esi=@(T) 100*10.^(-9.09718*(273.16./T-1)-3.56654*log10(273.16./T)+...
    0.876793*(1-T/273.16)+log10(6.1071));

%get current saturation ratio and saturation ratio gradient wrt pressure
[satfac,dsatfac]=dfrac(p,p0,satmax);

%get current liquid and ice fractions, gradients of these fractions wrt
%temperature, vapor mixing ratio, condensate mixing ratio, vapor pressure
conditions=getvapor(p,T,p0,T0,satmax,icefracFinal,TFinal,loadopt,fracload,Tunload);
[liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};

%linearly weighted latent heats
lF1=liqfrac*Lv+icefrac*Ls;

%Now use Equation 12 in supplemental information

B=rF1*lF1*(1/pF1+1/(p-pF1));

%last term in denominator in Eq.12
%Use a finite difference approximation to dp21/dT
V=B*satfac*(es(T+.01)-es(T-.01))/.02;

%middle term in denonimator in Eq. 12
W=Lf*rF*dliqfrac;

%last term in numerator in Eq. 12
A=B*es(T)*dsatfac;

%total heat capacity
Cp=cpd+cpv*rF1+cpl*rF*liqfrac+cpi*rF*icefrac; 

%the lapse rate
dTdp=((T*R0+rF1*lF1)/(p-pF1)-A)/(Cp+W+V);

