Aerosol Optics MERRA

[mjiacono]

Pull request to merge components for deriving aerosol optical properties based on MERRA aerosol optics data.

[RobertPincus]

Hi @mjiacono Thanks for this. Two requests.

1. Can you please change the target of the pull request to the new feature-merra-aerosol-optics branch? I will make some changes (adding CI, putting reference files outside the repo) before merging to develop.
2. Can you please add the ability to automatically run the tests and verify the results? That means adding a tests target to your Makefile that runs the code and a compare-to-reference.py file that exits with an error code if the results produced by the code don't match the input files.

Thanks - Robert

[mjiacono]

1. Target changed to the new branch
2. I'll work on this

[RobertPincus]

@mjiacono How are users to know what integers map to what aerosol type? Is this documented somewhere? In any case we should provide some "introspection" i.e. a way to use variables names (MERRA_BLACK_CARBON) instead of integers.

If we were free to make our own mapping it could open up some efficiencies, like grouping together all the tables where the properties depend on relative humidity.

[RobertPincus]

@mjiacono Thanks for including variable bnd_limits_wavenumber in the netcdf files. Still shy from the recent cloud physics update, can you confirm a) that the SW data are in strictly increasing order of wavenumber and b) that the bnd_limits_wavenumber accurately reflects the data in the fields?

[mjiacono]

@RobertPincus The aerosol types are defined in the ty_aerosol_optics class. I added comments to identify the aero_type integer values in the aerosol_optics function. We can consider converting aero_types from integers to character names, but wouldn’t it be less efficient (in terms of memory) to store character string variable names over an (ncol,nlay) array rather than integers?

[mjiacono]

@RobertPincus The SW aerosol optics coefficient tables are in the proper increasing wavenumber order for RRTMGP in agreement with the values in bnd_limits_wavenumber.

[RobertPincus]

@RobertPincus The aerosol types are defined in the ty_aerosol_optics class. I added comments to identify the aero_type integer values in the aerosol_optics function. We can consider converting aero_types from integers to character names, but wouldn’t it be less efficient (in terms of memory) to store character string variable names over an (ncol,nlay) array rather than integers?

@mjiacono I can see two opportunities here. For clarity, what if you move the private definitions of the mapping between integer and aerosol type to somewhere users can see them? I.e. in the module header, not inside the type definition:

    ! MERRA aerosol type 
    integer, parameter, public :: merra_dust_aero = 1

