Settings file is written in YAML format and that brings many benefits: easy to write, clear to read, self-explanatory names, flexible and powerful. The concepts are organized in blocks that are translated to YAML thanks to the fix indentation (we suggest 4 white-spaces). So, for example, the first level defines the different modules:

                    

input:
    # Here settings related with input module
    segment: # It defines a description of the input segment
        x: 2 # It defines the maximum size of x dimension of the segment to 2.
output:
    # Settings linked with the output
    # ...
retrieval:
    # Definition of inversion strategy
    # ...
...
            

When grasp is called, the first step is to read settings and the second is to prepare the environment for the settings defined in the main structures. If the settings are not valid, an explicit error message will be printed. Please read the first line of it carefully to understand the error.

The settings parameters can also be defined by the command line. In this case, after the first argument (settings file name), extra settings parameters can be defined with the syntax key=value where key is the parameter name in "dot syntax". For example, in GRASP, it is equivalent to be called with the argument input.driver=sdata or to have defined in the settings file the following content:

input:
    driver: sdata
                

It is important to define clearly how the relative paths have to be defined in the settings files. Relative paths are always relative to the file that defines it. In the next section, the reader will learn about how to include other settings file inside of a settings file, but this rule will stay valid: Relative paths are defined from settings file that define it. In the case of usage of command line, relative paths are relative to the current working directory. In the case of absolute paths, all this complexity disappears but the results are less portables.

                    

output:
    segment:
        function: hdf
        stream: "GRASP_Banizoumbou_20080101_20080331_2x2+3286+1376.hdf"
    tile:
        function: [ ascii, hdf ] 
        stream: [ "GRASP_Banizoumbou_20080101_20080331.txt",
                "GRASP_Banizoumbou_20080101_20080331.hdf" ]
                

                    

output:
    segment:
        function: hdf
        stream: "GRASP_Banizoumbou_{tile_from(%Y%m%d)}_{tile_to(%Y%m%d)}_{segment_nx}x{segment_ny}+{segment_corner_column(4)}+{segment_corner_row(4)}.hdf"
    tile:
        function: [ ascii, hdf ] 
        stream: [ "GRASP_Banizoumbou_{tile_from(%Y%m%d)}_{tile_to(%Y%m%d)}.txt",
                "GRASP_Banizoumbou_{tile_from(%Y%m%d)}_{tile_to(%Y%m%d)}.hdf" ]
                

The ensemble of characteristics available to be retrieved or simulated by GRASP is open to user selection in the settings file inside retrieval.constraints.characteristic section. An example of the general structure of them is showed below:

                
retreival:
    constraints:
        characteristic[1]:
            type: characteristic_name
            retrieved: true
            mode[1]:
                initial_guess:
                    value:                          [0.0, 0.0, 0.0, 0.0, ...]
                    min:                            [0.0, 0.0, 0.0, 0.0, ...]
                    max:                            [0.0, 0.0, 0.0, 0.0, ...]
                    index_of_wavelength_involved:   [0.0, 0.0, 0.0, 0.0, ...]
                single_pixel:
                    smoothness_constraints:
                        difference_order: 0.0
                        lagrange_multiplier: 0.0
                multi_pixel:
                    smoothness_constraints:
                        derivative_order_of_X_variability:    0.0
                        lagrange_multiplier_of_X_variability: 0.0
                        derivative_order_of_Y_variability:    0.0
                        lagrange_multiplier_of_Y_variability: 0.0
                        derivative_order_of_T_variability:    0.0
                        lagrange_multiplier_of_T_variability: 0.0
                    
            mode[2]:
            
                ...
            
            mode[3]:
            
                ...
                
            ...
                
                

The number assigned at each characteristic has no relevance at all while coherence is maintained. The retrieved field can be set to 'true' if the corresponding characteristic is going to be a retrieved parameter in the inversion; or to 'false' if it is just a fixed value for the forward simulation. The number of modes included in each characteristic depends on the nature of each one. For characteristics representing optical or mycrophisical aerosol properties (size distribution or refractive index for example), the number of modes corresponds to the number of aerosol modes selected for the retrieval/simulation (tipically one or two if fine/coarse distinction is made). For characteristics representing surface properties three modes will be needed to describe the associated model; except for polarization that only one is nedeed. The fields present in initial_guess part contain one element for each wavelength following the structure provided in the SDATA file in the case of optical wavelength dependent characteristics; one element for each bin for size distribution related characteristics; and one element for other mycrophysical magnitudes or non-wavelength dependent optical characteristics. In the former case, index_of_wavelength_involved should be filled with zeros for all the corresponding bins. The rest of the elements included in single_pixel or multi_pixel parts are always formed by one single element.

