[BUG/ISSUE] Metrics.py showing 0.000 year lifetime of CH3CCl3 wrt to OH

[Twize]

Description of the problem

Hi, I recently ran a 2-month C48 simulation using GCHP v13.3.2 using the default settings. The run completed fine with no issues and all of the output files (including the Metrics ones) were produced. I ran metrics.py, and I get the following:

============================================================================== GEOS-Chem FULL-CHEMISTRY SIMULATION METRICS

Simulation start : 20170101 000000z Simulation end : 20170301 000000z ==============================================================================

Mass-weighted mean OH concentration = 10.67896387472 x 10^5 molec cm-3

CH3CCl3 lifetime w/r/t tropospheric OH = 0.0000 years

CH4 lifetime w/r/t tropospheric OH = 10.5778 years

The other two values seem reasonable, and I compared with a GEOS-Chem Classic v13.3.2 run over a similar time frame, and obtained close values for the mass-weighted mean OH and CH4 lifetimes. I also checked the 'LossOHbyMCFcolumnTrop' values inside the Metrics NetCDF files using ncdump, and compared them with those from GEOS-Chem Classic and get relatively similar values (same orders of magnitude). I ran a different 1-month test, and get a non-zero value for the CH3CCl3 lifetime there when I run Metrics.py, so I am not quite sure of the cause here.

Is this likely a bug/calculation error in the Metrics.py code? Or might this point to a bigger issue with my simulation?

Files

• Run log: gchp.log
• Example metrics output file: GEOSChem.Metrics.20170215_2200z.nc4.zip

Not sure what other files may be useful, let me know if there are others that I can upload.

Software versions

• Compilers (Intel or GNU, and version): gcc v9.2.0
• NetCDF version: netcdf-fortran v4.6.0
• MPI version: openMPI v4.0.3
• ESMF version: ESMF v8.3.1

[lizziel]

Hi @Twize, it is hard to say what the issue is given the information. The best way to debug this is to open the python script metrics.py and follow along the code on how the CH3CCl3 lifetime is computed from the data. You can add python prints to check values at various stages to see where a zero is introduced. This should give insight about what the problem is.

[Twize]

Hi @lizziel, thanks, I will give this a go tomorrow and see if I can figure out where things might be going wrong.

[Twize]

Hi @lizziel, I found the source of the issue. I'm not sure if it has to do with my Anaconda environment (Python 3.7.13) and/or Xarray version (0.20.2), or if there is something missing from the metrics.py code, but the values in the Metrics dataset had dtype np.float32 after being read into the code instead of the intended np.float64 (which is specified in the comments in the function lifetimes_wrt_oh()). This caused the np.nansum on line 245 to exceed the maximum precision of np.float32 which is 1E38, causing the result of the sum to be infinity.

Changing line 245 of metrics.py from: oh_loss_by_mcf = np.nansum(ds["LossOHbyMCFcolumnTrop"].values) to: oh_loss_by_mcf = np.nansum(ds["LossOHbyMCFcolumnTrop"].values.astype('float64'))

produces the expected result, and yields a reasonable value of ~6.24 years for for the lifetime of CH3CCl3 w.r.t. to OH.

[lizziel]

Hi @Twize, thanks for reporting back on what the issue was. It is strange it is not preserving the original type of the data. I tried to reproduce the issue using xarray 0.20.2 but the data type is being correctly read as float64. To be safe we can add the astype command to the standard version.

[Twize]

Yeah, hard to say what exactly in my Python environment might be causing this behaviour (possibly Numpy? I'm using v1.21.4), as I have run metrics.py many times on a different server with a different Anaconda environment and never had this issue. Might be worth explicitly specifying astype('float64') in the standard version like you mentioned just to be safe.

I should add that I am using Intel Python3 as my Anaconda base on Scinet Niagara, so it may also have something to do with that but not sure.

[lizziel]

Aha, I just realized my test used GC-Classic and not GCHP, hence my output was REAL8. GCHP diagnostics are limited by the NASA MAPL library to only output as REAL4, hence your dataset includes floats rather than doubles. I will make an issue on GEOS-Chem repo about this.

[Twize]

Glad you were able to figure out the source of the issue!