Merge lp:~maddevelopers/mg5amcnlo/virtual_polyfit_EW into lp:~maddevelopers/mg5amcnlo/3.0.2

[Rikkert Frederix]

Replaced the grids to determine the approximate virtual by a polynomial fit using the polfit.f code (by Shampine, Davenport and Huddleston). It's still based on doing this one dimension at a time, but the polynomial fit is smoother than the grids, since the latter required rather coarse binning and had fixed bins which could not be adjusted to follow the importance sampling.

Also needed to replace the FKSParamReader with a fortran90 module to be able to use it correctly in the mint_module.f90

[marco zaro (marco-zaro)]

Thanks Rikkert,
the code changes look correct, for what I can understand of the polyfit module
I have tried p p > t t~ [QCD QED], results are compatible within the uncertainty (also split-order by split-order). The Polyfit run is generally (much) faster:
If I require 1 per mil accuracy for that process, without polyfit it took 18 minutes on my macbook, while with polyfit less than 8, so great job!!

review: Approve

=== modified file 'Template/NLO/Cards/FKS_params.dat'
--- Template/NLO/Cards/FKS_params.dat	2017-08-04 19:52:15 +0000
+++ Template/NLO/Cards/FKS_params.dat	2019-09-12 11:14:33 +0000
@@ -153,6 +153,16 @@
 0.005d0
 ! Default :: 0.005d0
 !
+! Use a polynomial fit for the virtual corrections. Gets a better accuracy,
+! but uses more memmory and disk space.
+#UsePolyVirtual
+.True.
+! Default :: .True.
+!
+! ==========================================================================
+! Combination of matrix elements
+! ==========================================================================
+!
 ! For fixed order calculations only (ignored for (N)LO+PS runs): This
 ! parameter determines if parton flavour configurations with idential
 ! matrix elements should be separated at the level of the analysis or
@@ -170,11 +180,5 @@
 ! Default :: .false.
 !
 ! ==========================================================================
-! Arbitrary numerical parameters used in the FKS formalism
-! ==========================================================================
-!
-! To be implemented by the FKS authors
-!
-! ==========================================================================
 ! End of FKS_params.dat file
 ! ==========================================================================
 
=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2017-06-23 14:14:03 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2019-09-12 11:14:33 +0000
@@ -4,12 +4,11 @@
 c that calls the OLP and returns the virtual weights. For convenience
 c also the born_wgt is passed to this subroutine.
 c
+      use FKSParams
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
       include 'born_nhel.inc'
-c general MadFKS parameters
-      include 'FKSParams.inc'
       double precision pi, zero,mone
       parameter (pi=3.1415926535897932385d0)
       parameter (zero=0d0)
 
=== modified file 'Template/NLO/SubProcesses/BinothLHA_OLP.f'
--- Template/NLO/SubProcesses/BinothLHA_OLP.f	2015-08-13 09:26:37 +0000
+++ Template/NLO/SubProcesses/BinothLHA_OLP.f	2019-09-12 11:14:33 +0000
@@ -6,11 +6,11 @@
 c
 C************************************************************************
 c
+      use FKSParams
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
       include "Binoth_proc.inc"
-      include "FKSParams.inc"
       double precision pi
       parameter (pi=3.1415926535897932385d0)
       double precision pin(0:3,nexternal-1),p(0:4,nexternal-1)
 
