%Concatenate all NCEP known files into a reanalysis file
%in GFDL format.
%Written by Max Kaznady, winter 2004/2005
function NCEP2GFDL(output)
%first we must download the reanalysis data
! ./get_reanalysis.csh
nc = netcdf(output, 'clobber')
if isempty(nc), error('## Bad netcdf operation.'), end
nc.description = ncchar('Reanalysis file for use with CDS from NCEP/NCAR data in GFDL format.');
nc.description = 'Reanalysis file for use with CDS from NCEP/NCAR data in GFDL format.';
nc('lat') = 73;
nc('lon') = 144;
nc('time') = 12;
nc('level') = 17;
nc{'lat'} = ncfloat('lat');
nc{'lon'} = ncfloat('lon');
nc{'time'} = ncfloat('time');
nc{'level'} = ncfloat('level');
nc{'lat'}(:) = -90:2.5:90;
nc{'lon'}(:) = 0:2.5:357.5; %360? 357.5 before...
nc{'time'}(:) = [0 744 1416 2160 2880 3624 4344 5088 5832 6552 7296 8016]; %no idea - just copied
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
nc_t_ref=netcdf('air.mon.ltm.nc');
nc{'t_ref'} = ncfloat('time', 'lat', 'lon');
if isempty(nc{'t_ref'}), error(' ## Bad vardef operation.'), end
superstition=nc_t_ref{'air'}(:,:,:);
nc{'t_ref'}(:,:,:)=superstition(:,:,:);
nc{'t_ref'}(:)=flipdim(nc{'t_ref'}(:),2);
nc{'t_ref'}.missing_value = ncfloat(nc_t_ref{'air'}.missing_value);
nc{'t_ref'}.missing_value(:) = nc_t_ref{'air'}.missing_value(:);
nc{'t_ref'}.add_offset = ncfloat(nc_t_ref{'air'}.add_offset);
nc{'t_ref'}.add_offset(:) = nc_t_ref{'air'}.add_offset(:);
nc{'t_ref'}.scale_factor = ncfloat(nc_t_ref{'air'}.scale_factor);
nc{'t_ref'}.scale_factor(:) = nc_t_ref{'air'}.scale_factor(:);
nc{'t_ref'}.units = ncchar('degK');
nc{'t_ref'}.units = 'degK';
%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{'t_ref'}(i,:,:)=nc{'t_ref'}(i,:,:)+273.16;
end
%ucomp - zonal wind component
nc_ucomp=netcdf('uwnd.mon.ltm.nc');
nc{'ucomp'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'ucomp'}), error(' ## Bad vardef operation.'), end
superstition=nc_ucomp{'uwnd'}(:,:,:,:);
nc{'ucomp'}(:,:,:,:)=superstition(:,:,:,:);
nc{'ucomp'}(:)=flipdim(nc{'ucomp'}(:),3);
nc{'ucomp'}.missing_value = ncfloat(nc_ucomp{'uwnd'}.missing_value);
nc{'ucomp'}.missing_value(:) = nc_ucomp{'uwnd'}.missing_value(:);
nc{'ucomp'}.add_offset = ncfloat(nc_ucomp{'uwnd'}.add_offset);
nc{'ucomp'}.add_offset(:) = nc_ucomp{'uwnd'}.add_offset(:);
nc{'ucomp'}.scale_factor = ncfloat(nc_ucomp{'uwnd'}.scale_factor);
nc{'ucomp'}.scale_factor(:) = nc_ucomp{'uwnd'}.scale_factor(:);
nc{'ucomp'}.units = ncchar('m/s');
nc{'ucomp'}.units = 'm/s';
%vcomp - meridional wind component
nc_vcomp=netcdf('vwnd.mon.ltm.nc');
nc{'vcomp'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'vcomp'}), error(' ## Bad vardef operation.'), end
superstition=nc_vcomp{'vwnd'}(:,:,:,:);
nc{'vcomp'}(:,:,:,:)=superstition(:,:,:,:);
nc{'vcomp'}(:)=flipdim(nc{'vcomp'}(:),3);
nc{'vcomp'}.missing_value = ncfloat(nc_vcomp{'vwnd'}.missing_value);
nc{'vcomp'}.missing_value(:) = nc_vcomp{'vwnd'}.missing_value(:);
nc{'vcomp'}.add_offset = ncfloat(nc_vcomp{'vwnd'}.add_offset);
nc{'vcomp'}.add_offset(:) = nc_vcomp{'vwnd'}.add_offset(:);
nc{'vcomp'}.scale_factor = ncfloat(nc_vcomp{'vwnd'}.scale_factor);
nc{'vcomp'}.scale_factor(:) = nc_vcomp{'vwnd'}.scale_factor(:);
nc{'vcomp'}.units = ncchar('m/s');
nc{'vcomp'}.units = 'm/s';
%temp - air temperature (pressure-based)
nc_temp=netcdf('temp.mon.ltm.nc');
nc{'temp'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'temp'}), error(' ## Bad vardef operation.'), end
superstition=nc_temp{'air'}(:,:,:,:);
nc{'temp'}(:,:,:,:)=superstition(:,:,:,:);
nc{'temp'}(:)=flipdim(nc{'temp'}(:),3);
nc{'temp'}.missing_value = ncfloat(nc_temp{'air'}.missing_value);
nc{'temp'}.missing_value(:) = nc_temp{'air'}.missing_value(:);
nc{'temp'}.add_offset = ncfloat(nc_temp{'air'}.add_offset);
nc{'temp'}.add_offset(:) = nc_temp{'air'}.add_offset(:);
nc{'temp'}.scale_factor = ncfloat(nc_temp{'air'}.scale_factor);
nc{'temp'}.scale_factor(:) = nc_temp{'air'}.scale_factor(:);
nc{'temp'}.units = ncchar('degK');
nc{'temp'}.units = 'degK';
%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{'temp'}(i,:,:)=nc{'temp'}(i,:,:)+273.16;
end
%precip - total water (surface data) - NCEP/NCAR - take GFDL's data and ncrcat it?
nc_precip=netcdf('prate.mon.ltm.nc');
%nc_precip=netcdf('TEMP7.nc');
nc{'precip'} = ncfloat('time', 'lat', 'lon');
if isempty(nc{'precip'}), error(' ## Bad vardef operation.'), end
%interpolate data
for i=1:12,
nc{'precip'}(i, :, :)=interpolate_data_xlongitude_ylatitude(nc_precip, nc, nc_precip{'prate'}(i,:,:),nc{'precip'}(i,:,:));
end
%superstition=nc_precip{'prate'}(:,:,:);
%superstition=nc_precip{'precip'}(:,:,:);
%nc{'precip'}(:,:,:)=superstition(:,:,:);
%nc{'precip'}(:)=flipdim(nc{'precip'}(:),2);
nc{'precip'}.missing_value = ncfloat(nc_precip{'prate'}.missing_value);
nc{'precip'}.missing_value(:) = nc_precip{'prate'}.missing_value(:);
%nc{'precip'}.missing_value = ncfloat(nc_precip{'precip'}.missing_value);
%nc{'precip'}.missing_value(:) = nc_precip{'precip'}.missing_value(:);
nc{'precip'}.add_offset = ncfloat(nc_precip{'prate'}.add_offset);
nc{'precip'}.add_offset(:) = nc_precip{'prate'}.add_offset(:);
nc{'precip'}.scale_factor = ncfloat(nc_precip{'prate'}.scale_factor);
nc{'precip'}.scale_factor(:) = nc_precip{'prate'}.scale_factor(:);
%nc{'precip'}.add_offset = ncfloat(nc_precip{'precip'}.add_offset);
%nc{'precip'}.add_offset(:) = nc_precip{'precip'}.add_offset(:);
%nc{'precip'}.scale_factor = ncfloat(nc_precip{'precip'}.scale_factor);
%nc{'precip'}.scale_factor(:) = nc_precip{'precip'}.scale_factor(:);
nc{'precip'}.units = ncchar('kg/m2/s');
nc{'precip'}.units = 'kg/m2/s';
%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{'precip'}(i,:,:)=nc{'precip'}(i,:,:).*86400;
%kg/m^2/s=(86400000kg/m^2)*(mm/day)
%end
%slp - sea-level pressure for m_map and slp plots
%NOTE: relies on ps, so I need that too
%nc_precip=netcdf('pr_wtr.mon.ltm.nc');
nc_slp=netcdf('slp.mon.ltm.nc');
nc{'slp'} = ncfloat('time', 'lat', 'lon');
if isempty(nc{'slp'}), error(' ## Bad vardef operation.'), end
superstition=nc_slp{'slp'}(:,:,:);
nc{'slp'}(:,:,:)=superstition(:,:,:);
nc{'slp'}(:)=flipdim(nc{'slp'}(:),2);
nc{'slp'}.missing_value = ncfloat(nc_slp{'slp'}.missing_value);
nc{'slp'}.missing_value(:) = nc_slp{'slp'}.missing_value(:);
nc{'slp'}.add_offset = ncfloat(nc_slp{'slp'}.add_offset);
nc{'slp'}.add_offset(:) = nc_slp{'slp'}.add_offset(:);
nc{'slp'}.scale_factor = ncfloat(nc_slp{'slp'}.scale_factor);
nc{'slp'}.scale_factor(:) = nc_slp{'slp'}.scale_factor(:);
nc{'slp'}.units = ncchar('hPa');
nc{'slp'}.units = 'hPa';
%ps - surface pressure, needed by slp
nc_ps=netcdf('pres.mon.ltm.nc');
nc{'ps'} = ncfloat('time', 'lat', 'lon');
if isempty(nc{'ps'}), error(' ## Bad vardef operation.'), end
superstition=nc_ps{'pres'}(:,:,:);
nc{'ps'}(:,:,:)=superstition(:,:,:);
nc{'ps'}(:)=flipdim(nc{'ps'}(:),2);
nc{'ps'}.missing_value = ncfloat(nc_ps{'pres'}.missing_value);
nc{'ps'}.missing_value(:) = nc_ps{'pres'}.missing_value(:);
nc{'ps'}.add_offset = ncfloat(nc_ps{'pres'}.add_offset);
nc{'ps'}.add_offset(:) = nc_ps{'pres'}.add_offset(:);
nc{'ps'}.scale_factor = ncfloat(nc_ps{'pres'}.scale_factor);
nc{'ps'}.scale_factor(:) = nc_ps{'pres'}.scale_factor(:);
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=netcdf('hgt.mon.ltm.nc');
nc{'hght'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'hght'}), error(' ## Bad vardef operation.'), end
superstition=nc_hght{'hgt'}(:,:,:,:);
nc{'hght'}(:,:,:,:)=superstition(:,:,:,:);
nc{'hght'}(:)=flipdim(nc{'hght'}(:),3);
for i=1:12,
for level=1:17,
ind=find(nc{'level'}(level)*100>nc{'ps'}(i, :, :));
if ~isempty(ind)
temp_om=squeeze(nc{'hght'}(i, level, :, :));
temp_om(ind)=nan;
nc{'hght'}(i, level, :, :)=temp_om;
end
end
end
nc{'hght'}.missing_value = ncfloat(nc_hght{'hgt'}.missing_value);
nc{'hght'}.missing_value(:) = nc_hght{'hgt'}.missing_value(:);
nc{'hght'}.add_offset = ncfloat(nc_hght{'hgt'}.add_offset);
nc{'hght'}.add_offset(:) = nc_hght{'hgt'}.add_offset(:);
nc{'hght'}.scale_factor = ncfloat(nc_hght{'hgt'}.scale_factor);
nc{'hght'}.scale_factor(:) = nc_hght{'hgt'}.scale_factor(:);
nc{'hght'}.units = ncchar('m');
nc{'hght'}.units = 'm';
%omega - vertical velocity
nc_omega=netcdf('omega.mon.ltm.nc');
nc{'omega'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'omega'}), error(' ## Bad vardef operation.'), end
superstition=nc_omega{'omega'}(:,:,:,:);
%insert NaNs for the levels here - no data available:
for i=13:17,
superstition(:,i,:,:)=nan(12, 1, 73, 144);
end
nc{'omega'}(:,:,:,:)=superstition(:,:,:,:);
nc{'omega'}(:)=flipdim(nc{'omega'}(:),3);
%apply masking
for i=1:12,
for level=1:17,
ind=find(nc{'level'}(level)*100>nc{'ps'}(i, :, :));
if ~isempty(ind)
temp_om=squeeze(nc{'omega'}(i, level, :, :));
temp_om(ind)=nan;
nc{'omega'}(i, level, :, :)=temp_om;
end
end
end
nc{'omega'}.missing_value = ncfloat(nc_omega{'omega'}.missing_value);
nc{'omega'}.missing_value(:) = nc_omega{'omega'}.missing_value(:);
nc{'omega'}.add_offset = ncfloat(nc_omega{'omega'}.add_offset);
nc{'omega'}.add_offset(:) = nc_omega{'omega'}.add_offset(:);
nc{'omega'}.scale_factor = ncfloat(nc_omega{'omega'}.scale_factor);
nc{'omega'}.scale_factor(:) = nc_omega{'omega'}.scale_factor(:);
nc{'omega'}.units = ncchar('pascals/s');
nc{'omega'}.units = 'pascals/s';%for each month:
%for i=1:12,
% nc_pres{'ps'}(find(nc_pres{'ps'}(i,:,:)==nc_pres{'ps'}.missing_value(:)))=nan;
%end
%rh - relative humidity
nc_rh=netcdf('rhum.mon.ltm.nc');
nc{'rh'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'rh'}), error(' ## Bad vardef operation.'), end
superstition=nc_rh{'rhum'}(:,:,:,:);
%superstition=superstition+nc_rh{'rhum'}.add_offset(:); %data offset
%insert NaNs for the levels here - no data available:
for i=9:17,
superstition(:,i,:,:)=nan(12, 1, 73, 144);
end
nc{'rh'}(:,:,:,:)=superstition(:,:,:,:);
nc{'rh'}(:)=flipdim(nc{'rh'}(:),3);
nc{'rh'}.units = ncchar('%');
nc{'rh'}.units = '%';
nc{'rh'}.missing_value = ncfloat(nc_rh{'rhum'}.missing_value);
nc{'rh'}.missing_value(:) = nc_rh{'rhum'}.missing_value(:);
nc{'rh'}.add_offset = ncfloat(nc_rh{'rhum'}.add_offset);
nc{'rh'}.add_offset(:) = nc_rh{'rhum'}.add_offset(:);
nc{'rh'}.scale_factor = ncfloat(nc_rh{'rhum'}.scale_factor);
nc{'rh'}.scale_factor(:) = nc_rh{'rhum'}.scale_factor(:);
%sphum - specific humidity
nc_sphum=netcdf('shum.mon.ltm.nc');
nc{'sphum'} = ncfloat('time', 'level', 'lat', 'lon');
if isempty(nc{'sphum'}), error(' ## Bad vardef operation.'), end
superstition=nc_sphum{'shum'}(:,:,:,:);
%insert NaNs for the levels here - no data available:
for i=9:17,
superstition(:,i,:,:)=nan(12, 1, 73, 144);
end
nc{'sphum'}(:,:,:,:)=superstition(:,:,:,:);
nc{'sphum'}(:)=flipdim(nc{'sphum'}(:),3);
nc{'sphum'}.missing_value = ncfloat(nc_sphum{'shum'}.missing_value);
nc{'sphum'}.missing_value(:) = nc_sphum{'shum'}.missing_value(:);
for i=1:12
nc{'sphum'}(i,:,:,:)=nc{'sphum'}(i,:,:,:)./1000;
end
nc{'sphum'}.units = ncchar('kg/kg');
nc{'sphum'}.units = 'kg/kg';
nc{'sphum'}.add_offset = ncfloat(nc_sphum{'shum'}.add_offset);
nc{'sphum'}.add_offset(:) = nc_sphum{'shum'}.add_offset(:);
nc{'sphum'}.scale_factor = ncfloat(nc_sphum{'shum'}.scale_factor);
nc{'sphum'}.scale_factor(:) = nc_sphum{'shum'}.scale_factor(:);
%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);
%and now we must cleanup the reanalysis files
! ./cleanup_reanalysis.csh
