[BUG/ISSUE] Error in GCHP with GFAS

[thng0526]

Hi developers,

I am a postgraduate student from The Chinese University of Hong Kong. I am currently trying to run GCHP version 13.2.1 with GFAS emission inventories. However, errors are encountered. The error message guide me to HEMCO.log but no error message was found in the log file. The related messages are listed below. You may also find the attached log file at the end of this report.

From 347168_print_out.log,

>===============================================================================
> GEOS-Chem ERROR: Error encountered in "HCO_Run"!
> -> at HCOI_GC_Run (in module GeosCore/hco_interface_gc_mod.F90)
>
>THIS ERROR ORIGINATED IN HEMCO!  Please check the HEMCO log file for
additional error messages!
>===============================================================================
>
>===============================================================================
>GEOS-Chem ERROR: Error encountered in "HCOI_GC_Run"!
> -> at Emissions_Run (in module GeosCore/emissions_mod.F90)
>===============================================================================
>
>===============================================================================
>GEOS-Chem ERROR: Error encountered in "HCO_Run"!
> -> at HCOI_GC_Run (in module GeosCore/hco_interface_gc_mod.F90)
>
>THIS ERROR ORIGINATED IN HEMCO!  Please check the HEMCO log file for
additional error messages!
>===============================================================================
>
>===============================================================================
>GEOS-Chem ERROR: Error encountered in "HCOI_GC_Run"!
> -> at Emissions_Run (in module GeosCore/emissions_mod.F90)
>===============================================================================

From HEMCO.log, the tail of the file:

>Use inorganic iodine emissions (extension module)
> HOI: HOI          81
> I2: I2          89
>TIMEZONES (i.e. OFFSETS FROM UTC) WERE READ FROM A FILE

The log and config files: 347168_error.log 347168_print_out.log ExtData.rc.txt HEMCO_Config.rc.txt HEMCO_Diagn.rc.txt HEMCO.log

I checked previous bug/issue and found similar problem with HEMCO standalone in 2019, which is closed already. https://github.com/geoschem/HEMCO/issues/25 Could you give me some advice on the issue.

Cheers, Leo NG

[LiamBindle]

The ExtData.rc and HEMCO_Config.rc updates look correct to me. @lizziel or @yantosca do you have any ideas?

For reference, the top of the stack trace is gchp_chunk_mod.F90:1191 (calling EMISSIONS_RUN) and it crashed in the first time step.

[yantosca]

@thng0526, thanks for writing. It seems like it got past convection OK. Maybe it is choking on a missing file for emissions?

You can also turn on VERBOSE: 3 in HEMCO_Config.rc and that should give you more debug output to the HEMCO.log. That might help you figure out if it is a file that is missing.

[lizziel]

This issue of not printing the error to the log file may now be fixed in 13.3 with this update: https://github.com/geoschem/HEMCO/pull/111. To get around that issue try changing this line at the top of HEMCO_Config.rc:

Logfile: HEMCO.log

to this:

Logfile: *

[thng0526]

Hi @yantosca and @lizziel. Thanks for your suggestions. The error message is now displayed and shown below.

Evaluate field GFAS_ACET

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! HEMCO ERROR: Error getting vertical index at location 1 48 : GFAS_ACET ERROR LOCATION: GET_CURRENT_EMISSIONS ERROR LOCATION: HCO_CalcEmis (HCO_CALC_MOD.F90) ERROR LOCATION: HCO_RUN (hco_driver_mod.F90)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

The full log file is also attached below. 349723_print_out.log

Could you give me some advice on where I went wrong? Thanks.

[lizziel]

Hi @thng0526, I'm glad you were able to get the error message now. I'll see if I can reproduce your issue and will investigate.

[lizziel]

I was able to reproduce this. The issue appears to be in retrieving the upper limit vertical level from the scale factors. HEMCO_Config.rc uses xyL=1:scal300 for GFAS, but upper limit scal300 is returning a very small number (1e-31) in my test run. I am not sure if anyone has tried to use this feature of HEMCO in GCHP before since GFAS is disabled by default. I'll see if I can trace where in HEMCO the problem is getting the scale value.

[thng0526]

It is great to heard that you can reproduce the issue. Thanks for letting me know your progress @lizziel. Please let me know if anything that I can help.

[lizziel]

I checked the file and it looks like the issue is because the scale factor variable mami is all missing values.

Scale factor applied in HEMCO_Config.rc:

#==============================================================================
# --- GFAS scale factors ---
#==============================================================================
(((GFAS
300 GFAS_EMITL $ROOT/GFAS/v2018-09/$YYYY/GFAS_$YYYY$MM.nc mami 2003-2021/1-12/1-31/0 C xy\
 m 1
)))GFAS

Scale factor file source defined in your ExtData.rc:

GFAS_EMITL m       N Y F%y4-%m2-01T00:00:00 none none mami          ./HcoDir/GFAS/v2018-0\
9/%y4/GFAS_%y4%m2.nc

Scale factor attributes and values in a GFAS file:

        float mami(time, lat, lon) ;
                mami:_FillValue = -1.e-31f ;
                mami:units = "m" ;
                mami:long_name = "Mean altitude of maximum injection" ;
                mami:missing_value = -1.e-31f ;
 mami =
  _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, 
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, 
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, etc

I am looking into why there are missing values in the file. Stay tuned.

[yantosca]

Tagging @kilicomu

[kilicomu]

@lizziel I've just taken a look at what seems to be the relevant GFAS file for this breaking simulation - GFAS_201001.nc. I counted missing mami values and there are definitely sensible-looking values there.

These values correspond to emission locations as you can see from the following plots, one of mami and the other of cofire:

