ISORROPIA II

Overview

Brief description

The ISORROPIA II package performs aerosol thermodynamical equilibrium. It partitions nitrate (HNO3 and NIT) and ammonia (NH3 and NH4) between the gas and aerosol phases. Inputs to the partitioning routine include temperature and RH. ISORROPIA II has significant benefits over previous implementations of ISORROPIA, especially for partitioning of nitrate at low RH.

From Pye et al [2009]:

ISORROPIA II [Fountoukis and Nenes, 2007] is implemented in GEOS-Chem to compute gas-aerosol equilibrium partitioning of nitric acid and ammonia. Particles in this study are not size-resolved; however, they can be generally assumed to represent PM2.5 since formation of sulfate-nitrate-ammonium on coarse mode sea salt and dust is excluded. Submicrometer-sized particles are likely to reach gas-aerosol equilibrium on time-scales less than the 1 hour computational time step used here [Meng and Seinfeld, 1996].

Sodium and chloride from accumulation mode sea salt are considered in the gas-aerosol equilibrium along with sulfate, nitrate, and ammonium. Calcium, magnesium, and potassium concentrations are not considered in the present study because of the issues with dust emissions previously mentioned. All inorganic aerosols are assumed to exist on the upper, metastable branch of the hygroscopic hysterisis curve. Although this assumption may not hold at higher altitudes in the free troposphere [Wang et al., 2008], since the focus of this study is mainly on surface-level concentrations, where humidities reach high values on a daily basis, the metastable assumption is acceptable.

Authors and collaborators

• Thanos Nenes (Georgia Tech) -- Principal Investigator
• Havala O. T. Pye (Caltech)

--Bob Y. 12:48, 2 March 2010 (EST)

Implementation notes

ISORROPIA II has been implemented into GEOS-Chem v8-03-01 and later versions.

Code structure

ISORROPIA II consists of the following files:

