Create variables to track non-passive species, so that they can be excluded from the Rosenbrock error norm (in int/rosenbrock*.f90 integrators only)

[yantosca]

Overview

This is the companion PR to #66 and #144. In this PR we create variables that can be used to exclude dummy species (i.e. those added to chemical reactions for the purpose of obtaining diagnostic information) from the Rosenbrock error norm computation. This should prevent the issue where integration results are changed depending on the number of dummy species in the mechanism (cf geoschem/geos-chem#3175).

Technical description

Two new module-level variables have been added to the ROOT_Global module (for F90 only). These are:

1. NonPassiveSpc_Count: The number of non-passive species
2. NonPassiveSpc_Indices: The KPP species indices for non-passive species

These optional variables (which are stored in ROOT_Global) are constructed when the Initialize routine is called.

• If Initialize is called with the PassiveSpc_ATOL_Threshold optional argument, then all species that have an absolute tolerance (ATOL) greater than or equal to PassiveSpc_ATOL_Threshold will be denoted as passive species. In this case, NonPassiveSpc_Count will be equal to NVAR - the number of passive species and NonPassiveSpc_Indices will contain the index list of only the non-passive species.

• If Initialize is called without any optional arguments, then NonPassiveSpc_Count will be equal to NVAR and NonPassiveSpc_Indices will contain the index list of all variable species.

We have also updated the algorithm in the Ros_ErrorNorm function (in each of the int/rosenbrock*.f90 files) to compute the error norm on the basis of non-passive species. We have also optimized the function to facilitate vectorization and to remove unnecessary ELSE blocks, which can be a bottleneck:

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
                               AbsTol, RelTol, VectorTol )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~> Computes the "scaled norm" of the error vector Yerr
!~~~>
!~~~> Now uses separate loops for scalar and vector tolerances,
!~~~> as this facilitates loop vectorization for better performance.
!~~~>
!~~~> Also uses NonPassiveSpc_Count and NonPassiveSpc_Indices (constructed
!~~~> in Initialize), so that we can exclude "passive" (e.g. prod/loss)
!~~~> species from the error norm computation.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IMPLICIT NONE

! Input arguments
   KPP_REAL, INTENT(IN) :: Y(N), Ynew(N), Yerr(N), AbsTol(N), RelTol(N)
   LOGICAL, INTENT(IN) ::  VectorTol
! Local variables
   KPP_REAL :: Err, Scale, Ymax
   INTEGER  :: I, IDX

   Err = ZERO

   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   !~~~> Vector Tolerances (per-species AbsTol & RelTol)
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IF ( VectorTol ) THEN

      DO IDX = 1, NonPassiveSpc_Count
         I     = NonPassiveSpc_Indices(IDX)
         Ymax  = MAX( ABS( Y(I) ), ABS( Ynew(I) ) )
         Scale = AbsTol(I) + RelTol(I) * Ymax
         Err   = Err + ( Yerr(I) / Scale )**2
      ENDDO

      ! Normalize the error norm by the number of non-dummy species
      ! and prevent it from getting smaller than 1e-10
      Err           = SQRT( Err / NonPassiveSpc_Count )
      ros_ErrorNorm = MAX( Err, 1.0d-10 )
      RETURN

   ENDIF

   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   !~~~> Scalar Tolerance (same AbsTol & RelTol for all species)
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO IDX = 1, NonPassiveSpc_Count
      I     = NonPassiveSpc_Indices(IDX)
      Ymax  = MAX( ABS( Y(I) ), ABS( Ynew(I) ) )
      Scale = AbsTol(1) + RelTol(1) * Ymax
      Err   = Err + ( Yerr(I) / Scale )**2
   ENDDO
   
   ! Normalize the error norm by the number of non-dummy species
   ! and prevent it from getting smaller than 1e-10
   Err           = SQRT( Err / NonPassiveSpc_Count )
   ros_ErrorNorm = MAX( Err, 1.0d-10 )

  END FUNCTION ros_ErrorNorm

