Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest file contains nans for the terrain in some cells with fine grid-spacing


New member
Hello all,

I am testing meshes that have a central region with spacings of around 1km, and my problem is that, for some of the cells in the central fine-resolution region, the terrain ('ter') value given by init_atmosphere when computing the file is NaN. However, the file is generated without throwing errors, and I only realised this was happening because using that file to generate the initial conditions produced a segmentation fault due to some zgrig(k, iCell) being NaNs.

The issue

In this picture, the terrain ('ter') variable of the file is plotted (zooming over the central fine-resolution region). There are 7 cells that have invalid (NaN) terrain values, colored dark grey.

I attach a zip folder ( with the input regional file; the namelist.init_atmosphere and streams.init_atmsphere files I used; and the output file that has been generated. The WPS_GEOG files are used were downloaded from the offical webpage, and from a message here in the forum that pointed to some extra files (landuse_30s and modis_landuse_20class_30s I think).

My theory

I have checked the source code and it seems to me that the HGT interpolation in mpas_init_atm_static.F does not check whether there are cells that can end up without terrain information. As I understand, the interpolation procedure consists of assigning to each MPAS cell the mean terrain of the topo_gmted2010_30s grid-points contained in that MPAS cell. However it does not take into account that it is possible that some MPAS cells, if they have fine grid-spacings, do not contain any of the topo_gmted2010_30s grid-points. Those cells end up with NaN values and no error is detected. I think that there should be a at least test for that situation, and it would be nice to implement a solution/algorithm that avoids this issue.

For now, in my source code, I added this snippet :

do iCell = 1,nCells

!MGB Test that no cells are empty of terrain data
if (nhs(iCell) == 0) then
call mpas_log_write('--- Cell $i was not filled by any topo', intArgs=(/iCell/))
call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Invalid topography dataset: at least one cell has NaN terrain ', messageType=MPAS_LOG_ERR)
call mpas_log_write(' This happens because the grid-spacing is finer than in the reference topo dataset.', messageType=MPAS_LOG_ERR)
call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write('Use finer reference topo or modify the code.', messageType=MPAS_LOG_CRIT)
end if

ter(iCell) = ter(iCell) / real(nhs(iCell))
end do
call mpas_log_write('--- end interpolate TER')

In fact, this also happens for other variables like the vegetation fraction (greenfrac) variable and soil category (isltyp), and the incorrect value is 0 instead of NaN, which is even worse because the next steps of the process will probably not fail but will be using very wrong data. So, I think tests should be added after all the interpolations that can potentially have this issue.


A solution I guess it would be to create a higher-resolution topography/landuse... option for the WPS input files, but I do believe something more elegant can be done that does not invlove using up much more storage space.

What do you suggest? Has someone found a work-around this issue? I believe some studies have used MPAS with resolutions around 1km, so I assume they managed to solve this, but I could not find anything here in the forum.

Any help will be appreciated,


    1.7 MB · Views: 0