[RegCNET] Possible small bug in postproc.f

[Travis O'Brien]

Dear RegCM Users,

I have come across an issue in postproc.f that can sporadically cause 
segmentation faults (and similar memory access errors): in the 
subroutine CALCVD(), which calculates the horizontal vorticity (named 
'VOR' in the netcdf files) and divergence (DIV), the algorithm for 
computing these quantities attempts to access parts of memory outside of 
the bounds of the U and V input arrays (see below for a description).

This is problematic because sometimes the subroutine can access the 
out-of-bound memory and get whatever value is stored there (which often 
seems to be 0, but I don't think that's guaranteed and may depend on the 
compiler).  However, sometimes if a person is unlucky, the memory is 
outside of the bounds of the program, which can cause the operating 
system to kill the program without any warning or error message other 
than the cryptic "Segmentation Fault" error.

I don't believe that this matters much on the data analysis side of the 
fence, because the DIV and VOR will only be affected by this at grid 
cells that are on the edge of the domain, and I think that these values 
shouldn't even be specified since these quantities are derivatives and 
so shouldn't have values at the domain edges.  Does anyone have any 
thoughts on this?

I propose to fix the issue by changing the do loops so that they 
terminate at nx1-1 and ny1-1, and then setting the edges of the arrays 
to be whatever the missing_value attribute is set to be in the netcdf 
files.  My implementation of this fix is attached as a patch to the 
version of postproc.f that is supplied in the current version of RegCM.

Cheers,
Travis O'Brien

Description of the issue--
in postproc1.param, nx1 and ny1 are defined to be the same as nx and ny, 
which is the size of the U and V arrays in postproc.f:
parameter (ny1=ny,nx1=nx)

in postproc.f, there is a do loop that goes from 1 to nx1/ny1, which are 
the upper bounds of the U and V arrays.  A piece of code tries to access 
i+1 and j+1, which are beyond the bounds of U and V when i=nx1 or j=ny1:
subroutine CALCVD(....)
    .
    .
real   u(nx,ny,nz)
    .
    .
do j=1,ny1
do i=1,nx1
    .
  u2=u(i+1,j  ,k)/dmap(i+1,j  )
    .
enddo
enddo

-- 
Travis A. O'Brien
Graduate Student Researcher
Earth and Planetary Science Dept.
UC Santa Cruz

tobrien at ucsc.edu
(831) 459-3504

-------------- next part --------------
A non-text attachment was scrubbed...
Name: postproc.f.patch
Type: text/x-patch
Size: 691 bytes
Desc: not available
URL: <https://lists.ictp.it/pipermail/regcnet/attachments/20080129/81fa5a0e/attachment-0002.bin>