NCL Website header
NCL Home > Documentation > Manuals > Reference

Extending the NCL function and procedure set

If you have a Fortran or C function or procedure that you'd like to call from NCL, you can do it by "wrapping" this function a number of different ways.

The wrapper file is a C program that describes all the arguments and passes them back and forth between the function/procedure you want to call, and the NCL script that is calling it.

There are multiple ways you can generate this C wrapper file and create the resultant shared object:

  • Using the WRAPIT script - this method only works with Fortran and Fortran 90 routines. It generates the wrapper file, compiles all the necessary code, and creates the shared object file. This is the easiest method to use, but doesn't allow you to have control over the interface.

  • Using the wrapit77 application - this method is similar to using WRAPIT, but you need to compile the code yourself. This method allows you to modify the wrapper interface if you want. It requires knowing how to generate shared objects on your system.

  • Writing your own wrapper as a built-in function
  • Writing your own wrapper as a shared object - these two methods work for wrapping both C and Fortran routines. They are a bit more complicated and require a fair amount of C and Fortran knowledge. Also, if the routine you're wrapping is a Fortran routine, they require that you know how to call Fortran from C. These methods allow you to have full control over the interface to the routine you are wrapping. The difference between the two is that one allows you to create a new NCL executable so that you don't need to load anything to access your new function, and the other method requires that you load a shared object.
Be sure to visit the additional examples section and the special considerations section for trouble shooting.

Using the WRAPIT script

This method only works for Fortran 77 or 90 code.

Step 1
Fortran 77 - Write code and add special wrapper text
Fortran 90 - Write code and write separate special wrapper text
Commercial library routine - Write separate special wrapper text
Step 2
Fortran 77 - Run WRAPIT to compile the external Fortran code(s)
Fortran 90 - Run WRAPIT to compile the external Fortran code(s) and the separate wrapper text
Commercial library routine - Run WRAPIT to compile the separate wrapper text and link in the commercial library
Step 3
Call the shared object from an NCL script.

Step 1

Fortran 77 - Write code and add special wrapper text:

If you have a Fortran 77 routine, use the special interface delimiters "C NCLFORTSTART" and "C NCLEND" to bracket the argument declarations right within the Fortran code.

For example, assume you have a Fortran code called "ex01.f":

C NCLFORTSTART
      subroutine cquad (a, b, c, nq, x, quad)
      real x(nq),quad(nq)
C NCLEND
C
C  Calculate quadratic polynomial values.
C
      do 10 i=1,nq
        quad(i) = a*x(i)**2 + b*x(i) + c
   10 continue

      return
      end

C NCLFORTSTART
      function arcln (numpnt, pointx, pointy)
      dimension pointx(numpnt),pointy(numpnt)
C NCLEND

C
C  Calculate arc lengths.
C
      if (numpnt .lt. 2) then
        print *, 'arcln: number of points must be at least 2'
        stop
      endif
      arcln = 0.
      do 10 i=2,numpnt
        pdist = sqrt((pointx(i)-pointx(i-1))**2 + 
     +                        (pointy(i)-pointy(i-1))**2)
        arcln = arcln + pdist
   10 continue

      return
      end

Only the argument variables should be included between the delimiters. If a particular variable is not typed, then the usual Fortran rules will apply for typing. Additionally, only the Fortran subroutine(s) actually called from NCL require the interface delimiters. The rest of the Fortran code may include other subroutines without delimiters.

Fortran 90 - Write code and write special wrapper text:

If you have a Fortran 90 file, then you need to create a separate "stub" file that contains the "C NCLFORTSTART" and "C NCLEND" delimiters. For example, assume the above "cquad" routine was in a Fortran 90 file called "cquad.f90":

subroutine cquad(a,b,c,nq,x,quad)
  implicit none
  integer, intent(in)  ::nq
  real,    intent(in)  ::a,b,c,x(nq)
  real,    intent(out) ::quad(nq)
  integer              ::i

  quad = a*x**2+b*x+c
  return
end subroutine cquad

Create a separate file, called something like "cquad90.stub" that contains nothing more than the following lines:

C NCLFORTSTART
      subroutine cquad(a,b,c,nq,x,quad)
      real a,b,c
      integer nq
      dimension x(nq),quad(nq)
C NCLEND

Commercial library routine - Write special wrapper text:

If you have a commercial library routine you want to call from NCL, use the same method as with the Fortran 90 routine to create a separate "stub" file.

For example, assume you want to call the IMSL routine "rline" to fit a line to a set of data points via least-squares. The "rline" arguments are rline(nobs,x,y,b0,b1,stat) where "nobs" is the number of observations, "x" and "y" are the data vectors, "b0" and "b1" are the intercept and slope, and "stat" is a vector of length 12 containing assorted statistics. The Fortran stub file you create (call it "rline.stub") would look like:

C NCLFORTSTART
      subroutine rline (n,x,y,b0,b1,stat)
      integer n                                ! explicit typing NOT required
      real    x(n), y(n), b0, b1, stat(12)
C NCLEND

Step 2

Fortran 77 - Run WRAPIT to compile the external code(s):

Run WRAPIT on your code to generate a shared object. It will have the same name as your Fortran file, except with ".so" appended instead of ".f" or ".f90. For the Fortran 77 example:

   WRAPIT ex01.f

This should create a file called "ex01.so".

Fortran 90 - Run WRAPIT to compile the external code(s) and the separate wrapper text:

Using the "cquad90.stub" file you created in the previous step, and assuming "cquad.f90" is the Fortran 90 file that contains the "cquad" subroutine, then you would run WRAPIT as follows:

   WRAPIT cquad90.stub cquad.f90

This should create a file called "cquad90.so".

Commercial library routine - Run WRAPIT to compile the separate wrapper text and link in the commercial library:

To build a shared object that includes a call to one or more commercial library routines, you must include the list of commercial libraries on WRAPIT command line so the shared object will be properly linked:

   WRAPIT -l imsl_mp rline.stub

This should create a file called "rline.so".

Note for all types of routines:

If WRAPIT does not appear to work, then you can modify the WRAPIT script directly ($NCARG_ROOT/bin/WRAPIT) and change the paths to your appropriate local compilers. WRAPIT does not contain paths to any libraries, so they may have to be added depending upon your system.

Step 3 - Call the shared object from an NCL script

Before you can call external subroutines from your NCL script, you need to load the shared object you created in Step 2. There are two ways to load a shared object:

  1. Use the external statement
  2. Use the NCL_DEF_LIB_DIR environment variable
Loading the shared object using external:

At the top of your NCL script and before the "begin" statement, add an "external" statement to load the shared object, followed by a name you want to give the shared object, followed by the path to the shared object in double quotes.

The name of the shared object is arbitrary, but by convention, is capitalized. If the shared object is in a different directory, be sure to include the path to it (you can use relative or absolute paths). If the shared object is in the same directory as your script, be sure to put a "./" in front of it.

Now, to call any one of your functions or procedures from your NCL script, precede it with two colons ("::") and the name you gave the shared object. The example below is from the Fortran 77 example, and assumes the shared object "ex01.so" is in the same directory as your script:

external EX01 "./ex01.so"
 
begin
;
; Calculate three values of a quadratic equation
;
   nump = 3
   x    = (/ -1., 0.0, 1.0 /)
   qval = new(nump,float)              
   EX01::cquad(-1., 2., 3., nump, x, qval) ; Call the NCL version of
                                           ; your Fortran subroutine.
   print("Polynomial value = " + qval)     ; Should be (/0,3,4/)
 
