[BUG/ISSUE] Possible bug in hco_geotools

[sdeastham]

While running GCHPctm with full debug flags, I came across an unexpected error just after

HEMCO/src/Core/hco_geotools_mod.F90

Line 1315 in c1c21cc

IF ( PRESENT( FldName ) ) THEN

. It seems that HCO_EvalFld is being called with HcoState%Grid%PBLHEIGHT%Val, but that array is still null (as set during the state initialization). Since HCO_EvalFld expects an already-allocated array, this results in an error when certain debug flags are on, and presumably some kind of memory corruption when they're not. Questions:

1. Am I missing the place where HcoState%Grid%PBLHEIGHT%Val gets allocated with non-zero size? In hco_state_mod's initialization the PBLHEIGHT allocation is with nx=0 and ny=0 (

   HEMCO/src/Core/hco_state_mod.F90

   Line 286 in 665c4a1

   CALL HCO_ArrInit ( HcoState%Grid%PBLHEIGHT, 0, 0, RC )

   ), which just results in a null pointer.
2. If I'm not missing this, should we be putting a call to HCO_ArrAssert( HcoState%Grid%PBLHEIGHT... immediately after the IF ( PRESENT( FldName ) ) THEN line? This would mirror the approach taken for the following options (IF ( PRESENT( PBLM) ) and so on).

[sdeastham]

Edit: I realize that it's possible that FldName will not be found, in which case perhaps we need to perform the HCO_EvalFld request with a dummy array and then only run HCO_ArrAssert if FOUND == .TRUE.?

[lizziel]

