copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
***s/r posiz3  -- Precomputing for vertical interpolation in vertint
*

      subroutine posiz3 (p,h,href,hint,ktop,kbot,ng,nk,nkref) 19
      implicit none
*
      integer ng,nk,nkref
      real p(ng,nk,2),href(ng,nkref),hint(ng,nk),
     $     h(ng,0:nkref+2),ktop(ng,nk),kbot(ng,nk)
*
* arguments:
* p(i,k,1)   o interpolation indexes for "vertint"
* p(i,k,2)   o precomputed values for interpolator "vertint"
* h          o monotonically increasing reference coordinates
* href       i reference coordinates
* hint       i interpolation locations
* ni,nj      i horizontal dimensions for all vectors
* nk         i vertical dimension for p and hint
* nkref      i vertical dimension for h and href
*
*  note: p and h should be directly passed to "vertint"
*
#include "vinterpo.cdk"
#include "nbcpu.cdk"
*
      integer i,k,kz,kc,kk
      real*8 capz,e1
*-------------------------------------------------------------------
*
      if (nk.gt.1) then
         if (hint(1,1).gt.hint(1,2)) then
            write (6,900)
            stop
         endif  
         if (href(1,1).gt.href(1,2)) then
            write (6,901)
            stop
         endif  
      endif
*
* Vertical halo for "vertint". Physically useless
* but yet imperative to allow for over indexing.
*
         do k=1,nkref
         do i=1,ng
            h(i,k) = href(i,k)
         end do
         end do
         do i=1,ng
            h(i,0)       = -href(i,2)
            h(i,nkref+1) = 1.1*href(i,nkref)
            h(i,nkref+2) = 1.2*href(i,nkref)
         end do
*
* Evaluate interpolation locations in terms of fractional index of href
* and precompute h2**2 / (h1+h2+h3) * capz * (capz-1) (Thomas and Cote)
*
      do k=1,nk        
         do i=1,ng
            ktop(i,k) = nkref
            kbot(i,k) = 1
         end do
         do kk=1,nkref/2
            do i=1,ng
               kc = ktop(i,k)+kbot(i,k)
               kc = ishft(kc,-1)
               if (h(i,kc).gt.hint(i,k)) then
                  ktop(i,k) = kc
               else
                  kbot(i,k) = kc
               endif
            end do
         end do
         do i=1,ng
            if (hint(i,k).ge.h(i,nkref)) kbot(i,k) = nkref
            if (hint(i,k).lt.h(i,    1)) kbot(i,k) = 1
         end do
         do i=1,ng
            kz = kbot(i,k) 
            capz=   max(0., ( hint(i,k   ) - h(i,kz) ) ) / 
     $                      (    h(i,kz+1) - h(i,kz) )
            if (kz.eq.nkref) capz = 0.
            p(i,k,1) = float(kz) + capz
            e1   = capz - 1.
            if ((kz.lt.2).or.(kz.ge.nkref-1)) capz = 0.
            p(i,k,2)= capz * e1 * ( h(i,kz+1) - h(i,kz  ) )**2./
     $                            ( h(i,kz+2) - h(i,kz-1) )
         end do
      end do
*
      if (v_interp.eq.'NEAREST') then  
         do k=1,nk        
         do i=1,ng
            p(i,k,1) = nint(p(i,k,1))
            p(i,k,2) = 0.
         end do
         end do
      endif
*
 900  format (/' Vertical interpolation locations must be presented in'/
     $         ' a monotonically increasing order -- POSIZ --'/)
 901  format (/' Vertical reference coordinates (along with data)'/
     $         ' must be presented in a monotonically increasing '/
     $         ' order -- POSIZ --'/)
*
*-------------------------------------------------------------------
       return
       end