(I think you should add "merra" to be explicit and avoid name clashes.

We could be fancier and make some data structure that users can inquire the number of aerosol types and their names, but maybe that's not gilding the lily.

With this in place we can think about combining some of the 7x3 tables.

[RobertPincus]

@mjiacono Note that the "continuous integration in a box" is failing because the software containers are now out of sync with the CI; that's my problem not yours.

[mjiacono]

@RobertPincus The aerosol types are defined in the ty_aerosol_optics class. I added comments to identify the aero_type integer values in the aerosol_optics function. We can consider converting aero_types from integers to character names, but wouldn’t it be less efficient (in terms of memory) to store character string variable names over an (ncol,nlay) array rather than integers?

@mjiacono I can see two opportunities here. For clarity, what if you move the private definitions of the mapping between integer and aerosol type to somewhere users can see them? I.e. in the module header, not inside the type definition:

    ! MERRA aerosol type 
    integer, parameter, public :: merra_dust_aero = 1

(I think you should add "merra" to be explicit and avoid name clashes.

We could be fancier and make some data structure that users can inquire the number of aerosol types and their names, but maybe that's not gilding the lily.

With this in place we can think about combining some of the 7x3 tables.

@RobertPincus MERRA aerosol types renamed and declarations moved before the ty_aerosol_optics class within the module header.

[RobertPincus]

@mjiacono There are three tables per aerosol type, one each for extinction, single-scattering albedo, and asymmetry parameter. This makes 21 tables total - a lot to keep track of.

Since each property for a given type depends on band, band and size, or band, size, and RH. Better to have a single table for each aerosol type, where the last index is 1:3. Define (private) integer parameters ext_index, ssa_index, and asy_index. Replace references to e.g. tau_salt_table, ssa_salt_table, asy_salt_table with salt_table(:, :, :, ext_index), salt_table(:, :, :,ssa_index), salt_table(:, :, :, asy_index) Or make the first index correspond to the quantity, return the three interpolated values as a size-3 vector, and unpack those into the output.

It would also be worth changing the netCDF file to have 7 tables rather than 21.

[RobertPincus]

@mjiacono You could go ahead and synchronize your repo with the main RTE-RRTMGP repo with the "sync fork" button just above the list of code. This will make your continuous integration stop failing for the wrong reasons.

[mjiacono]

@mjiacono You could go ahead and synchronize your repo with the main RTE-RRTMGP repo with the "sync fork" button just above the list of code. This will make your continuous integration stop failing for the wrong reasons.

@RobertPincus My fork of the rte-rrtmgp repo has been synchronized with the main repo for several days. The "sync fork" button reports: "This branch is not behind the upstream earth-system-radiation:main"

[RobertPincus]

@mjiacono You could go ahead and synchronize your repo with the main RTE-RRTMGP repo with the "sync fork" button just above the list of code. This will make your continuous integration stop failing for the wrong reasons.

@RobertPincus My fork of the rte-rrtmgp repo has been synchronized with the main repo for several days. The "sync fork" button reports: "This branch is not behind the upstream earth-system-radiation:main"

??? main in the RTE repo was updated at 11 am this morning. This is not visible?

[mjiacono]

@mjiacono You could go ahead and synchronize your repo with the main RTE-RRTMGP repo with the "sync fork" button just above the list of code. This will make your continuous integration stop failing for the wrong reasons.

@RobertPincus My fork of the rte-rrtmgp repo has been synchronized with the main repo for several days. The "sync fork" button reports: "This branch is not behind the upstream earth-system-radiation:main"

??? main in the RTE repo was updated at 11 am this morning. This is not visible?

Your change to main became visible to the fork a few minutes ago. It's synced again.

[mjiacono]

@mjiacono There are three tables per aerosol type, one each for extinction, single-scattering albedo, and asymmetry parameter. This makes 21 tables total - a lot to keep track of.

Since each property for a given type depends on band, band and size, or band, size, and RH. Better to have a single table for each aerosol type, where the last index is 1:3. Define (private) integer parameters ext_index, ssa_index, and asy_index. Replace references to e.g. tau_salt_table, ssa_salt_table, asy_salt_table with salt_table(:, :, :, ext_index), salt_table(:, :, :,ssa_index), salt_table(:, :, :, asy_index) Or make the first index correspond to the quantity, return the three interpolated values as a size-3 vector, and unpack those into the output.

It would also be worth changing the netCDF file to have 7 tables rather than 21.

@RobertPincus I completed changes to the code and data files to consolidate the aerosol tables.

[RobertPincus]

@mjiacono The most recent round of changes look very good, I hope also to you. Seems to me there are a countably small number of changes to be made:

1. Using check_extents instead of repeated calls to size
2. Encoding the RH discretization in the file rather than in code. (Revisiting the way the discretization is a good idea but optional.)
3. Encapsulating the interpolation and simplifying a) so as not to produce intermediate values of taussa and taussa*g, and b) to assign the interpolated values directly to the output fields to avoid unneeded copies. I see you've started on this.

Otherwise things seem pretty close, thanks very much.

[mjiacono]

@mjiacono The most recent round of changes look very good, I hope also to you. Seems to me there are a countably small number of changes to be made:

1. Using `check_extents` instead of repeated calls to `size`

2. Encoding the RH discretization in the file rather than in code. (Revisiting the way the discretization is a good idea but optional.)

3. [Encapsulating the interpolation ](https://github.com/earth-system-radiation/rte-rrtmgp/pull/202#discussion_r1069681206)  and simplifying a) so as not to produce intermediate values of tau_ssa and tau_ssa*g, and b) to assign the interpolated values directly to the output fields to avoid unneeded copies. I see you've started on this.

Otherwise things seem pretty close, thanks very much.

@RobertPincus (1) and (2) and (3) are completed.

[RobertPincus]

@mjiacono The most recent round of changes look very good, I hope also to you. Seems to me there are a countably small number of changes to be made:

1. Using check_extents instead of repeated calls to size
2. Encoding the RH discretization in the file rather than in code. (Revisiting the way the discretization is a good idea but optional.)
3. Encapsulating the interpolation and simplifying a) so as not to produce intermediate values of tau_ssa and tau_ssa*g, and b) to assign the interpolated values directly to the output fields to avoid unneeded copies. I see you've started on this.

Otherwise things seem pretty close, thanks very much.

@mjiacono

On point 1, see lines 273 and following.

On point 3, a) there's still lots of repeated code in compute_all_from_table - can you encapsulate this in a subroutine, as suggested? And b) you're still going from tau, ssa, g to tau, taussa, taussa*g and back again unnecessarily?

[mjiacono]

@mjiacono The most recent round of changes look very good, I hope also to you. Seems to me there are a countably small number of changes to be made:

