wiki:Documentation/UserGuide/ExtractPixel

Version 6 (modified by dgoll, 4 years ago) (diff)

--

Extract a single pixel from a large restart file

Author: M. McGrath
Last revision: S. Luyssaert (2020/02/28)
Last check: 2020/04/20, D. Goll

Objective

This item should help to: (1) find the latitude and longitude of the pixel that crashed, (2) find the neighbouring pixels, and (3) make restart files containing a single pixel starting from a restart file covering a much larger domain.

Background

Sometimes it happens that you are running a simulation over a large spatial area and you have hundreds or thousands of grid points. If the simulation crashes and you know what pixel is causing the problem, it's nice to be able to use the previous restart file to debug. How do you convert a restart file with thousands of pixels into a single pixel, though? This far we have been using two methods: (1) re-run the model for a 3 x 3 pixels test case and (2) extract the restart information from a restart file covering a much larger domain.

Find the latitude and longitude of the pixel that crashed

The first problem is of course to know which pixel is causing the crash. The pixel number is easy to find: it is ji (sechiba) or ipts (stomate) in the code. If your code does not have the following lines of code then paste them right at the start of sechiba_init. This subroutine is only called the first time step so you will have to look in the out_orchidee file of the first time step of your run to find the latitude and longitude of the pixel.

! Debug
! It is good to leave this in here.  It is only written out once,
! it takes almost no time, and it's necessary for identifying a 
! problem pixel.
DO ipts=1,kjpindex
   WRITE(numout,'(A,I6,10F20.10)') 'pixel number to lat/lon: ',ipts,lalo(ipts,1:2)
ENDDO
!-

Compile the model and run until it crashes. If you know the ipts on which the model crashed you can now derive the latitude and longitude of pixel that is causing the crash.

Find the neighbouring pixels

Often the fasted way to reproduce a crash is to rerun the whole run but for 3x3 test case on 1 or 2 processors. A 3x3 test case is more likley to reproduce the crash than a 1 pixel test case because some files need to be interpolated at the start of the simulation. Interpolation of a 3x3 test case is more similar to the interpolation of a large spatial domain than the interpolation of a single pixel. A 5x5 test case would probably be even better but of course slower and more difficult to debug.

Setting up a 3x3 test case requires that you know the neighbouring pixels. If you are using a regular grid, the neighbours can be easily derived from the grid size, if you are using a zoomed grid, the neighbouring pixels can be derived from the "pixel number to lat/lon" information written to the out_orchidee file or you can use

ncdump -h XXX_stomate_history.nc 

If rerunning the whole simulation would take too much time your best option is to make restart files containing a single pixel starting from a restart file covering a much larger domain.

Make restart files containing a single pixel starting from a restart file covering a much larger domain

The tool to use is ncks, which stands for NetCDF Kitchen Sink. In English, the expression "I brought everything but the kitchen sink" means that I brought absolutely everything with me. ncks is a bit similar; it's a tool that does almost everything NetCDF related, according to the authors.

The first problem is figuring out which pixel you need. Sometimes you know the latitude and longitude, which, in theory, ncks can use as input to extract the pixel. However, in the ORCHIDEE restart files, the latitude and longitude are not stored as dimensions, which means ncks cannot do anything with this information. The restart files store the latitude and longitude as arrays, NAV_LAT and NAV_LON. Assuming that our problem pixel corresponds to 37.5E longitude and 55.5N latitude, we need to find out which pixel this corresponds to. This can be done with two commands.

ncks -H -C -v nav_lat EUROPE_19581231_stomate_rest.nc | grep 55.5

This commands prints out a bunch of lines, looking something like this:

y[19] x[0] nav_lat[988]=55.5 
y[19] x[1] nav_lat[989]=55.5 
y[19] x[2] nav_lat[990]=55.5 
y[19] x[3] nav_lat[991]=55.5 
y[19] x[4] nav_lat[992]=55.5 
y[19] x[5] nav_lat[993]=55.5 
y[19] x[6] nav_lat[994]=55.5 

We can see here that the y value of all these lines is equal to 19. Therefore, grid point 19 in the y direction must correspond to a latitude of 55.5 degrees. We can use a similar command for the longitude.

ncks -H -C -v nav_lon EUROPE_19581231_stomate_rest.nc | grep 37.5

From this, we see that the x value is always constant, and equal to 49.

y[0] x[49] nav_lon[49]=37.5 
y[1] x[49] nav_lon[101]=37.5 
y[2] x[49] nav_lon[153]=37.5 
y[3] x[49] nav_lon[205]=37.5 
y[4] x[49] nav_lon[257]=37.5 
y[5] x[49] nav_lon[309]=37.5 
y[6] x[49] nav_lon[361]=37.5 

Therefore, we know that (x,y)=(49,19) corresponds to the pixel that we want. With this information, it's just a matter of extracting the correct pixel into our new restart file.

ncks -a -O -d y,19 -d x,49 EUROPE_19671231_stomate_rest.nc new_stomate_restart.nc

If you want to extract a rectangular 2 by 2 domain use

ncks -a -O -d y,19,20 -d x,49,50 EUROPE_19671231_stomate_rest.nc new_stomate_restart.nc

Tell libIGCM where to find the restart files by completing the config.card. Do this for SRF, SBG and OOL

#D-- SRF - SECHIBA
# Use this section to start from restart file for SRF if OverRule=n
# WriteFrequency is not used any more for ORCHIDEE
# Output files are now managed in sechiba.card
[SRF]
WriteFrequency=""
Restart=n
RestartDate=...
RestartJobName=...
RestartPath=...