Bugfix/vertical regridding

[nicholasbalasus]

Name and Institution (Required)

Name: Nick Balasus
Institution: Harvard University

Confirm you have reviewed the following documentation

• 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 geoschem/geos-chem#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]

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 geoschem/geos-chem#1916). WIll dig further.

[yantosca]

When applying PR geoschem/geos-chem#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 #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

daf54d0