[BUG/ISSUE] Potential div-by-zero error in emissions_mod (routine MMR_Compute_Flux)

[yantosca]

What institution are you from?

GCST

Overview

While working on GCHP, I discovered a potential issue in the routine MMR_Compute_Flux, which is used for certain passive species in the TransportTracers simulation. The code could potentially lead to a div-by-zero error if the Total_Area variable is zero in this code block:

geos-chem/GeosCore/emissions_mod.F90

Lines 488 to 550 in 9ab6931

	!=======================================================================
	! Compute the surface flux needed to restore the total burden
	!=======================================================================
	! Loop over species
	DO N = 1, nAdvect
	! Point to the Species Database entry for species N
	SpcInfo => State_Chm%SpcData(N)%Info
	! Only do calculation if this is an MMR tracer
	IF ( TRIM(SpcInfo%Name) == 'GlobEmis90dayTracer' .OR. &
	TRIM(SpcInfo%Name) == 'NHEmis90dayTracer' .OR. &
	TRIM(SpcInfo%Name) == 'SHEmis90dayTracer' ) THEN
	! Initialize
	Total_Spc = 0.0e+0_fp
	Total_Area = 0.0e+0_fp
	Mask = 1.0e+0_fp
	! Loop over grid boxes
	DO L = 1, State_Grid%NZ
	DO J = 1, State_Grid%NY
	DO I = 1, State_Grid%NX
	! Compute mol of Tracer needed to achieve the desired value
	Total_Spc = Total_Spc + &
	( GlobalBurden - State_Chm%Species(I,J,L,N)) * &
	(State_Met%AIRNUMDEN(I,J,L)/ AVO) * State_Met%AIRVOL(I,J,1)
	! To distribute it uniformly on the surface, compute the total
	! area [m2]
	IF ( L == 1 ) THEN
	! Latitude of grid box
	YMID = State_Grid%YMid(I,J)
	! Define mask if needed
	IF (TRIM(SpcInfo%Name)=='NHEmis90dayTracer') THEN
	IF ( YMID < 0.0 ) MASK(I,J) = 0.0e+0_fp
	ELSE IF (TRIM(SpcInfo%Name)=='SHEmis90dayTracer') THEN
	IF ( YMID >= 0.0 ) MASK(I,J) = 0.0e+0_fp
	ENDIF
	Total_Area = Total_Area + State_Grid%Area_M2(I,J) * MASK(I,J)
	ENDIF
	ENDDO
	ENDDO
	ENDDO
	! Compute flux [mol/m2]
	Flux(:,:) = ( Total_Spc / Total_Area ) * MASK(:,:)
	! Update species concentrations [mol/mol]
	Spc(:,:,1,N) = Spc(:,:,1,N) + Flux(:,:) * &
	AVO / ( State_Met%BXHEIGHT(:,:,1) * State_Met%AIRNUMDEN(:,:,1) )
	ENDIF ! MMR tracer
	! Free pointers
	SpcInfo => NULL()
	ENDDO

Total_Area is multiplied by MASK here, which can be zero for certain passive species depending on latitude.

Solution

We can rearrange the code so that we avoid the div-by-zero condition as such:

    ! Compute the surface flux needed to restore the total burden
    !========================================================================

    ! Loop over species
    DO N = 1, nAdvect

       ! Point to the Species Database entry for species N
       SpcInfo => State_Chm%SpcData(N)%Info

       ! Only do calculation if this is an MMR tracer
       IF ( TRIM(SpcInfo%Name) == 'GlobEmis90dayTracer'   .OR.               &
            TRIM(SpcInfo%Name) == 'NHEmis90dayTracer'     .OR.               &
            TRIM(SpcInfo%Name) == 'SHEmis90dayTracer'   ) THEN

          ! Initialize
          Total_Spc   = 0.0_fp
          Total_Area  = 0.0_fp
          Mask        = 1.0_fp
          Flux        = 0.0_fp

          !------------------------------------------------------------------
          ! To distribute it uniformly on the surface, 
          ! compute the total area [m2]
          !------------------------------------------------------------------
          DO J = 1, State_Grid%NY
          DO I = 1, State_Grid%NX

             ! Latitude of grid box
             YMID = State_Grid%YMid(I,J)

             ! Define mask if needed
             IF (TRIM(SpcInfo%Name) == 'NHEmis90dayTracer') THEN
                IF ( YMID <  0.0_fp ) MASK(I,J) = 0.0_fp
             ELSE IF (TRIM(SpcInfo%Name) == 'SHEmis90dayTracer') THEN
                IF ( YMID >= 0.0_fp ) MASK(I,J) = 0.0_fp
             ENDIF

             ! Accumulate the total area
             Total_Area = Total_Area + ( State_Grid%Area_M2(I,J) * MASK(I,J) )
          ENDDO
          ENDDO
 
          !------------------------------------------------------------------
          ! We only need to update concentrations when Total_Area > 0.
          ! If Total_Area == 0, then we need to skip the computation of
          ! Flux, or else we will incur a div-by-zero condition.
          !
          !   -- Bob Yantosca (11 Apr 2022)
          !------------------------------------------------------------------
          IF ( Total_Area > 0.0_fp ) THEN 

             ! Compute mol of Tracer needed to achieve the desired value
             DO L = 1, State_Grid%NZ
             DO J = 1, State_Grid%NY
             DO I = 1, State_Grid%NX
                Total_Spc = Total_Spc                                        &
                          + ( GlobalBurden - State_Chm%Species(I,J,L,N)  )   &
                          * ( State_Met%AIRNUMDEN(I,J,L) / AVO           )   &
                          * State_Met%AIRVOL(I,J,1)
             ENDDO
             ENDDO
             ENDDO 

             ! Compute flux [mol/m2]
             Flux(:,:)    = ( Total_Spc / Total_Area )

             ! Update species concentrations [mol/mol]
             Spc(:,:,1,N) = Spc(:,:,1,N)                                     &
                          + ( Flux(:,:) * AVO            )                   &
                          / ( State_Met%BXHEIGHT(:,:,1)                      &
                          *   State_Met%AIRNUMDEN(:,:,1) )
          ENDIF

       ENDIF ! MMR tracer

       ! Free pointers
       SpcInfo => NULL()

    ENDDO

Also note, I believe the multiplication of Flux by MASK(:,:) is now redundant, as the mask has already been multiplied by total area. So when MASK=0, Total_Area=0 as well.

[yantosca]

NOTE: I think this issue does not affect GEOS-Chem Classic, because the domain is global. But it is possible that one core might have all grid boxes where the mask is zero, and thus, a division-by-zero error will happen on that core. Recall that in the standard cubed-sphere grid (without stretching or rotation), one of the faces is completely in the southern hemisphere, and another face is completely in the northern hemisphere.

[yantosca]

Actually I think we should leave this line as-is:

             ! Compute flux [mol/m2]
             Flux(:,:)    = ( Total_Spc / Total_Area ) * MASK(:,:)

[yantosca]

Just realized that for GCHP, we do not yet use the 90 day passive species (GlobEmis90DayTracer, NHEmis90DayTracer, SHEmis90DayTracer). So this is not as critical of a fix, but it should go into the code at some point (probably for 14.0.0).

[yantosca]

Further discussion will be in PR #1216

[stale]

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 Stale bot from closing this issue.

[msulprizio]

This fix has now been pushed to dev and will be included in 14.0.0.