% function parabolicw
function [W,z,dt]=parabolicw(MAXW,topkm)
% always compute to 7 km
dz=8;                                      % delta z
nstep=6*1000/dz;  
z=1000:dz:7000;

% hard code for 5 km or 7 km tops
if (topkm==7) 
factor=9;
root=7;
else 
factor=6.25;
root=6;
end;

topm=topkm*1000;
W(1)=0;

for i=2:nstep+1;
% kludge for top of 5 km profile
if ((z(i)>3500) & (topkm==5))
W(i)=-1*(MAXW-.5)*(z(i)/1000-2)*(z(i)/1000-5)/2.25+.5;
else 
W(i)=-1*(MAXW-.5)*(z(i)/1000-1)*(z(i)/1000-root)/factor+.5;
end
if (W(i)<=0) 
W(i)=0.0; 
end
if (z(i)>topm) 
W(i)=0.0; 
end;
if (W(i)>0)
dt(i)=dz/W(i);
else 
dt(i)=0;
end;

end; % for loop

