Scaling issue in HEMCO for extension emissions like fire and soil NOx

[Ruijun-Dang]

I attempted to scale the GFED NH3 emissions within a specified region following the HEMCO guideline. It was successful only when using intel19 compiler, but it failed if switching to gfortran (commonly used) or intel23. This issue arose because our Cannon no longer supports intel19 after this recent maintenance.

I am using HEMCO v3.4.0, GC v13.4.1, and gfortran v10.2.0 and tried to locate the problem. I found that HEMCO didn't read the ScaleField I added on Line 213 of HEMCO_Config.rc (attached). This happened because HEMCO didn't execute Line 872, which is responsible for obtaining the ScaleField. Strangely, HEMCO incorrectly determined that PRESENT(SpcScal_char) is FALSE, despite the existence of parameter SpcScal_char. I am unsure about how is this issue related to the different compliers and how to fix this especially for gfortran.

HEMCO_Config.txt

[yantosca]

Thanks @Ruijun-Dang. I might try to build a test code to test only the subroutine if the various arguments are passed. Stay tuned.

[Ruijun-Dang]

Hi Bob, Thanks for checking! TC1 is a variable that I added earlier for the new soil NOx scheme, and in my current code, TC1 is already commented out. It’s strange that the code you copied differs from mine.

[yantosca]

Hi @Ruijun-Dang. I've put in some debug printout in HEMCO. It does seem to be finding the name of the scale field to apply:

----- 111
   2 Scaling_CO      1.050000
  25 Scaling_NAP      0.000275
----- 111
  14 ScaleField_NH3 GFED_NH3_SCALE

I also put in some debug output in this routine that gets called from hcox_gfed_mod.F90

      ! Check for masking
       CALL HCOX_SCALE( HcoState, SpcArr, TRIM(Inst%SpcScalFldNme(N)), RC )
       IF ( RC /= HCO_SUCCESS ) THEN
           CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
           RETURN
       ENDIF
       stop

this calls the routine HCOX_SCALE_dp2d in Extensions/hcox_tools_mod.F90:

    IF ( TRIM(SCALENAME) /= TRIM(HCOX_NOSCALE) ) THEN

       ! Get mask field
       CALL HCO_EvalFld ( HcoState, TRIM(SCALENAME), SCAL, RC )
       IF ( RC /= HCO_SUCCESS ) THEN
           CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
           RETURN
       ENDIF

       print*, 'sum arr  before : ', sum(arr)
       print*, 'sum scal before : ', sum(scal)

       DO J = 1, HcoState%NY
       do I = 1, HcoState%NX
          if ( arr(i,j) > 0.0 ) then
             print*, '==> ', I, J, arr(i,j), scal(i,J), arr(i,j)*scal(i,j)
          endif
       enddo
       enddo

       ! Set array to zero outside of mask region
       Arr = Arr * SCAL

       print*, 'sum arr after   : ', sum(arr)

    ENDIF

    ! Return w/ success
    RC = HCO_SUCCESS

when I do this, I get output such as:

 ==>            1          86   3.2377348827640749E-013   1.0000000000000000        3.2377348827640749E-013
 ==>            2          86   5.1102400276166482E-013   1.0000000000000000        5.1102400276166482E-013
 ==>            3          86   1.7700078009439567E-013   1.0000000000000000        1.7700078009439567E-013
 ==>            4          86   3.6869739146839539E-013   1.0000000000000000        3.6869739146839539E-013
 ==>            5          86   2.2345554195592751E-015   1.0000000000000000        2.2345554195592751E-015
etc. everywhere

which makes sense given the definition of GFED_NH3_SCALE:

111 GFED_NH3_SCALE  1                                                         -             2000/1/1/0             C  xy unitless  * 666 1 1 

If I change the GFED_NH3_SCALE to 2, I get the same output. Looking into this now.

[yantosca]

Hi Bob, Thanks for checking! TC1 is a variable that I added earlier for the new soil NOx scheme, and in my current code, TC1 is already commented out. It’s strange that the code you copied differs from mine.

