Home > atmlab > demos > asg_demo.m

asg_demo

PURPOSE ^

ASG_DEMO A simple demonstration how to use ASG

SYNOPSIS ^

function [G,FS]=asg_demo(file1,file2, varargin )

DESCRIPTION ^

 ASG_DEMO   A simple demonstration how to use ASG

    Gives an example on basic usage of ASG

    The example generates two atmospheric states
    including clouds.
    In the example radar data from Cloudsat
    is used to create cloud structure.
    ECMWF and Fascod climatological
    data is used to create other atmospheric 
    parameters.
    A lot of default settings are made in this
    demo, which may be re-considered depending
    on what output data should be used for

 FORMAT [G]=asg_demo([ dim,instrument ])

 OUT     G    structure with atmospheric variables on gformat data
         FS   structure specifying some settings

 IN      file1 cloudsat level 2-b geoprof file 
         file2 cloudsat ecmwf-aux file 
 OPT     sensor  
              the choice of instrument determines
              which gas species to include in G
              can be odin (default) or amsu 
         dim  
              desired dimensionality of data
              can be 1 (pressure) or 3(pressure,latitude,longitude)
              default=3 
         leg  
              the leg of Cloudsat-orbit data will be taken from
              1,2(default),or 3 
         use_lats 
              the center latitudes 
              
         lat_res 
              latitude resolution to use
              default=0.25 degrees 
         dlat  
              defines how large latitude sections that will be considered
              default=4 degrees

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

asg_demo.m

SOURCE CODE ^

