wiki:FpCtbtoOverview

Version 10 (modified by dearn, 9 years ago) (diff)

--

CTBTO Project Wiki Page

Update: 2015-03-18

This project is carried out in collaboration by the following partners:

The aim of the project is to enhance the FLEXPART modeling system for CTBTO operations, and to contribute the enhancements into the mainstream FLEXPART distribution for the benefit of others. This should enable CTBTO to benefit from new science and techniques used in future FLEXPART versions, while insuring that their own enhancements remain part of the distribution.

The scope of the project includes the following functionality areas:

  1. Allow GRIB input data to be pre-processed before FLEXPART runs
  2. Allow reading of both NCEP and ECMWF GRIB files using the same executable
  3. Allow dynamic allocation of simulation grid and particles
  4. Add configuration option for outputting same units in backward runs as in forward runs
  5. Support reading ECMWF GRIB files in various sub-formats
  6. Add Makefiles for ifort

Work is currently ongoing on the second functionality-area above. More specifically, this entails:

  • Unifying the GRIB reading routines in FLEXPART so that one executable may use either NCEP or ECMWF met files.
  • Building a simple, initial regression testing system by which modifications to FLEXPART can be validated for consistency in output and performance
  • Setting up a private developmental repository and ticketing system for the project in a way that developmental releases (a CTBTO-developmental branch) can be made accessible to the mainstream flexpart.eu community for evaluation and integration into the main branch
  • Planning and providing a presentation on the project at the FLEXPART Developers meeting to be held in Vienna during EGU week.

Updates on the project will be described in this public wikipage and will also be introduced in the ticketing system for easy follow up from FLEXPART developers and the community.

Update: 2015-06-30

To perform the unification in the GRIB reading routines in FLEXPART to achieve a single executable, the following steps have taken place:

  1. Documentation and test cases have been added into the source code within their respective directories. The testing environment, FLEXtest, is constituted of a set of simple modules that can be independently called but that are meant to be in a high level structure. http://flexpart.eu:8088/svn/CTBTOdevel/trunk/flexpart-testing/doc/ for readily available mark down documentation. An example of the functioning and outcome of the testing is added below.
  1. The FLEXPART version in the trunk (as of date 16 February 2015) has been used to program the unification of routines. The code has modified as follows:
  • renaming of duplicated routines so that both versions could be used in a single build
  • adding a routine to detect the format to be used (detectformat.f90)
  1. Added to the distribution a fairly simple and documented tested python-based environment documented python-based testing environment.

The code, documentation and environment may be obtained in http://flexpart.eu:8088/svn/CTBTOdevel/trunk/


EXAMPLE ON USAGE OF THE TEST ENVIONMENT '

The file FlexpartErrors/test/test_FlexpartErrors.py contains a variety of test cases which can serve as examples on how to use this class. However, the following is intended to be a more user-friendly introduction. Be sure fist to set your $PYTHONPATH correctly. If you have put the testing environment in directory /home/mememe/flexpart-testing you would set your $PYTHONPATH in the bash shell as follows

export PYTHONPATH=/home/mememe/flexpart-testing:$PYTHONPATH

The data used in these examples come from a tiny, nested outgrid of 4 rows, 3 columns, 5 levels, 2 species, 6 releases and 1 age class in the mother domain. In the nest, there are 10 rows and 5 columns of grid points. The control model comes from a simulation using ECMWF met data, and the test run comes from the same simulation using NCEP met data. There are three timestamps:

['20140919010000', '20140919020000', '20140919030000']

The test data here is available in the unittest_data directory of the code distribution. As a first example, let's get a grid that represents the difference (test minus control) between the two datasets, under various conditions. First, we create the FlexpartOutput objects as follows. You will need to change the value of DATA_DIR to the location of your distribution.

import flexread.FlexpartOutput as flexout
import FlexpartErrors as flexerr

DATA_DIR='/u1/uaf/morton/flexpart-testing/unittest_data'

