[BUG/ISSUE] Typo in CloudHet #906
Comments
|
Thanks Viral. Also tagging @msl3v on this. We can try to get this into 13.3.0 and benchmark. |
yantosca
added
the
topic: Chemical Mechanisms
This update addresses the issue from geoschem/geos-chem #906. The MAX( ff, 1.0e30 ) expression in CloudHet, CloudHet1R, and CloudHet2R functions should be MIN( ff, 1.0e30 ). Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
|
This has been fixed with commit c6adf4b. Awaiting PR review and merge. |
|
This has now been merged into the 13.3.0 development branch and tagged with 13.3.0-alpha.4. |
|
Hi @viral211 and @cdholmes. I merged this PR as 13.3.0-alpha.4 but @msulprizio and I noticed that it caused fullchem sims to choke on the 2nd chemistry timestep with a could not integrate error: I am wondering if that MIN value is even needed. In my new hetrates code I now use a SafeDiv function to set the result at 1e30 if the division cannot be done. So we may be able to get away with this: ! Ratio of volume inside to outside cloud
! ff has a range [0,+inf], so cap it at 1e30
ff = SafeDiv( H%CldFr, H%ClearFr, 1.0e+30_dp ) |
|
Also, this was not caught in integration tests because the integration tests only do 1 chemistry timestep, and stopped before the timestep of the error. We should probably extend integration tests to 1hr of simulation runtime. |
|
Omitting the the MIN line, as you have done, should be fine. Does that fix the integration error? |
|
It does not fix the integration error, I just treid. The original code in the old gckpp_HetRates.F90 was: geos-chem/KPP/fullchem/gckpp_HetRates.F90 Lines 1276 to 1294 in ba58da4 If you print out some of the values of ff, xx, and kHet, there are sometimes negative xx and negative kHet: I would think that you wouldn't want to have negative reaction rates. Should these be flushed to zero? |
|
@cdholmes, flushing negative xx to zero gets us past the integration error: !------------------------------------------------------------------------
! Grid-average loss frequency
!
! EXACT expression for entrainment-limited uptake
!------------------------------------------------------------------------
! Ratio (in cloud) of heterogeneous loss to detrainment, s/s
kk = kI * tauc
! Ratio of volume inside to outside cloud
! ff has a range [0,+inf], so cap it at 1e30
ff = SafeDiv( H%CldFr, H%ClearFr, 1.0e+30_dp )
!ff = MAX( ff, 1.0e+30_dp )
ff = MIN( ff, 1.0e+30_dp )
! Ratio of mass inside to outside cloud
! xx has range [0,+inf], but ff is capped at 1e30, so this shouldn't overflow
xx = ( ff - kk - 1.0_dp ) / 2.0_dp + &
SQRT( 1e0_dp + ff**2 + kk**2 + 2*ff**2 + 2*kk**2 - 2*ff*kk ) / 2.0_dp
! Prevent negative xx and negative kHet (bmy, 9/27/21)
xx = MAX( xx, 0.0_dp )
! Overall heterogeneous loss rate, grid average, 1/s
! kHet = kI * xx / ( 1d0 + xx )
! Since the expression ( xx / (1+xx) ) may behave badly when xx>>1,
! use the equivalent 1 / (1 + 1/x) with an upper bound on 1/x
kHet = kI / ( 1.0_dp + SafeDiv( 1.0_dp, xx, 1.0e+30_dp ) )
! Overall loss rate in a particular reaction branch, 1/s
kHet = kHet * branch |
|
|
Hi @yantosca. The calculation of
I wonder if that's he reason for negative values of xx. |
|
@cdholmes, here is printout from a few locations where xx and khet are zero: I have a larger log file (75 MB) if you'd liek to see the whole thing. |
|
@viral211 thanks. I was sure that I had fixed that in the new code but I guess it fell thru the cracks. Sorry for the bother. |
|
It looks like the bug in the xx definition was causing the negative values. I checked a few examples from the log excerpt that Bob sent. Can you confirm that the negatives are gone when xx is defined correctly? |
|
I was able to run a 1-day fullchem_benchmark run and I didn't get the integration error. With the negatives the run bombed out before 00:30. So I think we are OK now. |
|
I ran an internal 1-month benchmark to evaluate the impact of these changes (tagged as Changes are considerably large, as suspected by Viral in the original post above. I will close out this issue for now since this will officially be benchmarked in 13.3.0, but @viral211, @cdholmes, or others, please feel free to comment here on any of the changes to model output. |
|
Thanks @msulprizio. I was expecting a general increase in ozone because NOx loss on clouds is now much slower. The decrease in ozone over the southern ocean was unexpected. It is caused by the increase in Br, but I am not sure what drives it. I spoke with Xuan Wang and we suspect that it is because the loss of Br through the HOBr+S(IV) reaction in clouds is now slower (along with other HOBr cloud reactions) and more HOBr is available to react with HBr/Br-/Cl- on aerosols and release radicals. Xuan is looking into this further. |
There is a typo in CloudHet in gckpp_Hetrates.F90. The line
geos-chem/KPP/fullchem/gckpp_HetRates.F90
Line 1282 in ba58da4
should be:
ff = min( ff, 1.0e30_dp )I have run it by Chris Holmes @cdholmes. I think this would change the benchmark output significantly.
The text was updated successfully, but these errors were encountered: