function gcontour_gmt(z, v, x, y, psfile, sc_mess, xlabel, ylabel, ...
         xtickint, ytickint, xannotint, yannotint, title, text1, ...
         x_inches, y_inches, tension, coast, landscape, ...
         vertical_scale, rectangular_grid, append, more_code, ident)

[nrows ncols ] = size(z);

% Do some error checking.

lenx = length(x);
if lenx ~= ncols
	err = 1;
	stm = [ 'The x axis has length ' num2str(lenx) ...
              ' but there are ' num2str(ncols) ' columns'];
	disp(stm)
	return
end

leny = length(y);
if leny ~= nrows
	err = 1;
	stm = [ 'The y axis has length ' num2str(leny) ...
              ' but there are ' num2str(nrows) ' rows'];
	disp(stm)
	return
end

% Write out the (x, y, z) array of data to be contoured.  Note
% that missing values are assumed to be marked with a NaN and are
% not written out because psmask does not recognise them as missing
% values.

aa = [];
for j = 1:ncols
   for i = 1:nrows
      if ~isnan(z(i, j))
         aa = [ aa; x(j), y(i), z(i, j) ];
      end
   end
end
data_file = [ 'gmt' num2str(ident) '.dat' ];
str = [ 'delete ' data_file ];
eval(str);
str = [ 'save ' data_file ' aa /ascii' ];
eval(str);

% Write out the array of contour intervals and whether the contours
% should be annotated or not.

lenv = length(v);
con_file = [ 'gmt' num2str(ident) '.con' ];
%delete gmt.con
str = [ 'delete ' con_file ];
eval(str);
for i = 1:lenv
	fprintf(con_file, '%14.8e  a\n', v(i));
end

% Write out the array of colours between each contour interval.

cpt_file = [ 'gmt' num2str(ident) '.cpt' ];
%delete gmt.cpt
str = [ 'delete ' cpt_file ];
eval(str);
colour_int = 250/(lenv - 1);
for i = 1:lenv-1
	ss = round(colour_int*i);
	fprintf(cpt_file, '%14.8e %3.0f ', v(i), ss);
	fprintf(cpt_file, '%3.0f %3.0f ', ss, ss);
	fprintf(cpt_file, '%14.8e %3.0f ', v(i+1), ss);
	fprintf(cpt_file, '%3.0f %3.0f\n', ss, ss);
end

% Calculate variables required for seeting gmt options.  Note that
% this assumes that the input data is on a rectangular grid.  The
% information is used by psmask to make the mask if there are missing
% values and xyz2grd if there are no missing values.


x_inc = x(2) - x(1);
y_inc = y(2) - y(1);
west = x(1);
east = x(lenx);
south = y(1);
north = y(leny);

if append
   a_string = ' >>';
else
   a_string = ' >';
end

if more_code
   k_string = ' -K';
else
   k_string = '';
end

if landscape
   p_string = '';
else
   p_string = ' -P';
end

if vertical_scale
   h_string = '';
   scx_inches = x_inches + 0.25;
   scy_inches = y_inches/2;
%   sclen_inches = y_inches/2;
%   scwid_inches = y_inches/5;
   scwid_inches = min([ (11.0 - x_inches)/4 2 ]);
   sclen_inches = min([ 5*scwid_inches y_inches ]);
else
   h_string = 'h';
   scx_inches = y_inches/2;
   scy_inches = 10.25;
%   sclen_inches = y_inches/2;
%   scwid_inches = y_inches/5;
   scwid_inches = min([ (9.75 - y_inches)/2 x_inches/5 ]);
   sclen_inches = min([ 5*scwid_inches x_inches ]);
end

dx_inches = 0.05;
dy_inches = 0.05;
xscale = x_inches/(east - west);
yscale = y_inches/(north - south);
dx = dx_inches/xscale;
dy = dy_inches/yscale;
i1_string = [ ' -I' num2str(x_inc) '/' num2str(y_inc) ];
i2_string = [ ' -I' num2str(dx) '/' num2str(dy) ];
j_string = [ ' -Jx' num2str(xscale) '/' num2str(yscale) ];
r_string = [ ' -R' num2str(west) '/' num2str(east) '/' ...
              num2str(south) '/' num2str(north) ];
xpos_text = west;
ypos_text = south - 0.5/yscale;

b_string = [ ' -Bf' num2str(xtickint) 'a' num2str(xannotint) ':"' ...
         xlabel '":/f' num2str(ytickint) 'a' num2str(yannotint) ':"' ...
         ylabel '":/:."' title '":WeSn' ];