Thanks for reporting this. I think you are right. The target of pointer HcoState%Grid%PBLHEIGHT%Val must be allocated. Most of the HcoState%Grid members have their Val pointers point to either a State_Met or a State_Grid array so no problems there. HcoState%Grid%PEDGE%Val is different and is explicitly allocated. PSFC and ZSFC are allocated during the call to HCO_ArrAssert (this is somewhat hidden, which isn't ideal). But when HcoState%Grid%PBLHEIGHT%Val is used in HCO_EvalFld it is not allocated.

Like you say, using HCO_ArrAssert prior to setting HcoState%Grid%PBLHEIGHT%Val in HCO_EvalFld is needed. That subroutine is called in the second part of the if...else... block, but not in the first which is where Hco_EvalFld is called. Like as done for PSFC and ZSFC, it should be moved to before the if...else block to check the target dimensions and allocate if not yet associated.

I don't think it is needed to use a dummy array in HCO_EvalFld, but do think adding in a found argument is needed. If not found then HcoState%Grid%PBLHEIGHT%Val will not have changed and it will move to the next block to set it. If found, the target of HcoState%Grid%PBLHEIGHT%Val will have already been allocated and it can therefore be filled without issue.

It looks like this variable only impacts calculating the vertical dilution factor applied to distribute emissions into multiple vertical levels. But in GEOS-Chem the State_Met PBL height is passed so this fix shouldn't cause any diffs.

[lizziel]

We can put this bug fix in HEMCO 3.0.0 (for use in GEOS-Chem 13.0.0).

[jimmielin]

I think the PBLHEIGHT is asserted elsewhere, namely when the mixing height is passed in from State_Chm to HEMCO. I was actually wondering about this behavior too, and had to replicate this sequence of calls in HEMCO_CESM for HEMCO to work:

        HcoState%NX = my_IM
        HcoState%NY = my_JM
        HcoState%NZ = LM

        ! Pass Ap, Bp values, units [Pa], [unitless]
        ! later remove masterproc
        call HCO_VertGrid_Define(HcoState%Config,                &
                                 zGrid = HcoState%Grid%zGrid,    &
                                 nz    = HcoState%NZ,            &
                                 Ap    = Ap,                     &
                                 Bp    = Bp,                     &
                                 RC    = HMRC)
        ASSERT_(HMRC==HCO_SUCCESS)

        ! Point to grid variables
        HcoState%Grid%XMID%Val         => XMid   (my_IS:my_IE, my_JS:my_JE)
        HcoState%Grid%YMID%Val         => YMid   (my_IS:my_IE, my_JS:my_JE)
        HcoState%Grid%XEdge%Val        => XEdge  (my_IS:my_IE, my_JS:my_JE)
        HcoState%Grid%YEdge%Val        => YEdge  (my_IS:my_IE, my_JS:my_JE)
        HcoState%Grid%YSin%Val         => YSin   (my_IS:my_IE, my_JS:my_JE)
        HcoState%Grid%AREA_M2%Val      => AREA_M2(my_IS:my_IE, my_JS:my_JE)

        ! ...
        call HCO_CalcVertGrid(HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, HMRC)

        ! Pass boundary layer height to HEMCO (PBLm = PBL mixing height) [m]
        call HCO_SetPBLm(HcoState, FldName='PBL_HEIGHT', PBLM=State_HCO_PBLH, &
                         DefVal=1000.0_hp, & ! default value
                         RC=HMRC)

The call to hco_geotools_mod's HCO_SetPBLm will assert PBLHEIGHT correctly:

             ! Make sure size is ok
             CALL HCO_ArrAssert( HcoState%Grid%PBLHEIGHT, &
                                 HcoState%NX, HcoState%NY, RC )
             IF ( RC /= HCO_SUCCESS ) RETURN

[sdeastham]

Hi Haipeng, That’s the exact problem though – if a FldName is passed to HCO_SetPBLm (as is the case in the snippet you gave), then the first check in HCO_SetPBLm passes HcoState%Grid%PBLHeight%Val to a routine which expects an already-allocated array. This occurs first-thing in HCO_SetPBLm before and without running HCO_ArrAssert, resulting in the error I reported. Regards, Seb

[lizziel]

See my earlier response. Also, please try to follow the link before responding rather than straight from email. If responding via email then it includes the thread history which can get kind of long.

[jimmielin]

Hi Haipeng, That’s the exact problem though – if a FldName is passed to HCO_SetPBLm (as is the case in the snippet you gave), then the first check in HCO_SetPBLm passes HcoState%Grid%PBLHeight%Val to a routine which expects an already-allocated array. This occurs first-thing in HCO_SetPBLm before and without running HCO_ArrAssert, resulting in the error I reported. Regards, Seb

Thanks for clarifying! Indeed this looks like an oversight. I haven't read PBL_HEIGHT through HEMCO before (does GEOS-Chem - the CTM - do it using this field name? Presumably it reads PBLH from met data and passes that to HEMCO, not through this field name?); if this field is not found, it HCO_EvalFld will return without validating NX and NY, so we're all safe.

I agree that the array needs to be allocated further up the chain, probably in state initialization or in calculating the vertical grid. Thanks for reporting this!

If responding via email then it includes the thread history which can get kind of long.

<off-topic>Indeed, it's unbelievable that GitHub hasn't fixed this. It's an annoying bug.</off-topic>

[lizziel]

Yes, GEOS-Chem passes the data from State_Met to HEMCO to set PBLHEIGHT if FldName is not found. But in theory using data read in by HEMCO should work. Here's the relevant bit of code that passes the field to HEMCO.
https://github.com/geoschem/geos-chem/blob/c4b68994f9720f246f38a76489bff4913cf489ed/GeosCore/hcoi_gc_main_mod.F90#L2414-L2433

It's actually a bit circular since we now read met data via HEMCO. But the name of the PBL height read in and stored in HEMCO (PBLH) is not the same as the FldName passed to the subroutine to look for it in HEMCO (PBL_HEIGHT).

[jimmielin]

I wonder if we want to accept two ways of reading in the PBL height. Does this affect HEMCO standalone? For all other models there are other met fields that need to be passed into HEMCO anyway. If we don't have a reason to keep the PBL_HEIGHT field functionality, we can probably remove it since it's a hidden feature that needs maintenance.

[sdeastham]

Oh, wow - so if I understand this correctly (quite possible that I don't...) then the PBL_HEIGHT check actually "deliberately" fails in GEOS-Chem? So the error I'm seeing is semi-benign - Fortran is being passed a null pointer when it wants an allocated array, which is bad, but in reality doesn't matter because we don't actually use the passed variable. In either case though I guess the safest thing would be to move the HCO_ArrAssert to the start of HCO_SetPBLm - since every single possible option needs it anyway, and the routine actually causes the model to stop if none of the possible options are satisfied.

Note: I would cautiously second Haipeng's suggestion although I probably lack some understanding!

[lizziel]

The subroutine HCO_SetPBLm is not actually called in HEMCO, only in the HEMCO-GEOS-Chem interface, so I don't think this impacts HEMCO standalone at all. Also, the hard-coded PBL_HEIGHT field name is not in HCO_SetPBLm. That subroutine accepts any FldName. The string passed to it should be customized in the specific interface between HEMCO and a different model, such as is done for GEOS-Chem.

Using PBL_HEIGHT in the GEOS-Chem interface is probably historical baggage but we'd need to look more into the git history to know. It is optional so could be removed if the history reveals it's never going to be needed for GEOS-Chem. That would certainly make the intention of the call more understandable within GEOS-Chem.

[jimmielin]

Also, the hard-coded PBL_HEIGHT field name is not in HCO_SetPBLm. That subroutine accepts any FldName. The string passed to it should be customized in the specific interface between HEMCO and a different model, such as is done for GEOS-Chem.

Thanks for clarifying! Learn something everyday... I actually replicated this verbatim from the HEMCO-GEOS-Chem interface without understanding that the field can be arbitrarily named (and also not used). My understanding that for SetPblm:

  SUBROUTINE HCO_SetPBLm( HcoState, FldName, PBLM, DefVal, RC )
    TYPE(HCO_State),  POINTER                  :: HcoState    ! HEMCO state object
    CHARACTER(LEN=*), OPTIONAL, INTENT(IN   )  :: FldName     ! field name
    REAL(hp),         OPTIONAL, POINTER        :: PBLM(:,:)   ! pbl mixing height
    REAL(hp),         OPTIONAL, INTENT(IN   )  :: DefVal      ! default value

Probably only either FldName or PBLM should be present in the subroutine call. Otherwise, there may be a surprise "hidden feature" triggered by, for example, someone naming one of their HEMCO fields as the specified FldName, short-circuiting the input of PBL mixing height from the external model, which uses the PBLM attribute.

The implications of this is unclear; would there be a situation where HEMCO would want to use a different mixing height than the rest of the model?

[lizziel]

I added a fix for this to GEOS-Chem 12.8.1. It will be ported to the HEMCO repository when it gets updated for compatibility with 12.8.1 after that release.