Error interpolating SST data using multiprocessing with METGRID

[greent12]

I am trying to initialize WRF with GFS 0.25 degree data and using an external SST source (RTG 1/12 deg data). The GFS and SST data can be obtained here:

GFS: Index of /pub/data/nccf/com/gfs/prod
SST: Index of /pub/data/nccf/com/nsst/prod

I've made one change to the METGRID table for the SST variable. The original entry was:

========================================
name=SST
interp_option=sixteen_pt+four_pt
fill_missing=0.
missing_value=-1.E30
flag_in_output=FLAG_SST
========================================

The new entry is:

========================================
name=SST
interp_option=sixteen_pt+four_pt+average_4pt
masked = land
interp_mask = LANDSEA(1)
fill_missing=0.
missing_value=-1.E30
flag_in_output=FLAG_SST
========================================

When I run metgrid with one processor, the program runs to completion. However, when I run the program with multiple processors, I get a segmentation fault.
The stack trace points to the sixteen_pt interpolation function in the metgrid source code. Additionally, when I change the interpolation option, even to the most basic "nearest_neighbor" setting, the segmentation fault still occurs when using more than one processor. If I change the 'interp_mask' from LANDSEA(1) to LANDMASK(1) in order to use the mask from geogrid, the program works with any number of processors.

The SST data I am using does not come with a LANDSEA mask, and I am aware that the LANDSEA mask I am asking for is from the GFS data. I noticed that 'interp_module.f90' in the metgrid source code declares the 'array' and 'masked_array' with the same dimensions. In this case, the 'mask_array' would correspond to the dimensions of the GFS data, and the 'array' would correspond to the SST data (different dimensions than the GFS data). Is it possible that this is causing some issues with indexing out of bounds of one of these arrays?

I am attaching my namelist and METGRID.TBL files. The variable tables I used for ungribing the GFS and SST data are the 'Vtable.GFS' and 'Vtable.SST' files in the WPS/ungrib/Variable_tables directory. I've been able to reproduced this error on different systems, so i'm hoping this will be easily reproducible.

The SST data in the above link is only avaialbe for the past two days, so adjustments to my namelist for dates will be needed.

I am using WPS v4.2

 

[Ming Chen]

I am suspicious that the problem is attributed to LANDSEA mask. LANDMASK is derived from WPS static data, while LANDSEA is from the external foricng data. For your case, you use RTG SST as input, and I believe that you should create a land-sea mask specifically for RTG SST. Let's suppose the land-sea mask field for RTG data is named as LANDSEA-RTG, then in etgrid.tbl, we need to add the entry like:

name=SST
interp_option=sixteen_pt+four_pt
masked = land
fill_missing=0.
missing_value=-1.E30 ==> This is the value in RTG that is assumed to represent missing data
interp_mask=LANDSEA-RTG(maskval). ==> RTG data points will not be used in interpolation if the corresponding point in LANDSEA-RTG is maskval
flag_in_output=FLAG_SST

Please try and let me know whether it work.

 

[greent12]

Hi @Ming Chen ,

Thanks for your response. I have created a landsea mask for the 1/12 degree RTG SST data. I did this by doing a nearest neighbor interpolation of the water category in the modis 30s landuse data using cyclic BCs in longitude. I have attached a figure of what that mask looks (RTG_landmask.png). This mask now has the same dimensions as the RTG SST data.

I then wrote this mask out to intermediate file format. I've attached this file as well. I believe the intermediate file is written out correctly, as the 'rd_intermediate.exe' utility runs successfully and prints out the information consistent with what I wrote to it. The prefix for my intermediate files are "LANDSEA_RTG" and the field name inside the intermediate file is "SSTLSEA".

Metgrid successfully ran, however, inspection of the SST field shows a large swath of SSTs across the ocean in my domain (SST.png) that got filled with missing values (0.). Oddly, the other water points seemed to have been interpolated/assigned values correctly. Inspection of what was written out by metgrid for the SSTLSEA mask (SSTLSEA.png) shows that the mask is being correctly represented in my domain. Any idea of what may be going on here?

Note: I've expanded my domain compared to the original namelist I attached and changed to a lat-lon projection, so I am attaching the new namelist and my new METGRID.TBL file. Note: I did try this with a lambert projection as well, and the same problem persists, so it seems to be independent of projection used.

 

[Ming Chen]

Hi,

