[ferret_users] Problem with external function

[James Caveen <james_caveen@xxxxxxxxxx>]

Dear ferreters,

I have been trying to solve a problem with an external function for
quite some time now and I finally decide to turn to the Ferret
Community.  I realize that  this message is somewhat long but here I go
anyway.

 
I have written a Ferret  external function to reconstruct a 3D grid that
is stored using compression by gathering as described in Section 8.2 of
the CF Convention :
(http://www.cgd.ucar.edu/cms/eaton/cf-metadata/CF-current.html#gath)

and everything works fine when I use one NetCDF file.  However, if I try
to open a second file, Ferret's limit for irregular axis is exceeded;
my data is stored in a 295000 element vector with uneven spacing. 

Following is a ncdum -h of one of my files:


netcdf COMP_TEST.ROM_TEMP_copie1 {
dimensions:
        XGULF = 236 ;
        YGULF = 150 ;
        ZGULF = 73 ;
        ZGULFedges = 74 ;
        TGULF = UNLIMITED ; // (4 currently)
        IWET = 295000 ;
variables:
        float XGULF(XGULF) ;
                XGULF:point_spacing = "even" ;
        float YGULF(YGULF) ;
                YGULF:point_spacing = "even" ;
        float ZGULF(ZGULF) ;
                ZGULF:positive = "down" ;
                ZGULF:point_spacing = "uneven" ;
                ZGULF:edges = "ZGULFedges" ;
        float ZGULFedges(ZGULFedges) ;
        float TGULF(TGULF) ;
                TGULF:units = "HOURS since 1901-01-15 00:00:00" ;
                TGULF:time_origin = "15-JAN-1901" ;
        int IWET(IWET) ;
                IWET:compress = "XGULF YGULF ZGULF" ;
                IWET:point_spacing = "uneven" ;
        float TEMP(TGULF, IWET) ;
                TEMP:fill_value = 9.96921e+36f ;
                TEMP:missing_value = 9.96921e+36f ;
                TEMP:long_name = "TEMP" ;
                TEMP:units = "celcius" ;

and here is an example script where I use my nc_rar() external function:

set memory/size=500
use COMP_TEST.ROM_TEMP.nc
let a = x[gx=xgulf]*0 + y[gy=ygulf]*0 + z[gz=zgulf]*0 + t[gt=tgulf]*0
let b = nc_rar(a,temp])
shade b[l=1,k=1]

        

So here I go with my Question 1:
All my files use the same IWET axis.  Is there a way that I can tell
Ferret that the IWET axis of file two is identical to the IWET axis of
file one and thus force Ferret not to allocate extra space for the
second IWET axis?


Being somewhat naive, I attempted to fool Ferret by pretending the IWET
axis is evenly spaced as follows:

 
        int IWET(IWET) ;
                IWET:compress = "XGULF YGULF ZGULF" ;
                IWET:point_spacing = "even" ;

Now I can open two,three,... NetCDF files without encountering the
memory limit problem but my external function doesn't work anymore (all
indices returned by ef_get_coordinates() are wrong).  

So, here is my second question:
Could someone briefly explain to me what pre-treatment Ferret does to
axis variables before passing them on to an external function?


Finally, being somewhat stubborn, I decided to convert my IWET axis into
a regular variable (WET) and to define an evenly spaced axis XWET which
is common to all my variables:   

netcdf COMP_TEST.ROM_TEMP_copie1 {
dimensions:
        XGULF = 236 ;
        YGULF = 150 ;
        ZGULF = 73 ;
        ZGULFedges = 74 ;
        TGULF = UNLIMITED ; // (4 currently)
        XWET = 295000 ;
variables:
        float XGULF(XGULF) ;
                XGULF:point_spacing = "even" ;
        float YGULF(YGULF) ;
                YGULF:point_spacing = "even" ;
        float ZGULF(ZGULF) ;
                ZGULF:positive = "down" ;
                ZGULF:point_spacing = "uneven" ;
                ZGULF:edges = "ZGULFedges" ;
        float ZGULFedges(ZGULFedges) ;
        float TGULF(TGULF) ;
                TGULF:units = "HOURS since 1901-01-15 00:00:00" ;
                TGULF:time_origin = "15-JAN-1901" ;
        int XWET(XWET) ;
                IWET:compress = "XGULF YGULF ZGULF" ;
                IWET:point_spacing = "even" ;
        float TEMP(TGULF, XWET) ;
                TEMP:fill_value = 9.96921e+36f ;
                TEMP:missing_value = 9.96921e+36f ;
                TEMP:long_name = "TEMP" ;
                TEMP:units = "celcius" ;
        float WET(XWET) ;
                WET:fill_value = 9.96921e+36f ;
                WET:missing_value = 9.96921e+36f ;
                WET:long_name = "Indices points mouilles" ;


I then modified my external function so as to supply to it the WET
variable which contains the wet-cell indices:

set memory/size=500
use COMP_TEST.ROM_TEMP.nc
let a = x[gx=xgulf]*0 + y[gy=ygulf]*0 + z[gz=zgulf]*0 + t[gt=tgulf]*0
let b = nc_rar(a,temp,wet)
shade b[l=1,k=1]


The problem is that now inside my nc_rar() function all variables (temp
and wet) have zero values everywhere.  Everything else inside my
function appears to be right: arg_lo_ss, arg_hi_ss, etc. so I obviously
have an addressing problem or something similar. 

So here is my final question:

Has anybody ever encountered this behaviour with external functions and
managed to solve it?


I realize that Ferret does not support the CF Convention for NetCDF
files but the disk space savings resulting from compression by gathering
are very interesting.  So, I intend to investigate this avenue further
for our own needs.  Any tip or advice to help me solve this problem
would be greatly appreciated.




James



-- 
James Caveen
Analyste-programmeur
Service des technologies de l'information
Université du Québec à Rimouski
310 des Ursulines
Rimouski (Québec) G5L 3A1
Tel: (418) 723-1986 poste 1295
Fax: (418) 724-1842