if (imodify==3)

disp('Please specify a starting velocity model...')
disp('Choose from the following:')
disp('return= begin from an arbitrary initial model or');
disp('1=      to begin from existing velocity model file')
init=input('--->');
if (sum(init)==0); init=0; end;
if (init==1);
  iread=0;
  while (iread~=1)
    disp('Velocity file must have an extension ".vmod"')
    modfile=input('Type name of velocity file without the extension >','s');
    fid=fopen([modfile '.vmod'],'r');
    if (fid<0)
      disp('Sorry file does not exist');
    else
      fclose(fid);
      eval(['load ' modfile '.vmod']);
      iread=1;
    end
  end
  if (sum(size(str2num(modfile(1))))>0)
    eval(['data=X' modfile ';']);
  else
    eval(['data=' modfile ';']);
  end
  vees=data(:,2)';
  zees=data(:,1)';
  numlayer=length(zees)-1;
  disp([modfile ' has ' num2str(numlayer) ' layers.']);
  disp('~~~~~~~~~~~~~~~~~~~~');
else
  disp('Enter number of velocity layers you wish to model (default is 3)');
  numlayer=input('--->');
  if (sum(numlayer)==0); numlayer=3; end
  zees=[0:numlayer].*5;
  zees(2)=zees(2)./2;
  vees=1000+113.3333*zees;
end
end

eval(fig);
vm=vees; zm=zees;

if (imodify~=3)
  if (imodify==1);
    vees=vees_save;
    zees=zees_save;
    subplot(121); hold off;
  end
  subplot(121); 
  plot(vees,zees,'-o')
  vt=vees+0.02*vmax; zt=zees-0.02.*ymax;
  for k=1:numlayer+1;
    text(vt(k),zt(k),num2str(k));
  end

  grid on
  axis([0 vmax -0.1 ymax]);
  xlabel ('Velocity (km/s)'); ylabel ('Depth (m)')
  set(gca, 'ydir', 'reverse')
  set(gca,'xtick',[0:dv:vmax],'xticklabel',[0:dv:vmax]/1000)
  [i1]=input('Enter 1st node to modify  >');
  [i2]=input('Enter last node to modify  >');
  i1=max([i1 1]); i2=min([i2 numlayer+1]);
  for i=i1:i2;
    disp(['Click to mark velocity node:  ' num2str(i)])
    [vm(i),zm(i)]=ginput(1);
    subplot(121); hold on;
    plot(vm(i),zm(i),'*')
  end;
  for i=i1:i2;
      if (i>1)
      zm(i)=max([zm(i) zm(i-1)+1e-4.*ymax]);	%subsequent layers must be deeper
    end 
    if (i<numlayer+1);
      zm(i)=min([zm(i) zm(i+1)-1e-4.*ymax]);	%subsequent layers must be deeper
    end    
  end
    
end

vees_save=vees; zees_save=zees;
vees=vm;  zees=zm;
zees(1)=0;
subplot(121); hold off
plot(vees,zees,'-o'); hold on;
vt=vees+0.02*vmax; zt=zees-0.02.*ymax;
for k=1:numlayer+1;
  text(vt(k),zt(k),num2str(k));
end
grid on
axis([0 vmax -0.1 ymax]);
xlabel ('Velocity (km/s)'); ylabel ('Depth (m)')
set(gca, 'ydir', 'reverse')
set(gca, 'xtick',[0:dv:vmax],'xticklabel',[0:dv:vmax]/1000)
hold on;