View source for Aerosol thermodynamical equilibrium

== RPMARES ==

=== Bug in sulfate_mod.f caused by switch to RPMARES ===

'''''Lyatt Jaegle (jaegle@atmos.washington.edu) wrote:'''''

:I've been trying to run [[GEOS-Chem versions under development#v8-01-01|v8-01-01]] in an offline aerosol simulation. It seems to crash in <tt>sulfate_mod.f</tt>, at line 1892

    CALL GET_GNO3( I, J, L, HNO3_kg )

:GET_GNO3 is in <tt>isorropia_mod.f</tt>, which is no longer called (I saw [[#Reverting to RPMARES due to ISORROPIA giving unrealistic values|the discussion on the Wiki page]]), so I think that GAS_HNO3(I,J,L) is not defined.

:What is the best way of fixing this?  Should we simply set HNO3_vv to the monthly mean value of HNO3 read in GLOBAL_HNO3_MOD and bypass isorropia?

'''''Bob Yantosca (yantosca@seas.harvard.edu) wrote:'''''

:The simple thing seems to be something like this:

     USE GLOBAL_HNO3_MOD, ONLY : GET_HNO3_UGM3
     USE DAO_MOD,         ONLY : AIRVOL
 
     ...
 
     IF ( ITS_A_FULLCHEM_SIM() ) THEN
 
        ! Convert HNO3 [v/v] to [equivalents]
        HNO3_vv = STT(I,J,L,IDTHNO3)
        HNO3_eq = HNO3_vv * AD(I,J,L) / ( 28.97d0 / 63d0 ) / 63.d-3
 
     ELSE
 
        !%%% Get gas-phase HNO3 from ISORROPIA code
        !%%%CALL GET_GNO3( I, J, L, HNO3_kg )
               
        ! Get HNO3 in ug/m3, then multiply by volume in m3
        ! and 1e-6 kg/ug to get HNO3 in kg
        HNO3_kg = GET_HNO3_UGM3( I, J, L ) * AIRVOL(I,J,L) * 1e-6
 
        ! Convert HNO3 [kg] first to [v/v]
        HNO3_vv = HNO3_kg * ( 28.97d0 / 63d0 )   / AD(I,J,L)
 
        ! Then convert HNO3 [kg] to [equivalents]
        HNO3_eq = HNO3_kg / 63d-3
 
     ENDIF

:You might want to doublecheck the unit conversion but I think that's right.

--[[User:Bmy|Bob Y.]] 10:38, 26 June 2008 (EDT)

=== Run dies in RPMARES unexpectedly ===

Colette L. Heald <heald@atmos.colostate.edu> reported a run dying silently in RPMARES.

'''''Philippe Le Sager (plesager@seas.harvard.edu) wrote:'''''
:Again aerosol chemistry is the problem. RPMARES algorithm is not handling every possible cases. I found that a negative concentration for H+ can be returned when solving a cubic equation. I made a general fix for it even though it may have been a case of underflow propagation when I caught one instance.

:We previously take care of an overflow case when dealing with ionic activity coefficients, now with your restart file I came across an underflow issue. Instead of too small divisor, I got too large ones basically. The problem can happen either in the divisions or a little bit later when computing the A0 and A1 coefficients, i.e., at lines denoted with a @ here:

    DO NNN = 1, 50 ! loop count reduced 04/09/2001 by FSB
    NITR = NNN
 
    ! set up equilibrium constants including activities
    ! solve the system for hplus first then sulfate & nitrate
  @ RK2SA = K2SA * GAMAS2 * GAMAS2 / (GAMAS1 * GAMAS1 * GAMAS1)
  @ RKNA = KNA / ( GAMANA * GAMANA )
    RKNWET = RKNA * WH2O
    T21 = ZSO4 - MNH4
    T221 = ZSO4 + T21
 
    ! set up coefficients for cubic
    A2 = RK2SA + RKNWET - T21
  @ A1 = RK2SA * RKNWET - T21 * ( RK2SA + RKNWET )
   & - RK2SA * ZSO4 - RKNA * TNO3
  @ A0 = - (T21 * RK2SA * RKNWET
    & + RK2SA * RKNWET * ZSO4 + RK2SA * RKNA * TNO3 )
 
    CALL CUBIC ( A2, A1, A0, NR, CRUTES )

: So, I tried to make a general fix that handle both cases of small and large GAMAxx and should be compiler/machine independent, with a new subroutine in <tt>error_mod.f</tt> and a new <tt>rpmares_mod.f</tt>. You can grab them at: ftp://ftp.as.harvard.edu/pub/exchange/phs/GEOS-Chem-UPDATE/

--[[User:Plesager|phs]] 16:57, 6 June 2008 (EDT)

=== Fix for floating point error in RPMARES ===

In [[GEOS-Chem_versions_under_development#v8-01-01|GEOS-Chem v8-01-01]], we have added a fix into routine RPMARES (in rpmares_mod.f) to avoid a floating point error that can occur under some conditions.  

'''''Philippe Le Sager (plesager@seas.harvard.edu) wrote:'''''
:In RPMARES, the hydrogen ion (HPLUS) obtained as the root of a cubic equation is so large (286 moles/kg) that it triggers a large binary activity coefficient (LGAMA0) and then a very small mean molal ionic activity (GAMA=10**TRM=2.44e-159) when input into ACTCOF. That all seems correct, but this is so close to zero (~2e-159) that an error occurs when re-injected into RPMARES to compute the coefficients of the cubic equation for H+ at the next step.

:We basically need to shunt the cubic equation for zero ionic activity coefficients to solve for H+, and then sulfate and nitrate.

'''''Rokjin Park (rjpark@snu.ac.kr) replied:'''''
:I agree with Philippe that the case with extremely low NH3 cases we don’t need to solve the cubic equation.
:RPMARES should simply return the original values to avoid the NaN error.

'''''Bob Yantosca (yantosca@seas.harvard.edu) replied:'''''
:FYI, Philippe and I have put in a fix in routine RPMARES:

       CALL ACTCOF ( CAT, AN, GAMS, MOLNU, PHIBAR )
 
       GAMANA = GAMS( 1, 2 )
       GAMAS1 = GAMS( 1, 1 )
       GAMAS2 = GAMS( 1, 3 )
       GAMAAN = GAMS( 2, 2 )
 
       !------------------------------------------------------------
       ! Add robustness: now check if GAMANA or GAMAS1 is too small
       ! for the division in RKNA and RK2SA. If they are, return w/
       ! original values: basically replicate the procedure used
       ! after the current DO-loop in case of no-convergence
       ! (phs, bmy, rjp, 4/10/08)
       !--------------------------------------------------------------
       IF ( ( ABS( GAMANA ) < EPS ) .OR.
    &       ( ABS( GAMAS1 ) < EPS ) ) THEN
 
          ! Reset to original values
          ANH4  = TNH4 * MWNH4
          GNH3  = FLOOR
          GNO3  = GNO3_IN
          ANO3  = ANO3_IN
          ASO4  = TSO4 * MWSO4
          AHSO4 = FLOOR
 
          ! Update water
          CALL AWATER ( IRH, TSO4, TNH4, TNO3, AH2O )
 
          ! Exit this subroutine
          RETURN
       ENDIF

:so then if either of the GAMANA or GAMAS1 variables that get returned from ACTCOF is zero (or smaller than a minimum value) then we don't continue the iteration but then reset to original values and exit the subroutine.  This definitely fixed the problem.

--[[User:Bmy|Bmy]] 16:21, 11 April 2008 (EDT)

=== Reverting to RPMARES due to ISORROPIA giving unrealistic values ===

Due to the [[#01-Apr-2008|ISORROPIA bug at low RH]] discussed below, we will revert to the RPMARES aerosol thermodynamic equilibrium code in [[GEOS-Chem_versions_under_development#v8-01-01|GEOS-Chem v8-01-01]].

'''''Daniel Jacob (djacob@fas.harvard.edu) wrote:'''''
:Thanks! The HNO3 simulation is doing much better but I still wonder if the high values in the free troposphere (where the RH is low) are caused by the ISORROPIA bug. We do need to replace ISORROPIA with RPMARES in the standard model until Havala gives us the corrected ISORROPIA version.  Daniel 

'''''Bob Yantosca (yantosca@seas.harvard.edu) replied:'''''
:Hi all,

:I will put RPMARES into the std code.  I will do 2 1-yr benchmark simulations with the GEOS-5 v8-01-01 std code:

:* v8-01-01-geos5 w/ ISORROPIA (for clean comparison between v8-01-01 and previous models)
:* v8-01-01-geos5 w/ RPMARES (for clean comparison between RPMARES & ISORROPIA w/ the same version)
 
:RPMARES will be used in all NRT-ARCTAS runs for model dates 4/1/08 and beyond.

:Thanks and stay tuned,
:Bob Y.

--[[User:Bmy|Bmy]] 13:16, 4 April 2008 (EDT)

== ISORROPIA ==

=== 01-Apr-2008 ===

The GEOS-5 met fields [[GEOS-5 issues|incorrectly reported the relative humidity (RH) field as fraction instead of percent.]]  Therefore...

'''''Bob Yantosca (yantosca@seas.harvard.edu) wrote:'''''

:... GEOS-Chem was expecting to read in RH in percent (0..100) but was actually reading in fraction (0..1).  Then when we got to the places where RH was used in the code (e.g. ISORROPIA, drydep_mod.f, sulfate_mod.f), we were dividing that by 100 again.  So then this would further reduce the maximum RH from order 1 to order 0.01.  

:One of the symptoms that we saw was that HNO3 (in many locations, though not all) was completely wiped out in the PBL.  

:Would you expect a reduction in HNO3 by ISORROPIA if it were given 100x too small RH?  Do you know if ISORROPIA would give junk results in this case?

'''''Havala Olson Taylor Pye (havala@caltech.edu) replied:'''''

:The current version of ISORROPIA in the code definitely has problems at low RH and this is one of my primary motivations for updating it. I found that under certain conditions, at low RH, ISORROPIA would convert all nitrate (HNO3 + NIT) to the aerosol phase which leads to significant overpredictions in NIT. In the attached plot, below an RH of about 0.30 (30%) all nitrate was converted to the aerosol phase (the total amounts of all the species and T are the same for each data point). So by multiplying the RH by 100 to correct it, you probably moved from a region where ISORROPIA makes no sense to a region where is actually gives good predictions. Hope this clears things up. I've been working on ISORROPIAII and it gives predictions very similar to RPMARES (as far as where nitrate concentrations are high and low) and very different from ISORROPIA currently in GEOS-Chem. Let me know if you want any more information.

:Havala

--[[User:Bmy|Bmy]] 13:15, 4 April 2008 (EDT)