Wet deposition: Difference between revisions
| Line 65: | Line 65: | ||
=== Resolve very high tracer concentrations in MERRA wet deposition === | === Resolve very high tracer concentrations in MERRA wet deposition === | ||
'''''This fix | '''''This fix will be validated with [[GEOS-Chem_v11-01#v11-01d|1-month benchmark v11-01d]].''''' | ||
'''''Viral Shah wrote:''''' | '''''Viral Shah wrote:''''' | ||
Revision as of 20:20, 21 October 2015
This page describes the current wet deposition scheme used in GEOS-Chem.
Overview
Original algorithm
The Harvard Atmospheric Chemistry Modeling Group developed a wet deposition scheme (including scavenging of soluble tracer in convective updrafts, as well as rainout and washout of soluble tracers) for the GMI model. This scheme was then implemented into GEOS-Chem. Jacob et al [2000] describes the algorithm in full. This scheme is also described in Liu et al [2001].
Further updates
A number of updates have been added since the implementation of the original wet scavenging algorithm. These include:
Allow both washout and rainout when precipitation forms
NOTE: This update was incorporated into GEOS-Chem v9-01-01.
Qiaoqiao Wang wrote:
- When there is new formation of precipitation in lower layer k, rainout will be applied to the whole precipitation area: max(Fk,Fk+1), considering the contribution of precipitation formation overhead. This will overestimate rainout effect when Fk+1 is much larger than Fk. Therefore, we now apply rainout effect to precipitation area Fk and washout effect to the area: max (0, Fk+1-Fk) in the same grid box.
Updates for aerosol scavenging efficiency
This update was tested in the 1-month benchmark simulation v9-01-03e and approved on 02 Feb 2012.
Qiaoqiao Wang wrote:
- The bulk below-cloud scavenging parameterization of Dana and Hales (k = 0.1P, where P is the precipitation rate mm h-1) used in the standard GEOS-Chem model integrates scavenging efficiencies over typical aerosol size distributions. This overestimates scavenging as it does not account for the preferential removal of the very fine and coarse particles over the course of the precipitation event, shifting the aerosol size distribution toward the more scavenging-resistant accumulation mode that accounts for most of aerosol mass. Now we use the below-cloud scavenging coefficients (k = a Pb constructed by Feng (2007, 2009)) integrated over accumulation mode for most aerosols and over coarse mode for coarse dust and sea salt.
--Bob Y. 15:50, 11 January 2011 (EST)
Updates for MERRA met fields
In GEOS-Chem v9-01-01, we have implemented an improved wet deposition scheme which uses the precipitation fields directly from the MERRA reanalysis product.
The wet deposition algorithm for GEOS-5 is relatively unchanged, except that we allow both rainout and washout to form within a grid box simultaneously. The wet deposition code has been partitioned into several new subroutines within wetdep_mod.f in order to allow for better compatibility with the wet deposition scheme used for the MERRA met fields.
--Bob Y. 14:12, 9 March 2011 (EST)
Add scavenging by snow
This update was tested in the 1-month benchmark simulation v9-01-03e and approved on 02 Feb 2012.
Qiaoqiao Wang wrote:
- For in-cloud scavenging by rain droplets, we assume 100% of water-soluble aerosols are rained out. But in the case of snow, only dust and hydrophobic BC are considered to be IN and then could be rained out. Note that HNO3 is also assumed to be rained out by snow as it forms a monolayer in ice crystal. The below-cloud scavenging coefficients are also higher for snow than for rain droplets.
--Bob Y. 14:14, 9 March 2011 (EST)
Impaction scavenging for hydrophobic BC and homogeneous IN removal
These updates were validated with 1-month benchmark simulation v11-01b and 1-year benchmark simulation v11-01b-Run0. This version was approved on 19 Aug 2015.
From Wang et al. [2014]:
- We modify the scheme by (1) scavenging hydrophobic aerosol (hydrophobic BC and OC) in convective updrafts, since this would take place by impaction [Ekman et al., 2004], and (2) scavenging water-soluble aerosol from cold clouds by homogeneous freezing of solution droplets at T < 237 K [Friedman et al., 2011]. We conducted 222Rn-210Pb simulation to test the general model representation of aerosol deposition and found a lifetime of tropospheric 210Pb aerosol against deposition of 8.6 days.
--Melissa Sulprizio 11:13, 26 November 2013 (EST)
Validation
See Liu et al [2001].
Previous issues that have been resolved
Resolve very high tracer concentrations in MERRA wet deposition
This fix will be validated with 1-month benchmark v11-01d.
Viral Shah wrote:
- I am running v10.01, nested-grid NA full-chemistry w/ SOA simulation with the GEOS-FP met fields, and have found instances of extremely high concentrations of certain tracers that develop all of a sudden. I traced this back to a bug in the DO_MERRA_CONVECTION routine at the following:
! Check if the tracer is an aerosol or not
IF ( AER == .TRUE. ) THEN
!---------------------------------------------------------
! Washout of aerosol tracers
! This is modeled as a kinetic process
!---------------------------------------------------------
! Define ALPHA, the fraction of raindrops that
! re-evaporate when falling from (I,J,L+1) to (I,J,L)
ALPHA = ( REEVAPCN(K) * AD(K) )
& / ( PDOWN(K+1) * AREA_M2 * 10e+0_fp )
! ALPHA2 is the fraction of the rained-out aerosols
! that gets resuspended in grid box (I,J,L)
ALPHA2 = 0.5e+0_fp * ALPHA
- Here, the value of ALPHA, which should be less than or equal to 1, is sometimes much higher than 1, and for the attached example was >10^5. Since ALPHA is used to calculate resuspension of precipitating mass from above, these high values of ALPHA lead to a gain in tracer mass at the lower levels in excess of what is coming from above. I suggest we replace the bold font code with the following piece of code which is used elsewhere in the DO_MERRA_CONVECTION routine.
! %%%% CASE 1 %%%%
! Partial re-evaporation. Less precip is leaving
! the grid box then entered from above.
IF ( PDOWN(K+1) > PDOWN(K) .AND.
& PDOWN(K) > TINYNUM ) THEN
! Define ALPHA, the fraction of raindrops that
! re-evaporate when falling from grid box
! (I,J,L+1) to (I,J,L)
ALPHA = ( REEVAPCN(K) * AD(K) )
& / ( PDOWN(K+1) * AREA_M2 * 10e+0_fp )
! For safety
ALPHA = MIN( ALPHA, 1e+0_fp )
! ALPHA2 is the fraction of the rained-out aerosols
! that gets resuspended in grid box (I,J,L)
ALPHA2 = 0.5e+0_fp * ALPHA
ENDIF
! %%%% CASE 2 %%%%
! Total re-evaporation. Precip entered from above,
! but no precip is leaving grid box (ALPHA = 2 so
! that ALPHA2 = 1)
IF ( PDOWN(K) < TINYNUM ) THEN
ALPHA2 = 1e+0_fp
ENDIF
- A second bug I noticed is in the lines right below the first one. The current code is as follows:
! GAINED is the rained out aerosol coming down from
! grid box (I,J,L+1) that will evaporate and re-enter
! the atmosphere in the gas phase in grid box (I,J,L).
GAINED = T0_SUM * ALPHA2
! Amount of aerosol lost to washout in grid box
WETLOSS = Q(K,IC) * BMASS(K) / TCVV *
& WASHFRAC - GAINED
! LOST is the rained out aerosol coming down from
! grid box (I,J,L+1) that will remain in the liquid
! phase in grid box (I,J,L) and will NOT re-evaporate.
LOST = T0_SUM - GAINED
! Update tracer concentration (V. Shah, mps, 5/20/15)
Q(K,IC) = Q(K,IC) - WETLOSS * TCVV / BMASS(K)
- In this case, the amount of tracer GAINED is (mistakenly) not added to Q before WETLOSS is calculated. The fix for the bold faced line should be as follows:
WETLOSS = ( Q(K,IC) * BMASS(K) / TCVV + GAINED ) *
& WASHFRAC - GAINED
- This will be consistent with the code (below) for the 'non-aerosol' case
! MASS_WASH is the total amount of non-aerosol tracer
! that is available for washout in grid box (I,J,L).
! It consists of the mass in the precipitating
! part of box (I,J,L), plus the previously rained-out
! tracer coming down from grid box (I,J,L+1).
! (Eq. 15, Jacob et al, 2000).
MASS_WASH = ( F_WASHOUT * Q(K,IC) ) * BMASS(K) /
& TCVV + T0_SUM
! WETLOSS is the amount of tracer mass in
! grid box (I,J,L) that is lost to washout.
! (Eq. 16, Jacob et al, 2000)
WETLOSS = MASS_WASH * WASHFRAC - T0_SUM
--Melissa Sulprizio (talk) 16:43, 14 September 2015 (UTC)
Bugs in MERRA wet deposition
These fixes were included in the GEOS-Chem v10-01 public release.
Carey Friedman wrote:
- I think I've found two bugs in the MERRA wet deposition code that potentially affect all MERRA users but have particular relevance for Hg/POPs users.
- In the POPs model, which I am currently developing (i.e., not the public version), wet deposition mass is passed to depo_pops_mod.F via the routines ADD_POPG_WD and ADD_POPP_WD in convection_mod.F and wetscav_mod.F. This is analogous to what happens in the mercury model: wet dep mass is passed to depo_mercury_mod.F via ADD_Hg2_WD and ADD_HgP_WD. For POPs (and I suspect for Hg, since I literally copied the Hg code, though I haven't tested Hg explicitly), the mass that is passed to the depo module is less than the mass that is archived in the wet dep diagnostics (AD38 and AD39) when MERRA is used. They are equal using GEOS5. For Hg, I think this means mass is lost going from the air to the ocean.
- I've found two areas in the MERRA wet dep code that seem to be causing mass to be lost. The first, in convection_mod.F, is in the DO_MERRA_CONVECTION subroutine for washout of non-aerosol tracers below the cloud base. The current code looks like this:
! MASS_NOWASH is the amount of non-aerosol tracer in
! grid box (I,J,L) that is NOT available for washout.
MASS_NOWASH = ( 1d0 - F_WASHOUT ) * Q(K,IC)
! MASS_WASH is the total amount of non-aerosol tracer
! that is available for washout in grid box (I,J,L).
! It consists of the mass in the precipitating
! part of box (I,J,L), plus the previously rained-out
! tracer coming down from grid box (I,J,L+1).
! (Eq. 15, Jacob et al, 2000).
MASS_WASH = ( F_WASHOUT * Q(K,IC) ) + T0_SUM
! WETLOSS is the amount of tracer mass in
! grid box (I,J,L) that is lost to washout.
! (Eq. 16, Jacob et al, 2000)
WETLOSS = MASS_WASH * WASHFRAC - T0_SUM
- When I change the last equation to
WETLOSS = (MASS_WASH - T0_SUM)* WASHFRAC
- then the mass passed to depo_pops_mod.F and the AD38 diagnostic become equal. This seems to make sense physically based on the comments, but I'm hoping someone more familiar with the convection process can confirm this. If this indeed fixes a bug, then there is analogous code for washout of non-aerosol tracers in DO_WASHOUT_ONLY in wetscav_mod.F that needs to be changed too.
- The second, in wetscav_mod.F, is in the DO_WASHOUT_ONLY subroutine. The following lines, when commented, allow the mass passed to depo_pops_mod.F and AD39 to be equal:
#if defined( MERRA ) || defined( GEOS_57 )
!%%% BUG FIX: Only apply this section of when running
!%%% GEOS-Chem with MERRA met fields. (bmy, hamos, 5/26/11)
!
! Check if WASHFRAC = NaN or WASHFRAC < 0.1 %
!
! If WASHFRAC = NaN, then DSTT = NaN and SAFETY trips because
! tracer concentrations must be finite.WASHFRAC = NaN when F = 0.
! When less than 0.1% of a soluble tracer is available for washout
! DSTT < 0 and SAFETY trips. (Helen Amos, 20100928)
IF ( IT_IS_NAN( WASHFRAC ) ) THEN
CYCLE
ELSEIF ( WASHFRAC < 1D-3 ) THEN
CYCLE
ENDIF
#endif
- If I change this code to
IF ( IT_IS_NAN( WASHFRAC ) ) THEN
CYCLE
ELSEIF ( WASHFRAC < 1D-20 ) THEN
CYCLE
ENDIF
- rather than commenting it, then this reduces the amount of mass that is lost, but the % lost slowly increases over time such that is becomes significant over the course of several years. This is a problem for POPs, some of which run for >100 years.
- Incidentally, POPs also passes dry dep to depo_pops_mod.F the same way Hg does and there are no problems there going between GEOS5 and MERRA.
Helen Amos wrote:
- Remembering back 4-5 years ago when we did this... the fix is only needed for MERRA because the problem was never encountered with GEOS-5 during testing. We never figured out why MERRA met fields had the problem and GEOS-5 didn't, other than the problem seemed to be on the data assimilation side. 1D-3 was selected as a lower limit based on a series of tests on what trapped negative DSTT values and what didn't. This isn't pretty or exact, so if you come up with a better solution, I definitely welcome it!
- DSTT was tested with the Hg and Pb-Be-Rn simulations.
Carey Friedman wrote:
- I would recommend further investigation, since as Helen said, at one point 1d-3 was deemed the lowest cutoff before you started tripping SAFETY. So I don't want to recommend 1d-20, or to comment the code (even though it seems to work for POPs) only to cause people's runs to crash.
- My best recommendation would be to first fix the multiplication by WASHFRAC issue in DO_MERRA_CONVECTION (convection_mod.F) and DO_WASHOUT_ONLY (wetscav_mod.F) and then re-test Hg and Pb-Be-Rn to see if the trap to keep SAFETY from tripping is still needed at all. The WASHFRAC multiplication bugs will make WETLOSS less than it should be, and that could be contributing to the DSTT < 0 problem.
--Melissa Sulprizio 15:37, 2 February 2015 (EST)
Viral Shah wrote:
- I have been also examining these bugs in the wetscav and convection routines, and looking at how they affect the Hg simulation, in particular. My findings/recommendations are as follows:
- 1) For the WASHOUT_ONLY routine (in WETSCAV_MOD):
- The condition requiring that the WASHFRAC>1D-3 or any number greater than 0 for washout to occur, prevents the partial resuspension of the dissolved mass falling from above, and leads to an overestimate in the wet deposited mass. This is consistent with Shaojie's results of an increase in the atmospheric mass of Hg0 following this fix. This also leads to a lower Hg2 mass in the lower troposphere. I recommend that the condition on WASHFRAC should be removed.
- 2) For the DO_MERRA_CONVECTION routine (CONVECTION_MOD):
- I think, the bug here is very similar to the first bug. A condition on F(K,IC)>0 prevents the archiving of deposited mass in the wet deposition diagnostic, but it is recorded in the mass passed from the air to the ocean. After comparing with Jacob et al (2000), I think the current code for calculating washout follows the documented algorithm well, and does not seem to be causing the bug. It would be great to see if this fixes the issue in Carey's POP simulation.
- 3) In the course of reviewing the code, I found two other bugs in DO_MERRA_CONVECTION. First, the code is missing a statement that subtracts the wet deposited amount from the atmospheric mass. And, second, there seems to be an inconsistency in units with how the variable T0_SUM was used.
- I have attached a document describing this is more detail, and the fixes I propose for these bugs. I am also attaching the following link to the benchmark plots, comparing the results with the original code and the code with my proposed fixes.
--Melissa Sulprizio (talk) 20:30, 20 May 2015 (UTC)
Bug fixes for scavenging by co-condensation
These updates were validated with the 1-month benchmark simulation v10-01f and approved on Approved 13 Jan 2015.
Duncan Fairlie wrote:
- One detail that has got me puzzled concerns scavenging by co-condensation:
- Formulation under Mari et al. (2000) (for H2O2) for fraction in solid phase looks like
- Whereas in the code wetdep scheme it looks like
- Where W is the cloud LWC. I.e. The square root term is inverted in the 2 formulations.
- The code changes apply to subroutine COMPUTE_F, and concern the parameters CONV_H2O2 and CONV_NH3. To be consistent with the formulation in Mari et al. (2000), CONV_H2O2 for example should be defined as
- And similarly for CONV_NH3. I.e. The square-root term should be inverted w.r.t. what is currently in the code.
- Replace lines:
! CONV_H2O2 = 0.6 * SQRT( 1.9 ), used for the ice to gas ratio for H2O2
! 0.6 is ( sticking coeff H2O2 / sticking coeff water )
! 1.9 is ( molecular weight H2O2 / molecular weight water )
REAL*8, PARAMETER :: CONV_H2O2 = 8.27042925126d-1
! CONV_NH3 = 0.6 * SQRT( 0.9 ), used for the ice to gas ratio for NH3
! 0.6 is ( sticking coeff NH3 / sticking coeff water )
! 0.9 is ( molecular weight NH3 / molecular weight water )
REAL*8, PARAMETER :: CONV_NH3 = 5.69209978831d-1
- With
! CONV_H2O2 = 0.6 * SQRT( 0.53 ), used for the ice to gas ratio for H2O2
! 0.6 is ( sticking coeff H2O2 / sticking coeff water )
! 0.53 is ( molecular weight water / molecular weight H2O2 )
REAL*8, PARAMETER :: CONV_H2O2 = 4.36564d-1
! CONV_NH3 = 0.6 * SQRT( 1.06 ), used for the ice to gas ratio for NH3
! 0.6 is ( sticking coeff NH3 / sticking coeff water )
! 1.06 is ( molecular weight water / molecular weight NH3 )
REAL*8, PARAMETER :: CONV_NH3 = 6.17395d-1
- That should correct the issue. These changes should also apply to the RAINOUT subroutine, in addition to COMPUTE_F.
- This also applies to Equ. 9 in the documentation of wet scavenging: Jacob et al., 2000, Harvard wet deposition scheme for GMI. The term sqrt(MH2O2/MH2O) in Equ. 9 should read sqrt(MH2O/MH2O2). With associated changes to the text following.
--Melissa Sulprizio 14:56, 19 December 2014 (EST)
Prevent div-by-zero in WASHFRAC routines
These fixes were tested with 1-month benchmark GEOS-Chem v9-02p and approved on 13 Sep 2013.
In the following routines in GeosCore/wetscav_mod.F:
- WASHFRAC_FINE_AEROSOL
- WASHFRAC_COARSE_AEROSOL
- WASHFRAC_SIZE_AEROSOL
- WASHFRAC_HNO3
we now make sure to only compute the WASHFRAC quantity when F > 0. In other words, wherever we see an expression of the form:
WASHFRAC = F *(1d0 - EXP( -K_WASH *
& (PP / F*3.6d4 )**0.61d0 * DT / 3.6d3 ))
we now bracket this with an IF statement, such as this:
IF ( F > 0d0 ) THEN
WASHFRAC = F *(1d0 - EXP( -K_WASH *
& (PP / F*3.6d4 )**0.61d0 * DT / 3.6d3 ))
ELSE
WASHFRAC = 0d0
ENDIF
Because F multiplies the entire expression, WASHFRAC would become zero whenever F=0. But the IF statement also avoids placing F=0 into the denominator of a division, which would result in a div-by-zero condition.
This error condition was caught by the GEOS-Chem unit tester when compiling with option FPE=yes.
--Bob Y. 12:19, 12 September 2013 (EDT)
Bug in Henry's constant
This update was tested in the 1-month benchmark simulation v9-02l and approved on 26 Jun 2013. This update is included in Adjoint v35.
Fabien Paulot wrote:
- I think there is a small error in the definition of the Henry's constant, which is used in the wet deposition routine
- The universal gas constant is defined as
R = 8.32d-2 L atm K-1 mol-1
- But I think it should be 8.314/101 or
R = 8.2057d-2 L atm K-1 mol-1
- This error (if error there is) seems to go back a very long way (Jacob 2000)..
- This is small but this will decrease the ratio of liquid to gas by 1.6%
--Melissa Payer 17:44, 15 May 2013 (EDT)
Anomalous PRECTOT values on 22 July 2010
This update was tested in the 1-month benchmark simulation v9-01-03g and approved on 27 Feb 2012.
Several users have reported their 2° x 2.5° GEOS-5 simulations dying in the wet deposition module on 22 July 2010. We have traced this to an error in the GEOS-5 met fields. The PREACC variable in routine MAKE_QQ of wetscav_mod.F (which holds the GEOS-5 met field PRECTOT) contains the total precipitation at the ground. This should always be larger than the PRECON variable (which holds the GEOS-5 met field PRECCON), which is the convective fraction.
However, on at 06 GMT on 22 July 2010, the PREACC array contains some anomalous values (of the order 1e-20). This results in a situation where PREACC is larger than PRECON. Because the fraction of the precipitating box (FRAC) computed according to:
FRAC = ( PREACC(I,J) - PRECON(I,J) ) / PREACC(I,J)
then, this can lead to an (unphysical) large negative value for FRAC.
We recommend the following kludge. Add the following code in routine MAKE_QQ:
!==============================================================
! %%%%% FOR ALL OTHER MET FIELDS EXCEPT MERRA %%%%%
!==============================================================
IF ( PREACC(I,J) > 0d0 ) THEN
! Large scale or convective fraction of precipitation
IF ( LS ) THEN
FRAC = ( PREACC(I,J) - PRECON(I,J) ) / PREACC(I,J)
ELSE
FRAC = PRECON(I,J) / PREACC(I,J)
ENDIF
!##################################################################
!### KLUDGE: On July 22, 2010, there is an error where PREACC
!### is set to a very small number (but nonzero). This causes
!### FRAC to blow up to near infinity and PDOWN to be a very
!### large number. For now, limit FRAC to be between 0 and 1.
!###
IF ( FRAC > 1d0 ) FRAC = 1d0
IF ( FRAC < 0d0 ) FRAC = 0d0
!##################################################################
This will always make sure that the precipitating fraction of the grid box is between 0 and 1. Here is an example:
---> DATE: 2010/07/22 GMT: 06:00 X-HRS: 6.000
- Found all 33 A-3 met fields for 2010/07/22 07:30
- Found all 5 I-6 met fields for 2010/07/22 12:00
- LINOZ_CHEM3: Doing LINOZ stratospheric chemistry
### FRAC(99,75): -6.807457425009224E+016
### RESET FRAC(99,75): 0.000000000000000E+000
NOTE: Simulations done with the MERRA met fields are not affected because we use a different wet scavenging algorithm.
--Bob Y. 17:06, 1 December 2011 (EST)
Washout fix for non-aerosol species
A detailed description of this fix can be found in the Appendix of Amos et al. (2012, ACP).
This update was tested in the 1-month benchmark simulation v9-01-02f and approved on 02 Aug 2011.
Prior to GEOS-Chem v9-01-02, this bug was causing the wet deposition algorithm to underestimate washout of highly soluble gases:
Helen Amos wrote:
- In section 2.2.2 of Jacob et al. (2000), the first bullet point states:
If fi,L > Fi/f, scavenging is limited by mass transfer and we follow exactly the same procedure as for aerosols and HNO3 (section 2.2.1).
- Prior to GEOS-Chem v9-01-02, WASHFRAC_LIQ_GAS was not executing washout of non-aerosol tracers as dictated by this statement.
- Washout proceeds in two ways: (1) Henry's law equilibrium and (2) limited by mass transfer. Aerosols, HNO3, and highly soluble gases should proceed through washout obeying the mass transfer limitation. Moderately and sparingly soluble gases should proceed through washout obeying Henry's law equilibrium. The purpose of WASHFRAC_LIQ_GAS is to calculate WASHFRAC for non-aerosol species, where WASHFRAC is the fraction of a tracer available for washout. WASHFRAC is calculated using Henry's law if fi,L <= Fi/f (Eqn 15-16, Jacob et al. 2000), and WASHFRAC is calculated using mass transfer if fi,L > Fi/f (Eqn 14, Jacob et al. 2000). The algorithms in WASHOUT are written such that if you specify that a tracer proceed through washout according to Henry's law, the code expects WASHFRAC to be calculated using Henry's law (Eqn 15-16, Jacob et al. 2000). Similarly, if you specify that a tracer proceed through washout limited by mass transfer, then the code expects WASHFRAC to be calculated using mass transfer (Eqn 14). Prior to GEOS-Chem v9-01-02, every soluble gas (except HNO3) was arbitrarily forced to proceed through washout according to Henry's law, even when fi,L > Fi/f and WASHFRAC had been calculated from mass transfer equations. In short, there was a mismatch in information between WASHFRAC_LIQ_GAS and WASHOUT. This mismatch caused an underestimate in the washout of highly soluble gases because too much tracer was being re-evaporated (negative term in WETLOSS, where WETLOSS is the amount of tracer lost in a grid box).
- A logical is now passed as an input/output parameter to WASHFRAC_LIQ_GAS and this logical will specify whether a non-aerosol species should proceed through washout according to Henry's law equilibrium or limited by mass transfer. The code is now consistent with Section 2.2.2 of Jacob et al. (2000).
This fix was originally tested in the 1-month benchmark simulation v9-01-02b. A second bug in WASHFRAC_LIQ_GAS was discovered, invalidating v9-01-02b (see discussion below). The complete fix for WASHFRAC_LIQ_GAS has been added to the standard code and approved in the 1-month benchmark simulation v9-01-02f.
Update 7/28/2011: The original algorithm we tested in 1-month benchmark v9-01-02b was incorrect. Prior to v9-01-02b, routine WASHFRAC_LIQ_GAS looked like this:
IF (WASHFRAC > WASHFRAC_F_14) THEN WASHFRAC = WASHFRAC_F_14
In benchmark simulation v9-01-02b we had modified this IF statement (incorrectly) to:
IF ( WASHFRAC > WASHFRAC_F_14 ) THEN WASHFRAC = WASHFRAC_F_14 AER = TRUE. ENDIF
However, this IF statement should have been:
IF ( WASHFRAC > WASHFRAC_F_14 ) THEN WASHFRAC = F * WASHFRAC_F_14 AER = TRUE. ENDIF
This fix has been added to the standard code and was approved in the v9-01-02f benchmark. For more information, please refer to the Appendix of Amos et al. (2012, ACP).
--Helen Amos 03:57, 13 Aug 2011 (EST)
Negative tracer in routine WETDEP because of negative RH
See this post: GEOS-5 issues#Small negative RH value in 20060206.a6.2x25 file
Fixes are available at ftp://ftp.as.harvard.edu/pub/geos-chem/patches/v8-01-01.
--phs 16:31, 6 June 2008 (EDT)
Negative tracer in routine WETDEP
Dylan Millet wrote:
- I'm having a run die consistently at the same time (October 1, 2005; first time step of the month) in large-scale wetdep, with an STT element < 0.
- Platform: Linux cluster
- Threads: 8
- Version: v7-4-13 out of the box.
- GEOS4, 4x5, 30L, full chemistry
- IFORT 10.1
- In Section 6 (No Downward Precip) of wetscav_mod.f, subroutine safety is getting called.
WETDEP - STT < 0 at 1 1 29 for tracer 7 in area 6
- (First of all it seems odd to do wetdep for L=29, this is 63 km up). Have you seen anything like this? I ran for the whole year starting Jan 1 successfully until this point.
- ... By the way, the problem persists when I turn off chemistry altogether.
Philippe Le Sager replied:
- I used your restart file and the same input.geos (w/ chemistry on and off). My code went thru without problem. I tried both Sun Studio and Ifort 9 compilers, and the later on two different machines (altix and ceres). I used v7-04-13 and v8-01-01. I never reproduced your error.
- We just got the new Ifort 10, and tried it too. I run v8-01-01 without an error. But when I tried v7-04-13, I finally reproduced your error, with the exact same negative values!
- In other words: the bug happens with IFort 10 and v7-04-13 only.
- Also, have a look at this recent development. This is not the reason for your bug (I tried v8 w/ ifort 10 and isorropia -like v7-04-13- and it did not crash), but using RPMARES instead of Isorropia may be a way to fix it.
- ... More about the Ifort 10 / v7-04-13 issue. When I wanted to debug with TotalView, I could not reproduce the bug anymore.... because I simply suppress any optimization. So, I did more test and found that if the default -O2 optimization is used, GEOS-Chem crashes. But it works fine with -O1. It is hard to tell what happens, since only the emissions step is done between reading the restart file and the crash.
- Bob and I will further test Ifort 10 for optimization on our machines. Maybe we will find something... For the time being, you may have to switch to -O1, at least for the run that crashes. You will find the optimization flag at the beginning of the Makefile.ifort.
Long story short: This appears to be an optimization issue with IFORT 10 and v7-04-13. Upgrading to GEOS-Chem v8-01-01 should solve this problem.
--Bmy 10:38, 17 April 2008 (EDT)
Negative tracer in routine WETDEP #2
If you find find negative values caused by the wet deposition (particularly on data date 22 July 2010), then here is a workaround:
Mark Parrington wrote
- I'm getting negative values in the wet deposition and am having some difficulty in identifying a possible reason. I've experienced this problem in 3 different versions of the model (v8-02-01, v8-02-04, and v8-03-02) and have tried turning off the chemistry and changing the optimization. Has anyone else had a similar problem?
Claire Carouge replied:
- I don't know if this would help, but a crash on negative or NaN values in wet deposition doesn't always mean that there is a problem in wet deposition.
- GEOS-Chem has the subroutine CHECK_STT in GeosCore/tracer_mod.f. This function is used at some points in the code to check the concentration values (STT array) and stop the code if there is any negative, NaN, or Inf values. In particular, it's used in the wet deposition, but not before this point.
- I think the first thing to do is to identify were the bad values are created. You can insert calls to the CHECK_STT subroutine at different places in main.f (e.g. after transport, after PBL, after emissions, etc.) and see where the code stops. Then you can refine by adding the calls to CHECK_STT into the problematic subroutine and so on.
Mark Parrington followed up:
- I traced the bug to excessive evaporation in a subroutine called MAKE_QQ in wetscav_mod.f - which gives very small, ~1.0E-22, values of the accumulated precipitation (PREACC) leading to very large negative values of the convective fraction of precipitation (FRAC). To get around this I hard-wired the value of FRAC in MAKE_QQ to be 0.67 for the large scale fraction or else 0.33, if the value of PREACC is less than 1.0E-10 (around line number 474 in wetscav_mod.f v9-01-01):
! Large scale or convective fraction of precipitation
IF ( LS ) THEN
FRAC = ( PREACC(I,J) - PRECON(I,J) ) / PREACC(I,J)
! mp hack for -ve wet deposition
if ( preacc(i,j) .lt. 1.0E-10 ) frac = 0.67
ELSE
FRAC = PRECON(I,J) / PREACC(I,J)
! mp hack for -ve wet deposition
if ( preacc(i,j) .lt. 1.0E-10 ) frac = 0.33
ENDIF
- I was successfully able to run the model up to the end of October this way. I found this issue in v8-02-04 and v9-01-01 but only for the 2x2.5 simulation.
--Bob Y. 12:07, 28 April 2011 (EDT)
Outstanding issues
Low tropospheric 210Pb lifetime against deposition in v11-01b
Hongyu Liu wrote:
- The tropospheric 210Pb lifetime against deposition of 6.6 days (G-C/GEOS-FP) is the shortest among G-C versions. Switching from GEOS-5 to GEOS-FP met fields (v9-02r) reduced 210Pb lifetime from 9.3 to 7.7 days. See: http://wiki.seas.harvard.edu/geos-chem/index.php/Rn-Pb-Be_simulation#Budget_of_Pb210. If v11-01b had been driven by GEOS-5, the 210Pb lifetime may very well be ~8 days. Wang, Q. et al. [2014] used G-C/GEOS-5 and "find a lifetime of tropospheric 210Pb aerosol against deposition of 8.6 days, as compared to a best estimate of 9 days constrained by observations [Liu et al., 2001]." So it appears that the short 210Pb lifetime (6.6 days) in v11-01b is largely due to switching to GEOS-FP. Wet deposition needs to be re-evaluated and tuned for GEOS-FP.
- Bo Zhang has done a series of G-C v11-01b simulations of 210Pb with GEOS-FP (2013 - year for benchmark) and GEOS-5 (2009 - year used in Wang et al. 2014)...We were able to reproduce the 210Pb lifetime (GEOS-5, 2009) reported in Wang Q. et al. (2014) using the standard v11-01b code, where scavenging due to cold cloud is not applied to 210Pb. Switching from GEOS-5 to FP reduces 210Pb lifetime by about 2 days. (But note that G-C/GEOS5 uses MOISTQ to approximate 3-D precipitation formation rates while G-C/FP takes the precipitation formation rates from the FP archive.) Bo is currently evaluating these simulations with 210Pb surface concentration / flux & UT/LS observations.
Bo Zhang and Hongyu Liu provided the following quick fix recommendation on 10/8/15 that will be implemented in v11-01d:
- Recommendation for a quick fix: For GEOS-FP (and MERRA2), turn off cold cloud scavenging & change the condensed water content (CWC) back to 1.5x10-3kg/m3. The resulting 210Pb aerosol RT against wetdep (7.8 days) would be close to that in GEOS-5 (with cold cloud scaveging & CWC=1x10-3kg/m3, Q. Wang et al. 2014).
--Lizzie Lundgren (talk) 18:23, 13 October 2015 (UTC)
Values of F_PRIME greater than 1 in WETDEP
In GEOS-Chem v8-03-01, the value of the variable F_PRIME may end up being slightly greater than 1 under some circumstances. This can cause a negative tracer condition to occur in subroutine WETDEP of wetscav_mod.f.
Kelley Wells wrote:
- The following error message is continually produced in the 03 to 06 UTC time period of 29 July 2010 (2° x 2.5°, GEOS-5).
WETDEP - STT < 0 at 73 65 11 for tracer 1 in area 4
- We've determined that the issue is due to the variable F_PRIME being slightly greater than 1.0 at this point, thus resulting in F > 1.0 and finally RAINFRAC > 1.0. But, it's currently a mystery as to why F_PRIME is greater than 1.0, since both variables that go into its calculation (Q and K_RAIN) have positive values.
Dylan Millet wrote:
- The value of F_PRIME is only very marginally >1, so that
print*,f_prime gives 1.00000000
- but
print*,f_prime - 1d0 gives something on the order of 1E-12
- As best we can tell this is what gives the negative tracer value which causes the run to crash. Inserting a line of code setting:
F = MAX( F_PRIME, FTOP ) !### Kludge, don't let F be greater than 1 F = MIN( F, 1d0 )
- Seems to have fixed the issue. Perhaps we can be more careful and use a more elaborate IF test such as:
IF( ( ( F - 1d0 ) > 0d0 ) ) .and. ( ( F - 1d0 ) < 1d-6 ) ) THEN
F = 1d0
ENDIF
- which would catch this type of numerical instability error but not risk masking some other problem.
This type of error can be attributed to numerical instability. You can test by recompiling with the DEBUG=yes option, which will turn off all optimization. You might find that this error disappears when the optimization is turned off. However, turning off all optimization is only recommended for debugging purposes, as code that is not optimized will run much more slowly. Therefore, we recommend adding an IF statement (as shown above) to make sure that the value of F never exceeds 1.
NOTE: In GEOS-Chem v9-01-01 the wet deposition code was rewritten such that subroutine WETDEP was broken up into smaller subroutines, for better compatibility with the new MERRA wet scavenging algorithm. We are unsure if this condition exists in v9-01-01 as of this writing.
--Bob Y. 09:24, 11 March 2011 (EST)
References
- Amos, H. M., D. J. Jacob, C. D. Holmes, J. A. Fisher, Q. Wang, R. M. Yantosca, E. S. Corbitt, E. Galarneau, A. P. Rutter, M. S. Gustin, A. Steffen, J. J. Schauer, J. A. Graydon, V. L. St. Louis, R. W. Talbot, E. S. Edgerton, Y. Zhang, and E. M. Sunderland, Gas-Particle Partitioning of Atmopsheric Hg(II) and Its Effect on Global Mercury Deposition, Atmos. Chem. Phys.,12,591-603, 2012. PDF
- Dana, M.T., and J.M. Hales, Statistical aspects of the washout of polydisperse aerosols, Atmos. Environ, 10, 45-50, 1976.
- Domine, F., and E. Thibert, Mechanism of incorporation of trace gases in ice grown from the gas phase, Geophys. Res. Lett., 23, 3627-3630, 1996.
- Giorgi, F., and W.L. Chameides, Rainout lifetimes of highly soluble aerosols as inferred from simulations with a general circulation model, J. Geophys. Res., 91, 14,367-14,376, 1986.
- Jacob, D.J., Heterogeneous chemistry and tropospheric ozone, Atmos. Environ., 34, 2131-2159, 2000. PDF
- Jacob, D.J. H. Liu, C.Mari, and R.M. Yantosca, Harvard wet deposition scheme for GMI, Harvard University Atmospheric Chemistry Modeling Group, revised March 2000. PDF
- Levine, S.Z., and S.E. Schwartz, In-cloud and below-cloud scavenging of nitric acid vapor, Atmos. Environ., 16, 1725-1734, 1982.
- Liu, H., D.J. Jacob, I. Bey, and R.M. Yantosca, Constraints from 210Pb and 7Be on wet deposition and transport in a global three-dimensional chemical tracer model driven by assimilated meteorological fields, J. Geophys. Res., 106, 12,109-12,128, 2001. PDF
- Mari, C., D.J. Jacob, and P. Bechtold, Transport and scavenging of soluble gases in a deep convective cloud, J. Geophys. Res., 105, 22,255-22,267, 2000. PDF
- Selin, N.E. and D.J. Jacob, Seasonal and spatial patterns of mercury wet deposition in the United States: North American vs. intercontinental sources, Atmospheric Environment, 42, 5193-5204, 2008. PDF
- Wang, Q., D.J. Jacob, J.A. Fisher, J. Mao, E.M. Leibensperger, C.C. Carouge, P. Le Sager, Y. Kondo, J.L. Jimenez, M.J. Cubison, and S.J. Doherty, Sources of carbonaceous aerosols and deposited black carbon in the Arctic in winter-spring: implications for radiative forcing, Atmos. Chem. Phys., 11, 12,453-12,473, 2011. PDF
- Wang, Q., D.J. Jacob,J.R Spackman, A.E. Perring, J.P. Schwarz, N. Moteki, E.A. Marais, C. Ge, J. Wang and S.R.H. Barrett, Global budget and radiative forcing of black carbon aerosol: constraints from pole-to-pole (HIPPO) observations across the Pacific, J. Geophys. Res., 119, 195-206, 2014. PDF