=== removed file 'Template/NLO/SubProcesses/FKSParamReader.f'
--- Template/NLO/SubProcesses/FKSParamReader.f	2017-09-01 10:34:24 +0000
+++ Template/NLO/SubProcesses/FKSParamReader.f	1970-01-01 00:00:00 +0000
@@ -1,252 +0,0 @@
-      subroutine FKSParamReader(filename, printParam, force) 
-
-      implicit none
-
-      logical HasReadOnce, force
-      data HasReadOnce/.FALSE./
-
-
-      character(*) filename
-      CHARACTER*64 buff, buff2, mode
-      include "FKSParams.inc"
-      include "orders.inc"
-      
-      integer i,j
-
-      logical printParam, couldRead, paramPrinted
-      data paramPrinted/.FALSE./
-      couldRead=.False.
-!     Default parameters
-
-      if (HasReadOnce.and..not.force) then
-          goto 901
-      endif
-
-!     Make sure to have default parameters if not set 
-!     in the FKSParams.dat card (if it is an old one for instance)
-      call DefaultFKSParam()
-
-      open(68, file=fileName, err=676, action='READ')
-      do
-         read(68,*,end=999) buff
-         if(index(buff,'#').eq.1) then
-
-           if (buff .eq. '#IRPoleCheckThreshold') then
-             read(68,*,end=999) IRPoleCheckThreshold
-             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
-               stop 'IRPoleCheckThreshold must be >= -1.0d0.'
-             endif 
-
-           elseif (buff .eq. '#PrecisionVirtualAtRunTime') then
-             read(68,*,end=999) PrecisionVirtualAtRunTime
-             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
-               stop 'PrecisionVirtualAtRunTime must be >= -1.0d0.'
-             endif 
-
-           else if (buff .eq. '#NHelForMCoverHels') then
-             read(68,*,end=999) NHelForMCoverHels
-             if (NHelForMCoverHels .lt. -1) then
-               stop 'NHelForMCoverHels must be >= -1.'
-             endif
-           else if (buff .eq. '#QCD^2==') then
-             read(68,*,end=999) QCD_squared_selected
-             if (QCD_squared_selected .lt. -1) then
-               stop 'QCD_squared_selected must be >= -1.'
-             endif 
-           else if (buff .eq. '#QED^2==') then
-             read(68,*,end=999) QED_squared_selected
-             if (QED_squared_selected .lt. -1) then
-               stop 'QED_squared_selected must be >= -1.'
-             endif 
-           else if (buff .eq. '#VirtualFraction') then
-             read(68,*,end=999) Virt_fraction
-             if (Virt_fraction .lt. 0 .or. virt_fraction .gt.1) then
-                stop 'VirtualFraction should be a fraction'/
-     $               /' between 0 and 1'
-             endif 
-           else if (buff .eq. '#MinVirtualFraction') then
-             read(68,*,end=999) Min_Virt_fraction
-             if (min_virt_fraction .lt. 0 .or. min_virt_fraction .gt.1)
-     $            then
-                stop 'VirtualFraction should be a fraction'/
-     $               /' between 0 and 1'
-             endif 
-           else if (buff .eq. '#SeparateFlavourConfigurations') then
-             read(68,*,end=999) separate_flavour_configs
-
-           else if (buff .eq. '#VetoedContributionTypes') then
-             read(68,*,end=999) VetoedContributionTypes(0)
-             if (VetoedContributionTypes(0) .lt. 0.or.
-     &           VetoedContributionTypes(0) .gt. maxContribsSelected) then
-                write(*,*) 'VetoedContributionTypes length should be '/
-     &               /'>= 0 and <=',maxContribsSelected
-                stop 'Format error in FKS_params.dat.'
-             endif
-             read(68,*,end=999) (VetoedContributionTypes(I),I=1,
-     $          VetoedContributionTypes(0))
-             do I=1,VetoedContributionTypes(0)
-               if (VetoedContributionTypes(I).lt.1.or.
-     $             VetoedContributionTypes(I).gt.maxContribType) then
-                write(*,*) 'VetoedContributionTypes must be >=1 and '/
-     &               /'<=',maxContribType
-                stop 'Format error in FKS_params.dat.'
-               endif
-             enddo
-             do I=VetoedContributionTypes(0)+1,maxContribsSelected
-               VetoedContributionTypes(I)=-1
-             enddo
-
-           else if (buff .eq. '#SelectedContributionTypes') then
-             read(68,*,end=999) SelectedContributionTypes(0)
-             if (SelectedContributionTypes(0) .lt. 0 .or.
-     &           SelectedContributionTypes(0) .gt. maxContribsSelected) then
-                write(*,*) 'SelectedContributionTypes length should be '/
-     &               /'>= 0 and <=',maxContribsSelected
-                stop 'Format error in FKS_params.dat.'
-             endif
-             read(68,*,end=999) (SelectedContributionTypes(I),I=1,
-     $          SelectedContributionTypes(0))
-             do I=1,SelectedContributionTypes(0)
-               if (SelectedContributionTypes(I).lt.1.or.
-     $             SelectedContributionTypes(I).gt.maxContribType) then
-                write(*,*) 'SelectedContributionTypes must be >=1 and '/
-     &               /'<=',maxContribType
-                stop 'Format error in FKS_params.dat.'
-               endif
-             enddo
-             do I=SelectedContributionTypes(0)+1,maxContribsSelected
-               SelectedContributionTypes(I)=-1
-             enddo
-
-           else if (buff .eq. '#SelectedCouplingOrders') then
-             read(68,*,end=999) SelectedCouplingOrders(1,0)
-             if (SelectedCouplingOrders(1,0) .lt. 0 .or.
-     &           SelectedCouplingOrders(1,0) .gt. maxCouplingsSelected) then
-                write(*,*) 'SelectedCouplingOrders length should be >='/
-     &               /' 0 and <=',maxCouplingsSelected
-                stop 'Format error in FKS_params.dat.'                
-             endif
-             do j = 2, maxCouplingTypes
-               SelectedCouplingOrders(j,0) = SelectedCouplingOrders(1,0)
-             enddo
-             do j = 1, SelectedCouplingOrders(1,0)
-               read(68,*,end=999) (SelectedCouplingOrders(i,j),i=1,
-     $            nsplitorders)
-               do i=nsplitorders+1,maxCouplingTypes
-                 SelectedCouplingOrders(i,j)=-1
-               enddo
-             enddo
-
-           else
-             write(*,*) 'The parameter name ',buff(2:),
-     &' is not reckognized.'
-             stop 'Format error in FKS_params.dat.'
-           endif
-         endif
-      enddo
-  999 continue
-      couldRead=.True.
-      goto 998      
-
-  676 continue
-      write(*,*) 'ERROR :: MadFKS parameter file ',fileName,
-     &' could not be found or is malformed. Please specify it.'
-      stop 
-C     Below is the code if one desires to let the code continue with
-C     a non existing or malformed parameter file
-      write(*,*) 'WARNING :: The file ',fileName,' could not be ',
-     & ' open or did not contain the necessary information. The ',
-     & ' default MadFKS parameters will be used.'
-      call DefaultFKSParam()
-      goto 998
-
-  998 continue
-
-      if(printParam.and..not.paramPrinted) then
-      write(*,*)
-     & '==============================================================='
-      if (couldRead) then      
-        write(*,*) 'INFO: MadFKS read these parameters from '
-     &,filename
-      else
-        write(*,*) 'INFO: MadFKS used the default parameters.'
-      endif
-      write(*,*)
-     & '==============================================================='
-      write(*,*) ' > IRPoleCheckThreshold      = ',IRPoleCheckThreshold
-      write(*,*) ' > PrecisionVirtualAtRunTime = '
-     $     ,PrecisionVirtualAtRunTime
-      if (SelectedContributionTypes(0).gt.0) then
-        write(*,*) ' > SelectedContributionTypes = ',
-     &   (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
-      else
-        write(*,*) ' > SelectedContributionTypes = All'
-      endif
-      if (VetoedContributionTypes(0).gt.0) then
-        write(*,*) ' > VetoedContributionTypes   = ',
-     &    (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
-      else
-        write(*,*) ' > VetoedContributionTypes   = None'
-      endif
-      if (QCD_squared_selected.eq.-1) then
-        write(*,*) ' > QCD_squared_selected      = All'
-      else
-        write(*,*) ' > QCD_squared_selected      = ',QCD_squared_selected
-      endif
-      if (QED_squared_selected.eq.-1) then
-        write(*,*) ' > QED_squared_selected      = All'
-      else
-        write(*,*) ' > QED_squared_selected      = ',QED_squared_selected
-      endif
-      if (SelectedCouplingOrders(1,0).gt.0) then
-        do j=1,SelectedCouplingOrders(1,0)
-          write(*,*) ' > SelectedCouplingOrders(',j,') = ',
-     &      (SelectedCouplingOrders(i,j),i=1,nsplitorders)
-        enddo
-      else
-        write(*,*) ' > SelectedCouplingOrders    = All'
-      endif
-      write(*,*) ' > NHelForMCoverHels         = ',NHelForMCoverHels
-      write(*,*) ' > VirtualFraction           = ',Virt_fraction
-      write(*,*) ' > MinVirtualFraction        = ',Min_virt_fraction
-      write(*,*) ' > SeparateFlavourConfigs    = ',separate_flavour_configs
-      write(*,*)
-     & '==============================================================='
-      paramPrinted=.TRUE.
-      endif
-
-      close(68)
-      HasReadOnce=.TRUE.
-  901 continue
-      end
-
-      subroutine DefaultFKSParam() 
-
-      implicit none
-      integer i,j
-      include "FKSParams.inc"
-
-      IRPoleCheckThreshold=1.0d-5
-      NHelForMCoverHels=5
-      PrecisionVirtualAtRunTime=1d-3
-      Virt_fraction=1d0
-      QED_squared_selected=-1
-      QCD_squared_selected=-1
-      Min_virt_fraction=0.005d0
-      separate_flavour_configs=.false.
-      IncludeBornContributions=.true.
-      SelectedContributionTypes(0)=0
-      VetoedContributionTypes(0)=0
-      do i=1, maxContribsSelected
-        SelectedContributionTypes(I)=-1
-        VetoedContributionTypes(I)=-1
-      enddo
-      do j=1,maxCouplingTypes
-        SelectedCouplingOrders(j,0) = 0
-      enddo
-      do j=1,maxCouplingsSelected
-        do i=1,maxCouplingTypes
-          SelectedCouplingOrders(i,j) = -1
-        enddo
-      enddo
-      end
 
=== added file 'Template/NLO/SubProcesses/FKSParams.f90'
--- Template/NLO/SubProcesses/FKSParams.f90	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/FKSParams.f90	2019-09-12 11:14:33 +0000
@@ -0,0 +1,248 @@
+!====================================================================
+!
+!  Define common block with all general parameters used by MadFKS 
+!  See their definitions in the file FKS_params.dat.
+!
+!====================================================================
+module FKSParams
+  character(len=64), parameter ::  paramFileName='FKS_params.dat'
+  integer,parameter :: maxContribsSelected=100, &
+                       maxCouplingsSelected=100, &
+                       maxContribType=15, &
+                       maxCouplingTypes=20
+  real*8 :: IRPoleCheckThreshold,Virt_fraction, PrecisionVirtualAtRunTime,Min_virt_fraction
+  integer  :: NHelForMCoverHels,VetoedContributionTypes(0:maxContribsSelected), &
+              SelectedContributionTypes(0:maxContribsSelected),QED_squared_selected, &
+              SelectedCouplingOrders(maxCouplingTypes,0:maxCouplingsSelected), &
+              QCD_squared_selected
+  logical :: separate_flavour_configs,IncludeBornContributions,use_poly_virtual
+
+contains
+
+  subroutine FKSParamReader(filename, printParam, force)
+    ! Reads the file 'filename' and sets the parameters found in that file.
+    implicit none
+    logical, save :: HasReadOnce=.False.,paramPrinted=.false.
+    logical :: force,couldRead,printParam
+    character(*) :: filename
+    CHARACTER(len=64) :: buff, buff2, mode
+    include "orders.inc"
+    integer :: i,j
+    couldRead=.False.
+    if (HasReadOnce.and..not.force) then
+       goto 901
+    endif
+! Make sure to have default parameters if not set in the FKSParams.dat card
+! (if it is an old one for instance)
+    call DefaultFKSParam()
+! Overwrite the default parameters from file:
+    open(68, file=fileName, err=676, action='READ')
+    do
+       read(68,*,end=999) buff
+       if(index(buff,'#').eq.1) then
+          if (buff .eq. '#IRPoleCheckThreshold') then
+             read(68,*,end=999) IRPoleCheckThreshold
+             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
+                stop 'IRPoleCheckThreshold must be >= -1.0d0.'
+             endif
+          elseif (buff .eq. '#PrecisionVirtualAtRunTime') then
+             read(68,*,end=999) PrecisionVirtualAtRunTime
+             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
+                stop 'PrecisionVirtualAtRunTime must be >= -1.0d0.'
+             endif
+          else if (buff .eq. '#NHelForMCoverHels') then
+             read(68,*,end=999) NHelForMCoverHels
+             if (NHelForMCoverHels .lt. -1) then
+                stop 'NHelForMCoverHels must be >= -1.'
+             endif
+          else if (buff .eq. '#QCD^2==') then
+             read(68,*,end=999) QCD_squared_selected
+             if (QCD_squared_selected .lt. -1) then
+                stop 'QCD_squared_selected must be >= -1.'
+             endif
+          else if (buff .eq. '#QED^2==') then
+             read(68,*,end=999) QED_squared_selected
+             if (QED_squared_selected .lt. -1) then
+                stop 'QED_squared_selected must be >= -1.'
+             endif
+          else if (buff .eq. '#VirtualFraction') then
+             read(68,*,end=999) Virt_fraction
+             if (Virt_fraction .lt. 0 .or. virt_fraction .gt.1) then
+                stop 'VirtualFraction should be a fraction between 0 and 1'
+             endif
+          else if (buff .eq. '#MinVirtualFraction') then
+             read(68,*,end=999) Min_Virt_fraction
+             if (min_virt_fraction .lt. 0 .or. min_virt_fraction .gt.1) then
+                stop 'VirtualFraction should be a fraction between 0 and 1'
+             endif
+          else if (buff .eq. '#SeparateFlavourConfigurations') then
+             read(68,*,end=999) separate_flavour_configs
+          else if (buff .eq. '#UsePolyVirtual') then
+             read(68,*,end=999) use_poly_virtual
+          else if (buff .eq. '#VetoedContributionTypes') then
+             read(68,*,end=999) VetoedContributionTypes(0)
+             if (VetoedContributionTypes(0) .lt. 0.or. &
+                  VetoedContributionTypes(0) .gt. maxContribsSelected) then
+                write(*,*) 'VetoedContributionTypes length should be >= 0 and <=', &
+                     maxContribsSelected
+                stop 'Format error in FKS_params.dat.'
+             endif
+             read(68,*,end=999) (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
+             do I=1,VetoedContributionTypes(0)
+                if (VetoedContributionTypes(I).lt.1.or. &
+                     VetoedContributionTypes(I).gt.maxContribType) then
+                   write(*,*) 'VetoedContributionTypes must be >=1 and <=',maxContribType
+                   stop 'Format error in FKS_params.dat.'
+                endif
+             enddo
+             do I=VetoedContributionTypes(0)+1,maxContribsSelected
+                VetoedContributionTypes(I)=-1
+             enddo
+
+          else if (buff .eq. '#SelectedContributionTypes') then
+             read(68,*,end=999) SelectedContributionTypes(0)
+             if (SelectedContributionTypes(0) .lt. 0 .or. &
+                  SelectedContributionTypes(0) .gt. maxContribsSelected) then
+                write(*,*) 'SelectedContributionTypes length should be >= 0 and <=', &
+                     maxContribsSelected
+                stop 'Format error in FKS_params.dat.'
+             endif
+             read(68,*,end=999) (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
+             do I=1,SelectedContributionTypes(0)
+                if (SelectedContributionTypes(I).lt.1.or. &
+                     SelectedContributionTypes(I).gt.maxContribType) then
+                   write(*,*) 'SelectedContributionTypes must be >=1 and <=',maxContribType
+                   stop 'Format error in FKS_params.dat.'
+                endif
+             enddo
+             do I=SelectedContributionTypes(0)+1,maxContribsSelected
+                SelectedContributionTypes(I)=-1
+             enddo
+          else if (buff .eq. '#SelectedCouplingOrders') then
+             read(68,*,end=999) SelectedCouplingOrders(1,0)
+             if (SelectedCouplingOrders(1,0) .lt. 0 .or. &
+                  SelectedCouplingOrders(1,0) .gt. maxCouplingsSelected) then
+                write(*,*) 'SelectedCouplingOrders length should be >= 0 and <=', &
+                     maxCouplingsSelected
+                stop 'Format error in FKS_params.dat.'                
+             endif
+             do j = 2, maxCouplingTypes
+                SelectedCouplingOrders(j,0) = SelectedCouplingOrders(1,0)
+             enddo
+             do j = 1, SelectedCouplingOrders(1,0)
+                read(68,*,end=999) (SelectedCouplingOrders(i,j),i=1,nsplitorders)
+                do i=nsplitorders+1,maxCouplingTypes
+                   SelectedCouplingOrders(i,j)=-1
+                enddo
+             enddo
+          else
+             write(*,*) 'The parameter name ',buff(2:),'is not reckognized.'
+             stop 'Format error in FKS_params.dat.'
+          endif
+       endif
+    enddo
+999 continue
+    couldRead=.True.
+    goto 998      
+
+676 continue
+    write(*,*) 'ERROR :: MadFKS parameter file ',fileName, &
+         ' could not be found or is malformed. Please specify it.'
+    stop 
+    !   Below is the code if one desires to let the code continue with
+    !   a non existing or malformed parameter file
+    write(*,*) 'WARNING :: The file ',fileName,' could not be ', &
+         ' open or did not contain the necessary information. The ', &
+         ' default MadFKS parameters will be used.'
+    call DefaultFKSParam()
+    goto 998
+998 continue
+
+    if(printParam.and..not.paramPrinted) then
+       write(*,*) &
+            '==============================================================='
+       if (couldRead) then      
+          write(*,*) 'INFO: MadFKS read these parameters from ',filename
+       else
+          write(*,*) 'INFO: MadFKS used the default parameters.'
+       endif
+       write(*,*) &
+            '==============================================================='
+       write(*,*) ' > IRPoleCheckThreshold      = ',IRPoleCheckThreshold
+       write(*,*) ' > PrecisionVirtualAtRunTime = ',PrecisionVirtualAtRunTime
+       if (SelectedContributionTypes(0).gt.0) then
+          write(*,*) ' > SelectedContributionTypes = ', &
+               (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
+       else
+          write(*,*) ' > SelectedContributionTypes = All'
+       endif
+       if (VetoedContributionTypes(0).gt.0) then
+          write(*,*) ' > VetoedContributionTypes   = ', &
+               (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
+       else
+          write(*,*) ' > VetoedContributionTypes   = None'
+       endif
+       if (QCD_squared_selected.eq.-1) then
+          write(*,*) ' > QCD_squared_selected      = All'
+       else
+          write(*,*) ' > QCD_squared_selected      = ',QCD_squared_selected
+       endif
+       if (QED_squared_selected.eq.-1) then
+          write(*,*) ' > QED_squared_selected      = All'
+       else
+          write(*,*) ' > QED_squared_selected      = ',QED_squared_selected
+       endif
+       if (SelectedCouplingOrders(1,0).gt.0) then
+          do j=1,SelectedCouplingOrders(1,0)
+             write(*,*) ' > SelectedCouplingOrders(',j,') = ', &
+                  (SelectedCouplingOrders(i,j),i=1,nsplitorders)
+          enddo
+       else
+          write(*,*) ' > SelectedCouplingOrders    = All'
+       endif
+       write(*,*) ' > NHelForMCoverHels         = ',NHelForMCoverHels
+       write(*,*) ' > VirtualFraction           = ',Virt_fraction
+       write(*,*) ' > MinVirtualFraction        = ',Min_virt_fraction
+       write(*,*) ' > SeparateFlavourConfigs    = ',separate_flavour_configs
+       write(*,*) ' > UsePolyVirtual            = ',use_poly_virtual
+       write(*,*) &
+            '==============================================================='
+       paramPrinted=.TRUE.
+    endif
+
+    close(68)
+    HasReadOnce=.TRUE.
+901 continue
+  end subroutine FKSParamReader
+
+  subroutine DefaultFKSParam() 
+    ! Sets the default parameters
+    implicit none
+    integer i,j
+    IRPoleCheckThreshold=1.0d-5
+    NHelForMCoverHels=5
+    PrecisionVirtualAtRunTime=1d-3
+    Virt_fraction=1d0
+    QED_squared_selected=-1
+    QCD_squared_selected=-1
+    Min_virt_fraction=0.005d0
+    separate_flavour_configs=.false.
+    use_poly_virtual=.true.
+    IncludeBornContributions=.true.
+    SelectedContributionTypes(0)=0
+    VetoedContributionTypes(0)=0
+    do i=1, maxContribsSelected
+       SelectedContributionTypes(I)=-1
+       VetoedContributionTypes(I)=-1
+    enddo
+    do j=1,maxCouplingTypes
+       SelectedCouplingOrders(j,0) = 0
+    enddo
+    do j=1,maxCouplingsSelected
+       do i=1,maxCouplingTypes
+          SelectedCouplingOrders(i,j) = -1
+       enddo
+    enddo
+  end subroutine DefaultFKSParam
+
+end module FKSParams
 
=== removed file 'Template/NLO/SubProcesses/FKSParams.inc'
--- Template/NLO/SubProcesses/FKSParams.inc	2017-08-04 19:52:15 +0000
+++ Template/NLO/SubProcesses/FKSParams.inc	1970-01-01 00:00:00 +0000
@@ -1,31 +0,0 @@
-!====================================================================
-!
-!  Define common block with all general parameters used by MadFKS 
-!  See their definitions in the file FKS_params.dat.
-!
-!====================================================================
-
-      character*64 paramFileName
-      parameter ( paramFileName='FKS_params.dat')
-
-	  integer maxContribsSelected, maxCouplingsSelected
-	  parameter( maxContribsSelected=100, maxCouplingsSelected=100 )
-	  integer maxContribType, maxCouplingTypes
-	  parameter ( maxContribType=15, maxCouplingTypes=20 )
-
-      real*8 IRPoleCheckThreshold,Virt_fraction,
-     &	     PrecisionVirtualAtRunTime,Min_virt_fraction
-
-      integer NHelForMCoverHels, VetoedContributionTypes(0:maxContribsSelected),
-     &        SelectedContributionTypes(0:maxContribsSelected),
-     &        QED_squared_selected, QCD_squared_selected,
-     &        SelectedCouplingOrders(maxCouplingTypes,0:maxCouplingsSelected)
-
-      logical separate_flavour_configs, IncludeBornContributions
-
-      common /FKSPARAMS/IRPoleCheckThreshold,virt_fraction,
-     &       PrecisionVirtualAtRunTime,min_virt_fraction,
-     &       NHelForMCoverHels,separate_flavour_configs,
-     &       QCD_squared_selected, QED_squared_selected,
-     &       VetoedContributionTypes, SelectedContributionTypes,
-     &       SelectedCouplingOrders
 
=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2018-12-21 15:27:32 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2019-09-12 11:14:33 +0000
@@ -3,6 +3,7 @@
 c     This is the driver for the whole calulation
 c**************************************************************************
       use mint_module
+      use FKSParams
       implicit none
 C
 C     CONSTANTS
@@ -64,7 +65,6 @@
       common /to_polecheck/force_polecheck, polecheck_passed
       integer ret_code_ml
       common /to_ret_code/ret_code_ml
-      include 'FKSParams.inc'
       
 C-----
 C  BEGIN CODE
 
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2019-09-12 09:03:04 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2019-09-12 11:14:33 +0000
@@ -4,6 +4,7 @@
 c**************************************************************************
       use extra_weights
       use mint_module
+      use FKSParams
       implicit none
 C
 C     CONSTANTS
@@ -66,9 +67,6 @@
       include "timing_variables.inc"
       real*4 tOther, tTot
 
-c general MadFKS parameters
-      include "FKSParams.inc"
-
 c applgrid
       integer iappl
       common /for_applgrid/ iappl
 
=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2019-09-12 09:03:04 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2019-09-12 11:14:33 +0000
@@ -4,6 +4,7 @@
 c**************************************************************************
       use extra_weights
       use mint_module
+      use FKSParams
       implicit none
 C
 C     CONSTANTS
@@ -74,9 +75,6 @@
       include "timing_variables.inc"
       real*4 tOther, tTot
 
-c general MadFKS parameters
-      include "FKSParams.inc"
-
       double precision deravg,derstd,dermax,xi_i_fks_ev_der_max
      &     ,y_ij_fks_ev_der_max
       integer ntot_granny,derntot,ncase(0:6)
 
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2019-09-12 09:03:04 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2019-09-12 11:14:33 +0000
@@ -1456,6 +1456,7 @@
 c        contribution
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
@@ -1463,7 +1464,6 @@
       include 'coupl.inc'
       include 'fks_info.inc'
       include 'q_es.inc'
-      include 'FKSParams.inc'
       include 'orders.inc'
       integer type,i,j
       logical foundIt,foundOrders
@@ -1651,6 +1651,7 @@
       use weight_lines
       use extra_weights
       use mint_module
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
@@ -1658,7 +1659,6 @@
       include 'timing_variables.inc'
       include 'genps.inc'
       include 'orders.inc'
-      include 'FKSParams.inc'
       integer orders(nsplitorders)
       integer i,j,k,iamp,icontr_orig
       logical virt_found
@@ -2011,11 +2011,11 @@
 c wgts() array to include the weights.
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'timing_variables.inc'
-      include 'FKSParams.inc'
       include 'genps.inc'
       integer i,kr,kf,iwgt_save,dd
       double precision xlum(maxscales),dlum,pi,mu2_r(maxscales),c_mu2_r
@@ -2094,11 +2094,11 @@
 c computations (ickkw.eq.-1).
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'timing_variables.inc'
-      include 'FKSParams.inc'
       include 'genps.inc'
       integer i,ks,kh,iwgt_save
       double precision xlum(maxscales),dlum,pi,mu2_r(maxscales)
@@ -2200,11 +2200,11 @@
 c wgts() array to include the weights.
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'timing_variables.inc'
-      include 'FKSParams.inc'
       include 'genps.inc'
       integer n,izero,i,nn
       parameter (izero=0)
@@ -5877,12 +5877,22 @@
             tOLP=tOLP+(tAfter-tBefore)
             virtual_over_born=virt_wgt/born_wgt
             if (ickkw.ne.-1) then
-               virt_wgt=virt_wgt-average_virtual(0,ichan)*born_wgt
+               if (use_poly_virtual) then
+                  virt_wgt=virt_wgt-polyfit(0)*born_wgt
+               else
+                  virt_wgt=virt_wgt-average_virtual(0,ichan)*born_wgt
+               endif
                do iamp=1,amp_split_size
                   if (amp_split_virt(iamp).eq.0d0) cycle
+                  if (use_poly_virtual) then
+                     amp_split_virt(iamp)=amp_split_virt(iamp)-
+     $                    polyfit(iamp)
+     $                    *amp_split_born_for_virt(iamp)
+               else
                   amp_split_virt(iamp)=amp_split_virt(iamp)-
      $                 average_virtual(iamp,ichan)
      $                 *amp_split_born_for_virt(iamp)
+                  endif
                enddo
             endif
             if (abrv.ne.'virt') then
@@ -5902,12 +5912,21 @@
      $        amp_split_virt_save(1:amp_split_size)
       endif
       if (abrv(1:4).ne.'virt' .and. ickkw.ne.-1) then
-         avv_wgt=average_virtual(0,ichan)*born_wgt
-         do iamp=1, amp_split_size
-            if (amp_split_born_for_virt(iamp).eq.0d0) cycle
-            amp_split_avv(iamp)= average_virtual(iamp,ichan)
-     $           *amp_split_born_for_virt(iamp)
-         enddo
+         if (use_poly_virtual) then
+            avv_wgt=polyfit(0)*born_wgt
+            do iamp=1, amp_split_size
+               if (amp_split_born_for_virt(iamp).eq.0d0) cycle
+               amp_split_avv(iamp)= polyfit(iamp)
+     $              *amp_split_born_for_virt(iamp)
+            enddo
+         else
+            avv_wgt=average_virtual(0,ichan)*born_wgt
+            do iamp=1, amp_split_size
+               if (amp_split_born_for_virt(iamp).eq.0d0) cycle
+               amp_split_avv(iamp)= average_virtual(iamp,ichan)
+     $              *amp_split_born_for_virt(iamp)
+            enddo
+         endif
       endif
 
 c eq.(MadFKS.C.13)
 
=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
--- Template/NLO/SubProcesses/makefile_fks_dir	2018-12-21 15:27:32 +0000
+++ Template/NLO/SubProcesses/makefile_fks_dir	2019-09-12 11:14:33 +0000
@@ -26,16 +26,16 @@
 
 # Files for all executables
 FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) $(patsubst	\
-	%.f,%.o,$(wildcard matrix_*.f)) real_me_chooser.o		\
+	%.f,%.o,$(wildcard matrix_*.f)) FKSParams.o real_me_chooser.o	\
 	chooser_functions.o genps_fks.o setcuts.o setscales.o		\
 	veto_xsec.o $(patsubst %.f,%.o,$(wildcard b_sf_???.f))		\
 	born.o sborn_sf.o extra_cnt_wrapper.o $(patsubst		\
 	%.f,%.o,$(wildcard born_cnt_*.f)) fks_Sij.o			\
 	$(fastjetfortran_madfks) fks_singular.o montecarlocounter.o	\
 	reweight_xsec.o boostwdir2.o initcluster.o cluster.o		\
-	splitorders_stuff.o reweight.o get_color.o FKSParamReader.o	\
+	splitorders_stuff.o reweight.o get_color.o 			\
 	iproc_map.o MC_integer.o $(reweight_xsec_events_pdf_dummy)	\
-	$(applgrid_interface) weight_lines.o mint_module.o
+	$(applgrid_interface) weight_lines.o mint_module.o polfit.o
 
 # Files needed for mintFO & mintMC
 RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o 		\
 
=== modified file 'Template/NLO/SubProcesses/mint_module.f90'
--- Template/NLO/SubProcesses/mint_module.f90	2019-08-15 13:57:46 +0000
+++ Template/NLO/SubProcesses/mint_module.f90	2019-09-12 11:14:33 +0000
@@ -64,8 +64,8 @@
 ! nintegrals>6 : virtual and born order by order
 !
 
-
 module mint_module
+  use FKSParams ! contains use_poly_virtual
   implicit none
   integer, parameter, private :: nintervals=32    ! max number of intervals in the integration grids
   integer, parameter, public  :: ndimmax=60       ! max number of dimensions of the integral
@@ -92,7 +92,7 @@
   integer, dimension(maxchannels), public :: iconfigs
   double precision, public :: accuracy,min_virt_fraction_mint,wgt_mult
   double precision, dimension(0:n_ave_virt,maxchannels), public :: average_virtual
-  double precision, dimension(0:n_ave_virt), public :: virt_wgt_mint,born_wgt_mint
+  double precision, dimension(0:n_ave_virt), public :: virt_wgt_mint,born_wgt_mint,polyfit
   double precision, dimension(maxchannels), public :: virtual_fraction
   double precision, dimension(nintegrals,0:maxchannels), public :: ans,unc
   logical :: only_virt,new_point,pass_cuts_check
@@ -340,9 +340,13 @@
        ! overwrite xgrid with the new xgrid
        if (regridded(kchan)) xgrid(1:nint_used,1:ndim,kchan)=xgrid_new(1:nint_used,1:ndim)
     enddo
-    do k_ord_virt=0,n_ord_virt
-       call regrid_ave_virt(k_ord_virt)
-    enddo
+    if (use_poly_virtual) then
+       call do_polyfit()
+    else
+       do k_ord_virt=0,n_ord_virt
+          call regrid_ave_virt(k_ord_virt)
+       enddo
+    endif
 ! Regrid the MC over integers (used for the MC over FKS dirs)
     call regrid_MC_integer
   end subroutine update_integration_grids
@@ -664,10 +668,18 @@
           isix=2*k_ord_virt+6
        endif
        if (f(ithree).ne.0d0) then
-          f(isix)=f(isix)/virtual_fraction(ichan)
-          virtual=(f(ithree)+average_virtual(k_ord_virt,ichan)*f(isix))*virtual_fraction(ichan)
-          born=f(isix)*virtual_fraction(ichan)
-          call fill_ave_virt(x,k_ord_virt,virtual,born)
+          born=f(isix)
+          ! virt_wgt_mint=(virtual-average_virtual*born)/virtual_fraction. Compensate:
+          if (use_poly_virtual) then
+             virtual=f(ithree)*virtual_fraction(ichan)+ &
+                  polyfit(k_ord_virt)*f(isix)
+             call add_point_polyfit(ichan,k_ord_virt,x(1:ndim-3), &
+                  virtual/born,born/wgt_mult)
+          else
+             virtual=f(ithree)*virtual_fraction(ichan)+ &
+                  average_virtual(k_ord_virt,ichan)*f(isix)
+             call fill_ave_virt(x,k_ord_virt,virtual,born)
+          endif
        else
           f(isix)=0d0
        endif
@@ -801,7 +813,11 @@
        if(imode.eq.0) nhits(icell(kdim),kdim,ichan)=nhits(icell(kdim),kdim,ichan)+1
     enddo
     do k_ord_virt=0,n_ord_virt
-       call get_ave_virt(x,k_ord_virt)
+       if (use_poly_virtual) then
+          call get_polyfit(ichan,k_ord_virt,x(1:ndim-3),polyfit(k_ord_virt))
+       else
+          call get_ave_virt(x,k_ord_virt)
+       endif
     enddo
   end subroutine get_random_x
   
@@ -988,14 +1004,18 @@
 
   subroutine reset_mint_grids
     implicit none
-    integer :: kdim,kchan,kint,k_ord_virt
+    integer :: kdim,kchan,kint
     do kint=0,nint_used
        xgrid(kint,1:ndim,1:nchans)=dble(kint)/nint_used
     enddo
     nhits(1:nint_used,1:ndim,1:nchans)=0
     regridded(1:nchans)=.true.
     nhits_in_grids(1:nchans)=0
-    call init_ave_virt
+    if (use_poly_virtual) then
+       call init_polyfit(ndim-3,nchans,n_ord_virt,1000)
+    else
+       call init_ave_virt
+    endif
     ans_chan(0:nchans)=0d0
     if (double_events) then
        ! when double events, start with the very first channel only. For the
@@ -1032,11 +1052,13 @@
              write (12,*) 'MAX',(ymax(j,i,kchan),i=1,ndim)
           enddo
        endif
-       do j=1,nintervals_virt
-          do k=0,n_ord_virt
-             write (12,*) 'AVE',(ave_virt(j,i,k,kchan),i=1,ndim)
+       if (.not.use_poly_virtual) then
+          do j=1,nintervals_virt
+             do k=0,n_ord_virt
+                write (12,*) 'AVE',(ave_virt(j,i,k,kchan),i=1,ndim)
+             enddo
           enddo
-       enddo
+       endif
        if (imode.ge.1) then
           write (12,*) 'MAX',ymax_virt(kchan)
        endif
@@ -1046,6 +1068,7 @@
        write (12,*) 'AVE',virtual_fraction(kchan),average_virtual(0,kchan)
     enddo
     write (12,*) 'IDE',(ifold(i),i=1,ndim)
+    if (use_poly_virtual) call save_polyfit(12)
     close (12)
   end subroutine write_grids_to_file
   
@@ -1053,6 +1076,7 @@
 ! Read the MINT integration grids from file
     implicit none
     integer :: i,j,k,kchan,idum
+    integer,dimension(maxchannels) :: points
     character(len=3) :: dummy
     open (unit=12, file='mint_grids',status='old')
     ans(1,0)=0d0
@@ -1066,11 +1090,13 @@
              read (12,*) dummy,(ymax(j,i,kchan),i=1,ndim)
           enddo
        endif
-       do j=1,nintervals_virt
-          do k=0,n_ord_virt
-             read (12,*) dummy,(ave_virt(j,i,k,kchan),i=1,ndim)
+       if (.not.use_poly_virtual) then
+          do j=1,nintervals_virt
+             do k=0,n_ord_virt
+                read (12,*) dummy,(ave_virt(j,i,k,kchan),i=1,ndim)
+             enddo
           enddo
-       enddo
+       endif
        if (imode.ge.2) then
           read (12,*) dummy,ymax_virt(kchan)
        endif
@@ -1083,6 +1109,18 @@
     enddo
     read (12,*) dummy,(ifold(i),i=1,ndim)
     unc(1,0)=sqrt(unc(1,0))
+    ! polyfit stuff:
+    if (use_poly_virtual) then
+       do kchan=1,nchans
+          read (12,*) dummy,points(kchan)
+       enddo
+       do kchan=1,nchans
+          backspace(12)
+       enddo
+       call init_polyfit(ndim-3,nchans,n_ord_virt,maxval(points(1:nchans)))
+       call restore_polyfit(12)
+       call do_polyfit()
+    endif
     close (12)
 ! check for zero cross-section: if restoring grids corresponding to
 ! sigma=0, just terminate the run
 
=== added file 'Template/NLO/SubProcesses/polfit.f'
--- Template/NLO/SubProcesses/polfit.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/polfit.f	2019-09-12 11:14:33 +0000
@@ -0,0 +1,720 @@
+      subroutine init_polyfit(ndims,nchans,k_ord_virt,npoints)
+! wrapper subroutine that saves all information that's needed for calls
+! to the polfit() and pvalue() subroutines. Since I was too lazy to
+! write a proper module, use fortran entry statements to jump to the
+! right locations in this subroutine.
+      implicit none
+      integer ndims,nchans,ndim,nchan,maxpoint,maxdeg,ierr,i,ichan,ic
+     $     ,iunit,npoints,n_ord_virt,k_ord_virt,ko
+      integer,allocatable,dimension(:) :: n
+      integer,allocatable,dimension(:,:,:) :: ndeg
+      real*8,allocatable,dimension(:) :: r,yp,a
+      real*8,allocatable,dimension(:,:,:) :: y,w,x2d,temp3
+      real*8,allocatable,dimension(:,:,:,:) :: a2d
+      logical,allocatable,dimension(:) :: valid_ord_virt
+      parameter (maxdeg=20)
+      real*8 absXS,virt,ave_fun,fun_at_x,eps,sum_w,yfit,x(*)
+      logical fit_done,verbose
+      parameter (verbose=.false.)
+      character(len=3) :: dummy
+      save
+      maxpoint=npoints
+      nchan=nchans
+      ndim=ndims
+      n_ord_virt=k_ord_virt
+      if (allocated(ndeg)) deallocate(ndeg)
+      allocate(ndeg(ndim,0:n_ord_virt,nchan))
+      if (allocated(n)) deallocate(n)
+      allocate(n(nchan))
+      if (allocated(a)) deallocate(a)
+      allocate(a(maxpoint*3+maxdeg*3+3))
+      if (allocated(y)) deallocate(y)
+      allocate(y(maxpoint,0:n_ord_virt,nchan))
+      if (allocated(w)) deallocate(w)
+      allocate(w(maxpoint,0:n_ord_virt,nchan))
+      if (allocated(r)) deallocate(r)
+      allocate(r(maxpoint))
+      if (allocated(yp)) deallocate(yp)
+      allocate(yp(maxdeg))
+      if (allocated(x2d)) deallocate(x2d)
+      allocate(x2d(maxpoint,ndim,nchan))
+      if (allocated(a2d)) deallocate(a2d)
+      allocate(a2d(maxdeg*3+3,ndim,0:n_ord_virt,nchan))
+      if (allocated(valid_ord_virt)) deallocate(valid_ord_virt)
+      allocate(valid_ord_virt(0:n_ord_virt))
+      n(1:nchan)=0
+      fit_done=.false.
+      valid_ord_virt(0:n_ord_virt)=.false.
+      return
+
+      entry add_point_polyfit(ichan,k_ord_virt,x,virt,absXS)
+      ! only increment 'n' and 'x2d' if it's the zeroth "order" in the
+      ! virtual/born ratio
+      valid_ord_virt(k_ord_virt)=.true.
+      if (k_ord_virt.eq.0) then
+         if (n(ichan).eq.maxpoint) then
+! increase the size of the allocated variables:
+            maxpoint=maxpoint+1000
+            allocate(temp3(maxpoint,0:n_ord_virt,nchan))
+            do i=1,nchan
+               temp3(1:n(i),0:n_ord_virt,i)=y(1:n(i),0:n_ord_virt,i)
+            enddo
+            call move_alloc(temp3,y)
+            allocate(temp3(maxpoint,0:n_ord_virt,nchan))
+            do i=1,nchan
+               temp3(1:n(i),0:n_ord_virt,i)=w(1:n(i),0:n_ord_virt,i)
+            enddo
+            call move_alloc(temp3,w)
+            deallocate(r)
+            allocate(r(maxpoint))
+            allocate(temp3(maxpoint,ndim,nchan))
+            do i=1,nchan
+               temp3(1:n(i),1:ndim,i)=x2d(1:n(i),1:ndim,i)
+            enddo
+            call move_alloc(temp3,x2d)
+            deallocate(a)
+            allocate(a(maxpoint*3+maxdeg*3+3))
+         endif
+         n(ichan)=n(ichan)+1
+! add the point to the saved information
+         do i=1,ndim
+            x2d(n(ichan),i,ichan)=x(i)
+         enddo
+! initialise all to zero
+         y(n(ichan),0:n_ord_virt,ichan)=0d0
+         w(n(ichan),0:n_ord_virt,ichan)=0d0
+      endif
+      y(n(ichan),k_ord_virt,ichan)=virt
+      w(n(ichan),k_ord_virt,ichan)=sqrt(abs(absXS))
+      return
+
+      entry do_polyfit()
+      ! perform the fit for each dimension and each channel
+      do ic=1,nchan
+         do ko=0,n_ord_virt
+            if (.not.valid_ord_virt(ko)) cycle
+            do i=1,ndim
+               eps=-1d0
+               call polfit(n(ic),x2d(1,i,ic),y(1,ko,ic),w(1,ko,ic)
+     $                          ,maxdeg,ndeg(i,ko,ic),eps,r,ierr,a)
+               ! this is the information used by pvalue() to interpolate
+               ! later:
+               a2d(1:maxdeg*3+3,i,ko,ic)=a(1:maxdeg*3+3)
+               ! print something in the logs:
+               if (verbose) write (*,'(a,i3,i3,i3,i8,i3,i3,f10.4)')
+     $              'polyfit -',ic,ko,i,n(ic),ndeg(i,ko,ic),ierr,eps
+            enddo
+         enddo
+      enddo
+      fit_done=.true.
+      return
+
+      entry get_polyfit(ichan,k_ord_virt,x,fun_at_x)
+      ! if we haven't do a fit to the data, just return zero. This is
+      ! typically the case only during the very first iteration.
+      if ((.not.fit_done) .or. (.not.valid_ord_virt(k_ord_virt))) then
+         fun_at_x=0d0
+         return
+      endif
+      ! Do the interpolation. Since we have several dimensions, we use
+      ! the following:
+      !       f(x) = average + sum_i (f_i(x_i)-average)
+      ! where 'i' loops over all dimensions.  In other words, the
+      ! information from each of the separate dimensions should only add
+      ! (or subtract) something from the average value. To obtain the
+      ! average, simply call pvalue() to return the zeroth-order
+      ! polynomial.
+      fun_at_x=0d0
+      do i=1,ndim
+         a(1:maxdeg*3+3)=a2d(1:maxdeg*3+3,i,k_ord_virt,ichan)
+         call pvalue(ndeg(i,k_ord_virt,ichan),0,x(i),yfit,yp,a)
+         call pvalue(0,0,x(i),ave_fun,yp,a)
+         fun_at_x=fun_at_x+(yfit-ave_fun)+ave_fun/dble(ndim)
+      enddo
+      return
+
+      entry save_polyfit(iunit)
+      do ic=1,nchan
+         write (iunit,*) 'POL',n(ic)
+      enddo
+      do ic=1,nchan
+         do i=1,n(ic)
+            write (iunit,*) 'POL',x2d(i,1:ndim,ic),
+     &           y(i,0:n_ord_virt,ic),w(i,0:n_ord_virt,ic)
+         enddo
+      enddo
+      return
+
+      entry restore_polyfit(iunit)
+      do ic=1,nchan
+         read (iunit,*) dummy,n(ic)
+      enddo
+      do ic=1,nchan
+         do i=1,n(ic)
+            read (iunit,*) dummy,x2d(i,1:ndim,ic)
+     &           ,y(i,0:n_ord_virt,ic),w(i,0:n_ord_virt,ic)
+         enddo
+      enddo
+      do ic=1,nchan
+         do ko=0,n_ord_virt
+            do i=1,n(ic)
+               if (y(i,ko,ic).ne.0d0) valid_ord_virt(ko)=.true.
+            enddo
+         enddo
+      enddo
+      if (verbose) write (*,*) 'polyfit - valid orders virt'
+     $     ,valid_ord_virt(0:n_ord_virt)
+      return
+      end
+
+
+
+
+
+*DECK POLFIT      
+      SUBROUTINE POLFIT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A)
+c
+c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+cRF: make double precision and replace DO-loops and arirthmic
+c    IF-statements to comply more modern standards. Since we made it
+c    double precision, the 'extended, partial-double' precision part has
+c    become useless. Could remove it, but there is no real reason for
+c    that (should only be marginally faster).
+c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+c
+C***BEGIN PROLOGUE  POLFIT
+C***PURPOSE  Fit discrete data in a least squares sense by polynomials
+C            in one variable.
+C***LIBRARY   SLATEC
+C***CATEGORY  K1A1A2
+C***TYPE      SINGLE PRECISION (POLFIT-S, DPOLFT-D)
+C***KEYWORDS  CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
+C***AUTHOR  Shampine, L. F., (SNLA)
+C           Davenport, S. M., (SNLA)
+C           Huddleston, R. E., (SNLL)
+C***DESCRIPTION
+C
+C     Abstract
+C
+C     Given a collection of points X(I) and a set of values Y(I) which
+C     correspond to some function or measurement at each of the X(I),
+C     subroutine  POLFIT  computes the weighted least-squares polynomial
+C     fits of all degrees up to some degree either specified by the user
+C     or determined by the routine.  The fits thus obtained are in
+C     orthogonal polynomial form.  Subroutine  PVALUE  may then be
+C     called to evaluate the fitted polynomials and any of their
+C     derivatives at any point.  The subroutine  PCOEF  may be used to
+C     express the polynomial fits as powers of (X-C) for any specified
+C     point C.
+C
+C     The parameters for  POLFIT  are
+C
+C     Input --
+C         N -      the number of data points.  The arrays X, Y and W
+C                  must be dimensioned at least  N  (N .GE. 1).
+C         X -      array of values of the independent variable.  These
+C                  values may appear in any order and need not all be
+C                  distinct.
+C         Y -      array of corresponding function values.
+C         W -      array of positive values to be used as weights.  If
+C                  W(1) is negative,  POLFIT  will set all the weights
+C                  to 1.0, which means unweighted least squares error
+C                  will be minimized.  To minimize relative error, the
+C                  user should set the weights to:  W(I) = 1.0/Y(I)**2,
+C                  I = 1,...,N .
+C         MAXDEG - maximum degree to be allowed for polynomial fit.
+C                  MAXDEG  may be any non-negative integer less than  N.
+C                  Note -- MAXDEG  cannot be equal to  N-1  when a
+C                  statistical test is to be used for degree selection,
+C                  i.e., when input value of  EPS  is negative.
+C         EPS -    specifies the criterion to be used in determining
+C                  the degree of fit to be computed.
+C                  (1)  If  EPS  is input negative,  POLFIT  chooses the
+C                       degree based on a statistical F test of
+C                       significance.  One of three possible
+C                       significance levels will be used:  .01, .05 or
+C                       .10.  If  EPS=-1.0 , the routine will
+C                       automatically select one of these levels based
+C                       on the number of data points and the maximum
+C                       degree to be considered.  If  EPS  is input as
+C                       -.01, -.05, or -.10, a significance level of
+C                       .01, .05, or .10, respectively, will be used.
+C                  (2)  If  EPS  is set to 0.,  POLFIT  computes the
+C                       polynomials of degrees 0 through  MAXDEG .
+C                  (3)  If  EPS  is input positive,  EPS  is the RMS
+C                       error tolerance which must be satisfied by the
+C                       fitted polynomial.  POLFIT  will increase the
+C                       degree of fit until this criterion is met or
+C                       until the maximum degree is reached.
+C
+C     Output --
+C         NDEG -   degree of the highest degree fit computed.
+C         EPS -    RMS error of the polynomial of degree  NDEG .
+C         R -      vector of dimension at least NDEG containing values
+C                  of the fit of degree  NDEG  at each of the  X(I) .
+C                  Except when the statistical test is used, these
+C                  values are more accurate than results from subroutine
+C                  PVALUE  normally are.
+C         IERR -   error flag with the following possible values.
+C             1 -- indicates normal execution, i.e., either
+C                  (1)  the input value of  EPS  was negative, and the
+C                       computed polynomial fit of degree  NDEG
+C                       satisfies the specified F test, or
+C                  (2)  the input value of  EPS  was 0., and the fits of
+C                       all degrees up to  MAXDEG  are complete, or
+C                  (3)  the input value of  EPS  was positive, and the
+C                       polynomial of degree  NDEG  satisfies the RMS
+C                       error requirement.
+C             2 -- invalid input parameter.  At least one of the input
+C                  parameters has an illegal value and must be corrected
+C                  before  POLFIT  can proceed.  Valid input results
+C                  when the following restrictions are observed
+C                       N .GE. 1
+C                       0 .LE. MAXDEG .LE. N-1  for  EPS .GE. 0.
+C                       0 .LE. MAXDEG .LE. N-2  for  EPS .LT. 0.
+C                       W(1)=-1.0  or  W(I) .GT. 0., I=1,...,N .
+C             3 -- cannot satisfy the RMS error requirement with a
+C                  polynomial of degree no greater than  MAXDEG .  Best
+C                  fit found is of degree  MAXDEG .
+C             4 -- cannot satisfy the test for significance using
+C                  current value of  MAXDEG .  Statistically, the
+C                  best fit found is of order  NORD .  (In this case,
+C                  NDEG will have one of the values:  MAXDEG-2,
+C                  MAXDEG-1, or MAXDEG).  Using a higher value of
+C                  MAXDEG  may result in passing the test.
+C         A -      work and output array having at least 3N+3MAXDEG+3
+C                  locations
+C
+C     Note - POLFIT  calculates all fits of degrees up to and including
+C            NDEG .  Any or all of these fits can be evaluated or
+C            expressed as powers of (X-C) using  PVALUE  and  PCOEF
+C            after just one call to  POLFIT .
+C
+C***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
+C                 Curve fitting by polynomials in one variable, Report
+C                 SLA-74-0270, Sandia Laboratories, June 1974.
+C***ROUTINES CALLED  PVALUE, XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   740601  DATE WRITTEN
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890531  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C   920527  Corrected erroneous statements in DESCRIPTION.  (WRB)
+C***  END PROLOGUE  POLFIT
+      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
+      DOUBLE PRECISION TEMD1,TEMD2
+      DIMENSION X(*), Y(*), W(*), R(*), A(*), YP(0)
+      DIMENSION CO(4,3)
+      SAVE CO
+      DATA  CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2),
+     1      CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3),
+     2  CO(4,3)/-13.086850,-2.4648165,-3.3846535,-1.2973162,
+     3          -3.3381146,-1.7812271,-3.2578406,-1.6589279,
+     4          -1.6282703,-1.3152745,-3.2640179,-1.9829776/
+C***  FIRST EXECUTABLE STATEMENT  POLFIT
+      M = ABS(N)
+      IF (M .EQ. 0) GO TO 30
+      IF (MAXDEG .LT. 0) GO TO 30
+      A(1) = MAXDEG
+      MOP1 = MAXDEG + 1
+      IF (M .LT. MOP1) GO TO 30
+      IF (EPS .LT. 0.0  .AND.  M .EQ. MOP1) GO TO 30
+      XM = M
+      ETST = EPS*EPS*XM
+      IF (W(1) .LT. 0.0) GO TO 2
+      DO I = 1,M
+        IF (W(I) .LE. 0.0) GO TO 30
+      enddo
+      GO TO 4
+ 2    DO I = 1,M
+         W(I) = 1.0
+      enddo
+ 4    IF (EPS .GE. 0.0) GO TO 8
+C
+C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR
+C CHOOSING DEGREE OF POLYNOMIAL FIT
+C
+      IF (EPS .GT. (-.55)) GO TO 5
+      IDEGF = M - MAXDEG - 1
+      KSIG = 1
+      IF (IDEGF .LT. 10) KSIG = 2
+      IF (IDEGF .LT. 5) KSIG = 3
+      GO TO 8
+ 5    KSIG = 1
+      IF (EPS .LT. (-.03)) KSIG = 2
+      IF (EPS .LT. (-.07)) KSIG = 3
+C
+C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING
+C
+ 8    K1 = MAXDEG + 1
+      K2 = K1 + MAXDEG
+      K3 = K2 + MAXDEG + 2
+      K4 = K3 + M
+      K5 = K4 + M
+      DO I = 2,K4
+         A(I) = 0.0
+      enddo
+      W11 = 0.0
+      IF (N .LT. 0) GO TO 11
+C
+C UNCONSTRAINED CASE
+C
+      DO I = 1,M
+        K4PI = K4 + I
+        A(K4PI) = 1.0
+        W11 = W11 + W(I)
+      enddo
+      GO TO 13
+C
+C CONSTRAINED CASE
+C
+ 11   DO I = 1,M
+        K4PI = K4 + I
+        W11 = W11 + W(I)*A(K4PI)**2
+      enddo
+C
+C COMPUTE FIT OF DEGREE ZERO
+C
+ 13   TEMD1 = 0.0D0
+      DO I = 1,M
+        K4PI = K4 + I
+        TEMD1 = TEMD1 + DBLE(W(I))*DBLE(Y(I))*DBLE(A(K4PI))
+      enddo
+      TEMD1 = TEMD1/DBLE(W11)
+      A(K2+1) = TEMD1
+      SIGJ = 0.0
+      DO I = 1,M
+        K4PI = K4 + I
+        K5PI = K5 + I
+        TEMD2 = TEMD1*DBLE(A(K4PI))
+        R(I) = TEMD2
+        A(K5PI) = TEMD2 - DBLE(R(I))
+        SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
+      enddo
+      J = 0
+C
+C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION
+C
+      if (eps.lt.0d0) then
+         goto 24
+      elseif (eps.eq.0d0) then
+         goto 26
+      else
+         goto 27
+      endif
+      
+C
+C INCREMENT DEGREE
+C
+ 16   J = J + 1
+      JP1 = J + 1
+      K1PJ = K1 + J
+      K2PJ = K2 + J
+      SIGJM1 = SIGJ
+C
+C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1
+C
+      IF (J .GT. 1) A(K1PJ) = W11/W1
+C
+C COMPUTE NEW A COEFFICIENT
+C
+      TEMD1 = 0.0D0
+      DO I = 1,M
+        K4PI = K4 + I
+        TEMD2 = A(K4PI)
+        TEMD1 = TEMD1 + DBLE(X(I))*DBLE(W(I))*TEMD2*TEMD2
+      enddo
+      A(JP1) = TEMD1/DBLE(W11)
+C
+C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS
+C
+      W1 = W11
+      W11 = 0.0
+      DO I = 1,M
+        K3PI = K3 + I
+        K4PI = K4 + I
+        TEMP = A(K3PI)
+        A(K3PI) = A(K4PI)
+        A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP
+        W11 = W11 + W(I)*A(K4PI)**2
+      enddo
+C
+C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE
+C PRECISION
+C
+      TEMD1 = 0.0D0
+      DO I = 1,M
+        K4PI = K4 + I
+        K5PI = K5 + I
+        TEMD2 = DBLE(W(I))*DBLE((Y(I)-R(I))-A(K5PI))*DBLE(A(K4PI))
+        TEMD1 = TEMD1 + TEMD2
+      enddo
+      TEMD1 = TEMD1/DBLE(W11)
+      A(K2PJ+1) = TEMD1
+C
+C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND
+C ACCUMULATE SUM OF SQUARES OF ERRORS.  THE POLYNOMIAL EVALUATIONS ARE
+C COMPUTED AND STORED IN EXTENDED PRECISION.  FOR THE I-TH DATA POINT,
+C THE MOST SIGNIFICANT BITS ARE STORED IN  R(I) , AND THE LEAST
+C SIGNIFICANT BITS ARE IN  A(K5PI) .
+C
+      SIGJ = 0.0
+      DO I = 1,M
+        K4PI = K4 + I
+        K5PI = K5 + I
+        TEMD2 = DBLE(R(I)) + DBLE(A(K5PI)) + TEMD1*DBLE(A(K4PI))
+        R(I) = TEMD2
+        A(K5PI) = TEMD2 - DBLE(R(I))
+        SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
+      enddo
+C
+C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE
+C MAXDEG  HAS BEEN REACHED
+C
+      if (eps.lt.0d0) then
+         goto 23
+      elseif (eps.eq.0d0) then
+         goto 26
+      else
+         goto 27
+      endif
+C
+C COMPUTE F STATISTICS  (INPUT EPS .LT. 0.)
+C
+ 23   IF (SIGJ .EQ. 0.0) GO TO 29
+      DEGF = M - J - 1
+      DEN = (CO(4,KSIG)*DEGF + 1.0)*DEGF
+      FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN
+      FCRIT = FCRIT*FCRIT
+      F = (SIGJM1 - SIGJ)*DEGF/SIGJ
+      IF (F .LT. FCRIT) GO TO 25
+C
+C POLYNOMIAL OF DEGREE J SATISFIES F TEST
+C
+ 24   SIGPAS = SIGJ
+      JPAS = J
+      NFAIL = 0
+      IF (MAXDEG .EQ. J) GO TO 32
+      GO TO 16
+C
+C POLYNOMIAL OF DEGREE J FAILS F TEST.  IF THERE HAVE BEEN THREE
+C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND.
+C
+ 25   NFAIL = NFAIL + 1
+      IF (NFAIL .GE. 3) GO TO 29
+      IF (MAXDEG .EQ. J) GO TO 32
+      GO TO 16
+C
+C RAISE THE DEGREE IF DEGREE  MAXDEG  HAS NOT YET BEEN REACHED  (INPUT
+C EPS = 0.)
+C
+ 26   IF (MAXDEG .EQ. J) GO TO 28
+      GO TO 16
+C
+C SEE IF RMS ERROR CRITERION IS SATISFIED  (INPUT EPS .GT. 0.)
+C
+ 27   IF (SIGJ .LE. ETST) GO TO 28
+      IF (MAXDEG .EQ. J) GO TO 31
+      GO TO 16
+C
+C RETURNS
+C
+ 28   IERR = 1
+      NDEG = J
+      SIG = SIGJ
+      GO TO 33
+ 29   IERR = 1
+      NDEG = JPAS
+      SIG = SIGPAS
+      GO TO 33
+ 30   IERR = 2
+      write (*,*) 'SLATEC', 'POLFIT', 'INVALID INPUT PARAMETER.', 2,
+     +     1
+      stop 1
+      GO TO 37
+ 31   IERR = 3
+      NDEG = MAXDEG
+      SIG = SIGJ
+      GO TO 33
+ 32   IERR = 4
+      NDEG = JPAS
+      SIG = SIGPAS
+C
+ 33   A(K3) = NDEG
+C
+C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT
+C ALL THE DATA POINTS IF  R  DOES NOT ALREADY CONTAIN THESE VALUES
+C
+      IF(EPS .GE. 0.0  .OR.  NDEG .EQ. MAXDEG) GO TO 36
+      NDER = 0
+      DO I = 1,M
+        CALL PVALUE (NDEG,NDER,X(I),R(I),YP,A)
+      enddo
+ 36   EPS = SQRT(SIG/XM)
+ 37   RETURN
+      END
+
+
+*DECK PVALUE
+      SUBROUTINE PVALUE (L, NDER, X, YFIT, YP, A)
+c
+cRF: make double precision and replace DO-loops to comply with more
+c    modern standards
+c
+C***BEGIN PROLOGUE  PVALUE
+C***PURPOSE  Use the coefficients generated by POLFIT to evaluate the
+C            polynomial fit of degree L, along with the first NDER of
+C            its derivatives, at a specified point.
+C***LIBRARY   SLATEC
+C***CATEGORY  K6
+C***TYPE      SINGLE PRECISION (PVALUE-S, DP1VLU-D)
+C***KEYWORDS  CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION
+C***AUTHOR  Shampine, L. F., (SNLA)
+C           Davenport, S. M., (SNLA)
+C***DESCRIPTION
+C
+C     Written by L. F. Shampine and S. M. Davenport.
+C
+C     Abstract
+C
+C     The subroutine  PVALUE  uses the coefficients generated by  POLFIT
+C     to evaluate the polynomial fit of degree  L , along with the first
+C     NDER  of its derivatives, at a specified point.  Computationally
+C     stable recurrence relations are used to perform this task.
+C
+C     The parameters for  PVALUE  are
+C
+C     Input --
+C         L -      the degree of polynomial to be evaluated.  L  may be
+C                  any non-negative integer which is less than or equal
+C                  to  NDEG , the highest degree polynomial provided
+C                  by  POLFIT .
+C         NDER -   the number of derivatives to be evaluated.  NDER
+C                  may be 0 or any positive value.  If NDER is less
+C                  than 0, it will be treated as 0.
+C         X -      the argument at which the polynomial and its
+C                  derivatives are to be evaluated.
+C         A -      work and output array containing values from last
+C                  call to  POLFIT .
+C
+C     Output --
+C         YFIT -   value of the fitting polynomial of degree  L  at  X
+C         YP -     array containing the first through  NDER  derivatives
+C                  of the polynomial of degree  L .  YP  must be
+C                  dimensioned at least  NDER  in the calling program.
+C
+C***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
+C                 Curve fitting by polynomials in one variable, Report
+C                 SLA-74-0270, Sandia Laboratories, June 1974.
+C***ROUTINES CALLED  XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   740601  DATE WRITTEN
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890531  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***  END PROLOGUE  PVALUE
+      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
+      DIMENSION YP(*),A(*)
+      CHARACTER*8 XERN1, XERN2
+C***FIRST EXECUTABLE STATEMENT  PVALUE
+      IF (L .LT. 0) GO TO 12
+      NDO = MAX(NDER,0)
+      NDO = MIN(NDO,L)
+      MAXORD = A(1) + 0.5
+      K1 = MAXORD + 1
+      K2 = K1 + MAXORD
+      K3 = K2 + MAXORD + 2
+      NORD = A(K3) + 0.5
+      IF (L .GT. NORD) GO TO 11
+      K4 = K3 + L + 1
+      IF (NDER .LT. 1) GO TO 2
+      DO I = 1,NDER
+         YP(I) = 0.0
+      enddo
+ 2    IF (L .GE. 2) GO TO 4
+      IF (L .EQ. 1) GO TO 3
+C
+C L IS 0
+C
+      VAL = A(K2+1)
+      GO TO 10
+C
+C L IS 1
+C
+ 3    CC = A(K2+2)
+      VAL = A(K2+1) + (X-A(2))*CC
+      IF (NDER .GE. 1) YP(1) = CC
+      GO TO 10
+C
+C L IS GREATER THAN 1
+C
+ 4    NDP1 = NDO + 1
+      K3P1 = K3 + 1
+      K4P1 = K4 + 1
+      LP1 = L + 1
+      LM1 = L - 1
+      ILO = K3 + 3
+      IUP = K4 + NDP1
+      DO I = ILO,IUP
+         A(I) = 0.0
+      enddo
+      DIF = X - A(LP1)
+      KC = K2 + LP1
+      A(K4P1) = A(KC)
+      A(K3P1) = A(KC-1) + DIF*A(K4P1)
+      A(K3+2) = A(K4P1)
+C
+C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES
+C
+      DO I = 1,LM1
+        IN = L - I
+        INP1 = IN + 1
+        K1I = K1 + INP1
+        IC = K2 + IN
+        DIF = X - A(INP1)
+        VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1)
+        IF (NDO .LE. 0) GO TO 8
+        DO N = 1,NDO
+          K3PN = K3P1 + N
+          K4PN = K4P1 + N
+          YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN)
+       enddo
+C
+C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS
+C
+        DO N = 1,NDO
+          K3PN = K3P1 + N
+          K4PN = K4P1 + N
+          A(K4PN) = A(K3PN)
+          A(K3PN) = YP(N)
+        enddo
+ 8      A(K4P1) = A(K3P1)
+        A(K3P1) = VAL
+      enddo
+C
+C NORMAL RETURN OR ABORT DUE TO ERROR
+C
+ 10   YFIT = VAL
+      RETURN
+C
+ 11   WRITE (XERN1, '(I8)') L
+      WRITE (XERN2, '(I8)') NORD
+      write (*,*) 'SLATEC', 'PVALUE',
+     *   'THE ORDER OF POLYNOMIAL EVALUATION, L = ' // XERN1 //
+     *   ' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 //
+     *     ', COMPUTED BY POLFIT -- EXECUTION TERMINATED.', 8, 2
+      stop 1
+      RETURN
+C
+ 12   write (*,*)  'SLATEC', 'PVALUE',
+     +    'INVALID INPUT PARAMETER.  ORDER OF POLYNOMIAL EVALUATION ' //
+     +     'REQUESTED IS NEGATIVE -- EXECUTION TERMINATED.', 2, 2
+      stop 1
+      RETURN
+      END
+      
 
