country mask file not working for masking emissions in gchp

[alix-lamotte]

Your name

Alix Lamotte

Your affiliation

University of York

What happened? What did you expect to happen?

Hello,

We [@amy1916] are using GEOSChem v14.2.3 (with the addition of the fix for KPP integration failures geoschem/geos-chem#2006) for a C48S3 simulation.

we are trying to set up a run with masked anthropogenic emissions over France using the country mask (from the countrymask_0.1x0.1.nc file). We modified the HEMCO_config.rc and ExtData.rc files as suggested here ( https://hemco.readthedocs.io/en/latest/hco-ref-guide/more-examples.html#scale-or-zero-emissions-with-a-rectangular-mask ). We followed the steps as best as we could, but we're not certain it has worked correctly.

However, we are now running into issues using the country mask file [as detailed below].

Thank you for your help in advance,
Best wishes
Alix

What are the steps to reproduce the bug?

We modified the HEMCO_config.rc and ExtData.rc files as suggested here ( https://hemco.readthedocs.io/en/latest/hco-ref-guide/more-examples.html#scale-or-zero-emissions-with-a-rectangular-mask ).
HEMCO_Config.rc
HEMCO_Config.txt

ExtData.rc
ExtData.txt

scalefactor.txt
scalefactor.txt

After many attempts, we still weren't able to get a correct run.

The latests error messages were on the countrymask_0.1x0.1.nc : unsupported subclass (not real32) for units of attribute named _FillValue in CountryID in /users/wxw515/scratch/masked_emission/countrymask_0.1x0.1.nc
The file is downloaded from the Hemco read the docs, so we are unsure why this is causing problems. In an attempt to fix this, we modified the countrymask file using python and xarray to transform 'CountryID' from float64 to int32.
But the run still didn't work (but we haven't got a clear error message).

The attached files are :
the first log file
Real32_error.log

screenshot of : ncdump -h countrymask_0.1x0.1.nc (the modified version)

the second log file
Real32_error_attempted_fix.log

Please attach any relevant configuration and log files.

No response

What GCHP version were you using?

14.2.3

What environment were you running GCHP on?

Local cluster

What compiler and version were you using?

gcc 11.2.0

What MPI library and version were you using?

OpenMPI 4.1.1

Will you be addressing this bug yourself?

Yes, but I will need some help

Additional information

No response

[yantosca]

Thanks for writing @alix-lamotte and @amy1916. I think what is happening is that GCHP is choking on the unsigned integer fields of the country file. The MAPL I/O library (which GCHP uses) assumes either float or double inputs.

What you can do is to edit the file with the netCDF operators to change the uint to a float, such as:

cp countrymask_0.1x0.1.nc bad.countrymask.0.1x0.1.nc
ncap2 -s 'CountryID=CountryID.convert(NC_FLOAT);cell_area=cell_area.convert(NC_FLOAT)' countrymask_0.1x0.1.nc out.nc
mv out.nc countrymask_0.1x0.1.nc

and then see if that works for GCHP. If it does we can make a similar fix on our end as well.

[yantosca]

Also tagging @lizziel and @yidant for reference

[lizziel]

I wonder if you need to add a _FillValue attribute to CountryID. The error countrymask_0.1x0.1.nc : unsupported subclass (not real32) for units of attribute named _FillValue in CountryID is curious because it looks like there is no _FillValue in the file for that variable. It might be that MAPL requires it but has misleading error messaging when it is not found.

[yantosca]

Thanks @lizziel, that may also be true. I thought MAPL only needed long_name and units but who knows.

[lizziel]

I just noticed your second try has errors in ExtDataGridCompMod.F90. Errors in ExtData can be further diagnosed by turning on ExtData debug prints and checking log allPEs.log. See the GCHP debugging page in ReadTheDocs for information on how to do that. https://gchp.readthedocs.io/en/stable/user-guide/debugging.html

[lizziel]

Where did you get the country mask from? I am looking at the file we have at Harvard and is different from yours:

~/gcgrid/gcdata/ExtData/HEMCO/MASKS/v2014-07 $ ncdump countrymask_0.1x0.1.nc  -h
netcdf countrymask_0.1x0.1 {
dimensions:
        lon = 3600 ;
        lat = 1800 ;
        time = UNLIMITED ; // (1 currently)
variables:
        double lon(lon) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "lon" ;
                lon:units = "degrees_east" ;
                lon:axis = "X" ;
        double lat(lat) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "lat" ;
                lat:units = "degrees_north" ;
                lat:axis = "Y" ;
        double time(time) ;
                time:standard_name = "time" ;
                time:long_name = "time" ;
                time:units = "days since 2000-01-01 00:00:00" ;
                time:calendar = "standard" ;
        int CountryID(time, lat, lon) ;
                CountryID:long_name = "country mask" ;
                CountryID:units = "count" ;
                CountryID:_FillValue = -999 ;
                CountryID:missing_value = -999 ;
        double cell_area(lat, lon) ;
                cell_area:standard_name = "area" ;
                cell_area:long_name = "area of grid cell" ;
                cell_area:units = "m2" ;

// global attributes:
                :CDI = "Climate Data Interface version 1.6.4 (http://code.zmaw.de/projects/cdi)" ;
                :Conventions = "CF-1.4" ;
                :history = "Mon Apr 27 16:43:15 2015: cdo merge countrymask_0.1x0.1.nc area.nc tmp.nc\n",
                        "Mon Apr 27 16:43:04 2015: cdo gridarea countrymask_0.1x0.1.nc area.nc" ;
                :CDO = "Climate Data Operators version 1.6.4 (http://code.zmaw.de/projects/cdo)" ;

[yantosca]

Thanks for your patience @alix-lamotte and @amy1916. @lizziel and I were just looking into this. It turns out that there is some special logic in HEMCO (when used in a non-ESMF environment like GC-Classic) to read the country mask data from disk. This logic isn't present in the MAPL library, which handles all I/O for GCHP.

But a way forward would be to create individual mask files for the countries that you are interested in, and then specify them in ExtData.rc and HEMCO_Config.rc like you would any other mask file. In GCPy 1.5.0 we have an example file that does this specifically. See: https://github.com/geoschem/gcpy/blob/main/gcpy/examples/hemco/make_mask_file.py.

[yantosca]

I added an issue to remind the GCST to update the HEMCO documentation accordingly: geoschem/HEMCO#283

[alix-lamotte]

Hello, thank you for your answers !

@lizziel the countrymask file that we got is the one from this page (https://hemco.readthedocs.io/en/stable/hco-ref-guide/more-examples.html#mask-file-examples)

@yantosca we were able to get a successful run with your solution of creating an individual mask file for each country.

Thank you very much for your answers.

[yantosca]

Thanks @alix-lamotte and @amy1916. Maybe that mask file is out of date. I will make sure that we replace it wtih the good one. I think what may have happened is that we re-generated the mask file to be in the proper format for GCHP but maybe we didn't update the example file.

I'll go ahead and close out this issue. Please let us know if you run into any other problems.

[amy1916]

Hi @yantosca. Thanks for all your help on this issue.

Also, I wonder if it would be possible to update the HEMCO readthedocs scale factor / mask example with details on what to add into the ExtData.rc file when using GCHP, alongside HEMCO_Config.rc?

I think I got tripped up on not adding the scale factor / mask lines to ExtData.rc whilst doing a similar masked emission simulation previously in GCHP but figured it out by looking at previous GitHub issues etc (Apologies if this is already somewhere and I've missed it!).

Thanks again, Amy.