<p><b>mpetersen@lanl.gov</b> 2011-02-07 14:48:34 -0700 (Mon, 07 Feb 2011)</p><p>BRANCH COMMIT: Additions to module_spline_interpolation.F.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F        2011-02-07 20:17:24 UTC (rev 713)
+++ branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F        2011-02-07 21:48:34 UTC (rev 714)
@@ -6,7 +6,8 @@
 
   public ::   CubicSplineCoefficients, InterpolateCubicSpline, &amp;
     InterpolateColumnCubicSpline, IntegrateCubicSpline, &amp;
-    IntegrateColumnCubicSpline
+    IntegrateColumnCubicSpline, InterpolateColumnLinear, &amp;
+    TestInterpolateColumn
 
   contains
 
@@ -340,5 +341,133 @@
 
  end subroutine IntegrateColumnCubicSpline
 
+
+ subroutine InterpolateColumnLinear( &amp;
+               phiIn,phiOut, &amp;
+               zIn,zOut, &amp;
+               kmaxIn,kmaxOut)
+
+! !DESCRIPTION:
+!  Given kmaxIn pairs (zIn,phiIn), compute the linear interpolation
+!  phiOut at each zOut.
+!  This subroutine assumes that both zIn and zOut are monotonically
+!  increasing.  This should occur naturally if the z coordinates are
+!  computed from layer thicknesses in the calling routine.
+!  It also assumes that all values of zOut are within the first
+!  and last values of zIn.
+
+! !INPUT PARAMETERS:
+
+  real (kind=RKIND), dimension(kmaxIn), intent(in) :: &amp;
+    zIn,         &amp;! layer thickness, input grid
+    phiIn         ! interpolation variable, input grid
+
+  real (kind=RKIND), dimension(kmaxOut), intent(in) :: &amp;
+    zOut          ! layer thickness, output grid
+
+  integer, intent(in) :: &amp;
+    KmaxIn,      &amp;! number of layers, input grid
+    KmaxOut       ! number of layers, output grid
+
+! !OUTPUT PARAMETERS:
+
+  real (kind=RKIND), dimension(kmaxOut), intent(out) :: &amp;
+    phiOut        ! interpolation variable, output grid
+
+!-----------------------------------------------------------------------
+!
+!  local variables
+!
+!-----------------------------------------------------------------------
+
+  integer :: &amp;
+    kIn, kOut ! counters
+
+  kOut = 1
+
+  kInLoop: do kIn = 1,kmaxIn-1
+
+    do while(zOut(kOut) &lt; zIn(kIn+1)) 
+
+      phiOut(kOut) = phiIn(kIn)  &amp;
+        + (phiIn(kIn+1)-phiIn(kIn)) &amp;
+         /(zIn(kIn+1)  -zIn(kIn)  ) &amp;
+         *(zOut(kOut)  -zIn(kIn)  )
+
+      kOut = kOut + 1
+
+      if (kOut&gt;kmaxOut) exit kInLoop
+
+    enddo
+  
+  enddo kInLoop
+
+  end subroutine InterpolateColumnLinear
+
+
+subroutine TestInterpolateColumn
+
+!  Test function to show how to operate the cubic spline subroutines
+
+  integer, parameter :: &amp;
+    kmaxIn = 10
+  real (kind=RKIND), dimension(kmaxIn) :: &amp;
+    phiIn, zIn, phi2DerIn
+
+  integer, parameter :: &amp;
+    kmaxOut = 100
+  real (kind=RKIND), dimension(kmaxOut) :: &amp;
+    phiOut, zOut
+
+  integer :: &amp;
+    k
+
+!-----------------------------------------------------------------------
+!
+!  Create zIn, PhiIn, zOut
+!
+!-----------------------------------------------------------------------
+
+   do k=1,kmaxIn
+      zIn(k) = k-4
+      ! power function:
+      !phiIn(k) = zIn(k)**2
+      ! trig function:
+      phiIn(k) = sin(zIn(k)/2)
+   enddo
+
+   do k=1,kmaxOut
+      zOut(k) = zIn(1) + k/(kmaxOut+1.0)*(zIn(kmaxIn)-zIn(1))
+   enddo
+
+!-----------------------------------------------------------------------
+!
+!  Interpolate
+!
+!-----------------------------------------------------------------------
+
+   ! First, compute second derivative values at each node, phi2DerIn.
+   call CubicSplineCoefficients(zIn,phiIn,kmaxIn, &amp;
+      2.0e30, 2.0e30, phi2DerIn)
+
+   ! Compute interpolated values phiOut.
+   call InterpolateColumnCubicSpline(phiIn, phiOut, &amp;
+       phi2DerIn, zIn, zOut, kmaxIn, kmaxOut)
+
+!-----------------------------------------------------------------------
+!
+!  Output
+!
+!-----------------------------------------------------------------------
+
+! The following output can be copied directly into Matlab
+  print '(a,10f8.4,a)', 'zIn = [',zIn,'];'
+  print '(a,10f8.4,a)', 'phiIn = [',phiIn,'];'
+  print '(a,100f8.4,a)', 'zOut = [',zOut,'];'
+  print '(a,100f8.4,a)', 'phiOut = [',phiOut,'];'
+  print *, &quot;plot(zIn,phiIn,'-*r',zOut,phiOut,'x')&quot;
+
+end subroutine TestInterpolateColumn
+
 end module spline_interpolation
 

</font>
</pre>