CONTROL_OUTPUT_DIR = DATA_DIR + '/flexout_forward_tiny_complex_withnest'
TEST_OUTPUT_DIR = DATA_DIR + '/ncepflexout_forward_tiny_complex_withnest'

# Create the two objects used for containing the FLEXPART output
# By default, this will point to the mother nest of the domains.
control_output = flexout.FlexpartOutput(output_dir=CONTROL_OUTPUT_DIR)
test_output = flexout.FlexpartOutput(output_dir=TEST_OUTPUT_DIR)

# Create the error object
error_obj = flexerr.FlexpartErrors(control=control_output,
                                   test=test_output)

# Get the diff_grid from default parameters
diff_grid = error_obj.get_diff_grid()

print diff_grid

In the above, we create control_output and test_output objects, then create an error_obj object by initialising with these two objects. With no arguments to the diff_grid method, we use default values, giving us a horizontal slice at level 1, for a roughly middle timestamp

[[ 0.          0.          0.        ]
 [ 0.         -0.09765793  0.        ]
 [ 0.          0.18723705  0.        ]
 [ 0.          0.          0.        ]]

Exercising a little more specificity, in the next case we get the diff grid from last timestamp, 20140919030000, level 2, species 1, release 5, as follows

diff_grid = error_obj.get_diff_grid(timestamp='20140919030000',
                                    level=2,
                                    species=1,
                                    release=5
                                     )

print diff_grid
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [ -2.90811108e-02  -6.62608682e-05   0.00000000e+00]
 [ -1.18361338e-01  -8.55133458e-05   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00]]

Up to now, we have used the get_diff_grid() method to illustrate how to access different components of the output files. In most cases, the user will not even use get_diff_grid() unless he or she wants more detail. This FlexpartErrors class is used primarly to calculate the error between two grids. The available methods are described in http://flexpart.eu:8088/svn/CTBTOdevel/trunk/flexpart-testing/doc

In the following example we go through the calculation of errors on different portions of the output grids for the nest

import flexread.FlexpartOutput as flexout
import FlexpartErrors as flexerr

DATA_DIR='/u1/uaf/morton/flexpart-testing/unittest_data'

CONTROL_OUTPUT_DIR = DATA_DIR + '/flexout_forward_tiny_complex_withnest'
TEST_OUTPUT_DIR = DATA_DIR + '/ncepflexout_forward_tiny_complex_withnest'

# Create the two objects used for containing the FLEXPART output
# By default, this will point to the mother nest of the domains.
control_output = flexout.FlexpartOutput(output_dir=CONTROL_OUTPUT_DIR,
                                        nest=True)
test_output = flexout.FlexpartOutput(output_dir=TEST_OUTPUT_DIR,
                                     nest=True)

# Create the error object
error_obj = flexerr.FlexpartErrors(control=control_output,
                                   test=test_output)

# Calculate mean bias for the default time and horizontal slice of the next
bias = error_obj.mean_bias()
print 'Mean bias for default horizontal slice: ' + str(bias)

# Calculate the max error (max magnitude, whether positive or negative) of
# the entire time series and entire volume for species 2, release 1
maxerr = error_obj.max_error(species=2, release=1,
                             volume=True, timeseries=True)
print 'Max error for full timeseries and volume, species 2, release 1: ' + \
        str(maxerr)

# Calculate the RMSE of last timestep, level 2, for species 1, release 2
rmse = error_obj.rmse(timestamp='20140919030000', level=2,
                        species=1, release=2)
print 'RMSE for level 2, species 1, release 2, last timestamp: ' + \
        str(rmse)
Mean bias for default horizontal slice: 0.0449703330779
Max error for full timeseries and volume, species 2, release 1: 3.79051153362
RMSE for level 2, species 1, release 2, last timestamp: 0.146961582495

FlexpartErrorsTemporaryDocumentation?

hosted by ZAMG