;
; Calculate an arc length.
;
   xc = (/ 0., 1., 2. /)
   yc = (/ 0., 1., 0. /)
   arclen = EX01::arcln(nump,xc,yc)     ; Call the NCL version of
                                        ; your Fortran function.
   print("Arc length = " + arclen)      ;  should be 2.82843
end

Loading the shared object using an environment variable:

Create a directory on the machine you wish to run NCL, and put your shared object(s) in this directory. Set the environment variable NCL_DEF_LIB_DIR to this path. For example:

setenv NCL_DEF_LIB_DIR /home/haley/shared_objects

NCL will recognize the path given by the NCL_DEF_LIB_DIR environment variable as another place to look for shared objects. You can include multiple directory paths by separating them with colons:

setenv NCL_DEF_LIB_DIR /home/shea/shared_objects/:/home/haley/shared_objects

Now you can call the shared object just like a built-in NCL function:

begin
;
; Calculate three values of a quadratic equation
;
   nump = 3
   x    = (/ -1., 0.0, 1.0 /)
   qval = new(nump,float)              
   cquad(-1., 2., 3., nump, x, qval)       ; Call the NCL version of
                                           ; your Fortran subroutine.
   print("Polynomial value = " + qval)     ; Should be (/0,3,4/)
 
;
; Calculate an arc length.
;
   xc = (/ 0., 1., 2. /)
   yc = (/ 0., 1., 0. /)
   arclen = arcln(nump,xc,yc)           ; Call the NCL version of
                                        ; your Fortran function.
   print("Arc length = " + arclen)      ;  should be 2.82843
end

Be sure to see the section on special considerations for trouble shooting.

WRAPIT options

There are some command line options you can use with WRAPIT that are useful for debugging. These options must appear on the WRAPIT command line before the Fortran file names:
-d
Turns on array bounds, turns off optimization, displays some debug information, and prevents file clean up.

-h or -help
Gives you LOTS of information about WRAPIT, including information about other options.

-in
Use the Intel compiler.

-l <libname>
Passes a library name to the linker.

-L <libpath>
Passes a directory path to the linker.

-lf
Use the Lahey compiler.

-n <so name>
Assigns a name to the created shared object.

-pg
Use the Portland compiler.

-q32
Specifies 32-bit IRIX or AIX bit precision (the default is to use 64-bit precision, if available).

-r8
Promotes Fortran floats of real*4 to real*8 (if available).

Using the wrapit77 application

This section is for more advanced users who want to have more control over the interface to the routine they are calling.

Basically, the steps below are some of the ones that the WRAPIT script does for you automatically (see the previous section), but are broken down into individual steps here so you can intervene and make changes where you want.

Generating your own wrapper involves using a utility application called wrapit77. This utility facilitates generating a C source code wrapper for Fortran 77 functions and subroutines. Below is a step by step description of how to use wrapit77.

Starting with a Fortran procedure and ending with a functionally equivalent NCL procedure is a five-step process:

Step 1
Create code blocks that describe the Fortran procedure and its arguments. These code blocks are similar to Fortran 90 interface blocks (although the blocks created here must use Fortran 77 syntax), and C function prototypes.
Step 2
Create an NCL wrapper function by running the program wrapit77 that is supplied with the NCAR Graphics Software (version 4.1 or later). wrapit77 uses the code blocks in Step 1 as input.
Step 3
Create object modules (".o" files) for the NCL wrapper function generated in Step 2 and for the original Fortran source.
Step 4
Create a dynamic shared object from the object files created in Step 3 using the UNIX ld command. A dynamic shared object is an object module whose entries are loaded on an as-needed basis.
Step 5
Inform NCL where to locate the dynamic shared object created in Step 4.
Consider each of the above steps in detail:

Step 1 - describe your procedure

A wrapit interface block is a sequence of Fortran 77 statements that specifies a procedure and its arguments. These statements are bracketed by Fortran comment lines:

C NCLFORTSTART 

at the start and

C NCLEND 

at the end. These comment lines are known as wrapit interface block delimiters. Any blanks in the above two lines are ignored and the lines are case-insensitive.

The following is an example of a wrapit interface block for a Fortran function that finds the arc length of a polyline:

C NCLFORTSTART 
      FUNCTION ARCLN(NUMPNT, POINTX, POINTY)
      DIMENSION POINTX(NUMPNT), POINTY(NUMPNT)
C NCLEND

The statements between wrapit interface block delimiters should contain only a procedure statement and non-executable declarations describing its arguments. Any other statements (except comment statements) in a wrapit interface block will cause wrapit77 to exit with an error. Unless specifically typed, all variables are typed according to Fortran default typing.

Fortran 90 codes can use wrapit77 to generate an NCL wrapper function if you can provide a wrapit interface block for it, i.e. provided that you can specify the Fortran 90 procedure interface with Fortran 77 syntax. In some cases this will not be possible, for example if you are using pointers or structures as arguments, if the Fortran 90 procedure can have missing arguments, if the procedure is recursive, or if you are using keyword arguments. Many times, however, it is easy to convert Fortran 90 syntax to Fortran 77. For example, the Fortran 90 statements:

      CHARACTER (LEN=*)::NAME
      REAL,INTENT(IN)::NUMBER

could be specified as

      CHARACTER*(*) NAME
      REAL NUMBER

in the wrapit interface block.

Note: wrapit77 scans its input file for wrapit interface blocks, so if you can easily arrange your Fortran codes so that the wrapit interface blocks are embedded in them, it may be more convenient for you to use your original Fortran source file as input to wrapit77, rather than use a separate file containing only wrapit interface blocks. See example 2 in this section for details on this.

Step 2 - create an NCL wrapper function by running wrapit77

wrapit77 scans its input for wrapit interface blocks and produces an NCL wrapper function as output. Any source code in the input to wrapit77 that lies outside of a wrapit interface block is ignored, so the input to wrapit77 can simply be a sequence of wrapit interface blocks, or it can be complete Fortran source with embedded wrapit interface blocks.

For example, if your input wrapit interface block is in file wrapper_input, then executing

wrapit77 < wrapper_input >! wrapper_W.c

will create the NCL wrapper function wrapper_W.c. See the examples in this section for sample uses of wrapit77. You can name your NCL wrapper function anything you want, but the usual convention is to append the "_W" followed by ".c" (since the wrapper function is C code).

At this point, you can modify the "wrapper_W.c" file and change the interface however you see fit. For example, if your Fortran routine requires that several work arrays be created and passed to it, you can do this from the C wrapper routine, hence removing the need to create the work arrays inside the NCL script.

Step 3 - create object modules

Object modules (".o" files) need to be created for your original Fortran source and for the NCL wrapper function created in Step 2. The NCL wrapper function is C source and contains some special include files. You should use the nhlcc command to create the ".o" file for the NCL wrapper function, since it will correctly locate the include files. For example, if you have named your NCL wrapper function wrapper_W.c, then execute:

nhlcc -c wrapper_W.c

to create wrapper_W.o. For consistency with the use of nhlcc, you may want to use the command nhlf77 to create the ".o" file for your Fortran source, but any way is fine.

Step 4 - create a dynamic shared object module

