CMAQ5.5 Negative or Zero [H+] concentration specified in HLCONST

Hello.
I’m using the released version of CMAQv5.5. Unfortunately, I’m dealing with this error message:
Negative or Zero [H+] concentration specified in HLCONST
PM3EXIT: date&time specified as 0
Date&time specified as 0

I don’t know how to solve this error. Please, let me know if you need more information
Attached are the CCTM script and the log file.
run_cctm_Dominio1_FONDECYT.csh (39.1 KB)
CTM_LOG_005_v55_intel_Dominio01_Fondecyt_cb6r5_ae7_aq_m3dry_20190111.txt (89.3 KB)

Hi Ernesto,
Please recompile in debug mode, rerun the model, and report back whether the abort occurs in the same place. If it crashes, please provide the stack trace.

What landuse scheme and land-surface model did your WRF simulation use?

Dear @cgnolte . Thanks for your quick response.
I recompiled in debug mode and rerun the model, but it didn’t work. Unfortunately, I don’t understand what’s the error in the log files attached.

About the landuse scheme, I used Thermal Diffusion (option 1 in namelist.input for WRF simulation)

Please, if you help me solve this issue.

Best
run_cctm_Dominio1_FONDECYT_debug.csh (39.9 KB)
CCTM_36005431.txt (45.1 KB)
CTM_LOG_000.v55_intel_Dominio01_Fondecyt_cb6r5m_ae7_aq_m3dry_20190111.txt (106.5 KB)

The stack trace at the end of the CCTM_36005431.txt file indicates you are getting a divide by 0 error in m3dry. This is not my area of expertise, but there is likely an inconsistency between the land-water mask and the leaf area index, where some land grid cell (LWMASK=1) has a 0 value for LAI.

To verify whether that is the problem, you could run in a debugger, or just add a write statement prior to line 666 in m3dry.F, something like:

IF ( LAICR .EQ. 0.0 .OR. DIF0(L) .EQ. 0.0 ) THEN
   WRITE(*,*) C,R, LAICR, DIF0(L)
END IF

As for how to fix this issue, if it is possible to rerun WRF you might do so with a different land-surface scheme, such as PX. You can try modifying m3dry to avoid these error conditions, but from reading the WRF documentation it seems the 5-layer thermal diffusion scheme is simplistic, with soil temperature only and without soil moisture, so I am not sure whether this can work with CMAQ.

Dear @cgnolte
Thanks again for your reply.
I added the statement you sent


but it crashes again at the same line (see the error in log file).

I don’t know why this version of CMAQ has this kind of issue. I have used the same land-surface scheme with version CMAQ5.3 and it worked well. I suspect some mistake in the new m3dry.F file

The underlying issue appears to be that you have grid cells where the MCIP (WRF) LWMASK field equals 1 (i.e. land) but the LAI field (and likely VEGF field) equals 0. This is not a situation expected in the m3dry code. To determine where and why this happens, you’d probably have to dig into the WRF preprocessing and the WRF LSM code to determine how LWMASK, LAI, and VEGF are handled.

In previous versions of the code, the if/then block starting on line 516 defined the “water” case as

IF ( ( NINT(GRID_DATA%LWMASK( c,r )) .EQ. 0 ) .OR. ( vegcr .EQ. 0.0 ) ) THEN  ! water

In CMAQv55, this block was updated to only rely on LWMASK to determine whether a grid cell should be treated as water (in which case LAI isn’t needed) or land (in which case non-zero LAI is needed):

 IF ( NINT(GRID_DATA%LWMASK( c,r )) .EQ. 0 ) THEN  ! water

We would have to talk to the developer why this change was made. In any case, if my assumption is right that you have (unrealistic, or at least unexpected by m3dry) grid cells with LWMASK = 1 but LAI = VEGF = 0, the previous code would have protected you because VEGF = 0 would have branched these grid cells to the water case where LAI is not needed but the new code will branch these cells to the land case where the LAI = 0 value will cause a problem.

3 Likes

Dear @hogrefe.christian
Thanks for your reply for this post.
I changed the line 516 as you mentioned and it worked. So, now I’m confused about what statement must be. At least, the issue was identified and the solution is shown for other users.
Best

Hello Ernesto,

to help us better understand why you encountered a situation where LWMASK was 1 but VEGF and LAI were 0, could you please share which landuse scheme you used in your WRF simulations and whether your wrfout files contained variables VEGFRA and/or LAI?

Thanks,

Christian

Dear @hogrefe.christian
The landuse scheme is “Thermal Diffusion”.
In this link you can download one of the wrfout files to view the VEGFRA and LAI variables and their information. Hope this could help you to understand this issue.

Hello Ernesto,

thank you for posting the wrfout file. It seems your landuse scheme (not land-surface model) was MODIFIED_IGBP_MODIS_NOAH.

