[BUG/ISSUE] xyL inconsistent diagnostic

[barronh]

Describe the bug

There is an inconsistency between the diagnostic interpretation of xyL mass that changes by time step. The diagnostic is initialized with the xy mass split evenly over assigned layers, but subsequent time steps have the mass duplicated. I'm somewhat confused about which is right.

Background

I'm trying to use the xyL=1:N option described in the users guide.

HEMCO v2 ... for 2D fields it is legal to define a range of levels, in which case the emissions are uniformly distributed across these levels (maintaining the original total emissions).

I was expecting the sum of xy emissions to be conserved, but what I get is duplication of the 2D mass in each layer for all timesteps except for the first. If duplication is intended, perhaps, the wiki could be updated to say "(duplicating the original total emissions in each layer)." Which ever way is intended, the diagnostic should be consistent.

When checking my implementation (xyL=1:3), daily diagnostics showed 291.7% of input emissions. When I output hourly, the first has 100% of input emission and all subsequent hours have 300%.

To Reproduce

I ran with a real emission file and found the issue, but have boiled it down to something small and repeatable.

I'm using v2.2 that comes with GCv12.8.0, but in the offline mode. I am running with a single input file for a single variable NO2. To make sure the test case is easy, I've disabled all scaling factors and I am just using an output grid of equal resolution but larger in scope to not lose mass on the edges (0.1 degree x 0.1 degree). My config file has xy and xyL derived from a 2D NO2 variable. The xyz uses a 3D NO23D variable that duplicates mass in each layer. My diagnostics include NO2_XY_TOTAL, NO2_XYZ_TOTAL, and NO2_XYL_TOTAL.

The XY and XYZ diagnostics are consistently 100% and 100% of their input variables. On the first time step, the XYL diagnostic is 100% original mass. On all subsequent time steps, the XYL is 300% of the xy input mass.

Attached is a tar.gz with my setup and instructions (REAMDE) to reproduce.
xyLtest.tar.gz

Versions Etc.

HEMCO: v2.2.0 compiled stand alone from GEOS-Chem v12.8.0 distribution
Compiler: intel 19.0.1
Libraries: netcdf-4.7.3, hdf5-1.10.5
Cluster environment (uname -a): Linux atmos1 3.10.0-1062.12.1.el7.x86_64 #1 SMP Thu Dec 12 06:44:49 EST 2019 x86_64 x86_64 x86_64 GNU/Linux

[barronh]

Update

I added a bunch of print statements and I find that the DilFact on line 935 of hco_calc_mod.F[1] is 1/3 in the first hour, but 1 in all the subsequent hours. I tracked it down to line 2596-2598.[2] In standalone mode, the box heights are zero. There is a commented out warning for when this will produce odd results.

There is also an option "Vertical weights" on line 409 of hco_state_mod.F90, which can be set to disable "accurate weights."

Recommendations

I recommend two simple updates.

1. If I am not reading it wrong, the default dilution factor on line 2598 should be DilFact = 1.0_hp / REAL(UppLL-LowLL+1,hp).
2. I'd recommend finding a way to add a warning (not every I/J/K), but once.

Remaining Question

That leaves me with a question. Could I have added met fields to get more representative results?

[1] https://github.com/geoschem/HEMCO/blob/master/src/Core/hco_calc_mod.F90#L935
[2] https://github.com/geoschem/HEMCO/blob/master/src/Core/hco_calc_mod.F90#L2597
[3] https://github.com/geoschem/HEMCO/blob/master/src/Core/hco_state_mod.F90#L409

[lizziel]

@barronh, thanks for looking more into this. If I am understanding correctly there is a separate issue that setting DilFact to 1.0 (and issuing a warning message) should occur only if box height is available. If box height is zero in standalone mode then box height should be deemed unavailable, and the correct code in the ELSE should be executed:

! Compute 'accurate' dilution factor if boxheights are available
IF ( HcoState%Options%VertWeight .AND. &
     ASSOCIATED( HcoState%Grid%BXHEIGHT_M%Val ) ) THEN

   ! Height of grid box of interest (in m)
   dh = HcoState%Grid%BXHEIGHT_M%Val(I,J,L)

.... more code ....

   ! compute dilution factor: the new flux should emit the same mass per
   ! volume, i.e. flux_total/column_total = flux_level/column_level
   ! --> flux_level = fluxtotal * column_level / column_total.
   IF ( h2 > h1 ) THEN
      DilFact = dh / ( h2 - h1 )
   ELSE
      !write(*,*) 'Warning: GetDilFact h2 not greater than h1!!! ',LowLL,UppLL,L,EmisL1,EmisL2,EmisL1Unit,EmisL2Unit
      DilFact = 1.0_hp
   ENDIF

! Approxiate dilution factor otherwise
ELSE

   DilFact = 1.0_hp / REAL(UppLL-LowLL+1,hp)
ENDIF

Is this your understanding too?

[lizziel]

To clarify further, by "should" I mean HEMCO would need to be changed to do this. This looks like a bug.

[barronh]

I think I follow your point. Ideally, the code should have noted the unavailability and used a uniform distribution.

In addition, is there a case where h2 is not greater than h1 and the routine should return 1?

[lizziel]

is there a case where h2 is not greater than h1 and the routine should return 1?

That's a good question. I'll think more about the potential cases, although it probably safest to add that even if it should never happen.

[lizziel]

I looked more into a fix for this bug. The problem is HcoState%Grid%BXHEIGHT_M%Val is associated within subroutine HCO_CalcVertGrid even when box height is not available. This is because there is a call to HCO_ArrAssert for box height during which HcoState%Grid%BXHEIGHT_M and HcoState%Grid%BXHEIGHT_M%Val are both initialized if not already associated. The call to HCO_ArrAssert occurs before any checks for whether box height is actually available.

The dilution factor then ends up being wrong because logic in GetDilFact assumes box height is available and non-zero if HcoState%Grid%BXHEIGHT_M%Val is associated. However, the aforementioned code allows for the possibility that the value is associated but zero. This results in a zero-value delta h calculation between levels, and a resultant dilution factor of 1.0.

HCO_CalcVertGrid is called within ExtState_SetFields which itself is called after the first call HCO_Run. Hence the issue is only seen after the first timestep.

I added an update to nullify HcoState%Grid%BXHEIGHT_M%Val within HCO_CalcVertGrid if box height is unavailable. This fixes the erroneous duplication of emissions across all levels if using the xyL option with 2D emissions. The emissions are now uniformly distributed such that total emissions are preserved.

For the case where box height is available, I removed the automatic setting of dilution factor to 1.0 if h2 not greater than h1. If this occurs the model now exits with an error. I do not think there is a valid reason for h2 to be <= h1 when box height is available and criteria are met to call GetDilFact. However, if this ends up introducing a bug for a test case outside of our benchmark runs then we will fix it as we get user reports.

Thank you @barronh for reporting this and providing a test run directory and data. The provided materials were very helpful for reproducing the issue and honing in on a fix quickly. This update will go into HEMCO 3.0.0 and I will post a commit hash when it is available.

[barronh]

I am so grateful that you tracked this down and fixed it! Thank you!

[lizziel]

No problem! The fix is in commit a8f775f. You may notice the incorrect module name is used in the error handling message. That is fixed in a subsequent commit.