[Bugfix] Numerical error when running larger number of Jacobian elements

[sabourbaray]

Name and Institution (Required)

Name: Sabour Baray Institution: Environment and Climate Change Canada

Description of your issue or question

Running an experiment with a larger number of state vector elements (~800) produced a checkerboard pattern in the sensitivities beginning at state vector element 373, where some alternating elements had sensitivities of 0 (see figure below).

The issue arises because the Jacobian model runs are not perturbing the requested state vector element correctly, resulting in the same concentration field as the base-case (Jacobian 0000) run. This is due to a numerical error in the comparison in GeosCore/global_ch4_mod.F90:

https://github.com/geoschem/geos-chem/blob/4722f288e90291ba904222f4bbe4fc216d17c34a/GeosCore/global_ch4_mod.F90#L520C1-L535C16

      !------------------------------------------------------------
      ! Perturb emissions for analytical inversion
      !------------------------------------------------------------
      IF ( Input_Opt%AnalyticalInv ) THEN

         ! Only apply emission perturbation to current state vector
         ! element number
         IF ( Input_Opt%StateVectorElement .GT. 0 ) THEN
            IF ( STATE_VECTOR(I,J) == Input_Opt%StateVectorElement ) THEN
               State_Chm%CH4_EMIS(I,J,1) = State_Chm%CH4_EMIS(I,J,1) * Input_Opt%PerturbEmis
               !Print*, 'Analytical Inversion: Scaled state vector element ', &
               !        Input_Opt%StateVectorElement, ' by ', &
               !        Input_Opt%PerturbEmis
            ENDIF
         ENDIF
      ENDIF

After adding diagnostic print statements and running a test case without parallelization in this block of code, we can see there is a small numerical error.

Entered Analytical Inversion block
Current StateVectorElement:         381
No match at indices (          19 ,          22 )
STATE_VECTOR(I,J): 383.000030517578, Input_Opt%StateVectorElement: 381
Entered Analytical Inversion block
Current StateVectorElement:         381
No match at indices (          20 ,          22 )
STATE_VECTOR(I,J): 382.000000000000, Input_Opt%StateVectorElement: 381
Entered Analytical Inversion block
Current StateVectorElement:         381
No match at indices (          21 ,          22 )
STATE_VECTOR(I,J): 380.999969482422, Input_Opt%StateVectorElement: 381
Entered Analytical Inversion block
Current StateVectorElement:         381
No match at indices (          22 ,          22 )
STATE_VECTOR(I,J): 380.000000000000, Input_Opt%StateVectorElement: 381

In this example we can see the state vector element for this run 381 was skipped and not perturbed because the comparison to 380.999969482422 was not true.

A temporary solution to this is to add a tolerance, although this is probably not the preferred approach:

      !------------------------------------------------------------
      ! Perturb emissions for analytical inversion
      !------------------------------------------------------------
      IF ( Input_Opt%AnalyticalInv ) THEN
         ! Only apply emission perturbation to current state vector
         ! element number
         IF ( Input_Opt%StateVectorElement .GT. 0 ) THEN
            IF (ABS(STATE_VECTOR(I,J) - Input_Opt%StateVectorElement) < tolerance ) THEN
               State_Chm%CH4_EMIS(I,J,1) = State_Chm%CH4_EMIS(I,J,1) * Input_Opt%PerturbEmis
               !Print*, 'Analytical Inversion: Scaled state vector element ', &
               !        Input_Opt%StateVectorElement, ' by ', &
               !        Input_Opt%PerturbEmis
            ENDIF
         ENDIF
      ENDIF

I think the original issue in how HCO_GC_EvalFld interacts with the state vector built by the IMI, where there may be conversions between data types that introduces this small numerical error:

! =================================================================
! Get fields for CH4 analytical inversions if needed
! =================================================================
IF ( Input_Opt%AnalyticalInv ) THEN

   ! Evaluate the state vector field from HEMCO
   CALL HCO_GC_EvalFld( Input_Opt, State_Grid, 'CH4_STATE_VECTOR', &
                     STATE_VECTOR, RC)

[msulprizio]

Thanks @sabourbaray! We will try to fix this in GEOS-Chem 14.2.1.

I am tagging @hannahnesser who is also encountering this bug.

[hannahnesser]

Yes, this is very similar to what I find. My problems occur beginning around state vector element ~133--so I wonder a bit if it's limited to large numbers of state vector elements.

@sabourbaray--could you share your state vector file? I wonder if there are similarities in our files that are triggering this bug. It's odd that other inversions haven't encountered it.

[sabourbaray]

@hannahnesser sure, see attached for the state vector: StateVector.zip

[msulprizio]

I have added a fix in this pull request: https://github.com/geoschem/geos-chem/pull/1976.

The issue is that STATE_VECTOR is defined as REAL(FP), while Input_Opt%StateVectorElement is an INTEGER. My solution was to round STATE_VECTOR(I,J) to the nearest integer when comparing the two.

I tested this with @hannahnesser's setup and it now perturbs Jacobian runs that were previously identical to the prior simulation (i.e. where no emissions perturbation was occurring).

[github-actions[bot]]

This issue has been automatically marked as stale because it has not had recent activity. If there are no updates within 7 days it will be closed. You can add the "never stale" tag to prevent the issue from closing this issue.

[github-actions[bot]]

Closing due to inactivity