
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  map.m : to manipulate kmt- maps for MOM 2
%
% NOTE : This matlab script manipulates your kmt.dta
%        and your delta.kmt.'imt'x'jmt'x'km'.h file,
%        if it's able to get them. Unless there is a 
%        undo function, be careful what you do.
%
%        To obtain the information of the models axis the user
%        has to supply the file grid.dta (in 32 bit ieee format).
%
%        To view another (different) kmt-map you will need the files 
%        kmt.dta and grid.dta with an additional suffix.
%        These files will not be change.
%
%        Please copy this script to your $HOME/matlab 
%        directory. Then you are able to manipulate the 
%        kmt-map related files in a directory of your choice.
%       
%     All the comments and complains to ceden@ifm.uni-kiel.de
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


fprintf('\n\n ===== map.m ===== \n\n');

clear
clear global

% USER INPUT

 machine_format = 'ieee-be';

% Choose the form the following for machine_format
%
%        'native'      or 'n' - local machine format - the default
%        'ieee-le'     or 'l' - IEEE floating point with little-endian
%                               byte ordering
%        'ieee-be'     or 'b' - IEEE floating point with big-endian
%                               byte ordering
%        'vaxd'        or 'd' - VAX D floating point and VAX ordering
%        'vaxg'        or 'g' - VAX G floating point and VAX ordering
%        'cray'        or 'c' - Cray floating point with big-endian
%                               byte ordering
%        'ieee-le.l64' or 'a' - IEEE floating point with little-endian
%                               byte ordering and 64 bit long data type
%                               (currently supported on: alpha, sgi, ibm_rs,
%                                hp300, vax, and cray platforms only)
% END USER INPUT

%--------------------------------------------------------
% get the kmt-map from unformatted fortran written files
%--------------------------------------------------------

fkmt=fopen('kmt.dta','r',machine_format);

if fkmt == -1 
 fprintf(' Error: the file kmt.dta isn''t in working directory\n');
 error(' ==> map.m');
end;

fortran_format1=fread(fkmt,1,'long');

% a first check if the proper machine_format was chosen
if fortran_format1 ~=80 & fortran_format1~=172
 fprintf('\n         The variable machine_format seems to be not correct!\n');
 fprintf('         change this in the USER INPUT section directly.\n');
 fprintf('\n');
 error(' ==> map.m');
end;

if fortran_format1 == 80
 % read AWI security check
 AWI_ID = fread(fkmt,80,'char');
 comstr=setstr(AWI_ID)';
 fprintf(' AWI security check from kmt.dta : \n  %s\n',comstr(1:50));
 fprintf(' please check, if this seems correct\n');
 fread(fkmt,2,'long');
 stamp=fread(fkmt,32,'char');
 iotext=fread(fkmt,80,'char');
 expnam=fread(fkmt,60,'char');
 fread(fkmt,1,'long');
else
 stamp=fread(fkmt,32,'char');
 iotext=fread(fkmt,80,'char');
 expnam=fread(fkmt,60,'char');
 fread(fkmt,1,'long');
end;

fprintf(' stamp = %s\n',setstr(stamp));
fprintf(' expnam = %s\n',setstr(expnam));

fread(fkmt,1,'long');
imt=fread(fkmt,1,'long');
jmt=fread(fkmt,1,'long');
km=fread(fkmt,1,'long');
fread(fkmt,1,'long');

fread(fkmt,1,'long');
stamp2=fread(fkmt,32,'char');
iotext2=fread(fkmt,80,'char');
expnam2=fread(fkmt,60,'char');
fread(fkmt,1,'long');

fread(fkmt,1,'long');
kmt=fread(fkmt,[imt jmt],'long')';
fread(fkmt,1,'long');

fclose(fkmt);
clear fkmt

%-----------------------------------------------
% get the axis
%-----------------------------------------------

fgrid=fopen('grid.dta','r',machine_format);

if fgrid == -1 
 fprintf(' Error: the file grid.dta isn''t in working directory\n');
 error(' ==> map.m');
end;

if fread(fgrid,1,'long') == 80
 % read AWI security check
 comstr = fread(fgrid,80,'char');
 fprintf(' AWI security check from grid.dta :\n  %s \n',comstr(1:50));
 fprintf(' please check, if this seems correct\n');
 fread(fgrid,2,'long');
end;

imt2=fread(fgrid,1,'long');
jmt2=fread(fgrid,1,'long');
km2=fread(fgrid,1,'long');
fread(fgrid,1,'long');

% another security check
if imt2 ~= imt | jmt2 ~= jmt | km2 ~= km 
 error(' grid.dta does not match file kmt.dta');
end;

xu=ones(imt,4);
yu=ones(jmt,4);
zw=ones(km,4);

fread(fgrid,1,'long');
xu(:,4)=fread(fgrid,imt,'float');  % dxtdeg
yu(:,4)=fread(fgrid,jmt,'float');  % dytdeg
xu(:,2)=fread(fgrid,imt,'float');  % dxudeg
yu(:,2)=fread(fgrid,jmt,'float');  % dyudeg
zw(:,4)=fread(fgrid,km,'float');   % dzt
fread(fgrid,1,'float');            % dzw(0)   
zw(:,2)=fread(fgrid,km,'float');   % dzw(1:km)
xu(:,3)=fread(fgrid,imt,'float');  % xt
xu(:,1)=fread(fgrid,imt,'float');  % xu
yu(:,3)=fread(fgrid,jmt,'float');  % yt
yu(:,1)=fread(fgrid,jmt,'float');  % yu
zw(:,3)=fread(fgrid,km,'float');   % zt
zw(:,1)=fread(fgrid,km,'float');   % zw
fread(fgrid,1,'long');

fclose(fgrid);
clear fgrid

zt=[0; zw(:,1)/100];
clear zw

%-----------------------------------------------
% some initialisations of variables and arrays
%-----------------------------------------------

% now construct delta.kmt.imtxjmtxkm.h name

delta_name=sprintf('delta.kmt.%ix%ix%i.h',imt,jmt,km);

% check if file is available

checkfile=delta_name;
fid=fopen(checkfile,'r');
if fid == -1 
 fprintf([' Warning: the file ' checkfile ' isn''t in working directory\n']);
 fprintf(' map.m will generate a new one.');
 fid=fopen(checkfile,'w');
 fwrite(fid,'c this file is generated by map.m !!!\n','char');
 fwrite(fid,'c this seems to be not correct \n','char');
end;
fclose(fid);

% format-string to write in delta.kmt.imtxjmtxkm.h

form_del='      kmt( %3i, %3i) =    %2i\n';

% Zoom variables

js=1;
je=jmt;
is=1;
ie=imt;

dzoom = 1;

% default increment value for manipulating kmt- map

dkmt = 1;

% Axis for figures