Files in ISOROPIA/ subdirectory (prior to GEOS-Chem v11-02 (aka 12.0.0)
-----------------------------------------------------------------------
Makefile                 Makefile for ISORROPIA II code
isoropiaIIcode.F         Source code file with ISORROPIA II routines
isrpia.inc               ISOROPIA II header file with common blocks

Files in ISOROPIA/ subdirectory (in GEOS-Chem v11-02 (aka 12.0.0) and later versions)
-------------------------------------------------------------------------------------
Makefile                 Makefile for ISORROPIA II code
isoropiaII_main_mod.F90  Source code file with ISORROPIA II routines

Files in GeosCore/ subdirectory:
--------------------------------
isoropiaII_mod.f         "Interface" code between GEOS-Chem and ISORROPIA II

--Bob Yantosca (talk) 20:14, 6 June 2018 (UTC)

Additional Documentation

Havala O. T. Pye wrote:

1. Documentation (including a user manual) for ISORROPIAII can be found on the ISORROPIA website.
2. isoropiaIIcode.f is essentially ISOFWD.FOR and ISOCOM.FOR of the ISORROPIA II box model pasted together.
3. For more information on ISORROPIA II, see
   • Fountoukis, C. and Nenes, A.: ISORROPIA II: a computationally efficient thermodynamic equilibrium model for K+–Ca2+–Mg2+–NH4+–Na+–SO42−–NO3−–Cl−–H2O aerosols, Atmos. Chem. Phys., 7, 4639-4659, 2007. pdf
4. The implementation by Pye et al. 2009 JGR did not include Ca, K, or Mg since dust emissions were not used.

--Bob Y. 12:28, 17 March 2010 (EDT)

Validation

See Pye et al [2009]

References

1. Fountoukis, C., and A. Nenes (2007), ISORROPIA II: A computationally efficient thermodynamic equilibrium model for K+-Ca²⁺-Mg²⁺-NH₄⁺-Na⁺-SO₄^2-NO₃^--Cl-H₂O aerosols, Atmos. Chem. Phys., 7(17), 4639-4659.
2. Meng, Z. Y., and J. H. Seinfeld, Time scales to achieve atmospheric gas-aerosol equilibrium for volatile species, Atmos. Environ., 30(16), 28892900, doi:10.1016/1352-2310(95)00493-9., 1996.
3. Pye, H. O. T., H. Liao, S. Wu, L. J. Mickley, D. J.Jacob, D. K. Henze, and J. H. Seinfeld, Effect of changes in climate and emissions on future sulfate-nitrate-ammonium aerosol levels in the United States, J. Geophys. Res., 114, D01205, 2009. PDF
4. Wang, J., A. A. Hoffman, R. J. Park, D. Jacob, and S. T. Martin, Global distribution of solid and aqueous sulfate aerosols: Effect of the hysteresis of particle phase transitions, J. Geophys. Res., 113, D11206, doi:10.1029/2007JD009367, 2008. PDF

--Bob Y. 16:00, 22 February 2010 (EST)

Previous issues now resolved

Investigating persistent noise observed in ISORROPIA output

ISORROPIA was updated to version 2.2 in GEOS-Chem 12.3.0, which was released on 01 Apr 2019.

Seb Eastham wrote:

I’ve done some minor detective work at my end. Short version: the “noise” is likely the result of extremely non-linear behavior in certain areas of the parameter domain. This non-linear behavior is at least partially due to incomplete convergence of the solution rather than true non-linearity in the physical system. This analysis does not explore the implications of the finding for global aerosols, but my initial experiments suggest that this occurs at low sulfur concentrations and at all temperatures and humidities. Full explanation follows below.

Experiments:

I developed a stand-alone implementation of ISORROPIA to look at the partitioning behavior under conditions typical to the upper troposphere, where we seem to get the most noise. The plot below, which I call ISORROPIA Analysis, shows the results of parameter sweeps with ISORROPIA under the following conditions:

• T = 255 K
• RH = 30%
• Total NH3 (NH4 + NH3) fixed at 30 nmol/m3
• Total SO4 varying from 0 to 1 nmol/m3
• Total NO3 varying from 3 to 6 nmol/m3

This is obviously a contrived set of experiments, but it shows what I suspect is the root of this problematic behavior. The sweep was performed in 4 environments:

1. Top-left: Using the current GEOS-Chem implementation of ISORROPIA v2.0, by overriding the inputs to ISORROPIA in GEOS-Chem. The pixelation in this plot is the result of using a much smaller number of samples, and is not related to ISORROPIA. This plot is purely to show that the ‘in-GEOS-Chem’ implementation of ISORROPIA gives identical results to the standalone implementation (used for all other plots), and will not be mentioned again.
2. Top-right: Results from standalone implementation of ISORROPIA (details in the final block of text, below). This exactly reproduced the GEOS-Chem behavior.
3. Bottom left: The same standalone implementation, upgraded to ISORROPIA v2.2 using Shannon’s fixes.
4. Bottom right: The original implementation, but with a much lower tolerance and higher iteration count when calculating activities (EPSACT increased from 0.05 to 10^-6, NSWEEP increased from 4 to 20).

The main takeaways are:

1. The solution is extremely non-linear in this region. An increase of 0.01 nmol/m3 NO3 could result in +0.02 nmol/m3 NH3 if SO4 = 0, or -0.03 nmol/m3 NH3 if SO4 = 0.45
2. The solution is not fully converged under the standard GEOS-Chem settings for ISORROPIA, resulting in additional non-linearity not present in the fully converged solution. For example, there is an entire domain of solutions in the standard implementation which does not exist when the tolerance is reduced (note the region around SO4 = 0.4 nmol/m3). This also significantly and artificially increases the number of distinct regions in which the gradient of (in this case gaseous NH3 concentration) varies with a minute change in concentration, such as might occur through a photochemical response to changes in ozone or AOD elsewhere in the column.
3. Even with the enhanced convergence criteria, some “noisy” regions remain. However, the change in the solution pattern suggests that this non-linearity is not entirely physical.

To further explore this issue, I’ve created the second plot shown below (ISORROPIA Analysis Single), which shows broader non-convergence issues with increased precursor concentrations and higher temperatures. The conditions here are:

• T = 280 K
• RH = 30%
• Total SO4 fixed at 30 nmol/m3
• Total NH3 varying from 20 to 40 nmol/m3
• Total NO3 varying from 20 to 40 nmol/m3

The left hand plot shows the results with standard convergence criteria (EPSACT = 0.05, NSWEEP = 4), compared to results for enhanced criteria (EPSACT = 10^-6, NSWEEP = 20) on the right. These plots are with ISORROPIA v2.0. As you can see, the solution becomes extremely non-linear in the (TNH3 > 30 and TNO3 > 30) region. The improved convergence criteria do not result in a smoother solution, just a different one. This plot shows that the sign of the change in gas-phase NH3 could flip depending on a small difference in background conditions.

Finally, the full extent of the issue is revealed by estimating the sensitivity of gas-phase NH3 with respect to total NO3 (ISORROPIA Sensitivity), calculated using a slice of the results from the ISORROPIA Analysis Single dataset.

This shows the rate of change of the gas-phase NH3 concentration with respect to a change in total NO3, plotted as a function of this total. If two adjacent grid cells contain ~30 nmol/m3 of NO3 but initially differ in NO3 loading by 1 nmol/m3, they could have sensitivities of opposite sign.

A final note: editing the code to increase the accuracy of the root finding method can be done in other ways, including some of the v2.2 fixes and an increase in the initial number of subdivisions. Some combinations of these approaches were found to mitigate some of the nonlinearity in this region, but always with limited success. If these non-linearities are purely or largely the result of incomplete solution convergence, a dedicated investigation into enhancing or replacing certain core numerical methods in ISORROPIA may help.

Conclusions

1. The noise observed here is likely the result of non-linearity of the partitioning solution generated by ISORROPIA.
2. This non-linearity can result in grid cells with nearly identical conditions having opposite sensitivities with respect to a small change in precursor concentration.
3. Decreasing the solution tolerance and increasing the number of iterations mitigates the issue but does not solve it.
4. A dedicated investigation into the underlying numerical techniques may yield a solution.

Methods

See this section for more information about the ISORROPIA software updates.

--Bob Yantosca (talk) 20:03, 25 January 2017 (UTC)

Bug fixes for ISORROPIA II stable mode

These two fixes were included in v11-02e (approved 24 Mar 2018).

Jingyuan Shao (UW) wrote:

In isoropiaII_mod.F (line 558) when calculate aerosol PH by isorropia, the value would be too high when we choose stable mode. After I checked the code, I think code should be changed as follows :

        IF ( AERLIQ(8) < 1e-32_fp ) THEN
           ! Aerosol is dry so HPLUSTEMP and PH_SAV are undefined
           ! We force HPLUSTEMP to 1d20 and PH_SAV to -999e+0_fp.
           ! (hotp, ccc, 12/18/09)
           HPLUSTEMP       = 1e+20_fp
           PH_SAV(I,J,L)   = -999e+0_fp
        ELSE
           HPLUSTEMP       = AERLIQ(1) / AERLIQ(8) * 1e+3_fp/18e+0_fp

to

        IF ( AERLIQ(8) < 1e-18_fp ) THEN
           ! Aerosol is dry so HPLUSTEMP and PH_SAV are undefined
           ! We force HPLUSTEMP to 1d20 and PH_SAV to -999e+0_fp.
           ! (hotp, ccc, 12/18/09)
           HPLUSTEMP       = 1e+20_fp
           PH_SAV(I,J,L)   = -999e+0_fp
        ELSE
           HPLUSTEMP       = AERLIQ(1) / AERLIQ(8) * 1e+3_fp/18e+0_fp

or

        IF ( OTHER(1) .eq. 1.d0 ) THEN
           ! Aerosol is dry so HPLUSTEMP and PH_SAV are undefined
           ! We force HPLUSTEMP to 1d20 and PH_SAV to -999e+0_fp.
           ! (hotp, ccc, 12/18/09)
           HPLUSTEMP       = 1e+20_fp
           PH_SAV(I,J,L)   = -999e+0_fp
        ELSE
           HPLUSTEMP       = AERLIQ(1) / AERLIQ(8) * 1e+3_fp/18e+0_fp

the default is metastable mode, there is no influence if we choose metastable mode because metastable means supersaturation state, water always exist on aerosol if we choose metastable mode.

Shaojie Song (Harvard) wrote:

We have identified and fixed coding errors in ISORROPIA which significantly affect aerosol water pH calculations when the stable state for aerosol phase is applied. These errors exist in source codes for the cases G2 and O2, and also affect the pH results for the cases G1 and O1. The standard ISORROPIA codes fail to account for the partitioning of ammonia between the gas and aqueous phases. More information can be seen in our recent ACPD paper: https://www.atmos-chem-phys-discuss.net/acp-2018-6/.

In order to resolve this issue, in isoropiaIIcode.F of GEOS-Chem v11-01, both line 20401 and line 25829 (which consist of the same codes):

     PSI4 = MIN(PSI5+PSI6,CHI4)

should be replaced with

     IF(W(2).GT.TINY) THEN       ! Accounts for NH3 evaporation
        BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
        CC   = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4
        DD   = MAX(BB*BB-4.d0*CC,ZERO)  ! Patch proposed by Uma Shankar, 19/11/01
        PSI4 =0.5d0*(-BB - SQRT(DD))
     ELSE
        PSI4 = TINY
     ENDIF

The same coding errors probably exist in all of the previous versions of GEOS-Chem because they exist in the versions from 1.5 to 2.1 of ISORROPIA. Below is an example showing the average aerosol pH by GEOS-Chem v10-01 during Jan 1st-Jan 15th, 2015 (left: standard model; right: revised model). A big difference can be found over some land areas such as USA, Europe, and East Asia, which is mainly because their relatively higher concentrations of ammonia. The marked value represents the pH of Beijing, which is reduced from about 7 (neutral) to about 4 (moderately acidic). Note that there is no influence when using metastable mode (the default).

--Melissa Sulprizio (talk) 18:48, 30 January 2018 (UTC)

Implement ISORROPIA v2.0 as a Fortran module

This update was included in v11-02a and approved on 12 May 2017.

Seb Eastham has repackaged the existing ISORROPIA v2.0 code as a Fortran module. This finally allows us to remove files ISOROPIA/isoropiaII_code.F and ISOROPIA/isrpia.inc, which have proven to be problematic for some newer compilers (e.g. GNU Fortran v5 and v6). He writes:

1. I’ve adapted ISORROPIA v2.0 into a Fortran module, converting COMMON blocks and BLOCK DATA statements into module-level arrays. This can safely replace the current implementation in GEOS-Chem.
   • NOTE: This update should now avoid the THREADPRIVATE error in GNU Fortran v5 and v6.
2. All REAL variables are declared using a custom KIND parameter. Implicit variables are forbidden (using IMPLICIT NONE at the module level) and all variables are declared explicitly. This allows precision to be specified easily.
3. The number of OpenMP THREADPRIVATE variables has also been significantly reduced to minimize overhead and improve code readability.
4. The module can be switched between the legacy version (v2.0) and the updated version (v2.2) by (un)commenting the line #define ISORROPIA_V22 at the top of the module file.
5. When compiled as ISORROPIA v2.0, the file exactly reproduces results produced with GEOS-Chem (zero differences w/r/t prior model versions).
6. I’ve put together a wrapper program which can take the isorropiaII_main_mod.F file exactly as written, compile it into a standalone program, run a parameter sweep and store the output into a NetCDF file. Contact me if this might be of use to you.

Update 21 Feb 2017: The ISORROPIA module now allows you to compile GEOS-Chem with GNU Fortran v4.9 and higher.

Update 26 May 2017: We have brought the ISORROPIA module update into the GEOS-Chem v11-01 public release code, in the v11-01-ISORROPIA-Gfortran branch. If you would like to use Gfortran to compile the GEOS-Chem v11-01 public-release code, we recommend that you check out the code as follows:

 git clone https://bitbucket.org/gcst/geos-chem Code.v11-01
 cd Code.v11-01
 git checkout v11-01-ISORROPIA-Gfortran

--Bob Yantosca (talk) 18:30, 26 May 2017 (UTC)

Optimization and/or parallelization issues in ISORROPIA II

During the benchmarking process for GEOS-Chem v9-01-02, we discovered that ISORROPIA does not produce exactly the same output on single processor compared to multi-processor:

Daven Henze wrote:

There is an unidentified OpenMP error which contributes to differences shown [in the v9-01-02a 1-month benchmark], as confirmed by independent tests of consecutive evaluations of a single executable compiled with and without parallelization. The impacts of this OpenMP error are typically less than a fraction of a percent, and negligible in terms of absolute value, with the exception of NIT and NH3. This is an issue to investigate further.

Bob Yantosca wrote:

Here is a quick update on the issue in ISORROPIA.

1. If you compile with –O0 (all optimization turned off), you get identical results when using ISORROPIA as single vs. multi processor.
2. When you compile with –O1, then you get different results when using ISORROPIA as single vs. multi processor. This inclines me to believe that the ISORROPIA code is sensitive to optimization.
3. It appears that issue is occurring at the computation of the GNH3, GHNO3, and GHCL variables in function FUNCG5A within file isoropiaIIcode.f. These variables are computed as:

     GNH3      = MAX(CHI4 - PSI4, TINY)              ! Gas NH3
     GHNO3     = MAX(CHI5 - PSI5, TINY)              ! Gas HNO3
     GHCL      = MAX(CHI6 - PSI6, TINY)              ! Gas HCl

Here the CHI* and PSI* are common-block variables and TINY is a small number (1d-20, but I used 1d-28 for testing). The CHI* and PSI* values are held threadprivate for the OpenMP parallelization. It looks as if there is some kind of roundoff error at the limit of precision that is creeping in here that could be caused by the optimization. For example, at grid box (37,21,11) every once and a while we get this roundoff error:

   with single processor
   ### CHI5 :   4.326986976485327E-008
   ### PSI5 :   4.312588684524357E-008
   ### GHNO3:   1.439829196097068E-010

   with multi processor
   ### CHI5 :   4.326986976485327E-008
   ### PSI5 :   4.312588684524356E-008
   ### GHNO3:   1.439829196097134E-010

with similar roundoff errors (not shown here) happening in the computation of GNH3.

GHNO3 gets saved back into the tracer array STT(I,J,L,IDTHNO3) and GNH3 gets saved back into STT(I,J,L,IDTNH3). So that is why these we see these very small (and seemingly random) differences in HNO3, NH3, and NIT, whereas the other tracers are unaffected.

I don’t have a good sense of why precisely this is occurring, other than to say that ISORROPIA uses a lot of internal common blocks and other obsolete coding practices. The ISORROPIA code is extremely confusing to follow. My guess is that either something is not being held private for the parallelization, or that variables in common blocks are not being initialized properly.

Solution

Shannon Capps wrote:

Please see this attached report that documents my replication of the issue with both the GEOS-Chem and ANISORROPIA versions of ISORROPIA II in a box model....I would propose adding a flag (-fp-model source) for the compilation of ISORROPIA only. Add this to the file ISOROPIA/Makefile as follows:

   isoropiaIIcode.o:  isoropiaIIcode.F  isrpia.inc
           $(F90) $(R8) -fp-model source -c $<

Compiling the GEOS-Chem implementation of ISORROPIA with this flag causes the box model to produce results that "diff" finds identical in all combinations (e.g., m0 v. s0, m3 v. s0, etc.).

The -fp-source option of the Intel Fortran Compiler tells the compiler to only do "safe" optimizations (i.e. nothing that would affect the numerical precision of the output). For some reason, ISORROPIA II is very sensitive to the compiler optimization.

Bob Yantosca has done some additional 3-day test runs to assess the impact of using the -fp-model source option. The results are:

Run #	Machine	# CPUs	Wall time [mm:ss]	Scalability	Mean OH [1e5 molec/cm3]
1	terra-03-c.as.harvard.edu	4	18:10	3.5935	13.0384629976933
2	terra-04-c.as.harvard.edu	4	09:58	3.8123	13.0384629976933
3	terra-03-c.as.harvard.edu	4	16:38	3.7546	13.0384629976933
4	terra-04-c.as.harvard.edu	4	11:00	3.5090	13.0384629976933
5	terra-11-c.as.harvard.edu	8	09:29	7.3364	13.0384629976933
6	terra-06-c.as.harvard.edu	8	06:18	7.1421	13.0384629976933
7	terra-11-c.as.harvard.edu	8	09:10	7.5530	13.0384629976933
8	terra-01-c.as.harvard.edu	8	05:52	7.4549	13.0384629976933

All of the output files from runs 1-8 were 100% binary identical to each other. Therefore, using -fp-model source seems to eliminate the random numerical noise generated by ISORROPIA.

Persistence of problem

Update 8/26/11: In GEOS-Chem v9-01-02 we shall apply the -fp-model source option to all source code files, not just ISORROPIA. This will help us to ensure the precision of the results as we make widespread changes for ESMF compatibility.

Update 9/1/11: Oops, this may not be solved after all. Runs longer than 3 days still manifest the problem. See this PDF presentation for more details.

--Bob Y. 15:32, 23 November 2011 (EST)

Minor fixes in isorropiaIIcode.f

Havala Pye wrote

There are 2 minor bug fixes for ISORROPIA II. There are two locations in file ISOROPIA/isorropiaIIcode.f where the variable MOLALR(4) should be MOLALR(9). These fixes affect how aerosol phase water is calculated. Since GEOS-Chem typically considers sodium and chloride, these subcases are likely not being used, but the code should be updated in case these subcases are used in the future.

This fix will be implemented in GEOS-Chem v9-01-02.

--Bob Y. 11:10, 1 March 2011 (EST)

Patches for v8-03-01

The following issues in ISORROPIA II have now been fixed as of August 2010:

1. Bug fix in ISORROPIA for offline aerosol
2. Patch to fix declaration of IONIC in ISORROPIA II
3. Bug fix for ND42 diagnostic and ISORROPIA II

--Bob Y. 14:44, 3 August 2010 (EDT)

Outstanding issues

Chemistry instability introduced by ISORROPIA v2.2

Katie Travis found unstable behavior in simulating remote HNO3 against ATom observations over the remote ocean in GEOS-Chem 12.3.0 and later. Seb Eastham suggests the following:

ISORROPIA is not recommended outside of a pretty narrow temperature and pressure range. There’s an optional flag in isorropiaII_mod.F called skip_if_p_and_t_are_out_of_range, which prevents ISORROPIA from running in grid boxes where P < 200 hPa or T < 250 K. This may not be relevant for your simulations, but it can help to reduce noise at higher altitudes. You could try enabling that (uncomment the line define SKIP_IF_P… in isorropiaII_mod.F, around line 108) and see if it helps. This should allow more reliable changes at lower altitude to propagate to higher altitude in a smoother fashion.

--Melissa Sulprizio (talk) 17:37, 14 January 2020 (UTC)

Parallelization issue for nested grid simulations

For those of you who are running one of the GEOS-Chem nested grid simulations, you may encounter a segmentation fault error. This error is caused by GEOS-Chem running up against the memory limits of the machine. One way to solve this problem is to compile ISORROPIA II for single processor only:

Sungshik Patrick Kim wrote:

I tracked down the memory bottlenecks in my [China nested-grid simulation]....The crashes then were caused by memory allocation in ISORROPIA (with all of the private variables that had to be copied), so after turning off the parallelization for the ISORROPIA driver routine, I can actually spin up the model!

To compile ISORROPIA II for single processor, while compiling everything else for multi-processor, you can add an extra comment character ! to the !$OMP commands in GeosCore/isoropiaII_mod.f, such as:

!!$OMP PARALLEL DO
!!$OMP+DEFAULT( SHARED )
!!$OMP+PRIVATE( I,    J,      L,       N,      WI,   WT,  GAS,  TEMPI )
!!$OMP+PRIVATE( RHI,  VOL,    TSO4,    TNH3,   TNA,  TCL, ANO3, GNO3  )
!!$OMP+PRIVATE( TCA,  TMG,    TK,      CNTRL,  SCASI                  )
!!$OMP+PRIVATE( TNO3, AERLIQ, AERSLD,  OTHER,  TNH4, TNIT             )
!!$OMP+PRIVATE( HPLUSTEMP,    NUM_SAV, DEN_SAV, HNO3_DEN              )
!!$OMP+SCHEDULE( DYNAMIC )

...

!!$OMP END PARALLEL DO

--Bob Y. 10:23, 8 April 2011 (EDT)

Aerosol pH

This issue is somewhat resolved.

Havala O. T. Pye wrote:

At the GEOS-Chem User's meeting (2009) Becky Alexander and her student Eric noted that ISORROPIA II sometimes returns a negative H+ concentration. In my tests, the negative values occurred in N. India and Western Russia area during limited times of the year. A quick fix has been implemented. A condition I tracked down was returning a negative [H+] when solving the HSO4 = H + SO4 equilibrium (CALCHS4 in isorropiaIIcode.f) when the actual answer should have been ~1e-27 mol/m3. (I plugged the same input values into a spreadsheet and it returned the same negative answer so it seems like a numerical precision issue.) For my case study, it came down to subtracting two numbers of very similar value (first 8 digits the same) which resulted in a negative. I put in a fix to reset H+ to 1d-30 in CALCHS4 if it goes negative.

--Bob Y. 10:29, 29 January 2010 (EST)