Dynamic shared objects are object modules whose entries are loaded on an as-needed basis. The UNIX command ld is used to create the dynamic shared objects (".so" files) used with NCL. The object modules (".o" files) you created in Step 3 will be used to create your dynamic shared object. For each of the supported systems (except for Crays, where wrapit77 is not currently supported), the requirements are different. Here is a summary of what you need to do on each system (where the NCL wrapper function object code is named fcode_W.o and the object for the original Fortran is named fcode.o) to create a dynamic shared object named fcode.so:

  • System running Linux RedHat, or DEC Alpha running Digital UNIX:
      nhlcc -c fcode_W.c
      nhlf77 -c fcode.f
      ld -shared -o fcode.so fcode_W.o fcode.o
    

  • HP 700 series running HPUX:
      nhlcc +z -c fcode_W.c
      nhlf77 +z -c fcode.f
      ld -b -o fcode.so fcode_W.o fcode.o 
    

  • SGI running IRIX64 (-64):
      nhlcc -c fcode_W.c
      nhlf77 -c fcode.f
      ld -64 -shared -o fcode.so fcode_W.o fcode.o
    

  • SGI running IRIX64 (-n32):
      nhlcc -c fcode_W.c
      nhlf77 -c fcode.f
      ld -n32 -shared -o fcode.so fcode_W.o fcode.o 
    

  • Sun running Solaris:
      nhlcc -c fcode_W.c
      nhlf77 -c fcode.f
      ld -o fcode.so fcode_W.o fcode.o -G
    

If the above doesn't cover your system, or it doesn't seem to work when you try it, then look at the $NCARG_ROOT/bin/WRAPIT script. This script was designed to work on many different systems with different compilers, and is fairly well comment so that you can follow it.

Step 5 - inform NCL where to locate your dynamic shared object

NCL gets its information about externally defined procedures from dynamic shared objects, as created in Step 4.

Once you have created your dynamic shared object containing the NCL wrapper function and the object for the original code, there are three ways to tell NCL where to look for the dynamic shared object:

  1. Use an NCL external statement.
  2. Use the system environment variable LD_LIBRARY_PATH.
  3. Use the NCL environment variable NCL_DEF_LIB_DIR.

Accessing dynamic shared objects using an external statement

In this case, you insert a statement of the form:
external SO_NAME "path_name"

at the beginning of your NCL script. In this statement SO_NAME indicates the name you want the dynamic shared object to be known by in the NCL program and path_name is a relative or full pathname of the dynamic shared object file you created. The final argument to the external statement is an NCL string variable, so the quotes are required. The case where path_name is simply the name of the dynamic shared object itself (and not a pathname) is special and is treated in the "Accessing dynamic shared objects using LD_LIBRARY_PATH" section.

Once the dynamic shared object has been accessed and named, you can call your new NCL function or procedure with an NCL statement of the form:

SO_NAME::fortran_name(arguments)

where SO_NAME is the name that you gave the dynamic shared object in your NCL external statement and fortran_name is the name of the Fortran function or subroutine that you specified in your wrapit interface block. fortran_name must be entered as either all uppercase, or all lowercase, regardless of what case the original Fortran procedure name was. Using an NCL name that is all uppercase or lowercase will produce the same results. For example, if your original Fortran procedure was named "CalcX", then you can use either "CALCX" or "calcx" to invoke that function from NCL, but you cannot use "CalcX" or any other mixed-case form of the name. Most NCL programmers tend to use lowercase for procedure names. See the examples in this section for complete examples of using dynamic shared objects and calling your new NCL procedures.

Accessing dynamic shared objects using LD_LIBRARY_PATH

In this case, you insert a statement of the form:
external SO_NAME "shared_object_name"

at the beginning of your NCL script. In this statement, shared_object_name is the name of your dynamic shared object file. NCL will look in the directories listed in the environment variable LD_LIBRARY_PATH for the specified dynamic shared object. If it does not find it there, an error will be reported. LD_LIBRARY_PATH is a system environment variable and is a colon-separated list of directory pathnames. In attempting to set this environment variable, make sure that it is not already set for use by other applications; if it is, just add the pathname of the dynamic shared object to it.

Once NCL has located the dynamic shared object, you can invoke the function or procedure in the same way that was discussed above under "Accessing dynamic shared objects using an external statement."

Accessing dynamic shared objects using NCL_DEF_LIB_DIR

If the NCL environment variable NCL_DEF_LIB_DIR is defined and is a valid directory pathname, NCL will look in that directory first for dynamic shared objects (you can have as many dynamic shared objects in this directory as you want). Any entry in any dynamic shared object in the directory specified by NCL_DEF_LIB_DIR will be loaded automatically as an NCL built-in. In particular this means that, in order to invoke your new function or procedure, you should not prepend SO_NAME:: to the function or procedure name as was required in either of the two other methods for accessing dynamic shared objects discussed above. NCL will try to load everything in the directory specified by NCL_DEF_LIB_DIR, so you will get warning messages about files in the directory that are not shared objects.

This approach can be more convenient than the other two approaches, but it lacks the advantage of having the SO_NAME:: prepended to externally defined procedures which makes their externally defined status quickly recognized from the syntax.

The environment variables NCL_DEF_LIB_DIR and LD_LIBRARY_PATH are independent. The dynamic shared objects in the directory specified by NCL_DEF_LIB_DIR are loaded before any dynamic shared object specified on an external statement, whether that dynamic shared object is specified by way of a relative or direct pathname or by way of a directory specified using the LD_LIBRARY_PATH environment variable.

Be sure to see the section on special considerations for trouble shooting.

Writing your own wrapper as a built-in function

NCL has been designed to allow users the ability to extend the built-in default function set. This functionality requires a fair amount of C and Fortran knowledge. In addition, if you intend to call a Fortran function, some specific knowledge of how to call Fortran from C is needed. With release 4.1 extending NCL to call FORTRAN subroutines and functions has been greatly simplified. At release 4.1 the external statement was added to support loading of specially created shared libraries. Also a utility called wrapit77 was added to generate all the C code necessary to build an NCL shared library.

If you want to build your own custom version of NCL which includes your functions and procedures, the following steps should be followed. However, if you want to enable the installed version of NCL to call your functions and procedures using the external statement, then follow the instructions for writing a wrapper; if you're calling a C function, and then follow the directions for building a shared object. If you're calling FORTRAN then follow all the steps for using wrapit77

The general approach for extending the NCL function set is as follows:

  • Write a function or procedure in C or Fortran.

  • Write a wrapper function that calls the internal NCL function NclGetArgValue to retrieve parameters from NCL and then calls the actual function. The wrapper function can also call the internal function NclReturnValue to pass a return value back to NCL if necessary.

  • Write a C function called NclAddUserFuncs and use SetArgTemplate with either NclRegisterFunc or NclRegisterProc to assign an NCL function name to the user function.

  • Compile and link files containing the actual function wrapper, the wrapper function, and NclAddUserFuncs with the library libncl.a to create a new executable of NCL.
The NCL binaries installed on your system should contain a library called "libncl.a". To determine if this library has been installed on your system, run NCL and enter the following command:
    ncl 0> system("ls " + ncargpath("lib"))
    
If libncl.a is not one of the libraries listed by the above command, contact your site administrator.

This section uses an example to demonstrate how to extend the function set. In this example we will write a function that evaluates a polynomial given a value and arrays of polynomial coefficients. The intended NCL function interface is:

function poly(
        a[*]:float,
        b[*]:float,
        c[*]:float,
        x[1]:float
)

Writing the function in C or Fortran

There are only a few restrictions necessary for C and Fortran functions. First and most important, the C function, Fortran function, or Fortran subroutine name must not conflict with any function names in the NCL library. Use the UNIX function "nm" to accomplish this. Second, for C, none of the data passed in from the wrapper can be freed. Finally, for Fortran, if I/O is being done you must make sure that the unit number being used is not already in use. It is also important to follow good programming practice by freeing any memory allocated by the function unless it is being passed back as a return value, and closing all files opened by the function. Also, only the internal functions specified on this page should be used.

The following is the C source for the function poly in a file called poly.c :