Looking at the LWMASK, LAI, and VEGFRA fields, the problem likely originates over Antarctica and the Atacama desert with zero LAI values. If you had used the PX LSM in WRF, it would have prevented the divide-by-zero issues in the CMAQ m3dry deposition calculations because it assigns non-zero but very small minimum vegetation fraction and LAI values even to the snow/ice and barren/sparsely vegetated MODIS landuse categories (see this piece of WRF code).

Hello Ernesto,

after some discussions with the developer, rather than reverting the IF statement on line 516 in CMAQv5.5 to the earlier version, it might also work to change lines 665 and 666 from

           rstom = MET_DATA%RS( c,r ) * dwat / dif0( l )
 &               + 1.0 / ( heff_ap / 3000.0 + 100.0 * meso( l ) ) / laicr

to

           rstom = MET_DATA%RS( c,r ) * dwat / dif0( l )
 &               + laicr / ( heff_ap / 3000.0 + 100.0 * meso( l ) )

If it works, it would be preferable since the code would then treat ozone deposition to grid cells in Antarctica (snow/ice land use) and desert areas as land-based deposition which would have higher rates than deposition to water.

1 Like

Hello,

I am using the released version of CMAQv5.5.
I also encountered this bug, so I changed lines 665 and 666 in m3dry.F.
It worked!

Then, I tried to compile the “large” version of CMAQv5.5 to run the ISAM version and ensure that variables can exceed 2048 to 16384.
However, I encountered the same bug again. Even after modifying lines 665 and 666 in m3dry.F, unfortunately, it didn’t work.

Are there any ways to resolve this problem?

It is unexpected that changing I/O API from the standard to the “large” version would have any impact on whether this portion of the m3dry code encounters this particular error.

Could you please compile in debug mode, post the log files, verify that your latest build includes the modified version of m3dry.F with the changes to lines 665 and 666 rather than the original file copied from the code repository (for debugging purposes, you could add some print statement to the modified file so that your CCTM log file confirms the use of the modified file), and share which WRF LSM and landuse scheme you were using? As discussed above, this error likely can also be avoided by using the PX LSM in WRF, which would also provide greater consistency with the dry deposition calculations in CMAQ.

Hi everyone,

I can solve this issue by reduce the core.

But I do not konw why it worked.

Best regards,

Haofan

1 Like

Hi @hogrefe.christian,

I ran into the same problem as described in this post with the LWMASK = 1 but VEG = 0 barren ground grid cells causing problems for hlconst.F. The land scheme in my WRF meteorology is the same as the original poster.

The code modification suggested to rstom did not solve the problem. Instead, the first suggestion of changing to IF ( ( NINT(GRID_DATA%LWMASK( c,r )) .EQ. 0 ) .OR. ( vegcr .EQ. 0.0 ) ) THEN ! water on line 516 of m3dry.F did work.

While this allows the model to successfully run, it makes an incorrect assumption about the land surface type for certain grid cells. In the domain that I’m running over, some of the barren land grid cells are inland and should not be lumped into the water category. Instead, I’m wondering if line 420 of m3dry.F can be modified from
IF ( ( NINT(GRID_DATA%LWMASK( c,r )) .NE. 0 ) .AND. ( vegcr .GT. 0.0 ) ) THEN ! land
to
IF ( NINT(GRID_DATA%LWMASK( c,r )) .NE. 0 ) THEN ! land
so that it also relies only on LWMASK. This would mean that the calculation of hplus and heff in m3dry have matching conditional statements. When I tested this in CMAQ, it worked but I wanted to check that this was an acceptable change to the model code before running the whole simulation period.

Thanks!
-Alicia

Hello Alicia,

which land-surface model (LSM) are you using in your WRF simulations?

The m3dry dry deposition code is designed for consistency with the PX LSM in WRF, and within that WRF PX LSM / CMAQ M3DRY framework, a case with LWMASK = 1 but VEGF = 0 (or LAI = 0) is not a valid combination. As mentioned in my initial responses to Ernesto, the WRF PX LSM assigns nominal non-zero VEGF and LAI values to all non-water land use categories.

Fundamentally, m3dry distinguishes between deposition to water (in which case the rc / bulk surface resistance term consists of only one pathway, i.e. the resistance of deposition to water) and deposition to land (in which case the rc term consists of parallel pathways to the non-vegetated soil, vegetated soil, cuticles, and stomata present in the grid cell) - see Clifton et al. (2023) and Galmarini et al. (2021) for some reviews of dry deposition schemes including m3dry.

