function conditions = getvapor(p,T,p0,T0,satmax,icefracFinal,TFinal,loadopt,fracload,Tunload)

g=9.8;
cpd=1005;
cpv=1870;
Lv=2.5e6;
Lf=3.34e5;
Ls=Lv+Lf;
R1=461.5;
R0=287;
ep=R0/R1;

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));

r0=ep*es(T0)/(p0-es(T0)); %cloud base water vapor mixing ratio
satfac=dfrac(p,p0,satmax); %saturation ratio

[liqfrac,dliqfrac]=lf_func(T,icefracFinal,TFinal);
icefrac=1-liqfrac;

pF1=es(T).*satfac; %vapor prsesure
rF1=ep.*pF1./(p-pF1); %vapor mixing ratio
rF=r0-rF1; %liquid + ice mixing ratio

lo_int=floor(loadopt);

switch lo_int
    case 0 %unload constant starting at Tunload
        rF(T<=Tunload)=rF(T<=Tunload)*fracload;
    case 1 %unload increasing linearly starting at Tunload. total unload at T1
        T1=(loadopt-lo_int)*1000;
        
        if Tunload==T1
            rF(T<=Tunload)=rF(T<=Tunload).*fracload;
        else
            rF(T<=Tunload & T>=T1)=rF(T<=Tunload & T>=T1).*(fracload-(Tunload-T(T<=Tunload & T>=T1))./(Tunload-T1)*fracload);
            rF(T<T1)=rF(T<T1)*fracload;
        end
    case 2 %unload with sudden increase at freezing level. Increase amount ffu
        ffu=loadopt-lo_int;
        rF(T<=Tunload)=rF(T<=Tunload)*fracload;
        if Tunload>273.15-5 && fracload>0
            rF(T<=273.15+TFinal)=rF(T<=273.15+TFinal)/fracload*max(fracload-ffu,0);
        end
    case 3 %unload based on pressure
        rF(p<=Tunload)=rF(p<=Tunload)*fracload;
    case 4 %unload based on pressure gradually
        T1=(loadopt-lo_int)*100000;
        
        if Tunload==T1
            rF(p<=Tunload)=rF(p<=Tunload).*fracload;
        else
            rF(p<=Tunload & p>=T1)=rF(p<=Tunload & p>=T1).*(1-(Tunload-p(p<=Tunload & p>=T1))./(Tunload-T1)*(1-fracload));
            rF(p<T1)=rF(p<T1)*fracload;
        end

    case 5 %unload based on condensate threshold
        rmax = (loadopt-lo_int)*100;
        rmax=rmax*1e-3;
        rF(rF>=rmax)=rmax;               
        inds=T<=-4+273.15 & T>-40+273.15;
        rF(inds)=(rF(inds)-fracload*rF(inds))/36.*(T(inds)-273.15+40)+fracload*rF(inds);
        rF(T<-40+273.15)=fracload*rF(T<-40+273.15);       
        
end

if rF<0
    error('Negative Condensate')
end
    
conditions={liqfrac icefrac dliqfrac rF1 rF pF1};