The list of available characteristics can be found below:

4.1.2.1. Initial guess through the algorithm

Understanding GRASP inversion procedure means understanding how the algorithm is starting from an initial guess and obtains a results array. This is an iterative procedure explained in the literature and introduced in Section 2.3, “GRASP Scientific Core algorithm”. The purpose of this section is to explain how initial guess is represented inside the code as an array which evolve in each iteration until geting the result array. The following diagram shows how initial guess is read from settings file and translated into an internal array in the code. This detail could look very technical and related with the development but understanding of some internal concepts of the code helps to understand how it works. The following diagram shows this transformation:

Figure 4.1. Translation of settings file into initial guess array

Once the initial guess is loaded, it is set as first array of characteristics to be retrieved. Retrieval process iterates over it until it gets the results. The results of GRASP retrieval is the array with the same shape as the initial guess, but containing the results of retrieval of these parameters to match them with SDATA file. Then, the GRASP forward model is called using this results array, GRASP obtains the rest of the derived results it provides. So, we can talk about two kinds of results: basic results and derived products. The following diagram shows this process:

Figure 4.2. Evolution of retrieved characteristics during GRASP processing

To know all the details about the products obtained by GRASP, please see Section 4.3.1, “The list of GRASP output parameters”.

The SDATA (sensor data) format is the original input data format of the GRASP code. It is a simple text format designed by the science team at the early stages of the development of the scientific code and it's the easiest way to create test data.

In the context of the GRASP project, the SDATA format has a number of pros:

It has also a number of cons:

Whatever the number and seriousness of the cons, one of the design objectives is to keep the code simple and flexible allowing the scientific community to play with the code. It does not make sense to implement a driver for a single user who wants to do some tests with GRASP. That's why this easy format is maintained by the developer team.

In the following description, the elements in fixed-width font are the snippets of content. The numeric values in these snippets (e.g. in 2 2 2 : NX NY NT) are given only as examples.

        

SDATA version 2.0
2   2   2  : NX NY NT

4       2008-01-04T13:15:00Z    70000.0   0   0   : NPIXELS  TIMESTAMP  HOBS_km  NSURF  IFGAS
1	    1		1	 3286	     1377   2.599	13.528       252.0	 100.0      6  0.443   0.490   0.565 ...
2	    1		1	 3287	     1377   2.657	13.528       242.0	 100.0      6  0.443   0.490   0.565 ...
1	    2		1	 3286	     1376   2.601	13.583       241.0	 100.0      6  0.443   0.490   0.565 ...
2	    2		1	 3287	     1376   2.658	13.583       239.0	 100.0      6  0.443   0.490   0.565 ...


4       2008-01-06T13:02:41Z    70000.0   0   0   : NPIXELS  TIMESTAMP  HOBS_km  NSURF  IFGAS
1	    1		1	 3286	     1377   2.599	13.528       252.0	 100.0      6  0.443   0.490   0.565 ...
2	    1		1	 3287	     1377   2.657	13.528       242.0	 100.0      6  0.443   0.490   0.565 ...
1	    2		1	 3286	     1376   2.601	13.583       241.0	 100.0      6  0.443   0.490   0.565 ...
2	    2		1	 3287	     1376   2.658	13.583       239.0	 100.0      6  0.443   0.490   0.565 ...
                

A SDATA file has a simple structure:

2 2 2 : NX NY NT

A CELL is a set of neighbouring pixels, that form the base of a SEGMENT. Each CELL has a HEADER and a number of PIXELs, supposedly acquired at the same time.

The CELL HEADER contains:

4 2008-01-04T13:15:00Z 70000.0 0 0 : NPIXELS TIMESTAMP HOBS NSURF IFGAS

Each line of data after the CELL HEADER represents exactly one pixel, with all its fields. The Table 4.5, “The PIXEL structure” describes the order and type of these fields. For the types, the Fortran notation is used: array types are described with the dimensions of arrays between parentheses, and the ordering is such that the first index increases faster. Indices start from 1, not from 0 like in C. For instance, when one reads real(nwl) for wavelengths, that means that one has to read a list of nwl real values that represent wavelengths.

The field pixel[ipix].meas[nwl].meas_type[nip] of pixel structure is a special code which defines the type of measure. The following table describes the valid codes and their interpretation:

This chapter main goal is to describe how the angles should be defined to be used inside of GRASP code. The universal spirit of GRASP, where many different instruments coexist (from satellite to ground based measurements), creates challenges to define a homogeneous way of defining the angles, keeping a unique geometry. As it is shown in the figure Figure 4.4, “Definition of GRASP geometry” GRASP angles are defined to be considered as "normal" for satellite reference.

