fp = (yhi - rG) / dist

            global2local_j = fp*j + fm*(j+1)

      _RL rG, fp, fm, dist, zlo, zhi

            fp = (zhi - rG) / dist

            global2local_k = fp*k + fm*(k+1)

      _RL rG, fp, fm, dist, xlo, xhi

            fp = (xhi - rG) / dist

            global2local_i = fp*i + fm*(i+1)

      _RL rG, fp, fm, dist, ylo, yhi

      z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t)

      b  = z1 + dble(fp)

         discri = z1*z1-dble(fp)*dble(fpa)

      subroutine CUBIC( t, f, fp, ta, fa, fpa, tlower, tupper )

            if (fp.lt.0.) t = tupper

            if (fp.ge.0.) t = tlower

         discri = dble(fp)/z1

            if (fp.lt.0.) t = tupper

            if (fp.ge.0.) t = tlower

         t    = t + fp*(ta-t)/sngl(b+discri)

      double precision    t, f, fp, ta, fa, fpa, tlower, tupper

         fp = DDOT( nn, dd, 1, gg, 1 )

         if (fp .gt. xpara2*dotdg) then

         call CUBIC( tact, fnew, fp, ta, fa, fpa, left, right )

         call CUBIC( tact, fnew, fp, ta, fa, fpa, left, right )

         fpa = fp

      double precision    fa, fpa, fp

      _RL x1,x2,x3,we,pe,fw,fp,pa,pb,pc,ax,ba,bb,t1,ca,cb,t2

        fp=pe-float(ip-1)

        pa = coef1(ip,iw-1,nn)*(1.-fp)+coef1(ip+1,iw-1,nn)*fp

        pb = coef1(ip,iw,  nn)*(1.-fp)+coef1(ip+1,iw,  nn)*fp

        pc = coef1(ip,iw+1,nn)*(1.-fp)+coef1(ip+1,iw+1,nn)*fp

        ba = coef2(ip,iw,  nn)*(1.-fp)+coef2(ip+1,iw,  nn)*fp

        bb = coef2(ip,iw+1,nn)*(1.-fp)+coef2(ip+1,iw+1,nn)*fp

        ca = coef3(ip,iw,  nn)*(1.-fp)+coef3(ip+1,iw,  nn)*fp

        cb = coef3(ip,iw+1,nn)*(1.-fp)+coef3(ip+1,iw+1,nn)*fp