Thanks @Ruijun-Dang. I was on a different branch where TC1 was still used. My bad

[yantosca]

Hi @Ruijun-Dang. So if you try adding this to the GFED section:

# Define a scale field for reducing NH3 over CONUS
111 GFED_NH3_SCALE   0.8/0.8/0.8/0.8/0.8/0.8/0.8/0.8/0.8/0.8/0.8/0.8 - 2000/1-12/1/0 C xy unitless * 1013 1 1

then this seems to work. When I used this, I get

               I         J.            scale factor                              array * scale factor
 ==>           53          43  0.80000001192092896        1.9367079757339347E-015
 ==>           54          43  0.80000001192092896        4.0430843390972404E-015
 ==>           55          43  0.80000001192092896        3.0206143354462638E-016
 ==>           56          43  0.80000001192092896        7.6552680354808959E-016
 ==>           59          43  0.80000001192092896        2.5269879428496747E-016
 ==>           60          43  0.80000001192092896        5.8963056592077322E-016
 ==>           61          43  0.80000001192092896        8.2212139130535175E-015
 ==>           62          43  0.80000001192092896        4.7806971259152073E-016
 ==>           63          43  0.80000001192092896        1.2267547981056939E-014
 ==>           64          43  0.80000001192092896        3.9574821949580880E-015
 ==>           65          43  0.80000001192092896        3.1229402009660097E-015
 ==>           66          43  0.80000001192092896        1.3602174779176424E-014
 ==>           67          43  0.80000001192092896        4.5447770528564785E-014
 ==>           68          43  0.80000001192092896        1.7423414083989194E-013
 ==>           69          43  0.80000001192092896        6.8743677775498441E-014
 ==>           73          43  0.80000001192092896        5.0654655802074555E-012
 ==>           74          43  0.80000001192092896        2.3496142055479046E-014
 ==>           75          43  0.80000001192092896        5.2863913415380469E-015
 ==>           78          43  0.80000001192092896        8.3277350500704566E-014
 ==>           79          43  0.80000001192092896        1.5355354493726849E-013
 ==>           80          43  0.80000001192092896        1.5282104608905001E-012
 ==>           81          43  0.80000001192092896        5.0939391700588453E-014
 ==>           83          43  0.80000001192092896        9.1297856510724852E-015
 ==>           84          43  0.80000001192092896        2.0050179536405719E-014
 ==>           87          43  0.80000001192092896        2.1914785345370648E-014
 ==>           88          43  0.80000001192092896        1.6253043073968589E-014
 ==>           89          43  0.80000001192092896        9.2924855246491672E-015
 ==>           90          43  0.80000001192092896        5.1564838150283649E-015
 ==>           91          43  0.80000001192092896        1.6699125161425314E-014

[yantosca]

You can also try with just an annual scale factor and see if it works. The annual one definitely works. Also you can refer to the mask directly (instead of to scale factor 666 NH3_reduction).

111 GFED_NH3_SCALE   0.8 - 2000/1/1/0 C xy unitless * 1013 1 1

Let me know if this works.

[Ruijun-Dang]

Hi @yantosca . I applied the annual scale factor as you suggested but it didn't work. The GFED_NH3 didn't change.

111 GFED_NH3_SCALE   0.8 - 2000/1/1/0 C xy unitless * 1013 1 1

And in my HEMCO.log, it looks like the HEMCO didn't get GFED_NH3_SCALE as a scale field...

   - Emit GFED species NH3     as model species NH3
      --> Will use scale factor:    1.00000000
      --> Will use scale field : none

[yantosca]

@Ruijun-Dang: Does the monthly scale factors work?

[Ruijun-Dang]

@yantosca Weirdly no. Both annual and monthly sf didn't work. Are you also using gfortran as the compiler?

[yantosca]

I was using the gfortran 10.2.0 compiler.

[yantosca]

This issue has been resolved by applying the updates in PR geoschem/HEMCO#218, as well as making some corrections to the HEMCO_Config.rc file.