CCCma2GFDL(input, output)

%Convert CCCma file into GFDL reanalysis format.
%TODO: global attributes
%written by Max Kaznady, summer 2005

function CCCma2GFDL(input, output)

nc = netcdf(output, 'clobber')
nc_orig = netcdf(input)
if isempty(nc), error('## Bad netcdf operation.'), end
nc.description = ncchar('Reanalysis file for use with CDS from CCCma data in GFDL format.');
nc.description = 'Reanalysis file for use with CDS from CCCma data in GFDL format.';
nc('lat') = length(nc_orig{'lat'}(:));
nc('lon') = length(nc_orig{'lon'}(:));
nc('time') = length(nc_orig{'time'}(:));
nc('level') = length(nc_orig{'z01'}(:));
nc{'lat'} = ncfloat('lat');
nc{'lon'} = ncfloat('lon');
nc{'time'} = ncfloat('time');
nc{'level'} = ncfloat('level');
nc{'lat'}(:) = nc_orig{'lat'}(:);
nc{'lon'}(:) = nc_orig{'lon'}(:); 
nc{'time'}(:) = nc_orig{'time'}(:);
nc{'level'}(:) = [1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10];

%t_ref - air temperature at 2m

%ucomp - zonal wind component

nc{'ucomp'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'ucomp'}), error(' ## Bad vardef operation.'), end
nc{'ucomp'}(:)=nc_orig{'U'}(:);
%nc{'ucomp'}(:)=flipdim(nc{'ucomp'}(:),2);
nc{'ucomp'}.long_name = ncchar(nc_orig{'U'}.long_name);
nc{'ucomp'}.long_name = nc_orig{'U'}.long_name(:);
nc{'ucomp'}.units = ncchar('m/s');
nc{'ucomp'}.units = 'm/s';

%vcomp - meridional wind component

nc{'vcomp'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'vcomp'}), error(' ## Bad vardef operation.'), end
nc{'vcomp'}(:)=nc_orig{'V'}(:);
%nc{'vcomp'}(:)=flipdim(nc{'vcomp'}(:),2);
nc{'vcomp'}.long_name = ncchar(nc_orig{'V'}.long_name);
nc{'vcomp'}.long_name = nc_orig{'V'}.long_name(:);
nc{'vcomp'}.units = ncchar('m/s');
nc{'vcomp'}.units = 'm/s';

%temp - air temperature (pressure-based)

nc{'temp'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'temp'}), error(' ## Bad vardef operation.'), end
nc{'temp'}(:)=nc_orig{'T'}(:);
%nc{'temp'}(:)=flipdim(nc{'temp'}(:),2);
nc{'temp'}.long_name = ncchar(nc_orig{'T'}.long_name);
nc{'temp'}.long_name = nc_orig{'T'}.long_name(:);
nc{'temp'}.units = ncfloat(nc_orig{'T'}.units);
nc{'temp'}.units(:) = nc_orig{'T'}.units(:);
nc{'temp'}.units = ncchar('degK');
nc{'temp'}.units = 'degK';

%precip - total water (surface data) - NCEP/NCAR - take GFDL's data and ncrcat it?

%slp - sea-level pressure for m_map and slp plots
%NOTE: relies on ps, so I need that too

%ps - surface pressure, needed by slp

nc{'ps'} = ncfloat('time', 'lat', 'lon');
if isempty(nc{'ps'}), error(' ## Bad vardef operation.'), end
for j=1:length(nc_orig{'time'}(:)),
    nc{'ps'}(j, :, :) = squeeze(nc_orig{'PS'}(j,:,:));
end
%nc{'ps'}(:)=flipdim(nc{'ps'}(:),1);
nc{'ps'}.long_name = ncchar(nc_orig{'PS'}.long_name);
nc{'ps'}.long_name = nc_orig{'PS'}.long_name(:);
nc{'ps'}.units = ncchar('pascals');
nc{'ps'}.units = 'pascals';
%need to change data only for lat and lon, preserve time
%don't know a shortcut for doing this, so:
for i=1:12
    nc{'ps'}(i,:,:)=nc{'ps'}(i,:,:).*100; 
end

%hght - geopotential height

nc{'hght'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'hght'}), error(' ## Bad vardef operation.'), end
nc{'hght'}(:)=nc_orig{'GZ'}(:)./9.81;
%nc{'hght'}(:)=flipdim(nc{'hght'}(:),2);
nc{'hght'}.long_name = ncchar(nc_orig{'GZ'}.long_name);
nc{'hght'}.long_name = nc_orig{'GZ'}.long_name(:);
nc{'hght'}.units = ncchar('m');
nc{'hght'}.units = 'm';

%omega - vertical velocity

nc{'omega'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'omega'}), error(' ## Bad vardef operation.'), end
nc{'omega'}(:)=nc_orig{'W'}(:).*(-9.81).*1.225; %atmospheric density is not 1
%nc{'omega'}(:)=flipdim(nc{'omega'}(:),2);
nc{'omega'}.long_name = ncchar(nc_orig{'W'}.long_name);
nc{'omega'}.long_name = nc_orig{'W'}.long_name(:);
nc{'omega'}.units = ncchar('pascals/s');
nc{'omega'}.units = 'pascals/s';

%rh - relative humidity

nc{'rh'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'rh'}), error(' ## Bad vardef operation.'), end
nc{'rh'}(:)=nc_orig{'RH'}(:).*100;
%nc{'rh'}(:)=flipdim(nc{'rh'}(:),2);
nc{'rh'}.long_name = ncchar(nc_orig{'RH'}.long_name);
nc{'rh'}.long_name = nc_orig{'RH'}.long_name(:);
nc{'rh'}.units = ncchar('%');
nc{'rh'}.units = '%';

%sphum - specific humidity

nc{'sphum'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'sphum'}), error(' ## Bad vardef operation.'), end
nc{'sphum'}(:)=nc_orig{'Q'}(:);
%nc{'sphum'}(:)=flipdim(nc{'sphum'}(:),2);
nc{'sphum'}.long_name = ncchar(nc_orig{'Q'}.long_name);
nc{'sphum'}.long_name = nc_orig{'Q'}.long_name(:);
nc{'sphum'}.units = ncchar('kg/kg');
nc{'sphum'}.units = 'kg/kg';

%land_mask - masking out data over open water

%  nc_land_mask=netcdf('land.nc');
%  nc{'land_mask'} = ncfloat('lat', 'lon');
%  if isempty(nc{'land_mask'}), error(' ## Bad vardef operation.'), end
%  superstition=nc_land_mask{'land'}(:);
%  nc{'land_mask'}(:,:)=superstition(:,:);
%  nc{'land_mask'}(:)=flipdim(nc{'land_mask'}(:),1);

%show file
nc
%save file
nc=close(nc);