Home > atmlab > arts > scenegen > asg_create_gfs.m

asg_create_gfs

PURPOSE ^

ASG_CREATE_GFS Creates atmopsheric states including clouds

SYNOPSIS ^

function [H,C] = asg_create_gfs(Gcs,grids_cs,C,workfolder)

DESCRIPTION ^

 ASG_CREATE_GFS   Creates atmopsheric states including clouds 
                  based on various input data 

    Creates atmospheric states including clouds, where the spatial
    structure of the clouds is taken from radar data.

    Reads in Cloudsat radar data from selected parts
    of the orbits based on the fields of *C*.
    
    Prepares data for the sub-function *asg_gen_state*
    where realisations of atmospheric states are performed,
    which can either be 1 or 3-d.

    FORMAT [H,C] = asg_create_gfs(Gcs,grids_cs,C,workfolder)

 OUT   H           Modified G format data
       C           Modified C data 
 IN    Gcs         Specification of clear sky fields, in gformat.
       grids_cs    Initial grids for clear sky fields. These grids are used
                   as put into Q before calling *asg_atmgrids*. Obtained
                   grids are used e.g. when calling *asg_rndmz* for clear 
                   sky fields. 
                   The input is given as an array of vectors, where element  
                   1 corresponds to Q.P_GRID, element 2 Q.LAT_GRID and 
                   element 3 Q.LON_GRID. 
       C           A cloud definition structure with fields
         DATA           name of datafile to be read in
         LEG            leg of cloudsat data 
         USE_LATS       latitude midpoints
         DLAT           latitude (scalar)
         LAT_RES        latitude resolution (scalar)
         LON_GRID       longitude grid 
         C.P_LIMS       pressure limits
         DIM            desired dimensionality of data
         DIMADD         (see asg_dimadd)
         PROPS
                    a psd and ssp defintiion structure 
                    C.PROPS(1) corresponds to the radar frequency             
                    C.PROPS(2) corresponds to the instrument frequencys 
                    with fields                    
             radar_back    single scattering properties data
             PSD    particle size distribution, only 'MH97' and 'Water' works
             method psd method, only 'gauss-laguerre' works
             x      particle sizes (see help gauss_laguerre)
             w      weights
             x_norm normalistaion factor 
             shape  particle shape, only 'sphere' works 


        Examle usage:                
        visit asg_demo to see an example how to use this function

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

asg_create_gfs.m

SOURCE CODE ^