Other modifications bundled into this PR:

• Added C-I test F90_ros_passivespc, which adds 3 passive species to the small_strato mechanism.
• Added new driver program drv/general_passivespc.f90, which executes CALL INITIALIZE( PassiveSpc_ATOL_Threshold = 1.0d25 ). This is used in the F90_ros_passivespc C-I test.
• Added KPP internal function F90_FunctionBeginNoArgsDecl (in src/code_f90.c). This is similar to F90_FunctionBegin, except that it does not write any subroutine arguments. This allows us to manually declare subroutine below any F90 USE statements.
• Added macro DefElmO in src/code.h, which defines a scalar F90 variable with the OPTIONAL attribute.
• Added GNU 14, 15, 16 compilers to the "Run C-I tests" GitHub action.
• Set Hacc = ZERO; and ErrOld = ZERO; in int/runge_kutta.c to avoid generating compiler warnings during C-I tests
• Updated ICNTRL and RCNTRL comment headers in the int/rosenbrock*.f90 integrator files.
• Updated ReadtheDocs documentation accordingly.

Related GitHub issues

• Closes [BUG/ISSUE] Prod/loss species are contributing to the Error Norm when they should be ignored #66
• Closes Adding dummy species into reactions (to track P/L rates) changes results (in KPP 3.3.0) #144
• See Adding inert reaction tags causes changes to chemistry in 14.6.3 geoschem/geos-chem#3175

Tagging @RolfSander @msl3v @jimmielin @obin1 @fsolmon @ujongebloed @mje1398

[RolfSander]

It's good to see that there is now a possibility to exclude dummy
species from the Rosenbrock error calculation! I do have a suggestion
though:

With the current implementation, you can use the optional argument only
for the four selected integrators. When you switch to another integrator
and forget to remove the optional argument again, the compiler will
produce an error. Of course, we could add the new optional arguments to
all integrators but that's a bit tedious.

My suggestion is that we put NonDummySpc_Count into ICNTRL(19). The
calculation of NonDummySpc_Indices could be done by KPP automatically at
the end of SUBROUTINE Initialize (which is generated by
GenerateInitialize() in gen.c).

Would that work?

[yantosca]

Thanks @RolfSander, I think that's a good idea. I'll try it.

[yantosca]

@RolfSander: I just noticed that there is the #DUMMYINDEX option that lets you set the indices of species that aren't in the mechanism to zero. I was thinking I should rename e.g. NonDumSpc_Count etc. so as not to confuse our nomenclature. Maybe I should refer to these as "diagnostic species" (e.g. NonDiagSpc_Count, etc)? Or maybe "actual" species?

[RolfSander]

NonDiagSpc_Count

Yes, we should avoid confusion. Unfortunately, diag occurs quite often
in the KPP code, meaning diagonal, not diagnostic.

ActualSpec_Count would be ok.

Here is yet another idea: During the discussion, @obin1 used the name
"P/L species". So maybe ProdLossSpc_Count?

[obin1]

I like ProdLossSpc_Count! This is a good name if we only use them for the production/loss diagnostics.

If we are looking for something more general, would PassiveSpc_Count work? I would imagine we want all passive species (not necessarily inert, but passive diagnostic species) to be excluded from the adaptive timestepping strategy.

[yantosca]

Thanks @RolfSander. I think NonPassive_SpcCount etc is more general. I'll try that.

[RolfSander]

Thanks, @yantosca! The name NonPassiveSpc_Count is fine. However, I
still see one problem: The new code is not at the end of SUBROUTINE Initialize but before the inlined code from F90_INIT.

As I am using F90_INIT to define my tolerances, this means that
NonPassiveSpc_Count will be calculated using potentially obsolete
ATOL values.

Would it be okay to move the new code after the inlined F90_INIT
block?

[yantosca]

Thanks @RolfSander. Yes, I can make that change. I'll also update the documentation as well. Stay tuned.

[RolfSander]

Thanks, @yantosca!, it looks good now!