x=1:imt;
y=1:jmt;
X=ones(jmt,1)*x;
Y=(y')*ones(1,imt);

% Graphical stuff

fig1=figure(1); clf; 
%set(fig1,'position',[10 10 600 600]);

fig2=figure(2); clf
set(fig2,'visible','off');
%set(fig2,'position',[10 10 500 500]);

fig3=figure(3); clf
set(fig3,'visible','off');

fig4=figure(4); clf
set(fig4,'visible','off');

cmap=jet;
[m n]=size(cmap);
cmap=cmap(m:-1:1,:);
clear m n

zoom off

% angles for 3D-view

angle1=0;
angle2=90;

% some logical variables

stop=1;
newplot=1;

fprintf([' writing changes to ' delta_name '\n']);
fnewh=fopen(delta_name,'a');
fprintf(fnewh,'c\n');
fprintf(fnewh,'c map.m : writing interactive changes \n');
fprintf('\n');

% initialize some arrays for undo - function

% the undo buffer for kmt-map
unkmt=kmt;
% the undo buffer for delta.kmt.xx.xx.xx.h file
unfh='c\n';
% the undo buffer for floodfill
unfhmap=[];

% initialize some variables for second kmt-map to compare with

sec_map=0;
xu_sec=[];
yu_sec=[];
zt_sec=[];
imt_sec=[];
jmt_sec=[];
km_sec=[];
kmt_sec=[];
x_sec=[];
y_sec=[];
X_sec=[];
Y_sec=[];
stamp_sec=[];
expnam_sec=[];

%--------------------------------------
% Start main loop 
%--------------------------------------  

fprintf('now starting main loop\n');
fprintf('Move cursor to graphical window, click on the buttons\n');
fprintf('and you will see what happens. \n');
fprintf(' - Never stop this script by ctrl c, because of the undo buffers!  -\n');
fprintf(' - If the script stopped with an error message,                    -\n');
fprintf(' - please remove kmt and delta file and start again.               -\n');

while stop == 1

 figure(fig1);
 if sec_map == 0
  set(fig3,'visible','off');
 else
  set(fig3,'visible','on');
 end;

 %-------------------------------------------
 % Contour the map, but only when necessary 
 %-------------------------------------------

 if newplot == 1

  % main plotting part
  % plot zmap not the kmt-map
  fprintf(' generating new map');

  if sec_map

   % get zmap of second kmt-map

   ij_back=[is ie js je];
   is=min(imt_sec,is);
   ie=min(imt_sec,ie);
   js=min(jmt_sec,js);
   je=min(jmt_sec,je);

   if je>js & ie>is

    % plot only if zoom area is valid

    zmap_sec=zeros(je-js+1,ie-is+1);
    for j=js:je
     for i=is:ie
      zmap_sec(j-js+1,i-is+1)=zt_sec(kmt_sec(j,i)+1);
     end;
    end;

   end;

   is=ij_back(1); ie=ij_back(2); js=ij_back(3); je=ij_back(4);
  end;

  zmap=zeros(je-js+1,ie-is+1);
  %for j=js:je
  % for i=is:ie
  %    zmap(j-js+1,i-is+1)=zt(kmt(j,i)+1);
  % end;
  %end;
  kmt_tmp=kmt(js:je,is:ie);
  for k=1:km
     fprintf('.');
     n_tmp=find( kmt_tmp == k);
     zmap(n_tmp)=ones(length(n_tmp),1)*zt(k+1);
  end;
  clear kmt_tmp
  fprintf(' done\n');

  % define color scale

  if exist('zmap_sec') == 1
   mm1=max([zmap(:); zmap_sec(:)]);
   mm2=min([zmap(:); zmap_sec(:)]);
  else
   mm1=max(zmap(:));
   mm2=min(zmap(:));
  end;

  ih = image(x(is:ie),y(js:je),floor((zmap-mm2)/(mm1-mm2)*63)+1);
  set(get(ih,'parent'),'ydir','normal');
  hold on
  if min(min(kmt(js:je,is:ie))) == 0
   contour(x(is:ie),y(js:je),kmt(js:je,is:ie),[9999 0.001],'k');
  end;
  %plot(X(js:je,is:ie),Y(js:je,is:ie),'k.','markersize',3);
  hold off
  xlabel('i -->')
  ylabel('j -->')
  colormap(cmap)
  ch=colorbar;
  set(ch,'YTick',[10 20 30 40 50 60]);
  set(ch,'YTickLabel', ...
      num2str([10:10:60]'/64*(mm1-mm2)+mm2,'%8.2f m')); 

  if sec_map

   % plot second kmt-map for comparision
 
   fig_back=gcf;
   figure(fig3);
   ij_back=[is ie js je];
   is=min(imt_sec,is);
   ie=min(imt_sec,ie);
   js=min(jmt_sec,js);
   je=min(jmt_sec,je);
   
   if je>js & ie>is

    % plot only if zoom area is valid
  
    ih=image(x_sec(is:ie),y_sec(js:je), ...
        floor((zmap_sec-mm2)/(mm1-mm2)*63)+1);
    set(get(ih,'parent'),'ydir','normal');
    hold on
    if min(min(kmt_sec(js:je,is:ie))) == 0
     contour(x_sec(is:ie),y_sec(js:je),kmt_sec(js:je,is:ie),[9999 0.001],'k');
    end;
    plot(X_sec(js:je,is:ie),Y_sec(js:je,is:ie),'k.','markersize',3);
    hold off
    xlabel('i -->')
    ylabel('j -->')
    colormap(cmap)
    ch=colorbar;
    set(ch,'YTick',[10 20 30 40 50 60]);
    set(ch,'YTickLabel', ...
      num2str([10:10:60]'/64*(mm1-mm2)+mm2,'%8.2f m')); 
    
    title(' second kmt-map')
   else
    clf
    title(' zoom area is not valid')
   end;

   figure(fig_back);
   is=ij_back(1); ie=ij_back(2); js=ij_back(3); je=ij_back(4);
   clear ij_back fig_back

  end;
  clear zmap zmap_sec
  newplot=0;
 end;

 %-----------------------
 % Wait for user input
 %-----------------------

 ui1=uicontrol('style','pushbutton');
 set(ui1,'units','normalized');
 set(ui1,'position',[0.005 0.94 0.1 0.05])
 set(ui1,'callback','choice=1;thisstop=0;');
 set(ui1,'string','Exit')

 ui2=uicontrol('style','pushbutton');
 set(ui2,'units','normalized');
 set(ui2,'position',[0.15 0.94 0.1 0.05])
 set(ui2,'callback','choice=2;thisstop=0;');
 set(ui2,'string','zoom in');

 ui3=uicontrol('style','pushbutton');
 set(ui3,'units','normalized');
 set(ui3,'position',[0.3 0.94 0.1 0.05])
 set(ui3,'callback','choice=3;thisstop=0;');
 set(ui3,'string','zoom out');

 ui4=uicontrol('style','pushbutton');
 set(ui4,'units','normalized');
 set(ui4,'position',[0.45 0.94 0.1 0.05])
 set(ui4,'callback','choice=4;thisstop=0;');
 set(ui4,'string','3-D View');

 ui5=uicontrol('style','pushbutton');
 set(ui5,'units','normalized');
 set(ui5,'position',[0.6 0.94 0.1 0.05])
 set(ui5,'callback','choice=5;thisstop=0;');
 set(ui5,'string','floodfill');

 ui6=uicontrol('style','pushbutton');
 set(ui6,'units','normalized');
 set(ui6,'position',[0.75 0.94 0.1 0.05])
 set(ui6,'callback','choice=6;thisstop=0;');
 set(ui6,'string','single pt.');

 ui7=uicontrol('style','pushbutton');
 set(ui7,'units','normalized');
 set(ui7,'position',[0.895 0.94 0.1 0.05])
 set(ui7,'callback','choice=7;thisstop=0;');
 set(ui7,'string','undo');

 ui8=uicontrol('style','text');
 set(ui8,'units','normalized');
 set(ui8,'position',[0.2 0.01 0.6 0.03])
 set(ui8,'string','');

 ui9=uicontrol('style','pushbutton');
 set(ui9,'units','normalized');
 set(ui9,'position',[0.005 0.2 0.09 0.05])
 set(ui9,'callback','choice=9;thisstop=0;');
 set(ui9,'string','potholes');

 ui10=uicontrol('style','pushbutton');
 set(ui10,'units','normalized');
 set(ui10,'position',[0.005 0.1 0.09 0.05])
 set(ui10,'callback','choice=10;thisstop=0;');
 set(ui10,'string','non adv.');

 ui11=uicontrol('style','pushbutton');
 set(ui11,'units','normalized');
 set(ui11,'position',[0.005 0.3 0.09 0.05])
 set(ui11,'callback','choice=11;thisstop=0;');
 set(ui11,'string','trenches');

 ui12=uicontrol('style','pushbutton');
 set(ui12,'units','normalized');
 set(ui12,'position',[0.005 0.4 0.09 0.05])
 set(ui12,'callback','choice=12;thisstop=0;');
 set(ui12,'string','killworth');

 ui13=uicontrol('style','pushbutton');
 set(ui13,'units','normalized');
 set(ui13,'position',[0.005 0.6 0.09 0.05])
 set(ui13,'callback','choice=13;thisstop=0;');
 set(ui13,'string','sec. kmt');

 if sec_map
  figure(fig3);
  ui14=uicontrol('style','text');
  set(ui14,'units','normalized');
  set(ui14,'position',[0.2 0.01 0.6 0.03])
  set(ui14,'string','');
  figure(fig1);
 end;

 ui15=uicontrol('style','pushbutton');
 set(ui15,'units','normalized');
 set(ui15,'position',[0.005 0.7 0.09 0.05])
 set(ui15,'callback','choice=14;thisstop=0;');
 set(ui15,'string','kmu map');

 ui16=uicontrol('style','pushbutton');
 set(ui16,'units','normalized');
 set(ui16,'position',[0.005 0.8 0.09 0.05])
 set(ui16,'callback','choice=16;thisstop=0;');
 set(ui16,'string','more pts');

 pointer=get(gcf,'pointer');
 set(gcf,'pointer','crosshair');

 choice = -1;
 thisstop=1;
 while thisstop == 1
     drawnow;

     % Give coordinates and depth of pointer selection

     p=get(gca,'Currentpoint');
     i=round(p(1,1));
     j=round(p(1,2));
     i=max(is,i);
     j=max(js,j);
     i=min(ie,i);
     j=min(je,j);
     lat=yu(j,3);
     lon=xu(i,3);
     while lon > 360
      lon=lon-360;
     end;
     if lon > 180
      lon=lon-360;
     end;
     if lon < 0
      set(ui8,'string', ...
       sprintf('kmt(%i,%i)=%i => %6.1f m (%6.2fN, %6.2fW)', ...
        i,j,kmt(j,i),zt(kmt(j,i)+1),lat,-lon));
     else
      set(ui8,'string', ...
       sprintf('kmt(%i,%i)=%i => %6.1f m (%6.2fN, %6.2fE)', ...
        i,j,kmt(j,i),zt(kmt(j,i)+1),lat,lon));
     end;

     if sec_map

      % Also in second map

      i=max(is,i);
      j=max(js,j);
      i=min([ie i imt_sec]);
      j=min([je j jmt_sec]);
      lat=yu_sec(j,3);
      lon=xu_sec(i,3);
      while lon > 360
       lon=lon-360;
      end;
      if lon > 180
       lon=lon-360;
      end;
      if lon < 0
       set(ui14,'string', ...
        sprintf('kmt(%i,%i)=%i => %6.1f m (%6.2fN, %6.2fW)', ...
         i,j,kmt_sec(j,i),zt_sec(kmt_sec(j,i)+1),lat,-lon));
      else
       set(ui14,'string', ...
        sprintf('kmt(%i,%i)=%i => %6.1f m (%6.2fN, %6.2fE)', ...
         i,j,kmt_sec(j,i),zt_sec(kmt_sec(j,i)+1),lat,lon));
      end;
     end;
 end; 

 set(gcf,'pointer',pointer);
 set(ui1,'visible','off');
 set(ui2,'visible','off');
 set(ui3,'visible','off');
 set(ui4,'visible','off');
 set(ui5,'visible','off');
 set(ui6,'visible','off');
 set(ui7,'visible','off');
 set(ui8,'visible','off');
 set(ui9,'visible','off');
 set(ui10,'visible','off');
 set(ui11,'visible','off');
 set(ui12,'visible','off');
 set(ui13,'visible','off');
 if sec_map
  set(ui14,'visible','off');
 end;    
 set(ui15,'visible','off');
 set(ui16,'visible','off');
 clear ui1 ui2 ui3 ui4 ui5 ui6 ui7 ui8 pointer 
 clear ui9 ui10 ui11 ui12 ui13 ui14 ui15 ui16

 %-------------------------------------
 % Now do what the user want
 %-------------------------------------

 if choice == 14

  %------------------------------------
  % show kmu-map in another window
  %------------------------------------


  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string',' kmu-map is in fourth graphical window');
  drawnow;

  set(fig4,'visible','on');
  figure(fig4);
  clf

  kmu=zeros(jmt,imt);
  for j=max(2,js):min(je,jmt-1)
   for  i=max(2,is):min( ie,imt-1)
    kmu(j,i) = min([kmt(j,i) kmt(j,i+1) kmt(j+1,i) kmt(j+1,i+1)]);
   end;
  end;

  zmap=zeros(je-js+1,ie-is+1);
  for j=js:je
   for i=is:ie
      zmap(j-js+1,i-is+1)=zt(kmu(j,i)+1);
   end;
  end;

  ih = image( x(is:ie),y(js:je),floor((zmap-mm2)/(mm1-mm2)*63)+1);
  set(get(ih,'parent'),'ydir','normal');
  hold on
  if min(min(kmu(js:je,is:ie))) == 0
   contour(x(is:ie),y(js:je),kmu(js:je,is:ie),[9999 0.001],'k');
  end;
  %plot(X(js:je,is:ie),Y(js:je,is:ie),'k.','markersize',3);
  hold off
  xlabel('i -->')
  ylabel('j -->')
  colormap(cmap)
  ch=colorbar;
  set(ch,'YTick',[10 20 30 40 50 60]);
  set(ch,'YTickLabel', ...
      num2str([10:10:60]'/64*(mm1-mm2)+mm2,'%8.2f m')); 

  ui3=uicontrol('style','pushbutton');
  set(ui3,'units','normalized');
  set(ui3,'position',[0.05 0.01 0.05 0.05])
  set(ui3,'callback','thisstop=0;');
  set(ui3,'string','Ok');

  ui8=uicontrol('style','text');
  set(ui8,'units','normalized');
  set(ui8,'position',[0.2 0.01 0.6 0.03])
  set(ui8,'string','');

  thisstop=1;
  while thisstop == 1

   drawnow;
   p=get(gca,'Currentpoint');
   i=round(p(1,1));
   j=round(p(1,2));
   i=max(is,i);
   j=max(js,j);
   i=min(ie,i);
   j=min(je,j);
   lat=yu(j,1);
   lon=xu(i,1);
   while lon > 360
      lon=lon-360;
   end;
   if lon > 180
      lon=lon-360;
   end;
   if lon < 0
      set(ui8,'string', ...
       sprintf('kmu(%i,%i)=%i => %6.1f m (%6.2fN, %6.2fW)', ...
        i,j,kmu(j,i),zt(kmu(j,i)+1),lat,-lon));
   else
      set(ui8,'string', ...
       sprintf('kmu(%i,%i)=%i => %6.1f m (%6.2fN, %6.2fE)', ...
        i,j,kmu(j,i),zt(kmu(j,i)+1),lat,lon));
   end;

  end;

  set(ui3,'visible','off');
  set(ui8,'visible','off');
  clear kmu zmap ui3 ui8

  figure(fig1);
  set(fig4,'visible','off');

  set(ui0,'visible','off');
  clear ui0

 elseif choice == 13

  %--------------------------------------
  % load another kmt-map to compare with
  %--------------------------------------


  ui1=uicontrol('style','pushbutton','units','normalized',...
    'position',[0.1 0.5 0.8 .05], 'callback','thisstop=0;','string', ...
    'give suffix of second kmt- and grid-file (click when ready)');
  ui2=uicontrol('style','edit','units','normalized',...
           'position',[0.3 0.45 0.4 .05], 'max',1, ...
           'string','   ');

  thisstop=1;
  while thisstop == 1
     drawnow;
  end;

  comstr=deblank(get(ui2,'string'));
  while length(comstr)>1 & comstr(1) == ' ' 
   comstr=comstr(2:length(comstr));
  end;

  set(ui1,'visible','off');
  set(ui2,'visible','off');
  clear ui1 ui2

  % check if files are in working directory

  all_there=1;

  checkfile=['kmt.dta.' comstr];
  fid=fopen(checkfile,'r');
  if fid == -1, all_there=0; fprintf([' file ' checkfile ' not found\n']); 
  else, fclose(fid); end;
  checkfile=['grid.dta.' comstr];
  fid=fopen(checkfile,'r');
  if fid == -1, all_there=0; fprintf([' file ' checkfile ' not found\n']); 
  else, fclose(fid); end;

  if all_there

   % get the kmt-map

   fkmt=fopen(['kmt.dta.' comstr],'r',machine_format);

   format_check = fread(fkmt,1,'long');

   % a first check if the proper machine_format was chosen
   if format_check ~=80 & format_check~=172
    fprintf('\n         The variable machine_format seems to be not correct\n');
    fprintf('         for this file !! \n');
    fprintf('\n');
    error(' ==> map.m');
   end;

   if format_check == 80
    % read AWI security check
    AWI_ID_sec = setstr(fread(fkmt,80,'char'))';
    fprintf([' AWI security check from kmt.dta.' comstr ':\n  %s\n'], ...
         AWI_ID_sec(1:50));
    fprintf(' please check, if this seems correct\n');
    fread(fkmt,2,'long');
   end;

   stamp_sec=fread(fkmt,32,'char');
   iotext_sec=fread(fkmt,80,'char');
   expnam_sec=fread(fkmt,60,'char');
   fread(fkmt,1,'long');

   fread(fkmt,1,'long');
   imt_sec=fread(fkmt,1,'long');
   jmt_sec=fread(fkmt,1,'long');
   km_sec=fread(fkmt,1,'long');
   fread(fkmt,1,'long');

   fread(fkmt,1,'long');
   stamp_sec=setstr(fread(fkmt,32,'char'))';
   fread(fkmt,80,'char');
   expnam_sec=setstr(fread(fkmt,60,'char'))';
   fread(fkmt,1,'long');

   fread(fkmt,1,'long');
   for j=1:jmt_sec
    kmt_sec(j,:)=fread(fkmt,imt_sec,'long')';
   end;
   fread(fkmt,1,'long');

   fprintf(' stamp of second kmt-map = %s\n',stamp_sec);
   fprintf(' expnam of second kmt-map = %s\n',expnam_sec);

   fclose(fkmt);
   clear fkmt

   % get the axis

   fgrid=fopen(['grid.dta.' comstr],'r',machine_format);

   if fread(fgrid,1,'long') == 80
    % read AWI security check
    AWI_ID_sec = setstr(fread(fgrid,80,'char'))';
    fprintf([' AWI security check from grid.dta.' comstr ':\n  %s\n'], ...
         AWI_ID_sec(1:50));
    fprintf(' please check, if this seems correct\n');
    fread(fgrid,2,'long');
   end;

   imt2=fread(fgrid,1,'long');
   jmt2=fread(fgrid,1,'long');
   km2=fread(fgrid,1,'long');
   fread(fgrid,1,'long');

   % another security check
   if imt2 ~= imt_sec | jmt2 ~= jmt_sec | km2 ~= km_sec 
    error([' grid.dta.' comstr ' does not match file kmt.dta.' comstr]);
   end;

   xu_sec=ones(imt_sec,4); 
   yu_sec=ones(jmt_sec,4);
   zw_sec=ones(km_sec,4);

   fread(fgrid,1,'long');
   xu_sec(:,4)=fread(fgrid,imt_sec,'float');  % dxtdeg
   yu_sec(:,4)=fread(fgrid,jmt_sec,'float');  % dytdeg
   xu_sec(:,2)=fread(fgrid,imt_sec,'float');  % dxudeg
   yu_sec(:,2)=fread(fgrid,jmt_sec,'float');  % dyudeg
   zw_sec(:,4)=fread(fgrid,km_sec,'float');   % dzt
   fread(fgrid,1,'float');            % dzw(0)   
   zw_sec(:,2)=fread(fgrid,km_sec,'float');   % dzw(1:km)
   xu_sec(:,3)=fread(fgrid,imt_sec,'float');  % xt
   xu_sec(:,1)=fread(fgrid,imt_sec,'float');  % xu
   yu_sec(:,3)=fread(fgrid,jmt_sec,'float');  % yt
   yu_sec(:,1)=fread(fgrid,jmt_sec,'float');  % yu
   zw_sec(:,3)=fread(fgrid,km_sec,'float');   % zt
   zw_sec(:,1)=fread(fgrid,km_sec,'float');   % zw
   fread(fgrid,1,'long');

   fclose(fgrid);
   clear fgrid

   zt_sec=[0; zw_sec(:,1)/100];
   clear zw_sec

   % Axis for figures

   x_sec=1:imt_sec;
   y_sec=1:jmt_sec;
   X_sec=ones(jmt_sec,1)*x_sec;
   Y_sec=(y_sec')*ones(1,imt_sec);

   % enable ploting of second map

   sec_map = 1;
   newplot=1;

  end;

  clear comstr all_there

 elseif choice == 12
  
  %----------------------------------
  % check for killworth criterium
  % search for topographic instabilities (based  on the  work of
  % Peter Killworth  ...  eqn 11 from Ocean Modeling, Nov 1987)
  %----------------------------------

  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string', 'check for killworth criterium ...');
  drawnow;
  fprintf(' check for killworth criterium (only in zoom area):\n');

  am=2.5e8;
  dtuv=3000;
  killcrit=zeros(jmt,imt);

  for j=max(2,js):min(je,jmt-1)

   cst=cos(yu(j,3)/180*pi);
   dyt=yu(j,4)*1852*60*100;
   dysq=dyt^2;

   for i=max(2,is):min(ie,imt-1)

     if kmt(j-1,i+1) ~= 0 & kmt(j,i+1) ~= 0 

       dxt=xu(i,4)*1852*60*100;
       dxsq = (dxt*cst)^2;
       clix = am*dtuv/dxsq;
       h1   = zt(kmt(j-1,i+1)+1)*100;
       h2   = zt(kmt(j,i+1)+1)*100;
       hx   = (8.0*h1*h2/(h1+h2)^2 + dxsq/dysq)/(4.0 + dxsq/dysq);

       if clix >= hx 
        fprintf('  criterium exceeded at kmt(%i,%i)=%i',i,j,kmt(j,i));
        fprintf(' clix=%f hx=%f\n',clix,hx);
       end;

       killcrit(j,i)=clix/hx*cst^2;
     end;
   end;
  end;
  fprintf(' done\n');

  set(ui0,'string', 'results are in second window ...');
  drawnow;

  figure(fig2);

  mm1_k=max(killcrit(:));
  mm2_k=min(killcrit(killcrit>0));
  n=find(killcrit == 0);
  killcrit(n)=ones(size(killcrit(n)))*mm2_k;
  ih=image(x(is:ie),y(js:je), ...
     floor((killcrit(js:je,is:ie)-mm2_k)/(mm1_k-mm2_k)*63)+1);
  set(get(ih,'parent'),'ydir','normal');
  hold on
  if min(min(kmt(js:je,is:ie))) == 0
   contour(x(is:ie),y(js:je),kmt(js:je,is:ie),[9999 0.001],'k');
  end;
  plot(X(js:je,is:ie),Y(js:je,is:ie),'k.','markersize',3);
  hold off
  xlabel('i -->')
  ylabel('j -->')
  colormap(cmap)
  ch=colorbar;
  set(ch,'YTick',[10 20 30 40 50 60]);
  set(ch,'YTickLabel', ...
      num2str([10:10:60]'/64*(mm1_k-mm2_k)+mm2_k,'%8.4f')); 

  title(' Killworth criterium')

  ui1=uicontrol('style','pushbutton');
  set(ui1,'units','normalized');
  set(ui1,'position',[0.05 0.05 0.05 .05]);     
  set(ui1,'string', 'ok');
  set(ui1,'callback','thisstop=0;');

  ui8=uicontrol('style','text');
  set(ui8,'units','normalized');
  set(ui8,'position',[0.2 0.01 0.6 0.03])
  set(ui8,'string','');

  thisstop=1;
  while thisstop == 1
   drawnow;
   p=get(gca,'Currentpoint');
   i=round(p(1,1));
   j=round(p(1,2));
   i=max(is,i);
   j=max(js,j);
   i=min(ie,i);
   j=min(je,j);
   set(ui8,'string', ...
   sprintf('kmt(%i,%i)=%i => killcrit=%6.4f ',i,j,kmt(j,i),killcrit(j,i)));
  end;

  set(ui1,'visible','off');
  set(ui8,'visible','off');
  clear ui1 ui8 mm1_k mm2_k
  clear killcrit
  
  figure(fig1);

  set(ui0,'visible','off');
  clear ui0

 elseif choice == 11

  %---------------------------------
  % check for corners and trenches
  %---------------------------------
 
  % mimic comes from j.dengg's subr tophole.F

%     Search for isolated corners and trenches, i.e. "t" grid boxes 
%     at any depth that can`t be influenced by advection
%     -------------------------------------------------------------

%      Examples: (x=land, .=sea)
%
%      1) x...x   2) .....   3) xxx..   4) x....
%         xx.xx      .....      xx..x      xx...
%         xx.xx      xx.xx      x..xx      x..xx
%         x...x      xxxxx      ..xxx      x....
%
%      Note: Case 4 is not covered by the present routine! However, it
%            is probably arguable if a diffusive flux should be 
%            possible through a gap like that.

%      This procedure mimics a check that will later be repeated for
%      the surface level in MOM (1)`s iperim.F.

  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string', 'check for isolated corners and trenches ...');
  drawnow;

  potholes=0;
  fprintf(' check for isolated corners and trenches (only in zoom area) :\n');

% Compute number of vertical levels on the "u" grid

  kmu=zeros(jmt,imt);
  for j=1: jmt-1
   for  i=1: imt-1
    kmu(j,i) = min([kmt(j,i) kmt(j,i+1) kmt(j+1,i) kmt(j+1,i+1)]);
   end;
  end; 

  for j=max(2,js): min(je, jmt-1)
    for i=max(2,is): min(ie,imt-1)
     k = kmt(j,i);
     if k ~= 0 & kmu(j,i) < k & kmu(j,i-1) < k & ...
          kmu(j-1,i) < k & kmu(j-1,i-1) < k

      fprintf('  found at kmt(%i,%i)=%i',i,j,kmt(j,i));

      tmp=0;
      if kmt(j-1,i) >= k tmp=tmp+1; end;
      if kmt(j+1,i) >= k tmp=tmp+1; end;
      if kmt(j,i+1) >= k tmp=tmp+1; end;
      if kmt(j,i-1) >= k tmp=tmp+1; end;

      if tmp == 1
       fprintf(' this is case 2 (just a pocket)');
      end;
 
      fprintf('\n');

      potholes=1;

     end;
   end;
  end;

  fprintf(' done\n');
  set(ui0,'visible','off');
  clear ui0

  if potholes

   ui0=uicontrol('style','pushbutton');
   set(ui0,'units','normalized');
   set(ui0,'position',[0.2 0.45 0.6 .1]);     
   set(ui0,'string', 'found some, look at text window');
   set(ui0,'callback','thisstop=0;');
 
   thisstop=1;
   while thisstop
    drawnow;
   end; 

   set(ui0,'visible','off');
   clear ui0

  end;

 elseif choice == 9

  %--------------------
  % check for potholes
  %--------------------

  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string', 'check for potholes ...');
  drawnow;

  potholes=0;
  fprintf(' check for potholes :\n');
  for i=2:imt-1
   for j=2:jmt-1
    if kmt(j,i)>kmt(j,i+1) & ...
       kmt(j,i)>kmt(j,i-1) & ...
       kmt(j,i)>kmt(j+1,i) & ...
       kmt(j,i)>kmt(j-1,i) 
     fprintf('  found pothole at kmt(%i,%i)=%i\n',i,j,kmt(j,i));
     potholes=1;
    end;
   end;
  end;
  fprintf(' done\n');
  set(ui0,'visible','off');
  clear ui0

  if potholes

   ui0=uicontrol('style','pushbutton');
   set(ui0,'units','normalized');
   set(ui0,'position',[0.2 0.45 0.6 .1]);     
   set(ui0,'string', 'found potholes, look at text window');
   set(ui0,'callback','thisstop=0;');
 
   thisstop=1;
   while thisstop
    drawnow;
   end; 

   set(ui0,'visible','off');
   clear ui0

  end;

 elseif choice == 10 

  %-------------------------------
  % check for non advective bays
  %-------------------------------

  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string', 'check for non advective bays');
  drawnow;

  potholes=0;
  fprintf(' check for non advective bays :\n');
  for i=2:imt-1
   for j=2:jmt-1
    kmu1=min( [kmt(j,i) kmt(j,i+1) kmt(j+1,i) kmt(j+1,i+1) ]);
    kmu2=min( [kmt(j-1,i) kmt(j-1,i+1) kmt(j,i) kmt(j,i+1) ]);
    kmu3=min( [kmt(j,i-1) kmt(j,i) kmt(j+1,i-1) kmt(j+1,i) ]);
    kmu4=min( [kmt(j-1,i-1) kmt(j-1,i) kmt(j,i-1) kmt(j,i) ]);

    if kmt(j,i)>0 & kmu1 == 0 & kmu2 == 0 & kmu3 == 0 & kmu4 == 0
     fprintf('  found non advective bay at kmt(%i,%i)=%i\n',i,j,kmt(j,i));
     potholes=1;
    end;
   end;
  end;
  fprintf(' done\n');
  set(ui0,'visible','off');
  clear ui0

  if potholes

   ui0=uicontrol('style','pushbutton');
   set(ui0,'units','normalized');
   set(ui0,'position',[0.2 0.45 0.6 .1]);     
   set(ui0,'string', 'found non advective bays, look at text window');
   set(ui0,'callback','thisstop=0;');
 
   thisstop=1;
   while thisstop
    drawnow;
   end; 

   set(ui0,'visible','off');
   clear ui0

  end;

 elseif choice == 16

  %--------------------------------
  % change depth on more points
  %--------------------------------

  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string', ...
    'click on points you want to change (middle butt. finishes)')

  % Wait for user input
  button=1;
  ti=1;
  clear xx yy ii jj
  hold on
  while button ~= 2
   [xx yy button]=ginput(1);
   if button ~=2
    ii(ti)=round(xx);
    jj(ti)=round(yy);
    if ii(ti)>=is & ii(ti)<=ie & jj(ti)>=js & jj(ti)<=je
     h_tmp=text(X(jj(ti),ii(ti)),Y(jj(ti),ii(ti)),'o');
     set(h_tmp,'HorizontalAlignment','center'); clear h_tmp
    end;
   end;
   ti=ti+1;
  end;
  hold off

  set(ui0,'visible','off');
  clear ui0

  isgood = 1;
  if exist('ii')
   for ti=1:length(ii)
    %if ii(ti)<is | ii(ti)>ie & jj(ti)<js & jj(ti)>je
    if ii(ti)<is & ii(ti)>ie & jj(ti)<js & jj(ti)>je
     isgood =0;
    end;
   end;
  else 
   isgood = 0;
  end;

  if isgood 

    tmp=kmt(jj(1),ii(1));

    ui1=uicontrol('style','slider');
    set(ui1,'units','normalized');
    set(ui1,'position',[0.05 0.2 0.05 0.6]);
    set(ui1,'callback', ...
     ['tmp=floor(get(ui1,''value'')*km);' ...
      'set(ui2,''string'',' ...
   'sprintf('' depth level %i => %6.1f m'',tmp,zt(tmp+1) ));']); 
    set(ui1,'value',0.5);

    ui2=uicontrol('style','text');
    set(ui2,'units','normalized');
    set(ui2,'position',[0.2 0.95 0.6 .05]);     
    set(ui2,'string', sprintf(' depth level %i => %6.1f m',tmp,zt(tmp+1)) );

    ui3=uicontrol('style','pushbutton');
    set(ui3,'units','normalized');
    set(ui3,'position',[0.2 0.0 0.2 .05]);     
    set(ui3,'string','accept');
    set(ui3,'callback','thisstop=0; thischoice=1;');
 
    ui4=uicontrol('style','pushbutton');
    set(ui4,'units','normalized');
    set(ui4,'position',[0.6 0.0 0.2 .05]);     
    set(ui4,'string','cancel');
    set(ui4,'callback','thisstop=0; thischoice=0;');

    thisstop=1;
    while thisstop == 1
     drawnow;
    end;

    set(ui1,'visible','off');
    set(ui2,'visible','off');
    set(ui3,'visible','off');
    set(ui4,'visible','off');
    clear ui1 ui2 ui3 ui4

    if thischoice == 1

      if tmp<0; tmp=0; end;
      if tmp>km; tmp=km; end;

      ui1=uicontrol('style','pushbutton','units','normalized',...
           'position',[0.2 0.5 0.6 .05], 'callback','thisstop=0;', ... 
           'string','some comment (click here when ready)');
      ui2=uicontrol('style','edit','units','normalized',...
           'position',[0.2 0.45 0.6 .05], 'max',1, ...
           'string','   ');

      thisstop=1;
      while thisstop == 1
       drawnow;
      end;

      comstr=get(ui2,'string');
      
      set(ui1,'visible','off');
      set(ui2,'visible','off');
      clear ui1 ui2

      % save old kmt-map
      unkmt=kmt;
      % flush text buffer for delta-file
      fprintf(fnewh,unfh);
      % flush binary buffer for delta file
      if ~isempty( unfhmap )
       [mm nn]=size( unfhmap );
       for mi=1:mm
        fprintf(fnewh,form_del,unfhmap(mi,1),unfhmap(mi,2),unfhmap(mi,3));
       end;
       fprintf(fnewh,'c map.m: done\n');
      end;
      unfhmap=[];

      for ti=1:length(ii)
       kmt(jj(ti),ii(ti))=tmp;
      end;

      % save change in text buffer for delta-file
      unfh=[sprintf('c map.m: change depth on more points\n')];
      unfh=[unfh sprintf([ 'c ' comstr '\n' ] ) ];
      clear comstr
      for ti=1:length(ii)
       unfh=[unfh sprintf(form_del,ii(ti),jj(ti),kmt(jj(ti),ii(ti)))];
      end;
      unfh=[unfh sprintf([ 'c done. \n' ] ) ];

      newplot=1;

    end;
  end;


 elseif choice == 6

  %--------------------------------
  % change depth of a single point
  %--------------------------------

  ui0=uicontrol('style','text');
  set(ui0,'units','normalized');
  set(ui0,'position',[0.2 0.95 0.6 .05]);     
  set(ui0,'string', ...
    'click on point on map, where you want to change the depth')

  % Wait for user input
  [xx yy button]=ginput(1);

  set(ui0,'visible','off');
  clear ui0

  ii=round(xx);
  jj=round(yy);

  if ii>=is & ii<=ie & jj>=js & jj<=je

    tmp=kmt(jj,ii);

    ui1=uicontrol('style','slider');
    set(ui1,'units','normalized');
    set(ui1,'position',[0.05 0.2 0.05 0.6]);
    set(ui1,'callback', ...
     ['tmp=floor(get(ui1,''value'')*km);' ...
      'set(ui2,''string'',' ...
   'sprintf('' previous value: kmt(%i,%i)=%i => %6.1f now %i => %6.1f'',' ...
      'ii,jj,kmt(jj,ii),zt(kmt(jj,ii)+1),tmp,zt(tmp+1) ));']); 
    set(ui1,'value',tmp/km);

    ui2=uicontrol('style','text');
    set(ui2,'units','normalized');
    set(ui2,'position',[0.2 0.95 0.6 .05]);     
    set(ui2,'string', ...
    sprintf('previous value: kmt(%i,%i)=%i => %6.1f m',ii,jj,tmp,zt(tmp+1)));

    ui3=uicontrol('style','pushbutton');
    set(ui3,'units','normalized');
    set(ui3,'position',[0.2 0.0 0.2 .05]);     
    set(ui3,'string','accept');
    set(ui3,'callback','thisstop=0; thischoice=1;');
 
    ui4=uicontrol('style','pushbutton');
    set(ui4,'units','normalized');
    set(ui4,'position',[0.6 0.0 0.2 .05]);     
    set(ui4,'string','cancel');
    set(ui4,'callback','thisstop=0; thischoice=0;');

    thisstop=1;
    while thisstop == 1
     drawnow;
    end;

    set(ui1,'visible','off');
    set(ui2,'visible','off');
    set(ui3,'visible','off');
    set(ui4,'visible','off');
    clear ui1 ui2 ui3 ui4

    if thischoice == 1

      if tmp<0; tmp=0; end;
      if tmp>km; tmp=km; end;

      ui1=uicontrol('style','pushbutton','units','normalized',...
           'position',[0.2 0.5 0.6 .05], 'callback','thisstop=0;', ... 
           'string','some comment (click here when ready)');
      ui2=uicontrol('style','edit','units','normalized',...
           'position',[0.2 0.45 0.6 .05], 'max',1, ...
           'string','   ');

      thisstop=1;
      while thisstop == 1
       drawnow;
      end;

      comstr=get(ui2,'string');
      
      set(ui1,'visible','off');
      set(ui2,'visible','off');
      clear ui1 ui2

      % save old kmt-map
      unkmt=kmt;
      % flush text buffer for delta-file
      fprintf(fnewh,unfh);
      % flush binary buffer for delta file
      if ~isempty( unfhmap )
       [mm nn]=size( unfhmap );
       for mi=1:mm
        fprintf(fnewh,form_del,unfhmap(mi,1),unfhmap(mi,2),unfhmap(mi,3));
       end;
       fprintf(fnewh,'c map.m: done\n');
      end;
      unfhmap=[];

      kmt(jj,ii)=tmp;
      % save change in text buffer for delta-file
      unfh=[sprintf('c map.m: change depth on a single point\n')];
      unfh=[unfh sprintf([ 'c ' comstr '\n' ] ) ];
      clear comstr
      unfh=[unfh sprintf(form_del,ii,jj,kmt(jj,ii))];
      newplot=1;

    end;
  end;

 elseif choice == 5

  %-----------------------------------------------
  %  floodfill of a closed region with landvalues
  %-----------------------------------------------

    ui0=uicontrol('style','text');
    set(ui0,'units','normalized');
    set(ui0,'position',[0.2 0.95 0.6 .05]);     
    set(ui0,'string','click on point on map, where you want to start')

    % Wait for user input
    [xx yy button]=ginput(1);

    set(ui0,'visible','off');
    clear ui0

    ii=round(xx);
    jj=round(yy);
    if ii>=is & ii<=ie & jj>=js & jj<=je

      ui1=uicontrol('style','pushbutton','units','normalized',...
        'position',[0.2 0.5 0.6 .05], 'callback','thisstop=0;', ... 
        'string','some comment (click here when ready)');
      ui2=uicontrol('style','edit','units','normalized',...
           'position',[0.2 0.45 0.6 .05], 'max',1, ...
           'string','   ');

      thisstop=1;
      while thisstop == 1
       drawnow;
      end;

      comstr=get(ui2,'string');
      
      set(ui1,'visible','off');
      set(ui2,'visible','off');
      clear ui1 ui2 

      ui0=uicontrol('style','text');
      set(ui0,'units','normalized');
      set(ui0,'position',[0.2 0.95 0.6 .05]);     
      set(ui0,'string','please wait a moment ...');
      set(ui0,'visible','on');
      drawnow;

      unkmt=kmt;
      fprintf(fnewh,unfh);
      if ~isempty( unfhmap )
       [mm nn]=size( unfhmap );
       for i=1:mm
        fprintf(fnewh,form_del,unfhmap(i,1),unfhmap(i,2),unfhmap(i,3)); 
       end;
       fprintf(fnewh,'c map.m: done\n');
      end;

      unfh=sprintf('c map.m: floodfill starting at i=%i, j=%i\n',ii,jj);
      unfh=[unfh sprintf(['c ' comstr '\n'])];
      clear comstr

      newplot=1;
      kmt(jj,ii)=-9;
      unfhmap=[ii jj 0];
      for i=(ii+1):ie
       if kmt(jj,i-1) == -9 & kmt(jj,i) > 0
        kmt(jj,i)=-9;
        unfhmap=[unfhmap; i jj 0];
       end;
      end;
      for i=(ii-1):-1:is
       if  kmt(jj,i+1) == -9 & kmt(jj,i) > 0
        kmt(jj,i)=-9;
        unfhmap=[unfhmap; i jj 0];
       end;
      end;
      percent=0;
      dpercent=1./(je-js+1)*100;
      fprintf(' floodfill started, %3i percent',percent);
      for j=(jj+1):je
        percent=percent+dpercent;
        fprintf(' \b\b\b\b\b\b\b\b\b\b\b\b%3i percent',round(percent));
        unfhmap_tmp=[];
        for i=(ii+1):(ie-1)
         if ( kmt(j-1,i-1) == -9 | kmt(j-1,i) == -9 | kmt(j-1,i+1) == -9  ...
            | kmt(j,i-1) == -9 ) & kmt(j,i) > 0
                             kmt(j,i)=-9;
            unfhmap_tmp=[unfhmap_tmp; i j 0];
         end;
        end; 
        unfhmap=[unfhmap;unfhmap_tmp];
        unfhmap_tmp=[];
        for i=ii:-1:(is+1)
         if ( kmt(j-1,i-1) == -9 | kmt(j-1,i) == -9 | kmt(j-1,i+1) == -9  ...
            | kmt(j,i+1) == -9 ) & kmt(j,i) > 0
                             kmt(j,i)=-9;
                unfhmap_tmp=[unfhmap_tmp; i j 0];
         end;
        end;
        unfhmap=[unfhmap;unfhmap_tmp];
      end;

      for j=(jj-1):-1:js
       percent=percent+dpercent;
       fprintf(' \b\b\b\b\b\b\b\b\b\b\b\b%3i percent',round(percent));
       unfhmap_tmp=[];
       for i=(ii+1):(ie-1)
        if ( kmt(j+1,i-1) == -9 | kmt(j+1,i) == -9 | kmt(j+1,i+1) == -9  ...
           | kmt(j,i-1) == -9 ) & kmt(j,i) > 0
                            kmt(j,i)=-9;
                unfhmap_tmp=[unfhmap_tmp; i j 0];
        end;
       end;
       unfhmap=[unfhmap;unfhmap_tmp];
       unfhmap_tmp=[];
       for i=ii:-1:(is+1)
        if ( kmt(j+1,i-1) == -9 | kmt(j+1,i) == -9 | kmt(j+1,i+1) == -9  ...
           | kmt(j,i+1) == -9 ) & kmt(j,i) > 0
                            kmt(j,i)=-9;
                unfhmap_tmp=[unfhmap_tmp; i j 0];
        end;
       end;
       unfhmap=[unfhmap;unfhmap_tmp];
      end;
      n=find( kmt == -9);
      kmt(n)=0*kmt(n);
      fprintf(' done\n');
    
      set(ui0,'visible','off');
      clear ui0

     end;    

  elseif choice == 7

    %--------------------------------------
    %  undo function
    %--------------------------------------

      kmt=unkmt;
      unfh='c oops \n';
      unfhmap=[];
      newplot=1;
      
  elseif choice == 1

    %--------------------------------
    % exit this m-script
    %--------------------------------
    
    fprintf('ok, exit\n');
    stop = 0;

  elseif choice == 4

    %---------------------------------------
    % enable the 3-D mode in another window
    %---------------------------------------

    ui0=uicontrol('style','text');
    set(ui0,'units','normalized');
    set(ui0,'position',[0.2 0.95 0.6 .05]);     
    set(ui0,'string',' use the second graphical window');

    % change window
    set(fig2,'visible','on');
    figure(fig2); clf

    % vertical axis in 3-D-mode is in depth
    je2=min(je+1,jmt);
    ie2=min(ie+1,imt);
    zmap=zeros((je2-js+1)*2,(ie2-is+1)*2);
    for j=js:je2
     for i=is:ie2
      zmap(((j-js+1)*2-1):((j-js+1)*2),((i-is+1)*2-1):((i-is+1)*2))= ...
        -zt(kmt(j,i)+1)*ones(2,2);
     end;
    end;

    % mesh the kmt- map

    surf([x(is:0.5:ie2) x(ie2)]-0.5,[y(js:0.5:je2) y(je2)]-0.5,zmap)
    clear ie2 je2
    view(angle1,angle2);
    title('3-D view')
    xlabel('i -->')
    ylabel('j -->')
    zlabel('[m]');
    
    ui1=uicontrol('style','slider');
    set(ui1,'units','normalized');
    set(ui1,'position',[0.2 0.01 0.6 0.04]);
    set(ui1,'callback', ...
      'angle1=get(ui1,''value'')*360-180;view(angle1,angle2);');
    set(ui1,'value',(angle1+180)/360);

    ui2=uicontrol('style','slider');
    set(ui2,'units','normalized');
    set(ui2,'position',[0.01 0.2 0.04 0.6]);
    set(ui2,'callback','angle2=get(ui2,''value'')*90;view(angle1,angle2);');
    set(ui2,'value',(angle2)/90);

    ui3=uicontrol('style','pushbutton');
    set(ui3,'units','normalized');
    set(ui3,'position',[0.05 0.01 0.05 0.05])
    set(ui3,'callback','stop=0;');
    set(ui3,'string','Quit');

    ui4=uicontrol('style','pushbutton');
    set(ui4,'units','normalized');
    set(ui4,'position',[0.9 0.01 0.05 0.05])
    set(ui4,'callback','printit=1;');
    set(ui4,'string','print');

    stop=1;
    printit=0;
    while stop == 1
     drawnow;

     if printit == 1

      printit = 0;

      ui5=uicontrol('style','pushbutton','units','normalized',...
        'position',[0.2 0.5 0.6 .05], 'callback','thisstop=0;', ... 
        'string','enter postscript filename (click here when ready)');
      ui6=uicontrol('style','edit','units','normalized',...
           'position',[0.2 0.45 0.6 .05], 'max',1, ...
           'string','   ');

      thisstop=1;
      while thisstop == 1
       drawnow;
      end;

      comstr=get(ui6,'string');
      
      set(ui5,'visible','off');
      set(ui6,'visible','off');
      clear ui5 ui6 

      if comstr == ''

       ui5=uicontrol('style','pushbutton','units','normalized',...
        'position',[0.2 0.5 0.6 .05], 'callback','thisstop=0;', ... 
        'string',' didn''t print ');

       thisstop=1;
       while thisstop == 1
        drawnow;
       end;

       set(ui5,'visible','off');
       clear ui5  
 
      else

       title(setstr(expnam)');
       orient portrait

       set(ui0,'visible','off');
       set(ui1,'visible','off');
       set(ui2,'visible','off');
       set(ui3,'visible','off');
       set(ui4,'visible','off');

       eval([' print -dpsc ' comstr ]);

       set(ui0,'visible','on');
       set(ui1,'visible','on');
       set(ui2,'visible','on');
       set(ui3,'visible','on');
       set(ui4,'visible','on');

       title('3-D View')

       ui5=uicontrol('style','pushbutton','units','normalized',...
        'position',[0.2 0.5 0.6 .05], 'callback','thisstop=0;', ... 
        'string',[' printed in file ' comstr ]);

       thisstop=1;
       while thisstop == 1
        drawnow;
       end;

       set(ui5,'visible','off');
       clear ui5  

      end;
      clear comstr

     end;
    end; 

    stop = 1;
    figure(fig1);
    clear zmap
    set(ui0,'visible','off');
    set(ui1,'visible','off');
    set(ui2,'visible','off');
    set(ui3,'visible','off');
    set(ui4,'visible','off');
    clear ui0 ui1 ui2 ui3 ui4

  elseif choice == 2
    
   %------------------------------------------
   % zoom in
   %------------------------------------------

   ui1=uicontrol('style','text');
   set(ui1,'units','normalized');
   set(ui1,'position',[0.2 0.95 0.6 .05]);
   set(ui1,'string','click to new north-western edge' );
   [xx1 yy1 button]=ginput(1);
   set(ui1,'string','click to new south-eastern edge' );
   [xx2 yy2 button]=ginput(1);
   set(ui1,'visible','off');
   clear ui1

   xx1=floor(xx1);
   yy1=floor(yy1);
   xx2=ceil(xx2);
   yy2=ceil(yy2);
   if ( xx1 < xx2 ) &  ( yy2 < yy1 )
    newplot=1;
    is=xx1;
    ie=xx2;
    js=yy2;
    je=yy1; 
    if is < 1; is =1; end;
    if js < 1; js =1; end;
    if ie > imt; ie =imt; end;
    if je > jmt; je =jmt; end;
   end;

  elseif choice == 3

   %---------------------------------------------
   % zoom out
   %---------------------------------------------

   stop = 1;
   ui1=uicontrol('visible','off');
   ui2=uicontrol('visible','off');

   set(ui1,'style','slider');
   set(ui1,'units','normalized');
   set(ui1,'position',[0.2 0.4 0.6 .05]);
   set(ui1,'callback', ...
     'set(ui2,''string'',[ num2str(get(ui1,''value'')*1000) '' %''] );');
   set(ui1,'value',dzoom/10);
   set(ui1,'visible','on');

   set(ui2,'style','pushbutton');
   set(ui2,'units','normalized');
   set(ui2,'position',[0.2 0.46 0.6 .1]);
   set(ui2,'callback','stop=0;');
   set(ui2,'string',[ num2str(get(ui1,'value')*1000) ' %'] );
   set(ui2,'visible','on');
 
   while stop == 1
    drawnow;
   end;
   set(ui1,'visible','off');
   set(ui2,'visible','off');
   tmp=dzoom;
   dzoom=floor( get(ui1,'value')*10 );
   clear ui1 ui2
   dzoom=min( ie-is, je-js)*dzoom;

   stop = 1;
   is=is-dzoom;
   js=js-dzoom;
   ie=ie+dzoom;
   je=je+dzoom;
   if is < 1; is =1; end;
   if js < 1; js =1; end;
   if ie > imt; ie =imt; end;
   if je > jmt; je =jmt; end;
   newplot=1;
   dzoom=tmp;
   clear tmp

  end; % if
end;

%--------------------------------
%  end of main loop: cleaning up
%--------------------------------

% flush the undo buffers for delta-file

fprintf(fnewh,unfh);
if ~isempty( unfhmap )
 [mm nn]=size( unfhmap );
 for mi=1:mm
  fprintf(fnewh,form_del,unfhmap(mi,1),unfhmap(mi,2),unfhmap(mi,3));
 end;
 fprintf(fnewh,'c map.m: done\n');
end;
unfhmap=[];

fprintf(fnewh,'c\nc\nc map.m: all done\nc\nc\nc\n');
fclose(fnewh);

%--------------------------------------------------------
% save new kmt.dta
%--------------------------------------------------------


fprintf('writing kmt.dta file \n');

fkmt=fopen('kmt.dta','w',machine_format);

if fortran_format1 == 80
 fwrite(fkmt,80,'long');
 fwrite(fkmt,AWI_ID,'char');
 fwrite(fkmt,80,'long');
end;

fwrite(fkmt,32+80+60,'long');
fwrite(fkmt,[stamp' iotext' expnam'],'char');
fwrite(fkmt,32+80+60,'long');

fwrite(fkmt,3*4,'long');
fwrite(fkmt,[imt jmt km],'long');
fwrite(fkmt,3*4,'long');

fwrite(fkmt,32+80+60,'long');
fwrite(fkmt,[stamp2' iotext2' expnam2'],'char');
fwrite(fkmt,32+80+60,'long');

fwrite(fkmt,imt*jmt*4,'long');
fwrite(fkmt,kmt','long');
fwrite(fkmt,imt*jmt*4,'long');

fclose(fkmt);
clear fkmt

if imt == 98 & jmt == 150 & km == 45

 fprintf('\n This kmt-array seems to be the FLAME coarse resolution \n');
 fprintf(' Atlantic configuration. So, writing also the kmt-array for \n');
 fprintf(' FLAME coarse resolution North Atlantic configuration, which \n');
 fprintf(' is a subset of this kmt-array, in file kmt_N_Atl.dta. \n');
 % fprintf(' (Without the AWI-security ID !! Change this if you want.) \n');

 fkmt=fopen('kmt_N_Atl.dta','w',machine_format);

 % uncomment this to write also the AWI_ID string
 %
 if fortran_format1 == 80
  fwrite(fkmt,80,'long');
  AWI_ID(1:16)='GRE   1 TO   1  ';
  fwrite(fkmt,AWI_ID,'char');
  fwrite(fkmt,80,'long');
 end;

 fwrite(fkmt,32+80+60,'long');
 fwrite(fkmt,[stamp' iotext' expnam'],'char');
 fwrite(fkmt,32+80+60,'long');

 fwrite(fkmt,3*4,'long');
 fwrite(fkmt,[87 89 km],'long');
 fwrite(fkmt,3*4,'long');

 fwrite(fkmt,32+80+60,'long');
 fwrite(fkmt,[stamp2' iotext2' expnam2'],'char');
 fwrite(fkmt,32+80+60,'long');

 % Get the kmt-array for North Atlantic

 N_kmt=kmt(62:150,1:87);

 % and insure the correct settings for the open boundaries

 for i=1:87
  N_kmt(2,i)=min( min(N_kmt(2,i),N_kmt(3,i)), min(N_kmt(2,i),N_kmt(4,i)));
 end;
 
 N_kmt(1,:)=N_kmt(2,:);
 N_kmt(3,:)=N_kmt(2,:);
 N_kmt(4,:)=N_kmt(3,:);

  % This is Arnes Code for obc_south:
  %  do i=1,imt
  %   kmt(i,2)=min(min(kmt(i,2),kmt(i,3)),min(kmt(i,2),kmt(i,4)))
  %   kmt(i,1)=kmt(i,2)
  %   kmt(i,3)=kmt(i,2)
  %   kmt(i,4)=kmt(i,3)
  %  enddo
 
 % solid northern, eastern and western walls 
  
 N_kmt(89,:)= zeros(1,87);
 N_kmt(:,1) = zeros(89,1);
 N_kmt(:,87)= zeros(89,1);
 

 fwrite(fkmt,87*89*4,'long');
 fwrite(fkmt,N_kmt','long');
 fwrite(fkmt,87*89*4,'long');

 clear N_kmt

 fclose(fkmt);
 
 if fortran_format1 == 80
   
  fkmt=fopen('kmt_N_Atl.dta','r',machine_format);
  fread(fkmt,1,'long');
  comstr=fread(fkmt,80,'char');
  fclose(fkmt);
  
  comstr=setstr(comstr)';
  
  fprintf(' new AWI_ID for this file : %s',comstr);
  
 end;

end;

fprintf('\nall done\n\n');

% All done, hopefully.

%-------------------------------------------------------------------------
%-------------------------------------------------------------------------

