Home > atmlab > demos > asg2y_demo.m

asg2y_demo

PURPOSE ^

ASG2Y_DEMO A simple demonstration how to use

SYNOPSIS ^

function [G,y,ydata,dy]=asg2y_demo(Q,varargin)

DESCRIPTION ^

 ASG2Y_DEMO   A simple demonstration how to use
    asg functions to simulate radiative transfer
    for generated atmospheric states (by asg_demo)

    Gives an example on basic usage of ASG2Y

    The example simulate radiative transfer 
    for a cloudy state generated by asg_demo 

    A lot of default settings are made in this
    demo, which may be re-considered depending
    on what data should be used for

 FORMAT [G,y,ydata,dy]=asg2y_demo( Q, [ instrument ])

 OUT     G      structure with atmospheric variables on gformat data
         y      As returned by *arts_y*.
         ydata  As returned by *arts_y*.
         dy     As returned by *arts_y*.        

 IN      Q       qarts setting structure

 OPT     sensor  choice of instrument 
                 can be odin (default) or amsu

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

asg2y_demo.m

SOURCE CODE ^

0001 % ASG2Y_DEMO   A simple demonstration how to use
0002 %    asg functions to simulate radiative transfer
0003 %    for generated atmospheric states (by asg_demo)
0004 %
0005 %    Gives an example on basic usage of ASG2Y
0006 %
0007 %    The example simulate radiative transfer
0008 %    for a cloudy state generated by asg_demo
0009 %
0010 %    A lot of default settings are made in this
0011 %    demo, which may be re-considered depending
0012 %    on what data should be used for
0013 %
0014 % FORMAT [G,y,ydata,dy]=asg2y_demo( Q, [ instrument ])
0015 %
0016 % OUT     G      structure with atmospheric variables on gformat data
0017 %         y      As returned by *arts_y*.
0018 %         ydata  As returned by *arts_y*.
0019 %         dy     As returned by *arts_y*.
0020 %
0021 % IN      Q       qarts setting structure
0022 %
0023 % OPT     sensor  choice of instrument
0024 %                 can be odin (default) or amsu
0025 
0026 % Created by Bengt Rydberg
0027 
0028 
0029 function [G,y,ydata,dy]=asg2y_demo(Q,varargin)
0030 %
0031 [sensor] = optargs( varargin, { 'odin' } );
0032 
0033 error( 'Can not be used presently. Not yet updated.' );
0034 
0035 
0036 %-----------------------------------
0037 %generate cloudy atmospheric states
0038 
0039 [G,FS]=asg_demo( 3,sensor );
0040 G=G{1};
0041 
0042 %-----------------------------------
0043 %SET Q-STUCTURE (for ARTS RT-simulations)
0044 
0045 if ~exist('Q')
0046   %in sub-function set_q a lot of default settings are made
0047   Q=asg_q(sensor);
0048 end
0049  
0050 
0051 %------------------------------------
0052 %ADD SSP OF ICE/WATER PARTICLES FOR THE RELEVANT INSTRUMENT FREQUENCIES
0053 FS.PROPS=asg_ssp(Q.F_GRID);
0054     
0055 ind=find(strncmp({G.NAME},'Particles',9));
0056 for j=ind
0057  G(j).PROPS=FS.PROPS(1).SSP(j-ind(1)+1);
0058 end
0059   
0060 %---------------------------------------
0061 %RUN RTE CALCULATIONS
0062  
0063 Q.P_GRID=G(1).GRID1;
0064 Q.LAT_GRID=G(1).GRID2;
0065 Q.LON_GRID=G(1).GRID3;
0066 
0067 %adjust sensor_pos
0068 if strcmp('amsu',sensor)
0069   
0070    %the sensor is placed right above the cloud
0071    Q.SENSOR_POS(:,2)=FS.USE_LATS(1);
0072   
0073    %randomize a surface emissivity
0074    emis=min(0.95+0.05*randn(length(Q.F_GRID),1),1);
0075    emis=max(0.8,emis);
0076    FS.SURF_EMIS=emis;
0077 
0078 end
0079   
0080 if strcmp('odin',sensor)
0081 
0082    %make sure the sensor looks into the cloud field
0083    if FS.LEG==1
0084       Q.SENSOR_POS(:,2)=-Q.SENSOR_POS(:,2)+FS.USE_LATS(1);
0085       Q.SENSOR_LOS(:,2)=180; 
0086    elseif FS.LEG==2
0087       Q.SENSOR_POS(:,2)=Q.SENSOR_POS(:,2)+FS.USE_LATS(1); 
0088       Q.SENSOR_LOS(:,2)=0; 
0089    else
0090       Q.SENSOR_POS(:,2)=-Q.SENSOR_POS(:,2)+FS.USE_LATS(1);
0091       Q.SENSOR_LOS(:,2)=180;
0092    end
0093   
0094 end
0095 
0096 
0097 %run arts simulations
0098 
0099 workfolder =create_tmpfolder;
0100 
0101 %don't consider water particles here/ should be updated
0102 ind=find(strncmp({G.NAME},'LWC',3));
0103 G=G(1:ind-1);
0104 
0105 [y,dy,ydata]=asg_gfs2i(G,Q,FS,sensor,workfolder)
0106 
0107 delete_tmpfolder(workfolder)
0108 
0109 
0110 
0111 
0112 
0113 %-sub-functions
0114 
0115 
0116 function Q=asg_q(instrument)
0117 
0118  Q = qarts;
0119  Q.OUTPUT_FILE_FORMAT='binary';
0120 
0121  %- Basic atmospheric variables
0122  %
0123  Q.ATMOSPHERE_DIM = 3;
0124  Q.P_GRID         = 100/16e3; % Unit here is pressure decade
0125  Q.LAT_GRID=[];
0126  Q.LON_GRID=[];
0127  Q.R_GEOID=6366700;
0128 
0129  %- Surface properties
0130  %
0131  if strcmp('amsu',instrument)
0132   se=0.6;
0133   Q.SURFACE_PROP_AGENDA = ...
0134   {'InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim, p_grid,lat_grid, lon_grid,rte_gp_p, rte_gp_lat, rte_gp_lon,t_field )', ...
0135   'Ignore(rte_los)',...
0136   ['NumericSet(surface_emissivity,',num2str(se),')'],'surfaceSimple'};
0137  end
0138 
0139  if strcmp('odin',instrument) 
0140 
0141    Q.SURFACE_PROP_AGENDA{1}=...
0142       ['InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim,'...
0143        'p_grid, lat_grid, lon_grid,rte_gp_p, rte_gp_lat,'...
0144        'rte_gp_lon,t_field)'];
0145    Q.SURFACE_PROP_AGENDA{2}='Ignore(rte_pos)';
0146    Q.SURFACE_PROP_AGENDA{3}='Ignore(rte_los)';
0147    Q.SURFACE_PROP_AGENDA{4}=['surfaceBlackbody(surface_los,'...
0148       'surface_rmatrix,surface_emission, f_grid, stokes_dim,surface_skin_t )'];
0149  
0150  end
0151 
0152  %- Spectroscopic stuff
0153  %
0154  Q.STOKES_DIM       = 1;
0155  Q.ABS_LINES_FORMAT = 'Arts';
0156  Q.ABSORPTION='OnTheFly';
0157  Q.ABS_LINESHAPE='Voigt_Kuntz6';
0158  Q.ABS_LINESHAPE_CUTOFF=-1;
0159  Q.ABS_LINESHAPE_FACTOR='quadratic';
0160 
0161 
0162  if strcmp('amsu',instrument)
0163   f_grid_file=fullfile( atmlab('ARTS_INCLUDES'),'amsu/amsub.f_grid_fast.xml');
0164   f_grid=xmlLoad(f_grid_file);
0165   Q.F_GRID           = f_grid;
0166   Q.ABS_LINES=fullfile(atmlab('ARTS_INCLUDES'),'amsu/amsub.hitran04_lines.xml');
0167   Q.ABS_MODELS={fullfile(atmlab('ARTS_INCLUDES'),'continua.arts')};
0168  end
0169 
0170  if strcmp('odin',instrument)
0171   Q.F_GRID           = [501.2e9 544.4e9];
0172   Q.ABS_LINES{1}=...
0173      fullfile(atmlab('ARTS_INCLUDES'),'odin/linefile.SM_AC2ab.xml'); 
0174   Q.ABS_LINES{2}=...
0175      fullfile(atmlab('ARTS_INCLUDES'),'odin/linefile.SM_AC1e.xml');
0176   Q.ABS_MODELS={fullfile(atmlab('ARTS_INCLUDES'),'odin/odinsmr.arts')};
0177  end
0178 
0179  %- RTE variables
0180  %
0181  Q.PPATH_STEP_AGENDA = { 'ppath_stepGeometric( ppath_step, atmosphere_dim, p_grid, lat_grid, lon_grid, z_field, r_geoid, z_surface, ppath_lmax)' };
0182 Q.PPATH_LMAX=1e3;
0183 Q.RTE_AGENDA= {'RteStd(iy,diy_dvmr,diy_dt,ppath,ppath_array,ppath_array_index,f_grid,stokes_dim,emission_agenda,abs_scalar_gas_agenda,rte_do_vmr_jacs,rte_do_t_jacs )'};
0184  Q.EMISSION_AGENDA={'emissionPlanck'};
0185  Q.IY_SPACE_AGENDA={'Ignore(rte_los)','Ignore(rte_pos)',...
0186               'MatrixCBR(iy,stokes_dim,f_grid)' };
0187  Q.J_DO=0;
0188 
0189  %- SENSOR variabels
0190  Q.ANTENNA_DIM=1;
0191  Q.IY_UNIT            = 'RJBT';
0192  Q.SENSOR_DO=0;
0193 
0194 
0195  %- Sensor POS and LOS
0196  %
0197  % Here hard coded
0198  %
0199  if strcmp('amsu',instrument)
0200    re=Q.R_GEOID;
0201    zplat = 850e3;
0202    xoff=[-7.5e3:2.5e3:7.5e3];
0203    yoff=[-7.5e3:2.5e3:7.5e3];
0204    xoff=repmat(xoff,length(xoff),1);
0205    yoff=repmat(yoff,1,length(yoff));
0206    xoff=xoff(:);
0207    yoff=yoff(:);
0208    s=sqrt(xoff.^2+yoff.^2);
0209    za=180-180/pi*atan(s/zplat);
0210    th=180/pi*atan(yoff./xoff);
0211    th(find(isnan(th)))=0;
0212    ivec=find(xoff==0 & yoff<0);
0213    th(ivec)=th(ivec)+360;
0214    ivec=find(yoff==0 & xoff<0);
0215    th(ivec)=180;
0216    ivec=find(xoff<0 & yoff>0);
0217    th(ivec)=th(ivec)+180;
0218    ivec=find(xoff<0 & yoff<0);
0219    th(ivec)=th(ivec)+180;
0220    ivec=find(xoff>0 & yoff<0);
0221    th(ivec)=th(ivec)+360;
0222    ivec=find(th>180);
0223    th(ivec)=th(ivec)-360;
0224    ppos=zeros(length(th),3);
0225    ppos(:,1)=re+zplat;
0226    plos=zeros(length(th),2);
0227    plos(:,1)=za;
0228    plos(:,2)=th;
0229    Q.SENSOR_POS=ppos;
0230    Q.SENSOR_LOS=plos;
0231  end
0232 
0233  if strcmp('odin',instrument)
0234    re=Q.R_GEOID;
0235    zplat = 600e3;
0236    tan=[-5e3 -2e3:5e2:-500 0:2.5e2:10e3 10.5e3:5e2:12e3 15e3]';
0237    tan=6e3;
0238    za = geomztan2za(re,zplat,tan);
0239    ppos=zeros(length(za),3);
0240    ppos(:,1)=re+zplat;
0241    ppos(:,2)=-21.5;
0242    plos=zeros(length(za),2);
0243    plos(:,1)=za;
0244    Q.SENSOR_POS=ppos;
0245    Q.SENSOR_LOS=plos;
0246  end
0247 
0248  %- MC variables
0249  %
0250  Q.CLOUDBOX.METHOD                      = 'MC';
0251  Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR       = 3;      
0252  Q.CLOUDBOX.METHOD_PRMTRS.MAX_TIME      = -1;
0253  Q.CLOUDBOX.METHOD_PRMTRS.MAX_ITER      = -1;
0254  Q.CLOUDBOX.METHOD_PRMTRS.Z_FIELD_IS_1D = 0;
0255  Q.CLOUDBOX.LIMITS                      = [];
0256  Q.CLOUDBOX.OPT_PROP_GAS_AGENDA={'ext_matInit','abs_vecInit',...
0257                                 'ext_matAddGas','abs_vecAddGas'};
0258  Q.CLOUDBOX.OPT_PROP_PART_AGENDA={'ext_matAddPart','abs_vecAddPart'};
0259  Q.CLOUDBOX.SPT_CALC_AGENDA={};
0260  Q.CLOUDBOX_DO=0;
0261 
0262 return
0263 
0264 %-------------------------------------------------------------------------
0265 
0266 function P=asg_ssp(f_grid)
0267  %properties for adding pnd fields
0268  alpha=0;Ni=10;xnorm=50e-6;
0269  [x_i,w_i]=gauss_laguerre(alpha,Ni,xnorm);
0270  Nj=3;xnormw=10e-6;
0271  [x_k,w_k]=gauss_laguerre(alpha,Nj,xnormw);
0272  %x_i will be the diameter of the particles considered
0273  %calculate ssp
0274  for i=1:2
0275 
0276   if i==1 
0277    F=create_ssp('sphere',f_grid,x_i/2);
0278    P(i).x=x_i;
0279    P(i).w=w_i;
0280    P(i).x_norm=xnorm;
0281    P(i).PSD='MH97';
0282   else
0283    T_grid=[273 283 293 303 313];
0284    for i1=1:length(f_grid)
0285      for j1=1:length(T_grid)
0286        rfr_index(i1,j1) = sqrt(eps_water_liebe93( f_grid(i1), T_grid(j1) ));
0287      end
0288    end
0289    F=create_ssp('sphere',f_grid,x_k,T_grid,rfr_index);
0290    P(i).x=x_k;
0291    P(i).w=w_k;
0292    P(i).x_norm=xnormw;
0293    P(i).PSD='Water';
0294   end
0295   P(i).SSP=F;
0296   P(i).method='gauss-laguerre';
0297   P(i).shape='sphere';  
0298  end
0299 
0300 return
0301 
0302 
0303 %-------------------------------------------------------
0304 function [y,dy,ydata]=asg_gfs2i(G,Q,FS,sensor,workfolder)
0305 
0306  %- Determine necessary size of cloud box and adjust it
0307  %
0308  [Q,G]=asg_set_cb(Q,G);
0309 
0310 
0311  %check if there is any liquid cloud
0312  %if not we don't need to include it in the simulations
0313  ind=find(strcmp({G.NAME},'LWC field'));
0314  if ~isempty(ind)
0315    ind=find(strcmp({G.NAME},'LWC field'));
0316    G(ind).PROPS=['liquidcloud-MPM93'];
0317    G(ind).DATA=G(ind).DATA*1e-3;
0318    G(ind).DATA_NAME='Volume mixing ratio';
0319    G(ind)= asg_regrid( asgD, G(ind), Q );
0320    G(ind).DATA(find(G(ind).DATA(:)>10e-3))=10e-3;
0321  else
0322    ind=setdiff(1:length(G),ind);
0323    G=G(ind);
0324  end
0325 
0326  %- Run ARTS
0327  %
0328  if strcmp(sensor,'amsu')
0329    
0330 
0331   %SURFACE EMISSIVITIES
0332   Se=[];
0333   for j=1:length(FS.SURF_EMIS)
0334     se=num2str(FS.SURF_EMIS(j));
0335     Se=[Se,',',se];
0336   end         
0337   Se=Se(2:end);
0338 
0339   Q.SURFACE_PROP_AGENDA= ...
0340     {['InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim,'...
0341       'p_grid, lat_grid, lon_grid,rte_gp_p, rte_gp_lat,'...
0342       'rte_gp_lon,t_field)'],...
0343       'Ignore(rte_pos)',...
0344       ['VectorCreate(surface_emissivity)'],... 
0345      ['VectorSet(surface_emissivity,[',Se,'])'],...
0346       ['surfaceFlatVaryingEmissivity(surface_los, surface_rmatrix,'...
0347       'surface_emission, f_grid, stokes_dim,atmosphere_dim, rte_los,'...
0348       'surface_skin_t,surface_emissivity)']};  
0349    Q = asg2q( asgD, G, Q, workfolder );
0350    Q.ABS_NLS=[];
0351    %create an include file
0352    fid=fopen(fullfile(workfolder,'include.arts'))
0353    if fid==-1
0354      fid=fopen(fullfile(workfolder,'include.arts'),'a');
0355      fprintf(fid,'Arts2 {\n'); 
0356      fprintf(fid,'#VectorSetExplicitly(abs_n2,[0.7808])\n');
0357      fprintf(fid,'VectorSet( abs_n2, [ 0.7808 ] )\n');
0358      fprintf(fid,'IndexSet( abs_p_interp_order,   5 )\n');
0359      fprintf(fid,'IndexSet( abs_t_interp_order,   7 )\n');
0360      fprintf(fid,'IndexSet( abs_nls_interp_order, 5 )}\n');
0361      fclose(fid);
0362    end
0363    Q.INCLUDES= {fullfile(workfolder,'include.arts')};
0364    Q                  = qarts_abstable( Q,9,100,[],[] );
0365    [y,ydata,dy] = arts_y( Q, workfolder );
0366 
0367  end
0368 
0369 
0370  if strcmp(sensor,'odin')
0371 
0372   for i= 1:length(Q.F_GRID)
0373     if 0
0374      if i==1
0375        cind=find(strncmp(lower({G.NAME}),'hno3',4));
0376        ind=setdiff(1:length(G),cind);
0377        G1=G(ind);
0378      else
0379        cind=find(strncmp(lower({G.NAME}),'clo',3));
0380        ind=setdiff(1:length(G),cind);
0381        G1=G(ind);
0382      end
0383    else
0384     G1=G;
0385    end
0386     
0387     Q1           = asg2q( asgD, G1, Q, workfolder );
0388     Q1.ABS_NLS    =[];
0389     Q1           = qarts_abstable( Q1,9,100,[],[] );
0390     
0391     %create an include file
0392     fid=fopen(fullfile(workfolder,'include.arts'))
0393     if fid==-1
0394      fid=fopen(fullfile(workfolder,'include.arts'),'a');
0395      fprintf(fid,'Arts2 {\n'); 
0396      fprintf(fid,'#VectorSetExplicitly(abs_n2,[0.7808])\n');
0397      fprintf(fid,'VectorSet( abs_n2, [ 0.7808 ] )\n');
0398      fprintf(fid,'IndexSet( abs_p_interp_order,   5 )\n');
0399      fprintf(fid,'IndexSet( abs_t_interp_order,   7 )\n');
0400      fprintf(fid,'IndexSet( abs_nls_interp_order, 5 )}\n');
0401      fclose(fid);
0402     end
0403     Q1.INCLUDES= {fullfile(workfolder,'include.arts')};
0404     Q1.ABS_LINES =Q.ABS_LINES{i};
0405     Q1.F_GRID    =Q.F_GRID(i);
0406     [y(:,i),ydata{i},dy(:,i)] = arts_y( Q1, workfolder );
0407   
0408   end
0409 
0410  end
0411 
0412 
0413 return
0414 
0415 %-----------------------------------------------------------------
0416 
0417 function [Q,G]=asg_set_cb(Q,G)
0418 
0419 iwc_ind = find( strcmp( lower({G.NAME}), 'iwc field' ) );
0420 [C1,C2,C3,C4,C5,C6]=deal([]);
0421 
0422 for i=1:size(G(iwc_ind).DATA,3)
0423   [c1,c2]=find(G(iwc_ind).DATA(:,:,i));
0424   if ~isempty(c1)
0425      C1=min([C1,vec2col(c1)']);
0426      C2=max([C2,vec2col(c1)']);
0427   end
0428   if ~isempty(c2)
0429      C3=min([C3,vec2col(c2)']);
0430      C4=max([C4,vec2col(c2)']);
0431   end
0432 end
0433 
0434 
0435 for i=1:size(G(iwc_ind).DATA,1)
0436   [c1,c2]=find(squeeze(G(iwc_ind).DATA(i,:,:)));
0437     if ~isempty(c2)
0438        C5=min([C5,vec2col(c2)']);
0439        C6=max([C6,vec2col(c2)']);
0440    end
0441 end
0442 %
0443 
0444 if ~isempty(C1) %if C1 is empty we don't have any cloud particles
0445     p1   = G(iwc_ind).GRID1(C1);
0446     p2   = G(iwc_ind).GRID1(C2);
0447     lat1 = G(iwc_ind).GRID2(C3);
0448     lat2 = G(iwc_ind).GRID2(C4);
0449     z    = p2z_cira86([p1 p2]',mean([lat1 lat2]),160);
0450     lon1 = G(iwc_ind).GRID3(C5);
0451     lon2 = G(iwc_ind).GRID3(C6);
0452     %
0453     Q.CLOUDBOX.LIMITS = [z(1)-2e3 z(2)+2e3 lat1 lat2 lon1 lon2];
0454     
0455 else
0456     %remove all particles in G
0457     ind  =find(strncmp({G.NAME},'Particles',9)==0);
0458     G    = G(ind);
0459 
0460 end

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