Home > atmlab > arts > scenegen > asg_cloudsat_ecmwf.m

asg_cloudsat_ecmwf

PURPOSE ^

ASG_CLOUDSAT_ECMWF reads in cloudsat ecmwf data

SYNOPSIS ^

function [G]=asg_cloudsat_ecmwf(data,leg,grids_cs,aux)

DESCRIPTION ^

ASG_CLOUDSAT_ECMWF reads in cloudsat ecmwf data

  reads in temperature, humidity, and ozone on
  gformat data. The ecmwf data from the cloudsat
  files only extends up to approximately 25 km.

  OUT G        gformat data 
  
  IN  data     name of datafiles
      leg      leg of cloudsat
      grids_cs grids
  OPT aux      auxilirary data (see gfin_cloudsat_ecmwf)

 example usage:
  leg=2;
  datadir='/home/bengt/CIWSIR/Dataset_gen/CloudSAT_data/Data/';
  data=[datadir,'2006166131201_00702_CS_ECMWF-AUX_GRANULE_P_R03_E00.hdf'];
  latitude=[-90 90];
  [G]=asg_cloudsat_ecmwf(data,leg,latitude,aux)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

asg_cloudsat_ecmwf.m

SOURCE CODE ^

0001 %ASG_CLOUDSAT_ECMWF reads in cloudsat ecmwf data
0002 %
0003 %  reads in temperature, humidity, and ozone on
0004 %  gformat data. The ecmwf data from the cloudsat
0005 %  files only extends up to approximately 25 km.
0006 %
0007 %  OUT G        gformat data
0008 %
0009 %  IN  data     name of datafiles
0010 %      leg      leg of cloudsat
0011 %      grids_cs grids
0012 %  OPT aux      auxilirary data (see gfin_cloudsat_ecmwf)
0013 %
0014 % example usage:
0015 %  leg=2;
0016 %  datadir='/home/bengt/CIWSIR/Dataset_gen/CloudSAT_data/Data/';
0017 %  data=[datadir,'2006166131201_00702_CS_ECMWF-AUX_GRANULE_P_R03_E00.hdf'];
0018 %  latitude=[-90 90];
0019 %  [G]=asg_cloudsat_ecmwf(data,leg,latitude,aux)
0020 
0021 function [G]=asg_cloudsat_ecmwf(data,leg,grids_cs,aux)
0022 
0023 
0024 
0025 if leg==2
0026 
0027     [G,D]=gfin_cloudsat_ecmwf([],[],data,aux,leg,[-90 90]);
0028 
0029 else
0030     %read in both leg 1 and 3
0031     
0032     [G1,D]=gfin_cloudsat_ecmwf([],[],data,aux,1,[-90 0]);
0033     [G2,D]=gfin_cloudsat_ecmwf([],[],data,aux,3,[0 90]);
0034     
0035     %merge leg 1 and 3
0036 
0037     for i=1:length(G1)
0038         G(i)=G1(i);
0039         G(i).DATA=[G2(i).DATA G1(i).DATA];
0040         G(i).GRID2=[vec2col(G2(i).GRID2)' vec2col(G1(i).GRID2)'];
0041     end
0042 
0043 end
0044 
0045 
0046 if any(strcmp(aux,'DEM_elevation'))
0047    dem=1;
0048    ind=find(strcmp({G(:).NAME},'DEM_elevation'));
0049    indv=find(G(ind).DATA==-9999);
0050    G(ind).DATA(indv)=0;
0051 end
0052 
0053 %get the ecmwf data on a pressure grid
0054 %G(1) holds pressure data on an altitude grid
0055 %G(2:end) holds atmospheric data on an altitude grid
0056 %
0057 %take the pressure grid closest to 0 degree latitude
0058 %as the pressure grid to interpolate the data on
0059 
0060 lat_ind=min( find( abs(G(1).GRID2)==min( abs(G(1).GRID2) ) ) );
0061 p_grid_ecmwf=G(1).DATA(:,lat_ind);
0062 
0063 %remove negative values of this grid
0064 p_grid_ecmwf=p_grid_ecmwf( find( p_grid_ecmwf>0 ) );
0065 
0066 %create a lower resolution lat_grid_ecmwf
0067 %lat_grid_ecmwf=[G(1).GRID2(1) G(1).GRID2(2:500:end-1)  G(1).GRID2(end)];
0068 for i=1:length(grids_cs{2})
0069     [m1,ind1]=min(abs(G(1).GRID2-grids_cs{2}(i)));
0070     lat_grid_ecmwf(i)= G(1).GRID2(ind1);
0071 end
0072 
0073 %now interpolate the data on this p_grid_ecmwf
0074 %and lat_grid_ecmwf
0075 %interp linearly and extrap to the closest value
0076   
0077 if dem
0078  ilen=length(G)-2;
0079 else
0080  ilen=length(G)-1;
0081 end
0082 
0083 for i=1:ilen
0084 
0085     G(i+1).GRID1=p_grid_ecmwf;
0086     G(i+1).GRID2=grids_cs{2};
0087     
0088     DATA=zeros(length(p_grid_ecmwf),length(lat_grid_ecmwf));
0089     DATA_extrap=zeros(length(p_grid_ecmwf),length(lat_grid_ecmwf));
0090     
0091     for j=1:length(lat_grid_ecmwf)
0092         ind_j=min(find(lat_grid_ecmwf(j)==G(1).GRID2));
0093         ind_p=find(G(1).DATA(:,ind_j)>0);
0094         
0095         DATA_extrap(:,j)=interp1(log(G(1).DATA(ind_p,ind_j)),...
0096                           G(i+1).DATA(ind_p,ind_j),...
0097                           log(p_grid_ecmwf),'nearest','extrap');
0098         
0099         DATA(:,j)=interp1(log(G(1).DATA(ind_p,ind_j)),...
0100                           G(i+1).DATA(ind_p,ind_j),...
0101                           log(p_grid_ecmwf));
0102     end
0103 
0104     DATA(find(isnan(DATA)))=DATA_extrap(find(isnan(DATA)));
0105     G(i+1).DATA=DATA;
0106 
0107 end
0108  
0109 if dem
0110  G=G(2:5);
0111  indg=find(G(4).GRID2(1:end-1)-G(4).GRID2(2:end)~=0);
0112  G(4).DATA=interp1(G(4).GRID2(indg),G(4).DATA(indg),...
0113                    grids_cs{2},'nearest','extrap');
0114  G(4).GRID2=grids_cs{2};
0115  %plot(grids_cs{2},de,G(4).GRID2,G(4).DATA,'r')
0116 else
0117  G=G(2:4);
0118 end
0119 
0120 G(1).NAME= 'Temperature field';
0121 G(1).DATA_UNIT='K';
0122 
0123 %now the units for humidity and ozone are specific humidity
0124 %we want them in vmr
0125 
0126 Ma=28.966; %molecular mass of dry air
0127 Mw=18.016; %molecular mass of water
0128 Mo=47.998; %molecular mass of ozone
0129 
0130 h2o_vmr=Ma/Mw*G(2).DATA.* (1./(1-G(2).DATA*(1+Ma/Mw)));
0131 
0132 %h2o_unk is a for calculating ozone vmr
0133 h2o_unk=(1-G(2).DATA)./(Ma/Mw*G(2).DATA+1-G(2).DATA);
0134 
0135 G(2).DATA=h2o_vmr; 
0136 G(2).NAME='Water vapour';
0137 G(2).DATA_NAME='Volume mixing ratio';
0138 
0139 G(3).DATA=G(3).DATA.*(Ma*h2o_unk+Mw*h2o_vmr)/Mo;
0140 G(3).NAME='Ozone';
0141 G(3).DATA_NAME='Volume mixing ratio';
0142 
0143 D.GRID1_NAME='Pressure';
0144 D.GRID1_UNIT='Pascal';
0145 
0146 
0147 
0148 G0 = gf_empty( 3 );
0149 field={G(:).NAME};
0150 ind=find(strcmp(field,{'Temperature field'}));
0151 H1=G0;
0152 H1=gf_set_fields(H1,'TYPE','atm_data','NAME',G(ind).NAME,...
0153                'SOURCE',G(ind).SOURCE,'DATA_NAME',G(ind).DATA_NAME,...
0154                'DATA_UNIT',G(ind).DATA_UNIT);
0155 H1=gf_set(H1,G(ind).DATA,[{G(ind).GRID1},{G(ind).GRID2'},{}]);     
0156 H1=gf_set_fields(H1,'GRID1_NAME','pressure',...
0157                'GRID1_UNIT','Pa','GRID2_NAME','latitude',...
0158                'GRID2_UNIT','degree');
0159 
0160 H1(2)=G0;
0161 ind=find(strcmp(field,{'Water vapour'}));
0162 H1(2)=G0;
0163 H1(2)=gf_set_fields(H1(2),'TYPE','atm_data','NAME','H2O',...
0164                'SOURCE',G(ind).SOURCE,'DATA_NAME','vmr',...
0165                'DATA_UNIT','1');
0166 H1(2)=gf_set(H1(2),G(ind).DATA,[{G(ind).GRID1},{G(ind).GRID2'},{}]);     
0167 H1(2)=gf_set_fields(H1(2),'GRID1_NAME','pressure',...
0168                'GRID1_UNIT','Pa','GRID2_NAME','latitude',...
0169                'GRID2_UNIT','degree');
0170 
0171 H1(3)=G0;
0172 ind=find(strcmp(field,{'Ozone'}));
0173 H1(3)=G0;
0174 H1(3)=gf_set_fields(H1(3),'TYPE','atm_data','NAME','O3',...
0175                'SOURCE',G(ind).SOURCE,'DATA_NAME','vmr',...
0176                'DATA_UNIT','1');
0177 H1(3)=gf_set(H1(3),G(ind).DATA,[{G(ind).GRID1},{G(ind).GRID2'},{}]);     
0178 H1(3)=gf_set_fields(H1(3),'GRID1_NAME','pressure',...
0179                'GRID1_UNIT','Pa','GRID2_NAME','latitude',...
0180                'GRID2_UNIT','degree');
0181         
0182 if length(field)>3
0183  for i=4:length(field)
0184    H1(i)=G0;
0185    if strcmp(field(i),'DEM_elevation')
0186       H1(i)=gf_set_fields(H1(i),'TYPE','atm_data','NAME',...
0187                'surface altitude','SOURCE',G(i).SOURCE,...
0188                'DATA_NAME',G(i).DATA_NAME,'DATA_UNIT','m');
0189       H1(i)=gf_set(H1(i),G(i).DATA',[{} {G(i).GRID2'} {} ]) ;
0190       H1(i)=gf_set_fields(H1(i),'GRID1_NAME','latitude',...
0191                'GRID1_UNIT','degree');
0192    else
0193       H1(i)=gf_set_fields(H1(i),'TYPE','atm_data','NAME',...
0194                G(i).NAME,'SOURCE',G(i).SOURCE,...
0195                'DATA_NAME',G(i).DATA_NAME,'DATA_UNIT',G(i).DATA_UNIT);
0196       H1(i)=gf_set(H1(i),G(i).DATA,[{G(i).GRID1} {G(i).GRID2} {G(i).GRID3} ]); 
0197    end
0198  end
0199 end 
0200 
0201 G=H1;
0202 
0203 % GFIN_CLOUDSAT_ECMWF reads in cloudsat data on gformat
0204 %
0205 %          G(1) will holds pressure data.
0206 %          G(2) will holds temperature data.
0207 %          G(3) will holds humidity data.
0208 %          G(4) will holds ozone data.
0209 %          If it is desired G(1+i) will holds
0210 %          auxilirary data
0211 %
0212 % FORMAT   G=gfin_cloudsat_ecmwf( D, G, data, aux, leg,latitude)
0213 %
0214 % IN       D        G format definition structure
0215 %                   can be empty
0216 %          G        original gformat data
0217 %                   can be empty
0218 %          data     name of a cloudsat file
0219 % OPT      leg      options 1,2, or 3 see cloudsat_read
0220 %                   for definitions
0221 %                   default is 2
0222 %          latitude latitude region to be read in [lat1 lat2]
0223 %                   default is [0 5]
0224 %          aux      auxilirary data,options can for instance be
0225 %                   {'DEM_elevation','Longitude','TAI_start','Profile_time'}
0226 %                   default are none
0227 %
0228 %example usage on how to read in data from
0229 %a cloudsat 2B-GEOPROF file
0230 %
0231 %datadir='/home/bengt/CIWSIR/Dataset_gen/CloudSAT_data/Data/';
0232 %data=[datadir,'2006166131201_00702_CS_ECMWF-AUX_GRANULE_P_R03_E00.hdf'];
0233 %leg=2;
0234 %latitude=[0 5];
0235 %[G,D]=gfin_cloudsat_ecmwf([],[],data,[],leg,latitude,[]);
0236 %
0237 %2007-12-04 created by Bengt Rydberg
0238 
0239 function [G,D]=gfin_cloudsat_ecmwf( D, G, data, aux, leg, latitude)
0240 
0241 
0242 if ~isempty(G) & isempty(D)
0243    error('If G is defined D must be defined')
0244 end
0245 
0246 %- Default values
0247 %
0248 if isempty(D)
0249   D.DIM        = 3;
0250   D.GRID1_NAME = 'Altitude';
0251   D.GRID1_UNIT = 'm';
0252   D.GRID2_NAME = 'Latitude';
0253   D.GRID2_UNIT = 'deg';
0254   D.GRID3_NAME = 'Longitude';
0255   D.GRID3_UNIT = 'deg';
0256   G            = asgG;
0257   G.DIMADD     = [];
0258 end
0259 
0260 fields= {'EC_height','Latitude','Pressure',...
0261                     'Temperature','Specific_humidity','Ozone'};
0262 
0263 %leg_DEFAULT      = 2;
0264 %latitude_DEFAULT = [0 5];
0265 %name_DEFAULT     = 'ECMWF_temperature';
0266 %type_DEFAULT     = 'K';
0267 %dims_DEFAULT     = [1 2];
0268 %aux_DEFAULT      = [];
0269 %set_defaults;
0270 
0271 [aux]= optargs({aux},{[]} );
0272 
0273 
0274 
0275 if D.DIM < 3
0276   error( 'D.DIM must be >= 3 to accommodate data of type cloudsat type.' );
0277 end
0278 
0279 for id = 1 : D.DIM
0280     gname = lower( D.(sprintf('GRID%d_NAME',id)) );
0281     if strncmp( gname, 'altitude', 8 )
0282       dims(1) = id;
0283     elseif strncmp( gname, 'latitude', 8 )
0284       dims(2) = id;
0285     elseif strncmp( gname, 'longitude', 8 )
0286       dims(3) = id;
0287     end
0288     
0289 end
0290 
0291 if any( dims == 0 )
0292     error( [ 'Setting of *dims* failed. Altitude, latitude or ',...
0293              'longitude was nof found among the grid names.' ] );
0294 end
0295 
0296 if ~isempty(aux)
0297    for i=1:length(aux)
0298        if isempty(find(strcmp(aux(i),fields)))
0299           fields{end+1}=aux{i};
0300        end
0301    end
0302 end
0303 
0304 %read data
0305 P = cloudsat_read(data,fields,'leg',leg,'lat',latitude);
0306 
0307 if ischar(P)
0308    error('error in reading the file')
0309 end
0310 
0311 for i=1:length(fields)
0312     if strcmp(fields{i},'EC_height')
0313        P.(fields{i})=flipud( vec2col( P.(fields{i}) ) );
0314     elseif strcmp(fields{i},'Latitude')
0315        P.(fields{i})=vec2col( P.(fields{i}) )';
0316     else
0317        P.(fields{i}) = flipud(P.(fields{i})');
0318     end
0319 end
0320 
0321 
0322 kvec=1:(length(fields)-2);
0323 
0324 H=G;
0325 for ik=kvec
0326     ind=2+ik;
0327     H(ik).NAME=fields{ind};
0328     H(ik).DATA_UNIT='?';
0329     H(ik).DATA_NAME=fields{ind};
0330     y=P.(fields{ind});
0331     
0332     %find out dimensionality of loaded data
0333     S=[size(y)==size(P.Latitude)];
0334     if S(1)==1 & S(2)==1 
0335        %data is a vector
0336        for id = 1 : length(dims)
0337            gname = sprintf( 'GRID%d', dims(id) );
0338            if id == 1
0339           [H(ik).(gname)] = [];
0340            elseif id == 2
0341           [H(ik).(gname)] = vec2col(P.Latitude)';
0342            elseif id == 3
0343               [H(ik).(gname)] = [];
0344            end
0345        end
0346        H(ik).DIMS=dims(2);
0347        H(ik).DATA=y; 
0348     elseif  S(1)==1 & S(2)==0
0349        %data is a scalar
0350        H(ik).GRID1=[];
0351        H(ik).GRID2=[];
0352        H(ik).GRID3=[]; 
0353        H(ik).DIMS=[]; 
0354        H(ik).DATA=y;   
0355     elseif  S(1)==0 & S(2)==1
0356        %data is a matrix
0357         for id = 1 : length(dims)
0358            gname = sprintf( 'GRID%d', dims(id) );
0359            if id == 1
0360               [H(ik).(gname)] = P.EC_height;
0361            elseif id == 2
0362           [H(ik).(gname)] = vec2col(P.Latitude)';
0363            elseif id == 3
0364               [H(ik).(gname)] = [];
0365            end
0366        end
0367        if dims(1)~=1 | dims(2)~=2 | dims(3)~=3
0368       %data is y(z,lat) change order if necessary
0369       grid_ind=[dims(1) dims(2)];
0370       if grid_ind(1)>grid_ind(2)
0371              y=y';
0372           end
0373        end
0374        H(ik).DATA=y;
0375        H(ik).DIMS=sort([dims(1) dims(2)]);
0376    else
0377        error('unknown dimensionality of data')
0378    end
0379 
0380    H(ik).SOURCE=data;
0381 end
0382  
0383 G=H;

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005