Bugfix/vertical regridding

[nicholasbalasus]

Name and Institution (Required)

Name: Nick Balasus Institution: Harvard University

Confirm you have reviewed the following documentation

• [x] Contributing guidelines

Describe the update

So far, this PR fixes the COLLAPSE routine, removes INFLATE, and removes some old GEOS-4 code.

I've tried to make it so this is how the code works:

• If your input is 72/73 (native GEOS-5) and your output is 47/48 (reduced GEOS-5), use COLLAPSE
• If your input is 102/103 (native GISS) and your output is 74/75 (reduced GISS), use COLLAPSE
• Otherwise, use MESSy

Expected changes

Previously, the COLLAPSE routine was not working as it should for model levels or edges. For example, for model levels, when it was trying to collapse two levels, it would actually just use the value of the first. When it would try to collapse four levels, it would average the bottom three. Now, it does pressure weighted averaging correctly when collapsing model levels and samples at model edges when collapsing model edges. This will impact all 47L simulations (meteorology, restart files, prescribed losses, etc.).

Related Github Issue(s)

So far, this will already solve #232. It's not there yet for https://github.com/geoschem/geos-chem/issues/1904 though which is why I mark it as a draft. Here is what still needs to be done that I need some guidance on, maybe from @jimmielin.

• Because I've removed the INFLATE behavior as Daniel suggested, I've run into a few issues, such as the default OH fields or CH4 restart file for 72 L being on 47 L. Ideally, this would be changed, but I'm guessing dropping this 47 L -> 72 L behavior would cause widespread issues. I'd like to call MESSy in these cases to do a better job of regridding than INFALTE used to do (just copying one reduced level to be two identical native levels), but the Restart file (for example) doesn't have the PS term that NC_GET_SIGMA_LEVELS needs. As is, MESSy is being called but is failing.
• As discussed, we need to be careful about how we regrid special units. I've accepted that for a lot of cases, we'll just have to print something like "WARNING: you are doing a pressure-weighted regridding which may not always be appropriate." But I actually think in most cases, we are almost OK. For restart files and for Cl, we regrid mol/mol dry using pressure coordinates. I think we can just convert these to mol/mol (if this is easy) before doing the pressure regridding. But for OH, it is in kg/m3, so we'll need to convert to kg/kg before moving on. So I am thinking something like this:

  IF (HEMCO config file units says unitless) AND (netCDF file units says mol/mol dry)
   convert to mol/mol, vertically regrid, convert back to mol/mol dry
  ELSE IF (we are in HEMCO units of kg/m3)
   convert to kg/kg, vertical regrid, convert back to kg/m3
  ELSE
   leave units as is, conduct vertical regridding

  The only limitation here is being able to do those unit conversions at non-regridded and regridded-resolution, which may be a big obstacle.

[jimmielin]

Thanks @nicholasbalasus. Here are some thoughts on some of the points:

Because I've removed the INFLATE behavior as Daniel suggested, I've run into a few issues, such as the default OH fields or CH4 restart file for 72 L being on 47 L. Ideally, this would be changed, but I'm guessing dropping this 47 L -> 72 L behavior would cause widespread issues. I'd like to call MESSy in these cases to do a better job of regridding than INFALTE used to do (just copying one reduced level to be two identical native levels), but the Restart file (for example) doesn't have the PS term that NC_GET_SIGMA_LEVELS needs. As is, MESSy is being called but is failing.

I think we could use something similar to how MODEL_WRF && MODEL_CESM handles it now, which is to hard-code the 72 level centers into hcoio_read_std_mod.F90 (maybe elsewhere would be better.) and do the regridding on that, to at least maintain backward compatibility with existing restart files. It might necessitate a major version change if we weren't, and I'd appreciate hearing from the GCST about how to proceed here.

[nicholasbalasus]

Thanks @jimmielin! I am thinking now that we shouldn't do any "magic" in HEMCO w.r.t. a different vertical regridding based on units (especially since HEMCO is "unitless" in many cases such as with restart files). I think the solution for the methane simulation is to start with OH fields that are in mol/mol instead of kg/m3 (but this is a methane simulation issue, not a HEMCO issue). Let me know if you agree.

If so, the only remaining issue would be 47 -> 72.

[toddmooring]