Figure 4.4. Definition of GRASP geometry

Since GRASP angles are defined using a satellite reference, it provokes some problems to define what we could call as an intuitive "ground based" reference system. That is why we are going to put special emphasis on the definition of the angles to these less intuitive applications. The intuitive reference for "ground based" measurements, in spherical geometry, is given as follows:

where the sub index "gb" makes reference to "ground based".

The conversion to the GRASP geometry is done as follows:

Here we propose some examples for better understanding of the process. Before defining the measurement angles introduced in the code, both "intuitive" and "GRASP", we need first to consider the instrument viewing angle for each scenario (θ_v, ϕ_v). The following figure and table will provide some examples of angles defined for the ground based applications.

Figure 4.5. Ground based angles definition example

Since sunphptometers are widely used with GRASP, the following table provides information specifically about these instruments, considering the instrument viewing angle for each scenario (θ_v, ϕ_v). They can be understood as the "movements of the motors". The process will be as follow:

instrument viewing angle -> angle in (intuitive) ground based -> angle in GRASP

In the case of nephelometer data angle definition, it is a bit different since only θ angle has to be defined, the rest of the angles will be ignored. To provide nephelometer data (phase matrix), the conversion to the GRASP geometry is done as follows:

Sunphotometers are widely used with GRASP. They take measurements of sky radiance and direct sun. Many inversion strategies can be used to retrieve sunphotometer data, but this section will explain how to define input data. The information that runs inside of GRASP has to be pre-processed in order to screen clouds, calibrate and normalize the data.

Following Section 4.2.1, “The SDATA format”, the direct sun measurements can be described as AOD or TOD defined in Table 4.6, “TYPES OF MEASUREMENTS” as measurements of type 11 or 12. At this point it is needed to take into account that if "ifgas" field is defined as 1 in the case of AOD, no gaseous absorption optical depth will be accounted, but in the case of TOD they will be subtracted. The gases also affect the radiance measurements but in lower magnitude. In the case of TOD + radiances with ifgas=1, the same model will be applied to all measurements. If AOD is used, some (minor) incongruences could come from the use of different models to calculate gases for AOD and for radiances.

Radiance measurements are defined with the constant MEAS_TYPE_I(41). Polarized measurements can be defined as Q,U (42, 43) or as polarization rate (44). It is also important to check how polarized data is going to be manipulated in the retrieval code based on inversion strategy defined in the settings file.

Note that all processing will be considering range corrected profile for one wavelength. The procedure for other profiles from different wavelengths are exactly the same. Range corrected profile implies that at least the background noise was subtracted and altitude corrections were applied to the raw signal, but if it's possible to consider all other corrections (electrical noise and overlap correction, dead time correction, gluing analog and photon-counting signals and all others that your system may have), you should apply them.

Step 1. Background noise subtraction and range correction.

Let B' be the estimation of the background noise. Usually B' is estimated as P(Z_B), accumulated and averaged around selected altitude Z_BZ, much higher than the maximum altitude of lidar extraction in step 2 (Z_max). For example with maximum altitude Z_max=15km, it is averaged around 30 km and accumulated for the whole period of lidar observation. The noise and range corrected signal will be:

