seg fault in gfdl_fv_sat_adj

[evankalina]

While using the latest HAFS develop with the FV3_HAFS_v0_gfdlmp_tedmf_nonsst physics suite and HYCOM coupling, I get a segmentation fault during the initialization phase of the forecast job on Orion:

forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
hafs_forecast.x 00000000044B24D9 Unknown Unknown Unknown
libpthread-2.17.s 00002B0BCC7E65D0 Unknown Unknown Unknown
hafs_forecast.x 000000000367A902 fv_sat_adj_mp_wqs 1248 gfdl_fv_sat_adj.F90
hafs_forecast.x 0000000003676B86 fv_sat_adj_mp_fv_ 658 gfdl_fv_sat_adj.F90
hafs_forecast.x 000000000367387D fv_sat_adj_mp_fv_ 336 gfdl_fv_sat_adj.F90
hafs_forecast.x 00000000031EFFED fdlmp_tedmf_nonss 166 ccpp_FV3_HAFS_v0_gfdlmp_tedmf_nonsst_fast_physics_cap.F90
hafs_forecast.x 0000000003136840 ccpp_static_api_m 481 ccpp_static_api.F90
hafs_forecast.x 000000000250D491 fv_mapz_mod_mp_la 842 fv_mapz.F90

The seg fault occurs during the calculation of saturation vapor pressure in gfdl_fv_sat_adj.F90, which uses a table lookup:

es = tablew (it) + (ap1 - it) * desw (it)

You might think that the index it is outside the bounds of tablew, but it seems unlikely based on how it is calculated:

do i = is, ie

ap1 = 10. * dim (ta (i), tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = tablew (it) + (ap1 - it) * desw (it)
wqsat (i) = es / (rvgas * ta (i) * den (i))
it = ap1 - 0.5
! finite diff, del_t = 0.1:
dqdt (i) = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta (i) * den (i))

enddo

tablew is declared as allocatable though, so if it hasn't been allocated and then we attempt to use it, it would trigger a seg fault. If tablew is allocated, it is given a length of 2621.

BTW, this failure is identical to one reported by @BinLiu-NOAA in ufs-weather-model issue #718 back in July 2021. He fixed the problem by updating to a newer version of UWM. However, I'm using more recent code (last updated 1/26/22), so unless the bug was fixed and then reintroduced, I do not think it is a version issue.

In Bin's post, he suggests that the appearance of this bug may be sensitive to processor layout. I'm using 20x12 and could try some different ones.

If you have any thoughts about this error, please post.

[BinLiu-NOAA]

@evankalina I was wondering if you are running one of the standard HAFS application level regression tests, or another different configuration or different test cases. Typically, this kind of errors indicate the model is experiencing numerical instability issues. Also, if you rerun the forecast, will it fail at the same model integration time step?

[evankalina]

@BinLiu-NOAA Thanks for these ideas. I am running the 6-km H2DC-style configuration. The call to run_hafs.py looks like this:

h2dc: hafsv0p2a with fgat+3denvar+enkf WITH 6KM DETERMINISTIC GRID

${PYTHON3} ./run_hafs.py -t ${dev} ${FULL_CYCLE} 00L HISTORY
config.EXPT=${EXPT} config.SUBEXPT=${expt}
config.run_gsi_vr=no config.run_gsi_vr_fgat=no config.run_gsi_vr_ens=no
config.run_gsi=yes config.run_fgat=yes config.run_envar=yes
config.run_ensda=no config.ENS_SIZE=40 config.run_enkf=no
gsi.online_satbias=no
gsi.use_bufr_nr=yes
config.NHRS=12
config.run_emcgraphics=no ${scrubopt}
../parm/hafsv0p2a_regional_da_AL_6km.conf
../parm/hafs_hycom.conf

The failure always occurs during the adiabatic solver initialization (na_init=1), quite a bit before the first timestep. While adding some print statements to the code near the failure point, I see that the exact location of the segmentation fault moves around by a few lines (i.e., it doesn't always happen at es = tablew (it) + (ap1 - it) * desw (it)). The values going into the equation for es look normal. This makes me think it is a resource issue, but changing the layout to 40x30 (instead of 20x12) did not help. I tried setting the stack size to both unlimited and 6 GB, but this did not change the result either.