[Bug] COLLAPSE subroutine is not mass conserving

[nicholasbalasus]

Name and Institution (Required)

Name: Nick Balasus
Institution: Harvard University

Confirm you have reviewed the following documentation

• Support guidelines
• User manuals
• Debugging GEOS-Chem and HEMCO errors
• Current and past Github issues

Description of your issue or question

We are currently using a 72-layer restart file as initial conditions for a 47-layer GEOS-Chem CH4 run. Using the BoundaryCondition collection, I am able to output how this restart file has been regridded with HEMCO at t = 0. This reveals that the mass has not been conserved.

I run GEOS-Chem using this bash script:

#!/bin/bash

# Environment file
source /n/home06/nbalasus/envs/gcclassic.rocky+gnu10.minimal.env

# Fucntion to run GCClassic (version = 14.2.1)
function run_gc {

    # (1) Clone, compile, and make a run directory

    dir=$1
    mkdir "${dir}"
    cd "${dir}"
    git clone https://github.com/geoschem/GCClassic.git
    cd GCClassic
    git checkout c1eb4a2 # most recent 14.2.1 commit as of 1 Aug 2023 @ 7:09 PM
    git submodule update --init --recursive
    cd run
    if [[ "${dir}" == *"72"* ]]; then
        lev="72"
        c="3\n1\n2\n1\n${dir}\ngc_2x25_merra2_${lev}L_CH4\nn\n" # CH4, MERRA2, 72 levels
    else
        lev="47"
        c="3\n1\n2\n2\n${dir}\ngc_2x25_merra2_${lev}L_CH4\nn\n" # CH4, MERRA2, 47 levels
    fi
    printf ${c} | ./createRunDir.sh
    cd "${dir}/gc_2x25_merra2_${lev}L_CH4/build"
    cmake ../CodeDir -DRUNDIR=..
    make -j
    make install
    cd "${dir}/gc_2x25_merra2_${lev}L_CH4/"

    # (2) Modify configuration files

    # Modify HISTORY.rc: 
    # - write CH4 every hour (instantaneous)
    # - write boundary conditions every 3 hours (instantaneous with CH4 + dry air)
    sed -i -e "s|'CH4',|#'CH4',|g" \
        -e "s|'Metrics',|#'Metrics',|g" \
        -e "s|#'BoundaryConditions',|'BoundaryConditions',|g" \
        -e "s|SpeciesConc.frequency:      00000100 000000|SpeciesConc.frequency:      00000000 010000|g" \
        -e "s|SpeciesConc.duration:       00000100 000000|SpeciesConc.duration:       00000000 010000|g" \
        -e "s|SpeciesConc.mode:           'time-averaged'|SpeciesConc.mode:           'instantaneous'|g" \
        -e "s|'SpeciesConcMND_?ALL?          ',|#'SpeciesConcMND_?ALL?          ',|g"  HISTORY.rc
    sed -i "/BoundaryConditions.fields:     'SpeciesBC_?ADV?             ',/a\ 'Met_AD',\n 'Met_DELP',\n 'Met_DELPDRY'," HISTORY.rc

    # Modify geoschem_config.yml
    # - modify start and end date
    sed -i -e "s|20190101|20170101|g" \
        -e "s|20190201|20170102|g" geoschem_config.yml

    # (3) Get the restart file
    cp /n/jacob_lab/Users/tmooring/GCClassic-correcttropopause_archdirs/gc_2x25_merra2_CH4_ct_CLSBCr_prod1/Restarts/GEOSChem.Restart.20170101_0000z.nc4 ./Restarts/

    # (4) Copy the run script and submit it
    cp runScriptSamples/operational_examples/harvard_cannon/geoschem.run .
    sed -i -e "s|--mem=15000|--mem=64000|g" geoschem.run
    sbatch geoschem.run
}

for dir in "/n/holyscratch01/jacob_lab/${USER}/conservation_test_72_levels" "/n/holyscratch01/jacob_lab/${USER}/conservation_test_47_levels"; do
    run_gc "$dir"
done

Then, using Python, we can see that the mass is not conserved (6.25 Tg have been created):

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np

mw_ch4 = 16.04
mw_air = 28.96
kg_to_g = 1e3
g_to_tg = 1e-12

with xr.open_dataset("/n/home06/nbalasus/holyscratch01/conservation_test_47_levels/gc_2x25_merra2_47L_CH4/OutputDir/GEOSChem.BoundaryConditions.20170101_0000z.nc4") as processed_47L:
    kg_dry_air = processed_47L["Met_AD"].isel(time=0)
    mol_dry_air = kg_dry_air * kg_to_g / mw_air
    mol_CH4_per_mol_dry_air = processed_47L["SpeciesBC_CH4"].isel(time=0)
    mol_CH4 = mol_CH4_per_mol_dry_air * mol_dry_air
    tg_CH4 = mol_CH4 * mw_ch4 * g_to_tg
    print(f"{tg_CH4.sum().values:.2f} Tg of CH4")
    # 5082.72 Tg of CH4

with xr.open_dataset("/n/home06/nbalasus/holyscratch01/conservation_test_72_levels/gc_2x25_merra2_72L_CH4/OutputDir/GEOSChem.BoundaryConditions.20170101_0000z.nc4") as processed_72L:
    kg_dry_air = processed_72L["Met_AD"].isel(time=0)
    mol_dry_air = kg_dry_air * kg_to_g / mw_air
    mol_CH4_per_mol_dry_air = processed_72L["SpeciesBC_CH4"].isel(time=0)
    mol_CH4 = mol_CH4_per_mol_dry_air * mol_dry_air
    tg_CH4 = mol_CH4 * mw_ch4 * g_to_tg
    print(f"{tg_CH4.sum().values:.2f} Tg of CH4")
    # 5076.47 Tg of CH4