@nicholasbalasus @jimmielin I think there might be a bug in your proposed changes to the DO loop at line 1300 (936) of the old (new) hco_interp_mod.F90. Basically in that loop you are trying to average NLEV levels together by summing them with appropriate weights.

Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) is first written into in the line above the DO loop, weighted by the amount that level InLev1 (of the fine grid) should contribute to the total. The purpose of the DO loop is to then properly include the contributions from the other NLEV-1 levels. So it seems to me that you should continue to index over 1, NLEV-1 rather than 1, NLEV. Or you could index over 2, NLEV.

Please let me know if I am somehow misunderstanding how this code is supposed to work.

[nicholasbalasus]

Good catch, thanks @toddmooring!

[nicholasbalasus]

This PR should be ready for review. I have tested it with this bash script: run-bugfix-sh.txt. Then, I use this python script: test-bugfix-py.txt. The result is that before the bug fix, the 47 and 72 layer simulations start with different masses of CH4. After the bug fix, the 47 and 72 layer simulations start with the same masses of CH4.

~~ Results before bugfix/vertical-regridding ~~ 72 L simulation (t=0) -> 5123.24 Tg of CH4 47 L simulation (t=0) -> 5129.15 Tg of CH4

~~ Results after bugfix/vertical-regridding ~~ 72 L simulation (t=0) -> 5123.24 Tg of CH4 47 L simulation (t=0) -> 5123.24 Tg of CH4

I plot vertical profiles that are natively 72L to show that the COLLAPSE bug is fixed (going from 72 to 47 L). You can see that before, when COLLAPSING two levels, the value of the previous level was being taken.

I also plot vertical profiles that are natively 47L to show that the INFLATE is working as it did before (going from 47 to 72 L). Note the vertical offset between Before HEMCO and After HEMCO (47L simulation) that is due to the default OH fields (OH/v2014-09/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc) having a different vertical coordinate than GEOS-Chem.

[nicholasbalasus]

This will change the output of any simulation that regrids from 72 L to 47 L (that is, all 47 L simulations because of the native resolution of the met fields). 72 L simulations should in theory not be impacted, as INFLATE behavior has remained (though I am naive w.r.t. non-CH4 simulations).

[yantosca]

Thanks @jimmielin and @nicholasbalasus. I'll merge this into 14.2.1 and test.

[yantosca]

Integration tests are running.

[yantosca]

@nicholasbalasus @jimmielin: The integration test for 05x0625_NA_47L_merra2_fullchem failed with this error:

HEMCO: Opening /n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/data/ExtData/HEMCO/AEIC2019/v2022-03/2019_monmean/AEIC_monmean_201907.0.5x0.625.36L.nc
 Invalid level standard name: ^_^@^@^@^@^@^@^@<9C>y<9D>f<FC>^?^@^@<A0>y<9D>f<FC>^?^@^@px<9D>f<FC>^?^@^@