(cofire is scaled by 10^-1 for illustrative purposes. Also, if anybody can tell me what the hell is going on with Panoply such that the legend tick values aren't coming through, that would be great...)

The GFAS dataset is relatively sparse, so I wouldn't expect much in the mami arrays. Also, from the notes on how the dataset is prepared:

Where there is no CO emission, i.e. no fire, mean altitude of maximum injection value is set to the output missing value

So I am specifically setting no fire to missing values! Is the simulation error a consequence of how a GCHP component is handling missing values, rather than a problem with the processed inventory? I'd be amazed if nobody had come across this as a problem using GC-Classic for the years it's been available.

[lizziel]

@kilicomu, thanks for looking into this. I should have taken a look at Panoply as well. I wonder if the issue is it is negative. I'll look more into it. It is indeed strange that no one would have reported an issue before this.

Re: Panoply, you need to manually adjust the tick formatting in Scale->TickFormat. Precision is one decimal place and not scientific notation by default. Try setting it to 1E instead.

[lizziel]

It looks like the issue is indeed a negative value. My run is failing here (I added the print to diagnose the point of failure):

    ! Upper level                                                                         
    IF ( isLevDct2 ) THEN
       EmisL = GetEmisL( HcoState, LevDct2, I, J )
       IF ( EmisL < 0.0_hp ) THEN
          print *, "EmisL < 0 for upper level in GetVertIndx: ",EmisL
          RC = HCO_FAIL
          RETURN
       ENDIF
       EmisLUnit = LevDct2_Unit
    ELSE

Log error message:

HEMCO: Leaving HCO_RUN (hco_driver_mod.F90) ( 1)
  --- Do convection now
 EmisL < 0 for upper level in GetVertIndx:   -1.0000013902544108E-031
 EmisL < 0 for upper level in GetVertIndx:   -1.0000002147600600E-031
 EmisL < 0 for upper level in GetVertIndx:   -1.0000012727049757E-031
 EmisL < 0 for upper level in GetVertIndx:   -1.0000026832981967E-031
 EmisL < 0 for upper level in GetVertIndx:   -9.9999892171627409E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999786377135835E-032
 EmisL < 0 for upper level in GetVertIndx:   -1.0000023306498915E-031
 EmisL < 0 for upper level in GetVertIndx:   -1.0000006849578003E-031
 EmisL < 0 for upper level in GetVertIndx:   -9.9999645317813737E-032
 EmisL < 0 for upper level in GetVertIndx:   -1.0000013902544108E-031
 EmisL < 0 for upper level in GetVertIndx:   -9.9999939191401442E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999903926570918E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999880416683901E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999833396909868E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999927436457934E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999633562870229E-032
 EmisL < 0 for upper level in GetVertIndx:   -1.0000005674083652E-031
 EmisL < 0 for upper level in GetVertIndx:   -1.0000026832981967E-031
 EmisL < 0 for upper level in GetVertIndx:   -9.9999821641966360E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999892171627409E-032
 EmisL < 0 for upper level in GetVertIndx:   -1.0000016253532810E-031
 EmisL < 0 for upper level in GetVertIndx:   -1.0000008025072354E-031
 EmisL < 0 for upper level in GetVertIndx:   -9.9999798132079344E-032
 EmisL < 0 for upper level in GetVertIndx:   -9.9999856906796885E-032
===============================================================================
GEOS-Chem ERROR: Error encountered in "HCO_Run"!
 -> at HCOI_GC_Run (in module GeosCore/hco_interface_gc_mod.F90)

THIS ERROR ORIGINATED IN HEMCO!  Please check the HEMCO log file for
additional error messages!
===============================================================================

[kilicomu]

Yeah, I am explicitly setting the missing value to -1E-31. Can't for the life of me remember why - for some reason I have it in my head that it at least used to be a defined HEMCO missing value?

I've had a look through some of the other HEMCO datasets and there doesn't seem to be a consistently used missing value. Do you have a preference for what I use?

[kilicomu]

I might be mistaking -1E-31 for 1E-31, just from glancing at the HEMCO source code. In any case, let me know what it ought to be and I can see about reworking the data.

[lizziel]

-1E-31 is the correct HCO_MISSVAL value so the files are correct. Looking more at the code it appears the problem is that GCHP regridding includes the missing values while HEMCO regridding (a2a) does not. The regridded values in GCHP can thus be slightly different from -1E-31 which makes the comparison with HCO_MISSVAL return false. This prevents assigning an injection height value 0.0 and instead allows the negative value to propagate through, ultimately causing the crash. Since HEMCO regridding ignores missing values this is an issue specific to GCHP and explains why we haven't had a report about the issue before.

Let's keep the values in the file the same. I'm going to bring up whether we can ignore missing values in MAPL regridding at my next meeting with GMAO.

@thng0526, I suggest that you not use the injection heights for GFAS until we sort this out. Even if you add code to bypass this issue, such as adding a tolerance in the comparison to HCO_MISSVAL, I suspect that incorporation of missing values in GCHP regridding will make the injection heights too low. You could try using a static injection height like for GFED, although I would not be able to advise on the accuracy of doing so.

[kilicomu]

Okay @lizziel, I'll keep things as they are for now. Give me a shout if you want me to look at anything else with this.

[thng0526]

Okay @lizziel, I will hold the GFAS simulations until the problem solves. Thanks for your help.

[lizziel]

After talking with GMAO this issue was identified as a bug in MAPL. I created a bug report here. We expect a fix in the version of MAPL included in GEOS-Chem 14.0.

@thng0526, you could try emitting at the surface only rather than use injection heights. @theloniuspunk may be able to comment on whether that is a valid approach for GFAS.

[theloniuspunk]

Yes, unless you are looking specifically at BB smoke and/or free tropopheric aerosol , emitting GFAS at the surface is a fine option. It doesn't change the results that much, and it may take care of your issues. The emitting between the surface and mami values as a default option doesn't make that much sense to me anyway (mami is the altitude of the maximum amount of injection, not the maximum altitude of injection, so the default GEOS-Chem emissions setting is inconsistent with the GFAS product anyways).

[sdeastham]

While reading this conversation, it occurred to me that our default implementation probably doesn't make sense (whether injecting everything at mami, spreading the injection evenly from the surface to apt, or something in between). Consider three cells which will be "averaged" during regridding:

Cell A: 100 kg emitted at an altitude of 20 m Cell B: 1 kg emitted at an altitude of 500 m Cell C: 0 kg emitted, no altitude specified

Assuming equal area for the three sub-cells, our current regridding would emit 101 kg of material at an altitude of 260 m, which is wrong. To get the correct answer we'd need to premultiply the emissions and altitude fields I think?

If we wanted to pull that off, it could be done by reading in a derived field of emissions * injection altitude, regridding that, then dividing by regridded emissions (probably in a HEMCO extension). In GCHP that might be possible with ExtData derived imports (if not, a lightweight interfacing GC would do it, which we have been talking about doing for a while). In GC-Classic I think an edit to the HEMCO interface code would be needed. Alternatively, we could pre-process the GFAS data by premultiplying the fields, but I instinctively despise preprocessors...

[theloniuspunk]

That's right, Seb.

When I've experimented with different GFAS implementations, I've had to regrid to my horizontal grid offline pre-HEMCO to do the regridding properly, for the reasons you stated here. I couldn't think of any way to do it correctly within HEMCO.

Related, since 2018 or 2019, the native GFAS includes an "APB" (altitude of plume bottom) along with MAMI and APT. What I ended up liking in my experiments was emitting evenly between APB and APT, though I thought of going 1 step further and weighting emissions towards MAMI (which always was in between APB and and APT in altitude, as one would expect), but I didn't go this far.

Note, none of the things I tried (surface only, surface-to-MAMI, MAMI only, and APB-to-APT) compared very well to the WE-CAN field campaign in terms of putting plumes at the right place, but possible that the NA region grid was still too coarse... or maybe GFAS just had the plume heights wrong. At this point, my sabbatical ended and I dropped things... unpublished. Ugh :/

[kilicomu]

@sdeastham Before going down that route, it's probably worth us exploring the additional injection height information that was added to the GFAS inventory in July 2018 - specifically the altitude of the plume bottom and explicit injection height. @theloniuspunk and I started looking into at some point but didn't get to the end of the investigation.

[kilicomu]

There are some GFAS files that I prepped for that available here if anybody finds themselves flush with time!

[theloniuspunk]

Yes, @kilicomu, This is what I was referring to in my post (this is Jeff). APB-APT seems like the most sensible option for post-summer-2018 sims. However, the issue of horizontal averaging that Seb brings up still applies. There's no way to do this averaging correctly in HEMCO, which is why I was doing it offline.

[kilicomu]

Yep, I'm with you. Just wanted to make sure that the extra fields in the GFAS data post July 2018 were known before going down the rabbit hole!

[lizziel]

This issue should be fixed following the upgrade to MAPL in GCHP 14.0. The MAPL update that fixes it is https://github.com/GEOS-ESM/MAPL/pull/1305. The new behavior is to read _FillValue per file per netCDF file and to use it to determine what values are missing. Those values are then reset to MAPL_UNDEF and subsequently properly handled during regridding.

I will keep this issue open pending testing with 14.0.

[kilicomu]

That's great, thanks for sorting that out.

[lizziel]

This issue will be fixed in 14.1 or sooner (possibly a 14.0.z version). See https://github.com/geoschem/geos-chem/pull/1255 and https://github.com/geoschem/MAPL/pull/21 if you would like to get a fix in sooner. The delay is because we are trying to release 14.0.0 prior to IGC which is in early June, and the fix has a component for GC-Classic that needs to be looked into more.

[lizziel]

This issue is fixed with https://github.com/geoschem/GCHP/issues/258 which will be included in 14.1.