function [ok]=  write_coastwatch( hdf_filename, year, day, sst_analysis, error_analysis, ...
                                  land_mask, ice_mask, processing_date, bad_val, max_val, min_val);

% Function to write daily output file from NOAA Hires Operational SST Analysis
% in coastwatch format.
% Two arrays, sst_analysis and error_analysis containing SST analysis and associated
% error analysis are stored, together with metadata as required by COASTWATCH HDF format
% description.
%
% Combined land/ice mask is put into the HDF file
%
% in the original land_mask: 	0 is land		(0x00)
%				2 is pacific ocean	(0x02)
%				3 is atalantic ocean	(0x03)
%				5 is indian ocean	(0x05)
%				6 is arctic		(0x06)
%				7 is antarctic		(0x07)
%
% in the original ice_mask: 	0 is water		(0x00)
%				1 is ice		(0x01)
%				2 is land		(0x02)
%
% to use cwrender, the new mask is set to new_mask=16*land_mask+4*ice_mask
% so the landmask is at bit 3 (0x08)
%				4  is ice		(0x04)
%				8  is land		(0x08)
%				32 is pacific ocean	(0x20)
%				48 is atalantic ocean	(0x30)
%				80 is indian ocean	(0x50)
%				96 is arctic		(0x60)
%				112 is antarctic 	(0x70)
%

% Work out days since 1 Jan 1970 as used by COASTWATCH.
date_origin=datenum(1970,1,1);
this_date=datenum(year,1,1)+day-1;
pass_date=this_date-date_origin;

