[FEATURE REQUEST] Add toggle to skip using EQUIVALENCE statements

[yantosca]

In gckpp_Global.F90, we have these statements

! VAR, FIX are chunks of array C
      EQUIVALENCE( C(1),VAR(1) )
      EQUIVALENCE( C(75),FIX(1) )

However, the Fortran EQUIVALENCE statement is not thread-safe, and thus cannot be used for multi-threaded processing (i.e. with OpenMP parallel DO loops).

For GEOS-Chem we have deleted these statements and just use COPY statements before the solver

       !=====================================================================
       ! Initialize the KPP "C" vector of species concentrations [molec/cm3]
       !=====================================================================
       DO N = 1, NSPEC
          SpcID = State_Chm%Map_KppSpc(N)
          C(N)  = 0.0_dp
          IF ( SpcId > 0 ) C(N) = State_Chm%Species(I,J,L,SpcID)
       ENDDO
       ...
       ! Set VAR and FIX arrays
       VAR(1:NVAR) = C(1:NVAR)
       FIX         = C(NVAR+1:NSPEC)

and after the solver

       ! Copy VAR and FIX back into C (mps, 2/24/16)
       C(1:NVAR)       = VAR(:)
       C(NVAR+1:NSPEC) = FIX(:)

I propose that we add a toggle that would allow us to skip adding the EQUIVALENCE statements for use-cases where KPP is done within a parallel loop.

[RolfSander]

Yes, let's get rid of the EQUIVALENCE statements. With the copy statements, however, twice as much memory is used. Also, we would create inconsistencies if the C array is used somewhere inside a solver instead of VAR or FIX. Maybe we can use POINTERs to replace EQUIVALENCE. I haven't tested it but I think it could work like this:

REAL, TARGET,  DIMENSION(NSPEC) :: C
REAL, POINTER, DIMENSION(:)     :: VAR, FIX

VAR => C(1:NVAR)
FIX => C(NVAR+1:NSPEC)

[RolfSander]

Here is a small test program that works fine for me:

PROGRAM pointer_target

  INTEGER, PARAMETER :: NSPEC = 10
  INTEGER, PARAMETER :: NVAR = 7

  REAL, TARGET,  DIMENSION(NSPEC) :: C
  REAL, POINTER, DIMENSION(:)     :: VAR, FIX

  VAR => C(1:NVAR)
  FIX => C(NVAR+1:NSPEC)

  ! define C, then check VAR and FIX:
  C = (/ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10. /)
  PRINT *, C
  PRINT *, VAR
  PRINT *, FIX

  ! edit the values:
  C(3)   = 33.
  VAR(2) = 22.
  FIX(2) = 99.
  PRINT *, C
  PRINT *, VAR
  PRINT *, FIX

END PROGRAM pointer_target

[yantosca]

Thanks Rolf. I will work on this as a separate PR. Stay tuned.

[yantosca]

Am wondering if we should put the

VAR=> C(1:NSPEC)
FIX => C(NVAR+1:NSPEC)

statements in the _Initialize.F90 file. We might also want to create a _Finalize.F90 file with a Finalize() routine that nullifies the pointers. For OpenMP parallelization, the pointers need to be assigned and nullified on the same loop iteration.

[RolfSander]

I agree that _Initialize.f90 is probably the right place. More precisely, the SUBROUTINE Initialize() inside this file. However, are we sure that this subroutine is always called first? Is it really mandatory? I cannot imagine any application that would skip it but I just want to make sure...

I also agree that the pointers should be nullified eventually. However, what do you mean by "same loop iteration"? As far as I can see, the pointers should be nullified at the end of the model run, not after each KPP call. But maybe I misunderstood you...

[yantosca]

So what I meant is that if you have an OpenMP parallel loop such as:

REAL(dp), TARGET  :: my_variable
REAL(dp), POINTER :: myPointer
!$OMP THREADPRIVATE ( myPointer )

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L ) 
DO L = 1, NZ
DO J = 1, NY
DO I = 1, NX
   ...
   myPointer => myVariable
   ...
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
myPointer => NULL()

Then this would lead to a segmentation fault because the pointer hasn't been nullified before the start of the next iteration. So what we'd need to do is this:

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L ) 
DO L = 1, NZ
DO J = 1, NY
DO I = 1, NX
   ...
   myPointer => myVariable
   ...

   myPointer => NULL()
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

So in other words, variables that are private to the loop don't exist outside of the loop, because the compiler will create as many copies of them on the stack as there are threads. So if you try to access the loop variable after the loop ends, the memory isn't there anymore.

[yantosca]

So if we have a place where we assign pointers VAR & FIX (say it's in routine Initialize) then we'd also need a way of nullifying the pointers. We could do it manually in the KPP driver program but then the user might not be aware that this needs to be done. Maybe that's the best approach.

[RolfSander]

OK, I think I understand. You want to create and nullify the pointer at every KPP call. In that case, however, _Initialize.f90 is not the right place. SUBROUTINE Initialize is only called once at the beginning of the model run, not at every KPP call. At least this is how we do it in MECCA. Is GEOS-Chem different?

To create and nullify the pointer at every KPP call, I think we would have to add the corresponding code to SUBROUTINE integrate in all integrators.

[yantosca]

That sounds good, to do it in the integrators.

Our Initialize routine for GEOS-Chem looks like this. Not sure if this is a modification done by @msl3v:

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Initialize - function to initialize concentrations
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Initialize ( )

  USE gckpp_Global

  INTEGER :: i
  REAL(kind=dp) :: x

  CFACTOR = 1.000000e+00_dp

  x = (0.)*CFACTOR
  DO i = 1, NVAR
    VAR(i) = x
  END DO

  x = (0.)*CFACTOR
  DO i = 1, NFIX
    FIX(i) = x
  END DO

! constant rate coefficients
! END constant rate coefficients

! INLINED initializations

! End INLINED initializations

END SUBROUTINE Initialize

[yantosca]

so we do call that on each (I,J,L) grid box since the vectors are private to the loop.

[RolfSander]

Your SUBROUTINE Initialize looks quite familiar, I don't see any modifications from @msl3v. Also, it makes sense to me that you call it for each grid box. However, do you really call it at every time step? I don't think this is the original intention for this subroutine.

Anyway, yes, let's put the pointer assignments and nullifications at the start and end, respectively, of SUBROUTINE integrate in the integrator files.

[msl3v]

I don't think I ever touched this.