copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
***s/r vertint2   -- Vertical interpolation routine
*

      subroutine vertint3 (fup,f,zd,href,ng,nk,nkref) 30
*
      implicit none
*
      integer ng,nk,nkref
      real    fup(ng,nk),f(ng,nkref),zd(ng,nk,2),href(ng,0:nkref+2)
*
* arguments:
* fup          o interpolated values
* f            i field to interpolate
* zd(i,j,k,1)  i interpolation indexes from "posiz"
* zd(i,j,k,2)  i precomputed values from "posiz"
* href         i monotonically increasing reference coordinates
* ng           i horizontal dimension for all vectors
* nk           i vertical dimension for fup and zd
* nkref        i vertical dimension for f and href
* 
#include "vinterpo.cdk"
#include "nbcpu.cdk"
*
      integer i, j, k, kz, kzm1, kzp1, kzp2
      real dum
      real*8  capz, e1, e2, e3, e4, x3mx2, x2mx1, f3mf2
*
*----------------------------------------------------------------------
*
      if ((v_interp.eq.'CUBIC_UQAM').or.(v_interp.eq.'NEAREST')) then  
* original UQAM
*
      do k = 1, nk
      do i = 1, ng
*
         kz   = int(zd(i,k,1))
         capz = zd(i,k,1) - float(kz)       
         e1=(1.0+2.0*capz) *(1.0-capz)*(1.0-capz)
         e2=(3.0-2.0*capz) *capz*capz
         e3=capz*(1.0-capz)*(1.0-capz)
         e4=(1.0-capz)*capz*capz
         if ((kz.lt.2).or.(kz.ge.nkref-1)) then
            e1= 1. - capz
            e2= capz
            e3= 0.
            e4= 0.
         endif
         kzm1 = max(1,kz - 1)
         kzp1 = min(nkref,kz + 1)
         kzp2 = min(nkref,kz + 2)
*
         x3mx2 =  href(i,kz+1) - href(i,kz)
         
         fup(i,k) = e1*f(i,kz  ) + e2*f(i,kzp1)   +
     $                e3*(f(i,kzp1)-f(i,kzm1))/
     $                      (href(i,kz+1)-href(i,kz-1))*x3mx2 -
     $                e4*(f(i,kzp2)-f(i,kz  ))/
     $                      (href(i,kz+2)-href(i,kz  ))*x3mx2
      end do
      end do
*
      elseif (v_interp.eq.'CUBIC_NEWTON') then
*     cubic based on Newton form of interpo. polynomials
*
      do k = 1, nk
      do i = 1, ng
*
         kz   = int(zd(i,k,1))
         capz = zd(i,k,1) - float(kz)
         e1   = 1.0 - capz
         kzm1 = max(1,kz - 1)
         kzp1 = min(nkref,kz + 1)
         kzp2 = min(nkref,kz + 2)
*
         x3mx2 =  href(i,kz+1) - href(i,kz  )
         x2mx1 =  href(i,kz  ) - href(i,kz-1)
         f3mf2 = (f(i,kzp1) - f(i,kz  )) / x3mx2
*
         fup(i,k) = e1 * f(i,kz  ) + capz * f(i,kzp1)   +

     $    ( f3mf2 - (f(i,kz  )-f(i,kzm1)) / x2mx1 ) /
     $      (href(i,kz+1)-href(i,kz-1) ) * zd(i,k,2) *
     $      (href(i,kz+2)-href(i,kz+1) + x3mx2*e1 ) +

     $    ( (f(i,kzp2)-f(i,kzp1))/(href(i,kz+2)-href(i,kz+1)) -
     $      f3mf2 )/ (href(i,kz+2)-href(i,kz  ) ) * zd(i,k,2) *
     $      ( x2mx1 +  x3mx2*capz )
*
      end do
      end do
*
      else
*
         write (6,101) v_interp
         stop
*
      endif
*
 101  format (/' -----  ABORT ABORT  in vertint -----'/
     $         ' Wrong choice of vertical interpolator: ',a/)
*-----------------------------------------------------------------------
      return
      end