% Rotate dataset from matlab y/x to standard x/y orientation.
A=fliplr(sst_analysis');

% Create COASTWATCH dataset.
% Create sst_analysis SDS
sd_id=hdfsd('start', hdf_filename,'DFACC_CREATE');
ds_name='sst_analysis';
ds_type='double';
ds_rank=ndims(A);
ds_dims=fliplr(size(A));

sds_id=hdfsd('create', sd_id, ds_name, ds_type, ds_rank, ds_dims);

% Set missing value to -999.
status=hdfsd('setattr', sds_id, 'missing_value', bad_val);

% Set range.
status=hdfsd('setrange', sds_id, max_val, min_val);

% Set label units format coord_system 
status=hdfsd('setdatastrs', sds_id, 'blended daily goes poes sst', 'celsius', 'decimal', 'global lat lon');

ds_start=zeros(1:ndims(A));
ds_stride=[];
ds_edges=fliplr(size(A));
stat1=hdfsd('writedata',sds_id, ds_start, ds_stride, ds_edges, A);
stat2=hdfsd('endaccess', sds_id);

% Create error_analysis SDS
% Rotate dataset from matlab y/x to standard x/y orientation.
A=fliplr(error_analysis');
ds_name='error_analysis';

sds_id=hdfsd('create', sd_id, ds_name, ds_type, ds_rank, ds_dims);

% Set missing value to -999.
status=hdfsd('setattr', sds_id, 'missing_value', bad_val);

% Set label units format coord_system 
status=hdfsd('setdatastrs', sds_id, 'sst error', 'celsius', 'decimal', 'global lat lon');

ds_start=zeros(1:ndims(A));
ds_stride=[];
ds_edges=fliplr(size(A));
stat1=hdfsd('writedata',sds_id, ds_start, ds_stride, ds_edges, A);
stat2=hdfsd('endaccess', sds_id);

% Create mask SDS
% combine land_mask and ice_mask by shifting land_maks left 4 bit
% because ice_mask takes only three least significant bits (0 for water, 1(4) for ice, 2(8) for land)

mask=16*uint8(land_mask)+4*uint8(ice_mask);

% Rotate dataset from matlab y/x to standard x/y orientation.

mask=fliplr(mask');
ds_name='graphics';
ds_type='uint8';

sds_id=hdfsd('create', sd_id, ds_name, ds_type, ds_rank, ds_dims);

% Set label units format coord_system 

status=hdfsd('setdatastrs', sds_id, 'land ice water mask', 'none', 'decimal', 'global lat lon');

ds_start=zeros(1:ndims(mask));
ds_stride=[];
ds_edges=fliplr(size(mask));
stat1=hdfsd('writedata',sds_id, ds_start, ds_stride, ds_edges, mask);
stat2=hdfsd('endaccess', sds_id);

% Write metadata to file.
sds_id=sd_id;
 
% Specify satellite datasets used.
newline=char(12);
%satellite_names=['noaa17', 10, 'noaa18', 10, 'goese', 10, 'goesw'];
satellite_names=['noaa18/metop-a goese/w'];
status=hdfsd('setattr', sds_id, 'satellite', satellite_names);

% Specify sensor used.
status=hdfsd('setattr', sds_id, 'sensor', ['AVHRR', 10, 'IMAGER']);

% Set Data Origin.
status=hdfsd('setattr', sds_id, 'origin', 'USDOC/NOAA/NESDIS COASTWATCH');
%status=hdfsd('setattr', sds_id, 'origin', 'USDOC/NOAA/NESDIS/COASTWATCH \n USDOC/NOAA/NESDIS/COASTWATCH');

% Set Data Source.
%status=hdfsd('setattr', sds_id, 'data_source', 'USDOC/NOAA/NESDIS COASTWATCH');
%status=hdfsd('setattr', sds_id, 'data_source', 'USDOC/NOAA/NESDIS/COASTWATCH \n USDOC/NOAA/NESDIS/COASTWATCH \n USDOC/NOAA/NESDIS/COASTWATCH');

% Set History.
history_text=['writecoastwatch Produced by NOAA Blended SST Analysis Processing on ' processing_date];
status=hdfsd('setattr', sds_id, 'history', history_text);

% First, set COASTWATCH HDF version number.
status=hdfsd('setattr', sds_id, 'cwhdf_version', '3.4');

% Data represent a composite
%status=hdfsd('setattr', sds_id, 'composite','true');

% Set pass date to days since 1 January 1970
%status=hdfsd('setattr', sds_id, 'pass_date', [int32(pass_date) int32(pass_date) int32(pass_date)]);
status=hdfsd('setattr', sds_id, 'pass_date', int32(pass_date));

% Data cover whole day.
%status=hdfsd('setattr', sds_id, 'start_time', [0 0 0]);
%status=hdfsd('setattr', sds_id, 'temporal_extent', [86399 86399 86399]);
status=hdfsd('setattr', sds_id, 'start_time', 0);
status=hdfsd('setattr', sds_id, 'temporal_extent', 86399);

% Both day and night data used in Analysis.
status=hdfsd('setattr', sds_id, 'pass_type', 'day');

polygon_latitude=[90 90 90 0 -90 -90 -90 0 90];
polygon_longitude=[-180 0 180 180 180 0 -180 -180 -180];
%et_affine=[0.0 -11094.626 11131.949 0.0 -20037508.343 9985163.163];
et_affine=[0.0 -11119.48742 11119.48742 0.0 -20015077.37124 10007538.68562];


% Set size. Data is global at 0.1 deg spatial resolution.
status=hdfsd('setattr', sds_id, 'rows', int32(1800));
status=hdfsd('setattr', sds_id, 'cols', int32(3600));
status=hdfsd('setattr', sds_id, 'projection_type', 'mapped');
status=hdfsd('setattr', sds_id, 'projection', 'equirectangular');
status=hdfsd('setattr', sds_id, 'gctp_sys', int32(17));
status=hdfsd('setattr', sds_id, 'gctp_zone', int32(0));
status=hdfsd('setattr', sds_id, 'gctp_parm', [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]);
status=hdfsd('setattr', sds_id, 'gctp_datum', int32(12));
%status=hdfsd('setattr', sds_id, 'polygon_latitude', polygon_latitude);
%status=hdfsd('setattr', sds_id, 'polygon_longitude', polygon_longitude);
status=hdfsd('setattr', sds_id, 'et_affine', et_affine);

ok=hdfsd('end', sd_id);