void poly_C(float *outval, float *a, float *b, float *c, float x, int n) {
        int i;

        for(i = 0; i < n ; i++) {
                outval[i] = (a[i] * x * x) + (b[i] * x) + c[i];
        }
}

Writing the wrapper function

The wrapper function must be written in C and use the function NclGetArgValue to retrieve the values of the parameters from NCL. The interface for the NclGetArgValue provides information as to the number of dimensions, dimension sizes, the missing values if any, and the type. The prototype for NclGetArgValues is located in "NclBuiltInSupport.h". The interface for NclGetArgValue is included below.

Important note: in NCL versions 6.0.0 and later, the variables that hold the dimension sizes (they usually have the word "dimsizes" in the name in the examples below) must be of type "ng_size_t" and not "int". This is because NCL was upgraded in V6.0.0 to allow dimension sizes to be integers or longs. On 64-bit systems, this allows you to have dimension sizes greater than 2 gigabytes.

If you are using version NCL 5.2.1 or earlier, then use "int" any place where you see "ng_size_t" below.

    void *NclGetArgValue(
            int arg_num,                        /* IN:  argument number to retrieve */
            int n_args,                         /* IN:  total number of arguments */
    	int* n_dims,                        /* OUT: number of dimensions of parameter */
            ng_size_t* dimsizes,                /* OUT: dimension sizes of parameter */
            NclScalar* missing,                 /* OUT: missing value if any */
            int *has_missing,                   /* OUT: whether or not a missing value is 
                                                    present in the parameter */
            NclBasicDataTypes *type,            /* OUT: data type of parameter */
            int access_type                     /* IN:  Either 0 for don't care, 2 for a read or 1 for a write*/
    );
    

The input argument arg_num is the number of the argument requested. The arguments are numbered from left to right (as they appear in NCL) starting at 0 through n-1 where n is the total number of arguments to the function. The input argument n_args is the total number of arguments for the NCL function. The final input parameter is access_type. Access_type tells NCL whether or not the parameter being accessed will be modified or not. This is necessary to support NCL's pass-by-value parameter passing scheme.

The actual data associated with the parameter is returned as a void pointer by the NclGetArgValue function which must be cast to the appropriate type. The parameters n_dims and dimsizes represent the rank of the parameter. The type parameter is set to one of the following basic data types. NclBasicDataTypes is defined in "NclDataDefs.h".

    typedef enum  {
    	NCL_none =      0,
    	NCL_char =      010,
    	NCL_byte =      011,
    	NCL_short =     020,
    	NCL_ushort =    021,
    	NCL_int =       040,
    	NCL_uint =      041,
    	NCL_float =     042,
    	NCL_long =      044,
    	NCL_ulong =     045,
    	NCL_int64 =     0100,
    	NCL_uint64 =    0101,
    	NCL_double =    0102,
    	NCL_string =    0200,
    	NCL_numeric =   01000,
    	NCL_enumeric =  02000,
    	NCL_snumeric =  04000,
    	NCL_logical =   010000,
    	NCL_obj =       020000,
    	NCL_list =      040000
    } NclBasicDataTypes;
    

The contents of the has_missing parameter is set to 1 if the parameter possibly contains missing values. There are some situations where it is possible to have has_missing set and not actually have a missing value in the data. This occurs when a subsection of a variable, which has missing values, is subscripted and the portion subscripted does not actually contain missing values. However, if the has_missing parameter is set to 0, then there are no missing values. The value of the missing value is returned in the union NclScalar. NclScalar is defined in "NclDataDefs.h". It is a union of all the primitive types in NCL.

    typedef union _NclScalar {
            double             doubleval;
            unsigned long long uint64val;
            long long          int64val;
            unsigned short     ushortval;
            unsigned int       uintval;
            unsigned long      ulongval;
            float              floatval;
            int                intval;
            long               longval;
            short              shortval;
            unsigned char      charval;
            string             stringval;
            byte               byteval;
            logical            logicalval;
            obj                objval;
    }NclScalar;
    

Finally, return values are passed back to the NCL environment via the NclReturnValue. NclReturnValue provides an interface for passing back data. It currently does not support passing back coordinate arrays or attributes, with the exception of missing values.

The prototype for the NclReturnValues function is:

    NhlErrorTypes NclReturnValue(
    void *value,                       /* Pointer to data to be returned */
    int n_dims,                        /* Number of dimensions of the data */
    ng_size_t* dimsizes,               /* Dimension sizes of the data */
    NclScalar* missing,                /* If non-NULL is the missing value for the data 
                                          otherwise NULL if no missing values in data */
    NclBasicDataTypes type,            /* Ncl basic type identifier */
    int copy_data                      /* True if data should be copied before returned */
    );
    
One major goal of designing the wrapper function should be to establish a scheme for handling missing values if the function you are adding does not already handle them. In the following example, missing values are filtered out and propagated through to the output.

The source for the wrapper function follows:

    #include <stdio.h>
    /*
    * The following are the required NCAR Graphics include files.
    * They should be located in ${NCARG_ROOT}/include
    */
    #include <ncarg/hlu/hlu.h>
    #include <ncarg/hlu/NresDB.h>
    #include <ncarg/ncl/defs.h>
    #include <ncarg/ncl/NclDataDefs.h>
    #include <ncarg/ncl/NclBuiltInSupport.h>
    
    /*
    * Declare C function version
    */
    extern void poly_C(float *outval, float *a, float *b, float *c, float x, int n);
     
    NhlErrorTypes poly_W( void )
    {
    
            float *a;
            int n_dims_a,has_missing_a;
            int dimsizes_a[NCL_MAX_DIMENSIONS];
            NclScalar missing_a;
            float *b;
            int n_dims_b,has_missing_b;
            ng_size_t dimsizes_b[NCL_MAX_DIMENSIONS];
            NclScalar missing_b;
            float *c;
            int n_dims_c,has_missing_c;
            int dimsizes_c[NCL_MAX_DIMENSIONS];
            NclScalar missing_c;
            float *x;
            NclScalar missing_x;
            float *output_val;
            NclScalar tmp_missing;
            int has_missing_x;
            ng_size_t i;
    
    /*
    * Retrieve parameters
    */
    
    /*
    * Note any of the pointer parameters can be set to NULL, which
    * implies you don't care about the its value. In this example
    * the type parameter is set to NULL because the function
    * is later registered to only accept floating point numbers.
    * Also the parameter x is later registered to be a single
    * dimension floating point value so dimension sizes and the
    * number of dimensions are not needed.
    */
            a = (float*)NclGetArgValue(
                    0,
                    4,
                    &n_dims_a,
                    dimsizes_a,
                    &missing_a,
                    &has_missing_a,
                    NULL,
                    2
            );
            b = (float*)NclGetArgValue(
                    1,
                    4,
                    &n_dims_b,
                    dimsizes_b,
                    &missing_b,
                    &has_missing_b,
                    NULL,
                    2
            );
            c = (float*)NclGetArgValue(
                    2,
                    4,
                    &n_dims_c,
                    dimsizes_c,
                    &missing_c,
                    &has_missing_c,
                    NULL,
                    2
            );
            x = (float*)NclGetArgValue(
                    3,
                    4,
                    NULL,
                    NULL,
                    &missing_x,
                    &has_missing_x,
                    NULL,
                    2
            );
    
    /*
    * This is the only dimension size check needed since the function
    * is registered to only accept single dimension parameters.
    */
            if((dimsizes_a[0] == dimsizes_b[0])&&(dimsizes_b[0] == dimsizes_c[0])) {
    /*
    * The case where x is a missing value should be handled. One option
    * is to print an error message and return NhlFATAL. However a better
    * option is to return an array of size n_dims_a filled with missing
    * values.
    */
                    if((has_missing_x) && (missing_x.floatval == *x)) {
                            output_val = (float*)malloc(sizeof(float)*dimsizes_a[0]);
                            for(i = 0; i < n_dims_a; i++) {
                                    output_val[i] = *x;
                            }
                            return(NclReturnValue(
                                    (void*)output_val,
                                    n_dims_a,
                                    dimsizes_a,
                                    &missing_x,
                                    NCL_float,
                                    0
                            ));
                    }
    /*
    * The following section allocates the output memory and calls the
    * poly function
    */
                    output_val = (float*)malloc(sizeof(float)*dimsizes_a[0]);
    
    /*
     * If dimsizes_a[0] is greater than 2 gigabytes, you will need to
     * make this argument a "long" on a 64-bit machine.
     */
                    poly_C(output_val,a,b,c,*x,(int)dimsizes_a[0]);
    /*
    * A scheme for handling missing values in the input. NCL's algebra
    * handles missing value by propagating the missing value of the left
    * most term in an expression to the result whenever any term contains
    * a missing value.
    */
                    if(has_missing_a || has_missing_b || has_missing_c) {
    /*
    * Chooses suitable missing value by checking from left to right for
    * missing values.
    */
                            if(has_missing_a) {
                                    tmp_missing.floatval == missing_a.floatval;
                            } else if (has_missing_b) {
                                    tmp_missing.floatval == missing_b.floatval;
                            } else {
                                    tmp_missing.floatval == missing_c.floatval;
                            }
    /*
    * Checks for missing values in the input and changes output value to missing
    * value if any of the terms contain missing values
    */
                            for(i = 0; i < dimsizes_a[0]; i++) {
                                    if((has_missing_a)&&(a[i] == missing_a.floatval)) {
                                            output_val[i] = tmp_missing.floatval;
                                    }
                                    if((has_missing_b)&&(b[i] == missing_b.floatval)) {
                                            output_val[i] = tmp_missing.floatval;
                                    }
                                    if((has_missing_c)&&(c[i] == missing_c.floatval)) {
                                            output_val[i] = tmp_missing.floatval;
                                    }
                            }
    /*
    * Returns value to NCL
    */
                            return(NclReturnValue(
                                    (void*)output_val,
                                    n_dims_a,
                                    dimsizes_a,
                                    &tmp_missing,
                                    NCL_float,
                                    0
                            ));
                    } else {
    /*
    * Returns value to NCL
    */
                            return(NclReturnValue(
                                    (void*)output_val,
                                    n_dims_a,
                                    dimsizes_a,
                                    NULL,
                                    NCL_float,
                                    0
                            ));
                    }
            } else {
                    NhlPError(NhlFATAL,NhlEUNKNOWN,"poly: the dimension sizes of parameters a, b and c must be identical");
                    return(NhlFATAL);
            }
    
    }
    

Writing NclAddUserFuncs to register the function with NCL

Note:If you intend to create a shared library for use with either NCL_DEF_LIBRARY_DIR or the external statement follow these instructions but name the function Init rather than NclAddUserFuncs. Also if you intend to call a FORTRAN function or subroutine the utility wrapit77 can be used in place of this step in the wrapper creation.

The next step in extending the NCL function set is to write a function called NclAddUserFuncs. This function takes no arguments and does not return anything. Its purpose is to allow the user to call NclRegisterFunc and NclRegisterProc to tell NCL the NCL string name of the function, what the function pointer to call is, how many parameters the function takes, and what the type and dimensionality of the parameters are.

There are three C functions that must be called to register a C wrapper function as as an NCL built-in function. For each function being registered in NclAddUserFuncs, the function NewArgs must be called to initialize a private array that will hold the descriptions of each parameter. NewArgs takes a single integer argument that is the number of arguments for the function to be registered. It returns a void pointer that must be used in subsequent calls to SetArgTemplate and NclRegisterFunc.

The next function called is SetArgTemplate. It must be called for each parameter to the function. SetArgTemplate takes the pointer returned by NewArgs, the current parameter number, the parameter's type, number of dimensions, and dimension sizes as input. SetArgTemplate does not return a value; it simply modifies the private array returned by NewArgs.

Once each parameter has been defined by SetArgTemplate, then either the function NclRegisterFunc or NclRegisterProc can be called depending on whether the C function being added is meant to be an NCL function or procedure. The prototypes for these four functions are located in the include file NclBuiltIns.h, which should be located in ${NCARG_ROOT}/include/ncarg/ncl:

    void *NewArgs(
    int  n                         /* Total number of parameters to function being registered */
    );
    
    void SetArgTemplate(
            void*   args,           /* Array returned by NewArgs */
            int     arg_num,        /* Argument number being defined */
            char*   type_name,      /* string representing type, set to NclANY if any 
                                       type is accepted */
            int     n_dims,         /* number of dimensions of argument, set to 
                                       NclANY if any number of dimensions supported */
            ng_size_t*    dimsizes    /* dimension size for each dimension, set to NclANY 
                               if all dimension sizes are supported */
    );
    
    extern void NclRegisterFunc(
            NclPubBuiltInProcWrapper thefuncptr,      /* The C wrapper function */
            void* args,                               /* The array returned by NewArgs and initialized
                                                         by successive calls to SetArgTemplate */
            char* fname,                              /* String name of function to be used by NCL */
            int n_args                                /* Total number of arguments */
    );
    

The following is the source for the function NclAddUserFuncs. More than one user function can be registered from within this function.

    
    #include <stdio.h>
    #include <ncarg/hlu/hlu.h>
    #include <ncarg/hlu/NresDB.h>
    #include <ncarg/ncl/defs.h>
    #include <ncarg/ncl/NclBuiltIns.h>
    
    /*
    * Declare wrapper function
    */
    extern NhlErrorTypes poly_W(
    #if NhlNeedProto
    void
    #endif
    );
    
    void NclAddUserFuncs
    #if NhlNeedProto
    (void)
    #else
    ()
    #endif
    {
            void *args;
            ng_size_t dimsizes[NCL_MAX_DIMENSIONS];
            int nargs = 0;
    /*
    * Create private argument array
    */
            args = NewArgs(4);
    /*
    * Configure first three parameters identically as single dimension float
    * arrays of any size
    */
            SetArgTemplate(args,nargs,"float",1,NclANY);nargs++;
            SetArgTemplate(args,nargs,"float",1,NclANY);nargs++;
            SetArgTemplate(args,nargs,"float",1,NclANY);nargs++;
    /*
    * Configure fourth parameter as single scalar floating point value
    */
            dimsizes[0] = 1;
            SetArgTemplate(args,nargs,"float",1,dimsizes);nargs++;
    /*
    * Register wrapper function pointer  and argument templates
    */
            NclRegisterFunc(poly_W,args,"poly",nargs);
    
            return;
    }
    
    /* 
    * NOTE: the following function stubs must be included with the 
    * source for function NclAddUserFuncs. Future releases of this
    * documentation will describe their uses. For now though NCL
    * WILL NOT compile WITHOUT them.
    */
    
    void NclAddUserFileFormats
    #if NhlNeedProto
    (void)
    #else
    ()
    #endif
    {
            return;
    }
    
    void NclAddUserHLUObjs
    #if        NhlNeedProto
    (void)
    #else
    ()
    #endif
    {
            return;
    }
    

Compiling and linking files

The directions for compiling a shared object can be found in the section titled Using wrapit77 to generate wrappers. This is the recommended method of extending NCL since a new version of NCL does not need to be created.

The final step is to compile and link the source files created with libncl.a and the other NCAR Graphics libraries. The simplest way to compile these functions into an executable is to use the nhlcc script from the UNIX prompt. nhlcc should be located in ${NCARG_ROOT}/bin. Nhlcc is a script that will link in the appropriate HLU, LLU, and other support libraries. The compile line for this example using the nhlcc command is:

    nhlcc -o ncl ${NCARG_ROOT}/lib/libncl.a AddUserFuncs.c poly_w.c poly.c -L${NCARG_ROOT}/lib -lnetcdf -ldf -ll
    
AddUserFuncs.c contains the function NclAddUserFuncs and the function stubs NclAddUserFileFormats and NclAddUserHLUObjs, poly_w.c contains the wrapper function, and poly.c contains the actual function source. The extra libraries on the end are for netCDF, HDF, and lex respectively, which are used by NCL. The following is a list of all libraries needed in the order they should linked in to compile line:
    libhlu.a		HLU library
    
    libncarg.a		NCAR Graphics LLU library
    
    libncarg_gks.a		NCAR Graphics GKS library
    
    libncarg_c.a		NCAR Graphics common code
    
    libl.a			System LEX library
    
    libnetcdf.a		NetCDF library
    
    libX11.a 		X11 library
    	
    libdf.a			HDF library
    	
    libF77.a 		NOTE: This library is part of the C-Fortran interface libraries it 
    			might be different for different architectures consult the Fortran 
    			programmers manual for the machine you are using.
    
    libI77.a		NOTE: This library is part of the C-Fortran interface libraries it 
    			might be different for different architectures consult the Fortran 
    			programmers manual for the machine you are using.
    
    libU77.a		NOTE: This library is part of the C-Fortran interface libraries it 
    			might be different for different architectures consult the Fortran 
    			programmers manual for the machine you are using.
    
    libm.a 			System Math library
    
    libc.a			C library
    

----------------------------------------------------------------------

Writing your own wrapper as a shared object

To create a wrapper and compile it as a shared object, the steps of writing a C or Fortran function or procedure and writing a wrapper function are the same as described in the "Writing your own wrapper as a built-in function" section.

However, instead of creating a C function called NclAddUserFuncs, you must create a C function name "Init" which contains the same calls as NclAddUserFuncs. Then you must compile each piece of source code and link them together as a shared library. The compilation step is described in the "Using the wrapit77 application" section.

Examples

This section contains examples illustrating how to incorporate sample Fortran and C codes into NCL. These examples show how to use both the WRAPIT and wrapit77 methods to generate the shared objects.

  • Example 1 -- Fortran subroutine and function
  • Example 2 -- Fortran embedded wrapit interface blocks
  • Example 3 -- Fortran subroutine from a commercial library
  • Example 4 -- Fortran subroutine with CHARACTER input and output arguments
  • Example 5 -- Fortran subroutine with a 2-dimensional array; printing
  • Example 6 -- C subroutine and function

Example 1 -- Fortran subroutine and function

Begin with the Fortran source (in file ex01.f):

      SUBROUTINE CQUAD (A, B, C, NQ, X, QUAD)
      REAL X(NQ),QUAD(NQ)
C
C  Calculate quadratic polynomial values.
C
      DO 10 I=1,NQ
        QUAD(I) = A*X(I)**2 + B*X(I) + C
   10 CONTINUE
C
      RETURN
      END
      FUNCTION ARCLN (NUMPNT, POINTX, POINTY)
      DIMENSION POINTX(NUMPNT),POINTY(NUMPNT)
C
C  Calculate arc lengths.
C
      IF (NUMPNT .LT. 2) THEN
        PRINT *, 'ARCLN: Number of points must be at least 2'
        STOP
      ENDIF
      ARCLN = 0.
      DO 10 I=2,NUMPNT
        PDIST = SQRT((POINTX(I)-POINTX(I-1))**2 + 
     +                        (POINTY(I)-POINTY(I-1))**2)
        ARCLN = ARCLN + PDIST
   10 CONTINUE
C
      RETURN
      END

This first example follows in detail the five-step process described above.

Step 1 - define the wrapit interface block.

Create a file ex01.stub that contains the following two wrapit interface blocks:

C NCLFORTSTART
      SUBROUTINE CQUAD (A,B,C,NQ,X,QUAD)
      REAL X(NQ),QUAD(NQ)
C NCLEND
 
C NCLFORTSTART
      FUNCTION ARCLN (NUMPNT,POINTX,POINTY)
      DIMENSION POINTX(NUMPNT),POINTY(NUMPNT)
C NCLEND

Step 2 - create an NCL wrapper function.

Execute the following:

wrapit77 < ex01.stub >! ex01_W.c

to create the NCL wrapper function ex01_W.c.

Step 3 - create the object modules (the ".o" files).

Execute the following:

nhlcc -c ex01_W.c
nhlf77 -c ex01.f

to create ex01_W.o and ex01.o (on an HP system you will need to add a "+z" flag to each of the above commands).

Step 4 - create a dynamic shared object module.

For this, use the appropriate ld command listed in the section above discussing this. For example, if you are running on a Sun Solaris machine, then you would execute:

  ld -o ex01.so ex01_W.o ex01.o -G

to create the dynamic shared object ex01.so.

Note: if you are on a system that WRAPIT recognizes, then steps 2-4 could be replaced with one call to WRAPIT:

  WRAPIT ex01.stub ex01.f

Step 5 - tell NCL where your dynamic shared object is.

The external statement in the following example script tells NCL where to look for the dynamic shared object you just created.

external EX01 "./ex01.so"
 
begin
;
; calculate three values of a quadratic equation
;
   nump = 3
   x = (/ -1., 0.0, 1.0 /)
   qval = new(nump,float)              
   EX01::cquad(-1, 2, 3, nump, x, qval) ; call the new NCL version of
                                        ; your original Fortran subroutine
   print("Polynomial value = " + qval)
 
;
; calculate an arc length.
;
   xc = (/ 0., 1., 2. /)
   yc = (/ 0., 1., 0. /)
   arclen = EX01::arcln(nump,xc,yc)     ; call the new NCL version of
                                        ; your original Fortran function
   print("Arc length = " + arclen)
 
end

If you submit the above script to the NCL interpreter, it produces the output:

opening: ./ex01.so
(0)     Polynomial value = 0
(1)     Polynomial value = 3
(2)     Polynomial value = 4
(0)     Arc length = 2.82843

The numbers in parentheses at the left in the above printout are an artifact of how the NCL print function works and are of no relevance in this example.

Example 2 -- Fortran embedded wrapit interface blocks

Instead of using a separate stub file as in example 1 above, you could have used the following code (in file ex02.f) as input to wrapit77:

C NCLFORTSTART
      SUBROUTINE CQUAD (A,B,C,NQ,X,QUAD)
      REAL X(NQ),QUAD(NQ)
C NCLEND
C
C  Calculate quadratic polynomial values.
C
      DO 10 I=1,NQ
        QUAD(I) = A*X(I)**2 + B*X(I) + C
   10 CONTINUE
C
      RETURN
      END

C NCLFORTSTART
      FUNCTION ARCLN (NUMPNT, POINTX, POINTY)
      DIMENSION POINTX(NUMPNT),POINTY(NUMPNT)
C NCLEND
C
C  Calculate arc lengths.
C
      IF (NUMPNT .LT. 2) THEN
        PRINT *, 'ARCLN: Number of points must be at least 2'
        STOP
      ENDIF
      ARCLN = 0.
      DO 10 I=2,NUMPNT
        PDIST = SQRT((POINTX(I)-POINTX(I-1))**2 +
     +                        (POINTY(I)-POINTY(I-1))**2)
        ARCLN = ARCLN + PDIST
   10 CONTINUE
C
      RETURN
      END

In this example, the wrapit interface blocks are embedded directly into the Fortran code, thus avoiding the need to create a separate stub file containing them. All that wrapit77 is looking for in its input is blocks delimited by comment lines containing NCLFORTSTART and NCLEND.

Now execute:

wrapit77 < ex02.f >! ex02_W.c

and proceed as in steps 3-5 of Example 1, or execute:

WRAPIT ex02.f

and proceed with step 5.

Example 3 -- Fortran subroutine from a commercial library

Suppose you are calling:

      SUBROUTINE LIBSUB(IARG1,RARG)

from a commercial library named libcommercial.a. Using the following wrapit interface block (in file libsub.stub):

C NCLFORTSTART
      SUBROUTINE LIBSUB(IARG,RARG)
      INTEGER IARG
      REAL RARG
C NCLEND

Create an NCL wrapper function and its object module by executing:

  wrapit77 < libsub.stub > libsub_W.c
  nhlcc -c libsub_W.c

To create a dynamic shared object named libsub.so, use the appropriate ld command by loading libsub_W.o and libcommercial.a. For example, if you are running on a Sun Solaris machine, you would execute:

  ld -o libsub.so libsub_W.o -l commercial -G

to create the dynamic shared object libsub.so.

To combine the above steps, you could also use WRAPIT:

  WRAPIT libsub.stub -l commercial

Example 4 -- Fortran subroutine with CHARACTER input and output arguments

Start with a Fortran subroutine that takes an input string and returns a number of letters based on the length of the input string:

      SUBROUTINE EX04 (STRIN,STROUT)
      CHARACTER*(*) STRIN
      CHARACTER*26  ABET,STROUT
      DATA ABET/'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
C
      IMX = MIN(LEN(STRIN),26)
      STROUT = ABET(1:IMX)
C
      RETURN
      END

You could use embedded wrapit interface blocks, as in example 2 above, or use the following wrapit interface block:

C NCLFORTSTART
      SUBROUTINE EX04 (STRIN,STROUT)
      CHARACTER*(*) STRIN
      CHARACTER*26  STROUT
C NCLEND

to create the NCL wrapper and dynamic shared object (named ex04.so). Passing the following NCL script to the NCL interpreter:

external EXAMPLE04_SO "./ex04.so"
 
begin
 
 cstr = new(26,character)        ;  create a character array of length 26
 EXAMPLE04_SO::ex04("fifteen letters",cstr)
 str = chartostring(cstr)
 print(str)
 
end

produces the output:

opening: ./ex04.so
 
 
Variable: str
Type: string
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes:   [1]
Coordinates: 
(0)     ABCDEFGHIJKLMNO

Example 5 -- subroutine with a 2-dimensional array; printing

Start with one subroutine that calculates a function of two variables and stores the results in a 2-dimensional array, and another subroutine that prints 2-dimensional arrays by rows.

      SUBROUTINE EX05(M,N,X,Y,FXY)
      REAL X(M),Y(N),FXY(M,N)
C
C  Calculate FXY(I,J) = 2*I+J
C
      DO 10 J=1,N
        DO 20 I=1,M
          FXY(I,J) = 2.*REAL(I) + REAL(J)
   20   CONTINUE
   10 CONTINUE
C
      RETURN
      END
      SUBROUTINE PRT2D(M,N,A)
      REAL A(M,N)
C
C  Print the array A by rows using an F6.1 format with
C  7 values per line.
C
      DO 10 J=1,N
        PRINT *,'Row',J,':'
        DO 20 I=1,M/7
          WRITE(6,500) (A(LL,J),LL=(I-1)*7+1,I*7)
  500     FORMAT(7F6.1)
   20   CONTINUE
        IF (MOD(M,7) .NE. 0) WRITE(6,500) (A(LL,J),LL=(M/7)*7+1,M)
        PRINT *,' '
   10 CONTINUE
C
      RETURN
      END

Use the following wrapit interface block:

C NCLFORTSTART
      SUBROUTINE EX05(M,N,X,Y,FXY)
      REAL X(M),Y(N),FXY(M,N)
C NCLEND
C NCLFORTSTART
      SUBROUTINE PRT2D(M,N,A)
      REAL A(M,N)
C NCLEND

to create the NCL wrapper function and the dynamic shared object (named ex05.so). Then the following NCL script:

external EX05 "./ex05.so"
 
begin
;
; calculate three values of a quadratic equation
;
   m = 11
   n = 3
   x = new(m,float)
   y = new(n,float)
   fxy = new((/n,m/),float)
   EX05::ex05(m,n,x,y,fxy)
   EX05::prt2d(m,n,fxy)
end

will create the 2-dimensional array fxy in a manner compatible with other NCL procedures. Passing the above NCL script to the NCL interpreter produces the output:

opening: ./ex05.so
 Row  1:
     3.0   5.0   7.0   9.0  11.0  13.0  15.0
    17.0  19.0  21.0  23.0
  
 Row  2:
     4.0   6.0   8.0  10.0  12.0  14.0  16.0
    18.0  20.0  22.0  24.0
  
 Row  3:
     5.0   7.0   9.0  11.0  13.0  15.0  17.0
    19.0  21.0  23.0  25.0

Example 6 -- C subroutine and function

This example is identical to Example 1, except it's done in C, and hence you can't use WRAPIT to create the shared object. Also, the C wrapper file generated in Example 1 was used and modified to clean up the interface. The main change was to remove the requirement that the NCL script has to pass in the length of the arrays.

Here's the main "cquad.c" program with the two functions:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void cquad (double a, double b, double c, int nq, double *x, double *quad)
{
  int i;
/*
 *  Calculate quadratic polynomial values.
 */
  for(i = 0; i < nq; i++) {
    quad[i] = a * (x[i] * x[i]) + b * x[i] + c;
  }
}

double arcln (int numpnt, double *pointx, double *pointy)
{
  int i;
  double arcln, pdist;

/*
 *  Calculate arc lengths.
 */
  if (numpnt < 2) {
    printf("arcln: Number of points must be at least 2");
    exit(EXIT_FAILURE);
  }
  
  arcln = 0.;
  for(i = 1; i < numpnt; i++) {
    pdist = sqrt(pow(pointx[i]-pointx[i-1],2) + pow(pointy[i]-pointy[i-1],2));
    arcln = arcln + pdist;
  }
  return(arcln);
}

Here's the modified "ex01_W.c" file:

#include <stdio.h>
/*
 * The following are the required NCAR Graphics include files.
 * They should be located in ${NCARG_ROOT}/include
 */
#include <ncarg/hlu/hlu.h>
#include <ncarg/hlu/NresDB.h>
#include <ncarg/ncl/defs.h>
#include <ncarg/ncl/NclDataDefs.h>
#include <ncarg/ncl/NclBuiltInSupport.h>
#include <ncarg/ncl/NclBuiltIns.h>

extern void cquad(double, double, double, int, double *, double *);
extern double arcln(int, double *, double *);

NhlErrorTypes cquad_W( void ) {
  double *A, *B, *C, *X, *QUAD;
  ng_size_t X_dimsizes[1], QUAD_dimsizes[1];
  ng_size_t i, nq;
  int inq;

/*
 * Retrieve arguments from NCL function call. It doesn't matter
 * what order you call them in.
 */ 
  QUAD = (double*) NclGetArgValue(
                                  4,
                                  5,
                                  NULL,
                                  QUAD_dimsizes,
                                  NULL,
                                  NULL,
                                  NULL,
                                  1);


  X = (double*) NclGetArgValue(
                               3,
                               5,
                               NULL,
                               X_dimsizes,
                               NULL,
                               NULL,
                               NULL,
                               1);


  C = (double*) NclGetArgValue(
                               2,
                               5,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               1);


  B = (double*) NclGetArgValue(
                               1,
                               5,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               1);


  A = (double*) NclGetArgValue(
                               0,
                               5,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               NULL,
                               1);

/*
 * Get length of X and check it against the length of QUAD.
 */
  nq = X_dimsizes[0];
  inq = (int)nq;       /* This must be < INT_MAX */
  if(QUAD_dimsizes[0] != nq) {
    NhlPError(NhlFATAL,NhlEUNKNOWN,"cquad: length of QUAD must be equal to the length of X");
    return(NhlFATAL);
  }

/*
 * Call C routine
 */
  cquad(*A,*B,*C,inq,X,QUAD);

/*
 * Since this is going to be an NCL procedure, we don't need to set
 * up return value here.
 */
  return(NhlNOERROR);
}

NhlErrorTypes arcln_W( void ) {
  double *POINTX, *POINTY;
  double ARCLN_ret;
  ng_size_t POINTX_dimsizes[1], POINTY_dimsizes[1];
  ng_size_t ARCLN_ret_dimsizes[1];
  ng_size_t i, npts;
  int inpts;

/*
 * Retrieve arguments from NCL function call. It doesn't matter
 * what order you call them in.
 */ 
  POINTX = (double*) NclGetArgValue(
                                    0,
                                    2,
                                    NULL,
                                    POINTX_dimsizes,
                                    NULL,
                                    NULL,
                                    NULL,
                                    1);


  POINTY = (double*) NclGetArgValue(
                                    1,
                                    2,
                                    NULL,
                                    POINTY_dimsizes,
                                    NULL,
                                    NULL,
                                    NULL,
                                    1);


  npts = POINTX_dimsizes[0];
  npts = (int)npts;              /* npts must be < INT_MAX */  
  if(POINTY_dimsizes[0] != npts) {
    NhlPError(NhlFATAL,NhlEUNKNOWN,"arcln: POINTX and POINTY must be the same length");
    return(NhlFATAL);
  }

/*
 * Call the C function.
 */
  ARCLN_ret = arcln(inpts,POINTX,POINTY);

/*
 * Set up return value, which is a double scalar in this case.
 */ 
  ARCLN_ret_dimsizes[0] = 1;
  return(NclReturnValue(&ARCLN_ret,1,ARCLN_ret_dimsizes,NULL,NCL_double,1));
}



void Init(void){
  void *args;
  ng_size_t dimsizes[NCL_MAX_DIMENSIONS];
  int nargs;


  nargs = 0;

/*
 * Register arcln function.
 */
  args = NewArgs(2);

  dimsizes[0] = -1;      /* -1 ===> Dimension can be any size. */

  SetArgTemplate(args,0,"double",1,dimsizes);nargs++;
  SetArgTemplate(args,1,"double",1,dimsizes);nargs++;

  NclRegisterFunc(arcln_W,args,"arcln",nargs);

/*
 * Register cquad function.
 */
  nargs = 0;
  args = NewArgs(5);

  dimsizes[0] = 1;
  SetArgTemplate(args,0,"double",1,dimsizes);nargs++;
  SetArgTemplate(args,1,"double",1,dimsizes);nargs++;
  SetArgTemplate(args,2,"double",1,dimsizes);nargs++;

  dimsizes[0] = -1;      /* -1 ===> Dimension can be any size. */
  SetArgTemplate(args,3,"double",1,dimsizes);nargs++;
  SetArgTemplate(args,4,"double",1,dimsizes);nargs++;

  NclRegisterProc(cquad_W,args,"cquad",nargs);

}

Compile and create the shared object as in the instructions above. For example, on a Sun system:

  nhlcc -c cquad.c
  nhlcc -c ex01_W.c
  ld -o cquad.so ex01_W.o cquad.o -G

Now that you have the "cquad.so" shared object, you can use the following NCL script to call your new functions:

external C "./cquad.so"
 
begin
;
; Calculate three values of a quadratic equation.
;
   nump = 3
   x    = (/ -1.d, 0.0d, 1.0d /)
   qval = new(nump,double)              
   C::cquad(-1, 2, 3, x, qval) ; Call the new NCL version of
                                        ; your original C subroutine.
   print("Polynomial value = " + qval)
 
;
; Calculate an arc length.
;
   xc = (/ 0.d, 1.d, 2.d /)
   yc = (/ 0.d, 1.d, 0.d /)
   arclen = C::arcln(xc,yc)     ; Call the new NCL version of
                                        ; your original C function.
   print("Arc length = " + arclen)
end

Special considerations

This section contains several things that you should know to avoid common problems.
Variable names:
Due to a bug in the parser, no variable between the delimiters can be named "data".

Array dimensioning:
For NCL arrays, the fastest-varying dimension is the rightmost, while for Fortran it is the leftmost dimension. Therefore, if XA is a Fortran array dimensioned idim x jdim, this array will be dimensioned jdim x idim in NCL. Also, Fortran array subscripts start at 1, whereas NCL array subscripts start at 0. Example 5 in this section illustrates these concepts.

Arrays of character strings:
Currently wrapit77 honors only non-dimensioned Fortran type CHARACTER variables. You cannot pass arrays of NCL strings to Fortran, nor can you pass Fortran CHARACTER arrays from Fortran back to NCL.

Passing strings from NCL to Fortran:
If you want to pass an NCL variable of type string to a Fortran procedure, then the argument to the Fortran procedure must be declared as CHARACTER*(*). See example 4 in this section.

Passing Fortran CHARACTER variables to NCL:
If you want to pass a Fortran CHARACTER variable back to NCL, then the Fortran argument must be a variable of type CHARACTER of fixed length, and the corresponding NCL variable must be a character array of the same length. If you want to use the NCL character array as an NCL string, you will need to use the NCL conversion function chartostring. See example 4 in this section.

Complex numbers:
NCL does not have a complex data type. If you want to bring complex numbers into NCL, you will have to do it by bringing in the real and imaginary parts as separate arrays. This will most likely require that you write an interface subroutine to your Fortran code that splits up the Fortran COMPLEX numbers into real and imaginary parts. Although you will not be able to do arithmetic on the complex numbers in NCL, you can still do analysis on the real and imaginary parts separately.

Procedure name conflicts:
If the procedure that you are incorporating into NCL has the same name as a currently existing built-in NCL procedure, NCL will choose its built-in and not your procedure. However, most UNIX ld commands recognize the -B symbolic flag, and using it when you create your dynamic shared object will force NCL to load your procedure in preference to its own built-ins. The -B symbolic can cause ld to report missing entries that in fact are ultimately not missing. It is probably best just to be careful to avoid defining a procedure with the same name as an NCL built-in.

NCL termination:
If a Fortran procedure that you have incorporated into NCL executes a STOP statement, or if a C function executes an exit statement, then the NCL interpreter will abort.

Unsupported Fortran 77 syntax in wrapit interface blocks:
Fortran COMMON blocks:
Fortran COMMON blocks are not allowed in a wrapit interface block. This would preclude your having adjustable arrays whose dimensions are passed in a COMMON block, or using COMMON to pass values for variables.
Fortran ENTRY statements:
There is no way to accommodate an ENTRY statement in a Fortran procedure.
Alternate return arguments:
Subroutines with alternate return arguments are not allowed.