Skip to content

Stratospheric-adjusted radiative forcing in GC-Classic (supersedes PR #2525)#3211

Open
yantosca wants to merge 29 commits intodev/14.8.0from
feature/strat-adjusted-rf
Open

Stratospheric-adjusted radiative forcing in GC-Classic (supersedes PR #2525)#3211
yantosca wants to merge 29 commits intodev/14.8.0from
feature/strat-adjusted-rf

Conversation

@yantosca
Copy link
Copy Markdown
Contributor

@yantosca yantosca commented Feb 20, 2026

Name and Institution

Name: Bob Yantosca (Harvard) on behalf of Connor Barker (UCL) and Eloise Marais (UCL)

NOTE: We have rebased the code from PR #2525 (by @cbarker211) atop GEOS-Chem 14.7.0. This PR now supersedes #2525 which will be closed.

Describe the update

This update extends the stratospheric adjustment code added by @sdeastham in GCHP (#2010) to GCClassic. The commit here simply modifies the main level script for GCClassic, and allows the dynamical heating rates to be read in via the HEMCO meteorology file. No changes are made to the underlying RRTMG source code beyond what Seb already has added to include stratospheric adjustment. For further details see the original pull request.

To include stratospheric adjustments of temperature relative to the baseline simulation in the radiative fluxes, follow these steps:

  • Run a baseline simulation with no pertubations.

    • Set the following in geoschem_config.yml (rrmtg must be activated):
      • fixed_dyn_heating: true
      • read_dyn_heating: false
      • Activate the DynHeating collection in HISTORY.rc (sample file provided below) - this outputs the dynamical heating rates.
  • Run a pertubation simulation.

    • Create a symbolic link in the run directory called Heating_Baseline to point to the DynHeating output from the baseline simulation:
      • e.g. ln -s /path/to/reference/rundir/OutputDir Heating_Baseline
    • Set the following in geoschem_config.yml (rrmtg must be activated):
      • fixed_dyn_heating: true
      • read_dyn_heating: true
    • Add the following to the HEMCO_Config.rc.gmao_metfields file - this loads the dynamical heating rates:
# --- DynHeating field ---
* DynHeating  ./Heating_Baseline/GEOSChem.DynHeating.$YYYY$MM$DD_$HH00z.nc4 DynHeating 2019/1-12/1-31/0-23 I xyz 1 * - 1 1

HISTORY.txt

As described in #2010, occasionally there will be a warning of "RRTMG FDH routine failed to converge for x of y columns on CPU 0". We find this happens sporadically for our simulations for 1/3312 columns. Failing to converge for >10 columns usually indicates an error somewhere.

The seasonally-evolving fixed dynamical heating approximation has also been extended to GCClassic, but is not yet tested (as stated by Seb in #2010). The fdh_to_toa flag can be added in geoschem_config.yml to extend the stratospheric adjustment calculations to the top of the model.

Expected changes

This code can be tested by simulating a single day with and without an increase of 100 ppmv CO2 (see #2010).

    - 20190701 - 20190702, no spinup, default restart file.
    - 4 simulations:
        - baseline   (co2_ppmv: 390.0, read_dyn_heating: false)
        - without_sa (co2_ppmv: 490.0, read_dyn_heating: false)
        - sa_strat   (co2_ppmv: 490.0, read_dyn_heating: true, fdh_to_toa: false)
        - sa_toa     (co2_ppmv: 490.0, read_dyn_heating: true, fdh_to_toa: true)
    - If successful, you should get the following results.  
        -           Sim  TOA RF / W m-2  Trop CO2 RF / W m-2  TOA RF - Trop RF / W m-2
         0        No SA            0.64                 1.25                     -0.61
         1  SA in strat            1.10                 1.24                     -0.09
         2    SA to TOA            1.16                 1.24                     -0.02
        - The expected changes are an increase in the TOA RF, and a small decrease in the Trop CO2 RF. 

References

Tagging @cbarker211 @eamarais

Katherine R. Travis and others added 8 commits February 11, 2026 11:33
…/10.1021/acsetair.5c00136. Includes SO2 oxidation in aerosol water by H2O2, O3, NO2, HCHO, and TMI-catalzyed O2. Addition of cloud HPLUS diagnostic since cloud pH diagnostic appears broken.

Signed-off-by: Katherine R. Travis <[email protected]>
…/10.1021/acsetair.5c00136. Includes SO2 oxidation in aerosol water by H2O2, O3, NO2, HCHO, and TMI-catalzyed O2. Addition of cloud HPLUS diagnostic since cloud pH diagnostic appears broken.

Signed-off-by: Katherine R. Travis <[email protected]>
…/10.1021/acsetair.5c00136. Includes SO2 oxidation in aerosol water by H2O2, O3, NO2, HCHO, and TMI-catalzyed O2. Addition of cloud HPLUS diagnostic since cloud pH diagnostic appears broken.

Signed-off-by: Katherine R. Travis <[email protected]>
--> Change 'uptake probability' to 'sticking coefficient' for consistency
--> Consolidate production trackers from 13 to 3 (PSO4MP for multiphase sulfate, and PHMSAQ and PHMSMP for HMS aqueous and multiphase production)
--> Improved computational efficiency by removing ELSE blocks.
--> Added _dp to constants
--> Add Ionic strength diagnostic to GCHP HISTORY template

Signed-off-by: Katherine R. Travis <[email protected]>
GeosCore/aerosol_mod.F90
GeosCore/planeflight_mod.F90
Headers/state_chm_mod.F90
KPP/fullchem/fullchem_HetStateFuncs.F90
- Renamed State_Chm%Isorrop* fields to State_Chm%Ate*, where "Ate"
  stands for "aerosol thermoydnamical equilibrium".  ISORROPIA is no
  longer used so this should avoid confusion.

run/GCClassic/HISTORY.rc.templates/HISTORY.rc.aerosol
run/GCClassic/HISTORY.rc.templates/HISTORY.rc.fullchem
run/GCHP/HISTORY.rc.templates/HISTORY.rc.fullchem
- Rename diagnostic fields from "Isorrop*" to "Ate*"

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/gckpp*F90
- Rebuilt with KPP 3.3.0

KPP/fullchem/CHANGELOG_fullchem.md
- Changed "[ TBD ]" to the standard "[Unreleased] - TBD"

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_HetStateFuncs.F90
- Removed id_DST{1..4} variables
- Added id_DSTbin{1..7} variables
- Updated the equation for DUST to use DSTbin{1..7} species, which
  were added in GEOS-Chem 14.7.0.

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_SulfurChemFuncs.F90
- We have removed the scavenging factor of 0.7 from the expression for
  the DUST variable in routine SET_SO2, based on a suggestion by
  Katie Travis (@ktravis213).
- Cosmetic changes (indentation, etc.

Signed-off-by: Bob Yantosca <[email protected]>
@yantosca yantosca added this to the 14.8.0 milestone Feb 20, 2026
@yantosca yantosca self-assigned this Feb 20, 2026
@yantosca yantosca added category: Feature Request New feature or request topic: Diagnostics Related to output diagnostic data topic: Aerosols Related to aerosol species and/or microphysics in GEOS-Chem labels Feb 20, 2026
@yantosca yantosca force-pushed the feature/strat-adjusted-rf branch from 0db327d to 7dc642f Compare February 20, 2026 23:28
@yantosca
Copy link
Copy Markdown
Contributor Author

I ran an integration test on the feature/strat-adjusted-rf branch. All GEOS-Chem Classic integration tests passed:

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

CodeDir   : 7ce5cb9 GEOS-Chem update: Merge PR #3138 (Diagnostic timestep update)
GEOS-Chem : 7dc642f80 Move array alloc/dealloc for strat-adjusted rf to routines in main.F90
HEMCO     : e23c43b HEMCO 3.12.0 release
Cloud-J   : 1dff6fe Merge pull request #31 from geoschem/release/8.0.2
HETP      : 2a99b24 Merge pull request #2 from geoschem/bugfix/initialize_local_variables

Using 24 OpenMP threads
Number of execution tests: 38

Submitted as SLURM job: 61416506
==============================================================================

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

All GEOS-Chem Classic tests were zero-diff w/r/t gcc.14.7.0, except:

  • gc_4x5_merra2_fullchem_APM
    • Parallel error (known issue)

yantosca added 12 commits March 2, 2026 15:12
KPP/custom/custom.eqn
KPP/fullchem/fullchem.eqn
- Updated the version number to 14.8.0
- Added citation Travis et al 2025 to HMS + OH + SO2 rxn

CHANGELOG.md
- Removed leftover line from merge conflict

Signed-off-by: Bob Yantosca <[email protected]>
run/shared/species_database.yml
- Removed several leftover dummy species

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_HetStateFuncs.F90
- Removed extraneous comma at end of line after id_DSTbin3
- Removed the 0.7 wet scavenging factor as per @ktravis213
KPP/custom/custom.eqn
KPP/fullchem/fullchem.eqn
- Added comments denoting additional multiphase sulfate rxns
  as coming from the Travis2025 paper

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_HetStateFuncs.F90
- Added the following species ID flags as arguments to routine
  fullchem_SetState_Het, calling Ind_ from within the subroutine
  (which is computationally inefficient):
  - id_DSTbin{1..7}, id_pFe, id_SO2, id_SO4
- Removed Ind_ from the USE State_Chm_Mod statement in the
  fullchem_SetStateHet routine

GeosCore/fullchem_mod.F90
- Added local variables id_DSTbin{1..7}, id_pFe, id_SO2, id_SO4,
  which are set in Init_FullChem adn passed as arguments to
  fullchem_SetStateHet

Signed-off-by: Bob Yantosca <[email protected]>
run/shared/species_database.yml
- Added metadata for PHMSAQ and PHMSMP dummy species

Signed-off-by: Bob Yantosca <[email protected]>
Headers/state_chm_mod.F90
- Bug fix: Rename the metadata for all former ISORROP* diagnostics
  to ATE*.  This indicates that ISORROPIA is no longer used, but
  we now use the term ATE (for aerosol thermodynamical equilibrium)
  instead.

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_RateLawFuncs.F90
- Removed k_ex1 variable, we only need k
- Renamed gammaKT to gamma
- Separated local varaiables from arguments
- Added comments
- Cosmetic changes (whitespace, indentation)

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_SetStateHet
- Add calls to SafeDiv to prevent div-by-zero conditions
- Do not call ReactoDiff_Corr when l_r is zero, as this will
  get placed into the denominator, which will cause an div-by-zero
- Removed the xradi variable, now use State_Chem%AeroRadi directly
- Removed e+0 from numerical constants
- Cosmetic changes (indentation, comments, line breaks, aesthetics)

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_RateLawFuncs.F90
- Renamed "CRITRH" parameter to "RH_35_PERCENT" for clarity
- Renamed "CRITRH2" parameter to "RH_50_PERDCENT" for clarity

CHANGELOG.md
- Updated accordingy

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_HetStateFuncs.F90
- Changed "hydroxide*3" to "hydroxide**3" in the expression for
  Fe_d_a.  This was a typo.

Signed-off-by: Bob Yantosca <[email protected]>
KPP/fullchem/fullchem_HetStateFuncs.F90
- In routine SetStateHet:
  - Removed a duplicate definition of K0, CR, pKa for SO2
  - Get Henry's law parameters from HENRY_K0 and HENRY_CR arrays
    for species SO2, H2O2, O3, and NO2.

Signed-off-by: Bob Yantosca <[email protected]>
yantosca added 5 commits April 8, 2026 16:15
KPP/fullchem/fullchem_HetStateFuncs.F90
- In routine SetStateHet:
  - Precompute the 1/ 297K term and store it in the Inv_T297 parameter
  - Precompute the 1 / TK term and store it in the Inv_TK parameter
  - Precompute the 298.15 / TK - 1 term and store it in the
    T298_over_TK_m1 variable
  - Use precomputed terms in equations
  - Renamed "Tk" variable to "TK" (all caps) for consistency
  - Cosmetic changes (indentation, comments, line breaks)

Signed-off-by: Bob Yantosca <[email protected]>
KPP/stubs/stub_fullchem_HetStateFuncs.F90
- Added species ID flags for DSTbin{1..7}, SO2, SO4, and pFe to the
  argument list of SetStateHet (to match fullchem_HetStateFuncs.F90)

KPP/fullchem/fullchem_HetStateFuncs.F90
- Removed an extraneous comment

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
GeosCore/apm_driv_mod.F90
- Bug fix: Added the extra output variable "IONIC" to the call
  to MACH_HETP_Main_15Cases.  This was causing APM simulations to
  fail at compile time.

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
GeosCore/apm_driv_mod.F90
- Removed unneeded variables (e.g. SCASI)from OpenMP loop
- Now save ionic strength output from MACH_HETP_Main_15Cases
  to the HETP_IONIC variable
- Zero several private loop variables

NOTE: APM may be broken, it seems to stop with a floating point error.
Ask the APM Team to investigate.

Signed-off-by: Bob Yantosca <[email protected]>
run/shared/download_data.yml
- Added remote & local data paths for a fullchem restart file
  containing the APM species

run/GCClassic/setupForRestarts.sh
run/shared/download_data.yml
- Modified to read the remote & local data paths to the APM restart file

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
@yantosca yantosca marked this pull request as ready for review April 22, 2026 15:05

END SUBROUTINE Print_Dry_Run_Warning
!EOC
#ifdef RRTMG
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason why this code is in main.F90 and not in an RRTMG module? It seems out of place here.

@yantosca
Copy link
Copy Markdown
Contributor Author

@msulprizio It's how it came to us. Maybe this can be moved.

Copy link
Copy Markdown
Contributor

@msulprizio msulprizio left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @yantosca. I really don't think it makes sense to have code specific for RRTMG in main.F90. I would recommend moving routines Init_RRTMG_Strat_Adj_RF and Cleanup_RRTMG_Strat_Adj_RF from Interfaces/GCClassic/main.F90 to GeosCore/rrtmg_rad_transfer_mod.F90. It will just need to be checked that this doesn't conflict with Interfaces/GCHP/gchp_chunk_mod.F90 or else block off with an #ifdef.

yantosca and others added 4 commits April 22, 2026 15:34
This merge brings PR # (Add multiphase sulfate in aerosol water,
following Travis et al 2025, by @ktravis213 and @yantosca) into
the GEOS-Chem 14.8.0 development stream.

PR #3202 implements multiphase (heterogeneous) sulfate chemistry
in KPP. This includes SO2 oxidation by NO2, O3, H2O2, TMI-catalyzed
O2, and HCHO. All oxidants but NO2 have a dependence on ionic
strength calculated in HETP. Diagnostics are included for the uptake
coefficient gamma in StateChm and production of sulfate from each
new pathway. There is a separate pull request for the addition of
the ionic strength diagnostic in HETP.

Signed-off-by: Bob Yantosca <[email protected]>
Interfaces/GCClassic/main.F90
- Consolidate RRTMG variable declarations into a single #ifdef block
- Moved the allocation and deallocation of DT_3D, DT_3D_UPDATE,
  and HR_3D out of the main routine and into routines
  "Init_RRTMG_Strat_Adj_RF" and "Cleanup_RRTMG_Strat_Adj_RF".
  These ensure that the allocation is done during init and the
  deallocation is done at finalization rather than on each timestep.
- Avoid using (:,:,:) in array assignments (it's not needed)

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
@yantosca yantosca force-pushed the feature/strat-adjusted-rf branch from 7dc642f to c82fff3 Compare April 23, 2026 14:54
@yantosca yantosca removed this from the 14.8.0 milestone Apr 28, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

category: Feature Request New feature or request topic: Aerosols Related to aerosol species and/or microphysics in GEOS-Chem topic: Diagnostics Related to output diagnostic data

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants