APMCTBC1/2 first-hour blow-up followed by stall in carbon_mod.F90

[Eva-bie]

Your name

Shujun

Your affiliation

Fudan University

What happened? What did you expect to happen?

I would like to report a reproducible issue in GCClassic 14.7.0 + APM. The problem appears very early in the simulation: two APM carbon species, APMCTBC1 and APMCTBC2, increase to unphysical values within the first simulation hour, from xx to xx, and the run later stalls around 2023-01-01 09:20 UTC before being killed by Slurm walltime. The final stack trace points to carbon_mod.F90.
Quantitative evidence from restart files

I checked the following restart files:

GEOSChem.Restart.20230101_0000z.nc4
SpeciesRst_APMCTBC1
min = 0.000000E+00
max = 2.095061E-10
sum = 1.029591E-07
SpeciesRst_APMCTBC2
min = 0.000000E+00
max = 4.280331E-12
sum = 6.545060E-09
GEOSChem.Restart.20230101_0100z.nc4
SpeciesRst_APMCTBC1
min = 5.492030E-32
max = 4.379038E+03
sum = 7.004929E+03
SpeciesRst_APMCTBC2
min = 2.669176E-48
max = 1.583638E+00
sum = 3.528099E+00
GEOSChem.Restart.20230101_0200z.nc4
SpeciesRst_APMCTBC1
max = 7.964653E+19
sum = 8.603694E+19
SpeciesRst_APMCTBC2
max = 2.880347E+16
sum = 3.150073E+16
GEOSChem.Restart.20230101_0300z.nc4
SpeciesRst_APMCTBC1
max = 8.373333E+34
sum = 1.072179E+35
SpeciesRst_APMCTBC2
max = 2.218090E+31
sum = 3.053532E+31

This suggests the problem begins during the first simulation hour (0000z -> 0100z).
Localized blow-up pattern

The maximum occurs in a small local region rather than globally. Around one problematic grid cell in the 0100z restart, I found:

SpeciesRst_APMCTBC1
[[9.67887276e-01 1.64448247e+01 8.31764840e-03]
 [1.81063254e+02 4.37903800e+03 1.79416608e-01]
 [4.11908349e-07 4.01659124e-03 1.40821048e-04]]
SpeciesRst_APMCTBC2
[[3.49970580e-04 5.94637398e-03 2.92674972e-06]
 [6.54791568e-02 1.58363788e+00 9.65649295e-05]
 [3.11566143e-10 2.38175896e-06 1.00794996e-07]]

This looks like a localized source/update problem rather than a domain-wide drift.

I also checked nearby meteorological fields in the restart and they do not look obviously pathological:

Met_DELPDRY is around ~15 hPa
Met_BXHEIGHT is around ~110–115 m

So this does not seem to be caused by an extremely small box mass or layer thickness.

The model later slows down dramatically around 2023-01-01 09:20 UTC and is eventually killed by Slurm. The final output is:

---> DATE: 2023/01/01  UTC: 09:20
srun: Job step aborted: Waiting up to 32 seconds for job step to finish.
slurmstepd: error: *** STEP ... CANCELLED ...
forrtl: error (78): process killed (SIGTERM)

gcclassic  ... carbon_mod_mp_zero... line 2559  carbon_mod.F90
gcclassic  ... carbon_mod_mp_soa ... line 2078  carbon_mod.F90
gcclassic  ... carbon_mod_mp_soa ... line 1812  carbon_mod.F90
gcclassic  ... carbon_mod_mp_chem... line 909   carbon_mod.F90
gcclassic  ... chemistry_mod...     line 576    chemistry_mod.F90
gcclassic  ... MAIN__               line 1529   main.F90

This makes me suspect that the early blow-up of APMCTBC1/APMCTBC2 is the upstream cause, and the later stall in carbon_mod.F90 is a downstream consequence.

