function y=f(Tc_new,Tc_old,wc,P,qvc_old,Te,qve,dt,lambda,r)
L=2.5e6;
Cp=1004;
g=9.8;
qvc_new=r*.622*es(Tc_new)/(P-es(Tc_new));
y=Tc_new-Tc_old+g/Cp*wc*dt+L/Cp*(qvc_new-qvc_old)...
 -lambda*(Te-Tc_old+L/Cp*(qve-qvc_old))*wc*dt;
