Home > atmlab > demos > asg_demo2.m

asg_demo2

PURPOSE ^

ASG_DEMO2 A simple demonstration how to use ASG

SYNOPSIS ^

function [G,FS]=asg_demo2(sensor,cs_file,ecmwf_file,leg,lat,lat_res,dlat)

DESCRIPTION ^

 ASG_DEMO2   A simple demonstration how to use ASG

    Gives an example on basic usage of ASG

    The example generates an atmospheric state
    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,FS]=asg_demo2(sensor,cs_file,ecmwf_file,leg,lat,lat_res,dlat)

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

 IN      
         sensor  the choice of instrument determines
              which gas species to include in G
              can be odin  or amsu 
      
         cs_file Cloudsat file from where to read data

         ecmwf_file Cloudsat aux-file from where to read data         

         leg  leg of cloudsat orbit to consider
              1=(-90,0);2=(-90,90);3=[0,90]

         lat  latitude of CloudSat data to consider 

         lat_res cloud lattude resolution to use

         dlat defines latitide/longitude extension of cloud

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

asg_demo2.m

SOURCE CODE ^

0001 % ASG_DEMO2   A simple demonstration how to use ASG
0002 %
0003 %    Gives an example on basic usage of ASG
0004 %
0005 %    The example generates an atmospheric state
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,FS]=asg_demo2(sensor,cs_file,ecmwf_file,leg,lat,lat_res,dlat)
0017 %
0018 % OUT     G    structure with atmospheric variables on gformat data
0019 %         FS   structure specifying some settings
0020 %
0021 % IN
0022 %         sensor  the choice of instrument determines
0023 %              which gas species to include in G
0024 %              can be odin  or amsu
0025 %
0026 %         cs_file Cloudsat file from where to read data
0027 %
0028 %         ecmwf_file Cloudsat aux-file from where to read data
0029 %
0030 %         leg  leg of cloudsat orbit to consider
0031 %              1=(-90,0);2=(-90,90);3=[0,90]
0032 %
0033 %         lat  latitude of CloudSat data to consider
0034 %
0035 %         lat_res cloud lattude resolution to use
0036 %
0037 %         dlat defines latitide/longitude extension of cloud
0038 %
0039 
0040 % Created by Bengt Rydberg
0041 function [G,FS]=asg_demo2(sensor,cs_file,ecmwf_file,leg,lat,lat_res,dlat)
0042 
0043 error( 'Can not be used presently. Not yet updated.' );
0044 
0045 %SET BASIC PROPERTIES
0046 
0047 %SET CLOUD VARIABLES
0048 %-
0049 
0050 %DEFINE FILE FROM WHERE TO READ DATA
0051 FS.DATA=cs_file;
0052 
0053 %DEFINE FROM WHAT LEG OF CLOUDSAT-ORBIT DATA WILL BE TAKEN FROM
0054 %Leg1=(-90,0);Leg2=(-90,90);Leg3=[0,90]
0055 FS.LEG=leg;
0056 
0057 %DEFINE FROM WHICH LATITUDES TO READ DATA
0058 FS.USE_LATS=lat;
0059 
0060 
0061 %DEFINE WHICH LATITUDE RESOLUTION TO USE
0062 FS.LAT_RES=lat_res; 
0063 
0064 %define pressure limits from where to consider Cloudsat data
0065 FS.P_LIMS=[900e2 80e2];
0066 dim=3;
0067 FS.DIM=dim;
0068 
0069 if dim==3
0070      
0071    %define the method to expand cloudsat data from 2-d to 3-d
0072    FS.DIMADD.METHOD='iaaft';
0073    
0074    %define how large latitude sections that will be considered
0075    FS.DLAT=dlat;
0076 
0077    %define longitude grid to use for cloud properties
0078    FS.LON_GRID=[-dlat/2:FS.LAT_RES:dlat/2];
0079 
0080 end
0081 
0082 %define properties for converting radar backscatter IWC/LWC fields
0083 f_grid=94e9;
0084 FS.PROPS=asg_ssp(f_grid);
0085 
0086 %END OF CLOUD VARIABLES SETTINGS
0087 
0088 %SET ATMOSPHERIC PROPERTIES (no cloud)
0089 %-
0090 
0091 %define from where to read data
0092 %the baseline here is that gas fields are read from fascod tropical
0093 %climatology(ftc)
0094 %however for H2O and O3 ECMWF data are here merged with Fascod data(ecmwfftc)
0095 %N2 is here fixed by giving it a numerical value
0096 %This in order to show some different options
0097 
0098 fascod     = fullfile( atmlab('ARTS_XMLDATA_PATH'),...
0099                            'atmosphere', 'fascod' );
0100 CS.FILE=ecmwf_file;
0101 
0102 %define from what leg of cloudsat-orbit data will be taken from
0103 CS.LEG=FS.LEG;
0104 
0105 CS.INSTRUMENT=sensor;
0106 
0107 
0108 %grids_cs will be the grids of the atmospheric variables (not for clouds)
0109 %however the grids will later be set tighter around the cloudbox
0110 %depending on the FS structure
0111 p_grid=z2p_simple([0:1e3:25e3 35e3 45e3 60e3]);
0112 CS.GRIDS{1}=p_grid; 
0113 CS.GRIDS{2}=sort([-90 90 FS.USE_LATS-10 FS.USE_LATS+10 ...
0114                  FS.USE_LATS-5 FS.USE_LATS+5 ...
0115                  FS.USE_LATS-dlat/2 FS.USE_LATS FS.USE_LATS+dlat/2]);
0116 
0117 if dim==3
0118   lon_grid=[FS.LON_GRID(1) 0 FS.LON_GRID(end)];
0119   CS.GRIDS{3}=unique(sort([-180 -10  -5 lon_grid 5 10 180]));
0120 end
0121 
0122 
0123 %mandatory fields to include
0124 
0125 fields=[{'temperature'},{'Altitude'},{'Surface altitude'}];
0126 atmfields(1).type='temperature';
0127 atmfields(1).source='ecmwfftc';
0128 atmfields(2).type='altitude';
0129 atmfields(2).source='';
0130 atmfields(3).type='surface altitude';
0131 atmfields(3).source='ecmwf';
0132 
0133 %define which gas types to include (depends on choice of instrument)
0134 
0135 if strcmp('odin',sensor)
0136 
0137    %gases to include for Odin
0138    gases=[{'H2O'},{'O3'},{'ClO'},{'N2O'},{'O2'},{'HNO3'},{'N2'}];
0139    %ABSORPTION MODELS
0140    abs_models{1}={'H2O','H2O-ForeignContStandardType',...
0141                   'H2O-SelfContStandardType'};
0142    abs_models{2}='O3';
0143    abs_models{3}='ClO';
0144    abs_models{4}='N2O';
0145    abs_models{5}={'O2','O2-PWR98'};
0146    abs_models{6}='HNO3';
0147    abs_models{7}='N2-SelfContMPM93';
0148     
0149 end 
0150  
0151 if strcmp('amsu',sensor)  
0152 
0153    %gases to include for amsu
0154    gases=[{'H2O'},{'O3'},{'O2'},{'N2'}];
0155 
0156    %absorption models amsu
0157    abs_models{1}='H2O-PWR98';
0158    abs_models{2}='O3';
0159    abs_models{3}='O2-PWR93';
0160    abs_models{4}='N2-SelfContStandardType';
0161 
0162 end
0163 
0164 
0165 len=length(atmfields);
0166 
0167 for i=1:length(gases) 
0168 
0169   atmfields(len+i).type='gas_species';
0170   atmfields(len+i).abs_models=abs_models{i};
0171   atmfields(len+i).gas=gases{i};
0172 
0173   if strcmp(gases{i},'H2O') | strcmp(gases{i},'O3')
0174 
0175      atmfields(len+i).source='ecmwfftc';
0176   
0177   elseif strcmp(gases{i},'N2')
0178 
0179      atmfields(len+i).source=0.7814;
0180   
0181   else
0182 
0183      atmfields(len+i).source='ftc'; 
0184    
0185   end
0186 
0187 end
0188 
0189 
0190 %Baseline for parameters describing how perturbations will be performed
0191 
0192 RND.FORMAT    = 'param'; 
0193 RND.SEPERABLE = 1;
0194 RND.CCO       = 0.01;               % Cut-off for correlation values
0195 RND.TYPE      = 'rel';              % Relative disturbances as default
0196 RND.DATALIMS  = [0];                % Do not allow negative values
0197 %
0198 RND.SI        = 0.1;                % 40% std. dev. as default
0199 %
0200 RND.CFUN1     = 'exp';              % Exp. correlation function for p-dim.
0201 RND.CL1       = [0.15 0.3 0.3]';    % Corr. length varies with altitude
0202 RND.CL1_GRID1 = [1100e2 10e2 1e-3];    
0203 %
0204 RND.CFUN2     = 'lin';              % Linear correlation function for lat-dim.
0205 RND.CL2       = 0.5;                % Corr. length 0.5 deg everywhere
0206 %
0207 RND.CFUN3     = 'lin';              % Linear correlation function for lat-dim.
0208 RND.CL3       = 0.5;                % Corr. length 0.5 deg everywhere
0209 
0210 %describe perturbations for each variable
0211 for i=1:length(atmfields)
0212   
0213   if strcmp(atmfields(i).type,'temperature')
0214 
0215      atmfields(i).RND=RND;
0216      atmfields(i).RND.TYPE='abs';
0217      atmfields(i).RND.SI=1;
0218    
0219   elseif strcmp(atmfields(i).gas,'H2O')
0220 
0221      atmfields(i).RND=RND;
0222      atmfields(i).RND.SI=0.1;
0223  
0224   elseif strcmp(atmfields(i).gas,'O3')
0225      
0226      atmfields(i).RND=RND;
0227      atmfields(i).RND.SI=0.2;
0228 
0229   elseif strcmp(atmfields(i).gas,'ClO') | strcmp(atmfields(i).gas,'N2O') | ...
0230          strcmp(atmfields(i).gas,'HNO3') 
0231 
0232      atmfields(i).RND=RND;
0233   
0234   else
0235      
0236      atmfields(i).RND=[];
0237  
0238   end
0239 
0240 end
0241  
0242 %Read basic data
0243 
0244 
0245 clim='subarctic-summer';
0246 clear G H
0247 for i=1:length(atmfields)
0248 
0249      %D=asgD;
0250      
0251      G(i)=asgG;
0252 
0253      %G = gf_empty( 3 )
0254 
0255 
0256 
0257     if strcmp(atmfields(i).type,'temperature')
0258 
0259        if strcmp(atmfields(i).source,'ftc') | ...
0260           strcmp(atmfields(i).source,'ecmwfftc') 
0261           %G1 = gfin_artsxml( asgD, [], ...
0262           %                   fullfile( fascod, [clim,'.t.xml'] ), ...
0263           %                   'Temperature field','t_field' );
0264           G1 = gf_artsxml(fullfile( fascod, [clim,'.t.xml'] ), ...
0265                              'Temperature field','t_field' );
0266           
0267  
0268        end
0269          
0270        if strcmp(atmfields(i).source,'ecmwfftc') 
0271        
0272           H=asg_cloudsat_ecmwf(CS.FILE,CS.LEG,CS.GRIDS,{'DEM_elevation'});
0273          
0274           for j=1:length(H)
0275               if strcmp(H(j).NAME,'Water vapour')          
0276                 H(j).NAME='H2O';
0277               end
0278               if strcmp(H(j).NAME,'Ozone')          
0279                 H(j).NAME='O3';
0280               end
0281           end
0282    
0283        end
0284 
0285      end
0286   
0287      if strcmp(atmfields(i).type,'altitude')
0288      
0289         %altitude field is here just initialised
0290         G(i).NAME            = 'Altitude field'; 
0291         G(i).DATA_NAME       = 'Altitude field';
0292         G(i).DATA_UNIT       = 'm';
0293         %G(i).SPCFC.P0       = 1013e2;
0294         %G(i).SPCFC.P0       = G(1).GRID1(1);
0295         G(i).SPCFC.P0        =CS.GRIDS{1}(1);
0296         grids=[{G(1).GRID1(1)} {0} {0} {0}]; 
0297         %G(i).SPCFC.Z0=gf_set( asgD, asgG,grids , 0, [], ...
0298         %                       'Altitude at P0', 'Altitude', 'm' );
0299         G(i).SPCFC.Z0=gf_set( asgG,0,grids,[],[]);
0300         G(i).SPCFC.Z0.DIMS=[];
0301         G(i).SPCFC.Z0.GRID1=[];
0302         G(i).SPCFC.Z0.GRID2=[];
0303         G(i).SPCFC.Z0.GRID3=[];
0304         G(i).SPCFC.Z0.GRID4=[];
0305 
0306      end 
0307        
0308      if strcmp(atmfields(i).type,'surface altitude')
0309           
0310         if strcmp(atmfields(i).source,'ecmwf')
0311           
0312            ind=find(strcmp('DEM_elevation',{H.NAME}));
0313            G(i)=H(ind);
0314            G(i).DATA_NAME='Surface altitude';
0315            G(i).NAME='Surface altitude';
0316            G(i).DIMADD.METHOD='expand';
0317                    
0318         end
0319        
0320      end
0321             
0322      
0323      if strcmp(atmfields(i).type,'gas_species')
0324 
0325         if strcmp(atmfields(i).source,'ftc') |  ...
0326            strcmp(atmfields(i).source,'ecmwfftc')               
0327            sourcefile= fullfile( fascod, [clim,'.',...
0328                                  atmfields(i).gas,'.xml'] );
0329            %G1= gfin_artsxml( asgD, [],sourcefile,atmfields(i).gas,...
0330            %                  'vmr_field' );
0331 
0332            G1= gf_artsxml( sourcefile,atmfields(i).gas,...
0333                             'vmr_field' );
0334 
0335            
0336            %G1 = gf_artsxml( fullfile( fascod, [clim,'.t.xml'] ), ...
0337            %                  'Temperature field','t_field' );
0338          
0339           
0340         end
0341             
0342         if strcmp(atmfields(i).source,'ecmwfftc')
0343 
0344            if ~exist('H')
0345 
0346               H=asg_cloudsat_ecmwf(CS.FILE,CS.LEG,CS.GRIDS,...
0347                                        {'DEM_elevation'});
0348            end
0349         end   
0350             
0351         if isnumeric(atmfields(i).source)
0352 
0353            grids=[{G(1).GRID1} {0} {0} {0}];
0354            data=atmfields(i).source*ones(length(grids{1}),1);
0355            %G(i)= gf_set( asgD, asgG, grids, data, [], ...
0356            %              atmfields(i).gas, 'Volume mixing ratio', '-' );
0357            %keyboard
0358            G(i)= gf_set( gf_empty( 1 ),data,grids,[],[]);
0359            G1= gf_set(  asgG,data,grids,[],[]);
0360   
0361          %G(i).SPCFC.Z0=gf_set( asgG,0,grids,[],[]);
0362        
0363         end
0364 
0365      end
0366         
0367      if strcmp(atmfields(i).source,'ftc') |  ...
0368         strcmp(atmfields(i).source,'ecmwfftc')   
0369 
0370          G(i).NAME=G1.NAME;
0371          G(i).DIMS=1;
0372          G(i).DATA=G1.DATA;
0373          G(i).DATA_NAME=G1.DATA_NAME;
0374          G(i).DATA_UNIT=G1.DATA_UNIT;
0375          G(i).GRID1=G1.GRID1;
0376          %G(i).GRID2=G1.GRID2;
0377          %G(i).GRID3=G1.GRID3;
0378          G(i).SOURCE=G1.SOURCE;
0379             
0380      end
0381      
0382      if strcmp(atmfields(i).source,'ecmwfftc')
0383      
0384         %merge data from ecmwf and fascod
0385         G(i)=merge_data2(asgD,G(i).NAME,G(i),H);
0386 
0387      end
0388  
0389  
0390     G(i).PROPS=atmfields(i).abs_models;
0391 
0392      
0393     G(i).RNDMZ=atmfields(i).RND;
0394 
0395 end
0396 
0397 
0398 keyboard
0399 
0400 
0401 %rename H2O field
0402 ind1=find(strcmp('H2O',{G.NAME}));
0403 G(ind1).NAME='Water vapour';
0404 G(ind1).SPCFC.FIXED_RH=0;  
0405 G(ind1).SPCFC.MOD_RH=0;  
0406 
0407 %-------------------------------------------------------------------------------
0408 %CREATE FULLSKY SCENARIOUS
0409 workfolder =create_tmpfolder;
0410 
0411 warning off
0412 [Gfs,lon] = asg_create_gfs(G,CS.GRIDS,FS,workfolder);
0413 delete_tmpfolder(workfolder)
0414 warning on
0415 
0416 FS.USE_LONS=lon;
0417 G=Gfs;
0418 %CHECK OUT IF THE CLOUDS ARE OVER LAND OR OCEAN
0419 %FS.L_OR_S=land_or_sea(FS.USE_LATS,FS.USE_LONS); %1 is over land
0420 
0421 
0422 
0423 %-------------------------------------------------------------------------------
0424 %sub-functions
0425 %
0426 
0427 function G=merge_data2(D,field,G,H)
0428  
0429     S=G;
0430     ind1=find(strcmp(field,{G.NAME}));
0431     ind2=find(strcmp(field,{H.NAME}));
0432     lims1{1}=[H(ind2).GRID1([1 end])];
0433     lims1{2}=[H(ind2).GRID2([1 end])];
0434     m1=min(find(H(ind2).GRID1(end)>G(ind1).GRID1));
0435     lims2{1}=[G(ind1).GRID1(m1) G(ind1).GRID1(end)];
0436     addpath /home/bengt/ARTS/atmlab_sep08/atmlab/gformat
0437     G(ind1)=gf_merge(asgD,H(ind2),lims1,asgD,G(ind1),lims2);
0438     rmpath /home/bengt/ARTS/atmlab_sep08/atmlab/gformat
0439     G(ind1).DIMADD=S(ind1).DIMADD;
0440 
0441 return
0442 
0443 %-------------------------------------------------------------------------------
0444 
0445 function P=asg_ssp(f_grid)
0446  %properties for adding pnd fields
0447  alpha=0;Ni=10;xnorm=50e-6;
0448  [x_i,w_i]=gauss_laguerre(alpha,Ni,xnorm);
0449  Nj=3;xnormw=10e-6;
0450  [x_k,w_k]=gauss_laguerre(alpha,Nj,xnormw);
0451  %x_i will be the diameter of the particles considered
0452  %calculate ssp
0453  for i=1:2
0454 
0455   if i==1 
0456    F=create_ssp('sphere',f_grid,x_i/2);
0457    P(i).x=x_i;
0458    P(i).w=w_i;
0459    P(i).x_norm=xnorm;
0460    P(i).PSD='MH97';
0461   else
0462    T_grid=[273 283 293 303 313];
0463    for i1=1:length(f_grid)
0464      for j1=1:length(T_grid)
0465        rfr_index(i1,j1) = sqrt(eps_water_liebe93( f_grid(i1), T_grid(j1) ));
0466      end
0467    end
0468    F=create_ssp('sphere',f_grid,x_k,T_grid,rfr_index);
0469    P(i).x=x_k;
0470    P(i).w=w_k;
0471    P(i).x_norm=xnormw;
0472    P(i).PSD='Water';
0473   end
0474   P(i).SSP=F;
0475   P(i).method='gauss-laguerre';
0476   P(i).shape='sphere';  
0477  end
0478 
0479 return
0480 
0481 
0482 %-------------------------------------------------------------------------------
0483 
0484 %Find out if we are over land or ocean
0485 function [l_or_s]=land_or_sea(lat_fs,lon_fs)
0486   [lat,lon,M]=land_sea_mask;
0487   for ip=1:length(lat_fs)
0488     if lon_fs(ip)<0
0489       l_data=round(360+lon_fs(ip));
0490     else
0491       l_data=round(lon_fs(ip));
0492     end
0493     lat_ind=find(round(lat_fs(ip))==lat);
0494     lon_ind=find(l_data==lon);
0495     l_or_s(ip)=M(lat_ind,lon_ind);
0496   end
0497 
0498 return

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