^A^@^@^@^@^@^@^@^G^@^@^@^@^@^@^@<C0>y<9D>f<FC>^?^@^@<9E><C0>+^A^@^@^@^@<E0><9C><9D>f<FC>^?^@^@4<90><9D>f<FC>^?^@^@,<90><9D>f<FC>^?^@^@X<91><9D>f<FC>^?^@^@ <86><9D>f<FC>^?^@^@`ق<C2>^@^@^@^@                ^G^@^@^@^@^@^@^@<80><F8><95><E6>^A^@^@^@kg/m2/s                        ^@^D^@^@^@^@^@^@^@^@^@^@^@^A^A^@^@^D
^@^@^@^@^@^@^@^A^@^@^@^@^@^@^@@^B^@^@i^A^@^@$^@^@^@^A^@^@^@^E^@^@^@^@^@^@^@<D0>^N^C<E7>^A^@^@^@^C^@^@^@
^@^@^@^@^@^@^@^@^@^@<F4> in /n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/data/ExtData/HEMCO/AEIC2019/v2022-03/2019_monmean/AEIC_monmean_201907.0.5x0.625.36L.nc

HEMCO ERROR: Cannot read sigma levels of /n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/data/ExtData/HEMCO/AEIC2019/v2022-03/2019_monmean/AEIC_monmean_201907.0.5x0.625.36L.nc

HEMCO ERROR: Error encountered in routine HCOIO_Read!

HEMCO ERROR: Error in HCOIO_DATAREAD called from HEMCO ReadList_Fill: AEIC19_MONMEAN_NO
 --> LOCATION: ReadList_Fill (HCO_ReadList_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!
===============================================================================

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x1533d3131b7f in ???
#1  0x1533d4335859 in sflush
        at ../.././libgfortran/io/unix.h:89
#2  0x5987cf in __hco_interface_gc_mod_MOD_hcoi_gc_run
===============================================================================
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!
===============================================================================

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x1533d3131b7f in ???
#1  0x1533d4335859 in sflush
        at ../.././libgfortran/io/unix.h:89
#2  0x5987cf in __hco_interface_gc_mod_MOD_hcoi_gc_run
#3  0x4e3fe7 in __emissions_mod_MOD_emissions_run
        at /n/home09/ryantosca/IT/14.2.r19/GCC_14.2.1_r19/CodeDir/src/GEOS-Chem/GeosCore/emissions_mod.F90:185
#4  0x407635 in geos_chem
        at /n/home09/ryantosca/IT/14.2.r19/GCC_14.2.1_r19/CodeDir/src/GEOS-Chem/Interfaces/GCClassic/main.F90:702
#5  0x40cf2e in main
        at /n/home09/ryantosca/IT/14.2.r19/GCC_14.2.1_r19/CodeDir/src/GEOS-Chem/Interfaces/GCClassic/main.F90:32
srun: error: holy2c12303: task 0: Segmentation fault (core dumped)

[nicholasbalasus]

Thanks, @yantosca. Looks like a 36L file, so it probably called MESSy and MESSy didn't work. I'll take a look at this tomorrow.

[yantosca]

@nicholasbalasus: The lev dimension of the HEMCO/AEIC2019/v2022-03/2019_monmean/AEIC_monmean_201907.0.5x0.625.36L.nc file is just integers but not eta levels.

 lev = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 
    21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36 ;

From the technical note HEMCO/AEIC2019/v2022-03/2019_monmean/AEIC_2019_technical_note.pdf:

Vertically, there are 36 layers of varying thickness. The data orientation is “positive”, such that the first layer is adjacent to the Earth’s surface and increasing number corresponds to increasing altitude. The layers are defined identically to the first 36 layers used in the NASA Global Modeling and Assimilation Office (GMAO) Modern Era Retrospective for Research and Analysis version 2 (MERRA-2). The Goddard Earth Observation System (GEOS) which produces MERRA-2 uses a hybrid vertical coordinate. The pressure at the lower edge of the grid cell at latitude index j, longitude index i, and layer k is defined as

p(i, j, k) = A(k) + B(k) × p_sfc (i, j)

where A(k) and B(k) are layer-specific constants given in Appendix A, and p_sfc is the surface pressure at the base of the column. Since there are 36 layers and the pressure must also be defined at the top of the uppermost layer, we provide the first 37 values for A and B.

So this is an edge case, where the data doesn't extend up to level 37 (~170 hPa), which makes sense, as this is probably higher than cruise altitude for typical airplanes. I think it was decided to truncate the data to save space since the horizontal resolution is 0.5 x 0.625.

[nicholasbalasus]

Added a fix to handle these files.

[yantosca]

Thanks @nicholasbalasus. Testing again w/ the fix.

[yantosca]

@nicholasbalasus @jimmielin: Some good news and bad news. The good news is that the fix for AEIC now allows the 47L fullchem integration test to pass:

gc_05x0625_NA_47L_merra2_fullchem...................Execute Simulation....PASS

but the carbon simulation integration test now fails:

gc_4x5_merra2_carbon................................Execute Simulation....FAIL

The log file says:

********************************************
* B e g i n   T i m e   S t e p p i n g !! *
********************************************

---> DATE: 2019/01/01  UTC: 00:00
 HEMCO already called for this timestep. Returning.
Carbon_Gases: Global OH is set to zero!

I think the global OH file used in the GCClassic carbon simulation is a 47-L file. Will test further and suggest a fix.

[nicholasbalasus]

Thanks @yantosca. A 47L -> 72L should be handled by this (and is in this example). But that is with the CH4 simulation...

I'm not sure if this is related, but I find for the CH4 simulation that when I save out BOH at t=0 (using the Boundary Condition collection), it is all equal to zero. Then at t=1 it looks right. The same is true for BCl. I wonder if there's a catch for this in the carbon simulation but not the methane simulation.

[nicholasbalasus]

^ But the behavior I suggested should not be relevant to this PR...

[yantosca]

Hi @nicholasbalasus. This might be an old-fashioned input error.

    IF ( Input_Opt%amIRoot ) THEN
       IF ( useGlobOH ) THEN
          IF ( useGlobOHv5 ) THEN
             WRITE( 6, 100 ) 'GLOBAL_OH_GCv5'
          ELSE
             WRITE( 6, 100 ) 'GLOBAL_OH'
          ENDIF
 100   FORMAT( 'Carbon_Gases: Using global OH oxidant field option: ', a )

       ELSE
          WRITE( 6, 110 )
 110      FORMAT( 'Carbon_Gases: Global OH is set to zero!' )
       ENDIF

    ENDIF

I bundled this PR with some updates for the carbon simulation (see https://github.com/geoschem/geos-chem/pull/1916). WIll dig further.

[yantosca]

When applying PR https://github.com/geoschem/geos-chem/pull/1923, the carbon integration test now passes, but the Hg integration tests now fail. (The TOMAS failures are known issues that will be resolved in 14.3.0).

==============================================================================
GEOS-Chem Classic: Execution Test Results

GCClassic #47fe2b7 GEOS-Chem submodule update: Merge PR #1916 (Fixes for CH4 in carbon sim)
GEOS-Chem #e5d91c03f Merge PR #1923 (Update carbon_gases_mod for consistency w/ PR #1916)
HEMCO     #daf54d0 Merge PR #235 (Bug fix for vertical regridding)

Using 24 OpenMP threads
Number of execution tests: 26

Submitted as SLURM job: 809695
==============================================================================

Execution tests:
------------------------------------------------------------------------------
gc_05x0625_NA_47L_merra2_CH4........................Execute Simulation....PASS
gc_05x0625_NA_47L_merra2_fullchem...................Execute Simulation....PASS
gc_4x5_47L_merra2_fullchem..........................Execute Simulation....PASS
gc_4x5_47L_merra2_fullchem_TOMAS15..................Execute Simulation....FAIL
gc_4x5_47L_merra2_fullchem_TOMAS40..................Execute Simulation....FAIL
gc_4x5_merra2_aerosol...............................Execute Simulation....PASS
gc_4x5_merra2_carbon................................Execute Simulation....PASS
gc_4x5_merra2_CH4...................................Execute Simulation....PASS
gc_4x5_merra2_CO2...................................Execute Simulation....PASS
gc_4x5_merra2_fullchem..............................Execute Simulation....PASS
gc_4x5_merra2_fullchem_aciduptake...................Execute Simulation....PASS
gc_4x5_merra2_fullchem_APM..........................Execute Simulation....PASS
gc_4x5_merra2_fullchem_benchmark....................Execute Simulation....PASS
gc_4x5_merra2_fullchem_complexSOA...................Execute Simulation....PASS
gc_4x5_merra2_fullchem_complexSOA_SVPOA.............Execute Simulation....PASS
gc_4x5_merra2_fullchem_LuoWd........................Execute Simulation....PASS
gc_4x5_merra2_fullchem_marinePOA....................Execute Simulation....PASS
gc_4x5_merra2_fullchem_RRTMG........................Execute Simulation....PASS
gc_4x5_merra2_Hg....................................Execute Simulation....FAIL
gc_4x5_merra2_metals................................Execute Simulation....PASS
gc_4x5_merra2_POPs_BaP..............................Execute Simulation....PASS
gc_4x5_merra2_tagCH4................................Execute Simulation....PASS
gc_4x5_merra2_tagCO.................................Execute Simulation....PASS
gc_4x5_merra2_tagO3.................................Execute Simulation....PASS
gc_4x5_merra2_TransportTracers......................Execute Simulation....PASS
gc_4x5_merra2_TransportTracers_LuoWd................Execute Simulation....PASS

Summary of test results:
------------------------------------------------------------------------------
Execution tests passed: 23
Execution tests failed: 3
Execution tests not yet completed: 0

[yantosca]

All GCHP integration tests pass.

==============================================================================
GCHP: Execution Test Results

GCClassic #a5676b8 GEOS-Chem submodule update: Merge PR #1916 (Fixes for CH4 in carbon sim)
GEOS-Chem #e5d91c03f Merge PR #1923 (Update carbon_gases_mod for consistency w/ PR #1916)
HEMCO     #daf54d0 Merge PR #235 (Bug fix for vertical regridding)

Number of execution tests: 5

Submitted as SLURM job: 809728
==============================================================================

Execution tests:
------------------------------------------------------------------------------
gchp_merra2_fullchem................................Execute Simulation....PASS
gchp_merra2_fullchem_benchmark......................Execute Simulation....PASS
gchp_merra2_fullchem_RRTMG..........................Execute Simulation....PASS
gchp_merra2_tagO3...................................Execute Simulation....PASS
gchp_merra2_TransportTracers........................Execute Simulation....PASS

Summary of test results:
------------------------------------------------------------------------------
Execution tests passed: 5
Execution tests failed: 0
Execution tests not yet completed: 0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  All execution tests passed!  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Furthermore all GCHP integration tests were identical to the previous integraton test for PR #1903 plus PR https://github.com/geoschem/hemco/pull/234:

Checking gchp_merra2_fullchem
   -> No differences in OutputDir
   -> No differences in Restarts

Checking gchp_merra2_fullchem_benchmark
   -> No differences in OutputDir
   -> No differences in Restarts

Checking gchp_merra2_fullchem_RRTMG
   -> No differences in OutputDir
   -> No differences in Restarts

Checking gchp_merra2_tagO3
   -> No differences in OutputDir
   -> No differences in Restarts

Checking gchp_merra2_TransportTracers
   -> No differences in OutputDir
   -> No differences in Restarts

[yantosca]

The error in the Hg simulation is:

HEMCO: Opening /n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/data/ExtData/HEMCO/MERCURY/v2018-04/ocean_fixed.nc

HEMCO ERROR:  ModelLev_Interpolate was called but MESSy should have been used: GLOBAL_OCEAN

HEMCO ERROR: ERROR 0
 --> LOCATION: ModelLev_Interpolate (hco_interp_mod.F90)

HEMCO ERROR: ERROR 16
 --> LOCATION: HCOIO_READ (HCOIO_READ_STD_MOD.F90)

Which is caused by the level index being = 1:

data:

 time = "2007-01-01", "2007-02-01", "2007-03-01", "2007-04-01", "2007-05-01", 
    "2007-06-01", "2007-07-01", "2007-08-01", "2007-09-01", "2007-10-01", 
    "2007-11-01", "2007-12-01" ;

 lat = -89, -86, -82, -78, -74, -70, -66, -62, -58, -54, -50, -46, -42, -38, 
    -34, -30, -26, -22, -18, -14, -10, -6, -2, 2, 6, 10, 14, 18, 22, 26, 30, 
    34, 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86, 89 ;

 lon = -180, -175, -170, -165, -160, -155, -150, -145, -140, -135, -130, 
    -125, -120, -115, -110, -105, -100, -95, -90, -85, -80, -75, -70, -65, 
    -60, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 
    20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 
    110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170, 175 ;

 lev = 1 ;
}

which should not be making it into the regridding code since there is only surface data.

[nicholasbalasus]

So is this the only xyz file in all of GCClassic with a z dimension of size 1? I'm wondering if we can adjust the HEMCO dimensions to just xy but I'm not really familiar with what that file is supposed to be.

[nicholasbalasus]

Though I'm also not sure how this got a IsModelLev value of True.

[nicholasbalasus]

Looks like the change would need to be made here...

https://github.com/geoschem/geos-chem/blob/4722f288e90291ba904222f4bbe4fc216d17c34a/GeosCore/mercury_mod.F90#L3267-L3277

If changing xyz to xy works, then we would just need to change OCEAN_CONC(:,:,1) to OCEAN_CONC(:,:)

What do you think @yantosca ?

[yantosca]

Thanks @nicholasbalasus. I can make a separate PR to fix the Hg code and HEMCO config files in the geoschem/geos-chem repo. Will let you know.

[nicholasbalasus]

Closing, as this was merged

https://github.com/geoschem/HEMCO/commit/daf54d07ed16d0948302070855468b6b27fdd9a5