function create_cira_xml


%=== Read southern hemisphere temperatures
%
fid = cira_open( pwd, 'sht.dat' );
[ps,SH] = cira_read( fid, 't' );
fclose(fid);


%=== Read northern hemisphere temperatures
%
fid = cira_open( pwd, 'nht.dat' );
[pn,NH] = cira_read( fid, 't' );
fclose(fid);

% pn and ps are identical!


%=== Combine the temperatures for SH and NH, and sort p
%
[TT,p,lat] = cira_combine( ps, SH, NH );
%
clear pn ps SH NH


%--- Fix bad temperature values
%
for im = 1:size(TT,1)
  for ilat = 1:size(TT,2)
    
    ind = find( TT(im,ilat,:) >= 999 );
    
    if ~isempty(ind)
      TT(im,ilat,ind) = TT(im,ilat,ind(end)+1);
    end  
  end
end


%--- Store "atmosphere" temperature files
%
A.grids     = { p, lat, 0 };
A.gridnames = { 'Pressure', 'Latitude', 'Longitude' };
A.dataname  = 'Temperature';
%
for im = 1:12

  filename = sprintf( 'cira86_month%d.t.xml', im );
  
  A.name = sprintf( 'CIRA86 temperature for month %d', im );
  A.data = squeeze( TT(im,:,:) )';
  
  xmlStore( filename, A, 'GriddedField3' );
  
end


%--- Create and store climatology
%
C          = make_climatology( p, lat, TT );
C.name     = 'CIRA86 temperature';
C.dataname = 'Temperature';
%
filename   = fullfile( fileparts(fileparts(pwd)), 'climatology', 'cira86', ...
                                                               'cira86.t.xml' );
xmlStore( filename, C, 'GriddedField4' );



%--- Convert to altitudes 
%
Z   = TT;
p0  = 350e2;
%
for im = 1:size(TT,1)
  for ilat = 1:size(TT,2)
    
    t = squeeze(TT(im,ilat,:));
    
    h2o = zeros( size(p) );
    
    ind = find( p > p0  &  t > 273.15 );
    if ~isempty(ind)
      h2o(ind) = 0.5 * e_eq_water( t(ind) ) ./ p(ind);
    end

    ind = find( p > p0  &  t <= 273.15 );
    if ~isempty(ind)
      h2o(ind) = 0.5 * e_eq_ice( t(ind) ) ./ p(ind);
    end

    Z(im,ilat,:) = pt2z( p, t, h2o, 1013e2, 0, lat(ilat) );

  end
end


%--- Store altitude files
%
for im = 1:12

  filename = sprintf( 'cira86_month%d.z.xml', im );
  
  A.data = squeeze( Z(im,:,:) )';
  
  xmlStore( filename, A, 'GriddedField3' );
  
end


%--- Create and store climatology
%
C          = make_climatology( p, lat, Z );
C.name     = 'CIRA86 geometrical altitudes';
C.dataname = 'Altitude';
%
filename   = fullfile( fileparts(fileparts(pwd)), 'climatology', 'cira86', ...
                                                               'cira86.z.xml' );
xmlStore( filename, C, 'GriddedField4' );


return
%----------------------------------------------------------------------------

% Copied from cira86.m in ami/file, with some modifications.

function fid = cira_open( folderpath, name )
  fid = fopen( fullfile(folderpath,name), 'r' );
  if fid < 0
    error(sprintf('The file %s could not be opened.',name));
  end
return



function [cp,A] = cira_read( fid, t_or_z )

  if t_or_z == 't'
    n = 20;
  elseif t_or_z == 'z'
    n = 19;
  else
    error('Unknown letter code');
  end

  % Allocate a tensor for data to read
  B = zeros(12,n,71);

  % Loop months and read
  for i = 1:12

    % Read 6 lines with text
    for j = 1:6
      fgetl(fid);
    end

    % Read the data for the month
    B(i,:,:) = fscanf( fid, '%f', [n 71] );

    % Read newline character
    fgetl(fid);

  end

  cp = squeeze( B(1,2,:) ) * 100;
  A  = B( :, n + (-16:0), : );

return


function [A,p,lat] = cira_combine( ps, SH, NH )

  % Data for +-80 are copied to be valid also for +-90

  A = zeros(12,35,71);
 
  lat = [-90 -80:5:80 90 ];
  
  A(:,1,:)     = SH(:,1,:);
  A(:,35,:)    = NH(:,17,:);
  A(:,2:18,:)  = SH;
  A(:,19:34,:) = NH(:,2:17,:);

  p   = ps( end:-1:1 );
  A   = A(:,:,end:-1:1);  

return


function C = make_climatology( p, lat, A )
  
  %- Data for doy 0 and 367 will be created to cover the complete year
  dd          = daysinmonth( 2007, 1:12 );  % Just some non leap-year
  doy         = [ 0 dd(1)/2 cumsum(dd(1:11))+dd(2:12)/2 367];
  
  C.grids     = { p, lat, 0, doy };
  C.gridnames = { 'Pressure', 'Latitude', 'Longitude', 'doy' };

  C.data     = zeros(71,35,1,14);
  
  for ilat = 1 : 35
    for d = 1 :14

      if d == 1
        C.data(:,ilat,1,d) = (A(1,ilat,:) + A(12,ilat,:) ) / 2;        
      elseif d == 14
        C.data(:,ilat,1,d) = C.data(:,ilat,1,1);
      else
        C.data(:,ilat,1,d) = A(d-1,ilat,:);
      end
    end
  end
    
return