I am currently testing whether turning off default ship emissions and/or default anthropogenic emissions changes this behavior. I can report those results later in the same issue, but I wanted to submit the core reproducible behavior first because the first-hour blow-up is already clear.

What are the steps to reproduce the bug?

Model: GEOS-Chem Classic 14.7.0
Mechanism: fullchem + APM
Meteorology: MERRA2
Grid: 4° × 5°
Run start time: 2023-01-01 00:00 UTC
Emissions: default model emissions only
Restart: using the restart file https
://drive.google.com/file/d/1ifKomOKCDgfsDcLdGD8cMOEpjpWqMkwA/view ?usp=drive_link

GC.log

HISTORY.txt

geoschem_config.yml
HEMCO_Config.txt

Please attach any relevant configuration and log files.

No response

What GEOS-Chem version were you using?

14.4.0

What environment were you running GEOS-Chem on?

Local cluster

What compiler and version were you using?

gcc-7.3.1

Will you be addressing this bug yourself?

No

In what configuration were you running GEOS-Chem?

GCClassic

What simulation were you running?

Full chemistry

As what resolution were you running GEOS-Chem?

2*2.5

What meterology fields did you use?

MERRA-2

Additional information

No response

[Eva-bie]

Small correction: although I selected 2° x 2.5° in the issue template dropdown by mistake, the actual failing case reported here is a 4° x 5° MERRA2 run.

[yantosca]

Tagging @ganluoasrc

[yantosca]

Thanks for writing @Eva-bie. I noticed you are using an old version of GCC (7.3.1). We have found that older versions of GCC sometimes can lead to numerical instability. We recommend users to use GCC 10 or higher just to be safe.

[Eva-bie]

Hi Bob, thanks for your suggestion!

I did a further code-level comparison and found a likely issue in the APM-coupled BC redistribution logic in GeosCore/carbon_mod.F90.

In v14.7.0 and v14.4.1 source, the relevant BC section originally behaved as if MASS(L) was not representing the total BC-bin mass consistently during redistribution to CTBC.

MASS(L)=Spc(N)%Conc(I,J,L)

when looping over BC bins. This seemed inconsistent with the later redistribution term:

OLD(L,K) * TC0(L) / MASS(L)

because in that form MASS(L) would not represent the total BC-bin mass at level L.

After checking the older v12.6.1 source, I found that the earlier implementation explicitly used the total BC-bin mass at level L:

MASS(L) = SUM(Spc(I,J,L,APMIDS%id_BCBIN1:IDTEMP))

This makes the redistribution term physically and mathematically more consistent, since TC0(L)/MASS(L) then represents the fractional contribution of the current BC bin to the total BC mass in that layer.

Based on this comparison, I tested a provisional modification in v14.7.0 so that MASS(L) is accumulated over all BC bins before being used in the redistribution. After recompiling and making sure the new executable was actually used, the abnormal blow-up of APMCTBC1/2 disappeared. The restart outputs from 0000z to 0300z are now stable and remain within reasonable magnitudes, instead of jumping to unrealistically large values in the first hour.

At the same time, I also noticed that the v14.4.1 implementation appears to be the same as v14.7.0 in this section. So I would like to ask whether, from the intended scientific/algorithmic assumption, MASS(L) here should indeed represent the total BC-bin mass at each level during redistribution.

In parallel, I am also updating the compiler toolchain as suggested, and that work is ongoing.

[yantosca]

Thanks @Eva-bie.

At the same time, I also noticed that the v14.4.1 implementation appears to be the same as v14.7.0 in this section. So I would like to ask whether, from the intended scientific/algorithmic assumption, MASS(L) here should indeed represent the total BC-bin mass at each level during redistribution.

That is a better question for @ganluoasrc, who is our APM developer.

[ganluoasrc]

Based on our coating assumption, it is MASS(L) = SUM(Spc(I,J,L,APMIDS%id_BCBIN1:IDTEMP)) as shown in earlier version. We use it in our group.