index = find(z == NaN);
missing_data = length(index);

% Write the shell script commands to gmt.sh

if append == 0
   delete gmt.sh
   str = '#!/bin/sh\n';
   fprintf('gmt.sh', str);
end
str = [ 'rm ' psfile '\n' ];
fprintf('gmt.sh', str);

if missing_data
   if rectangular_grid
      str = [ 'psmask ' data_file i1_string j_string ' -Dmisspath' ...
              r_string ' -K -O' p_string ' -V > /dev/null\n' ];
      fprintf('gmt.sh', str);


      str = [ 'grdmask misspath_0.xy' i2_string ' -Mmask.grd' ...
               r_string ' -FNaN/NaN/1 -V\n' ];
      fprintf('gmt.sh', str);

      str = [ 'i=1\n' ];
      fprintf('gmt.sh', str);
      str = [ 'while test -f misspath_$i.xy\n' ];
      fprintf('gmt.sh', str);
      str = [ 'do grdmask misspath_$i.xy' i2_string ...
              ' -Mmasktemp.grd' r_string ' -F1/NaN/NaN -V\n' ];
      fprintf('gmt.sh', str);
      str = [ '   grdmath mask.grd x masktemp.grd = mask.grd\n' ];
      fprintf('gmt.sh', str);
      str = [ '   i=`expr $i + 1`\n' ];
      fprintf('gmt.sh', str);
      str = [ 'done\n' ];
      fprintf('gmt.sh', str);
  
      str = ['surface ' data_file ' -Fgmt_out1.grd' i2_string ...
             r_string ' -T' num2str(tension) ' -V\n' ];
      fprintf('gmt.sh', str);
      str = [ 'grdmath gmt_out1.grd x mask.grd = gmt_out.grd\n'];
      fprintf('gmt.sh', str);
   else
      str = ['surface ' data_file ' -Fgmt_out.grd' i2_string ...
             r_string ' -T' num2str(tension) ' -V\n' ];
      fprintf('gmt.sh', str);
   end
else
   if rectangular_grid
      str = [ 'xyz2grd ' data_file ' -Ggmt_in.grd' i1_string ...
              r_string ' -V\n' ];
      fprintf('gmt.sh', str);
      str = [ 'grdsample gmt_in.grd -Ggmt_out.grd' i2_string ...
               r_string ' -V\n' ];
      fprintf('gmt.sh', str);
      str = [ 'rm gmt_in.grd\n' ];
      fprintf('gmt.sh', str);
   else
      str = ['surface ' data_file ' -Fgmt_out.grd' i2_string ...
             r_string  ' -T' num2str(tension) ' -V\n' ];
      fprintf('gmt.sh', str);
   end
end
str = [ 'grdimage gmt_out.grd -C' cpt_file j_string ' -K' p_string ...
        ' -V' a_string psfile '\n' ];
fprintf('gmt.sh', str);
str = [ 'grdcontour gmt_out.grd -C' con_file j_string ' -K -O' ...
        p_string ' -V >> ' psfile '\n' ];
fprintf('gmt.sh', str);
if coast
   str = [ 'pscoast' j_string ' -K -O' ...
        p_string r_string ' -V -W >> ' psfile '\n' ];
   fprintf('gmt.sh', str);
end
if scwid_inches > 0.5
   str = [ 'psscale -C' cpt_file ' -D' num2str(scx_inches) '/' ...
         num2str(scy_inches) '/' num2str(sclen_inches) '/' ...
         num2str(scwid_inches) h_string ' -K -L"' sc_mess '" -O'...
         p_string  ' -V >> ' psfile '\n' ];
   fprintf('gmt.sh', str);
end
str = [ 'pstext' j_string r_string ' -K -O' p_string ' -V <<END>> ' ...
        psfile '\n' num2str(xpos_text) ' ' num2str(ypos_text) ...
        ' 20 0.0 1 1 ' text1 '\n' ...
        'END\n' ];
fprintf('gmt.sh', str);

str = [ 'psbasemap' b_string j_string k_string r_string ' -O' ...
         p_string ' -U -V >> ' psfile '\n' ];
fprintf('gmt.sh', str);

str = [ 'rm ' con_file ' ' data_file ' ' cpt_file ' mask.grd ' ...
        'masktemp.grd gmt_out.grd gmt_out1.grd misspath_*.xy\n' ];
fprintf('gmt.sh', str);


!chmod u+x gmt.sh