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?

The first problem is of course to know which pixel is causing the crash. The pixel number is easy to find: it is ipts 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 lat and long 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
!-

The answer lies in the tool 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