S(Z_i) = (P(Z_i)-B')*Z_i²

Step 2. Altitude range selection and signal cropping.

The minimum Z_min and maximum Z_max altitudes are selected and the signal is cropped, so Z_min<Z_i<Z_max. The minimum altitude should be selected as low as possible, preferably in the region where the overlap correction could be correctly applied. The Z_max should be selected from the following considerations: maximum altitude where the noise levels of lidar measurement are acceptable and the amount of atmospheric aerosols is still noticeable.

In case if the lidar used is inclined The sounding zenith angle should be provided for the corresponding wavelenghts. The angle value should be placed in the position of the SZA corresponding to this wavelength. All SZA's for all vertical profile measurements should be the same. Vertically pointed lidar should have this value set to 0. Note: Keep in mind that retrieved aerosol vertical profiles will be retrieved for vertically projected altitudes, i.e. if your lidar is inclined, the maximum altitude of the profile will be equivalent to Z_max*cos(SZA).

Step 3. Backscatter of molecular profile

Backscatter vertical profile is calculated inside of GRASP based on standard atmosphere model for each lidar wavelength[1]:

β_mol(Z_i,λ) = N * σ(λ) * (P(Z_i)/P_SA) * (T_SA/T(Z_i))

        where:

Step 4. Reducing the number of points in profiles

GRASP/GARRLiC can use arbitrary altitude/range scale. However, to keep number of retrieved parameters reasonable and to fight higher noise contamination of lidar signals at higher altitudes, it is recommended to use a logarithmical altitude/range scale with N_Z points to represent aerosol profiles in the atmosphere. For that, we have to present all vectors (altitude/range vector and profiles of lidar signals) in logarithmically equidistant manner.

1. Move to logarithmic scale and find altitude/range step:

logarithmic scale: Z_i ^lg = lg(Z_i)

step in log scale: ΔZ = (Z_max ^lg - Z_min^lg)/N_Z

logarithmic altitude ranges (from h_k to h_k+1 ) for averaging data in logarithmically equidistant manner, k = 1 .. N_z: h_k = Z₀ ^lg + (k - 1) * ΔZ

2. Average the data profiles

A_k = (Σ_j=1ⁿ A_j(h_k, h_k+1) ) / n

        where:

After such procedure, the number of points in three main vectors reduced to N_Z points in logarithmically equidistant manner. These values should be placed in places corresponding to the zenith wieving (vza or θ_v) angles for the vavelenths corresponding to lidar or vertical measurements.

Hint: the altitude/range vectors for measurements at all wavelengths have to be the same.

Step 5. Profile normalization

Values of lidar signals (except for volume depolarization profiles) vary from instrument to instrument, from detector to detector, that is why GRASP/GARRLiC requires normalized lidar signal. For consistency with the molecular optical depth, the profile of the molecular backscatter inside the code has to be normalized as well. Normalized lidar and backscatter profiles:

A'_k = A_k / ∫_Zmin^Zmax A_k dZ

        where A represents profile of lidar signal or molecular backscatter.

Caution: integration have to be done using meters in altitudes.

At the end, for each wavelength you have to have normalized lidar profiles and altitude vector. Volume depolarization profiles don't need to be normalized, the only requirement for such observations is to be presented in the percentage range i.e. [1.0e-9, 100]

References: Anthony Bucholtz, Rayleigh-scattering calculations for the terrestrial atmosphere, Optical Society of America, 1995.

The output from GRASP is quite complex and strongly dependent on the settings used. As it was discussed in Section 4.1.2, “Retrieved characteristics”, there are two kinds of output products: direct and derived. Direct results are the values directly inverted in the initial guess array, then after obtaining them, forward model is called one more time to obtain the derived products. The list of direct and derived products obtained depends on the data and the inversion strategy selected (defined in the settings file). This chapter contains a complete list of products that can be obtained with GRASP, but it does not mean that all of them can be obtained with all inversion strategies.

In the internal GRASP output structures there are some products that are repeated between direct products structure (array with the same shape as initial guess) and derived products. The reason is that for some applications they are direct information, for others they are derived results. This depends on input data and inversion strategy, defined in the settings file.

Before defining the output, we are going to define the size of some arrays. This information is needed to understand the list of the products. For example, when a product such as AOD is defined as many times as wavelengths, it is because the output will be wavelength dependent (a value for each wavelength):

The following list includes all products that can be obtained using GRASP:

As it was explained in Section 4.1.2, “Retrieved characteristics” section, the retrieval algorithm works iteratively over an array of parameters (in its first definition it is called the initial guess), until it represents the best solution. This solution array has same shape as the initial guess (the same parameters and defined in the same position). Once it is obtained, a final call of the forward model with the resulting array provides a complete list of output products.

For some applications, it can be useful to use the forward model without inverting any data. It can be done easily in GRASP with the use of the setting parameter retrieval.convergence.stop_before_performing_retrieval=true. When no retrieval is performed, just one call of forward model is performed. If in the initial guess array the user has set an aerosol model, it will be used inside of the forward model obtaining therefore an entire output structure, with information in all fields.

This procedure can be used also to reprocess some data. If output parameters of a retrieval are stored, then they can be set as initial guess and then, running GRASP with the same settings, except for retrieval.convergence.stop_before_performing_retrieval=true the entire output can be obtained again. This procedure can be used to reprocess data with many objectives such as saving storage space (just save the output array of grasp and reprocess to obtain the rest, if it is needed) or obtaining extra products in the future.

The previous procedure can also help to simulate the input data. It is useful because a valid geometry is necessary in the input data to set an SDATA file. In this case, an aerosol model is set as an initial guess and the code works just for the forward run. Then, by using retrieval.debug.simulated_sdata_file parameter, the user can set a path to the simulated files. A SDATA file will be dump to that path where the geometry is the same and the measurements are filled with the output results of the forward run. Then, this SDATA file can be used to self-consistency tests, where synthetic data is retrieved.