[BUG/ISSUE] Potential error in gckpp_HetRates.F90: 1st order uptake of VOCs

I am not sure if this is a bug or a "feature", but I noticed what may be an error where we compute 1st order uptake of VOC species. If you look in routines HetLVOC (and in similar routines), there code such as:

!EOC
!------------------------------------------------------------------------------
!                  GEOS-Chem Global Chemical Transport Model                  !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: HetLVOC
!
! !DESCRIPTION: Sets the heterogenous chemistry rate for LVOC: condensation of
! low-volatility ISOPOOH oxidation products.
!\\
!\\
! !INTERFACE:
!
    FUNCTION HETLVOC( A, B ) RESULT( HET_LVOC )
!
! !INPUT PARAMETERS:
!
      ! Rate coefficients
      REAL(fp), INTENT(IN) :: A, B
!
! !RETURN VALUE:
!
      REAL(fp)             :: HET_LVOC
!
! !REMARKS:
!
! !REVISION HISTORY:
!  15 Jun 2017 - M. Sulprizio- Initial version based on calcrate.F from E.Marais
!  See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
      LOGICAL  :: DO_EDUCT
      INTEGER  :: N
      REAL(fp) :: XSTKCF, ADJUSTEDRATE

      ! Initialize
      HET_LVOC     = 0.0_fp
      ADJUSTEDRATE = 0.0_fp
      XSTKCF       = 0.0_fp

      ! Don't do PSC rate adjustment
      DO_EDUCT     = .FALSE.

      ! Loop over aerosol types
      DO N = 1, NAEROTYPE    ! NAEROTYPE=14

         ! Define gamma

         IF (N.eq.8 .and. RELHUM >= CRITRH) THEN    ! Note, CRITRH = 35%
            XSTKCF = B
         ENDIF

         IF (N.eq.13) THEN
            ! Calculate for stratospheric liquid aerosol
            ! Note that XSTKCF is actually a premultiplying
            ! factor in this case, including c-bar
            ADJUSTEDRATE = XAREA(N) * XSTKCF
         ELSE
            ! Reaction rate for surface of aerosol
            ADJUSTEDRATE=ARSL1K(XAREA(N),XRADI(N),XDENA,XSTKCF,XTEMP, &
                               (A**0.5_FP))
         ENDIF

         IF ( DO_EDUCT .and. N > 12 ) THEN
            ! PSC reaction - prevent excessive reaction rate
            IF (ADJUSTEDRATE.gt.(1.e+0_fp/HetMinLife)) THEN
               ADJUSTEDRATE = 1.e+0_fp/HetMinLife
            ENDIF
         ENDIF

         ! Add to overall reaction rate
         HET_LVOC = HET_LVOC + ADJUSTEDRATE
      END DO

    END FUNCTIOn HETLVOC

Following the execution of the loop:

Result HET_LVOC = 0 for aerosol types N=1 to N=7
Reaction probability xstkcf is set to the 2nd argument B only if N=8 (tropospheric sulfate) and if RH > 35%,
Uptake of VOC on aerosol type N=8 (trop sulfate) is computed and stored in HET_LVOC.
For aerosol types N > 8, xstckf is not reset to zero (but I think it should be).
a. For grid boxes where RH > 35%, The uptake of VOC by aerosol types N=9 thru N=14 is done.
b. But for grid boxes where RH <=35%, uptake of VOC is not computed.

Long story short, I think there is a missing statement to zero out xstkcf at the beginning of each N iteration:

      ! Loop over aerosol types
      DO N = 1, NAEROTYPE

         !##### NEED TO ADD THIS HERE! ######
         XSTCKF = 0.0_dp
         !################################

         ! Define gamma

         IF (N.eq.8 .and. RELHUM >= CRITRH) THEN    ! Note, CRITRH = 35%
            XSTKCF = B
          ENDIF

This will then skip computing the uptake of VOC on other aerosol types, regardless of the value of RH in the grid box.

@sdeastham, @msulprizio, @msl3v, @lizziel: any thoughts?

I am also attempting to replace the code above by this cleaner formulation. To match the behavior of the existing code, I have to do:

    !========================================================================
    ! Het1stOrderUptakeVOC begins here!
    !========================================================================

    ! Initialize
    rate  = 0.0_dp

    ! Only consider inorganic aqueous aerosols with RH > 35%.
    IF ( RELHUM >= CRITRH ) THEN

       ! Uptake by tropospheric sulfate (aerosol type 8)
       rate = rate + ARSL1K( XAREA(8 ), XRADI(8 ), XDENA, gamma, XTEMP, SrMw )

       ! Uptake by black carbon and organic carbon (aerosol types 9-10)
       rate = rate + ARSL1K( XAREA(9 ), XRADI(9 ), XDENA, gamma, XTEMP, SrMw )
       rate = rate + ARSL1K( XAREA(10), XRADI(10), XDENA, gamma, XTEMP, SrMw )

       ! Uptake by fine & coarse sea salt (aerosol types 11-12)
       rate = rate + ARSL1K( XAREA(11), XRADI(11), XDENA, gamma, XTEMP, SrMw )
       rate = rate + ARSL1K( XAREA(12), XRADI(12), XDENA, gamma, XTEMP, SrMw )

       ! Uptake by stratospheric sulfate (aerosol type 13)
       ! and by irregular ice cloud (aerosol type 14)
       rate = rate + XAREA(13) * gamma
       rate = rate + ARSL1K( XAREA(14), XRADI(14), XDENA, gamma, XTEMP, SrMw )

    ENDIF

but if there is a bug in the code from the prior comment and we should be setting the reaction proabability to zero after each iteration, then this all reduces to:

    IF ( RELHUM >= CRITRH ) THEN

       ! Uptake by tropospheric sulfate (aerosol type 8)
       rate = rate + ARSL1K( XAREA(8 ), XRADI(8 ), XDENA, gamma, XTEMP, SrMw )

    ENDIF

geoschem / geos-chem

[BUG/ISSUE] Potential error in gckpp_HetRates.F90: 1st order uptake of VOCs #668