d_ebalance.m

function dM=d_ebalance(Tsfc,T_fine,Vf,pres_fine,ea,sig,emissivity,...
    Cp,xLs,rho_air,De_h,B1,B2,tau,cc)
% cc=FOREST.cc;
% tau=FOREST.tau;
ems=emissivity.ems;
%analytical deriative of L
dL=-4*Tsfc^3*ems*sig*(Vf*ems - ems + cc*ems - Vf*cc*ems - ...
    cc*ems*tau + Vf*cc*ems*tau + 1);
%analytical derivative of T=latent+sensible given unstable
%conditions (Tsfc > Tair)
if Tsfc > T_fine
    dT = Cp*De_h*rho_air*((B1*(T_fine - Tsfc))/...
        (B2*(Tsfc - T_fine)^(1/2) + 1) - 1) + ...
        Cp*De_h*rho_air*(T_fine - Tsfc)*...
        (B1/(B2*(Tsfc - T_fine)^(1/2) + 1) - ...
        (B1*B2*(Tsfc - T_fine)^(1/2))/(2*(B2*(Tsfc - T_fine)^(1/2) ...
        + 1)^2)) + (311*De_h*rho_air*xLs*(ea - ...
        (12223*exp(((5613*Tsfc)/250 - 6.1328e3)/...
        (Tsfc - 3/5)))/20)*(B1/(B2*(Tsfc - T_fine)^(1/2) + 1) - ...
        (B1*B2*(Tsfc - T_fine)^(1/2))/(2*(B2*(Tsfc - T_fine)^(1/2) ...
        + 1)^2)))/(500*pres_fine) + ...
        (3801353*De_h*rho_air*xLs*exp(((5613*Tsfc)/250 - ...
        6.1328e3)/(Tsfc - 3/5))*(5613/(250*(Tsfc - 3/5)) - ...
        ((5613*Tsfc)/250 - 6.1328e3)/(Tsfc - 3/5)^2)*...
        ((B1*(T_fine - Tsfc))/(B2*(Tsfc - T_fine)^(1/2) + 1) - 1))...
        /(10000*pres_fine);
    %analytical derivative of T given stable conditions
elseif Tsfc < T_fine
    dT = (B1*Cp*De_h*rho_air*(T_fine - Tsfc))/...
        ((B1*(T_fine - Tsfc))/2 + 1)^3 - (Cp*De_h*rho_air)/...
        ((B1*(T_fine - Tsfc))/2 + 1)^2 - ...
        (3801353*De_h*rho_air*xLs*exp(((5613*Tsfc)/250 - ...
        6.1328e3)/(Tsfc - 3/5))*(5613/(250*(Tsfc - 3/5)) - ...
        ((5613*Tsfc)/250 - 6.1328e3)/(Tsfc - 3/5)^2))/...
        (10000*pres_fine*((B1*(T_fine - Tsfc))/2 + 1)^2) + ...
        (311*B1*De_h*rho_air*xLs*(ea - (12223*exp(((5613*Tsfc)/...
        250 - 6.1328e3)/(Tsfc - 3/5)))/20))/(500*pres_fine*...
        ((B1*(T_fine - Tsfc))/2 + 1)^3);
%     analytical derivative of T given neutral conditions
else
    dT=- Cp*De_h*rho_air - ...
        (3801353*De_h*rho_air*xLs*exp(((5613*Tsfc)/250 - 6.1328e3)/...
        (Tsfc - 3/5))*(5613/(250*(Tsfc - 3/5)) - ...
        ((5613*Tsfc)/250 - 6.1328e3)/(Tsfc - 3/5)^2))/(10000*pres_fine);
end
%total melt derivative
dM=dL+dT;%;+dG;
end