0001 % ASG_DEMO   A simple demonstration how to use ASG
0002 %
0003 %    Gives an example on basic usage of ASG
0004 %
0005 %    The example generates two atmospheric states
0006 %    including clouds.
0007 %    In the example radar data from Cloudsat
0008 %    is used to create cloud structure.
0009 %    ECMWF and Fascod climatological
0010 %    data is used to create other atmospheric
0011 %    parameters.
0012 %    A lot of default settings are made in this
0013 %    demo, which may be re-considered depending
0014 %    on what output data should be used for
0015 %
0016 % FORMAT [G]=asg_demo([ dim,instrument ])
0017 %
0018 % OUT     G    structure with atmospheric variables on gformat data
0019 %         FS   structure specifying some settings
0020 %
0021 % IN      file1 cloudsat level 2-b geoprof file
0022 %         file2 cloudsat ecmwf-aux file
0023 % OPT     sensor
0024 %              the choice of instrument determines
0025 %              which gas species to include in G
0026 %              can be odin (default) or amsu
0027 %         dim
0028 %              desired dimensionality of data
0029 %              can be 1 (pressure) or 3(pressure,latitude,longitude)
0030 %              default=3
0031 %         leg
0032 %              the leg of Cloudsat-orbit data will be taken from
0033 %              1,2(default),or 3
0034 %         use_lats
0035 %              the center latitudes
0036 %
0037 %         lat_res
0038 %              latitude resolution to use
0039 %              default=0.25 degrees
0040 %         dlat
0041 %              defines how large latitude sections that will be considered
0042 %              default=4 degrees
0043 %
0044 
0045 % Created by Bengt Rydberg
0046 
0047 function [G,FS]=asg_demo(file1,file2, varargin )
0048 %
0049 
0050 [sensor,dim,leg,use_lats,lat_res,dlat] = ...
0051       optargs( varargin, {'odin',3,2,13,0.25,4} );
0052 
0053 if sensor~='odin' & sensor~='amsu'
0054   error('sensor can only be odin or amsu')
0055 end
0056 
0057 if dim~=1 & dim~=3
0058  error('dim must be 1 or 3')
0059 end
0060 
0061 if ~isscalar(leg)
0062  error('leg must be a scalar')
0063 end
0064 
0065 if leg~=1 & leg~=2 & leg~=3
0066  error('leg must be 1,2, or 3')
0067 end
0068 
0069 if ~isvector(use_lats)
0070  error('use_lats must be a vector')
0071 end
0072 
0073 if leg==1
0074    if any(use_lats<-90) | any(use_lats>0)
0075       error('if leg=1 use_lats must be within (-90 - 0)')
0076    end
0077 elseif leg==2
0078    if any(use_lats<-90) | any(use_lats>90)
0079       error('if leg=2 use_lats must be within (-90 - 90)')
0080    end
0081 else
0082   if any(use_lats<0) | any(use_lats>90)
0083       error('if leg=3 use_lats must be within (0 - 90)')
0084   end
0085 end
0086 if ~isscalar(dlat)
0087  error('dlat must be a scalar')
0088 end
0089  
0090 
0091 %SET BASIC PROPERTIES
0092 
0093 %SET CLOUD VARIABLES
0094 %-
0095 
0096 %DEFINE FILE FROM WHERE TO READ DATA
0097 %toppath = fileparts( fileparts( which( 'atmlab_init' ) ) );
0098 %datapath=fullfile(toppath,'demos','data');
0099 %cloudsatfile='2006224035643_01541_CS_2B-GEOPROF_GRANULE_P_R04_E01.hdf';
0100 %FS.DATA=fullfile(datapath,cloudsatfile);
0101 FS.DATA=file1;
0102 
0103 %DEFINE FROM WHAT LEG OF CLOUDSAT-ORBIT DATA WILL BE TAKEN FROM
0104 %Leg1=(-90,0);Leg2=(-90,90);Leg3=[0,90]
0105 FS.LEG=leg;
0106 
0107 %DEFINE FROM WHICH LATITUDES TO READ DATA
0108 FS.USE_LATS=use_lats;
0109 
0110 %DEFINE WHICH LATITUDE RESOLUTION TO USE
0111 %data will be binned to this resolution
0112 FS.LAT_RES=lat_res; 
0113 
0114 %define pressure limits from where to consider Cloudsat data
0115 FS.P_LIMS=[900e2 80e2];
0116 
0117 FS.DIM=dim;
0118 
0119 if dim==3
0120      
0121    %define the method to expand cloudsat data from 2-d to 3-d
0122    FS.DIMADD.METHOD='iaaft';
0123    
0124    %define how large latitude sections that will be considered
0125    FS.DLAT=dlat;
0126 
0127    %define longitude grid to use for cloud properties
0128    FS.LON_GRID=[-1:FS.LAT_RES:1];
0129 
0130 end
0131 
0132 %define properties for converting radar backscatter IWC/LWC fields
0133 
0134 f_grid=94e9;
0135 lambda=constants('SPEED_OF_LIGHT')/f_grid;
0136 
0137 %first ice particles
0138 
0139 P(1).PSD='MH97';
0140 P(1).shape='sphere';
0141 P(1).method='gauss-laguerre';
0142 P(1).F_GRID=f_grid;
0143 %settings for gauss-laguerre
0144 alpha=0;Ni=10;xnorm=50e-6;
0145 [x_i,w_i]=gauss_laguerre(alpha,Ni,xnorm);
0146 %x_i will be the diameter of the particles considered
0147 P(1).x=x_i;
0148 P(1).w=w_i;
0149 P(1).x_norm=xnorm;
0150 P(1).T_grid=173:10:273;
0151 %calculate radar back-scattering per particle
0152 for i=1:length(P(1).T_grid)
0153      m=sqrt(eps_ice_liebe93( f_grid, P(1).T_grid(i) ));
0154      [y(:,i)]=mie_back(m,lambda,x_i);
0155 end
0156 P(1).radar_back=y;
0157 
0158 % water particles
0159 
0160 P(2).PSD='Water';
0161 P(2).shape='sphere';
0162 P(2).method='gauss-laguerre';
0163 P(2).F_GRID=f_grid;
0164 %settings for gauss-laguerre
0165 Nj=3;xnormw=10e-6;
0166 [x_k,w_k]=gauss_laguerre(alpha,Nj,xnormw);
0167 %x_k will be the diameter of the particles considered
0168 P(2).x=x_k;
0169 P(2).w=w_k;
0170 P(2).x_norm=xnormw;
0171  
0172 P(2).T_grid=[273:10:313];
0173 clear y
0174 %calculate radar back_scattering per particle
0175 for i=1:length(P(2).T_grid)
0176      m=sqrt(eps_water_liebe93( f_grid, P(2).T_grid(i) ));
0177      [y(:,i)]=mie_back(m,lambda,x_k);
0178 end
0179 P(2).radar_back=y;
0180 
0181 FS.PROPS=P;
0182 %END OF CLOUD VARIABLES SETTINGS
0183 
0184 
0185 
0186 %SET ATMOSPHERIC PROPERTIES (no cloud)
0187 %-
0188 
0189 %grids_cs will be the grids of the atmospheric variables (not for clouds)
0190 %however the grids will later be set tighter around the cloudbox
0191 %depending on the FS structure
0192 p_grid=z2p_simple([-1e3:1e3:25e3 30e3:5e3:70e3]);
0193 CS.GRIDS{1}=p_grid; 
0194 CS.GRIDS{2}=sort([-90 -50 -30:10:30 50 90 FS.USE_LATS]);
0195 
0196 if dim==3
0197   lon_grid=[FS.LON_GRID(1) 0 FS.LON_GRID(end)];
0198   CS.GRIDS{3}=unique(sort([-180 -150 -10  -5 lon_grid 5 10 150 180]));
0199 end
0200 
0201 
0202 %define from where to read data
0203 %the baseline here is that gas fields are read from fascod climatology (ftc)
0204 %however for H2O and O3 ECMWF data are here merged with Fascod data(ecmwfftc)
0205 %N2 is here fixed by giving it a numerical value
0206 %to show some different options
0207 
0208 fascod     = fullfile( atmlab('ARTS_XMLDATA_PATH'),...
0209                            'atmosphere', 'fascod' );
0210 clim='tropical';
0211 %atmdatadir=fullfile(fascod,clim);
0212 
0213 %ecmwffile='2006224035643_01541_CS_ECMWF-AUX_GRANULE_P_R04_E01.hdf';
0214 %CS.FILE=fullfile(datapath,ecmwffile);
0215 CS.FILE=file2;
0216 
0217 %define from what leg of cloudsat-orbit data will be taken from
0218 CS.LEG=FS.LEG;
0219 
0220 CS.INSTRUMENT=sensor;
0221 
0222 %mandatory fields to include
0223 
0224 fields=[{'temperature'},{'Altitude'},{'Surface altitude'}];
0225 atmfields(1).type='temperature';
0226 atmfields(1).source='ecmwfftc';
0227 atmfields(2).type='altitude';
0228 atmfields(2).source='';
0229 atmfields(3).type='surface altitude';
0230 atmfields(3).source='ecmwf';
0231 
0232 %define which gas types to include (depends on choice of instrument)
0233 
0234 if strcmp('odin',sensor)
0235 
0236    %gases to include for Odin
0237    gases=[{'H2O'},{'O3'},{'ClO'},{'N2O'},{'O2'},{'HNO3'},{'N2'}];
0238    %ABSORPTION MODELS
0239    abs_models{1}={'H2O','H2O-ForeignContStandardType',...
0240                   'H2O-SelfContStandardType'};
0241    abs_models{2}={'O3'};
0242    abs_models{3}={'ClO'};
0243    abs_models{4}={'N2O'};
0244    abs_models{5}={'O2','O2-PWR98'};
0245    abs_models{6}={'HNO3'};
0246    abs_models{7}={'N2-SelfContMPM93'};
0247 end
0248 
0249 if strcmp('smiles',sensor)
0250 
0251    %gases to include for Odin
0252    gases=[{'H2O'},{'O3'},{'ClO'},{'N2O'},{'O2'},{'HNO3'},{'N2'},{'HCl'}];
0253    %ABSORPTION MODELS
0254    abs_models{1}={'H2O','H2O-ForeignContStandardType',...
0255                   'H2O-SelfContStandardType'};
0256    abs_models{2}={'O3'};
0257    abs_models{3}={'ClO'};
0258    abs_models{4}={'N2O'};
0259    abs_models{5}={'O2','O2-PWR98'};
0260    abs_models{6}={'HNO3'};
0261    abs_models{7}={'N2-SelfContMPM93'};
0262    abs_models{8}={'HCl'};
0263 end
0264  
0265  
0266 if strcmp('amsu',sensor)  
0267 
0268    %gases to include for amsu
0269    gases=[{'H2O'},{'O3'},{'O2'},{'N2'}];
0270 
0271    %absorption models amsu
0272    abs_models{1}={'H2O-PWR98'};
0273    abs_models{2}={'O3'};
0274    abs_models{3}={'O2-PWR93'};
0275    abs_models{4}={'N2-SelfContStandardType'};
0276 
0277 end
0278 
0279 
0280 len=length(atmfields);
0281 
0282 for i=1:length(gases) 
0283 
0284   atmfields(len+i).type='gas_species';
0285   atmfields(len+i).abs_models=abs_models{i};
0286   atmfields(len+i).gas=gases{i};
0287 
0288   if strcmp(gases{i},'H2O') | strcmp(gases{i},'O3')
0289 
0290      atmfields(len+i).source='ecmwfftc';
0291   
0292   elseif strcmp(gases{i},'N2')
0293 
0294      atmfields(len+i).source=0.7814;
0295   
0296   else
0297 
0298      atmfields(len+i).source='ftc'; 
0299    
0300   end
0301 
0302 end
0303 
0304 
0305 %Baseline for parameters describing how perturbations will be performed
0306 
0307 RND.FORMAT    = 'param'; 
0308 RND.SEPERABLE = 1;
0309 RND.CCO       = 0.01;               % Cut-off for correlation values
0310 RND.TYPE      = 'rel';              % Relative disturbances as default
0311 RND.DATALIMS  = [0];                % Do not allow negative values
0312 %
0313 RND.SI        = 0.4;                % 40% std. dev. as default
0314 %
0315 RND.CFUN1     = 'exp';              % Exp. correlation function for p-dim.
0316 RND.CL1       = [0.15 0.3 0.3]';    % Corr. length varies with altitude
0317 RND.CL1_GRID1 = [1100e2 10e2 1e-3];    
0318 %
0319 RND.CFUN2     = 'lin';              % Linear correlation function for lat-dim.
0320 RND.CL2       = 0.5;                % Corr. length 0.5 deg everywhere
0321 %
0322 RND.CFUN3     = 'lin';              % Linear correlation function for lon-dim.
0323 RND.CL3       = 0.5;                % Corr. length 0.5 deg everywhere
0324 
0325 %describe perturbations for each variable
0326 for i=1:length(atmfields)
0327   
0328   if strcmp(atmfields(i).type,'temperature')
0329 
0330      atmfields(i).RND=RND;
0331      atmfields(i).RND.TYPE='abs';
0332      atmfields(i).RND.SI=1; %1 K noise added
0333       
0334 
0335   elseif strcmp(atmfields(i).gas,'H2O')
0336 
0337      atmfields(i).RND=RND;
0338      atmfields(i).RND.SI=0.1; 
0339  
0340   elseif strcmp(atmfields(i).gas,'O3')
0341      
0342      atmfields(i).RND=RND;
0343      atmfields(i).RND.SI=0.2;
0344 
0345   elseif strcmp(atmfields(i).gas,'ClO') | strcmp(atmfields(i).gas,'N2O') | ...
0346          strcmp(atmfields(i).gas,'HNO3') | strcmp(atmfields(i).gas,'HCl') 
0347 
0348      atmfields(i).RND=RND;
0349   
0350   else
0351      
0352      atmfields(i).RND=[];
0353  
0354   end
0355 
0356 end
0357  
0358 %CS.ATMFIELDS=atmfields;
0359 
0360 %end of clear sky variables settings
0361 %--------------------------------------------------
0362 
0363 %Read basic data
0364 
0365 
0366 clear G H G1
0367 
0368 G0 = gf_empty( 3 );
0369 G0.PROPS=[];
0370 G0.DIMADD.METHOD='expand';
0371 G0.RNDMZ=[];
0372 G0.SPCFC=[];
0373 G0.SURFACE=0;
0374 
0375 
0376 for i=1:length(atmfields)
0377      
0378     G(i)=G0;
0379 
0380     if strcmp(atmfields(i).type,'temperature')
0381 
0382        if strcmp(atmfields(i).source,'ftc') | ...
0383           strcmp(atmfields(i).source,'ecmwfftc') 
0384 
0385           G1 = gf_artsxml(fullfile( fascod, [clim,'.t.xml'] ), ...
0386                              'Temperature field','t_field' );
0387           G(i)=gf_set_fields(G(i),'TYPE',G1.TYPE,'NAME',G1.NAME,...
0388                'SOURCE',G1.SOURCE,'DATA_NAME',G1.DATA_NAME,...
0389                'DATA_UNIT',G1.DATA_UNIT);
0390           G(i)=gf_set(G(i),G1.DATA,[{G1.GRID1},{},{}]);     
0391           G(i)=gf_set_fields(G(i),'GRID1_NAME',G1.GRID1_NAME,...
0392                 'GRID1_UNIT',G1.GRID1_UNIT);
0393           G(i).RNDMZ=atmfields(i).RND;
0394        end
0395          
0396        if strcmp(atmfields(i).source,'ecmwfftc') 
0397        
0398           H=asg_cloudsat_ecmwf(CS.FILE,CS.LEG,CS.GRIDS,{'DEM_elevation'});
0399 
0400        end
0401 
0402      end
0403   
0404      if strcmp(atmfields(i).type,'altitude')
0405      
0406         %altitude field is here just initialised
0407         G(i)=gf_set_fields(G(i),'TYPE','atm_data','NAME','altitude field',...
0408              'DATA_NAME','altitude','DATA_UNIT','m');
0409         G(i).SPCFC.P0        = 1013e2;
0410         grids=[{ G(i).SPCFC.P0} {} {}]; 
0411         G(i).SPCFC.Z0=gf_set( G(i),0,grids);
0412         G(i).SPCFC.Z0=gf_set_fields(G(i).SPCFC.Z0,'NAME','Altitude at P0',...
0413                       'DATA_UNIT','m','GRID1_NAME','pressure',...
0414                       'GRID1_UNIT','Pa');      
0415      end 
0416        
0417      if strcmp(atmfields(i).type,'surface altitude')
0418           
0419         if strcmp(atmfields(i).source,'ecmwf')
0420   
0421            ind=find(strcmp({H.NAME},'surface altitude'));  
0422            G(i)=gf_set_fields(G(i),'TYPE',H(ind).TYPE,...
0423                'NAME',H(ind).NAME,'SOURCE',H(ind).SOURCE,...
0424                'DATA_NAME',H(ind).DATA_NAME,'DATA_UNIT',...
0425                 H(ind).DATA_UNIT);
0426            G(i)=gf_set(G(i),H(ind).DATA,[{H(ind).GRID1} {} {} ]); 
0427            G(i)=gf_set_fields(G(i),'GRID1_NAME',H(ind).GRID1_NAME,...
0428                'GRID1_UNIT',H(ind).GRID1_UNIT);
0429 
0430 
0431         end
0432        
0433      end
0434             
0435      
0436      if strcmp(atmfields(i).type,'gas_species')
0437 
0438         if strcmp(atmfields(i).source,'ftc') |  ...
0439            strcmp(atmfields(i).source,'ecmwfftc')               
0440            sourcefile= fullfile( fascod, [clim,'.',...
0441                                  atmfields(i).gas,'.xml'] );
0442           
0443            G1= gf_artsxml( sourcefile,atmfields(i).gas,'vmr_field');
0444            G(i)=gf_set_fields(G(i),'TYPE','atmdata','NAME',G1.NAME,...
0445                'DATA_NAME',G1.DATA_NAME,'DATA_UNIT','1');
0446            G(i)=gf_set(G(i),G1.DATA,[{G1.GRID1},{},{}]);     
0447            G(i)=gf_set_fields(G(i),'GRID1_NAME',G1.GRID1_NAME,...
0448                 'GRID1_UNIT',G1.GRID1_UNIT);
0449            G(i).PROPS=atmfields(i).abs_models;
0450            G(i).RNDMZ=atmfields(i).RND;
0451        
0452           
0453         end
0454             
0455         if strcmp(atmfields(i).source,'ecmwfftc')
0456 
0457            if ~exist('H')
0458 
0459               H=asg_cloudsat_ecmwf(CS.FILE,CS.LEG,CS.GRIDS,...
0460                                        {'DEM_elevation'});
0461 
0462            end
0463         end   
0464             
0465         if isnumeric(atmfields(i).source)
0466            
0467            G(i)=gf_set_fields(G(i),'NAME',atmfields(i).gas,...
0468                'SOURCE','set manually','DATA_NAME','Volume mixing ratio',...
0469                'DATA_UNIT','1');
0470            grids=[{G(1).GRID1} {} {}];
0471            data=atmfields(i).source*ones(length(grids{1}),1);
0472            G(i)= gf_set(  G(i),data,grids);
0473            G(i)=gf_set_fields(G(i),'GRID1_NAME',G1.GRID1_NAME,...
0474                 'GRID1_UNIT',G1.GRID1_UNIT);
0475            G(i).PROPS=atmfields(i).abs_models;
0476            G(i).RNDMZ=atmfields(i).RND;
0477        
0478         end
0479 
0480      end
0481 
0482      if strcmp(atmfields(i).source,'ecmwfftc')
0483      
0484         %merge data from ecmwf and fascod
0485         G(i)=merge_data2(G(i).NAME,G(i),H); 
0486      end
0487 
0488 end
0489 
0490 %H2O field
0491 ind1=find(strcmp('H2O',{G.NAME}));
0492 G(ind1).SPCFC.FIXED_RH=0;  
0493 G(ind1).SPCFC.MOD_RH=1;  
0494 
0495 
0496 %-------------------------------------------------------------------------------
0497 %CREATE FULLSKY SCENARIOUS
0498 workfolder =create_tmpfolder;
0499 
0500 warning off
0501 
0502 [Gfs,FS] = asg_create_gfs(G,CS.GRIDS,FS,workfolder);
0503 
0504 delete_tmpfolder(workfolder)
0505 warning on
0506 
0507 
0508 G=Gfs;
0509 
0510 
0511 
0512 %CHECK OUT IF THE CLOUDS ARE OVER LAND OR OCEAN
0513 %FS.L_OR_S=land_or_sea(FS.USE_LATS,FS.USE_LONS); %1 is over land
0514 
0515 
0516 
0517 %-------------------------------------------------------------------------------
0518 %sub-functions
0519 %
0520 
0521 function G=merge_data2(field,G,H)
0522 
0523     %merge 2-dim data (pressure x latitude (H))
0524     %with 1-dim data (pressure (G))
0525     %into 2 dim-data
0526     ind1=find(strcmp(field,{G.NAME}));
0527     ind2=find(strcmp(field,{H.NAME}));
0528     %find the min p_index of the 1-d data where
0529     %the 2-d data have no info
0530     m1=min(find(H(ind2).GRID1(end)>G(ind1).GRID1));
0531     %now merge the data
0532     data=G(ind1).DATA(m1:end)*ones(size(H(ind2).DATA,2),1)';
0533     %the new p_grid
0534     grid1=[H(ind2).GRID1',G(ind1).GRID1(m1:end)']';
0535     G(ind1).DIM=2;
0536     G(ind1)=gf_set(G(ind1),[H(ind2).DATA', data']',...
0537             [{grid1} {vec2col(H(ind2).GRID2)} ]);
0538     G(ind1)=gf_set_fields(G(ind1),'GRID1_NAME',H(ind2).GRID1_NAME,...
0539              'GRID1_UNIT',H(ind2).GRID1_UNIT,...
0540              'GRID2_NAME',H(ind2).GRID2_NAME,...
0541              'GRID2_UNIT',H(ind2).GRID2_UNIT);
0542  
0543 return
0544 
0545 %-------------------------------------------------------------------------------
0546 
0547 %Find out if we are over land or ocean
0548 function [l_or_s]=land_or_sea(lat_fs,lon_fs)
0549   [lat,lon,M]=land_sea_mask;
0550   for ip=1:length(lat_fs)
0551     if lon_fs(ip)<0
0552       l_data=round(360+lon_fs(ip));
0553     else
0554       l_data=round(lon_fs(ip));
0555     end
0556     lat_ind=find(round(lat_fs(ip))==lat);
0557     lon_ind=find(l_data==lon);
0558     l_or_s(ip)=M(lat_ind,lon_ind);
0559   end
0560 
0561 return

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