diff --git a/czspline/czhybrid_init.f90 b/czspline/czhybrid_init.f90
deleted file mode 100644
index 5a32f59..0000000
--- a/czspline/czhybrid_init.f90
+++ /dev/null
@@ -1,42 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Initialization of ezspline
-
-#include "czspline_handle_size.h"
-
-subroutine czhybrid_init2(handle, n1, n2, hspline, BCS1, BCS2, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1, n2
-  integer, intent(in) :: hspline(2)
-  integer, intent(in) :: BCS1(2), BCS2(2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  allocate(self % ptr)
-  call ezhybrid_init(self % ptr, n1, n2, hspline, ier, BCS1, BCS2)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czhybrid_init2
-
-subroutine czhybrid_init3(handle, n1, n2, n3, hspline, BCS1, BCS2, BCS3, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1, n2, n3
-  integer, intent(in) :: hspline(3)
-  integer, intent(in) :: BCS1(2), BCS2(2), BCS3(2)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  allocate(self % ptr)
-  call ezhybrid_init(self % ptr, n1, n2, n3, hspline, ier, BCS1, BCS2, BCS3)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czhybrid_init3
diff --git a/czspline/czspline_derivative.f90 b/czspline/czspline_derivative.f90
deleted file mode 100644
index e5b1430..0000000
--- a/czspline/czspline_derivative.f90
+++ /dev/null
@@ -1,167 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Compute derivatives
-
-#include "czspline_handle_size.h"
-
-! point
-subroutine czspline_derivative1(handle, i1, p1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1
-  real(fp), intent(in) :: p1
-  real(fp), intent(out) :: f
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, p1, f, ier)
-end subroutine czspline_derivative1
-
-! cloud
-subroutine czspline_derivative1_cloud(handle, i1, k, p1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1
-  integer, intent(in) :: k
-  real(fp), intent(in) :: p1(k)
-  real(fp), intent(out) :: f(k)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, k, p1, f, ier)
-end subroutine czspline_derivative1_cloud
-
-!array
-subroutine czspline_derivative1_array(handle, i1, k1, p1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1
-  integer, intent(in) :: k1
-  real(fp), intent(in) :: p1(k1)
-  real(fp), intent(out) :: f(k1)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, k1, p1, f, ier)
-end subroutine czspline_derivative1_array
-
-! point
-subroutine czspline_derivative2(handle, i1, i2, p1, p2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1, i2
-  real(fp), intent(in) :: p1, p2
-  real(fp), intent(out) :: f
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, i2, p1, p2, f, ier)
-end subroutine czspline_derivative2
-
-! cloud
-subroutine czspline_derivative2_cloud(handle, i1, i2, k, p1, p2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1, i2
-  integer, intent(in) :: k
-  real(fp), intent(in) :: p1(k), p2(k)
-  real(fp), intent(out) :: f(k)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, i2, k, p1, p2, f, ier)
-end subroutine czspline_derivative2_cloud
-
-!array
-subroutine czspline_derivative2_array(handle, i1, i2, k1, k2, p1, p2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1, i2
-  integer, intent(in) :: k1, k2
-  real(fp), intent(in) :: p1(k1), p2(k2)
-  real(fp), intent(out) :: f(k1, k2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, i2, k1, k2, p1, p2, f, ier)
-end subroutine czspline_derivative2_array
-
-! point
-subroutine czspline_derivative3(handle, i1, i2, i3, p1, p2, p3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1, i2, i3
-  real(fp), intent(in) :: p1, p2, p3
-  real(fp), intent(out) :: f
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, i2, i3, p1, p2, p3, f, ier)
-end subroutine czspline_derivative3
-
-! cloud
-subroutine czspline_derivative3_cloud(handle, i1, i2, i3, k, p1, p2, p3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1, i2, i3
-  integer, intent(in) :: k
-  real(fp), intent(in) :: p1(k), p2(k), p3(k)
-  real(fp), intent(out) :: f(k)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, i2, i3, k, p1, p2, p3, f, ier)
-end subroutine czspline_derivative3_cloud
-
-!array
-subroutine czspline_derivative3_array(handle, i1, i2, i3, k1, k2, k3, p1, p2, p3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: i1, i2, i3
-  integer, intent(in) :: k1, k2, k3
-  real(fp), intent(in) :: p1(k1), p2(k2), p3(k3)
-  real(fp), intent(out) :: f(k1, k2, k3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_derivative(self % ptr, i1, i2, i3, k1, k2, k3, p1, p2, p3, f, ier)
-end subroutine czspline_derivative3_array
diff --git a/czspline/czspline_free.f90 b/czspline/czspline_free.f90
deleted file mode 100644
index 61c7822..0000000
--- a/czspline/czspline_free.f90
+++ /dev/null
@@ -1,51 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Finalization of ezspline
-  
-#include "czspline_handle_size.h"
-
-subroutine czspline_free1(handle, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_free(self %ptr, ier)
-  deallocate(self %ptr)
-  handle = 0
-end subroutine czspline_free1
-
-subroutine czspline_free2(handle, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_free(self %ptr, ier)
-  deallocate(self %ptr)
-  handle = 0
-end subroutine czspline_free2
-
-subroutine czspline_free3(handle, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_free(self %ptr, ier)
-  deallocate(self %ptr)
-  handle = 0
-end subroutine czspline_free3
-
diff --git a/czspline/czspline_gradient.f90 b/czspline/czspline_gradient.f90
deleted file mode 100644
index 6b0a5c3..0000000
--- a/czspline/czspline_gradient.f90
+++ /dev/null
@@ -1,160 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Compute gradients
-  
-#include "czspline_handle_size.h"
-
-! point
-subroutine czspline_gradient1(handle, p1, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1
-  real(fp), intent(out) :: df
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, p1, df, ier)
-end subroutine czspline_gradient1
-
-! cloud
-subroutine czspline_gradient1_cloud(handle, k, p1, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k
-  real(fp), intent(in) :: p1(k)
-  real(fp), intent(out) :: df(k)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, k, p1, df, ier)
-end subroutine czspline_gradient1_cloud
-
-! array
-subroutine czspline_gradient1_array(handle, k1, p1, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1
-  real(fp), intent(in) :: p1(k1)
-  real(fp), intent(out) :: df(k1)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, k1, p1, df, ier)
-end subroutine czspline_gradient1_array
-
-! point
-subroutine czspline_gradient2(handle, p1, p2, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1, p2
-  real(fp), intent(out) :: df(2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, p1, p2, df, ier)
-end subroutine czspline_gradient2
-
-! cloud
-subroutine czspline_gradient2_cloud(handle, k, p1, p2, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k
-  real(fp), intent(in) :: p1(k), p2(k)
-  real(fp), intent(out) :: df(k, 2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, k, p1, p2, df, ier)
-end subroutine czspline_gradient2_cloud
-
-! array
-subroutine czspline_gradient2_array(handle, k1, k2, p1, p2, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1, k2
-  real(fp), intent(in) :: p1(k1), p2(k2)
-  real(fp), intent(out) :: df(k1, k2, 2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, k1, k2, p1, p2, df, ier)
-end subroutine czspline_gradient2_array
-
-! point
-subroutine czspline_gradient3(handle, p1, p2, p3, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1, p2, p3
-  real(fp), intent(out) :: df(3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, p1, p2, p3, df, ier)
-end subroutine czspline_gradient3
-
-! cloud
-subroutine czspline_gradient3_cloud(handle, k, p1, p2, p3, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k
-  real(fp), intent(in) :: p1(k), p2(k), p3(k)
-  real(fp), intent(out) :: df(k, 3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, k, p1, p2, p3, df, ier)
-end subroutine czspline_gradient3_cloud
-
-! array
-subroutine czspline_gradient3_array(handle, k1, k2, k3, p1, p2, p3, df, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1, k2, k3
-  real(fp), intent(in) :: p1(k1), p2(k2), p3(k3)
-  real(fp), intent(out) :: df(k1, k2, k3, 3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_gradient(self % ptr, k1, k2, k3, p1, p2, p3, df, ier)
-end subroutine czspline_gradient3_array
-
-
diff --git a/czspline/czspline_handle_size.h b/czspline/czspline_handle_size.h
deleted file mode 100644
index bcab5a4..0000000
--- a/czspline/czspline_handle_size.h
+++ /dev/null
@@ -1 +0,0 @@
-#define _ARRSZ 12 
diff --git a/czspline/czspline_init.f90 b/czspline/czspline_init.f90
deleted file mode 100644
index 0b50658..0000000
--- a/czspline/czspline_init.f90
+++ /dev/null
@@ -1,56 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Initialization of ezspline
-
-#include "czspline_handle_size.h"
-
-subroutine czspline_init1(handle, n1, BCS1, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1
-  integer, intent(in) :: BCS1(2)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  allocate(self % ptr)
-  call ezspline_init(self % ptr, n1, BCS1, ier)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czspline_init1
-
-subroutine czspline_init2(handle, n1, n2, BCS1, BCS2, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1, n2
-  integer, intent(in) :: BCS1(2), BCS2(2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  allocate(self % ptr)
-  call ezspline_init(self % ptr, n1, n2, BCS1, BCS2, ier)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czspline_init2
-
-subroutine czspline_init3(handle, n1, n2, n3, BCS1, BCS2, BCS3, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1, n2, n3
-  integer, intent(in) :: BCS1(2), BCS2(2), BCS3(2)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  allocate(self % ptr)
-  call ezspline_init(self % ptr, n1, n2, n3, BCS1, BCS2, BCS3, ier)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czspline_init3
diff --git a/czspline/czspline_interp.f90 b/czspline/czspline_interp.f90
deleted file mode 100644
index 7e71db6..0000000
--- a/czspline/czspline_interp.f90
+++ /dev/null
@@ -1,158 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Interpolate
-
-#include "czspline_handle_size.h"
-
-! point interpolation
-subroutine czspline_interp1(handle, p1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1
-  real(fp), intent(out) :: f
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, p1, f, ier)
-end subroutine czspline_interp1
-
-! cloud interpolation
-subroutine czspline_interp1_cloud(handle, k, p1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k
-  real(fp), intent(in) :: p1(k)
-  real(fp), intent(out) :: f(k)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, k, p1, f, ier)
-end subroutine czspline_interp1_cloud
-
-! array interpolation
-subroutine czspline_interp1_array(handle, k1, p1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1 
-  real(fp), intent(in) :: p1(k1)
-  real(fp), intent(out) :: f(k1)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, k1, p1, f, ier)
-end subroutine czspline_interp1_array
-
-! point interpolation
-subroutine czspline_interp2(handle, p1, p2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1, p2
-  real(fp), intent(out) :: f
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, p1, p2, f, ier)
-end subroutine czspline_interp2
-
-! cloud interpolation
-subroutine czspline_interp2_cloud(handle, k, p1, p2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k
-  real(fp), intent(in) :: p1(k), p2(k)
-  real(fp), intent(out) :: f(k)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, k, p1, p2, f, ier)
-end subroutine czspline_interp2_cloud
-
-! array interpolation
-subroutine czspline_interp2_array(handle, k1, k2, p1, p2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1, k2
-  real(fp), intent(in) :: p1(k1), p2(k2)
-  real(fp), intent(out) :: f(k1, k2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, k1, k2, p1, p2, f, ier)
-end subroutine czspline_interp2_array
-
-! point interpolation
-subroutine czspline_interp3(handle, p1, p2, p3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1, p2, p3
-  real(fp), intent(out) :: f
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, p1, p2, p3, f, ier)
-end subroutine czspline_interp3
-
-! cloud interpolation
-subroutine czspline_interp3_cloud(handle, k, p1, p2, p3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k
-  real(fp), intent(in) :: p1(k), p2(k), p3(k)
-  real(fp), intent(out) :: f(k)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, k, p1, p2, p3, f, ier)
-end subroutine czspline_interp3_cloud
-
-! array interpolation
-subroutine czspline_interp3_array(handle, k1, k2, k3, p1, p2, p3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1, k2, k3
-  real(fp), intent(in) :: p1(k1), p2(k2), p3(k3)
-  real(fp), intent(out) :: f(k1, k2 ,k3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_interp(self % ptr, k1, k2, k3, p1, p2, p3, f, ier)
-end subroutine czspline_interp3_array
diff --git a/czspline/czspline_isgridregular.f90 b/czspline/czspline_isgridregular.f90
deleted file mode 100644
index ce5218c..0000000
--- a/czspline/czspline_isgridregular.f90
+++ /dev/null
@@ -1,44 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Return error if original grid is not strictly increasing
-
-#include "czspline_handle_size.h"
-
-subroutine czspline_isgridregular1(handle, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_isgridregular(self % ptr, ier)
-end subroutine czspline_isgridregular1
-
-subroutine czspline_isgridregular2(handle, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_isgridregular(self % ptr, ier)
-end subroutine czspline_isgridregular2
-
-subroutine czspline_isgridregular3(handle, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_isgridregular(self % ptr, ier)
-end subroutine czspline_isgridregular3
diff --git a/czspline/czspline_isindomain.f90 b/czspline/czspline_isindomain.f90
deleted file mode 100644
index 15394ea..0000000
--- a/czspline/czspline_isindomain.f90
+++ /dev/null
@@ -1,148 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Return error if any of the target points lie outside the domain
-
-#include "czspline_handle_size.h"
-
-!point
-subroutine czspline_isindomain1(handle, p1, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, p1, ier)
-end subroutine czspline_isindomain1
-
-!cloud
-subroutine czspline_isindomain1_cloud(handle, k, p1, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: k
-  real(fp), intent(in) :: p1(k)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, k, p1, ier)
-end subroutine czspline_isindomain1_cloud
-
-!array
-subroutine czspline_isindomain1_array(handle, k1, p1, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: k1
-  real(fp), intent(in) :: p1(k1)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, k1, p1, ier)
-end subroutine czspline_isindomain1_array
-
-!point
-subroutine czspline_isindomain2(handle, p1, p2, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1, p2
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, p1, p2, ier)
-end subroutine czspline_isindomain2
-
-!cloud
-subroutine czspline_isindomain2_cloud(handle, k, p1, p2, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: k
-  real(fp), intent(in) :: p1(k), p2(k)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, k, p1, p2, ier)
-end subroutine czspline_isindomain2_cloud
-
-!array
-subroutine czspline_isindomain2_array(handle, k1, k2, p1, p2, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: k1, k2
-  real(fp), intent(in)  :: p1(k1), p2(k2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, k1, k2, p1, p2, ier)
-end subroutine czspline_isindomain2_array
-
-!point
-subroutine czspline_isindomain3(handle, p1, p2, p3, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in) :: p1, p2, p3
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, p1, p2, p3, ier)
-end subroutine czspline_isindomain3
-
-!cloud
-subroutine czspline_isindomain3_cloud(handle, k, p1, p2, p3, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: k
-  real(fp), intent(in) :: p1(k), p2(k), p3(k)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, k, p1, p2, p3, ier)
-end subroutine czspline_isindomain3_cloud
-
-!array
-subroutine czspline_isindomain3_array(handle, k1, k2, k3, p1, p2, p3, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: k1, k2, k3
-  real(fp), intent(in) :: p1(k1), p2(k2), p3(k3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_isindomain(self % ptr, k1, k2, k3, p1, p2, p3, ier)
-end subroutine czspline_isindomain3_array
diff --git a/czspline/czspline_load.f90 b/czspline/czspline_load.f90
deleted file mode 100644
index 2c5e3ee..0000000
--- a/czspline/czspline_load.f90
+++ /dev/null
@@ -1,53 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Load object from file
-
-#include "czspline_handle_size.h"
-
-subroutine czspline_load1(handle, filename, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  character(len=*), intent(in) :: filename
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  allocate(self % ptr)
-  call ezspline_load(self % ptr, filename, ier)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czspline_load1
-
-subroutine czspline_load2(handle, filename, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  character(len=*), intent(in) :: filename
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  allocate(self % ptr)
-  call ezspline_load(self % ptr, filename, ier, ' ')
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czspline_load2
-
-subroutine czspline_load3(handle, filename, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  character(len=*), intent(in) :: filename
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  allocate(self % ptr)
-  call ezspline_load(self % ptr, filename, ier)
-  handle = 0
-  handle = transfer(self, handle)
-end subroutine czspline_load3
diff --git a/czspline/czspline_pointer_types.f90 b/czspline/czspline_pointer_types.f90
deleted file mode 100644
index aa642c4..0000000
--- a/czspline/czspline_pointer_types.f90
+++ /dev/null
@@ -1,25 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Boiler plate code to generate pointer to types for 1, 2, and 3 dimensions
-
-module czspline_pointer_types
-
-  use ezspline_obj
-  use ezspline
-
-  type czspline1
-    type(ezspline1), pointer :: ptr => NULL()
-  end type czspline1
-
-  type czspline2
-    type(ezspline2), pointer :: ptr => NULL()
-  end type czspline2
-
-  type czspline3
-    type(ezspline3), pointer :: ptr => NULL()
-  end type czspline3
-
-end module czspline_pointer_types
diff --git a/czspline/czspline_save.f90 b/czspline/czspline_save.f90
deleted file mode 100644
index fb2fbb0..0000000
--- a/czspline/czspline_save.f90
+++ /dev/null
@@ -1,47 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Save object in file
-
-#include "czspline_handle_size.h"
-
-subroutine czspline_save1(handle, filename, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  character(len=*), intent(in) :: filename
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_save(self % ptr, filename, ier)
-end subroutine czspline_save1
-
-subroutine czspline_save2(handle, filename, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  character(len=*), intent(in) :: filename
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_save(self % ptr, filename, ier)
-end subroutine czspline_save2
-
-subroutine czspline_save3(handle, filename, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  character(len=*), intent(in) :: filename
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_save(self % ptr, filename, ier)
-end subroutine czspline_save3
diff --git a/czspline/czspline_setters.f90 b/czspline/czspline_setters.f90
deleted file mode 100644
index ca41f6b..0000000
--- a/czspline/czspline_setters.f90
+++ /dev/null
@@ -1,204 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Member setters
-
-#include "czspline_handle_size.h"
-
-! x1, x2, x3
-subroutine czspline_set_axes1(handle, k1, x1, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1
-  real(fp), intent(in) :: x1(k1)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % x1 = x1
-end subroutine czspline_set_axes1
-
-subroutine czspline_set_axes2(handle, k1, k2, x1, x2, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1, k2
-  real(fp), intent(in) :: x1(k1), x2(k2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % x1 = x1
-  self % ptr % x2 = x2
-end subroutine czspline_set_axes2
-
-subroutine czspline_set_axes3(handle, k1, k2, k3, x1, x2, x3, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: k1, k2, k3
-  real(fp), intent(in) :: x1(k1), x2(k2), x3(k3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % x1 = x1
-  self % ptr % x2 = x2
-  self % ptr % x3 = x3
-end subroutine czspline_set_axes3
-
-! isHermite
-subroutine czspline_set_ishermite1(handle, flag, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: flag
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % isHermite = flag
-end subroutine czspline_set_ishermite1
-
-subroutine czspline_set_ishermite2(handle, flag, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: flag
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % isHermite = flag
-end subroutine czspline_set_ishermite2
-
-subroutine czspline_set_ishermite3(handle, flag, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)  :: flag
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % isHermite = flag
-end subroutine czspline_set_ishermite3
-
-! Boundary types
-subroutine czspline_set_bctypes1(handle, bctype1, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)    :: bctype1(2)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % ibctype1 = bctype1
-end subroutine czspline_set_bctypes1
-
-subroutine czspline_set_bctypes2(handle, bctype1, bctype2, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)    :: bctype1(2), bctype2(2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % ibctype1 = bctype1
-  self % ptr % ibctype2 = bctype2
-end subroutine czspline_set_bctypes2
-
-subroutine czspline_set_bctypes3(handle, bctype1, bctype2, bctype3, ier)
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in)    :: bctype1(2), bctype2(2), bctype3(2)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % ibctype1 = bctype1
-  self % ptr % ibctype2 = bctype2
-  self % ptr % ibctype3 = bctype3
-end subroutine czspline_set_bctypes3
-
-! Boundary conditions bcval1min, bcval1max, ...
-subroutine czspline_set_bcvals1(handle, bcval1, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in)   :: bcval1(2)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % bcval1min = bcval1(1)
-  self % ptr % bcval1max = bcval1(2)
-end subroutine czspline_set_bcvals1
-
-subroutine czspline_set_bcvals2(handle, bcval1, bcval2, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in)   :: bcval1(2), bcval2(2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % bcval1min = bcval1(1)
-  self % ptr % bcval1max = bcval1(2)
-  self % ptr % bcval2min = bcval2(1)
-  self % ptr % bcval2max = bcval2(2)
-end subroutine czspline_set_bcvals2
-
-subroutine czspline_set_bcvals3(handle, bcval1, bcval2, bcval3, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  real(fp), intent(in)   :: bcval1(2), bcval2(2), bcval3(2)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  ier = 0
-  self % ptr % bcval1min = bcval1(1)
-  self % ptr % bcval1max = bcval1(2)
-  self % ptr % bcval2min = bcval2(1)
-  self % ptr % bcval2max = bcval2(2)
-  self % ptr % bcval3min = bcval3(1)
-  self % ptr % bcval3max = bcval3(2)
-end subroutine czspline_set_bcvals3
diff --git a/czspline/czspline_setup.f90 b/czspline/czspline_setup.f90
deleted file mode 100644
index 7b38676..0000000
--- a/czspline/czspline_setup.f90
+++ /dev/null
@@ -1,53 +0,0 @@
-!---------------------------------------------------------------------------
-! This code was developed at Tech-X (www.txcorp.com). It is free for any one
-! to use but comes with no warranty whatsoever. Use at your own risk. 
-! Thanks for reporting bugs to pletzer@txcorp.com. 
-!---------------------------------------------------------------------------
-! Compute spline coefficients
-
-#include "czspline_handle_size.h"
-
-subroutine czspline_setup1(handle, n1, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1
-  real(fp), intent(in) :: f(n1)
-  integer, intent(out) :: ier
-  type(czspline1) :: self
-  self = transfer(handle, self)
-  call ezspline_setup(self % ptr, f, ier)
-end subroutine czspline_setup1
-
-subroutine czspline_setup2(handle, n1, n2, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1, n2
-  real(fp), intent(in) :: f(n1, n2)
-  integer, intent(out) :: ier
-  type(czspline2) :: self
-  self = transfer(handle, self)
-  call ezspline_setup(self % ptr, f, ier)
-end subroutine czspline_setup2
-
-subroutine czspline_setup3(handle, n1, n2, n3, f, ier)
-  use precision_mod, only: fp
-  use ezspline_obj
-  use ezspline
-  use czspline_pointer_types
-  implicit none
-  integer, intent(inout) :: handle(_ARRSZ)
-  integer, intent(in) :: n1, n2, n3
-  real(fp), intent(in) :: f(n1, n2, n3)
-  integer, intent(out) :: ier
-  type(czspline3) :: self
-  self = transfer(handle, self)
-  call ezspline_setup(self % ptr, f, ier)
-end subroutine czspline_setup3
diff --git a/czspline/czspline_test.cxx b/czspline/czspline_test.cxx
deleted file mode 100644
index 88378cc..0000000
--- a/czspline/czspline_test.cxx
+++ /dev/null
@@ -1,392 +0,0 @@
-/**
- *---------------------------------------------------------------------------
- * This code was developed at Tech-X (www.txcorp.com). It is free for any one
- * to use but comes with no warranty whatsoever. Use at your own risk. 
- * Thanks for reporting bugs to pletzer@txcorp.com. 
- *---------------------------------------------------------------------------
- * 
- * Test C++ bindings to ezspline
- */
-#include "czspline_capi.h"
-
-// std includes
-#include <iostream>
-#include <vector>
-#include <cmath>
-#include <cassert>
-
-/** 
- * Test 1d interpolation
- */
-void test1d() {
-
-  // opaque handle
-  int handle[_ARRSZ];
- 
-  // error code
-  int ier = 0;
-  // boundary condition type
-  int bcs1[2]; 
-  // not-a-knot
-  bcs1[0] = 0; bcs1[1] = 0;
-  
-  // create axes
-  int n1 = 11;
-  double dx = 1.0 / double(n1 - 1);
-  std::vector<double> x1(n1);
-  for (int i = 0; i < n1; ++i) 
-    x1[i] = i * dx;
-  // constructor
-  F77NAME(czspline_init1)(handle, &n1, bcs1, &ier);
-  assert(ier == 0);
-  F77NAME(czspline_set_axes1)(handle, &n1, &x1[0], &ier);
-  assert(ier == 0);
-
-  // set original data
-  std::vector<double> f(n1); 
-  for (int i = 0; i < n1; ++i) 
-    f[i] = x1[i]*x1[i]*x1[i];
-  F77NAME(czspline_setup1)(handle, &n1, &f[0], &ier);
-  assert(ier == 0);
-  
-  // target axes
-  int m1 = 101;
-  double dy = 1.0 / double(m1 - 1);
-  std::vector<double> y1(m1);
-  for (int i = 0; i < m1; ++i) 
-    y1[i] = i * dy;
-
-  // interpolated values
-  std::vector<double> g(m1);
-  
-  // point interpolation
-  int m = m1/2;
-  double y = y1[m];
-  F77NAME(czspline_interp1)(handle, &y, &g[m], &ier);
-  assert(ier == 0);
-  double error_point = std::abs(g[m] - y*y*y);
-
-  // cloud interpolation
-  F77NAME(czspline_interp1_cloud)(handle, &m1, &y1[0], &g[0], &ier);
-  assert(ier == 0);
-  double error_cloud = 0.0;
-  for (int i = 0; i < m1; ++i)
-    error_cloud += std::abs(g[i] - y1[i]*y1[i]*y1[i]);
-  error_cloud /= double(m1);
-
-  // array interpolation
-  F77NAME(czspline_interp1_array)(handle, &m1, &y1[0], &g[0], &ier);
-  assert(ier == 0);
-  double error_array = 0.0;
-  for (int i = 0; i < m1; ++i)
-    error_array += std::abs(g[i] - y1[i]*y1[i]*y1[i]);
-  error_array /= double(m1);
-
-  // save to file
-  std::string fname = "test1d.nc";
-  F77NAME(czspline_save1)(handle, fname.c_str(), &ier, fname.size());
-  assert(ier == 0);
-
-  // load from file
-  int handle2[_ARRSZ];
-  F77NAME(czspline_load1)(handle2, fname.c_str(), &ier, fname.size());
-  assert(ier == 0);
-
-  // clean up
-  F77NAME(czspline_free1)(handle, &ier);
-  assert(ier == 0);
-  F77NAME(czspline_free1)(handle2, &ier);
-  assert(ier == 0);
-
-  std::cout << "test1d interpolation errors\n";
-  std::cout << "point: " << error_point << '\n';
-  std::cout << "cloud: " << error_cloud << '\n';
-  std::cout << "array: " << error_array << '\n';
-}
-
-/** 
- * Test 2d interpolation
- */
-void test2d() {
-
-  // opaque handle
-  int handle[_ARRSZ];
-
-  // error code
-  int ier = 0;
-  // boundary condition type
-  int bcs1[2]; 
-  // not-a-knot
-  bcs1[0] = 0; bcs1[1] = 0;
-  int bcs2[2]; 
-  // periodic
-  bcs2[0] = -1; bcs2[1] = -1;
-  
-  // create axes
-  int n1 = 11;
-  int n2 = 12;
-  double pi = 3.1415926535897931;
-  double dx1 = 1.0 / double(n1 - 1);
-  double dx2 = 2.0 * pi / double(n2 - 1);
-  std::vector<double> x1(n1);
-  for (int i = 0; i < n1; ++i) 
-    x1[i] = i * dx1;
-  std::vector<double> x2(n2);
-  for (int i = 0; i < n2; ++i) 
-    x2[i] = i * dx2;
-  // constructor
-  F77NAME(czspline_init2)(handle, &n1, &n2, bcs1, bcs2, &ier);
-  assert(ier == 0);
-  F77NAME(czspline_set_axes2)(handle, &n1, &n2, &x1[0], &x2[0], &ier);
-  assert(ier == 0);
-
-  // set original data
-  std::vector<double> f(n1*n2);
-  for (int i2 = 0; i2 < n2; ++i2) {
-    for (int i1 = 0; i1 < n1; ++i1) {
-      int i = i1 + n1 * i2;
-      f[i] = x1[i1]*x1[i1]*x1[i1] * cos(x2[i2]);
-    }
-  }
-  F77NAME(czspline_setup2)(handle, &n1, &n2, &f[0], &ier);
-  assert(ier == 0);
-  
-  // target axes
-  int m1 = 101;
-  int m2 = 102;
-  double dy1 = 1.0 / double(m1 - 1);
-  double dy2 = 2.0 * pi / double(m2 - 1);
-  std::vector<double> y1(m1);
-  for (int i = 0; i < m1; ++i) 
-    y1[i] = i * dy1;
-  std::vector<double> y2(m2);
-  for (int i = 0; i < m2; ++i) 
-    y2[i] = i * dy2;
-
-  // interpolated values
-  std::vector<double> g(m1 * m2);
-  
-  // point interpolation
-  int i1p = m1/2;
-  int i2p = m2/2;
-  double y1p = y1[i1p];
-  double y2p = y2[i2p];
-  int i = i1p + m1 * i2p;
-  F77NAME(czspline_interp2)(handle, &y1p, &y2p, &g[i], &ier);
-  assert(ier == 0);
-  double error_point = std::abs(g[i] - y1p*y1p*y1p * cos(y2p));
-
-  // cloud interpolation (m1 points)
-  i1p = 0;
-  i = i1p + i2p * m1;
-  std::vector<double> y2pp(m1, y2p);
-  F77NAME(czspline_interp2_cloud)(handle, &m1, &y1[i1p], &y2pp[0], 
-				  &g[i], &ier);
-  assert(ier == 0);
-  double error_cloud = 0.0;
-  for (int i1 = 0; i1 < m1; ++i1)
-    error_cloud += std::abs(g[i + i1] - y1[i1]*y1[i1]*y1[i1]*cos(y2[i2p]));
-  error_cloud /= double(m1);
-
-  // array interpolation
-  F77NAME(czspline_interp2_array)(handle, &m1, &m2, 
-				  &y1[0], &y2[0], &g[0], &ier);
-  assert(ier == 0);
-  double error_array = 0.0;
-  for (int i2 = 0; i2 < m2; ++i2) {
-    for (int i1 = 0; i1 < m1; ++i1) {
-      int i = i1 + i2 * m1;
-      error_array += std::abs(g[i] - y1[i1]*y1[i1]*y1[i1] * cos(y2[i2]));
-    }
-  }
-  error_array /= double(m1 * m2);
-
-  // save to file
-  std::string fname = "test2d.nc";
-  F77NAME(czspline_save2)(handle, fname.c_str(), &ier, fname.size());
-  assert(ier == 0);
-
-  // load from file
-  int handle2[_ARRSZ];
-  F77NAME(czspline_load2)(handle2, fname.c_str(), &ier, fname.size());
-  assert(ier == 0);
-
-  // clean up
-  F77NAME(czspline_free2)(handle, &ier);
-  assert(ier == 0);
-  F77NAME(czspline_free2)(handle2, &ier);
-  assert(ier == 0);
-
-  std::cout << "test2d interpolation errors\n";
-  std::cout << "point: " << error_point << '\n';
-  std::cout << "cloud: " << error_cloud << '\n';
-  std::cout << "array: " << error_array << '\n';
-}
-
-/** 
- * Test 3d interpolation
- */
-void test3d() {
-
-  // opaque handle
-  int handle[_ARRSZ];
-
-  // error code
-  int ier = 0;
-  // boundary condition type
-  int bcs1[2]; 
-  // not-a-knot
-  bcs1[0] = 0; bcs1[1] = 0;
-  int bcs2[2]; 
-  // periodic
-  bcs2[0] = -1; bcs2[1] = -1;
-  int bcs3[2];
-  // slope and 2nd derivative
-  bcs3[0] = +1; bcs3[1] = +2;
-
-  // create axes
-  int n1 = 11;
-  int n2 = 12;
-  int n3 = 13;
-  double pi = 3.1415926535897931;
-  double dx1 = 1.0 / double(n1 - 1);
-  double dx2 = 2.0 * pi / double(n2 - 1);
-  double dx3 = 2.0 * pi / double(n3 - 1);
-  std::vector<double> x1(n1);
-  for (int i = 0; i < n1; ++i) 
-    x1[i] = i * dx1;
-  std::vector<double> x2(n2);
-  for (int i = 0; i < n2; ++i) 
-    x2[i] = i * dx2;
-  std::vector<double> x3(n3);
-  for (int i = 0; i < n3; ++i) 
-    x3[i] = i * dx3;
-  // constructor
-  F77NAME(czspline_init3)(handle, &n1, &n2, &n3, 
-			  bcs1, bcs2, bcs3, &ier);
-  assert(ier == 0);
-  // boundary conditions (values for periodic and not-a-knot
-  // are not used but must still be passed to the setter)
-  double bcval1[2];
-  double bcval2[2];
-  double bcval3[2];
-  bcval3[0] = 1.0; // df/dx3 @ x3=0
-  bcval3[1] = 0.0; // d^2/dx3^2 @ x3=2*pi
-  F77NAME(czspline_set_bcvals3)(handle, bcval1, bcval2, bcval3, &ier);
-  assert(ier == 0);
-  // set axes (necessary unless (0,..1) for anything but periodic BCs
-  // or (0..2*pi) for periodic BCs
-  F77NAME(czspline_set_axes3)(handle, &n1, &n2, &n3, 
-			      &x1[0], &x2[0], &x3[0], &ier);
-  assert(ier == 0);
-
-  // set original data
-  std::vector<double> f(n1*n2*n3);
-  for (int i3 = 0; i3 < n3; ++i3) {
-    for (int i2 = 0; i2 < n2; ++i2) {
-      for (int i1 = 0; i1 < n1; ++i1) {
-	int i = i1 + n1*(i2 + n2*i3);
-	f[i] = x1[i1]*x1[i1]*x1[i1] * cos(x2[i2]) * sin(x3[i3]);
-      }
-    }
-  }
-  F77NAME(czspline_setup3)(handle, &n1, &n2, &n3, 
-			   &f[0], &ier);
-  assert(ier == 0);
-  
-  // target axes
-  int m1 = 101;
-  int m2 = 102;
-  int m3 = 103;
-  double dy1 = 1.0 / double(m1 - 1);
-  double dy2 = 2.0 * pi / double(m2 - 1);
-  double dy3 = 2.0 * pi / double(m3 - 1);
-  std::vector<double> y1(m1);
-  for (int i = 0; i < m1; ++i) 
-    y1[i] = i * dy1;
-  std::vector<double> y2(m2);
-  for (int i = 0; i < m2; ++i) 
-    y2[i] = i * dy2;
-  std::vector<double> y3(m3);
-  for (int i = 0; i < m3; ++i) 
-    y3[i] = i * dy3;
-
-  // interpolated values
-  std::vector<double> g(m1 * m2 * m3);
-  
-  // point interpolation
-  int i1p = m1/2;
-  int i2p = m2/2;
-  int i3p = m3/2;
-  double y1p = y1[i1p];
-  double y2p = y2[i2p];
-  double y3p = y3[i3p];
-  int i = i1p + m1*(i2p + m2*i3p);
-  F77NAME(czspline_interp3)(handle, &y1p, &y2p, &y3p, 
-			    &g[i], &ier);
-  assert(ier == 0);
-  double error_point = std::abs(g[i] - y1p*y1p*y1p * cos(y2p) * sin(y3p));
-
-  // cloud interpolation (m1 points along x1 axis)
-  i1p = 0;
-  i = i1p + m1*(i2p + m2*i3p);
-  std::vector<double> y2pp(m1, y2p);
-  std::vector<double> y3pp(m1, y3p);
-  F77NAME(czspline_interp3_cloud)(handle, &m1, 
-				  &y1[i1p], &y2pp[0], &y3pp[0],
-				  &g[i], &ier);
-  assert(ier == 0);
-  double error_cloud = 0.0;
-  for (int i1 = 0; i1 < m1; ++i1)
-    error_cloud += std::abs(g[i + i1] - 
-			    y1[i1]*y1[i1]*y1[i1]*cos(y2[i2p])*sin(y3[i3p]));
-  error_cloud /= double(m1);
-
-  // array interpolation
-  F77NAME(czspline_interp3_array)(handle, &m1, &m2, &m3,
-				  &y1[0], &y2[0], &y3[0],
-				  &g[0], &ier);
-  assert(ier == 0);
-  double error_array = 0.0;
-  for (int i3 = 0; i3 < m3; ++i3) {
-    for (int i2 = 0; i2 < m2; ++i2) {
-      for (int i1 = 0; i1 < m1; ++i1) {
-	int i = i1 + m1*(i2 + m2*i3);
-	error_array += std::abs(g[i] - 
-				y1[i1]*y1[i1]*y1[i1] * cos(y2[i2]) * sin(y3[i3]));
-      }
-    }
-  }
-  error_array /= double(m1 * m2 * m3);
-
-  // save to file
-  std::string fname = "test3d.nc";
-  F77NAME(czspline_save3)(handle, fname.c_str(), &ier, fname.size());
-  assert(ier == 0);
-
-  // load from file
-  int handle2[_ARRSZ];
-  F77NAME(czspline_load3)(handle2, fname.c_str(), &ier, fname.size());
-  assert(ier == 0);
-
-  // clean up
-  F77NAME(czspline_free3)(handle, &ier);
-  assert(ier == 0);
-  F77NAME(czspline_free3)(handle2, &ier);
-  assert(ier == 0);
-
-  std::cout << "test3d interpolation errors\n";
-  std::cout << "point: " << error_point << '\n';
-  std::cout << "cloud: " << error_cloud << '\n';
-  std::cout << "array: " << error_array << '\n';
-}
-
-/**
- * Main driver
- */
-int main() {
-    test1d();
-    test2d();
-    test3d();
-    return 0;
-}
diff --git a/czspline/f77name.h b/czspline/f77name.h
deleted file mode 100644
index fbd619c..0000000
--- a/czspline/f77name.h
+++ /dev/null
@@ -1,18 +0,0 @@
-#ifndef F77NAME_H  /*  F77NAME(c_name) make c_name fortran-callable */
-#define F77NAME_H  /*  dmc 9 Apr 1999 */
-
-#if __VMS || __IBM || __RS6000 || __HP || __ABS || defined(__xlC__)
-#define F77NAME(name) name
-#endif
-
-#if __CRAY
-#define F77NAME(name) name
-#endif
-
-/* fallthru:  append the "_", which most systems want. */
-
-#ifndef F77NAME
-#define F77NAME(name) name##_
-#endif
-
-#endif
diff --git a/ezspline/ezspline_cinterp.f90 b/ezspline/ezspline_cinterp.f90
index 22d6cf0..6cd6814 100644
--- a/ezspline/ezspline_cinterp.f90
+++ b/ezspline/ezspline_cinterp.f90
@@ -51,201 +51,168 @@ subroutine EZspline_cinterp(spline_o, k, p1, p2, p3, f, ier)
   ier = 0
   ifail=0
   if(spline_o%isReady /= 1) then
-     ier = 94
-     return
+    ier = 94
+    return
   end if
 
   if (spline_o%isLinear == 1) then
 
-     call vecpc3(ict, k, p1, p2, p3, k, f, &
-          & spline_o%n1, spline_o%x1pkg(1,1), &
-          & spline_o%n2, spline_o%x2pkg(1,1), &
-          & spline_o%n3, spline_o%x3pkg(1,1), &
-          & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
-          & iwarn,ifail)
+    call vecpc3(ict, k, p1, p2, p3, k, f, &
+         spline_o%n1, spline_o%x1pkg(1,1), &
+         spline_o%n2, spline_o%x2pkg(1,1), &
+         spline_o%n3, spline_o%x3pkg(1,1), &
+         spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
+         iwarn,ifail)
 
   else if (spline_o%isHybrid == 1) then
 
-     call vecintrp3d(ict, k, p1, p2, p3, k, f, &
-          & spline_o%n1, spline_o%x1pkg(1,1), &
-          & spline_o%n2, spline_o%x2pkg(1,1), &
-          & spline_o%n3, spline_o%x3pkg(1,1), &
-          & spline_o%hspline, &
-          & spline_o%fspl(1,1,1,1), size(spline_o%fspl,1), &
-          & size(spline_o%fspl,2), size(spline_o%fspl,3), &
-          & size(spline_o%fspl,4), iwarn, ifail)
-
-  else if (spline_o%isHermite ==0) then
-     !
-!!$     call vectricub(ict, k, p1, p2, p3, k, f, &
-!!$          & spline_o%n1,spline_o%x1pkg(1,1), &
-!!$          & spline_o%n2,spline_o%x2pkg(1,1), &
-!!$          & spline_o%n3,spline_o%x3pkg(1,1), &
-!!$          & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
-!!$          & iwarn, ifail)
-
-     call xlookup(k,p1,spline_o%n1,spline_o%x1pkg(1,1),2,ix,dxn,hx,hxi,ifail)
-     call xlookup(k,p2,spline_o%n2,spline_o%x2pkg(1,1),2,iy,dyn,hy,hyi,ifail)
-     call xlookup(k,p3,spline_o%n3,spline_o%x3pkg(1,1),2,iz,dzn,hz,hzi,ifail)
-
-!!$      call fvtricub(ict,k,k,f,ix,iy,iz,dxn,dyn,dzn,           &
-!!$     &   hx,hxi,hy,hyi,hz,hzi,spline_o%fspl(1,1,1,1), &
-!!$     &   spline_o%n1, spline_o%n2, spline_o%n3)
-
-
-     z36th=sixth*sixth
-     z216th=sixth*sixth*sixth
-     !
-     !  prepare useful parameters...
-     !
-     do v=1,k
-        j1=ix(v)
-        j2=iy(v)
-        j3=iz(v)
+    call vecintrp3d(ict, k, p1, p2, p3, k, f, &
+         spline_o%n1, spline_o%x1pkg(1,1), &
+         spline_o%n2, spline_o%x2pkg(1,1), &
+         spline_o%n3, spline_o%x3pkg(1,1), &
+         spline_o%hspline, &
+         spline_o%fspl(1,1,1,1), size(spline_o%fspl,1), &
+         size(spline_o%fspl,2), size(spline_o%fspl,3), &
+         size(spline_o%fspl,4), iwarn, ifail)
+
+  else if (spline_o%isHermite == 0) then
+
+    call xlookup(k,p1,spline_o%n1,spline_o%x1pkg(1,1),2,ix,dxn,hx,hxi,ifail)
+    call xlookup(k,p2,spline_o%n2,spline_o%x2pkg(1,1),2,iy,dyn,hy,hyi,ifail)
+    call xlookup(k,p3,spline_o%n3,spline_o%x3pkg(1,1),2,iz,dzn,hz,hzi,ifail)
+
+    z36th=sixth*sixth
+    z216th=sixth*sixth*sixth
+    !
+    !  prepare useful parameters...
+    !
+    do v=1,k
+      j1=ix(v)
+      j2=iy(v)
+      j3=iz(v)
+      !
+      !   ...in x direction
+      !
+      xp=dxn(v)
+      xpi=one-xp
+      xp2=xp*xp
+      xpi2=xpi*xpi
+      !
+      cx=xp*(xp2-ONE)
+      cxi=xpi*(xpi2-ONE)
+      hx2=hx(v)*hx(v)
+      !
+      !   ...and in y direction
+      !
+      yp=dyn(v)
+      ypi=ONE-yp
+      yp2=yp*yp
+      ypi2=ypi*ypi
+      !
+      cy=yp*(yp2-ONE)
+      cyi=ypi*(ypi2-ONE)
+      hy2=hy(v)*hy(v)
+      !
+      !   ...and in z direction
+      !
+      zp=dzn(v)
+      zpi=ONE-zp
+      zp2=zp*zp
+      zpi2=zpi*zpi
+      !
+      cz=zp*(zp2-ONE)
+      czi=zpi*(zpi2-ONE)
+      hz2=hz(v)*hz(v)
+      !
+      !  get desired values:
+      !
+      if(ict(1).eq.1) then
         !
-        !   ...in x direction
+        !  function value:
         !
-        xp=dxn(v)
-        xpi=one-xp
-        xp2=xp*xp
-        xpi2=xpi*xpi
+        somme=(                                                       &
+             zpi*(                                                    &
+             xpi*(ypi*spline_o%fspl(1,j1,j2,j3)  +yp*spline_o%fspl(1,j1,j2+1,j3))+            &
+             xp*(ypi*spline_o%fspl(1,j1+1,j2,j3)+yp*spline_o%fspl(1,j1+1,j2+1,j3)))          &
+             +zp*(                                                    &
+             xpi*(ypi*spline_o%fspl(1,j1,j2,j3+1)  +yp*spline_o%fspl(1,j1,j2+1,j3+1))+        &
+             xp*(ypi*spline_o%fspl(1,j1+1,j2,j3+1)+yp*spline_o%fspl(1,j1+1,j2+1,j3+1))))
         !
-!!$         if((ict(1).eq.1).or.(ict(3).eq.1).or.(ict(4).eq.1).or.         &
-!!$     &      (ict(6).eq.1).or.(ict(7).eq.1).or.(ict(10).eq.1)) then
-        cx=xp*(xp2-ONE)
-        cxi=xpi*(xpi2-ONE)
-        hx2=hx(v)*hx(v)
-!!$         end if
-!!$         if((ict(2).eq.1).or.(ict(8).eq.1).or.(ict(9).eq.1)) then
-!!$            cxd=THREE*xp2-ONE
-!!$            cxdi=-THREE*xpi2+ONE
-!!$         end if
+        somme=somme+sixth*hx2*(                                       &
+             zpi*(                                                    &
+             cxi*(ypi*spline_o%fspl(2,j1,j2,j3)  +yp*spline_o%fspl(2,j1,j2+1,j3))+            &
+             cx*(ypi*spline_o%fspl(2,j1+1,j2,j3)+yp*spline_o%fspl(2,j1+1,j2+1,j3)))          &
+             +zp*(                                                    &
+             cxi*(ypi*spline_o%fspl(2,j1,j2,j3+1)  +yp*spline_o%fspl(2,j1,j2+1,j3+1))+        &
+             cx*(ypi*spline_o%fspl(2,j1+1,j2,j3+1)+yp*spline_o%fspl(2,j1+1,j2+1,j3+1))))
         !
-        !   ...and in y direction
+        somme=somme+sixth*hy2*(                                       &
+             zpi*(                                                    &
+             xpi*(cyi*spline_o%fspl(3,j1,j2,j3)  +cy*spline_o%fspl(3,j1,j2+1,j3))+            &
+             xp*(cyi*spline_o%fspl(3,j1+1,j2,j3)+cy*spline_o%fspl(3,j1+1,j2+1,j3)))          &
+             +zp*(                                                    &
+             xpi*(cyi*spline_o%fspl(3,j1,j2,j3+1)  +cy*spline_o%fspl(3,j1,j2+1,j3+1))+        &
+             xp*(cyi*spline_o%fspl(3,j1+1,j2,j3+1)+cy*spline_o%fspl(3,j1+1,j2+1,j3+1))))
         !
-        yp=dyn(v)
-        ypi=ONE-yp
-        yp2=yp*yp
-        ypi2=ypi*ypi
+        somme=somme+sixth*hz2*(                                       &
+             czi*(                                                    &
+             xpi*(ypi*spline_o%fspl(4,j1,j2,j3)  +yp*spline_o%fspl(4,j1,j2+1,j3))+            &
+             xp*(ypi*spline_o%fspl(4,j1+1,j2,j3)+yp*spline_o%fspl(4,j1+1,j2+1,j3)))          &
+             +cz*(                                                    &
+             xpi*(ypi*spline_o%fspl(4,j1,j2,j3+1)  +yp*spline_o%fspl(4,j1,j2+1,j3+1))+        &
+             xp*(ypi*spline_o%fspl(4,j1+1,j2,j3+1)+yp*spline_o%fspl(4,j1+1,j2+1,j3+1))))
         !
-!!$         if((ict(1).eq.1).or.(ict(2).eq.1).or.(ict(4).eq.1).or.         &
-!!$     &      (ict(5).eq.1).or.(ict(7).eq.1).or.(ict(9).eq.1)) then
-        cy=yp*(yp2-ONE)
-        cyi=ypi*(ypi2-ONE)
-        hy2=hy(v)*hy(v)
-!!$         end if
-!!$         if((ict(3).eq.1).or.(ict(8).eq.1).or.(ict(10).eq.1)) then
-!!$            cyd=THREE*yp2-ONE
-!!$            cydi=-THREE*ypi2+ONE
-!!$         end if
-        !
-        !   ...and in z direction
-        !
-        zp=dzn(v)
-        zpi=ONE-zp
-        zp2=zp*zp
-        zpi2=zpi*zpi
-        !
-!!$         if((ict(1).eq.1).or.(ict(2).eq.1).or.(ict(3).eq.1).or.         &
-!!$     &      (ict(5).eq.1).or.(ict(6).eq.1).or.(ict(8).eq.1)) then
-        cz=zp*(zp2-ONE)
-        czi=zpi*(zpi2-ONE)
-        hz2=hz(v)*hz(v)
-!!$         end if
-!!$         if((ict(4).eq.1).or.(ict(9).eq.1).or.(ict(10).eq.1)) then
-!!$            czd=THREE*zp2-ONE
-!!$            czdi=-THREE*zpi2+ONE
-!!$         end if
-        !
-        !  get desired values:
-        !
-        if(ict(1).eq.1) then
-           !
-           !  function value:
-           !
-           somme=(                                                       &
-                &         zpi*(                                                    &
-                &           xpi*(ypi*spline_o%fspl(1,j1,j2,j3)  +yp*spline_o%fspl(1,j1,j2+1,j3))+            &
-                &            xp*(ypi*spline_o%fspl(1,j1+1,j2,j3)+yp*spline_o%fspl(1,j1+1,j2+1,j3)))          &
-                &        +zp*(                                                     &
-                &           xpi*(ypi*spline_o%fspl(1,j1,j2,j3+1)  +yp*spline_o%fspl(1,j1,j2+1,j3+1))+        &
-                &            xp*(ypi*spline_o%fspl(1,j1+1,j2,j3+1)+yp*spline_o%fspl(1,j1+1,j2+1,j3+1))))
-           !
-           somme=somme+sixth*hx2*(                                         &
-                &         zpi*(                                                    &
-                &           cxi*(ypi*spline_o%fspl(2,j1,j2,j3)  +yp*spline_o%fspl(2,j1,j2+1,j3))+            &
-                &            cx*(ypi*spline_o%fspl(2,j1+1,j2,j3)+yp*spline_o%fspl(2,j1+1,j2+1,j3)))          &
-                &        +zp*(                                                     &
-                &           cxi*(ypi*spline_o%fspl(2,j1,j2,j3+1)  +yp*spline_o%fspl(2,j1,j2+1,j3+1))+        &
-                &            cx*(ypi*spline_o%fspl(2,j1+1,j2,j3+1)+yp*spline_o%fspl(2,j1+1,j2+1,j3+1))))
+        somme=somme+z36th*hx2*hy2*(                                   &
+             zpi*(                                                    &
+             cxi*(cyi*spline_o%fspl(5,j1,j2,j3)  +cy*spline_o%fspl(5,j1,j2+1,j3))+            &
+             cx*(cyi*spline_o%fspl(5,j1+1,j2,j3)+cy*spline_o%fspl(5,j1+1,j2+1,j3)))          &
+             +zp*(                                                    &
+             cxi*(cyi*spline_o%fspl(5,j1,j2,j3+1)  +cy*spline_o%fspl(5,j1,j2+1,j3+1))+        &
+             cx*(cyi*spline_o%fspl(5,j1+1,j2,j3+1)+cy*spline_o%fspl(5,j1+1,j2+1,j3+1))))
            !
-           somme=somme+sixth*hy2*(                                         &
-                &         zpi*(                                                    &
-                &           xpi*(cyi*spline_o%fspl(3,j1,j2,j3)  +cy*spline_o%fspl(3,j1,j2+1,j3))+            &
-                &            xp*(cyi*spline_o%fspl(3,j1+1,j2,j3)+cy*spline_o%fspl(3,j1+1,j2+1,j3)))          &
-                &        +zp*(                                                     &
-                &           xpi*(cyi*spline_o%fspl(3,j1,j2,j3+1)  +cy*spline_o%fspl(3,j1,j2+1,j3+1))+        &
-                &            xp*(cyi*spline_o%fspl(3,j1+1,j2,j3+1)+cy*spline_o%fspl(3,j1+1,j2+1,j3+1))))
-           !
-           somme=somme+sixth*hz2*(                                         &
-                &         czi*(                                                    &
-                &           xpi*(ypi*spline_o%fspl(4,j1,j2,j3)  +yp*spline_o%fspl(4,j1,j2+1,j3))+            &
-                &            xp*(ypi*spline_o%fspl(4,j1+1,j2,j3)+yp*spline_o%fspl(4,j1+1,j2+1,j3)))          &
-                &        +cz*(                                                     &
-                &           xpi*(ypi*spline_o%fspl(4,j1,j2,j3+1)  +yp*spline_o%fspl(4,j1,j2+1,j3+1))+        &
-                &            xp*(ypi*spline_o%fspl(4,j1+1,j2,j3+1)+yp*spline_o%fspl(4,j1+1,j2+1,j3+1))))
-           !
-           somme=somme+z36th*hx2*hy2*(                                     &
-                &         zpi*(                                                    &
-                &           cxi*(cyi*spline_o%fspl(5,j1,j2,j3)  +cy*spline_o%fspl(5,j1,j2+1,j3))+            &
-                &            cx*(cyi*spline_o%fspl(5,j1+1,j2,j3)+cy*spline_o%fspl(5,j1+1,j2+1,j3)))          &
-                &        +zp*(                                                     &
-                &           cxi*(cyi*spline_o%fspl(5,j1,j2,j3+1)  +cy*spline_o%fspl(5,j1,j2+1,j3+1))+        &
-                &            cx*(cyi*spline_o%fspl(5,j1+1,j2,j3+1)+cy*spline_o%fspl(5,j1+1,j2+1,j3+1))))
-           !
-           somme=somme+z36th*hx2*hz2*(                                     &
-                &         czi*(                                                    &
-                &           cxi*(ypi*spline_o%fspl(6,j1,j2,j3)  +yp*spline_o%fspl(6,j1,j2+1,j3))+            &
-                &            cx*(ypi*spline_o%fspl(6,j1+1,j2,j3)+yp*spline_o%fspl(6,j1+1,j2+1,j3)))          &
-                &        +cz*(                                                     &
-                &           cxi*(ypi*spline_o%fspl(6,j1,j2,j3+1)  +yp*spline_o%fspl(6,j1,j2+1,j3+1))+        &
-                &            cx*(ypi*spline_o%fspl(6,j1+1,j2,j3+1)+yp*spline_o%fspl(6,j1+1,j2+1,j3+1))))
-           !
-           somme=somme+z36th*hy2*hz2*(                                     &
-                &         czi*(                                                    &
-                &           xpi*(cyi*spline_o%fspl(7,j1,j2,j3)  +cy*spline_o%fspl(7,j1,j2+1,j3))+            &
-                &            xp*(cyi*spline_o%fspl(7,j1+1,j2,j3)+cy*spline_o%fspl(7,j1+1,j2+1,j3)))          &
-                &        +cz*(                                                     &
-                &           xpi*(cyi*spline_o%fspl(7,j1,j2,j3+1)  +cy*spline_o%fspl(7,j1,j2+1,j3+1))+        &
-                &            xp*(cyi*spline_o%fspl(7,j1+1,j2,j3+1)+cy*spline_o%fspl(7,j1+1,j2+1,j3+1))))
-           !
-           somme=somme+z216th*hx2*hy2*hz2*(                                &
-                &         czi*(                                                    &
-                &           cxi*(cyi*spline_o%fspl(8,j1,j2,j3)  +cy*spline_o%fspl(8,j1,j2+1,j3))+            &
-                &            cx*(cyi*spline_o%fspl(8,j1+1,j2,j3)+cy*spline_o%fspl(8,j1+1,j2+1,j3)))          &
-                &        +cz*(                                                     &
-                &           cxi*(cyi*spline_o%fspl(8,j1,j2,j3+1)  +cy*spline_o%fspl(8,j1,j2+1,j3+1))+        &
-                &            cx*(cyi*spline_o%fspl(8,j1+1,j2,j3+1)+cy*spline_o%fspl(8,j1+1,j2+1,j3+1))))
-           !
-           f(v)=somme
-        else
-           write(6,*) 'ERROR ict(1) must be 1 (interpolation)'
-        end if
+        somme=somme+z36th*hx2*hz2*(                                   &
+             czi*(                                                    &
+             cxi*(ypi*spline_o%fspl(6,j1,j2,j3)  +yp*spline_o%fspl(6,j1,j2+1,j3))+            &
+             cx*(ypi*spline_o%fspl(6,j1+1,j2,j3)+yp*spline_o%fspl(6,j1+1,j2+1,j3)))          &
+             +cz*(                                                    &
+             cxi*(ypi*spline_o%fspl(6,j1,j2,j3+1)  +yp*spline_o%fspl(6,j1,j2+1,j3+1))+        &
+             cx*(ypi*spline_o%fspl(6,j1+1,j2,j3+1)+yp*spline_o%fspl(6,j1+1,j2+1,j3+1))))
+        !
+        somme=somme+z36th*hy2*hz2*(                                   &
+             czi*(                                                    &
+             xpi*(cyi*spline_o%fspl(7,j1,j2,j3)  +cy*spline_o%fspl(7,j1,j2+1,j3))+            &
+             xp*(cyi*spline_o%fspl(7,j1+1,j2,j3)+cy*spline_o%fspl(7,j1+1,j2+1,j3)))          &
+             +cz*(                                                    &
+             xpi*(cyi*spline_o%fspl(7,j1,j2,j3+1)  +cy*spline_o%fspl(7,j1,j2+1,j3+1))+        &
+             xp*(cyi*spline_o%fspl(7,j1+1,j2,j3+1)+cy*spline_o%fspl(7,j1+1,j2+1,j3+1))))
         !
+        somme=somme+z216th*hx2*hy2*hz2*(                              &
+             czi*(                                                    &
+             cxi*(cyi*spline_o%fspl(8,j1,j2,j3)  +cy*spline_o%fspl(8,j1,j2+1,j3))+            &
+             cx*(cyi*spline_o%fspl(8,j1+1,j2,j3)+cy*spline_o%fspl(8,j1+1,j2+1,j3)))          &
+             +cz*(                                                    &
+             cxi*(cyi*spline_o%fspl(8,j1,j2,j3+1)  +cy*spline_o%fspl(8,j1,j2+1,j3+1))+        &
+             cx*(cyi*spline_o%fspl(8,j1+1,j2,j3+1)+cy*spline_o%fspl(8,j1+1,j2+1,j3+1))))
         !
-     end do                             ! vector loop
+        f(v)=somme
+      else
+        write(6,*) 'ERROR ict(1) must be 1 (interpolation)'
+      end if
+      !
+      !
+    end do                             ! vector loop
 
 
-     if(ifail /= 0) ier = 97
+    if(ifail /= 0) ier = 97
 
   else
 
-     call vecherm3(ict, k, p1, p2, p3, k, f, &
-          & spline_o%n1, spline_o%x1pkg(1,1), &
-          & spline_o%n2, spline_o%x2pkg(1,1), &
-          & spline_o%n3, spline_o%x3pkg(1,1), &
-          & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
-          & iwarn,ifail)
+    call vecherm3(ict, k, p1, p2, p3, k, f, &
+         spline_o%n1, spline_o%x1pkg(1,1), &
+         spline_o%n2, spline_o%x2pkg(1,1), &
+         spline_o%n3, spline_o%x3pkg(1,1), &
+         spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
+         iwarn,ifail)
 
   end if
 
diff --git a/ezspline/ezspline_init.f90 b/ezspline/ezspline_init.f90
index 8d0d39a..647929b 100644
--- a/ezspline/ezspline_init.f90
+++ b/ezspline/ezspline_init.f90
@@ -18,10 +18,10 @@ subroutine EZspline_init1(spline_o, n1, BCS1, ier)
   ier = 0
 
   if(EZspline_allocated(spline_o)) then
-     ier = 100  ! allocated already
-     return
+    ier = 100  ! allocated already
+    return
   else
-     call EZspline_preInit(spline_o)
+    call EZspline_preInit(spline_o)
   end if
 
   spline_o%n1 = n1
@@ -40,27 +40,26 @@ subroutine EZspline_init1(spline_o, n1, BCS1, ier)
   if(iok /= 0) ier = 1
 
   do i = 1, 2
-
-     spline_o%bcval1min = 0.0_fp
-     spline_o%bcval1max = 0.0_fp
-     select case(BCS1(i))
-     case (-1)
-        spline_o%ibctype1(i) = -1
-     case (0)
-        spline_o%ibctype1(i) = 0
-     case (1)
-        spline_o%ibctype1(i) = 1
-     case (2)
-        spline_o%ibctype1(i) = 2
-     case default
-        ier = 2
-        spline_o%ibctype1(i) = 0
-     end select
-
+    spline_o%bcval1min = 0.0_fp
+    spline_o%bcval1max = 0.0_fp
+    select case(BCS1(i))
+    case (-1)
+      spline_o%ibctype1(i) = -1
+    case (0)
+      spline_o%ibctype1(i) = 0
+    case (1)
+      spline_o%ibctype1(i) = 1
+    case (2)
+      spline_o%ibctype1(i) = 2
+    case default
+      ier = 2
+      spline_o%ibctype1(i) = 0
+    end select
   end do
+
   if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) then
-     spline_o%ibctype1(1)=-1
-     spline_o%ibctype1(2)=-1
+    spline_o%ibctype1(1)=-1
+    spline_o%ibctype1(2)=-1
   end if
 
   !
@@ -70,7 +69,7 @@ subroutine EZspline_init1(spline_o, n1, BCS1, ier)
   if(BCS1(2)==-1) spline_o%x1max = ezspline_twopi
 
   spline_o%x1 = spline_o%x1min + (spline_o%x1max - spline_o%x1min)* &
-       & (/ (real(i-1,fp)/real(spline_o%n1-1, fp), i=1,spline_o%n1 ) /)
+       (/ (real(i-1,fp)/real(spline_o%n1-1, fp), i=1,spline_o%n1 ) /)
 
   spline_o%isReady = 0
 
@@ -78,8 +77,6 @@ subroutine EZspline_init1(spline_o, n1, BCS1, ier)
 end subroutine EZspline_init1
 
 
-
-
 subroutine EZspline_init2(spline_o, n1, n2, BCS1, BCS2, ier)
   use precision_mod, only: fp
   use EZspline_obj
@@ -188,12 +185,12 @@ subroutine EZspline_init2(spline_o, n1, n2, BCS1, BCS2, ier)
   if(BCS1(2)==-1) spline_o%x1max = ezspline_twopi
 
   spline_o%x1 = spline_o%x1min + (spline_o%x1max - spline_o%x1min)* &
-       & (/ (real(i-1,fp)/real(spline_o%n1-1, fp), i=1,spline_o%n1 ) /)
+       (/ (real(i-1,fp)/real(spline_o%n1-1, fp), i=1,spline_o%n1 ) /)
   spline_o%x2min = 0.0_fp
   spline_o%x2max = 1.0_fp
   if(BCS2(2)==-1) spline_o%x2max = ezspline_twopi
   spline_o%x2 = spline_o%x2min + (spline_o%x2max - spline_o%x2min)* &
-       & (/ (real(i-1,fp)/real(spline_o%n2-1, fp), i=1,spline_o%n2 ) /)
+       (/ (real(i-1,fp)/real(spline_o%n2-1, fp), i=1,spline_o%n2 ) /)
 
   spline_o%isReady = 0
 
@@ -342,18 +339,18 @@ subroutine EZspline_init3(spline_o, n1, n2, n3, BCS1, BCS2, BCS3, ier)
   if(BCS1(2)==-1) spline_o%x1max = ezspline_twopi
 
   spline_o%x1 = spline_o%x1min + (spline_o%x1max - spline_o%x1min)* &
-       & (/ (real(i-1,fp)/real(spline_o%n1-1, fp), i=1,spline_o%n1 ) /)
+       (/ (real(i-1,fp)/real(spline_o%n1-1, fp), i=1,spline_o%n1 ) /)
   spline_o%x2min = 0.0_fp
   spline_o%x2max = 1.0_fp
   if(BCS2(2)==-1) spline_o%x2max = ezspline_twopi
   spline_o%x2 = spline_o%x2min + (spline_o%x2max - spline_o%x2min)* &
-       & (/ (real(i-1,fp)/real(spline_o%n2-1, fp), i=1,spline_o%n2 ) /)
+       (/ (real(i-1,fp)/real(spline_o%n2-1, fp), i=1,spline_o%n2 ) /)
 
   spline_o%x3min = 0.0_fp
   spline_o%x3max = 1.0_fp
   if(BCS3(2)==-1) spline_o%x3max = ezspline_twopi
   spline_o%x3 = spline_o%x3min + (spline_o%x3max - spline_o%x3min)* &
-       & (/ (real(i-1,fp)/real(spline_o%n3-1, fp), i=1,spline_o%n3 ) /)
+       (/ (real(i-1,fp)/real(spline_o%n3-1, fp), i=1,spline_o%n3 ) /)
 
   spline_o%isReady = 0
 
diff --git a/ezspline/ezspline_load.f90 b/ezspline/ezspline_load.f90
index cbd069b..ec5a3a6 100644
--- a/ezspline/ezspline_load.f90
+++ b/ezspline/ezspline_load.f90
@@ -31,8 +31,6 @@ subroutine EZspline_load1(spline_o, filename, ier, spl_name)
 
   integer ncid, ifail
   integer dimlens(3)
-  character(len=2), parameter :: real8='R8'
-  character(len=3), parameter :: int='INT'
   integer n1, BCS1(2)
   real(fp), dimension(:), allocatable :: f
 
@@ -42,16 +40,16 @@ subroutine EZspline_load1(spline_o, filename, ier, spl_name)
   logical :: fullsv
 
   if(present(spl_name)) then
-     zpre = spl_name
+    zpre = spl_name
   else
-     zpre = ' '
+    zpre = ' '
   end if
 
   ier = 0
   call cdfOpn(ncid, filename, 'r') ! no error flag??
   if(ncid==0) then
-     ier=43
-     return
+    ier=43
+    return
   end if
 
   ! check if n1 is present; check if spline object has correct rank
@@ -164,8 +162,6 @@ subroutine EZspline_load2(spline_o, filename, ier, spl_name)
 
   integer ncid, ifail, in0, in1, in2
   integer dimlens(3)
-  character(len=2), parameter :: real8='R8'
-  character(len=3), parameter :: int='INT'
   integer n1, n2, BCS1(2), BCS2(2), hspline(2)
   real(fp), dimension(:,:), allocatable :: f
 
@@ -319,8 +315,6 @@ subroutine EZspline_load3(spline_o, filename, ier, spl_name)
 
   integer ncid, ifail, in0, in1, in2, in3
   integer dimlens(3)
-  character(len=2), parameter :: real8='R8'
-  character(len=3), parameter :: int='INT'
   integer n1, n2, n3, BCS1(2), BCS2(2), BCS3(2), hspline(3)
   real(fp), dimension(:,:,:), allocatable :: f
 
diff --git a/ezspline/ezspline_mod.f90 b/ezspline/ezspline_mod.f90
index d639c0a..02f26d3 100644
--- a/ezspline/ezspline_mod.f90
+++ b/ezspline/ezspline_mod.f90
@@ -3,9 +3,9 @@ module EZspline_type
 
   real(fp), parameter :: ezspline_twopi = 6.2831853071795865_fp
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! EZspline data types
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   type EZspline3
     !
@@ -828,7 +828,6 @@ end subroutine EZspline_setup1
   end interface
 
 
-
   interface EZspline_interp
     !
     ! Interpolation at grid point(s) p1, [p2, [p3]]. Result is returned in
@@ -1221,8 +1220,8 @@ end subroutine EZspline_save2
     subroutine EZspline_save1(spline_o,filename,ier,spl_name,fullsave)
       use EZspline_obj
       use EZcdf
-      type(EZspline1) :: spline_o
-      character(len=*) :: filename
+      type(EZspline1), intent(in) :: spline_o
+      character(len=*), intent(in) :: filename
       integer, intent(out) :: ier
 
       character(len=*), intent(in),optional :: spl_name
@@ -1253,8 +1252,7 @@ subroutine EZspline_load3(spline_o, filename, ier, spl_name)
       type(EZspline3) :: spline_o
       character(len=*) :: filename
       integer, intent(out) :: ier
-
-      character(len=*), intent(in),optional :: spl_name
+      character(len=*), intent(in), optional :: spl_name
     end subroutine EZspline_load3
 
     subroutine EZspline_load2(spline_o, filename, ier, spl_name)
@@ -1263,8 +1261,7 @@ subroutine EZspline_load2(spline_o, filename, ier, spl_name)
       type(EZspline2) :: spline_o
       character(len=*) :: filename
       integer, intent(out) :: ier
-
-      character(len=*), intent(in),optional :: spl_name
+      character(len=*), intent(in), optional :: spl_name
     end subroutine EZspline_load2
 
     subroutine EZspline_load1(spline_o, filename, ier, spl_name)
@@ -1273,8 +1270,7 @@ subroutine EZspline_load1(spline_o, filename, ier, spl_name)
       type(EZspline1) :: spline_o
       character(len=*) :: filename
       integer, intent(out) :: ier
-
-      character(len=*), intent(in),optional :: spl_name
+      character(len=*), intent(in), optional :: spl_name
     end subroutine EZspline_load1
 
   end interface
diff --git a/ezspline/ezspline_save.f90 b/ezspline/ezspline_save.f90
index 270ce7a..96165d6 100644
--- a/ezspline/ezspline_save.f90
+++ b/ezspline/ezspline_save.f90
@@ -11,13 +11,12 @@
 ! 1-D
 !
 
-subroutine EZspline_save1(spline_o, filename, ier, &
-     spl_name, fullsave)
+subroutine EZspline_save1(spline_o, filename, ier, spl_name, fullsave)
   use EZspline_obj
   use EZcdf
   implicit none
-  type(EZspline1) :: spline_o
-  character(len=*) :: filename
+  type(EZspline1),  intent(in) :: spline_o
+  character(len=*), intent(in) :: filename
   ! ier:
   ! 17=could not save spline object in file filename
   integer, intent(out) :: ier
@@ -25,11 +24,11 @@ subroutine EZspline_save1(spline_o, filename, ier, &
   !  if SPL_NAME is set, APPEND to file instead of creating new; prepend
   !  name to all NetCDF data items.  This allows one file to contain
   !  multiple items. (new dmc Mar 2006)
-  character(len=*), intent(in), OPTIONAL :: SPL_NAME  ! optional spline "name"
+  character(len=*), intent(in), optional :: spl_name  ! optional spline "name"
 
   !  if FULLSAVE is set .TRUE., save derived coefficients along with data
   !  this saves recalculating the coefficients when file is read (Mar 2006)
-  logical, intent(in), OPTIONAL :: FULLSAVE
+  logical, intent(in), optional :: fullsave
 
   integer ncid, ifail, in0, in1
   integer dimlens(3)
@@ -41,37 +40,36 @@ subroutine EZspline_save1(spline_o, filename, ier, &
 
   ier = 0
   if(spline_o%isReady /= 1) then
-     ier = 94
-     return
+    ier = 94
+    return
   end if
 
   in0 = size(spline_o%fspl,1)
   in1 = size(spline_o%fspl,2)
-
   if(present(spl_name)) then
-     call ezspline_spl_name_chk(spl_name,ier)
-     if(ier.ne.0) return
-     zpre=spl_name
+    call ezspline_spl_name_chk(spl_name,ier)
+    if(ier.ne.0) return
+    zpre=spl_name
   else
-     zpre=' '
+    zpre=' '
   end if
 
   if(present(fullsave)) then
-     fullsv=fullsave
+    fullsv=fullsave
   else
-     fullsv=.FALSE.
+    fullsv=.FALSE.
   end if
 
   if(zpre.eq.' ') then
-     call cdfOpn(ncid, filename, 'w') ! no error flag??
-     imodify=.FALSE.
+    call cdfOpn(ncid, filename, 'w') ! no error flag??
+    imodify=.FALSE.
   else
-     call cdfOpn(ncid, filename, 'v') ! Append if exists, else create new
-     imodify=.TRUE.
+    call cdfOpn(ncid, filename, 'v') ! Append if exists, else create new
+    imodify=.TRUE.
   end if
   if(ncid==0) then
-     ier = 39
-     return
+    ier = 39
+    return
   end if
 
   dimlens = (/1, 1, 1/) ! scalar
@@ -92,16 +90,16 @@ subroutine EZspline_save1(spline_o, filename, ier, &
   !
   dimlens = (/spline_o%n1, 1, 1/)
   if(.not.fullsv) then
-     call ezspline_defVar(ncid, trim(zpre)//'f', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'f', dimlens, real8, imodify, ier)
   else
-     dimlens = (/1, 1, 1/) ! scalar
-     call ezspline_defVar(ncid, trim(zpre)//'x1min', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x1max', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'ilin1', dimlens, int, imodify, ier)
-     dimlens = (/spline_o%n1, 4, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'x1pkg', dimlens, real8, imodify, ier)
-     dimlens = (/2, spline_o%n1, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'fspl', dimlens, real8, imodify, ier)
+    dimlens = (/1, 1, 1/) ! scalar
+    call ezspline_defVar(ncid, trim(zpre)//'x1min', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x1max', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'ilin1', dimlens, int, imodify, ier)
+    dimlens = (/spline_o%n1, 4, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'x1pkg', dimlens, real8, imodify, ier)
+    dimlens = (/2, spline_o%n1, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'fspl', dimlens, real8, imodify, ier)
   end if
   if(ier.ne.0) return
 
@@ -115,22 +113,22 @@ subroutine EZspline_save1(spline_o, filename, ier, &
   call cdfPutVar(ncid, trim(zpre)//'bcval1min', spline_o%bcval1min, ifail)
   call cdfPutVar(ncid, trim(zpre)//'bcval1max', spline_o%bcval1max, ifail)
   if(ifail/=0) then
-     ier=40
-     return
+    ier=40
+    return
   end if
 
   if(spline_o%isReady == 1) then
-     if(.not.fullsv) then
-        call cdfPutVar(ncid, trim(zpre)//'f', spline_o%fspl(1,:), ifail)
-     else
-        call cdfPutVar(ncid, trim(zpre)//'x1min', spline_o%x1min, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x1max', spline_o%x1max, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'ilin1', spline_o%ilin1, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x1pkg', spline_o%x1pkg, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'fspl', spline_o%fspl, ifail)
-     end if
+    if(.not.fullsv) then
+      call cdfPutVar(ncid, trim(zpre)//'f', spline_o%fspl(1,:), ifail)
+    else
+      call cdfPutVar(ncid, trim(zpre)//'x1min', spline_o%x1min, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x1max', spline_o%x1max, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'ilin1', spline_o%ilin1, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x1pkg', spline_o%x1pkg, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'fspl', spline_o%fspl, ifail)
+    end if
   else
-     ier = 93
+    ier = 93
   end if
 
   call cdfCls(ncid)
@@ -145,8 +143,7 @@ end subroutine EZspline_save1
 !! 2-D
 !!
 
-subroutine EZspline_save2(spline_o, filename, ier, &
-     spl_name, fullsave)
+subroutine EZspline_save2(spline_o, filename, ier, spl_name, fullsave)
   use EZspline_obj
   use EZcdf
   implicit none
@@ -159,11 +156,11 @@ subroutine EZspline_save2(spline_o, filename, ier, &
   !  if SPL_NAME is set, APPEND to file instead of creating new; prepend
   !  name to all NetCDF data items.  This allows one file to contain
   !  multiple items. (new dmc Mar 2006)
-  character(len=*), intent(in), OPTIONAL :: SPL_NAME  ! optional spline "name"
+  character(len=*), intent(in), optional :: spl_name  ! optional spline "name"
 
   !  if FULLSAVE is set .TRUE., save derived coefficients along with data
   !  this saves recalculating the coefficients when file is read (Mar 2006)
-  logical, intent(in), OPTIONAL :: FULLSAVE
+  logical, intent(in), optional :: fullsave
 
   integer ncid, ifail, in0, in1, in2
   integer dimlens(3)
@@ -175,8 +172,8 @@ subroutine EZspline_save2(spline_o, filename, ier, &
 
   ier = 0
   if(spline_o%isReady /= 1) then
-     ier = 94
-     return
+    ier = 94
+    return
   end if
 
   in0 = size(spline_o%fspl,1)
@@ -184,29 +181,29 @@ subroutine EZspline_save2(spline_o, filename, ier, &
   in2 = size(spline_o%fspl,3)
 
   if(present(spl_name)) then
-     call ezspline_spl_name_chk(spl_name,ier)
-     if(ier.ne.0) return
-     zpre=spl_name
+    call ezspline_spl_name_chk(spl_name,ier)
+    if(ier.ne.0) return
+    zpre=spl_name
   else
-     zpre=' '
+    zpre=' '
   end if
 
   if(present(fullsave)) then
-     fullsv=fullsave
+    fullsv=fullsave
   else
-     fullsv=.FALSE.
+    fullsv=.FALSE.
   end if
 
   if(zpre.eq.' ') then
-     call cdfOpn(ncid, filename, 'w') ! no error flag??
-     imodify=.FALSE.
+    call cdfOpn(ncid, filename, 'w') ! no error flag??
+    imodify=.FALSE.
   else
-     call cdfOpn(ncid, filename, 'v') ! Append if exists, else create new
-     imodify=.TRUE.
+    call cdfOpn(ncid, filename, 'v') ! Append if exists, else create new
+    imodify=.TRUE.
   end if
   if(ncid==0) then
-     ier = 39
-     return
+    ier = 39
+    return
   end if
 
   dimlens = (/1, 1, 1/) ! scalar
@@ -238,21 +235,21 @@ subroutine EZspline_save2(spline_o, filename, ier, &
   !
   dimlens = (/in1, in2, 1/)
   if(.not.fullsv) then
-     call ezspline_defVar(ncid, trim(zpre)//'f', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'f', dimlens, real8, imodify, ier)
   else
-     dimlens = (/1, 1, 1/) ! scalar
-     call ezspline_defVar(ncid, trim(zpre)//'x1min', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x1max', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'ilin1', dimlens, int, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x2min', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x2max', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'ilin2', dimlens, int, imodify, ier)
-     dimlens = (/spline_o%n1, 4, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'x1pkg', dimlens, real8, imodify, ier)
-     dimlens = (/spline_o%n2, 4, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'x2pkg', dimlens, real8, imodify, ier)
-     dimlens = (/in0, in1, in2/)
-     call ezspline_defVar(ncid, trim(zpre)//'fspl', dimlens, real8, imodify, ier)
+    dimlens = (/1, 1, 1/) ! scalar
+    call ezspline_defVar(ncid, trim(zpre)//'x1min', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x1max', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'ilin1', dimlens, int, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x2min', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x2max', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'ilin2', dimlens, int, imodify, ier)
+    dimlens = (/spline_o%n1, 4, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'x1pkg', dimlens, real8, imodify, ier)
+    dimlens = (/spline_o%n2, 4, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'x2pkg', dimlens, real8, imodify, ier)
+    dimlens = (/in0, in1, in2/)
+    call ezspline_defVar(ncid, trim(zpre)//'fspl', dimlens, real8, imodify, ier)
   end if
   if(ier.ne.0) return
 
@@ -274,26 +271,26 @@ subroutine EZspline_save2(spline_o, filename, ier, &
   call cdfPutVar(ncid, trim(zpre)//'bcval2min', spline_o%bcval2min, ifail)
   call cdfPutVar(ncid, trim(zpre)//'bcval2max', spline_o%bcval2max, ifail)
   if(ifail/=0) then
-     ier=40
-     return
+    ier=40
+    return
   end if
 
   if(spline_o%isReady == 1) then
-     if(.not.fullsv) then
-        call cdfPutVar(ncid, trim(zpre)//'f', spline_o%fspl(1,:,:), ifail)
-     else
-        call cdfPutVar(ncid, trim(zpre)//'x1min', spline_o%x1min, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x1max', spline_o%x1max, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'ilin1', spline_o%ilin1, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x2min', spline_o%x2min, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x2max', spline_o%x2max, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'ilin2', spline_o%ilin2, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x1pkg', spline_o%x1pkg, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x2pkg', spline_o%x2pkg, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'fspl', spline_o%fspl, ifail)
-     end if
+    if(.not.fullsv) then
+      call cdfPutVar(ncid, trim(zpre)//'f', spline_o%fspl(1,:,:), ifail)
+    else
+      call cdfPutVar(ncid, trim(zpre)//'x1min', spline_o%x1min, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x1max', spline_o%x1max, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'ilin1', spline_o%ilin1, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x2min', spline_o%x2min, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x2max', spline_o%x2max, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'ilin2', spline_o%ilin2, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x1pkg', spline_o%x1pkg, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x2pkg', spline_o%x2pkg, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'fspl', spline_o%fspl, ifail)
+    end if
   else
-     ier = 93
+    ier = 93
   end if
 
   call cdfCls(ncid)
@@ -308,8 +305,7 @@ end subroutine EZspline_save2
 !!! 3-D
 !!!
 
-subroutine EZspline_save3(spline_o, filename, ier, &
-     spl_name, fullsave)
+subroutine EZspline_save3(spline_o, filename, ier, spl_name, fullsave)
   use EZspline_obj
   use EZcdf
   implicit none
@@ -322,11 +318,11 @@ subroutine EZspline_save3(spline_o, filename, ier, &
   !  if SPL_NAME is set, APPEND to file instead of creating new; prepend
   !  name to all NetCDF data items.  This allows one file to contain
   !  multiple items. (new dmc Mar 2006)
-  character(len=*), intent(in), OPTIONAL :: SPL_NAME  ! optional spline "name"
+  character(len=*), intent(in), optional :: spl_name  ! optional spline "name"
 
   !  if FULLSAVE is set .TRUE., save derived coefficients along with data
   !  this saves recalculating the coefficients when file is read (Mar 2006)
-  logical, intent(in), OPTIONAL :: FULLSAVE
+  logical, intent(in), optional :: fullsave
 
   integer ncid, ifail, in0, in1, in2, in3
   integer dimlens(3)
@@ -338,8 +334,8 @@ subroutine EZspline_save3(spline_o, filename, ier, &
 
   ier = 0
   if(spline_o%isReady /= 1) then
-     ier = 94
-     return
+    ier = 94
+    return
   end if
 
   in0 = size(spline_o%fspl,1)
@@ -348,29 +344,29 @@ subroutine EZspline_save3(spline_o, filename, ier, &
   in3 = size(spline_o%fspl,4)
 
   if(present(spl_name)) then
-     call ezspline_spl_name_chk(spl_name,ier)
-     if(ier.ne.0) return
-     zpre=spl_name
+    call ezspline_spl_name_chk(spl_name,ier)
+    if(ier.ne.0) return
+    zpre=spl_name
   else
-     zpre=' '
+    zpre=' '
   end if
 
   if(present(fullsave)) then
-     fullsv=fullsave
+    fullsv=fullsave
   else
-     fullsv=.FALSE.
+    fullsv=.FALSE.
   end if
 
   if(zpre.eq.' ') then
-     call cdfOpn(ncid, filename, 'w') ! no error flag??
-     imodify=.FALSE.
+    call cdfOpn(ncid, filename, 'w') ! no error flag??
+    imodify=.FALSE.
   else
-     call cdfOpn(ncid, filename, 'v') ! Append if exists, else create new
-     imodify=.TRUE.
+    call cdfOpn(ncid, filename, 'v') ! Append if exists, else create new
+    imodify=.TRUE.
   end if
   if(ncid==0) then
-     ier = 39
-     return
+    ier = 39
+    return
   end if
 
   dimlens = (/1, 1, 1/) ! scalar
@@ -410,26 +406,26 @@ subroutine EZspline_save3(spline_o, filename, ier, &
   !
   dimlens = (/in1, in2, in3/)
   if(.not.fullsv) then
-     call ezspline_defVar(ncid, trim(zpre)//'f', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'f', dimlens, real8, imodify, ier)
   else
-     dimlens = (/1, 1, 1/) ! scalar
-     call ezspline_defVar(ncid, trim(zpre)//'x1min', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x1max', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'ilin1', dimlens, int, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x2min', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x2max', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'ilin2', dimlens, int, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x3min', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'x3max', dimlens, real8, imodify, ier)
-     call ezspline_defVar(ncid, trim(zpre)//'ilin3', dimlens, int, imodify, ier)
-     dimlens = (/spline_o%n1, 4, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'x1pkg', dimlens, real8, imodify, ier)
-     dimlens = (/spline_o%n2, 4, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'x2pkg', dimlens, real8, imodify, ier)
-     dimlens = (/spline_o%n3, 4, 1/)
-     call ezspline_defVar(ncid, trim(zpre)//'x3pkg', dimlens, real8, imodify, ier)
-     dimlens = (/in0*in1, in2, in3/)
-     call ezspline_defVar(ncid, trim(zpre)//'fspl', dimlens, real8, imodify, ier)
+    dimlens = (/1, 1, 1/) ! scalar
+    call ezspline_defVar(ncid, trim(zpre)//'x1min', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x1max', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'ilin1', dimlens, int, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x2min', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x2max', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'ilin2', dimlens, int, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x3min', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'x3max', dimlens, real8, imodify, ier)
+    call ezspline_defVar(ncid, trim(zpre)//'ilin3', dimlens, int, imodify, ier)
+    dimlens = (/spline_o%n1, 4, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'x1pkg', dimlens, real8, imodify, ier)
+    dimlens = (/spline_o%n2, 4, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'x2pkg', dimlens, real8, imodify, ier)
+    dimlens = (/spline_o%n3, 4, 1/)
+    call ezspline_defVar(ncid, trim(zpre)//'x3pkg', dimlens, real8, imodify, ier)
+    dimlens = (/in0*in1, in2, in3/)
+    call ezspline_defVar(ncid, trim(zpre)//'fspl', dimlens, real8, imodify, ier)
   end if
   if(ier.ne.0) return
 
@@ -457,31 +453,31 @@ subroutine EZspline_save3(spline_o, filename, ier, &
   call cdfPutVar(ncid, trim(zpre)//'bcval3min', spline_o%bcval3min, ifail)
   call cdfPutVar(ncid, trim(zpre)//'bcval3max', spline_o%bcval3max, ifail)
   if(ifail/=0) then
-     ier=40
-     return
+    ier=40
+    return
   end if
 
   if(spline_o%isReady == 1) then
-     if(.not.fullsv) then
-        call cdfPutVar(ncid, trim(zpre)//'f', spline_o%fspl(1,:,:,:), ifail)
-     else
-        call cdfPutVar(ncid, trim(zpre)//'x1min', spline_o%x1min, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x1max', spline_o%x1max, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'ilin1', spline_o%ilin1, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x2min', spline_o%x2min, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x2max', spline_o%x2max, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'ilin2', spline_o%ilin2, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x3min', spline_o%x3min, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x3max', spline_o%x3max, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'ilin3', spline_o%ilin3, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x1pkg', spline_o%x1pkg, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x2pkg', spline_o%x2pkg, ifail)
-        call cdfPutVar(ncid, trim(zpre)//'x3pkg', spline_o%x3pkg, ifail)
-        call ezspline_cdfput3(ncid, trim(zpre)//'fspl', spline_o%fspl, &
-             in0*in1, in2, in3, ifail)
-     end if
+    if(.not.fullsv) then
+      call cdfPutVar(ncid, trim(zpre)//'f', spline_o%fspl(1,:,:,:), ifail)
+    else
+      call cdfPutVar(ncid, trim(zpre)//'x1min', spline_o%x1min, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x1max', spline_o%x1max, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'ilin1', spline_o%ilin1, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x2min', spline_o%x2min, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x2max', spline_o%x2max, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'ilin2', spline_o%ilin2, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x3min', spline_o%x3min, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x3max', spline_o%x3max, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'ilin3', spline_o%ilin3, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x1pkg', spline_o%x1pkg, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x2pkg', spline_o%x2pkg, ifail)
+      call cdfPutVar(ncid, trim(zpre)//'x3pkg', spline_o%x3pkg, ifail)
+      call ezspline_cdfput3(ncid, trim(zpre)//'fspl', spline_o%fspl, &
+           in0*in1, in2, in3, ifail)
+    end if
   else
-     ier = 93
+    ier = 93
   end if
 
   call cdfCls(ncid)
@@ -506,19 +502,19 @@ subroutine ezspline_spl_name_chk(spl_name,ier)
   ilen = len(trim(spl_name))
 
   if(ilen.le.0) then
-     ier=50
+    ier=50
   else if(ilen.gt.20) then
-     ier=51
+    ier=51
   else
-     indx=index(zlegal,spl_name(1:1))
-     if(indx.le.10) then
-        ier=52
-     else
-        do ic=2,ilen
-           indx=index(zlegal,spl_name(ic:ic))
-           if(indx.le.0) ier=52
-        end do
-     end if
+    indx=index(zlegal,spl_name(1:1))
+    if(indx.le.10) then
+      ier=52
+    else
+      do ic=2,ilen
+        indx=index(zlegal,spl_name(ic:ic))
+        if(indx.le.0) ier=52
+      end do
+    end if
   end if
 
 end subroutine ezspline_spl_name_chk
@@ -546,42 +542,31 @@ subroutine ezspline_defVar(ncid, name, dimlens, ztype, imodify, ier)
   !-----------------------------
 
   if(.not.imodify) then
-
-     ! simply define the quantity...
-
-     call cdfDefVar(ncid, name, dimlens, ztype, ifail)
-     if(ifail.ne.0) ier=42
+    ! simply define the quantity...
+    call cdfDefVar(ncid, name, dimlens, ztype, ifail)
+    if(ifail.ne.0) ier=42
 
   else
-
-     ifail = nf_inq_varid(ncid, name, ivarid)
-     if(ifail.ne.0) then
-
-        ! assume item does not exist, so, just define it...
-
-        call cdfDefVar(ncid, name, dimlens, ztype, ifail)
-        if(ifail.ne.0) ier=42
-
-     else
-
-        ! item DOES exist; require dimension & type match...
-
-        call cdfInqVar(ncid, name, dimlens_loc, ztype_loc, ifail)
-        if(ztype.ne.ztype_loc) then
-           ier=53
-
-        else
-           do ii=1,3
-              id1=max(1,dimlens(ii))
-              id2=max(1,dimlens_loc(ii))
-              if(id1.ne.id2) ier=53
-           end do
-        end if
-
-        !  if ier= 0, no cdfDefVar is needed: object already defined
-        !  and attributes are correct.
-
-     end if
+    ifail = nf_inq_varid(ncid, name, ivarid)
+    if(ifail.ne.0) then
+      ! assume item does not exist, so, just define it...
+      call cdfDefVar(ncid, name, dimlens, ztype, ifail)
+      if(ifail.ne.0) ier=42
+    else
+      ! item DOES exist; require dimension & type match...
+      call cdfInqVar(ncid, name, dimlens_loc, ztype_loc, ifail)
+      if(ztype.ne.ztype_loc) then
+        ier=53
+      else
+        do ii=1,3
+          id1=max(1,dimlens(ii))
+          id2=max(1,dimlens_loc(ii))
+          if(id1.ne.id2) ier=53
+        end do
+      end if
+      !  if ier= 0, no cdfDefVar is needed: object already defined
+      !  and attributes are correct.
+    end if
   end if
 
 end subroutine ezspline_defVar
diff --git a/ezspline/ezspline_setup.f90 b/ezspline/ezspline_setup.f90
index ed83ebb..833e626 100644
--- a/ezspline/ezspline_setup.f90
+++ b/ezspline/ezspline_setup.f90
@@ -2,7 +2,7 @@ subroutine EZspline_setup1(spline_o, f, ier, exact_dim)
   use precision_mod, only: fp
   use EZspline_obj
   implicit none
-  type(EZspline1) spline_o
+  type(EZspline1) :: spline_o
   real(fp), dimension(:), intent(in) :: f
   ! ier:
   ! 0=ok
@@ -13,19 +13,19 @@ subroutine EZspline_setup1(spline_o, f, ier, exact_dim)
   !  dimensioning match between f and spline_o%fspl; default is F
 
   logical :: iexact
-  integer ifail
+  integer :: ifail
 
   integer :: ipx
-  integer iper, imsg, itol, inum, in1
-  real(fp) ztol, df1, df2
+  integer :: iper, imsg, itol, inum, in1
+  real(fp) :: ztol, df1, df2
 
   !-------------------------
   iexact=.FALSE.
   if(present(exact_dim)) iexact = exact_dim
 
   if( .not.EZspline_allocated(spline_o) ) then
-     ier = 98
-     return
+    ier = 98
+    return
   end if
 
   in1 = size(spline_o%fspl,2)
@@ -34,8 +34,8 @@ subroutine EZspline_setup1(spline_o, f, ier, exact_dim)
   if(size(f,1).lt.in1) return
 
   if(iexact) then
-     ier = 58
-     if(size(f,1).gt.in1) return
+    ier = 58
+    if(size(f,1).gt.in1) return
   end if
   ier = 0
 
@@ -51,90 +51,85 @@ subroutine EZspline_setup1(spline_o, f, ier, exact_dim)
   iper=0
   if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
   call genxpkg(spline_o%n1, spline_o%x1(1), spline_o%x1pkg(1,1),&
-       & iper,imsg,itol,ztol, spline_o%klookup1 ,ifail)
+       iper,imsg,itol,ztol, spline_o%klookup1 ,ifail)
   if(ifail/=0) ier=27
 
   spline_o%isReady = 0
 
-  spline_o%fspl(1, 1:in1) = &
-       &  f(1:in1)
+  spline_o%fspl(1,1:in1) = f(1:in1)
 
   if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
 
-     call mkspline(             &
-          &   spline_o%x1(1), spline_o%n1, &
-          &   spline_o%fspl(1,1), &
-          &   spline_o%ibctype1(1), spline_o%bcval1min, &
-          &   spline_o%ibctype1(2), spline_o%bcval1max, &
-          &   spline_o%ilin1, ifail)
+    call mkspline(spline_o%x1(1), spline_o%n1, spline_o%fspl(1,1), &
+         spline_o%ibctype1(1), spline_o%bcval1min, &
+         spline_o%ibctype1(2), spline_o%bcval1max, &
+         spline_o%ilin1, ifail)
 
-     if(ifail /= 0) then
-        ier = 98
-     else
-        spline_o%isReady = 1
-     end if
+    if(ifail /= 0) then
+      ier = 98
+    else
+      spline_o%isReady = 1
+    end if
 
   else if (spline_o%isLinear == 1) then
 
-     spline_o%ilin1=0
-     if(in1.gt.2) then
-        if(spline_o%x1pkg(3,4).eq.0.0_fp) spline_o%ilin1=1  ! evenly spaced grid
-     else
-        spline_o%ilin1=1
-     end if
+    spline_o%ilin1=0
+    if(in1.gt.2) then
+      if(spline_o%x1pkg(3,4).eq.0.0_fp) spline_o%ilin1=1  ! evenly spaced grid
+    else
+      spline_o%ilin1=1
+    end if
 
-     spline_o%isReady = 1   ! no coefficient setup necessary
+    spline_o%isReady = 1   ! no coefficient setup necessary
 
   else
-     !
-     ! Hermite polynomial coefficient setup based on Akima's method
-     ! (df/dx..etc not required)
-     !
-     ipx = 0
-     if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
-        ipx=1
-     else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
-          spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
-        ipx=99  ! an error...
-     else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
-        ipx=2
-        if(spline_o%ibctype1(1)==1) then
-           spline_o%fspl(2,1)=spline_o%bcval1min
-        else
-           ! hand implemented default BC
-           df1=(spline_o%fspl(1,2)-spline_o%fspl(1,1))/ &
-                (spline_o%x1(2)-spline_o%x1(1))
-           df2=(spline_o%fspl(1,3)-spline_o%fspl(1,2))/ &
-                (spline_o%x1(3)-spline_o%x1(2))
-           spline_o%fspl(2,1)=(3*df1-df2)/2
-        end if
-        inum=spline_o%n1
-        if(spline_o%ibctype1(2)==1) then
-           spline_o%fspl(2,inum)=spline_o%bcval1max
-        else
-           ! hand implemented default BC
-           df1=(spline_o%fspl(1,inum)-spline_o%fspl(1,inum-1))/ &
-                (spline_o%x1(inum)-spline_o%x1(inum-1))
-           df2=(spline_o%fspl(1,inum-1)-spline_o%fspl(1,inum-2))/ &
-                (spline_o%x1(inum-1)-spline_o%x1(inum-2))
-           spline_o%fspl(2,inum)=(3*df1-df2)/2
-        end if
-     end if
-     ifail = 0
-     call akherm1p(spline_o%x1(1), spline_o%n1, &
-          & spline_o%fspl(1,1), &
-          & spline_o%ilin1, &
-          & ipx, ifail)
-
-     if (ifail /=0 ) then
-        ier = 91
-     else
-        spline_o%isReady = 1
-     end if
-
+    !
+    ! Hermite polynomial coefficient setup based on Akima's method
+    ! (df/dx..etc not required)
+    !
+    ipx = 0
+    if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
+      ipx=1
+    else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
+         spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
+      ipx=99  ! an error...
+    else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
+      ipx=2
+      if(spline_o%ibctype1(1)==1) then
+        spline_o%fspl(2,1)=spline_o%bcval1min
+      else
+        ! hand implemented default BC
+        df1=(spline_o%fspl(1,2)-spline_o%fspl(1,1))/ &
+             (spline_o%x1(2)-spline_o%x1(1))
+        df2=(spline_o%fspl(1,3)-spline_o%fspl(1,2))/ &
+             (spline_o%x1(3)-spline_o%x1(2))
+        spline_o%fspl(2,1)=(3*df1-df2)/2
+      end if
+      inum=spline_o%n1
+      if(spline_o%ibctype1(2)==1) then
+        spline_o%fspl(2,inum)=spline_o%bcval1max
+      else
+        ! hand implemented default BC
+        df1=(spline_o%fspl(1,inum)-spline_o%fspl(1,inum-1))/ &
+             (spline_o%x1(inum)-spline_o%x1(inum-1))
+        df2=(spline_o%fspl(1,inum-1)-spline_o%fspl(1,inum-2))/ &
+             (spline_o%x1(inum-1)-spline_o%x1(inum-2))
+        spline_o%fspl(2,inum)=(3*df1-df2)/2
+      end if
+    end if
+    ifail = 0
+    call akherm1p(spline_o%x1(1), spline_o%n1, &
+         spline_o%fspl(1,1), &
+         spline_o%ilin1, &
+         ipx, ifail)
+    if (ifail /=0 ) then
+      ier = 91
+    else
+      spline_o%isReady = 1
+    end if
   end if
 
-
+  return
 end subroutine EZspline_setup1
 
 
@@ -198,18 +193,17 @@ subroutine EZspline_setup2(spline_o, f, ier, exact_dim)
   iper=0
   if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
   call genxpkg(spline_o%n1,spline_o%x1(1),spline_o%x1pkg(1,1),&
-       & iper,imsg,itol,ztol,spline_o%klookup1,ifail)
+       iper,imsg,itol,ztol,spline_o%klookup1,ifail)
   if(ifail/=0) ier=27
   iper=0
   if(spline_o%ibctype2(1)==-1 .OR. spline_o%ibctype2(2)==-1) iper=1
   call genxpkg(spline_o%n2,spline_o%x2(1),spline_o%x2pkg(1,1),&
-       & iper,imsg,itol,ztol,spline_o%klookup2,ifail)
+       iper,imsg,itol,ztol,spline_o%klookup2,ifail)
   if(ifail/=0) ier=27
 
   spline_o%isReady = 0
 
-  spline_o%fspl(1, 1:in1, 1:in2) = &
-       &  f(1:in1, 1:in2)
+  spline_o%fspl(1,1:in1,1:in2) = f(1:in1,1:in2)
 
   ! this fixes a VMS f90 compiler optimizer problem:
   if(ztol.eq.-1.2345d30) &
@@ -218,146 +212,145 @@ subroutine EZspline_setup2(spline_o, f, ier, exact_dim)
   if (spline_o%isHybrid == 1) then
 
      call mkintrp2d( &
-          &   spline_o%x1(1), spline_o%n1, &
-          &   spline_o%x2(1), spline_o%n2, &
-          &   spline_o%hspline, spline_o%fspl(1,1,1), &
-          &   in0,in1,in2, &
-          &   spline_o%ibctype1(1), spline_o%bcval1min(1), &
-          &   spline_o%ibctype1(2), spline_o%bcval1max(1), &
-          &   spline_o%ibctype2(1), spline_o%bcval2min(1), &
-          &   spline_o%ibctype2(2), spline_o%bcval2max(1), &
-          &   ifail)
+          spline_o%x1(1), spline_o%n1, &
+          spline_o%x2(1), spline_o%n2, &
+          spline_o%hspline, spline_o%fspl(1,1,1), &
+          in0,in1,in2, &
+          spline_o%ibctype1(1), spline_o%bcval1min(1), &
+          spline_o%ibctype1(2), spline_o%bcval1max(1), &
+          spline_o%ibctype2(1), spline_o%bcval2min(1), &
+          spline_o%ibctype2(2), spline_o%bcval2max(1), &
+          ifail)
 
      spline_o%ilin1 = 0  ! Hybrid does not compute this
      spline_o%ilin2 = 0  ! Hybrid does not compute this
 
      if(ifail /= 0) then
-        ier = 98
+       ier = 98
      else
-        spline_o%isReady = 1
+       spline_o%isReady = 1
      end if
 
   else if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
 
      call mkbicub(             &
-          &   spline_o%x1(1), spline_o%n1, &
-          &   spline_o%x2(1), spline_o%n2, &
-          &   spline_o%fspl(1,1,1), spline_o%n1, &
-          &   spline_o%ibctype1(1), spline_o%bcval1min(1), &
-          &   spline_o%ibctype1(2), spline_o%bcval1max(1), &
-          &   spline_o%ibctype2(1), spline_o%bcval2min(1), &
-          &   spline_o%ibctype2(2), spline_o%bcval2max(1), &
-          &   spline_o%ilin1, spline_o%ilin2, ifail)
+          spline_o%x1(1), spline_o%n1, &
+          spline_o%x2(1), spline_o%n2, &
+          spline_o%fspl(1,1,1), spline_o%n1, &
+          spline_o%ibctype1(1), spline_o%bcval1min(1), &
+          spline_o%ibctype1(2), spline_o%bcval1max(1), &
+          spline_o%ibctype2(1), spline_o%bcval2min(1), &
+          spline_o%ibctype2(2), spline_o%bcval2max(1), &
+          spline_o%ilin1, spline_o%ilin2, ifail)
 
      if(ifail /= 0) then
-        ier = 98
+       ier = 98
      else
-        spline_o%isReady = 1
+       spline_o%isReady = 1
      end if
 
   else if (spline_o%isLinear == 1) then
-
-     spline_o%ilin1=0
-     if(in1.gt.2) then
-        if(spline_o%x1pkg(3,4).eq.0.0_fp) spline_o%ilin1=1  ! evenly spaced grid
-     else
-        spline_o%ilin1=1
-     end if
-
-     spline_o%ilin2=0
-     if(in2.gt.2) then
-        if(spline_o%x2pkg(3,4).eq.0.0_fp) spline_o%ilin2=1  ! evenly spaced grid
-     else
-        spline_o%ilin2=1
-     end if
-
-     spline_o%isReady = 1   ! no coefficient setup necessary
+    spline_o%ilin1=0
+    if(in1.gt.2) then
+      if(spline_o%x1pkg(3,4).eq.0.0_fp) spline_o%ilin1=1  ! evenly spaced grid
+    else
+      spline_o%ilin1=1
+    end if
+
+    spline_o%ilin2=0
+    if(in2.gt.2) then
+      if(spline_o%x2pkg(3,4).eq.0.0_fp) spline_o%ilin2=1  ! evenly spaced grid
+    else
+      spline_o%ilin2=1
+    end if
+
+    spline_o%isReady = 1   ! no coefficient setup necessary
 
   else
-     !
-     ! Hermite polynomial coefficient setup based on Akima's method
-     ! (df/dx..etc not required)
-     !
-     ipx = 0
-     if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
-        ipx=1
-     else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
-          spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
-        ipx=99  ! an error...
-     else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
-        ipx=2
-        do jj=1,spline_o%n2
-           if(spline_o%ibctype1(1)==1) then
-              spline_o%fspl(2,1,jj)=spline_o%bcval1min(jj)
-           else
-              ! hand implemented default BC
-              df1=(spline_o%fspl(1,2,jj)-spline_o%fspl(1,1,jj))/ &
-                   (spline_o%x1(2)-spline_o%x1(1))
-              df2=(spline_o%fspl(1,3,jj)-spline_o%fspl(1,2,jj))/ &
-                   (spline_o%x1(3)-spline_o%x1(2))
-              spline_o%fspl(2,1,jj)=(3*df1-df2)/2
-           end if
-           inum=spline_o%n1
-           if(spline_o%ibctype1(2)==1) then
-              spline_o%fspl(2,inum,jj)=spline_o%bcval1max(jj)
-           else
-              ! hand implemented default BC
-              df1=(spline_o%fspl(1,inum,jj)-spline_o%fspl(1,inum-1,jj))/ &
-                   (spline_o%x1(inum)-spline_o%x1(inum-1))
-              df2=(spline_o%fspl(1,inum-1,jj)-spline_o%fspl(1,inum-2,jj))/ &
-                   (spline_o%x1(inum-1)-spline_o%x1(inum-2))
-              spline_o%fspl(2,inum,jj)=(3*df1-df2)/2
-           end if
-        end do
-     end if
-     ipy = 0
-     if (spline_o%ibctype2(1)==-1 .or. spline_o%ibctype2(2)==-1) then
-        ipy=1
-     else if (spline_o%ibctype2(1)<-1 .or. spline_o%ibctype2(1)>1 .or. &
-          spline_o%ibctype2(2)<-1 .or. spline_o%ibctype2(2)>1 ) then
-        ipy=99  ! an error...
-     else if (spline_o%ibctype2(1)==1 .or. spline_o%ibctype2(2)==1) then
-        ipy=2
-        do ii=1,spline_o%n1
-           if(spline_o%ibctype2(1)==1) then
-              spline_o%fspl(3,ii,1)=spline_o%bcval2min(ii)
-           else
-              ! hand implemented default BC
-              df1=(spline_o%fspl(1,ii,2)-spline_o%fspl(1,ii,1))/ &
-                   (spline_o%x2(2)-spline_o%x2(1))
-              df2=(spline_o%fspl(1,ii,3)-spline_o%fspl(1,ii,2))/ &
-                   (spline_o%x2(3)-spline_o%x2(2))
-              spline_o%fspl(3,ii,1)=(3*df1-df2)/2
-           end if
-           inum=spline_o%n2
-           if(spline_o%ibctype2(2)==1) then
-              spline_o%fspl(3,ii,inum)=spline_o%bcval2max(ii)
-           else
-              ! hand implemented default BC
-              df1=(spline_o%fspl(1,ii,inum)-spline_o%fspl(1,ii,inum-1))/ &
-                   (spline_o%x2(inum)-spline_o%x2(inum-1))
-              df2=(spline_o%fspl(1,ii,inum-1)-spline_o%fspl(1,ii,inum-2))/ &
-                   (spline_o%x2(inum-1)-spline_o%x2(inum-2))
-              spline_o%fspl(3,ii,inum)=(3*df1-df2)/2
-           end if
-        end do
-     end if
-     ifail = 0
-     call akherm2p(spline_o%x1(1), spline_o%n1, &
-          &         spline_o%x2(1), spline_o%n2, &
-          & spline_o%fspl(1,1,1), spline_o%n1, &
-          & spline_o%ilin1, spline_o%ilin2, &
-          & ipx,ipy, ifail)
-
-     if (ifail /=0 ) then
-        ier = 91
-     else
-        spline_o%isReady = 1
-     end if
+    !
+    ! Hermite polynomial coefficient setup based on Akima's method
+    ! (df/dx..etc not required)
+    !
+    ipx = 0
+    if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
+      ipx=1
+    else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
+         spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
+      ipx=99  ! an error...
+    else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
+      ipx=2
+      do jj=1,spline_o%n2
+        if(spline_o%ibctype1(1)==1) then
+          spline_o%fspl(2,1,jj)=spline_o%bcval1min(jj)
+        else
+          ! hand implemented default BC
+          df1=(spline_o%fspl(1,2,jj)-spline_o%fspl(1,1,jj))/ &
+               (spline_o%x1(2)-spline_o%x1(1))
+          df2=(spline_o%fspl(1,3,jj)-spline_o%fspl(1,2,jj))/ &
+               (spline_o%x1(3)-spline_o%x1(2))
+          spline_o%fspl(2,1,jj)=(3*df1-df2)/2
+        end if
+        inum=spline_o%n1
+        if(spline_o%ibctype1(2)==1) then
+          spline_o%fspl(2,inum,jj)=spline_o%bcval1max(jj)
+        else
+          ! hand implemented default BC
+          df1=(spline_o%fspl(1,inum,jj)-spline_o%fspl(1,inum-1,jj))/ &
+               (spline_o%x1(inum)-spline_o%x1(inum-1))
+          df2=(spline_o%fspl(1,inum-1,jj)-spline_o%fspl(1,inum-2,jj))/ &
+               (spline_o%x1(inum-1)-spline_o%x1(inum-2))
+          spline_o%fspl(2,inum,jj)=(3*df1-df2)/2
+        end if
+      end do
+    end if
+    ipy = 0
+    if (spline_o%ibctype2(1)==-1 .or. spline_o%ibctype2(2)==-1) then
+      ipy=1
+    else if (spline_o%ibctype2(1)<-1 .or. spline_o%ibctype2(1)>1 .or. &
+         spline_o%ibctype2(2)<-1 .or. spline_o%ibctype2(2)>1 ) then
+      ipy=99  ! an error...
+    else if (spline_o%ibctype2(1)==1 .or. spline_o%ibctype2(2)==1) then
+      ipy=2
+      do ii=1,spline_o%n1
+        if(spline_o%ibctype2(1)==1) then
+          spline_o%fspl(3,ii,1)=spline_o%bcval2min(ii)
+        else
+          ! hand implemented default BC
+          df1=(spline_o%fspl(1,ii,2)-spline_o%fspl(1,ii,1))/ &
+               (spline_o%x2(2)-spline_o%x2(1))
+          df2=(spline_o%fspl(1,ii,3)-spline_o%fspl(1,ii,2))/ &
+               (spline_o%x2(3)-spline_o%x2(2))
+          spline_o%fspl(3,ii,1)=(3*df1-df2)/2
+        end if
+        inum=spline_o%n2
+        if(spline_o%ibctype2(2)==1) then
+          spline_o%fspl(3,ii,inum)=spline_o%bcval2max(ii)
+        else
+          ! hand implemented default BC
+          df1=(spline_o%fspl(1,ii,inum)-spline_o%fspl(1,ii,inum-1))/ &
+               (spline_o%x2(inum)-spline_o%x2(inum-1))
+          df2=(spline_o%fspl(1,ii,inum-1)-spline_o%fspl(1,ii,inum-2))/ &
+               (spline_o%x2(inum-1)-spline_o%x2(inum-2))
+          spline_o%fspl(3,ii,inum)=(3*df1-df2)/2
+        end if
+      end do
+    end if
+    ifail = 0
+    call akherm2p(spline_o%x1(1), spline_o%n1, &
+         spline_o%x2(1), spline_o%n2, &
+         spline_o%fspl(1,1,1), spline_o%n1, &
+         spline_o%ilin1, spline_o%ilin2, &
+         ipx,ipy, ifail)
+
+    if (ifail /=0 ) then
+      ier = 91
+    else
+      spline_o%isReady = 1
+    end if
 
   end if
 
-
+  return
 end subroutine EZspline_setup2
 
 
@@ -640,5 +633,5 @@ subroutine EZspline_setup3(spline_o, f, ier, exact_dim)
 
   end if
 
-
+  return
 end subroutine EZspline_setup3
diff --git a/ezspline/qk_pspline.f90 b/ezspline/qk_pspline.f90
index 311c092..9243b0d 100644
--- a/ezspline/qk_pspline.f90
+++ b/ezspline/qk_pspline.f90
@@ -53,9 +53,9 @@ program qk_pspline
   character(len=8) :: zhdr
 
   real(fp) :: zsum,zorig,zdiff
-  real(fp),parameter :: zsumtol8 = 1.0d-4
-  real(fp),parameter :: zsumtol4 = 1.0d-3
-  real(fp),parameter :: epslon = 1.0d-34
+  real(fp),parameter :: zsumtol8 = 1.0e-4_fp
+  real(fp),parameter :: zsumtol4 = 1.0e-3_fp
+  real(fp),parameter :: epslon = 1.0e-34_fp
 
   !---------------------------------------------------------------------
   !  continuity test for values and derivatives...