0001 % ASG_CREATE_GFS   Creates atmopsheric states including clouds
0002 %                  based on various input data
0003 %
0004 %    Creates atmospheric states including clouds, where the spatial
0005 %    structure of the clouds is taken from radar data.
0006 %
0007 %    Reads in Cloudsat radar data from selected parts
0008 %    of the orbits based on the fields of *C*.
0009 %
0010 %    Prepares data for the sub-function *asg_gen_state*
0011 %    where realisations of atmospheric states are performed,
0012 %    which can either be 1 or 3-d.
0013 %
0014 %    FORMAT [H,C] = asg_create_gfs(Gcs,grids_cs,C,workfolder)
0015 %
0016 % OUT   H           Modified G format data
0017 %       C           Modified C data
0018 % IN    Gcs         Specification of clear sky fields, in gformat.
0019 %       grids_cs    Initial grids for clear sky fields. These grids are used
0020 %                   as put into Q before calling *asg_atmgrids*. Obtained
0021 %                   grids are used e.g. when calling *asg_rndmz* for clear
0022 %                   sky fields.
0023 %                   The input is given as an array of vectors, where element
0024 %                   1 corresponds to Q.P_GRID, element 2 Q.LAT_GRID and
0025 %                   element 3 Q.LON_GRID.
0026 %       C           A cloud definition structure with fields
0027 %         DATA           name of datafile to be read in
0028 %         LEG            leg of cloudsat data
0029 %         USE_LATS       latitude midpoints
0030 %         DLAT           latitude (scalar)
0031 %         LAT_RES        latitude resolution (scalar)
0032 %         LON_GRID       longitude grid
0033 %         C.P_LIMS       pressure limits
0034 %         DIM            desired dimensionality of data
0035 %         DIMADD         (see asg_dimadd)
0036 %         PROPS
0037 %                    a psd and ssp defintiion structure
0038 %                    C.PROPS(1) corresponds to the radar frequency
0039 %                    C.PROPS(2) corresponds to the instrument frequencys
0040 %                    with fields
0041 %             radar_back    single scattering properties data
0042 %             PSD    particle size distribution, only 'MH97' and 'Water' works
0043 %             method psd method, only 'gauss-laguerre' works
0044 %             x      particle sizes (see help gauss_laguerre)
0045 %             w      weights
0046 %             x_norm normalistaion factor
0047 %             shape  particle shape, only 'sphere' works
0048 %
0049 %
0050 %        Examle usage:
0051 %        visit asg_demo to see an example how to use this function
0052 
0053 
0054 %  2009-05-18 created by Bengt Rydberg
0055 
0056 function [H,C] = asg_create_gfs(Gcs,grids_cs,C,workfolder)
0057   
0058   
0059 %- Load CloudSat data
0060 %
0061 dim=C.DIM;
0062 C.AUX       = {'DEM_elevation','Longitude','TAI_start','Profile_time'};
0063 if dim==3
0064  lat_limits  = [ C.USE_LATS(1)-C.DLAT/2 C.USE_LATS(end)+C.DLAT/2 ];
0065 end
0066 if dim==1
0067   lat_limits = [C.USE_LATS(1)-C.LAT_RES C.USE_LATS(end)+C.LAT_RES]; 
0068 end
0069 
0070 %make sure we read in also just outside lat_limits
0071 lat_out     = [ lat_limits(1)-C.LAT_RES-0.1 lat_limits(end)+C.LAT_RES+0.1 ];
0072 
0073 [PATHSTR,name1,EXT,VERSN] = fileparts(C.DATA);
0074 if strcmp('.zip',EXT)
0075    unzip(C.DATA,workfolder)
0076    C.DATA=fullfile(workfolder,name1);
0077 end
0078 
0079 [Gdbz] = asg_cloudsat_dBZe([],C.DATA,C.LEG,lat_out,[],C.AUX,[]);
0080 
0081 %find out what longitudes and doys USE_LATS correspond to
0082 for li=1:length(C.USE_LATS)
0083  [lat_m,lat_ind(li)]=min( abs(Gdbz(2).GRID1-C.USE_LATS(li)));
0084 end
0085 lon_ind=find(strcmp({Gdbz.NAME},'Longitude'));
0086 C.USE_LONS=Gdbz(lon_ind).DATA(lat_ind);
0087 mjd_ind=find(strcmp({Gdbz.NAME},'mjd'));
0088 doy=mjd2doy(Gdbz(mjd_ind).DATA(lat_ind));
0089 C.MJD=Gdbz(mjd_ind).DATA(lat_ind);
0090 
0091 %- reduce the latitude resolution
0092 %  bin the data on a desired grid
0093 %
0094 grids{1} = Gdbz(1).GRID1;
0095 grids{1} = Gdbz(1).GRID1( find( Gdbz(1).GRID1>0 & Gdbz(1).GRID1<24e3 ) );
0096 
0097 if dim==3
0098    grids{2} = lat_limits(1) : C.LAT_RES : lat_limits(end);
0099 else
0100    grids{2}= sort([C.USE_LATS C.USE_LATS+C.LAT_RES C.USE_LATS-C.LAT_RES]);
0101 end
0102 
0103 
0104 Gdbz(1).DATA=10.^(Gdbz(1).DATA/10);
0105 [Gdbz]=asg_bin(Gdbz(1),grids,[1 2]);
0106 
0107 
0108 
0109 if dim==1
0110  for i=1:length(C.USE_LATS)
0111    l_ind(i)=find(Gdbz.GRID2==C.USE_LATS(i));
0112  end
0113  Gdbz.DATA=Gdbz.DATA(:,l_ind);
0114  Gdbz.GRID2=Gdbz.GRID2(l_ind);
0115 end
0116 
0117 %switch back to dBZ
0118 Gdbz.DATA=10*log10(Gdbz.DATA);
0119 Gdbz.DATA(find(Gdbz.DATA<=-50))=-50;
0120 
0121 
0122 % Switch to pressure grid.
0123 Gdbz.GRID1= z2p_cira86(Gdbz.GRID1,mean(lat_limits),doy(1) );
0124 Gdbz.GRID1_NAME='Pressure';
0125 Gdbz.GRID1_UNIT='Pa';
0126 
0127 
0128 %- Set Gdbz.GRID1
0129 %
0130 % Get out part of Gdbz that is inside the selected
0131 %C.P_LIMS and latitude limits .
0132 %
0133 grids_dbz{1} = Gdbz.GRID1;
0134 if dim==3
0135 
0136  limits=[fliplr(C.P_LIMS);lat_limits];
0137  [Gdbz]=asg_crop(Gdbz,limits,[1 2]);
0138 
0139 else
0140  
0141  limits=[fliplr(C.P_LIMS)];
0142  [Gdbz]=asg_crop(Gdbz,limits,[1]);
0143  if length(Gdbz.GRID2)<2
0144    Gdbz.DATA=Gdbz.DATA';
0145  end
0146 
0147 end
0148  
0149 %
0150 
0151 %add the psd and radar backscatter properties of the assumed particles
0152 %
0153 Gdbz.PROPS=C.PROPS;
0154 
0155 % Set DIMADD based on fields of C.
0156 %
0157 if dim==3
0158  grids_dbz{3} = C.LON_GRID;
0159  Gdbz.DIMADD=C.DIMADD;
0160 else
0161  Gdbz.DIMADD=[];
0162 end
0163 %
0164 
0165 %- Prepare data for each selected latitude
0166 %
0167 for ip = 1 : length(C.USE_LATS)
0168  
0169   Gcs1 = Gcs;
0170   %- Extract CloudSat data
0171   %
0172   if dim==3
0173 
0174    lat  = [C.USE_LATS(ip)-C.DLAT/2 C.USE_LATS(ip)+C.DLAT/2];
0175    limits=[fliplr(C.P_LIMS);lat];
0176    [G(ip)]=asg_crop(Gdbz,limits,[1 2]);
0177    grids_dbz{2} = G(ip).GRID2;
0178 
0179   else
0180 
0181    G(ip)=Gdbz;
0182    G(ip).DATA=Gdbz.DATA(:,ip);
0183    G(ip).GRID2= Gdbz.GRID2(ip);
0184 
0185   end 
0186 
0187   
0188   Q = qarts;
0189   Q.P_GRID=grids_cs{1};
0190 
0191   if dim==3
0192 
0193    lim1=round(C.USE_LATS(ip)-C.DLAT/2);
0194    lim2=round(C.USE_LATS(ip)+C.DLAT/2);
0195    if ip==1
0196       lat_grid=grids_cs{2};
0197    end
0198    %the lat grid around the cloudbox
0199    add_lat_grid=[ lim1-3 lim1-1 lim1-0.01 lim1 ...
0200                  C.USE_LATS(ip)-1 C.USE_LATS(ip) C.USE_LATS(ip)+1 ...
0201                  lim2 lim2+0.01 lim2+1 lim2+3];
0202    grids_cs{2}=unique(sort([lat_grid,add_lat_grid])); 
0203 
0204    Q.LAT_GRID=grids_cs{2};
0205    Q.LON_GRID=grids_cs{3}; 
0206    Q.ATMOSPHERE_DIM = 3;
0207    
0208   else
0209 
0210    Q.ATMOSPHERE_DIM = 1;
0211    Q.LAT_GRID=C.USE_LATS(ip);
0212    
0213   end 
0214  
0215   %- Generate a single scene
0216   %
0217   %see sub-function
0218 
0219   [H{ip}] = asg_gen_state( Gcs1, grids_cs, G(ip), ...
0220                           grids_dbz, Q,dim);
0221   
0222   
0223 end
0224 
0225 clear G
0226 G=H;
0227 
0228 
0229 function [G] = ...
0230         asg_gen_state(Gcs,grids_cs,Gdbz,grids_dbz,Q,dim)
0231 
0232 
0233 
0234  %- Create clear sky fields
0235  %
0236  if dim==3
0237     Gcs = asg_dimadd( Gcs, Q );
0238  end
0239 
0240  %regrid data
0241 
0242  Gcs = asg_regrid( Gcs, Q );
0243  Gcs = asg_rndmz( Gcs );
0244  Gcs = asg_hydrostat( Gcs, Q );
0245 
0246 
0247  %- Determine grids for dbz field
0248  %
0249  
0250  Q.P_GRID   = grids_dbz{1};
0251  if dim==3
0252     Q.LAT_GRID = grids_dbz{2};
0253     Q.LON_GRID = grids_dbz{3};
0254  end
0255  %
0256  
0257  %- Create final dbz field(s)
0258  %
0259  
0260 
0261  if dim==3
0262   Gdbz.DATA=10.^(Gdbz.DATA/10);
0263   Gdbz = asg_dimadd( Gdbz, Q );
0264   Gdbz = asg_regrid( Gdbz, Q );
0265   Gdbz.DATA=10*log10(Gdbz.DATA);
0266   %pad the data with "zeros" around the dbz field, to make sure interpolation
0267   %will cause no cloud "enlarging"
0268   Gdbz = asg_zeropad( Gdbz, Q);
0269   Gdbz.DATA([1:2 end-1:end],:,:)=-50;
0270   Gdbz.DATA(:,[1:2 end-1:end],:)=-50;
0271   Gdbz.DATA(:,:,[1:2 end-1:end])=-50;
0272  end
0273 
0274  
0275  %- Make re-gridding of fields so cloud and atmospheric variables match
0276  %
0277  %
0278  
0279  if dim==1
0280     Q.ATMOSPHERE_DIM=1;
0281  end
0282 
0283  Q.P_GRID=250/16e3;
0284  if dim==3
0285    Q.LAT_GRID=0.001;
0286    Q.LON_GRID=0.001;
0287  end
0288 
0289  %dummy fields
0290  Gdbz.RNDMZ=[];
0291  Gdbz.SPCFC=[];
0292  Gdbz.SURFACE=[];
0293  len=length(Gcs);
0294  Gcs(len+1)=Gdbz;
0295   
0296 
0297  
0298  [grid1a] = gf_get_grid( Gcs(1), 1 );
0299  [grid1b] = gf_get_grid( Gdbz, 1 );
0300  grid1=flipud(vec2col(sort(union(grid1a,grid1b))));
0301  grid1 = gridconvert( grid1, false, @log10, true );
0302  grid1 = gridthinning( grid1, Q.P_GRID );
0303  Q.P_GRID = gridconvert( grid1, true, @pow10 );
0304  
0305  if dim==3
0306   [grid1a] = gf_get_grid( Gcs(1), 2 );
0307   [grid1b] = gf_get_grid( Gdbz, 2 );
0308   grid1=vec2col(sort(union(grid1a,grid1b)));
0309   grid1 = gridthinning( grid1, Q.LAT_GRID );
0310   Q.LAT_GRID=grid1;
0311 
0312   [grid1a] = gf_get_grid( Gcs(1), 3 );
0313   [grid1b] = gf_get_grid( Gdbz, 3 );
0314   grid1=vec2col(sort(union(grid1a,grid1b)));
0315   grid1 = gridthinning( grid1, Q.LON_GRID );
0316   Q.LON_GRID=grid1;
0317  end
0318 
0319  ind=strcmp(lower({Gcs.NAME}),'radar_reflectivity');
0320  Gcs(ind).DATA=10.^(Gcs(ind).DATA/10);
0321  Gcs = asg_regrid( Gcs, Q );
0322  Gcs(ind).DATA=10*log10(Gcs(ind).DATA);
0323  Gcs = asg_hydrostat( Gcs, Q );  
0324 
0325  %extract iwc and particle number fields(pnd) fields
0326  Gcs = asg_dbz2pnd( Gcs, Q ,Gcs(end).PROPS);
0327 
0328  %- Modify water vapour to match cloud distribution
0329  %
0330  ind1=find(strcmp('h2o',lower({Gcs.NAME})));
0331  if  0 %Gcs(ind1).SPCFC.MOD_RH
0332 
0333    Gcs=asg_iwc_relhumid(Gcs,Q);
0334 
0335    Gcs = asg_hydrostat( Gcs, Q );
0336 
0337  end
0338 
0339 G=Gcs;
0340 
0341 return

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