chem: add scalar guards for p_hcl and p_ch3cl in add_anthropogenics t…

[AlphaBetaGammaChi]

Fixes a runtime SIGSEGV crash during emissions_driver when running chemical mechanisms that do not include anthropogenic chlorine (such as CRI-MECH with emiss_opt = 19 or emiss_opt = 20).

Cause
In chem/module_emissions_anthropogenics.F, assignments for p_hcl and p_ch3cl were placed in the generic anthropogenic emissions loop without param_first_scalar bounds checks.
For mechanisms where e_hcl and e_ch3cl are not present in the package, referencing chem(..., p_hcl) and emis_ant(..., p_e_hcl) causes an invalid memory reference and crashes on timestep 1:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace:
__module_emissions_driver_MOD_emissions_driver
chem_driver_
solve_interface_

[weiwangncar]

The regression results:

Test Type              | Expected  | Received |  Failed
= = = = = = = = = = = = = = = = = = = = = = = =  = = = =
Number of Tests        : 23           24
Number of Builds       : 60           57
Number of Simulations  : 158           150        0
Number of Comparisons  : 95           86        0

Failed Simulations are: 
None
Which comparisons are not bit-for-bit: 
None

[weiwangncar]

@dudhia Without this, the model can abort as it has been shown - using array indices outside do loops. If it needs to be that accurate, there are other ways to fix it.

[saneku]

I introduced

area=(dx/mapfac_mx(i,j))*(dy/mapfac_my(i,j))

And I could achievemass-balancee closure with it.

Using
area=dx*dx

does not provide mass-balance closure, as it is a less accurate area calculation.

[weiwangncar]

@saneku Can you move the area calculation to a do loop where the accurate area information is needed? For example area can be defined inside the i,j do loop that begins on line 455, for example. Or you can declare a local 2D array to store the area array to be used later?

[saneku]

I agree with @weiwangncar. The map-factor-aware area calculation should be preserved, but moved inside the i,j loops, where the indices are defined:

! For volcanic ash transport, vash variables are in mixing ratio here.
do j = jts, jte
do i = its, ite
if (erup_end(i,j).gt.0) then
area = (dx/mapfac_mx(i,j)) * (dy/mapfac_my(i,j))

    so2_mass = 1.5e4*3600.*1.e9/64./area

The area argument should then be removed from the preceding diagnostic message. Its format should also be changed from 2f15.3 to f15.3:
write(message,'(" ADJUSTED ASH HEIGHT: ",f15.3)') emiss_ash_height

[weiwangncar]

@AlphaBetaGammaChi Would you like to make this change in this PR?

[saneku]

@weiwangncar and @AlphaBetaGammaChi i have noticed that this correction is already implemented in
#2353

[weiwangncar]

@saneku Thanks for the alert. We will review PR-2353 and close this one afterwards.

[AlphaBetaGammaChi]

Hi @weiwangncar,

Thank you! Just to flag: PR #2353 addresses the emissions_driver.F area calculation, but the other half of this PR — the scalar guards for p_hcl and p_ch3cl in chem/module_emissions_anthropogenics.F — is a separate, independent bug fix that is not covered by #2353 I believe.

Without those guards, any CRI-MECH run (chem_opt=601 with emiss_opt=19 or emiss_opt=20) will still crash I believe with a SIGSEGV on timestep 1 due to unguarded array references for species not present in the mechanism.

Would it be possible to either:

1. Keep this PR open for the module_emissions_anthropogenics.F fix (I can revert the emissions_driver.F change so it doesn't conflict with Fixed the volc diags bounds crash #2353), or
2. Incorporate the scalar guards into Fixed the volc diags bounds crash #2353?

Happy to proceed either way!

[weiwangncar]

@AlphaBetaGammaChi I suggest to go with your proposal #1. Thanks for the reminder for the other fix.