I find this to be related to the COLLAPSE subroutine in HEMCO.

HEMCO/src/Core/hco_interp_mod.F90

Lines 1220 to 1311 in 477c7e8

	!------------------------------------------------------------------------------
	! Harmonized Emissions Component (HEMCO) !
	!------------------------------------------------------------------------------
	!BOP
	!
	! !IROUTINE: COLLAPSE
	!
	! !DESCRIPTION: Helper routine to collapse input levels onto the output grid.
	! The input data is weighted by the grid box thicknesses defined on top of
	! this module. The input parameter T determines the time slice to be considered,
	! and MET denotes the met field type of the input data (4 = GEOS-4 levels, GEOS-5
	! otherwise).
	!\\
	!\\
	! !INTERFACE:
	!
	SUBROUTINE COLLAPSE ( Lct, REGR_4D, OutLev, InLev1, NLEV, T, MET )
	!
	! !INPUT PARAMETERS:
	!
	REAL(sp), POINTER :: REGR_4D(:,:,:,:) ! 4D input data
	INTEGER, INTENT(IN) :: OutLev
	INTEGER, INTENT(IN) :: InLev1
	INTEGER, INTENT(IN) :: NLEV
	INTEGER, INTENT(IN) :: T
	INTEGER, INTENT(IN) :: MET ! 4=GEOS-4, 22=GISS E2.2, else GEOS-5
	!
	! !INPUT/OUTPUT PARAMETERS:
	!
	TYPE(ListCont), POINTER :: Lct ! HEMCO list container
	!
	! !REVISION HISTORY:
	! 30 Dec 2014 - C. Keller - Initial version
	! See https://github.com/geoschem/hemco for complete history
	!EOP
	!------------------------------------------------------------------------------
	!BOC
	INTEGER :: I, NZ, ILEV, TOPLEV
	REAL(hp) :: THICK
	REAL(hp), POINTER :: EDG(:)
	REAL(hp), ALLOCATABLE :: WGT(:)
	!=================================================================
	! COLLAPSE begins here
	!=================================================================
	! Init
	EDG => NULL()
	! Reset
	Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = 0.0_hp
	! Don't do anything if there are not enough levels in REGR_4D
	NZ = SIZE(REGR_4D,3)
	IF ( NZ < InLev1 ) RETURN
	! Get maximum level to be used for pressure thickness calculations.
	TOPLEV = InLev1 + ( NLEV-1 )
	! Get pointer to grid edges on the native input grid
	IF ( Met == 4 ) THEN
	EDG => G4_EDGE_NATIVE(InLev1:TOPLEV)
	ELSE IF ( Met == 22 ) THEN
	EDG => E102_EDGE_NATIVE(InLev1:TOPLEV)
	ELSE
	EDG => G5_EDGE_NATIVE(InLev1:TOPLEV)
	ENDIF
	! Thickness of output level
	THICK = EDG(1) - EDG(NLEV)
	! Get level weights
	ALLOCATE(WGT(NLEV))
	WGT = 0.0
	DO I = 1, NLEV-1
	WGT(I) = ( EDG(I) - EDG(I+1) ) / THICK
	ENDDO
	! Pass levels to output data, one after each other
	Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = REGR_4D(:,:,InLev1,T) * WGT(1)
	DO I = 1, NLEV-1
	ILEV = InLev1 + I
	IF ( NZ < ILEV ) EXIT
	Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) &
	+ ( REGR_4D(:,:,ILEV,T) * WGT(I+1) )
	ENDDO
	! Cleanup
	DEALLOCATE(WGT)
	EDG => NULL()
	END SUBROUTINE COLLAPSE

I can make these simple changes that account for the fact that we are regridding layers, not levels. For example, if InLev1 is equal to layer 37 (and should collapse 2 layers, 37 and 38), then the TOPLEV should be the upper edge of grid box 38 (which is layer 39).

sed -i -e "s|TOPLEV = InLev1 + ( NLEV-1 )|TOPLEV = InLev1 + NLEV|g" \
           -e "s|THICK = EDG(1) - EDG(NLEV)|THICK = EDG(1) - EDG(NLEV+1)|g" \
           -e "s|DO I = 1, NLEV-1|DO I = 1, NLEV|g" hco_interp_mod.F90

Before this, when two layers were supposed to be being collapse, only the lower level was being taken into account (and the first three when four layers were collapse). The reason I create an issue for this instead of a PR is the fact that it seems like you would want a different subroutine depending on if you are regridding layers (n=72 or 47) or levels (n=73 or 48) which I get the impression that HEMCO is supposed to do based on the codes and comments.

[nicholasbalasus]

Tagging @yantosca @jimmielin. I am happy to write a fix for this but would like your thoughts on the last portion.

[nicholasbalasus]

Also thanks to @toddmooring for flagging this!

[eastjames]

Issue #208 deals with the same subroutine and may be related

[yantosca]

Hi @eastjames @toddmooring @nicholasbalasus @jimmielin. I think you are correct, this seems to have been a bug in translation. We used to have the code in transfer_mod.F (going waaaaayyy back) but then this was brought into HEMCO. I would proceed with your fix after confirming that mass is conserved.

[yantosca]

Or it could have always been like that...

[nicholasbalasus]

Thank you, @yantosca. My changes does conserve mass now (but I've only checked this for one methane file). I'm going to go through the code a little more closely and do some more tests first and then will get back to you.

[nicholasbalasus]

Closing - addressed here daf54d0