end-of-line normalization
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.5.1 / src / eps09.f
index eadc639..eba1cd2 100644 (file)
-!****************************************************************************\r
-!\r
-!           EPS09.f\r
-!\r
-! An interface for the scale dependent nuclear modifications\r
-!       R_f^A(x,Q) = f_A(x,Q)/f_p(x,Q) \r
-! where f_A is the distribution of the parton flavour f for a PROTON in a\r
-! nucleus A, and f_p is the corresponding parton distribution in the \r
-! free proton.\r
-!  \r
-! When using this interface, please refer to:\r
-!  \r
-! K.J. Eskola, H. Paukkunen and C.A. Salgado,\r
-! "EPS09 - a New Generation of NLO and LO Nuclear Parton Distribution Functions,"\r
-! Published as JHEP04(2009) 065.\r
-! Eprint: arXiv:0902.4154 [hep-ph].\r
-!\r
-! Questions & comments to:\r
-!   hannu.paukkunen@phys.jyu.fi\r
-!   kari.eskola@phys.jyu.fi\r
-!   carlos.salgado@usc.es\r
-! \r
-! ***************************************************************************\r
-! Instructions:\r
-!\r
-! For given input values of\r
-!\r
-!     order: 1=LO, 2=NLO   ; integer\r
-!     pset : 1...31        ; integer\r
-!            1     = central fit\r
-!            2,3   = error sets S{+1}, S{-1}\r
-!            4,5   = error sets S{+2}, S{-2}\r
-!            ...   ...\r
-!            30,31 = error sets {S+15}, {S-15}\r
-!     A    : atomic number ; integer\r
-!     x    : Bjorken-x     ; double precision\r
-!     Q    : scale in GeV  ; double precision\r
-!\r
-! the command \r
-!\r
-!   Call EPS09(order, pset, A, x, Q, ruv, rdv, ru, rd, rs, rc, rb, rg)\r
-!\r
-! returns the bound proton nuclear corrections R_f^A(x,Q)\r
-! (in double precision) for\r
-!   \r
-!   ruv = up valence\r
-!   rdv = down valence\r
-!   ru  = up sea\r
-!   rd  = down sea\r
-!   rs  = strange\r
-!   rc  = charm\r
-!   rb  = bottom\r
-!   rg  = gluons\r
-!\r
-! The nuclear corrections for bound neutrons can be obtained\r
-! by the isospin symmetry, e.g. the total up quark distribution\r
-! per nucleon in a nucleus A with Z protons is\r
-!\r
-!  u_A(x,Q) =    Z/A * [ruv*uV_p(x,Q) + ru*uSea_p(x,Q)] +\r
-!            (A-Z)/A * [rdv*dV_p(x,Q) + rd*dSea_p(x,Q)]\r
-!\r
-! Note that the parametrization should only be applied at the\r
-! kinematical domain\r
-!\r
-!             1e-6 <= x <= 1\r
-!              1.3 <= Q <= 1000 GeV.\r
-!\r
-! No warning message is displayed if these limits are\r
-! exceeded, and outside these boundaries the modifications\r
-! are frozen to the boundary values, i.e\r
-!\r
-!   for Q > 1000, the modifications at Q=1000 are returned,\r
-!   for Q < 1.3,  the modifications at Q=1.3 are returned,\r
-!   for x > 1,    the modifications at x=1 are returned\r
-!   for x < 1e-6, the modifications at x=1e-6 are returned,\r
-!\r
-! The data used by the program for required order\r
-! and atomic number A, are stored in separate files\r
-!\r
-!   LO : EPS09LOR_A\r
-!   NLO: EPS09NLOR_A\r
-!\r
-! which must be located in the working directory.\r
-!\r
-! The error bands for absolute cross-sections and for\r
-! their nuclear ratios should be computed as explained\r
-! in Secs. 2.5 and 4 of arXiv:0902.4154 [hep-ph]. For\r
-! the absolute cross sections, both the errors in the\r
-! free-proton PDFs f_p(x,Q) and the errors in\r
-! the modifications R_f^A(x,Q) should be accounted for.\r
-! For the nuclear ratios, it is sufficient to account only\r
-! for the errors in the modifications R_f^A(x,Q).\r
-!\r
-! *********************************************************\r
-! *********************************************************\r
-\r
-\r
-\r
-      Subroutine EPS09(order, pset, AAA, xxx, QQQ,ruv, rdv, ru, rd, rs, rc, rb, rg)\r
-\r
-      Implicit none\r
-      Double precision :: ruv, rdv, ru, rd, rs, rc, rb, rg, QQQ, xxx\r
-      Double precision :: LSTEP, x, Q, Q2, allvalues(1:31,1:8,0:50,0:50)\r
-      Double precision :: x_i=0.000001, arg(4), fu(4), res, fg(3)\r
-      Double precision :: result(9), dummy\r
-      Double precision :: realQ, Q2min=1.69, Q2max=1000000.0, Qsteps=50.0\r
-      Double precision :: n_x, zero=0.0\r
-\r
-      Character (Len=50) filenimi\r
-\r
-      Integer :: xlinsteps=25, xlogsteps=25, startline, lineno\r
-      Integer :: k, p, t, Qpoint, xpoint, pset, iovar\r
-      Integer :: setnumber,j, A, openchannel, order, AAA\r
-      Integer :: psetlast = -10, Alast = -10, orderlast = -10\r
\r
-      character*512 dirpath,setpath\r
-\r
-      save Alast\r
-      save psetlast\r
-      save orderlast\r
-      save allvalues\r
-\r
-! *********************************************\r
-! Stop if the set specifications are wrong ones\r
-! *********************************************\r
-\r
-      If (order .NE. 1 .and. order .NE. 2) then\r
-      Write(*,*) 'Wrong order!'\r
-      Write(*,*) 'LO : order = 1'\r
-      Write(*,*) 'NLO: order = 2'\r
-      Stop\r
-      End If\r
-\r
-      If (pset .LT. 1 .or. pset .GT. 31) then\r
-      Write(*,*) 'Wrong set!'\r
-      Write(*,*) 'Central set: pset = 1'\r
-      Write(*,*) 'Error sets : pset = 2...31'\r
-      Stop\r
-      End If\r
-\r
-! ********************************\r
-! Make sure not to change any\r
-! specifications given by the user\r
-! ********************************\r
-\r
-      A  = AAA\r
-      x  = xxx\r
-      Q  = QQQ\r
-      Q2 = Q*Q \r
-\r
-! *******************************\r
-! Freeze x if it's < 10E-6 or > 1\r
-! *******************************\r
-\r
-      If (x .LT. x_i) Then\r
-      x = x_i\r
-      End If\r
-      If (x .GT. 1) Then\r
-      x = 1.0\r
-      End If\r
-\r
-! ************************************\r
-! Freeze Q^2 if it's < 1.69 or > 10E+6\r
-! ************************************\r
-\r
-      If (Q2 .LT. Q2min) Then\r
-      Q2 = Q2min\r
-      End If\r
-      If (Q2 .GT. Q2max) Then\r
-      Q2 = Q2max\r
-      End If\r
-\r
-! If the set specifications have been changed, read the tables again\r
-\r
-      If (A .NE. Alast .or. order .NE. orderlast) Then\r
-\r
-!      Write(*,*) 'Set changed!'\r
-\r
-! Read the table1\r
-\r
-      If (order .EQ. 1) then\r
-\r
-        If (A < 10) Then\r
-        Write(filenimi,'("EPS09LOR_", I1)'), A\r
-        Else If (A < 100) Then\r
-        Write(filenimi,'("EPS09LOR_", I2)'), A\r
-        Else If (A < 1000) Then\r
-        Write(filenimi,'("EPS09LOR_", I3)'), A\r
-        End If\r
-\r
-      Else\r
-\r
-        If (A < 10) Then\r
-        Write(filenimi,'("EPS09NLOR_", I1)'), A\r
-        Else If (A < 100) Then\r
-        Write(filenimi,'("EPS09NLOR_", I2)'), A\r
-        Else If (A < 1000) Then\r
-        Write(filenimi,'("EPS09NLOR_", I3)'), A\r
-        End If\r
-\r
-      End If\r
-\r
-      call getdirpath(dirpath) \r
-      setpath=dirpath(:LEN_TRIM(dirpath))//"/"//filenimi(:LEN_TRIM(filenimi))//".LHgrid"\r
-      print *,setpath(:LEN_TRIM(setpath))\r
-      OPEN (11, file = setpath(:LEN_TRIM(setpath)), status='OLD', IOSTAT=iovar)\r
-\r
-      If (iovar .NE. 0) Then\r
-      Write(*,*) 'Missing file: ',filenimi\r
-      stop\r
-      End If\r
-\r
-      Do setnumber = 1, 31\r
-\r
-      Do k = 0,50\r
-\r
-      Read(11,*) dummy\r
-\r
-      Do t = 0,50\r
-\r
-      Read(11,*) (allvalues(setnumber,p,k,t), p=1,8)\r
-\r
-      End Do\r
-      End Do\r
-\r
-      End Do\r
-\r
-      Close(11)\r
-\r
-      psetlast  = pset\r
-      Alast     = A\r
-      orderlast = order\r
-\r
-      End If\r
-\r
-! Find out the position in the loglog Q^2-grid\r
-\r
-      realQ  = Qsteps * (log(log(Q2)/log(Q2min)))/(log(log(Q2max)/log(Q2min)))\r
-      Qpoint = Aint(realQ)\r
-\r
-      If (Qpoint .LE. 0) Then\r
-         Qpoint = 1\r
-      End If\r
-      If (Qpoint .GE. Anint(Qsteps)-1) Then\r
-         Qpoint = Anint(Qsteps)-1\r
-      End If\r
-\r
-      LSTEP = (1.0/(xlogsteps)) * LOG(0.1/x_i)\r
-\r
-! *********************\r
-! Interpolate the grids \r
-! *********************\r
-\r
-      Do t=1,8\r
-\r
-! Find the position in the x-grid\r
-\r
-      If (x .LE. 0.1) then\r
-         n_x  = ((1.0/LSTEP) * Log(x/x_i))\r
-       xpoint = Aint(n_x)\r
-      Else\r
-       n_x    = ((x-0.1)*xlinsteps/(1.0-0.1) + xlogsteps)\r
-       xpoint = Aint(n_x)\r
-      End If\r
-\r
-      If (xpoint .LE. 0) Then\r
-        xpoint = 1\r
-      End If\r
-\r
-      If (t .EQ. 1 .or. t .EQ. 2) Then\r
-         If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then\r
-         xpoint = (xlinsteps+xlogsteps)-4\r
-         End If\r
-      End If\r
-\r
-      If (t .EQ. 3 .or. t .EQ. 4 .or. t .EQ. 5 .or. t .EQ. 6 .or. t .EQ. 7) Then\r
-        If (xpoint .GE. (xlinsteps+xlogsteps)-7) Then\r
-        xpoint = (xlinsteps+xlogsteps)-7\r
-        End If\r
-      End If\r
-       \r
-      If (t .EQ. 8) Then\r
-        If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then\r
-        xpoint = (xlinsteps+xlogsteps)-4\r
-        End If\r
-      End If\r
-\r
-      Do k = 1, 4\r
-      If (xpoint-2+k .LT. xlogsteps) Then\r
-      arg(k) = (x_i) * exp(LSTEP * (xpoint-2+k))\r
-      Else\r
-      arg(k) = 0.1 + (xpoint-2+k-xlogsteps) * (1-0.1)/xlinsteps\r
-      End If\r
-      End Do\r
\r
-      Do j=1,3\r
-\r
-      fu(1) = allvalues(pset,t,Qpoint-2+j,xpoint-1)\r
-      fu(2) = allvalues(pset,t,Qpoint-2+j,xpoint)\r
-      fu(3) = allvalues(pset,t,Qpoint-2+j,xpoint+1)\r
-      fu(4) = allvalues(pset,t,Qpoint-2+j,xpoint+2)\r
-      Call luovi(fu,arg,4,x,res)\r
-      fg(j) = res\r
-\r
-! *****************************************\r
-! *****************************************\r
-\r
-      End Do\r
-\r
-      arg(1) = Qpoint-1\r
-      arg(2) = Qpoint\r
-      arg(3) = Qpoint+1\r
-\r
-      Call luovi(fg,arg,3,realQ,res)\r
-  \r
-      result(t) = res\r
-\r
-      End Do\r
-\r
-      ruv = max(result(1),zero)\r
-      rdv = max(result(2),zero)\r
-      ru  = max(result(3),zero)\r
-      rd  = max(result(4),zero)\r
-      rs  = max(result(5),zero)\r
-      rc  = max(result(6),zero)\r
-      rb  = max(result(7),zero)\r
-      rg  = max(result(8),zero)\r
-\r
-200   Continue\r
-\r
-      End Subroutine EPS09\r
-\r
-! ********************************\r
-! Modified version of Cern Library\r
-! interpolation routine E100\r
-! ********************************\r
-\r
-      SUBROUTINE luovi(F,ARG,MMM,Z,SUM)\r
-\r
-      Implicit none\r
-      INTEGER :: MMM\r
-      Double precision  :: F(MMM), ARG(MMM), COF(MMM), SUM, Z\r
-      INTEGER :: M, MM, I, J, JNDEX, INDEX\r
-\r
-      MM = MIN(MMM, 20)\r
-      M = MM - 1\r
-      DO 1780 I= 1, MM\r
-      COF(I) = F(I)\r
- 1780 Continue\r
-      DO 1800 I= 1, M\r
-      DO 1790 J= I, M\r
-         JNDEX = MM - J\r
-         INDEX = JNDEX + I\r
-         COF(INDEX) = (COF(INDEX)-COF(INDEX-1))/(ARG(INDEX)-ARG(JNDEX))\r
- 1790 CONTINUE\r
- 1800 CONTINUE\r
-      SUM = COF(MM)\r
-      DO 1810 I= 1, M\r
-         INDEX = MM - I\r
-         SUM = (Z-ARG(INDEX))*SUM + COF(INDEX)\r
- 1810 CONTINUE\r
-\r
-      End SUBROUTINE luovi\r
-\r
+!****************************************************************************
+!
+!           EPS09.f
+!
+! An interface for the scale dependent nuclear modifications
+!       R_f^A(x,Q) = f_A(x,Q)/f_p(x,Q) 
+! where f_A is the distribution of the parton flavour f for a PROTON in a
+! nucleus A, and f_p is the corresponding parton distribution in the 
+! free proton.
+!  
+! When using this interface, please refer to:
+!  
+! K.J. Eskola, H. Paukkunen and C.A. Salgado,
+! "EPS09 - a New Generation of NLO and LO Nuclear Parton Distribution Functions,"
+! Published as JHEP04(2009) 065.
+! Eprint: arXiv:0902.4154 [hep-ph].
+!
+! Questions & comments to:
+!   hannu.paukkunen@phys.jyu.fi
+!   kari.eskola@phys.jyu.fi
+!   carlos.salgado@usc.es
+! 
+! ***************************************************************************
+! Instructions:
+!
+! For given input values of
+!
+!     order: 1=LO, 2=NLO   ; integer
+!     pset : 1...31        ; integer
+!            1     = central fit
+!            2,3   = error sets S{+1}, S{-1}
+!            4,5   = error sets S{+2}, S{-2}
+!            ...   ...
+!            30,31 = error sets {S+15}, {S-15}
+!     A    : atomic number ; integer
+!     x    : Bjorken-x     ; double precision
+!     Q    : scale in GeV  ; double precision
+!
+! the command 
+!
+!   Call EPS09(order, pset, A, x, Q, ruv, rdv, ru, rd, rs, rc, rb, rg)
+!
+! returns the bound proton nuclear corrections R_f^A(x,Q)
+! (in double precision) for
+!   
+!   ruv = up valence
+!   rdv = down valence
+!   ru  = up sea
+!   rd  = down sea
+!   rs  = strange
+!   rc  = charm
+!   rb  = bottom
+!   rg  = gluons
+!
+! The nuclear corrections for bound neutrons can be obtained
+! by the isospin symmetry, e.g. the total up quark distribution
+! per nucleon in a nucleus A with Z protons is
+!
+!  u_A(x,Q) =    Z/A * [ruv*uV_p(x,Q) + ru*uSea_p(x,Q)] +
+!            (A-Z)/A * [rdv*dV_p(x,Q) + rd*dSea_p(x,Q)]
+!
+! Note that the parametrization should only be applied at the
+! kinematical domain
+!
+!             1e-6 <= x <= 1
+!              1.3 <= Q <= 1000 GeV.
+!
+! No warning message is displayed if these limits are
+! exceeded, and outside these boundaries the modifications
+! are frozen to the boundary values, i.e
+!
+!   for Q > 1000, the modifications at Q=1000 are returned,
+!   for Q < 1.3,  the modifications at Q=1.3 are returned,
+!   for x > 1,    the modifications at x=1 are returned
+!   for x < 1e-6, the modifications at x=1e-6 are returned,
+!
+! The data used by the program for required order
+! and atomic number A, are stored in separate files
+!
+!   LO : EPS09LOR_A
+!   NLO: EPS09NLOR_A
+!
+! which must be located in the working directory.
+!
+! The error bands for absolute cross-sections and for
+! their nuclear ratios should be computed as explained
+! in Secs. 2.5 and 4 of arXiv:0902.4154 [hep-ph]. For
+! the absolute cross sections, both the errors in the
+! free-proton PDFs f_p(x,Q) and the errors in
+! the modifications R_f^A(x,Q) should be accounted for.
+! For the nuclear ratios, it is sufficient to account only
+! for the errors in the modifications R_f^A(x,Q).
+!
+! *********************************************************
+! *********************************************************
+
+
+
+      Subroutine EPS09(order, pset, AAA, xxx, QQQ,ruv, rdv, ru, rd, rs, rc, rb, rg)
+
+      Implicit none
+      Double precision :: ruv, rdv, ru, rd, rs, rc, rb, rg, QQQ, xxx
+      Double precision :: LSTEP, x, Q, Q2, allvalues(1:31,1:8,0:50,0:50)
+      Double precision :: x_i=0.000001, arg(4), fu(4), res, fg(3)
+      Double precision :: result(9), dummy
+      Double precision :: realQ, Q2min=1.69, Q2max=1000000.0, Qsteps=50.0
+      Double precision :: n_x, zero=0.0
+
+      Character (Len=50) filenimi
+
+      Integer :: xlinsteps=25, xlogsteps=25, startline, lineno
+      Integer :: k, p, t, Qpoint, xpoint, pset, iovar
+      Integer :: setnumber,j, A, openchannel, order, AAA
+      Integer :: psetlast = -10, Alast = -10, orderlast = -10
+      character*512 dirpath,setpath
+
+      save Alast
+      save psetlast
+      save orderlast
+      save allvalues
+
+! *********************************************
+! Stop if the set specifications are wrong ones
+! *********************************************
+
+      If (order .NE. 1 .and. order .NE. 2) then
+      Write(*,*) 'Wrong order!'
+      Write(*,*) 'LO : order = 1'
+      Write(*,*) 'NLO: order = 2'
+      Stop
+      End If
+
+      If (pset .LT. 1 .or. pset .GT. 31) then
+      Write(*,*) 'Wrong set!'
+      Write(*,*) 'Central set: pset = 1'
+      Write(*,*) 'Error sets : pset = 2...31'
+      Stop
+      End If
+
+! ********************************
+! Make sure not to change any
+! specifications given by the user
+! ********************************
+
+      A  = AAA
+      x  = xxx
+      Q  = QQQ
+      Q2 = Q*Q 
+
+! *******************************
+! Freeze x if it's < 10E-6 or > 1
+! *******************************
+
+      If (x .LT. x_i) Then
+      x = x_i
+      End If
+      If (x .GT. 1) Then
+      x = 1.0
+      End If
+
+! ************************************
+! Freeze Q^2 if it's < 1.69 or > 10E+6
+! ************************************
+
+      If (Q2 .LT. Q2min) Then
+      Q2 = Q2min
+      End If
+      If (Q2 .GT. Q2max) Then
+      Q2 = Q2max
+      End If
+
+! If the set specifications have been changed, read the tables again
+
+      If (A .NE. Alast .or. order .NE. orderlast) Then
+
+!      Write(*,*) 'Set changed!'
+
+! Read the table1
+
+      If (order .EQ. 1) then
+
+        If (A < 10) Then
+        Write(filenimi,'("EPS09LOR_", I1)'), A
+        Else If (A < 100) Then
+        Write(filenimi,'("EPS09LOR_", I2)'), A
+        Else If (A < 1000) Then
+        Write(filenimi,'("EPS09LOR_", I3)'), A
+        End If
+
+      Else
+
+        If (A < 10) Then
+        Write(filenimi,'("EPS09NLOR_", I1)'), A
+        Else If (A < 100) Then
+        Write(filenimi,'("EPS09NLOR_", I2)'), A
+        Else If (A < 1000) Then
+        Write(filenimi,'("EPS09NLOR_", I3)'), A
+        End If
+
+      End If
+
+      call getdirpath(dirpath) 
+      setpath=dirpath(:LEN_TRIM(dirpath))//"/"//filenimi(:LEN_TRIM(filenimi))//".LHgrid"
+      print *,setpath(:LEN_TRIM(setpath))
+      OPEN (11, file = setpath(:LEN_TRIM(setpath)), status='OLD', IOSTAT=iovar)
+
+      If (iovar .NE. 0) Then
+      Write(*,*) 'Missing file: ',filenimi
+      stop
+      End If
+
+      Do setnumber = 1, 31
+
+      Do k = 0,50
+
+      Read(11,*) dummy
+
+      Do t = 0,50
+
+      Read(11,*) (allvalues(setnumber,p,k,t), p=1,8)
+
+      End Do
+      End Do
+
+      End Do
+
+      Close(11)
+
+      psetlast  = pset
+      Alast     = A
+      orderlast = order
+
+      End If
+
+! Find out the position in the loglog Q^2-grid
+
+      realQ  = Qsteps * (log(log(Q2)/log(Q2min)))/(log(log(Q2max)/log(Q2min)))
+      Qpoint = Aint(realQ)
+
+      If (Qpoint .LE. 0) Then
+         Qpoint = 1
+      End If
+      If (Qpoint .GE. Anint(Qsteps)-1) Then
+         Qpoint = Anint(Qsteps)-1
+      End If
+
+      LSTEP = (1.0/(xlogsteps)) * LOG(0.1/x_i)
+
+! *********************
+! Interpolate the grids 
+! *********************
+
+      Do t=1,8
+
+! Find the position in the x-grid
+
+      If (x .LE. 0.1) then
+         n_x  = ((1.0/LSTEP) * Log(x/x_i))
+       xpoint = Aint(n_x)
+      Else
+       n_x    = ((x-0.1)*xlinsteps/(1.0-0.1) + xlogsteps)
+       xpoint = Aint(n_x)
+      End If
+
+      If (xpoint .LE. 0) Then
+        xpoint = 1
+      End If
+
+      If (t .EQ. 1 .or. t .EQ. 2) Then
+         If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then
+         xpoint = (xlinsteps+xlogsteps)-4
+         End If
+      End If
+
+      If (t .EQ. 3 .or. t .EQ. 4 .or. t .EQ. 5 .or. t .EQ. 6 .or. t .EQ. 7) Then
+        If (xpoint .GE. (xlinsteps+xlogsteps)-7) Then
+        xpoint = (xlinsteps+xlogsteps)-7
+        End If
+      End If
+       
+      If (t .EQ. 8) Then
+        If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then
+        xpoint = (xlinsteps+xlogsteps)-4
+        End If
+      End If
+
+      Do k = 1, 4
+      If (xpoint-2+k .LT. xlogsteps) Then
+      arg(k) = (x_i) * exp(LSTEP * (xpoint-2+k))
+      Else
+      arg(k) = 0.1 + (xpoint-2+k-xlogsteps) * (1-0.1)/xlinsteps
+      End If
+      End Do
+      Do j=1,3
+
+      fu(1) = allvalues(pset,t,Qpoint-2+j,xpoint-1)
+      fu(2) = allvalues(pset,t,Qpoint-2+j,xpoint)
+      fu(3) = allvalues(pset,t,Qpoint-2+j,xpoint+1)
+      fu(4) = allvalues(pset,t,Qpoint-2+j,xpoint+2)
+      Call luovi(fu,arg,4,x,res)
+      fg(j) = res
+
+! *****************************************
+! *****************************************
+
+      End Do
+
+      arg(1) = Qpoint-1
+      arg(2) = Qpoint
+      arg(3) = Qpoint+1
+
+      Call luovi(fg,arg,3,realQ,res)
+  
+      result(t) = res
+
+      End Do
+
+      ruv = max(result(1),zero)
+      rdv = max(result(2),zero)
+      ru  = max(result(3),zero)
+      rd  = max(result(4),zero)
+      rs  = max(result(5),zero)
+      rc  = max(result(6),zero)
+      rb  = max(result(7),zero)
+      rg  = max(result(8),zero)
+
+200   Continue
+
+      End Subroutine EPS09
+
+! ********************************
+! Modified version of Cern Library
+! interpolation routine E100
+! ********************************
+
+      SUBROUTINE luovi(F,ARG,MMM,Z,SUM)
+
+      Implicit none
+      INTEGER :: MMM
+      Double precision  :: F(MMM), ARG(MMM), COF(MMM), SUM, Z
+      INTEGER :: M, MM, I, J, JNDEX, INDEX
+
+      MM = MIN(MMM, 20)
+      M = MM - 1
+      DO 1780 I= 1, MM
+      COF(I) = F(I)
+ 1780 Continue
+      DO 1800 I= 1, M
+      DO 1790 J= I, M
+         JNDEX = MM - J
+         INDEX = JNDEX + I
+         COF(INDEX) = (COF(INDEX)-COF(INDEX-1))/(ARG(INDEX)-ARG(JNDEX))
+ 1790 CONTINUE
+ 1800 CONTINUE
+      SUM = COF(MM)
+      DO 1810 I= 1, M
+         INDEX = MM - I
+         SUM = (Z-ARG(INDEX))*SUM + COF(INDEX)
+ 1810 CONTINUE
+
+      End SUBROUTINE luovi
+