While the code changes you suggest would functionally work, sending the m3dry dry deposition calculation into the “land” branch even for grid cells where VEGF (and/or LAI) is zero violates the assumption that deposition to vegetated soil, leaf cuticles, and leaf stomata all occur in such “land” grid cells. If VEGF is zero in a “land” cell, I’m not clear what the computation of the vegetated soil, cuticular, and stomatal resistances and the derivation of rc from those components would actually represent in your case. It would depend on understanding how the WRF LSM calculated the variables used in the computation of these component resistances and then passed these through MCIP to CMAQ.

Conceptually, you could probably hack the m3dry code to handle land cells with zero VEGF and/or LAI in the same way heat and moisture flux calculations were handled for such cells in the WRF LSM code, possibly by adding branching for rc resistance calculations for “land” cells (based on LWMASK) that add a new case where VEGF is zero and computes rc as only deposition to non-vegetated soil and the existing case where VEGF is assumed to be non-zero and rc is computed from deposition to non-vegetated soil, vegetated soil, leaf cuticles, and leaf stomata. Doing so would likely involve some deep-dive into the code and require some testing of dry deposition velocities calculated for such grid cells.

1 Like

Hi Christian,

Thanks for the super detailed response!

The landuse scheme in the WRF simulations was MODIFIED_IGBP_MODIS with the NOAH LSM, which I think might have been what Ernesto was also using. I have no idea how LWMASK=1 VEG=0 grid cells were treated in WRF for things you mention such as heat flux.

In the Clifton et al. 2023 paper you included, based on equation 76 (if I’ve read it correctly) the m3dry calculation of surface resistance rc can be reduced down to rc = rg (ground resistance) in the case that fraction of vegetation cover is 0. In the m3dry.F module, the surface resistance rc is calculated in the code block that starts on line 516 with the conditional IF ( NINT(GRID_DATA%LWMASK( c,r )) .EQ. 0 ) THEN ! water and continues on line 576 with ELSE ! land. And following the code in this section, it seems like the calculation for surface resistance does reduce down to just ground resistance when LAI/VEGCR are 0. This is the original code, not what I’ve modified. So CMAQ v5.5 is already able to calculate surface resistance and then deposition velocity based only on the LWMASK setting.

The section that I wanted to modify starts on line 420. From what I understand of the code, this section determines canopy wetness and H+ concentration for vegetated land and wetbulb temperature for water. So modifying the conditional statement on line 420 from to rely on just the LWMASK should mean that all land grid cells, even the bare ground cells, have a canopy wetness and [H+] calculated in them. In the case that vegcr = 0, canopy wetness delta = 0. But I wasn’t sure I wasn’t understanding this section correctly, which is why I originally posted my question. Is my understanding of this section of the m3dry.F module correct?

Thanks!
-Alicia

Hello Alicia,

thanks for your response.

I had another look at the code. Your interpretation of the code section from lines 516 to 679 (specifically the “land” section starting on lines 576) is correct. While PX LSM and therefore CMAQ M3DRY do not expect any land cell where VEGF zero, and while rcut, rwetsfc, and rgndc are all calculated even though there is no vegetation, they ultimately don’t impact the rci calculation on lines 670 - 675 since vegcr is 0, so rci (inverted rc) is indeed just calculated from the deposition to soil (dry, wet, and/or snow-covered). My suggestion on Friday that additional coding might be needed in that section to handle this situation was not correct.

You are also correct that by modifying line 420 to ensure that the calculation of moisture parameters affecting the soil and cuticular resistances occurs for all land cells (even when VEGF is zero) will then allow the code to compute these resistances (for your scenario, most importantly the effect of moisture on soil resistance) in lines 516 - 579.

That said, there may be parts of the soil resistance calculations (e.g. the rgnd equations in lines 635 and 639) where the assumptions made in the WRF PX LSM / CMAQ M3DRY framework and the WRF NOAH LSM are inconsistent. Specifically, there may be inconsistencies about the depth of the first soil layer for which WRF soil moisture is calculated (variable SOIM1), which in turn would mean that the rgnd equation may need to look different for WRF NOAH LSM SOIM1 than WRF PX LSM SOIM1. This would be similar to this discussion about windblown dust and its dependency on soil moisture assumptions in the WRF LSM. Similarly, there is a comment on line 625 that “Canopy height is assumed to be 10 * z0 according to PX-LSM“. I do not know whether this assumption would be consistent with the assumptions in the NOAH LSM.

Ultimately, the code modification you envision should functionally work to have M3DRY compute dry deposition velocities for all grid cells for your met fields, but some aspects of these computations may be affected by inconsistencies in assumptions between WRF PX LSM and WRF NOAH LSM.

Hi Christian,

Your point about checking for inconsistencies between WRF PX LSM and WRF NOAH LSM absolutely makes sense. I’ll have to dig into this a bit more before deciding how to proceed. And thank you for confirming that the code change would at least functionally work.

Thanks for all your help!
-Alicia