LHAPDF/lhapdf5.5.1/src/wrapcteq6.f

   1 ! -*- F90 -*-
   2
   3
   4       subroutine CTEQ6evolve(x,Q,pdf)
   5       implicit real*8(a-h,o-z)
   6       include 'parmsetup.inc'
   7       character*16 name(nmxset)
   8       integer nmem(nmxset),ndef(nmxset),mmem
   9       common/NAME/name,nmem,ndef,mmem
  10       real*8 pdf(-6:6)
  11       integer nset
  12       Character Line*80
  13       PARAMETER (MXX = 96, MXQ = 20, MXF = 5, nhess = 40)
  14       PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX)
  15       Common                                                            &
  16      & / CtqPar1nhess /                                                 &
  17      &  Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset),                 &
  18      &  UPD(0:nhess,MXPQX,nmxset)                                       &
  19      & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset)                 &
  20      & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset)             &
  21      & / QCDtable /  Alambda(nmxset), Nfl(nmxset), Iorder(nmxset)       &
  22      & / Masstbl / Amass(6,nmxset)
  23       common/masses_LHA/cMass(nmxset),bMass(nmxset),tMass(nmxset)
  24       data pi / 3.141592653589793d0 /
  25       save
  26 !
  27       call getnset(iset)
  28       call getnmem(iset,imem)
  29 !
  30       U =         X * CtLhCtq6Pdf(imem,1,X,Q)
  31       D =         X * CtLhCtq6Pdf(imem,2,X,Q)
  32       USEA =      X * CtLhCtq6Pdf(imem,-1,X,Q)
  33       DSEA =      X * CtLhCtq6Pdf(imem,-2,X,Q)
  34       STR =       X * CtLhCtq6Pdf(imem,3,X,Q)
  35       CHM =       X * CtLhCtq6Pdf(imem,4,X,Q)
  36       BOT =       X * CtLhCtq6Pdf(imem,5,X,Q)
  37       GLU  =      X * CtLhCtq6Pdf(imem,0,X,Q)
  38 !
  39       pdf(0)  = glu
  40       pdf(1)  = d
  41       pdf(-1) = dsea
  42       pdf(2)  = u
  43       pdf(-2) = usea
  44       pdf(3)  = str
  45       pdf(-3) = str
  46       pdf(4)  = chm
  47       pdf(-4) = chm
  48       pdf(5)  = bot
  49       pdf(-5) = bot
  50       pdf(6)  = 0.0d0
  51       pdf(-6) = 0.0d0
  52
  53       return
  54 !
  55       entry CTEQ6read(nset)
  56
  57       call CtLhbldat1
  58       call CtLhbldat2
  59
  60       call LHct6set
  61                                             !*** nmem+1=number of member
  62       read(1,*)nmem(nset),ndef(nset)
  63       if(nmem(nset) .gt. nhess) then
  64          print *,'fatal error:  nmem=',nmem(nset),' > nhess=',nhess
  65          stop
  66       endif
  67       Read  (1, '(A)') Line
  68       Read  (1, '(A)') Line
  69       Read  (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
  70       Iorder(nset) = Nint(Dr)
  71       Nfl(nset) = Nint(Fl)
  72       Alambda(nset) = Al(nset)
  73
  74       cMass(nset) = Amass(4,nset)
  75       bMass(nset) = Amass(5,nset)
  76       tMass(nset) = Amass(6,nset)
  77
  78       Read  (1, '(A)') Line
  79       Read  (1, *) NX(nset),  NT(nset), NfMx(nset)
  80
  81       Read  (1, '(A)') Line
  82       Read  (1, *)                                                      &
  83      &  QINI(nset), QMAX(nset), (TV(I,nset), I =0, NT(nset))
  84
  85       Read  (1, '(A)') Line
  86       Read  (1, *) XMIN(nset), (XV(I,nset), I =0, NX(nset))
  87
  88       Do 11 Iq = 0, NT(nset)
  89          TV(Iq,nset) = Log(Log (TV(Iq,nset) /Al(nset)))
  90    11 Continue
  91 !
  92 !                  Since quark = anti-quark for nfl>2 at this stage,
  93 !                  we Read  out only the non-redundent data points
  94 !     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
  95
  96       Nblk = (NX(nset)+1) * (NT(nset)+1)
  97       Npts =  Nblk  * (NfMx(nset)+3)
  98
  99                                              !*** new version: allows nm
 100       do ihess = 0,nmem(nset)
 101
 102         Read  (1, '(A)') Line
 103         Read  (1, '(A)') Line
 104         Read  (1, *, IOSTAT=IRET) (UPD(ihess,I,nset), I=1,Npts)
 105
 106       enddo
 107       return
 108 !
 109       entry CTEQ6alfa(alfas,Qalfa)
 110       alfas = pi*CtLhALPI(Qalfa)
 111       return
 112 !
 113       entry CTEQ6init(Eorder,Q2fit)
 114
 115       return
 116 !
 117       entry CTEQ6pdf(mem)
 118 !        imem = mem
 119         call getnset(iset)
 120         call setnmem(iset,mem)
 121       return
 122 !
 123       END
 124
 125       subroutine LHct6set
 126       Implicit Double Precision (A-H,O-Z)
 127       common /ctq6co/ xlast, qlast, nxsave
 128       nxsave = -1000
 129       xlast = -2.
 130       qlast = -2.
 131       return
 132       END
 133
 134 !=======================================================================
 135       Function CtLhPartonX6 (imem,IPRTN, XX, QQ)
 136 !  Given the parton distribution function in the array U in
 137 !  COMMON / PEVLDT / , this routine interpolates to find
 138 !  the parton distribution at an arbitray point in x and q.
 139 !
 140       Implicit Double Precision (A-H,O-Z)
 141
 142       include 'parmsetup.inc'
 143       Parameter (MXX = 96, MXQ = 20, MXF = 5, nhess = 40)
 144       Parameter (MXQX= MXQ * MXX,   MXPQX = MXQX * (MXF+3))
 145
 146       Common                                                            &
 147      & / CtqPar1nhess /                                                 &
 148      &  Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset),                 &
 149      &  UPD(0:nhess,MXPQX,nmxset)                                       &
 150      & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset)                 &
 151      & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset)
 152
 153       common /ctq6co/ xlast, qlast, nxsave
 154         parameter(nqvec = 4)
 155
 156       Dimension fvec(4), fij(4)
 157       Dimension xvpow(0:mxx)
 158       Data OneP / 1.00001 /
 159                                 !**** choice of interpolation variable
 160       Data xpow / 0.3d0 /
 161       Save xvpow
 162       Data ixprint,iqprint/0,0/
 163       save ixprint,iqprint
 164
 165          save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
 166          save const1 , const2, const3, const4, const5, const6
 167          save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
 168
 169
 170       call getnset(iset)
 171       call getnmem(iset,imem)
 172
 173 ! store the powers used for interpolation on first call...
 174       if(nx(iset).ne. nxsave) then
 175          nxsave = nx(iset)
 176          xvpow(0) = 0.D0
 177          do i = 1, nx(iset)
 178             xvpow(i) = xv(i,iset)**xpow
 179          enddo
 180       endif
 181
 182       X = XX
 183       Q = QQ
 184
 185       if((x.lt.xmin(iset)).or.(x.gt.1.d0)) then
 186         ixprint=ixprint+1
 187         if(ixprint.lt.11) print 98,x
 188         if(ixprint.eq.10) print *,                                      &
 189      &  'more warning messages like the last suppressed.'
 190       endif
 191    98 format(' WARNING:  X=',e12.5,' OUT OF RANGE')
 192       if((q.lt.qini(iset)).or.(q.gt.qmax(iset))) then
 193         iqprint=iqprint+1
 194         if(iqprint.lt.11) print 99,q
 195         if(iqprint.eq.10) print *,                                      &
 196      &  'more warning messages like the last suppressed.'
 197       endif
 198    99 format(' WARNING:  Q=',e12.5,' OUT OF RANGE')
 199
 200 ! skip the initialization in x if same as in the previous call.
 201         if(x .eq. xlast) goto 100
 202         xlast = x
 203
 204 !      -------------    find lower end of interval containing x, i.e.,
 205 !                       get jx such that xv(jx) .le. x .le. xv(jx+1)...
 206       JLx = -1
 207       JU = Nx(iset)+1
 208    11 If (JU-JLx .GT. 1) Then
 209          JM = (JU+JLx) / 2
 210          If (X .Ge. XV(JM,iset)) Then
 211             JLx = JM
 212          Else
 213             JU = JM
 214          Endif
 215          Goto 11
 216       Endif
 217 !                     Ix    0   1   2      Jx  JLx         Nx-2     Nx
 218 !                           |---|---|---|...|---|-x-|---|...|---|---|
 219 !                     x     0  Xmin               x                 1
 220 !
 221       If     (JLx .LE. -1) Then
 222         Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX6 x=',x
 223         Stop
 224       ElseIf (JLx .Eq. 0) Then
 225          Jx = 0
 226       Elseif (JLx .LE. Nx(iset)-2) Then
 227
 228 !                For interior points, keep x in the middle, as shown abo
 229          Jx = JLx - 1
 230       Elseif (JLx.Eq.Nx(iset)-1 .or. x.LT.OneP) Then
 231
 232 !                  We tolerate a slight over-shoot of one (OneP=1.00001)
 233 !              perhaps due to roundoff or whatever, but not more than th
 234 !                                      Keep at least 4 points >= Jx
 235          Jx = JLx - 2
 236       Else
 237         Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX6 x=',x
 238         Stop
 239       Endif
 240 !          ---------- Note: JLx uniquely identifies the x-bin; Jx does n
 241
 242 !                       This is the variable to be interpolated in
 243       ss = x**xpow
 244
 245       If (JLx.Ge.2 .and. JLx.Le.Nx(iset)-2) Then
 246
 247 !     initiation work for "interior bins": store the lattice points in s
 248       svec1 = xvpow(jx)
 249       svec2 = xvpow(jx+1)
 250       svec3 = xvpow(jx+2)
 251       svec4 = xvpow(jx+3)
 252
 253       s12 = svec1 - svec2
 254       s13 = svec1 - svec3
 255       s23 = svec2 - svec3
 256       s24 = svec2 - svec4
 257       s34 = svec3 - svec4
 258
 259       sy2 = ss - svec2
 260       sy3 = ss - svec3
 261
 262 ! constants needed for interpolating in s at fixed t lattice points...
 263       const1 = s13/s23
 264       const2 = s12/s23
 265       const3 = s34/s23
 266       const4 = s24/s23
 267       s1213 = s12 + s13
 268       s2434 = s24 + s34
 269       sdet = s12*s34 - s1213*s2434
 270       tmp = sy2*sy3/sdet
 271       const5 = (s34*sy2-s2434*sy3)*tmp/s12
 272       const6 = (s1213*sy2-s12*sy3)*tmp/s34
 273
 274       EndIf
 275
 276   100      continue
 277
 278 ! skip the initialization in q if same as in the previous call.
 279         if(q .eq. qlast) goto 110
 280         qlast = q
 281
 282       tt = log(log(Q/Al(iset)))
 283
 284 !         --------------Now find lower end of interval containing Q, i.e
 285 !                          get jq such that qv(jq) .le. q .le. qv(jq+1).
 286       JLq = -1
 287       JU = NT(iset)+1
 288    12 If (JU-JLq .GT. 1) Then
 289          JM = (JU+JLq) / 2
 290          If (tt .GE. TV(JM,iset)) Then
 291             JLq = JM
 292          Else
 293             JU = JM
 294          Endif
 295          Goto 12
 296        Endif
 297
 298       If     (JLq .LE. 0) Then
 299          Jq = 0
 300       Elseif (JLq .LE. Nt(iset)-2) Then
 301 !                                  keep q in the middle, as shown above
 302          Jq = JLq - 1
 303       Else
 304 !                         JLq .GE. Nt-1 case:  Keep at least 4 points >=
 305         Jq = Nt(iset) - 3
 306
 307       Endif
 308 !                                   This is the interpolation variable i
 309
 310       If (JLq.GE.1 .and. JLq.LE.Nt(iset)-2) Then
 311 !                                        store the lattice points in t..
 312       tvec1 = Tv(jq,iset)
 313       tvec2 = Tv(jq+1,iset)
 314       tvec3 = Tv(jq+2,iset)
 315       tvec4 = Tv(jq+3,iset)
 316
 317       t12 = tvec1 - tvec2
 318       t13 = tvec1 - tvec3
 319       t23 = tvec2 - tvec3
 320       t24 = tvec2 - tvec4
 321       t34 = tvec3 - tvec4
 322
 323       ty2 = tt - tvec2
 324       ty3 = tt - tvec3
 325
 326       tmp1 = t12 + t13
 327       tmp2 = t24 + t34
 328
 329       tdet = t12*t34 - tmp1*tmp2
 330
 331       EndIf
 332
 333   110      continue
 334
 335 ! get the pdf function values at the lattice points...
 336 ! In this code, we store 8 flavors: u,ubar,d,dbar,s=sbar,c=cbar,b=bbar,g
 337 ! hence Iprtn=3,4,5 (s,c,b) are obtained from -3,-4,-5 (sbar,cbar,bbar)
 338
 339       If (Iprtn .GE. 3) Then
 340          Ip = - Iprtn
 341       Else
 342          Ip = Iprtn
 343       EndIf
 344       jtmp = ((Ip + NfMx(iset))*(NT(iset)+1)+(jq-1))*(NX(iset)+1)+jx+1
 345
 346       Do it = 1, nqvec
 347
 348          J1  = jtmp + it*(NX(iset)+1)
 349
 350        If (Jx .Eq. 0) Then
 351 !                          For the first 4 x points, interpolate x^2*f(x
 352 !                           This applies to the two lowest bins JLx = 0,
 353 !            We cannot put the JLx.eq.1 bin into the "interior" section
 354 !                           (as we do for q), since Upd(J1) is undefined
 355          fij(1) = 0
 356          fij(2) = Upd(imem,J1+1,iset) * XV(1,iset)**2
 357          fij(3) = Upd(imem,J1+2,iset) * XV(2,iset)**2
 358          fij(4) = Upd(imem,J1+3,iset) * XV(3,iset)**2
 359 !
 360 !                 Use CtLhPolint which allows x to be anywhere w.r.t. th
 361
 362          Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx)
 363
 364          If (x .GT. 0D0)  Fvec(it) =  Fx / x**2
 365 !                                              Pdf is undefined for x.eq
 366        ElseIf  (JLx .Eq. Nx(iset)-1) Then
 367 !                                                This is the highest x b
 368
 369 !** fix allow 4 consecutive elements with imem... mrw 19.9.2005
 370         fij(1) = Upd(imem,j1,iset)
 371         fij(2) = Upd(imem,j1+1,iset)
 372         fij(3) = Upd(imem,j1+2,iset)
 373         fij(4) = Upd(imem,j1+3,iset)
 374         Call CtLhPolint4 (XVpow(Nx(iset)-3), Fij(1), 4, ss, Fx, Dfx)
 375
 376         Fvec(it) = Fx
 377
 378        Else
 379
 380 ! for all interior points, use Jon's in-line function
 381 ! This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
 382 ! (This is cubic spline interpolation, as used by cteq; it was
 383 ! changed to polint in previous Durham releases (jcp).)
 384          sf2 = Upd(imem,J1+1,iset)
 385          sf3 = Upd(imem,J1+2,iset)
 386
 387          Fvec(it) = (const5*(Upd(imem,J1,iset)                          &
 388      &                             - sf2*const1 + sf3*const2)           &
 389      &             + const6*(Upd(imem,J1+3,iset)                        &
 390      &                             + sf2*const3 - sf3*const4)           &
 391      &             + sf2*sy3 - sf3*sy2) / s23
 392
 393        Endif
 394
 395       enddo
 396 !                                   We now have the four values Fvec(1:4
 397 !     interpolate in t...
 398
 399       If (JLq .LE. 0) Then
 400 !                         1st Q-bin, as well as extrapolation to lower Q
 401         Call CtLhPolint4(TV(0,iset), Fvec(1), 4, tt, ff, Dfq)
 402
 403       ElseIf (JLq .GE. Nt(iset)-1) Then
 404 !                         Last Q-bin, as well as extrapolation to higher
 405         Call CtLhPolint4(TV(Nt(iset)-3,iset), Fvec(1), 4, tt, ff, Dfq)
 406       Else
 407 !                         Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
 408 !       which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
 409 !                         the full range QV(0:Nt)  (in contrast to XV)
 410         tf2 = fvec(2)
 411         tf3 = fvec(3)
 412
 413         g1 = ( tf2*t13 - tf3*t12) / t23
 414         g4 = (-tf2*t34 + tf3*t24) / t23
 415
 416         h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12                      &
 417      &    +  (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
 418
 419         ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
 420       EndIf
 421
 422       CtLhPartonX6 = ff
 423
 424       Return
 425 !                                       ********************
 426       END
 427 !=======================================================================
 428       Function CtLhCtq6Pdf (imem,Iparton, X, Q)
 429       Implicit Double Precision (A-H,O-Z)
 430       include 'parmsetup.inc'
 431       Logical Warn
 432       Common                                                            &
 433      & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset)                 &
 434      & / QCDtable /  Alambda(nmxset), Nfl(nmxset), Iorder(nmxset)
 435
 436       Data Warn /.true./
 437       save Warn
 438
 439       call getnset(iset)
 440
 441       If (X .lt. 0D0 .or. X .gt. 1D0) Then
 442         Print *, 'X out of range in CtLhCtq6Pdf: ', X
 443         Stop
 444       Endif
 445       If (Q .lt. Alambda(iset)) Then
 446         Print *, 'Q out of range in CtLhCtq6Pdf: ', Q
 447         Stop
 448       Endif
 449
 450 ! added to force pdf = 0.0 at x=1.0 exactly - mrw
 451       if(x .eq. 1.0d0) then
 452           CtLhCtq6Pdf = 0.0d0
 453           return
 454       endif
 455 !
 456       If ((Iparton .lt. -NfMx(iset) .or. Iparton .gt. NfMx(iset))) Then
 457          If (Warn) Then
 458 !        put a warning for calling extra flavor.
 459              Warn = .false.
 460              Print *, 'Warning: Iparton out of range in CtLhCtq6Pdf: '  &
 461      &              , Iparton
 462          Endif
 463          CtLhCtq6Pdf = 0D0
 464          Return
 465       Endif
 466
 467       CtLhCtq6Pdf = CtLhPartonX6 (imem,Iparton, X, Q)
 468       if(CtLhCtq6Pdf.lt.0.D0)  CtLhCtq6Pdf = 0.D0
 469
 470       Return
 471
 472 !                             ********************
 473       END
 474       SUBROUTINE CtLhPOLINT4 (XA,YA,N,X,Y,DY)
 475 ! fast version of polint, valid only for N=4
 476 ! Have explicitly unrolled the loops.
 477
 478       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 479
 480       PARAMETER (NMAX=4)
 481       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 482
 483         if(n .ne. 4) then
 484            print *,'fatal CtLhPolint4 call',n
 485            stop
 486         endif
 487
 488       NS=1
 489       DIF=ABS(X-XA(1))
 490
 491         DIFT=ABS(X-XA(1))
 492         IF (DIFT.LT.DIF) THEN
 493           NS=1
 494           DIF=DIFT
 495         ENDIF
 496         C(1)=YA(1)
 497         D(1)=YA(1)
 498
 499         DIFT=ABS(X-XA(2))
 500         IF (DIFT.LT.DIF) THEN
 501           NS=2
 502           DIF=DIFT
 503         ENDIF
 504         C(2)=YA(2)
 505         D(2)=YA(2)
 506
 507         DIFT=ABS(X-XA(3))
 508         IF (DIFT.LT.DIF) THEN
 509           NS=3
 510           DIF=DIFT
 511         ENDIF
 512         C(3)=YA(3)
 513         D(3)=YA(3)
 514
 515         DIFT=ABS(X-XA(4))
 516         IF (DIFT.LT.DIF) THEN
 517           NS=4
 518           DIF=DIFT
 519         ENDIF
 520         C(4)=YA(4)
 521         D(4)=YA(4)
 522
 523
 524       Y=YA(NS)
 525       NS=NS-1
 526
 527
 528           HO=XA(1)-X
 529           HP=XA(2)-X
 530           W=C(2)-D(1)
 531           DEN=W/(HO-HP)
 532           D(1)=HP*DEN
 533           C(1)=HO*DEN
 534
 535
 536           HO=XA(2)-X
 537           HP=XA(3)-X
 538           W=C(3)-D(2)
 539           DEN=W/(HO-HP)
 540           D(2)=HP*DEN
 541           C(2)=HO*DEN
 542
 543
 544           HO=XA(3)-X
 545           HP=XA(4)-X
 546           W=C(4)-D(3)
 547           DEN=W/(HO-HP)
 548           D(3)=HP*DEN
 549           C(3)=HO*DEN
 550
 551         IF (2*NS.LT.3)THEN
 552           DY=C(NS+1)
 553         ELSE
 554           DY=D(NS)
 555           NS=NS-1
 556         ENDIF
 557         Y=Y+DY
 558
 559           HO=XA(1)-X
 560           HP=XA(3)-X
 561           W=C(2)-D(1)
 562           DEN=W/(HO-HP)
 563           D(1)=HP*DEN
 564           C(1)=HO*DEN
 565
 566           HO=XA(2)-X
 567           HP=XA(4)-X
 568           W=C(3)-D(2)
 569           DEN=W/(HO-HP)
 570           D(2)=HP*DEN
 571           C(2)=HO*DEN
 572
 573         IF (2*NS.LT.2)THEN
 574           DY=C(NS+1)
 575         ELSE
 576           DY=D(NS)
 577           NS=NS-1
 578         ENDIF
 579         Y=Y+DY
 580
 581           HO=XA(1)-X
 582           HP=XA(4)-X
 583           W=C(2)-D(1)
 584           DEN=W/(HO-HP)
 585           D(1)=HP*DEN
 586           C(1)=HO*DEN
 587
 588         IF (2*NS.LT.4-3)THEN
 589           DY=C(NS+1)
 590         ELSE
 591           DY=D(NS)
 592           NS=NS-1
 593         ENDIF
 594         Y=Y+DY
 595
 596       RETURN
 597       END
 598       SUBROUTINE CTLHPOLINT3 (XA,YA,N,X,Y,DY)
 599 ! fast version of polint, valid only for N=3
 600 ! Have explicitly unrolled the loops.
 601       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 602       PARAMETER (NMAX=3)
 603       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 604         if(n .ne. 3) then
 605            print *,'fatal CtLhPolint3 call',n
 606            stop
 607         endif
 608       NS=1
 609       DIF=ABS(X-XA(1))
 610         DIFT=ABS(X-XA(1))
 611         IF (DIFT.LT.DIF) THEN
 612           NS=1
 613           DIF=DIFT
 614         ENDIF
 615         C(1)=YA(1)
 616         D(1)=YA(1)
 617         DIFT=ABS(X-XA(2))
 618         IF (DIFT.LT.DIF) THEN
 619           NS=2
 620           DIF=DIFT
 621         ENDIF
 622         C(2)=YA(2)
 623         D(2)=YA(2)
 624         DIFT=ABS(X-XA(3))
 625         IF (DIFT.LT.DIF) THEN
 626           NS=3
 627           DIF=DIFT
 628         ENDIF
 629         C(3)=YA(3)
 630         D(3)=YA(3)
 631       Y=YA(NS)
 632       NS=NS-1
 633           HO=XA(1)-X
 634           HP=XA(2)-X
 635           W=C(2)-D(1)
 636           DEN=W/(HO-HP)
 637           D(1)=HP*DEN
 638           C(1)=HO*DEN
 639           HO=XA(2)-X
 640           HP=XA(3)-X
 641           W=C(3)-D(2)
 642           DEN=W/(HO-HP)
 643           D(2)=HP*DEN
 644           C(2)=HO*DEN
 645         IF (2*NS.LT.2)THEN
 646           DY=C(NS+1)
 647         ELSE
 648           DY=D(NS)
 649           NS=NS-1
 650         ENDIF
 651         Y=Y+DY
 652           HO=XA(1)-X
 653           HP=XA(3)-X
 654           W=C(2)-D(1)
 655           DEN=W/(HO-HP)
 656           D(1)=HP*DEN
 657           C(1)=HO*DEN
 658         IF (2*NS.LT.1)THEN
 659           DY=C(NS+1)
 660         ELSE
 661           DY=D(NS)
 662           NS=NS-1
 663         ENDIF
 664         Y=Y+DY
 665       RETURN
 666       END