Fixing CloudSat Precip flag surface elevation bug

[rodrigoguzman-lmd]

There is a bug in the CloudSat Precipitation Flag diagnostics. The computation of these flags is performed with information coming from the 2nd layer from the surface of the statistical CloudSat/Calipso grid (40 levels). However, the surface elevation was not taken into account for this calculation and the precipitation flags were wrong over land surfaces higher than 480m Above Sea Level.

This bug was discovered when implementing COSPv2.1 in the LMDZ model. The sample data provided with the standalone offline COSPv2 version does not allow to detect this bug because the 153 grid points of this dataset are over ocean. Indeed, performing the Python regression test provided with the offline COSPv2 after having corrected this CloudSat bug does not result in any difference between the diagnostics fields from the reference output file and the new output file ("All fields match"). We should consider adding land points to the sample data provided with COSPv2 in order to have a wider range of cases covered with this regression test.

All the lines that have been modified or added to fix this bug have a "!PREC_BUG" comment at the end of each of these lines. As a first step to discuss this provisional bug fix, the lines that have been replaced or removed from the code are left in the code as comments.

[dustinswales]

@rodrigoguzman-lmd
I looked over your changes and see no problem merging them onto the trunk. I'll just wait for a second PMC member (@alejandrobodas @RobertPincus) to chime in, but in the meantime go ahead and remove the commented out lines and all instances of "PREC_BUG".

[RobertPincus]

I will take a look today. Meanwhile can @rodrigoguzman-lmd<https://github.com/rodrigoguzman-lmd> turn his initial suggestion about adding land points to the regression tests into an issue on the Github site?

[RobertPincus]

Hi -

I'd like to see two changes, please. In two places some magic numbers have been introduced:

Line 397 of src/simulator/quickbeam/quickbeam.F90
cloudsat_preclvl_index(:) = 39 - floor( surfelev(:)/480. )

Less important is line 869 of driver/src/cosp2_test.f90
cloudsat_preclvl_index(:) = 39 - floor( cospstateIN%surfelev(:)/480. )

I realize that these are simply replacing a previous magic number in cosp_config.F90:
cloudsat_preclvl = 39

But it would be much better if the 39 and 480. were replaced with parameters or variables from the configuration type.

[alejandrobodas]

Hi @rodrigoguzman-lmd . Thanks for this, and apologies for being slow in looking in to this. Please see my review. I've used this capability for the first time, so I hope it's clear what I've done. I've made a few in-line modifications to the lines commented out, but only to learn how to use them, please delete all the comments as requested by @dustinswales . My main comment is regarding the magic numbers, which are not constants. I think we need to resolve this in a more general way.

[rodrigoguzman-lmd]

@dustinswales , @RobertPincus and @alejandrobodas
I have made some changes in the code and taken most of your comments into account. However, I can't replace the parameter "cloudsat_preclvl = 39" by a more general expression because changing this is beyond the scope of this bug fix.

[rodrigoguzman-lmd]

Some comments on the changes appearing in the output sample dataset from this bug fix and why having landpoints in the sample data is important.
Figs_pull_request_19.pdf
Python_scripts_pull_request_19.zip

[rodrigoguzman-lmd]

After having corrected the bug fix and when performing the regression test with the 153 points output data file, I get:

[rguzman@loholt2 driver]$ python test_cosp2Imp.py data/outputs/UKMO/cosp2_output_um.ref153.nc data/outputs/UKMO/cosp2_output_um.153.nc
############################################################################################
Treating relative differences less than 0.0010000000% as insignificant
All fields match

[rodrigoguzman-lmd]

