function create_msis_xml
  
p = 1.013000001 * z2p_simple( 0 : 1e3 : 140e3 )'; 
% Gives 1013 hPa at ground with some margin for numerical problems

load msisdata.mat
%
zMsis = 1e3 * zMsis;
NMsis = 1e6 * NMsis;

dd          = daysinmonth( 2007, 1:12 );  % Just some non leap-year
doy         = [ 0 cumsum(dd)];
doy(end)    = 367;

T.name      = 'MSIS90 temperature';
T.dataname  = 'Temperature';
T.grids     = { p, latMsis, 0, doy };
T.gridnames = { 'Pressure', 'Latitude', 'Longitude', 'doy' };
T.data      = zeros( length(p), 19, 1, 13 );

Z           = T;
Z.name      = 'MSIS90 altitude';
Z.dataname  = 'Altitude';

k = constants('BOLTZMANN_CONST');


for ilat = 1 : 19
  for d = 1 :13

    tthis = squeeze(TMsis(:,d,ilat));
    pthis = squeeze(NMsis(:,d,ilat)) .* tthis * k;

    T.data(:,ilat,1,d) = interpp( pthis, tthis, p );   
    Z.data(:,ilat,1,d) = interpp( pthis, zMsis, p );   
  
  end
end

xmlStore( 'msis90.t.xml', T, 'GriddedField4' );
xmlStore( 'msis90.z.xml', Z, 'GriddedField4' );