Skip to content

Fixing simplified chemistry - #10

Draft
noije wants to merge 4 commits into
plesager:add-hamm7from
noije:simchem_fixes
Draft

Fixing simplified chemistry#10
noije wants to merge 4 commits into
plesager:add-hamm7from
noije:simchem_fixes

Conversation

@noije

@noije noije commented Aug 24, 2026

Copy link
Copy Markdown
Collaborator

Fixes in simplified chemistry, described in #7 and #9. A switch has been introduced to allow for delayed mixing of gas-phase SO2 between the cloudy and clear-sky fractions.

@noije
noije requested a review from plesager August 26, 2026 08:30
ZFAQ_SO2 = 1._JPRD / (1._JPRD + 1._JPRD / (ZHCC_SO2_EFF * ZPNEB ) )
! -- Eq. (7.8) in Seinfeld and Pandis (Third Edition), i.e. f_A / (1 + f_A)
! -- bug fix: ZPNEB in the original code by Remy and Bock replaced by ZCLW_VFRAC (=w_L in S&P)
ZFAQ_H2O2 = 1._JPRD / (1._JPRD + 1._JPRD / (ZHCC_H2O2 * ZCLW_VFRAC ) )

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I had a crash on this line caused by ZCLW_VFRAC going to 0 at some points. I think it is because when calculating the ZCLW_VFRAC:

ZCLW_VFRAC = 0._JPRB
IF (ZPNEB > 1.0E-12_JPRB) THEN
        ZCLW_VFRAC = PQLI(JL,JK) / ZPNEB * ZAIR_DENS / ZRHOLW
ENDIF

The liquid water content can be 0 even though cloud fraction is small. Thus I modified the IF to be:
IF (ZPNEB > 1.0E-12_JPRB .AND. PQLI(JL,JK) > 1.0E-20_JPRB) THEN
And also the below where there is the NOT:

      !IF (.NOT.(ZPNEB > 1.0E-12_JPRB )) THEN  ! 1.0E-12 is the default PNEB value
      IF (ZPNEB <= 1.0E-12_JPRB .OR. PQLI(JL,JK) <= 1.0E-20_JPRB ) THEN  ! 1.0E-12 is the default PNEB value

and it fixed the crash

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Thanks, Eemeli. I am a bit surprised because in my version the ZFAQ_* fractions are not calculated
IF ( ZCLW_VFRAC < 1.0E-10_JPRD ) THEN

Possibly, the product ZHCC * ZCLW_VFRAC is getting too small, but this can only happen if the dimensionless Henry coefficient gets small. At 298 K, the values are ~2 10^6 for H2O2, ~0.25 for O3 and ~4 10^5 for SO2 (effective value at pH=5). Then, there is a temperature dependence which can reduce the values, but this doesn't change the values by orders of magnitude (see, e.g., #7 for SO2).

Still, if division by a small number is the problem this can simply be solved by rewriting the fractions:

         ZFAQ_H2O2 =  ZHCC_H2O2 * ZCLW_VFRAC / ( 1._JPRD + ZHCC_H2O2  * ZCLW_VFRAC  )
         ZFAQ_O3   =  ZHCC_O3 * ZCLW_VFRAC / ( 1._JPRD + ZHCC_O3 * ZCLW_VFRAC  )
         ZFAQ_SO2  =  ZHCC_SO2_EFF * ZCLW_VFRAC /  ( 1._JPRD  + ZHCC_SO2_EFF * ZCLW_VFRAC  )

PS: I don't see anything like
!IF (.NOT.(ZPNEB > 1.0E-12_JPRB )) THEN ! 1.0E-12 is the default PNEB value
in my code. Am I missing something?

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Hmm.. It seems that I have some different version... I need to check still...

@noije
noije marked this pull request as draft August 27, 2026 11:16
@noije

noije commented Aug 28, 2026

Copy link
Copy Markdown
Collaborator Author

I am testing the new chemistry routine in EC-Earth4. After recompiling the whole code, a simulation I started yesterday has completed over 10 years and is still running.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants