+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ EZspline: an f90 interface to the PSPLINE routines supporting real*4 (r4) and real*8 (r8) precisions +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ A Pletzer Mon Apr 24 09:53:34 EDT 2000 (pletzer@pppl.gov) D McCune Tue Jan 18 2005 (dmccune@pppl.gov) D McCune 2006 (dmccune@pppl.gov) D McCune Tue Apr 24 2007 (dmccune@pppl.gov) ----------------------- 2007 update -- EZspline supports "hybrid" interpolation, using methods newly added to the pspline base library. 2d and 3d objects can now be designated for different methods of interpolation along each grid dimension. Use the new methods {EZhybrid_init} to initialize r4 or r8 2d or 3d hybrid interpolation objects. This method accepts an integer array argument "hspline" which has 2 elements or 3 elements for 2d or 3d objects respectively. The interpretation is: hspline(j) = -1 => zonal step function along dimension "j" (in this case the size of the data along this dimension is ONE LESS than the corresponding grid; the grid gives the step function zone boundaries. hspline(j) = 0 => piecewise linear along dimension "j" hspline(j) = 1 => Akima Hermite C1 segmented cubic along dimension "j". hspline(j) = 2 => C2 cubic Spline along dimension "j". Restriction: at present, Akima Hermite and Spline interpolation cannot be mixed in a single hybrid object. I.e. if hspline(j)=2 then hspline(k)=1 must not occur for k.ne.j. Boundary conditions can be specified for the end points of the dimensions having cubic interpolation. The treatment is the same as for "pure" Hermite or Spline objects. However, boundary condition array sizes are reduced by one along any zonal step function dimension (matching the treatment of the interpolation data itself). ----------------------- 2006 update -- EZspline_save supports two optional arguments: character*(*) :: spl_name ! name of spline being saved. ! alphanumeric, no imbedded blanks logical :: fullsave ! TRUE to save the spline coefficients If spl_name is used, it allows multiple spline objects to be saved in a single file; the spl_name value must be given again when the file is read using EZspline_load. If the spl_name argument is omitted, the file contains but one spline object, as before. If fullsave=.TRUE. is set, the spline and all its coefficients are saved. If fullsave is omitted or .FALSE., only the data values are saved and the coefficients are recalculated at EZspline_load time. EZspline_load now supports the spl_name optional argument, to specify which spline object to read from a file containing multiple named spline objects. There is no fullsave argument; EZspline_load figures out the right thing to do based on file contents. ----------------------- 2005 update -- EZspline supports piecewise linear interpolation (linear, bilinear, trilinear for 1d 2d and 3d data respectively). Use EZlinear_Init instead of EZSpline_Init to set up a spline object for piecewise linear instead of spline or Hermite interpolation. Piecewise linear interpolation needs no boundary conditions; therefore, EZlinear_Init calls need no boundary condition type specification arguments. ----------------------- 1) AIM Provide an easy to use and flexible f90 interface to Doug McCune's PSPLINE routines. The main features of EZspline are: * provide a single interface to 1-d, 2-d or 3-d spline interpolation. * allow for a variety of boundary conditions (periodic, not-a-knot, 1st derivative or 2nd derivative imposed). * support both compact cubic splines and Akima Hermites. (In 2005, support for simple piecewise linear interpolation was added to EZspline). (In 2007, support for hybrid objects-- Spline along some dimensions, piecewise linear and/or zonal step function along others-- was added to EZspline). * offers intuitive method for evaluating derivatives up to and including 2nd order (up to 1st order for Akima Hermite and Piecewise Linear). * use dynamic memory allocation. * implement error handling through a dedicated routine (EZspline_error). * use reasonable default parameters (uniform grid [0, 1] or [0, 2 pi] in the case of periodic boundary conditions). * interpolation objects can be dumped to or loaded from a netCDF file. 2) EZspline(n)_r(m) DATA TYPE Spline and Hermite coefficients are encapsulated into derived types bearing the name EZspline(n)_r(m) where "n" denotes the dimensionality or rank, 1, 2 or 3 and "m" the precision, either 4 or 8. The "r8" and "r4" precision are defined as integer, parameter :: ezspline_r8 = selected_real_kind(12,100) integer, parameter :: ezspline_r4 = selected_real_kind(6,37) in ezspline_mod.f90. The various interpolation methods are universally supported across types (1-d, 2-d or 3-d) and precisions through the use of interfaces. The EZspline(n)_r(m) types are accessed by importing the module use EZspline_obj a statement that should be added somewhere at the top of the routine/program (before any locally defined variables and before any IMPLICIT statement). A 2-D/real*8 interpolation object, for instance, is declared as type(EZspline2_r8) :: spl Members of the 'spl'object are accessed through the % operator. Typical members are spl%isHermite Set this to 1 if interpolation is Akima (spline is default) spl%x1 Grid in the 1st variable spl%x2 Grid in the 2nd variable spl%n1 Grid size in the 1st variable ... spl%bcval1min Boundary condition value on the left for the 1st variable spl%bcval1max Boundary condition value on the right for the 1st variable spl%bcval2min Boundary condition value on the left for the 2nd variable ... Some members such as x1, x2, x3..., bcval1min, bcval1max, bcval2min... can be explicitely set by the user. Others, such as n1, n2... and the boundary value type ibctype1, ibctype2... should only be set through a call to the constructor (EZspline_init). 3) EZspline's REFLECTIVE METHODS EZspline methods defined in the EZspline module. To access the the methods the following modules must be imported use EZspline_obj use EZspline somewhere at the top of the caller routine (before any local declaration and any IMPLICIT statement). The interpolation (spline or Akima) object is declared as follows: type(EZspline3_r8) :: spline_o ! 3-d cubic spline/real*8 object (somewhere at the top but below any USE and IMPLICIT declarations). All reflective methods are subroutines which take an interpolation object as 1st argument and return an error flag as last argument (0 is OK). It is strongly advised to to check the value of the returned flag (ier) after each method is invoked and map its value to an error message by calling EZspline_error(ier). call EZspline_xxxx(spline_o, ...., ier) ! generic EZspline method call call EZspline_error(ier) At the very least, the use of an EZspline interpolation involves the following steps: * initialization (EZspline_init or EZlinear_init or EZhybrid_init) * coefficient set-up (EZspline_setup) * interpolation proper (eg EZspline_interp) * freeing the memory (EZspline_free) We now go through every method in details 3.a) Initialization, memory allocation and boundary condition set-up subroutine EZspline_init(spline_o, n1, n2, n3, BCS1, BCS2, BCS3, ier) type(EZspline3_r8) spline_o integer, intent(in) :: n1, n2, n3 integer, intent(in) :: BCS1(2), BCS2(2), BCS3(2) integer, intent(out) :: ier Here BCS1, BCS2 and BCS3 are 2-element arrays describing the boundary conditions. Each element can be either: * -1 for periodic boundary condition * 0 for not a knot boundary condition * 1 if the 1st derivative is imposed * 2 if the 2nd derivative is imposed on grid n=1, 2 or 3 with the 1st and 2nd elements denoting the boundary conditions on the left and right respectively. Note-- Uniform grids [0,1] are assumed for x1, x2, x3 by default ([0, 2 pi] in the case of periodic boundary conditions). It is the responsibility of the user to explicitly set x1, x2 and x3 otherwise! The x1, x2 and x3 members should be set before EZspline_setup is called. To initialize an 3d object for piecewise linear interpolation, use instead: subroutine EZlinear_init(spline_o, n1, n2, n3, ier) type(EZspline3_r8) spline_o integer, intent(in) :: n1, n2, n3 integer, intent(out) :: ier In this case the default uniform grids always cover [0,1]. To initialize a 3d hybrid object, use instead: subroutine EZhybrid_init3_r8(spline_o, n1, n2, n3, hspline, ier, & BCS1, BCS2, BCS3) type(EZspline3_r8) spline_o integer, intent(in) :: n1, n2, n3 integer, intent(in) :: hspline(3) integer, intent(out) :: ier integer, intent(in), OPTIONAL :: BCS1(2), BCS2(2), BCS3(2) The integer array hspline(1:3) controls the interpolation to be used along each dimension, according to the code: hspline(j) = -1 => zonal step function along dimension "j" (in this case the size of the data along this dimension is ONE LESS than the corresponding grid; the grid gives the step function zone boundaries. hspline(j) = 0 => piecewise linear along dimension "j" hspline(j) = 1 => Akima Hermite C1 segmented cubic along dimension "j". hspline(j) = 2 => C2 cubic Spline along dimension "j". The boundary condition specifications BCS* are optional; if present they are interpreted as in EZspline_init. Boundary condition controls which prescribe a derivative value only have effect for the end points of dimensions along which cubic (Hermite or Spline) interpolation is applied. A positive boundary condition control for a dimension with piecewise linear or zonal step function interpolation will be reported as an error. For piewise linear or zonal dimensions, a periodic boundary condition can be specified, but, it only affects the default grid provided for that dimension. As with EZspline, default grids are provided that span [0:1], or [0:2pi] for dimensions with periodic boundary conditions specified. NOTE: at present, Spline and Hermite interpolation cannot be mixed, for technical reasons. An error will be reported if one of the hspline(1:3) elements is 1, and another is 2. 3.b) Spline/Akima Hermite/Hybrid object coefficients set-up Given the array f this routine computes all the necessary cubic coefficients for subsequent interpolation. subroutine EZspline_setup(spline_o, f, ier, exact_dim) type(EZspline3_r8) spline_o real(r8), dimension(:,:,:), intent(in) :: f integer, intent(out) :: ier logical, intent(in), OPTIONAL :: exact_dim Even for piecewise linear interpolation, this routine must be called-- this copies the data being interpolated into the object, and, calculates information needed to speed zone lookup during future interpolation calls. The optional argument controls array dimension size checking. Note that the expected size is the grid size (the n1, n2, or n3 arguments given in EZspline_init), or one less than the grid size in the case of dimensions of hybrid objects along which zonal step function interpolation is used. If exact_dim=.TRUE. is specified, all the dimensions of f(:,:,:) must match the corresponding non-coefficient dimensions of spline_o%fspl EXACTLY. If it is omitted or set .FALSE., the dimensions of f(:,:,:) can match or exceed the corresponding dimensions of spline_o%fspl; only the first (n1,n2,n3, or n1-1,n2-1,n3-1 as the case may be) points of f are copied into spline_o%fspl. 3.c) Interpolation methods The interpolation methods come in 3 flavors: point interpolation, cloud interpolation and array interpolation. The point interpolation operates on a single point (p1, p2, p3). The cloud interpolation operates on a list of unordered points, and is particularly useful for mapping a structured mesh onto an unstructured grid. The array interpolation operates on a structured grid. All interpolation routines share the same name but accept 3 types of arguments through the use of an interface: ! point interpolation: p1, p2, p3 and f are scalars call EZspline_interp(spline_o, p1, p2, p3, f, ier) ! ! cloud interpolation: p1, p2, p3 and f are rank-1 arrays of length n call EZspline_interp(spline_o, n, p1, p2, p3, f, ier) ! ! array interpolation: p1, p2, p3 are rank-1 arrays of ! length n1, n2 and n3 respectively. Here f has size (n1, n2, n3). call EZspline_interp(spline_o, n1, n2, n3, p1, p2, p3, f, ier) 3.d) Higher order derivatives These are achieved using EZspline_derivative with the first 3 integers denoting the derivative order. As above, the EZspline_derivative routines come in 3 flavors: point, cloud and array interpolation. ! ! evaluate the spline derivative ! d^{i1} d^{i2} d^{i3} f / d x1^{i1} d x2^{i2} d x2^{i3} ! subroutine EZspline_derivative(spline_o, i1, i2, i3, p1, p2, p3, f, ier) In some situations all derivatives may be required and it is then more efficient to call ! ! Return the gradient at point p1, p2, p3 in df(3)=(/fx, fy, fz/). ! subroutine EZspline_gradient(spline_o, p1, p2, p3, df, ier) 3.e) Control/diagnostics Because checking whether arguments are inside a given domain can be detrimental to efficiency, the responsibility to check this is left to the user. A return value of ier > 0 indicates failure of the test. ! ! control/diagnostics ! subroutine EZspline_isInDomain(spline_o, p1, p2, p3, ier) subroutine EZspline_isGridRegular(spline_o, ier) 3.f) persistence Persistence is the action of saving intermediate results for debugging and/or postprocessing purposes. Here, EZcdf_save saves all the data necessary to rebuild the cubic spline object, and EZcdf_load reads these data and recompute the cubic spline coefficients. ! ! save in or load from netCDF file -- ! note new arguments "spl_name" and "fullsave" which are OPTIONAL: ! see "2006 update" comments, above. ! subroutine EZspline_save(spline_o, filname, ier, & spl_name, fullsave) ! optional subroutine EZspline_load(spline_o, filname, ier, & spl_name) ! optional ! 3.g) modulo ! Map argument to (x1,2,3min, x1,2,3max) interval. This is useful to avoid ! an out-of-grid error when periodic boundary conditions are applied. This ! method has no effect when the boundary conditions are not periodic. ! subroutine EZspline_modulo3(spline_o, p1, p2, p3, ier) 3.g) free Every call to EZcdf_init should be followed by a call to EZspline_free to prevent memory leak. ! deallocate, reset etc ! call EZspline_free(spline_o, ier) Note that an error will be reported if an attempt is made to initialize an already-initialized spline object (which would be a sign of a memory leak). 4) EZspline's NON-REFLECTIVE METHODS A small number of methods do not take the interpolation object as first argument. Hence, these methods can be called even after the object has been freed. 4.a) Error handling The meaning a an error code is printed using EZspline_error. This routine is a look up table printing error/warning messages. ! ! error handle routine ! call EZspline_error(ier) 4.b) Save data in netCDF file format Save data in netCDF file 'filename'. To save an EZspline1,2,3 object use EZspline_save method. subroutine EZspline_2NetCDF_array3(n1, n2, n3, x1, x2, x3, f, filename, ier) 5) EXAMPLE Here is a full 3-d spline interpolation example taking periodic boundary conditions on the first two variables, and not a knot boundary condition on the last variable. Accordingly, the bcval1,2,3min/bcval1,2,3max values need not be set. program spline_test ! example of ezspline calls ! ------------------------- ! On PPPL cluster machine use the following command ot compile: ! On Alpha-Linux: ! f90 -assume no2underscores -I/usr/ntcc/mod -o spline_test spline_test.f90 -L/usr/ntcc/lib -lpspline -lezcdf -L/usr/local/lib -lnetcdf ! On Alpha-OSF1: ! f90 -I/usr/ntcc/mod -o spline_test spline_test.f90 -L/usr/ntcc/lib -lpspline -lezcdf -L/usr/local/lib -lnetcdf ! On Solaris: ! f90 -M/usr/ntcc/mod -dalign -o spline_test spline_test.f90 -L/usr/ntcc/lib -lpspline -lezcdf -L/usr/local/lib -lnetcdf ! On Linux (Lahey-Fujitsu): ! lf95 -I /usr/ntcc/lff95/mod -o spline_test spline_test.f90 -L/usr/ntcc/lff95/lib -lpspline -lezcdf -L/usr/local/lib -lnetcdf use EZspline_obj ! import the modules use EZspline implicit none integer, parameter :: r8 = selected_real_kind(12,100) ! real*8 kind real(r8), parameter :: twopi = 6.2831853071795862320_r8 real(r8), dimension(:), allocatable :: x1, x2, x3 ! independent ! variables real(r8), dimension(:,:,:), allocatable :: f ! function values integer n1, n2, n3, ier, BCS1(2), BCS2(2), BCS3(2) type(EZspline3_r8) :: spline_o ! 3-d EZspline object real(r8) p1, p2, p3, fp integer m, i1, i2, i3, i real(r8), dimension(:), allocatable :: y1, y2, y3, fy integer k1, k2, k3 real(r8), dimension(:), allocatable :: z1, z2, z3 real(r8), dimension(:,:,:), allocatable :: fz n1 = 11 n2 = 6 n3 = 21 allocate(x1(n1), x2(n2), x3(n3), f(n1,n2,n3)) BCS1 = (/ -1, -1 /) ! periodic BCS2 = (/ -1, -1 /) ! periodic BCS3 = (/ 0, 0 /) ! not a knot ! initialize/allocate memory print *,'initializing...' call EZspline_init(spline_o, n1, n2, n3, BCS1, BCS2, BCS3, ier) call EZspline_error(ier) !spline_o%x1 = ... ! necessary if spline_o%x1 not in [0, 2 pi] !spline_o%x2 = ... ! necessary if spline_o%x2 not in [0, 2 pi] ! set explicitely spline_o%x3 spline_o%x3 = -1._r8 + 2._r8*(/ ( (real(i-1,r8)/real(n3-1,r8))**2, i=1, n3 ) /) ! need to set explicitly the following if boundary conditions ! are anything but not-a-knot or periodic, i.e. BCS(n) has ! element /= -1, 0. ! spline_o%bcval1min = ... ; spline_o%bcval1max = ... ! spline_o%bcval2min = ... ; spline_o%bcval2max = ... ! spline_o%bcval3min = ... ; spline_o%bcval3max = ... ! compute cubic spline coefficients do i3=1, n3 do i2 = 1, n2 do i1 = 1, n1 f(i1,i2,i3) = (cos(twopi*i1/(n1-1))* & & sin(twopi*i2/(n2-1))**2)*exp(spline_o%x3(i3)) enddo enddo enddo print *,'setting up...' call EZspline_setup(spline_o, f, ier) call EZspline_error(ier) ! save object call EZspline_save(spline_o, "spline.nc", ier) call EZspline_error(ier) ! point interpolation p1 = twopi/2._r8 p2 = 0._r8 p3 = -0.5_r8 print *,'point interpolation...' call EZspline_interp(spline_o, p1, p2, p3, fp, ier) call EZspline_error(ier) ! cloud interpolation m = 21 allocate(y1(m), y2(m), y3(m), fy(m)) y1 = spline_o%x1min + (spline_o%x1max-spline_o%x1min)* & & (/ ( real(i-1,r8)/real(m-1,r8), i=1, m ) /) y2 = spline_o%x2min + (spline_o%x2max-spline_o%x2min)* & & (/ ( real(i-1,r8)/real(m-1,r8), i=1, m ) /) y3 = spline_o%x3min + (spline_o%x3max-spline_o%x3min)* & & (/ ( real(i-1,r8)/real(m-1,r8), i=1, m ) /) print *,'cloud interpolation...' call EZspline_interp(spline_o, m, y1, y2, y3, fy, ier) call EZspline_error(ier) ! array interpolation k1 = 5 k2 = 4 k3 = 3 allocate(z1(k1), z2(k2), z3(k3), fz(k1,k2,k3)) z1 = spline_o%x1min + (spline_o%x1max-spline_o%x1min)* & & (/ ( real(i-1,r8)/real(k1-1,r8), i=1, k1 ) /) z2 = spline_o%x2min + (spline_o%x2max-spline_o%x2min)* & & (/ ( real(i-1,r8)/real(k2-1,r8), i=1, k2 ) /) z3 = spline_o%x3min + (spline_o%x3max-spline_o%x3min)* & & (/ ( real(i-1,r8)/real(k3-1,r8), i=1, k3 ) /) print *,'array interpolation...' call EZspline_interp(spline_o, k1, k2, k3, z1, z2, z3, fz, ier) call EZspline_error(ier) ! clean up and free up memory print *,'cleaning up' call Ezspline_free(spline_o, ier) call EZspline_error(ier) deallocate(x1, x2, x3, f) deallocate(y1, y2, y3, fy) deallocate(z1, z2, z3, fz) end program spline_test