Fixing simplified chemistry - #10
Conversation
| 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 ) ) |
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
Hmm.. It seems that I have some different version... I need to check still...
|
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. |
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.