=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2019-08-28 08:02:11 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2019-09-12 11:14:33 +0000
@@ -2376,9 +2376,12 @@
         # latter is only the only line that contains integers.
         for j,fs in enumerate([files_mint_grids,files_MC_integer]):
             linesoffiles=[]
+            polyfit_data=[]
             for f in fs:
                 with open(f,'r+') as fi:
-                    linesoffiles.append(fi.readlines())
+                    data=fi.readlines()
+                    linesoffiles.append([ dat for dat in data if 'POL' not in dat.split()[0] ])
+                    polyfit_data.append([ dat for dat in data if 'POL'     in dat.split()[0] ])
             to_write=[]
             for rowgrp in zip(*linesoffiles):
                 action=list(set([row.strip().split()[0] for row in rowgrp])) # list(set()) structure to remove duplicants
@@ -2414,8 +2417,30 @@
                             write_string.append(int(sum(floatgrp)/len(floatgrp)))
                 else:
                     raise aMCatNLOError('Unknown action for combining grids: %s' % action[0])
-                    
                 to_write.append(action[0] + " " + (" ".join(str(ws) for ws in write_string)) + "\n")
+
+            if polyfit_data:
+                # special for data regarding virtuals. Need to simply append
+                # all the data, but skipping doubles.
+                for i,onefile in enumerate(polyfit_data):
+                    # Get the number of channels, and the number of PS points per channel
+                    data_amount_in_file=[int(oneline.split()[1]) for oneline in onefile if len(oneline.split())==2]
+                    if i==0:
+                        filtered_list=[ [] for i in range(len(data_amount_in_file)) ]
+                    start=len(data_amount_in_file)
+                    for channel,channel_size in enumerate(data_amount_in_file):
+                        end=start+channel_size
+                        for data_channel in onefile[start:end]:
+                            if data_channel not in filtered_list[channel]:
+                                filtered_list[channel].append(data_channel)
+                        start=end
+                # The amount of data in each file (per channel):
+                for channel in filtered_list:
+                    to_write.append("POL " + str(len(channel)) + "\n")
+                # All the data:
+                for ch in filtered_list:
+                    for dat in ch:
+                        to_write.append(dat)
             # write the data over the master location
             if j==0:
                 with open(pjoin(location,'mint_grids'),'w') as f:
 