1. Using check_extents instead of repeated calls to size
2. Encoding the RH discretization in the file rather than in code. (Revisiting the way the discretization is a good idea but optional.)
3. Encapsulating the interpolation and simplifying a) so as not to produce intermediate values of tau_ssa and tau_ssa*g, and b) to assign the interpolated values directly to the output fields to avoid unneeded copies. I see you've started on this.

Otherwise things seem pretty close, thanks very much.

@mjiacono

On point 1, see lines 273 and following.

On point 3, a) there's still lots of repeated code in compute_all_from_table - can you encapsulate this in a subroutine, as suggested? And b) you're still going from tau, ssa, g to tau, tau_ssa, tau_ssa*g and back again unnecessarily?

@RobertPincus On point 3, a) of the seven tables, there are three different sizes, so you would need three subroutines to do what you're suggesting. I don't see the advantage in terms of amount of code or clarity; b) it's not clear what you're after here; the intermediate values are used to derive the final optics for either 1-stream or 2-stream cases.

[RobertPincus]

@mjiacono Thanks for the explanation, I see the point about needing tau and taussa as products in order to fill in either _1scl or 2str optical properties. For g, however, this shouldn't be necessary, since g isn't used at all in the _1scl optical properties. Perhaps you can interpolate tau, taussa, and g?

Please also introduce private integer parameters TAU = 1, SSA = 2, G = 3 and replace the hard-coded 1, 2, 3 in the first indexes of the tables used in interp_all_from_table with these parameters.

If you aren't going to implement interp_aero_table_2d and interp_aero_table_3d please remove them (indeed please remove any unused code e.g. at line 133)

[mjiacono]

@mjiacono Thanks for the explanation, I see the point about needing tau and tau_ssa as products in order to fill in either _1scl or 2str optical properties. For g, however, this shouldn't be necessary, since g isn't used at all in the _1scl optical properties. Perhaps you can interpolate tau, tau_ssa, and g?

Please also introduce private integer parameters TAU = 1, SSA = 2, G = 3 and replace the hard-coded 1, 2, 3 in the first indexes of the tables used in interp_all_from_table with these parameters.

If you aren't going to implement interp_aero_table_2d and interp_aero_table_3d please remove them (indeed please remove any unused code e.g. at line 133)

@RobertPincus I addressed these points and added the linear interpolation function. Unless you have other suggestions I'll move on to compiling and testing.

[RobertPincus]

@mjiacono The code itself looks ready to merge, thanks for all your efforts.

Can you add an automated test? Ideally this would be Makefile targets for tests and check as in the examples and test directory. make in your aerosol_optics/test directory will build your code and test harness; make tests will run the tests; make check will compare the values to known good output, which you can put in the repo for now.

I'll ensure that these tests are run every time a change is made to the code.

[RobertPincus]

Note that having new code exercise tests is one of the thing we ask contributors to do.

[mjiacono]

@mjiacono The code itself looks ready to merge, thanks for all your efforts.

Can you add an automated test? Ideally this would be Makefile targets for tests and check as in the examples and test directory. make in your aerosol_optics/test directory will build your code and test harness; make tests will run the tests; make check will compare the values to known good output, which you can put in the repo for now.

I'll ensure that these tests are run every time a change is made to the code.

@RobertPincus I am working on setting up a test for aerosol optics, but I haven't used the newer testing infrastructure before and at the moment it isn't building for me. One problem is that I can't access the auxiliary input files for the RFMIP clear sky test because we can't use ftp anymore. I'll try to download them some other way.

[mjiacono]

@RobertPincus I updated my fork to include modifications to the all-sky test to include aerosol_optics. The /testing directory in my fork is now obsolete and can be disregarded. Note that I left the cloud testing in place, which was set up to run in two out of every three profiles. I set up the aerosols to run in every other profile, so the results reflect all combinations: clear-only, clear-cloud, clear-aerosol, and clear-cloud-aerosol. I did not revise any of the documentation to reflect the aerosol extension code or testing. A plot showing the clear-aerosol and clear-only results is attached here.

[RobertPincus]

@mjiacono Love your approach to testing the aerosols in all combinations of clouds. Thanks also for adding the figure. Could you please git rm your testing directory? I'll move the netcdf file with the updated all-sky results to the continuous integration staging area before merging.

This seems done to me, thanks very much, I'll try to get this merged next week.

[mjiacono]

@mjiacono Love your approach to testing the aerosols in all combinations of clouds. Thanks also for adding the figure. Could you please git rm your testing directory? I'll move the netcdf file with the updated all-sky results to the continuous integration staging area before merging.

This seems done to me, thanks very much, I'll try to get this merged next week.

@RobertPincus Glad you're satisfied with the results. I removed the testing directory. Let me know if you spot anything else.