Thank you for uploading the plots. Let's try one more option:

please delete the following section from your METGRID.TBL:

===========================================
name=SSTLSEA; from_input=LANDSEA_RTG
interp_option=nearest_neighbor
fill_missing=-1.
fill_lev=200100:LANDMASK(1)
===========================================

Then rerun metgrid.exe. Let me know how it works.

 

[greent12]

Hi Ming,

I tried this and still results in the swath of missing SSTs across the gulf.

However, I think I did found the fix for this. In my METGIRD.TBL I used 'interp_land_mask = SSTLSEA(1)' instead of 'interp_mask = SSTLSEA(1)'.

Old entry:
========================================
name=SST
interp_option=sixteen_pt+four_pt+average_4pt
masked = land
interp_mask = SSTLSEA(1)
fill_missing=0
missing_value=-1.E30
flag_in_output=FLAG_SST
========================================

New entry:
========================================
name=SST
interp_option=sixteen_pt+four_pt+average_4pt
masked=land
interp_land_mask = SSTLSEA(1)
fill_missing=0
missing_value=-1.E30
flag_in_output=FLAG_SST
========================================

I've attached what the SST field from metgrid looks like. (SST_correct.png). The field looks blocky because i've asked for 0.25 degrees grid spacing and also asked panoply not to interpolate.

My question now is, what is the correct usage of the 'interp_mask', and 'interp_land_mask', and why would changing this entry in the table change the results from metgrid? It is true that the mask that I created is intended to be the mask for 'land', but when should you use 'interp_mask'.

Thanks,
Tyler

 

[Ming Chen]

Tyler,

Many thanks for the update, which saves me lots of time to look into this issue.

Regarding your question of interp_land_mask and interp_mask, my understanding is that, INTERP_LAND_MASK specifies the name of the field to be used as an interpolation mask when interpolating to water points, while INTERP_MASK gives the name of the field to be used as an interpolation mask for interpolation to any points (i.e., not distingquish land and water points).

For your case, we only need to interpolate SST to water point, which is why INTERP_LAND_MASK should be used.

 

[greent12]

Hi Ming,

That makes sense to me. However, I see entries in the default 'METGRID.TBL.ARW', for example, 'SEAICE' that use 'interp_mask' and 'masked=land'. If it is true that 'interp_land_mask' should be used when interpolating to water points using the 'masked=land' key, why would this field (and many others in the default table) not use 'interp_land_mask'?

 

[greent12]

Hi Ming,

I just wanted to provide you an update on some findings.

From a few previous posts, I noted how there was a large swath of missing SST values over the gulf when trying to use 'interp_mask' and 'masked=land'. I tracked down the reason that this was happening. This was an artifact caused by not adding enough decimal points to my 'dx' when writing out the intermediate file for the SST landmask. This is a 1/12 degree dataset, making the dx ~ 0.083333 (repeating). In the routine 'get_interp_masks', the conditional statement at the end, which tries to determine if it is a global dataset or not, was not being met due to lack of precision: abs(4320 x 0.083333 - 360) = 0.00144 which is not less than the 0.0001 threshold. Because of this, the three gridpoint pad was not being applied to the mask array. When this mask was being passed in to the interpolation routine, the dimensions did not match what the routine was expecting (three gridpoint pad on either side in longitude dimension), causing weird behavior with indexing the mask. The dimensions of the mask are expected to be the same as the array being passed in, which in this case is the 1/12 SST data that was correctly determined to be a global dataset, and therefore had the three gridpoint pad. I added some more precision to my dx when writting my intermediate file, which fixed the issue.

In regards to using 'interp_land_mask", I noticed that the return code from the 'storage_get_field' call inside the 'interp_met_field' routine was 1 (could not find the field?). Because of this, the calls to the interpolation routines were not passing in the mask as an optional argument, so the mask was not even used. This did not manifest itself in the plot above called 'SST_correct.png' because the SST field has values over land. So even though I was intending for the interpolation routines to not use data I deemed to be over land in the SST dataset, the fact that it was using those values for the interpolation wouldnt cause the interpolated value to look suspicious. If the SST field would have had missing values over land, then this would have given weird values near the coastlines.

As I think I have solved my problem for good here, I am still curious as to what the purpose of having the options 'interp_mask' and 'interp_(land/sea)_mask' in the metgrid table is? I'm just curious to know so that in the future I can use these keywords appropriately.

Thanks!