[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

• [x] Support guidelines
• [x] User manuals
• [x] Debugging GEOS-Chem and HEMCO errors
• [x] 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. https://github.com/geoschem/HEMCO/blob/477c7e820c253c24e61d7c53df1911cd0e896abf/src/Core/hco_interp_mod.F90#L1220-L1311

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 https://github.com/geoschem/HEMCO/commit/daf54d07ed16d0948302070855468b6b27fdd9a5