When performing the regression test with the 300 points output data file (see Issue #20 ):

[rguzman@loholt2 driver]$ python test_cosp2Imp.py data/outputs/UKMO/cosp2_output_um.ref300.nc data/outputs/UKMO/cosp2_output_um.300.nc
############################################################################################
Treating relative differences less than 0.0010000000% as insignificant
ptcloudsatflag0: 12.33 % of values differ, relative range: -1.00e+00 to -5.00e-02
ptcloudsatflag3: 0.67 % of values differ, relative range: -2.00e-01 to -2.00e-01
cloudsatpia: 0.67 % of values differ, relative range: -3.35e-01 to -3.35e-01
All other fields match.

[rodrigoguzman-lmd]

So we can only see the changes of this bug fix when enough land grid points are included in the sample dataset.
In Fig.1 we can see on the top panels the surface elevation ("orography" in the input sample dataset) for the 153 points case (left) and the 300 points case (right). In the middle panels of Fig.1, the difference between the CALIPSO lidar simulator altitude of opacity (z_opaque) computed with respect to the Sea Level (ASL) and respect to the Surface Elevation (ASE) shows that, as expected, no difference appears between these 2 variables when opaque clouds have been detected in the grid point (missing values when there are no opaque clouds at all in the grid point) over ocean (left). The difference of altitude between these 2 variables only appears over land and corresponds, as expected, to the Surface Elevation altitude where opaque clouds exist (right). The bottom panels of Fig.1 show flags 0,1,2 and 3 of the CloudSat Precip diagnostics for the 153 points case (left) and the 300 points case (right).

[rodrigoguzman-lmd]

Fig.2 shows the Surface Elevation for the 2 case scenarios again on the top panels (153 points case on the left column, 300 points case on the right column). Then, middle panels show the differences between the new output (with the bug fix) and the reference output for flags 0, 1, 2 and 3 of the CloudSat Precip diagnostics (new minus reference) , and bottom panels show the corresponding difference for the CloudSat PIA diagnostic. We can see that, as notified by the regression test, no differences appear for the 153 points case and some differences appear for the 300 points case.

[rodrigoguzman-lmd]

However, what the regression test is telling us is unconsistent with what we can see on Fig.2 where differences also exist beteween the new and the reference fields for "precipflag1" and "precipflag2" (right middle panel). This is because the regression test does not correctly account for differences between the two datasets when the reference value is equal to 0. Indeed, if we look at the «  test_cosp2Imp.py » python script, the problem comes from line 99 :

    diff[numpy.where(var_ref == 0.0)] = 0.0

As no relative difference can be computed when « var_ref » equals 0 (dividing by 0), the difference is set to 0. Nevertheless, if the value of the diagnostic is 0 in the reference data but is different from 0 in the new data (as is the case for "precipflag1" and "precipflag2") which is being tested, then the script will miss notifying such difference. In order to avoid missing these differences when performing the regression test, I suggest adding the following 2 lines (1 comment line) at line 100 of the «  test_cosp2Imp.py » script :

    # Relative difference = 1 if var_ref = 0 and var != var_ref
    diff[(var_ref == 0.0) & (var != var_ref)] = 1.0

[rodrigoguzman-lmd]

After adding these 2 lines of code in « test_cosp2Imp.py », the regression test gives :

[rguzman@loholt2 driver]$ python test_cosp2Imp.py data/outputs/UKMO/cosp2_output_um.ref300.nc data/outputs/UKMO/cosp2_output_um.300.nc
############################################################################################
Treating relative differences less than 0.0010000000% as insignificant
ptcloudsatflag0: 13.00 % of values differ, relative range: -1.00e+00 to 1.00e+00
ptcloudsatflag1: 5.00 % of values differ, relative range: 1.00e+00 to 1.00e+00
ptcloudsatflag2: 3.33 % of values differ, relative range: 1.00e+00 to 1.00e+00
ptcloudsatflag3: 7.67 % of values differ, relative range: -2.00e-01 to 1.00e+00
cloudsatpia: 16.00 % of values differ, relative range: -3.35e-01 to 1.00e+00
All other fields match.

Which is now consistent with what Fig.2 shows.

[rodrigoguzman-lmd]

Finally, I think a more thorough evaluation of this bug fix in COSPv2 CloudSat Precip diagnostics should be done with at least one GCM control run. That way we will be able to see how realistic the precipitation patterns and frequencies are at global scale in order to be sure that there are no bugs left in the code and that this bug fix is satisfying. We have performed such run with the IPSL coupled model (Atmosphere = LMDZ model) and we would be happy to compare our results to any other climate model willing to test this new CloudSat Precip diagnostics. @dustinswales , @RobertPincus and @alejandrobodas , any clues on this query? Thank you for your reviews of the code I've submitted and I hope we will be able to merge this bug fix to the master branch soon.

[RobertPincus]

@dustinswales @alejandrobodas, @rodrigoguzman-lmd notes that replacing the magic "39" is difficult. Do either of you know what this number refers to, and/or how it might be parameterized or determined from cosp_init or elsewhere?

[alejandrobodas]

@RobertPincus , I am not familiar with this number, it seems to be used in the calculation of the cloudsat precip diagnostics in quickbeam. If it is just the number of levels in the output grid minus 1 (40-1), then it should be easy to replace. However, if it's linked to a specific height, then the calculation will depend on the definition of the output grid. Unless @dustinswales can propose a simple solution, I would agree with @rodrigoguzman-lmd and treat it as a separate issue.

[dustinswales]

@rodrigoguzman-lmd
I will work on trying this in CESM2, I promise...
I see you got rid of the 39 and 480 magic numbers and replaced them with parameters from cosp_config.F90. This is perfect, the only other thing I think we should do is define cloudsat_preclvl as a parameter in cosp_config.F90. Then set it during cosp_init(), cloudsat_preclvl=Nlvgrid-1.

Nlvgrid is provided by the host-model and is set during initialization, so any dependencies on Nlvgrid should be set alongside during initialization.
If at some point someone decides to changes this statistical grid, I highly doubt this(?), but they will have to address how to compute the index for the 480-960m level, but that's not on us here.

[roj222]

Hi all, Would someone remind me who wrote this code and (from my perspective) is responsible for it ? CloudSat can't actually see anything 1-bin above the surface (due to clutter) and the operational algorithms (of which I am aware) make a precipitation determination using bins further from the surface.   I am surprised that the simulator would not take this into account. More generally, land may have a fixed height in any GCM grid cell but this is not true for the observations.   Meaning that in a given area the size of a GCM-grid cell, CloudSat will (at times) be looking for precipitation at different altitudes.  I would also certainly expect to find precipitation differences between high-terrain and neighboring low-lands.  So I am not sure how one should account for this in observational-model comparisons. Given the code situation, it seems likely to me that this has not been thought through and it may not be a good idea to be using this simulator over land.   At least, I think there needs to be a discussion here, and the limitations made clear to perspective users. Cheers, Roj.

…

On 3/20/19 9:52 AM, dustinswales wrote: @rodrigoguzman-lmd <https://github.com/rodrigoguzman-lmd> I will work on trying this in CESM2, I promise... I see you got rid of the 39 and 480 magic numbers and replaced them with parameters from cosp_config.F90. This is perfect, the only other thing I think we should do is define cloudsat_preclvl as a parameter in cosp_config.F90. Then set it during cosp_init(), cloudsat_preclvl=Nlvgrid-1. Nlvgrid is provided by the host-model and is set during initialization, so any dependencies on Nlvgrid should be set alongside during initialization. If at some point someone decides to changes this statistical grid, I highly doubt this(?), but they will have to address how to compute the index for the 480-960m level, but that's not on us here. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#19 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALoqHZb6tDvWvjTFY7LxKlDd09ikpJ1Nks5vYmc-gaJpZM4aldl7>.

[alejandrobodas]

Hi @roj222 , I believe this code was introduced in pull request #12 , and the diagnostics seem to be documented in this paper:
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017JD028213

[alejandrobodas]

Hi all, this request seems to be stalled. My personal view is that, if these diagnostics are not reliable over land, then we should mask them. I think we should still make the changes proposed in this pull request, just to avoid the re-introduction of the bug if at some point the diagnostics are calculated over land.
I would like to hear your opinions on this.

Cheers,
Alejandro

[rodrigoguzman-lmd]

Ok.

[dustinswales]

@rodrigoguzman-lmd
We (@jenkayco ) are in the process of testing these diagnostics in CESM2 and encountering some issues over land points. COSP2 is accessed in CESM2(CAM) as an external, so it's simplest to just point to your developmental branch to test these changes. I'm actively debugging this, so rest assured I will cleanup when it's working. The changes I'm testing are from within Jen's COSP1.4 code, which is working fine. Just trying to see where the differences are originating in the COSP2 implementation.
A question for you. Are you running a version of your code w/ COSP2 using these Cloudsat diagnostics? If so, do the diagnostics look sensible over land?

[rodrigoguzman-lmd]

@dustinswales @jenkayco , We did runs with the IPSL (LMDZ) coupled model with COSPv2 inline having my temporary debugged version of the CloudSat simulator precip diagnostics. Attached to this message you will find some figures with comments regarding this land point issue.
IPSL_COSPv2_one_year_figures.pdf

[dustinswales]

@rodrigoguzman-lmd
Thanks for these plots.
These changes are working in CESM2, sorry for hijacking you pull request. There is a small change that we added regarding the level used over the land points, which is now one layer higher over land.
I will merge these changes onto the trunk w/ tag v2.1.4.