=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2019-08-15 13:09:08 +0000
+++ madgraph/iolibs/export_fks.py	2019-09-12 11:14:33 +0000
@@ -624,8 +624,7 @@
                      'FKS_params.dat',
                      'initial_states_map.dat',
                      'OLE_order.olc',
-                     'FKSParams.inc',
-                     'FKSParamReader.f',
+                     'FKSParams.f90',
                      'cuts.inc',
                      'unlops.inc',
                      'pythia_unlops.f',
@@ -688,7 +687,8 @@
                      'randinit',
                      'sudakov.inc',
                      'maxconfigs.inc',
-                     'timing_variables.inc']
+                     'timing_variables.inc',
+                     'polfit.f']
 
         for file in linkfiles:
             ln('../' + file , '.')
@@ -1149,24 +1149,26 @@
 
         amp_split_size=len(amp_split_orders)
 
-        text = 'C The orders to be integrated for the Born and at NLO\n'
+        text = '! The orders to be integrated for the Born and at NLO\n'
         text += 'integer nsplitorders\n'
         text += 'parameter (nsplitorders=%d)\n' % len(split_orders)
         text += 'character*3 ordernames(nsplitorders)\n'
         text += 'data ordernames / %s /\n' % ', '.join(['"%3s"' % o for o in split_orders])
         text += 'integer born_orders(nsplitorders), nlo_orders(nsplitorders)\n'
-        text += 'C the order of the coupling orders is %s\n' % ', '.join(split_orders)
+        text += '! the order of the coupling orders is %s\n' % ', '.join(split_orders)
         text += 'data born_orders / %s /\n' % ', '.join([str(max_born_orders[o]) for o in split_orders])
         text += 'data nlo_orders / %s /\n' % ', '.join([str(max_nlo_orders[o]) for o in split_orders])
-        text += 'C The position of the QCD /QED orders in the array\n'
+        text += '! The position of the QCD /QED orders in the array\n'
         text += 'integer qcd_pos, qed_pos\n'
-        text += 'C if = -1, then it is not in the split_orders\n'
+        text += '! if = -1, then it is not in the split_orders\n'
         text += 'parameter (qcd_pos = %d)\n' % qcd_pos
         text += 'parameter (qed_pos = %d)\n' % qed_pos
-        text += 'C this is to keep track of the various coupling combinations entering each ME\n'
+        text += '! this is to keep track of the various \n'
+        text += '! coupling combinations entering each ME\n'
         text += 'integer amp_split_size, amp_split_size_born\n'
         text += 'parameter (amp_split_size = %d)\n' % amp_split_size
-        text += 'parameter (amp_split_size_born = %d) ! the first entries in amp_split are for the born\n' % amp_split_size_born
+        text += '! the first entries in the next line in amp_split are for the born \n'
+        text += 'parameter (amp_split_size_born = %d)\n' % amp_split_size_born
         text += 'double precision amp_split(amp_split_size)\n'
         text += 'double complex amp_split_cnt(amp_split_size,2,nsplitorders)\n'
         text += 'common /to_amp_split/amp_split, amp_split_cnt\n'