C     --=========================================--
C     Main subroutine
C     --=========================================--

      SUBROUTINE MG5_1_SLOOPMATRIX(P_USER,ANS)
C     
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C     
C     Returns amplitude squared summed/avg over colors
C     and helicities for the point in phase space P(0:3,NEXTERNAL)
C     and external lines W(0:6,NEXTERNAL)
C     
C     Process: u u~ > u u~ [ virt = QCD ] @1
C     
      IMPLICIT NONE
      INCLUDE 'coef_specs.inc'
      INCLUDE 'polynomial_specs.inc'
C     
C     USER CUSTOMIZABLE OPTIONS
C     
C     This parameter is designed for the check timing command of MG5.
C      It skips the loop reduction.
      LOGICAL SKIPLOOPEVAL
      PARAMETER (SKIPLOOPEVAL=.FALSE.)
C     For timing checks. Stops the code after having only initialized
C      its arrays from the external data files 
      LOGICAL BOOTANDSTOP
      PARAMETER (BOOTANDSTOP=.FALSE.)
C     
C     CONSTANTS
C     
      CHARACTER*512 PARAMFNAME,HELCONFIGFNAME,LOOPFILTERFNAME
      CHARACTER*512 COLORNUMFNAME,COLORDENOMFNAME, HELFILTERFNAME
      CHARACTER*512 PROC_PREFIX
      PARAMETER ( PARAMFNAME='MadLoopParams.dat')
      PARAMETER ( HELCONFIGFNAME='HelConfigs.dat')
      PARAMETER ( LOOPFILTERFNAME='LoopFilter.dat')
      PARAMETER ( HELFILTERFNAME='HelFilter.dat')
      PARAMETER ( COLORNUMFNAME='ColorNumFactors.dat')
      PARAMETER ( COLORDENOMFNAME='ColorDenomFactors.dat')
      PARAMETER ( PROC_PREFIX='MG5_1_')

      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=2)
      INTEGER    NLOOPS, NLOOPGROUPS, NCTAMPS
      PARAMETER (NLOOPS=22, NLOOPGROUPS=13, NCTAMPS=54)
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=76)
      INTEGER    NCOLORROWS
      PARAMETER (NCOLORROWS=NLOOPAMPS)
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NWAVEFUNCS,NLOOPWAVEFUNCS
      PARAMETER (NWAVEFUNCS=8,NLOOPWAVEFUNCS=50)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)
      REAL*8     ZERO
      PARAMETER (ZERO=0D0)
      REAL*16     MP__ZERO
      PARAMETER (MP__ZERO=0E0_16)
      COMPLEX*16 IMAG1
      PARAMETER (IMAG1=(0D0,1D0))
C     These are constants related to the split orders
      INCLUDE 'nsqso_born.inc'
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=2)
      INTEGER ANS_DIMENSION
      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
      INTEGER NSQSOXNLG
      PARAMETER (NSQSOXNLG=NSQUAREDSO*NLOOPGROUPS)
      INTEGER NSQUAREDSOP1
      PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
C     The total number of loop reduction libraries
C     At present, there are only CutTools,PJFry++,IREGI,Golem95
      INTEGER NLOOPLIB
      PARAMETER (NLOOPLIB=4)
C     Only CutTools provides QP
      INTEGER QP_NLOOPLIB
      PARAMETER (QP_NLOOPLIB=1)
      INTEGER MAXSTABILITYLENGTH
      DATA MAXSTABILITYLENGTH/20/
      COMMON/MG5_1_STABILITY_TESTS/MAXSTABILITYLENGTH

C     
C     ARGUMENTS 
C     
      REAL*8 P_USER(0:3,NEXTERNAL)
C     
C     The zeroth component of the second dimension is the result
C      summed over all
C     contributing split orders. The zeroth component of the first one
C      is the Born.
C     Notice that the upper bound of the second integer is not number
C      of squared orders
C     combination for the loops but the maximum between this number
C      for the Born
C     contributions and the loop ones. There are some cases for which
C      the Born contrib.
C     has squared split order contributions than the loop does. For
C      example
C     
C     generate u u~ > d d~ QCD^2<=2 QED^2<=99 [virt=QCD]
C     
C     It is however somehow academical. This is why ANS_DIMENSION is
C      not just NSQSO but rather MAX(NSQSO,NSQSO_BORN)
C     
      REAL*8 ANS(0:3,0:ANS_DIMENSION)
C     
C     LOCAL VARIABLES 
C     
      INTEGER I,J,K,H,HEL_MULT,I_QP_LIB,DUMMY

      CHARACTER*512 PARAMFN,HELCONFIGFN,LOOPFILTERFN,COLORNUMFN
     $ ,COLORDENOMFN,HELFILTERFN
      CHARACTER*512 TMP
      SAVE PARAMFN
      SAVE HELCONFIGFN
      SAVE LOOPFILTERFN
      SAVE COLORNUMFN
      SAVE COLORDENOMFN
      SAVE HELFILTERFN

      INTEGER CTMODEINIT_BU
      REAL*8 MLSTABTHRES_BU
      INTEGER NEWHELREF
      LOGICAL HEL_INCONSISTENT
      REAL*8 P(0:3,NEXTERNAL)
C     DP_RES STORES THE DOUBLE PRECISION RESULT OBTAINED FROM
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
C     THE STAB_STAGE COUNTER I CORRESPONDANCE GOES AS FOLLOWS
C     I=1 -> ORIGINAL PS, CTMODE=1
C     I=2 -> ORIGINAL PS, CTMODE=2, (ONLY WITH CTMODERUN=-1)
C     I=3 -> PS WITH ROTATION 1, CTMODE=1, (ONLY WITH CTMODERUN=-2)
C     I=4 -> PS WITH ROTATION 2, CTMODE=1, (ONLY WITH CTMODERUN=-3)
C     I=5 -> POSSIBLY MORE EVALUATION METHODS IN THE FUTURE, MAX IS
C      MAXSTABILITYLENGTH
C     IF UNSTABLE IT GOES TO THE SAME PATTERN BUT STAB_INDEX IS THEN
C      I+20.
      LOGICAL EVAL_DONE(MAXSTABILITYLENGTH)
      LOGICAL DOING_QP_EVALS
      INTEGER STAB_INDEX,BASIC_CT_MODE


      REAL*8 ACC(0:NSQUAREDSO)
      REAL*8 DP_RES(3,0:NSQUAREDSO,MAXSTABILITYLENGTH)
C     QP_RES STORES THE QUADRUPLE PRECISION RESULT OBTAINED FROM
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
      REAL*8 QP_RES(3,0:NSQUAREDSO,MAXSTABILITYLENGTH)
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      INTEGER NATTEMPTS
      DATA NATTEMPTS/0/
      DATA IC/NEXTERNAL*1/
      REAL*8 HELSAVED(3,NCOMB)
      INTEGER ITEMP
      LOGICAL LTEMP
      REAL*8 BORNBUFF(0:NSQSO_BORN),TMPR
      REAL*8 BUFFR(3,0:NSQUAREDSO),BUFFR_BIS(3,0:NSQUAREDSO),TEMP(0:3
     $ ,0:NSQUAREDSO),TEMP1(0:NSQUAREDSO)
      REAL*8 TEMP2
      COMPLEX*16 COEFS(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE)
      COMPLEX*16 CFTOT
      LOGICAL FOUNDHELFILTER,FOUNDLOOPFILTER
      DATA FOUNDHELFILTER/.TRUE./
      DATA FOUNDLOOPFILTER/.TRUE./
      LOGICAL LOOPFILTERBUFF(NSQUAREDSO,NLOOPGROUPS)
      DATA ((LOOPFILTERBUFF(J,I),J=1,NSQUAREDSO),I=1,NLOOPGROUPS)
     $ /NSQSOXNLG*.FALSE./

      LOGICAL AUTOMATIC_TIR_CACHE_CLEARING
      DATA AUTOMATIC_TIR_CACHE_CLEARING/.TRUE./
      COMMON/MG5_1_RUNTIME_OPTIONS/AUTOMATIC_TIR_CACHE_CLEARING

      INTEGER IDEN
      DATA IDEN/36/
      INTEGER HELAVGFACTOR
      DATA HELAVGFACTOR/4/
      LOGICAL DONEHELDOUBLECHECK
      DATA DONEHELDOUBLECHECK/.FALSE./
      INTEGER NEPS
      DATA NEPS/0/
C     Below are variables to bypass the checkphase and insure
C      stability check to take place
      LOGICAL OLD_CHECKPHASE, OLD_HELDOUBLECHECKED
      INTEGER OLD_GOODHEL(NCOMB)
      LOGICAL OLD_GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/MG5_1_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
C     
C     FUNCTIONS
C     
      INTEGER MG5_1_TIRCACHE_INDEX
      INTEGER MG5_1_ML5SOINDEX_FOR_BORN_AMP
      INTEGER MG5_1_ML5SOINDEX_FOR_LOOP_AMP
      INTEGER MG5_1_ML5SQSOINDEX
      INTEGER MG5_1_ISSAME
      LOGICAL MG5_1_ISZERO
      INTEGER SET_RET_CODE_U
C     
C     GLOBAL VARIABLES
C     
      INCLUDE 'coupl.inc'
      INCLUDE 'mp_coupl.inc'
      INCLUDE 'MadLoopParams.inc'


      LOGICAL CHOSEN_SO_CONFIGS(NSQUAREDSO)
      DATA CHOSEN_SO_CONFIGS/.TRUE./
      COMMON/MG5_1_CHOSEN_LOOP_SQSO/CHOSEN_SO_CONFIGS

      INTEGER N_DP_EVAL, N_QP_EVAL
      DATA N_DP_EVAL/1/
      DATA N_QP_EVAL/1/
      COMMON/MG5_1_N_EVALS/N_DP_EVAL,N_QP_EVAL

      LOGICAL CHECKPHASE
      DATA CHECKPHASE/.TRUE./
      LOGICAL HELDOUBLECHECKED
      DATA HELDOUBLECHECKED/.FALSE./
      COMMON/MG5_1_INIT/CHECKPHASE, HELDOUBLECHECKED
      INTEGER NTRY
      DATA NTRY/0/
      REAL*8 REF
      DATA REF/0.0D0/

      LOGICAL MP_DONE
      DATA MP_DONE/.FALSE./
      COMMON/MG5_1_MP_DONE/MP_DONE
C     A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE
C      AVAILABLE OR NOT
      LOGICAL LOOPLIBS_AVAILABLE(4)
      DATA LOOPLIBS_AVAILABLE/.TRUE.,.FALSE.,.TRUE.,.TRUE./
      COMMON/MG5_1_LOOPLIBS_AV/ LOOPLIBS_AVAILABLE
C     A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS
C      AVAILABLE OR NOT IN THE LOOPLIBS
C     PJFry++ and Golem95 do not support direction test
      LOGICAL LOOPLIBS_DIRECTEST(4)
      DATA LOOPLIBS_DIRECTEST /.TRUE.,.TRUE.,.TRUE.,.TRUE./

C     PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT
C      MODIFIED FOR THE PURPOSE OF THE STABILITY TEST
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
C      USED ANYWHERE ELSE
      REAL*8 PS(0:3,NEXTERNAL)
      COMMON/MG5_1_PSPOINT/PS
C     AGAIN BELOW, MP_PS IS THE FIXED (POSSIBLY IMPROVED) MP PS POINT
C      AND MP_P IS THE ONE WHICH CAN BE MODIFIED (I.E. ROTATED ETC.)
C      FOR STABILITY PURPOSE
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
      COMMON/MG5_1_MP_PSPOINT/MP_PS,MP_P

      REAL*8 LSCALE
      INTEGER CTMODE
      COMMON/MG5_1_CT/LSCALE,CTMODE
      LOGICAL MP_PS_SET
      DATA MP_PS_SET/.FALSE./

C     The parameter below sets the convention for the helicity filter
C     For a given helicity, the attached integer 'i' means
C     'i' in ]-inf;-HELOFFSET[ -> Helicity is equal, up to a sign, to
C      helicity number abs(i+HELOFFSET)
C     'i' == -HELOFFSET        -> Helicity is analytically zero
C     'i' in ]-HELOFFSET,inf[  -> Helicity is contributing with weight
C      'i'. If it is zero, it is skipped.
C     Typically, the hel_offset is 10000
      INTEGER HELOFFSET
      DATA HELOFFSET/10000/
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      COMMON/MG5_1_FILTERS/GOODAMP,GOODHEL,HELOFFSET

      INTEGER HELPICKED
      DATA HELPICKED/-1/
      COMMON/MG5_1_HELCHOICE/HELPICKED
      INTEGER USERHEL
      DATA USERHEL/-1/
      COMMON/MG5_1_USERCHOICE/USERHEL

C     This integer can be accessed by an external user to set its
C      target squared split order.
C     If set to a value different than -1, the code will try to avoid
C      computing anything which
C     does not contribute to contributions of squared split orders
C      SQSO_TARGET and below.
      INTEGER SQSO_TARGET
      DATA SQSO_TARGET/-1/
      COMMON/MG5_1_SOCHOICE/SQSO_TARGET
C     The following logical are used to broadcast the fact that the
C      target 'required' CT and
C     loop split orders contributions have been reached already and
C      the rest can be skipped.
      LOGICAL UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE
     $ ,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE
     $ ,CTCALL_REQ_SO_DONE,FILTER_SO
      DATA UVCT_REQ_SO_DONE/.FALSE./
      DATA MP_UVCT_REQ_SO_DONE/.FALSE./
      DATA CT_REQ_SO_DONE/.FALSE./
      DATA MP_CT_REQ_SO_DONE/.FALSE./
      DATA LOOP_REQ_SO_DONE/.FALSE./
      DATA MP_LOOP_REQ_SO_DONE/.FALSE./
      DATA CTCALL_REQ_SO_DONE/.FALSE./
      DATA FILTER_SO/.FALSE./
      COMMON/MG5_1_SO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE
     $ ,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_S
     $ O_DONE,CTCALL_REQ_SO_DONE,FILTER_SO

C     Allows to forbid the zero helicity double check, no matter the
C      value in MadLoopParams.dat
C     This can be accessed with the SET_FORBID_HEL_DOUBLECHECK
C      subroutine of MadLoopCommons.dat
      LOGICAL FORBID_HEL_DOUBLECHECK
      COMMON/FORBID_HEL_DOUBLECHECK/FORBID_HEL_DOUBLECHECK

      INTEGER I_SO
      DATA I_SO/1/
      COMMON/MG5_1_I_SO/I_SO
      INTEGER I_LIB
      DATA I_LIB/1/
      COMMON/MG5_1_I_LIB/I_LIB
C     UTIL NOW, ONLY CUTTOOLS PROVIDE QP
      LOGICAL QP_TOOLS_AVAILABLE
      DATA QP_TOOLS_AVAILABLE/.FALSE./
      INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
      COMMON/MG5_1_LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS

      COMPLEX*16 AMP(NBORNAMPS)
      COMMON/MG5_1_AMPS/AMP
      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/MG5_1_W/W

      COMPLEX*32 MPW(20,NWAVEFUNCS)
      COMMON/MG5_1_MP_W/MPW

      COMPLEX*16 WL(MAXLWFSIZE,0:LOOPMAXCOEFS-1,MAXLWFSIZE,0:NLOOPWAVEF
     $ UNCS)
      COMPLEX*16 PL(0:3,0:NLOOPWAVEFUNCS)
      COMMON/MG5_1_WL/WL,PL

      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
      COMMON/MG5_1_LCOEFS/LOOPCOEFS



      COMPLEX*16 AMPL(3,NCTAMPS)
      COMMON/MG5_1_AMPL/AMPL

      COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
      LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
      COMMON/MG5_1_LOOPRES/LOOPRES,S

      INTEGER CF_D(NCOLORROWS,NBORNAMPS)
      INTEGER CF_N(NCOLORROWS,NBORNAMPS)
      COMMON/MG5_1_CF/CF_D,CF_N

      INTEGER HELC(NEXTERNAL,NCOMB)
      COMMON/MG5_1_HELCONFIGS/HELC

      REAL*8 PREC,USER_STAB_PREC
      DATA USER_STAB_PREC/-1.0D0/
      COMMON/MG5_1_USER_STAB_PREC/USER_STAB_PREC

C     Return codes H,T,U correspond to the hundreds, tens and units
C     building returncode, i.e.
C     RETURNCODE=100*RET_CODE_H+10*RET_CODE_T+RET_CODE_U

      INTEGER RET_CODE_H,RET_CODE_T,RET_CODE_U
      REAL*8 ACCURACY(0:NSQUAREDSO)
      DATA (ACCURACY(I),I=0,NSQUAREDSO)/NSQUAREDSOP1*1.0D0/
      DATA RET_CODE_H,RET_CODE_T,RET_CODE_U/1,1,0/
      COMMON/MG5_1_ACC/ACCURACY,RET_CODE_H,RET_CODE_T,RET_CODE_U

      LOGICAL MP_DONE_ONCE
      DATA MP_DONE_ONCE/.FALSE./
      COMMON/MG5_1_MP_DONE_ONCE/MP_DONE_ONCE

      CHARACTER(512) MLPATH
      COMMON/MLPATH/MLPATH

C     This variable controls the general initialization which is
C      *common* between all MadLoop SubProcesses.
C     For example setting the MadLoopPath or reading the ML runtime
C      parameters.
      LOGICAL ML_INIT
      COMMON/ML_INIT/ML_INIT

C     This variable controls the *local* initialization of this
C      particular SubProcess.
C     For example, the reading of the filters must be done independentl
C     y by each SubProcess.
      LOGICAL LOCAL_ML_INIT
      DATA LOCAL_ML_INIT/.TRUE./

C     ----------
C     BEGIN CODE
C     ----------

      IF(ML_INIT) THEN
        ML_INIT = .FALSE.
        CALL PRINT_MADLOOP_BANNER()
        TMP = 'auto'
        CALL SETMADLOOPPATH(TMP)
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
        CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.)
        IF (FORBID_HEL_DOUBLECHECK) THEN
          DOUBLECHECKHELICITYFILTER = .FALSE.
        ENDIF

C       Make sure that NROTATIONS_QP and NROTATIONS_DP are set to zero
C        if AUTOMATIC_TIR_CACHE_CLEARING is disabled.
        IF(.NOT.AUTOMATIC_TIR_CACHE_CLEARING) THEN
          IF(NROTATIONS_DP.NE.0.OR.NROTATIONS_QP.NE.0) THEN
            WRITE(*,*) '##INFO: AUTOMATIC_TIR_CACHE_CLEARING i'
     $       //'s disabled, so MadLoop automatically resets NROTATIONS'
     $       //'_DP and NROTATIONS_QP to 0.'
            NROTATIONS_QP=0
            NROTATIONS_DP=0
          ENDIF
        ENDIF
      ENDIF

      IF (LOCAL_ML_INIT) THEN
        LOCAL_ML_INIT = .FALSE.
        QP_TOOLS_AVAILABLE=.FALSE.
        INDEX_QP_TOOLS(1:QP_NLOOPLIB+1)=0
C       SKIP THE ONES THAT NOT AVAILABLE
        J=1
        DO I=1,NLOOPLIB
          IF(MLREDUCTIONLIB(J).EQ.0)EXIT
          IF(.NOT.LOOPLIBS_AVAILABLE(MLREDUCTIONLIB(J)))THEN
            MLREDUCTIONLIB(J:NLOOPLIB-1)=MLREDUCTIONLIB(J+1:NLOOPLIB)
            MLREDUCTIONLIB(NLOOPLIB)=0
          ELSE
            J=J+1
          ENDIF
        ENDDO
        IF(MLREDUCTIONLIB(1).EQ.0)THEN
          STOP 'No available loop reduction lib is provided. Make sur'
     $     //'e MLReductionLib is correct.'
        ENDIF
        J=0
        DO I=1,NLOOPLIB
          IF(MLREDUCTIONLIB(I).EQ.1)THEN
            J=J+1
            IF(.NOT.QP_TOOLS_AVAILABLE)QP_TOOLS_AVAILABLE=.TRUE.
            INDEX_QP_TOOLS(J)=I
          ENDIF
        ENDDO

C       Setup the file paths
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
        CALL JOINPATH(MLPATH,PROC_PREFIX,TMP)
        CALL JOINPATH(TMP,HELCONFIGFNAME,HELCONFIGFN)
        CALL JOINPATH(TMP,LOOPFILTERFNAME,LOOPFILTERFN)
        CALL JOINPATH(TMP,COLORNUMFNAME,COLORNUMFN)
        CALL JOINPATH(TMP,COLORDENOMFNAME,COLORDENOMFN)
        CALL JOINPATH(TMP,HELFILTERFNAME,HELFILTERFN)

        CALL MG5_1_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)

        OPEN(1, FILE=COLORNUMFN, ERR=104, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NCOLORROWS
          READ(1,*,END=105) (CF_N(I,J),J=1,NBORNAMPS)
        ENDDO
        GOTO 105
 104    CONTINUE
        STOP 'Color factors could not be initialized from fil'
     $   //'e MG5_1_ColorNumFactors.dat. File not found'
 105    CONTINUE
        CLOSE(1)
        OPEN(1, FILE=COLORDENOMFN, ERR=106, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NCOLORROWS
          READ(1,*,END=107) (CF_D(I,J),J=1,NBORNAMPS)
        ENDDO
        GOTO 107
 106    CONTINUE
        STOP 'Color factors could not be initialized from fil'
     $   //'e MG5_1_ColorDenomFactors.dat. File not found'
 107    CONTINUE
        CLOSE(1)
        OPEN(1, FILE=HELCONFIGFN, ERR=108, STATUS='OLD',              
     $       ACTION='READ')
        DO H=1,NCOMB
          READ(1,*,END=109) (HELC(I,H),I=1,NEXTERNAL)
        ENDDO
        GOTO 109
 108    CONTINUE
        STOP 'Color helictiy configurations could not be initialize'
     $   //'d from file MG5_1_HelConfigs.dat. File not found'
 109    CONTINUE
        CLOSE(1)

C       SETUP OF THE COMMON STARTING EXTERNAL LOOP WAVEFUNCTION
C       IT IS ALSO PS POINT INDEPENDENT, SO IT CAN BE DONE HERE.
        DO I=0,3
          PL(I,0)=(0.0D0,0.0D0)
        ENDDO
        DO I=1,MAXLWFSIZE
          DO J=0,LOOPMAXCOEFS-1
            DO K=1,MAXLWFSIZE
              IF(I.EQ.K.AND.J.EQ.0) THEN
                WL(I,J,K,0)=(1.0D0,0.0D0)
              ELSE
                WL(I,J,K,0)=(0.0D0,0.0D0)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
        IF(BOOTANDSTOP) THEN
          WRITE(*,*) '##Stopped by user request.'
          STOP
        ENDIF
      ENDIF

      IF(NTRY.EQ.0) THEN
        HELDOUBLECHECKED=(.NOT.DOUBLECHECKHELICITYFILTER).OR.(HELICITYF
     $   ILTERLEVEL.EQ.0)
        OPEN(1, FILE=LOOPFILTERFN, ERR=100, STATUS='OLD',          
     $    ACTION='READ')
        DO J=1,NLOOPGROUPS
          READ(1,*,END=101) (GOODAMP(I,J),I=1,NSQUAREDSO)
        ENDDO
        GOTO 101
 100    CONTINUE
        FOUNDLOOPFILTER=.FALSE.
        DO J=1,NLOOPGROUPS
          DO I=1,NSQUAREDSO
            GOODAMP(I,J)=(.NOT.USELOOPFILTER)
          ENDDO
        ENDDO
 101    CONTINUE
        CLOSE(1)

        IF (HELICITYFILTERLEVEL.EQ.0) THEN
          FOUNDHELFILTER=.TRUE.
          DO J=1,NCOMB
            GOODHEL(J)=1
          ENDDO
          GOTO 122
        ENDIF
        OPEN(1, FILE=HELFILTERFN, ERR=102, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NCOMB
          READ(1,*,END=103) GOODHEL(I)
        ENDDO
        GOTO 103
 102    CONTINUE
        FOUNDHELFILTER=.FALSE.
        DO J=1,NCOMB
          GOODHEL(J)=1
        ENDDO
 103    CONTINUE
        CLOSE(1)
        IF (HELICITYFILTERLEVEL.EQ.1) THEN
C         We must make sure to remove the matching-helicity optimisatio
C         n, as requested by the user.
          DO J=1,NCOMB
            IF ((GOODHEL(J).GT.1).OR.(GOODHEL(J).LT.-HELOFFSET)) THEN
              GOODHEL(J)=1
            ENDIF
          ENDDO
        ENDIF
 122    CONTINUE
      ENDIF

C     First compute the borns, it will store them in ANS(0,I)
C     It is left untouched for the rest of MadLoop evaluation.
C     Notice that the squared split order index I does NOT
C     correspond to the same ordering of J for the loop ME 
C     results stored in ANS(K,J), with K in [1-3].The ordering 
C     of each can be obtained with ML5SOINDEX_FOR_SQUARED_ORDERS 
C     and SQSOINDEX_FROM_ORDERS for the loop ME and born ME 
C     respectively. For this to work, we assume that there is 
C     always more squared split orders in the loop ME than in the
C     born ME, which is practically always true. In any case, only
C     the split_order summed value I=0 is used in ML5 code.
      DO I=0,NSQSO_BORN
        BORNBUFF(I)=0.0D0
      ENDDO
      CALL MG5_1_SMATRIXHEL_SPLITORDERS(P_USER,USERHEL,BORNBUFF(0))
      DO I=0,NSQSO_BORN
        ANS(0,I)=BORNBUFF(I)
      ENDDO

C     We set here the reference to the born summed over all split
C      orders
      REF=0.0D0
      DO I=1,NSQSO_BORN
        REF=REF+ANS(0,I)
      ENDDO

      MP_DONE=.FALSE.
      MP_DONE_ONCE=.FALSE.
      MP_PS_SET=.FALSE.
      STAB_INDEX=0
      DOING_QP_EVALS=.FALSE.
      EVAL_DONE(1)=.TRUE.
      DO I=2,MAXSTABILITYLENGTH
        EVAL_DONE(I)=.FALSE.
      ENDDO

      IF(.NOT.BYPASS_CHECK) THEN
        NTRY=NTRY+1
      ENDIF

      IF (USER_STAB_PREC.GT.0.0D0) THEN
        MLSTABTHRES_BU=MLSTABTHRES
        MLSTABTHRES=USER_STAB_PREC
C       In the initialization, I cannot perform stability test and
C        therefore guarantee any precision
        CTMODEINIT_BU=CTMODEINIT
C       So either one choses quad precision directly
C       CTMODEINIT=4
C       Or, because this is very slow, we keep the orignal value. The
C        accuracy returned is -1 and tells the MC that he should not
C        trust the evaluation for checks.
        CTMODEINIT=CTMODEINIT_BU
      ENDIF

      IF(DONEHELDOUBLECHECK.AND.(.NOT.HELDOUBLECHECKED)) THEN
        HELDOUBLECHECKED=.TRUE.
        DONEHELDOUBLECHECK=.FALSE.
      ENDIF

      CHECKPHASE=(NTRY.LE.CHECKCYCLE).AND.(((.NOT.FOUNDLOOPFILTER
     $ ).AND.USELOOPFILTER).OR.(.NOT.FOUNDHELFILTER))

      IF (WRITEOUTFILTERS) THEN
        IF ((HELICITYFILTERLEVEL.NE.0).AND.(.NOT. CHECKPHASE).AND.(.NOT
     $   .FOUNDHELFILTER)) THEN
          OPEN(1, FILE=HELFILTERFN, ERR=110, STATUS='NEW',ACTION='WRIT'
     $     //'E')
          DO I=1,NCOMB
            WRITE(1,*) GOODHEL(I)
          ENDDO
 110      CONTINUE
          CLOSE(1)
          FOUNDHELFILTER=.TRUE.
        ENDIF

        IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDLOOPFILTER).AND.USELOOPFI
     $   LTER) THEN
          OPEN(1, FILE=LOOPFILTERFN, ERR=111, STATUS='NEW',ACTION='WRI'
     $     //'TE')
          DO J=1,NLOOPGROUPS
            WRITE(1,*) (GOODAMP(I,J),I=1,NSQUAREDSO)
          ENDDO
 111      CONTINUE
          CLOSE(1)
          FOUNDLOOPFILTER=.TRUE.
        ENDIF
      ENDIF

      IF (BYPASS_CHECK) THEN
        OLD_CHECKPHASE = CHECKPHASE
        OLD_HELDOUBLECHECKED = HELDOUBLECHECKED
        CHECKPHASE = .FALSE.
        HELDOUBLECHECKED = .TRUE.
        DO I=1,NCOMB
          OLD_GOODHEL(I)=GOODHEL(I)
          GOODHEL(I)=1
        ENDDO
        DO I=1,NSQUAREDSO
          DO J=1,NLOOPGROUPS
            OLD_GOODAMP(I,J)=GOODAMP(I,J)
            GOODAMP(I,J)=.TRUE.
          ENDDO
        ENDDO
      ENDIF

      IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
        HELPICKED=1
        CTMODE=CTMODEINIT
      ELSE
        IF (USERHEL.NE.-1) THEN
          IF(GOODHEL(USERHEL).EQ.-HELOFFSET) THEN
            DO I=0,NSQUAREDSO
              ANS(1,I)=0.0D0
              ANS(2,I)=0.0D0
              ANS(3,I)=0.0D0
            ENDDO
            GOTO 9999
          ENDIF
        ENDIF
        HELPICKED=USERHEL
        IF (CTMODERUN.NE.-1) THEN
          CTMODE=CTMODERUN
        ELSE
          CTMODE=1
        ENDIF
      ENDIF

      DO I=1,NEXTERNAL
        DO J=0,3
          PS(J,I)=P_USER(J,I)
        ENDDO
      ENDDO


      IF (IMPROVEPSPOINT.GE.0) THEN
C       Make the input PS more precise (exact onshell and energy-moment
C       um conservation)
        CALL MG5_1_IMPROVE_PS_POINT_PRECISION(PS)
      ENDIF

      DO I=1,NEXTERNAL
        DO J=0,3
          P(J,I)=PS(J,I)
        ENDDO
      ENDDO

      DO K=1, 3
        DO I=0,NSQUAREDSO
          BUFFR(K,I)=0.0D0
        ENDDO
        DO I=1,NCTAMPS
          AMPL(K,I)=(0.0D0,0.0D0)
        ENDDO
      ENDDO

C     Start by using the first available loop reduction library and qp
C      library.
      I_LIB=1
      I_QP_LIB=1

      GOTO 208
C     MadLoop jumps to this label during stability checks when it
C      recomputes a rotated PS point
 200  CONTINUE
C     For the computation of a rotated version of this PS point we
C      must reset the TIR cache since this changes the definition of
C      the loop denominators.
      CALL MG5_1_CLEAR_TIR_CACHE()
 208  CONTINUE

C     MadLoop jumps to this label during initialization when it goes
C      to the computation of the next helicity.
 205  CONTINUE

      IF (.NOT.MP_PS_SET.AND.(CTMODE.EQ.0.OR.CTMODE.GE.4)) THEN
        CALL MG5_1_SET_MP_PS(P_USER)
        MP_PS_SET = .TRUE.
      ENDIF

      LSCALE=DSQRT(ABS((P(0,1)+P(0,2))**2-(P(1,1)+P(1,2))**2-(P(2,1)
     $ +P(2,2))**2-(P(3,1)+P(3,2))**2))

      CTCALL_REQ_SO_DONE=.FALSE.
      FILTER_SO = (.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.(SQSO_TARG
     $ ET.NE.-1)

      DO I=1,NLOOPGROUPS
        DO J=0,LOOPMAXCOEFS-1
          DO K=1,NSQUAREDSO
            LOOPCOEFS(J,K,I)=(0.0D0,0.0D0)
          ENDDO
        ENDDO
      ENDDO

      DO I=1,NLOOPGROUPS
        DO J=1,3
          DO K=1,NSQUAREDSO
            LOOPRES(J,K,I)=(0.0D0,0.0D0)
          ENDDO
        ENDDO
      ENDDO

      DO K=1,3
        DO I=0,NSQUAREDSO
          ANS(K,I)=0.0D0
        ENDDO
      ENDDO

C     Check if we directly go to multiple precision
      IF (CTMODE.GE.4) THEN
        IF (.NOT.MP_DONE) THEN
          CALL MG5_1_MP_COMPUTE_LOOP_COEFS(MP_P,BUFFR_BIS)
C         It should be safe to directly set MP_DONE to true already
C          here. But maybe I overlooked something.
          MP_DONE=.TRUE.
        ENDIF
C       Even if MP_DONE is .TRUE. we should anyway skip the
C       double precision evaluation as it as already been
C       computed in quadruple precision.
        GOTO 300
      ENDIF

      DO H=1,NCOMB
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
     $   .NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHE
     $   L(H).NE.0)))) THEN
          DO I=1,NEXTERNAL
            NHEL(I)=HELC(I,H)
          ENDDO

          UVCT_REQ_SO_DONE=.FALSE.
          CT_REQ_SO_DONE=.FALSE.
          LOOP_REQ_SO_DONE=.FALSE.

          IF (.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.HELPICKED.EQ.
     $     -1) THEN
            HEL_MULT=GOODHEL(H)
          ELSE
            HEL_MULT=1
          ENDIF


C         Helas calls for the born amplitudes and counterterms
C          associated to given loops
          CALL IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1))
          CALL OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
          CALL OXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3))
          CALL IXXXXX(P(0,4),ZERO,NHEL(4),-1*IC(4),W(1,4))
          CALL FFV1P0_3(W(1,1),W(1,2),GC_5,ZERO,ZERO,W(1,5))
C         Amplitude(s) for born diagram with ID 1
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),GC_5,AMP(1))
          CALL FFV1P0_3(W(1,1),W(1,3),GC_5,ZERO,ZERO,W(1,6))
C         Amplitude(s) for born diagram with ID 2
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),GC_5,AMP(2))
          CALL FFV1P0_3(W(1,4),W(1,3),GC_5,ZERO,ZERO,W(1,7))
C         Counter-term amplitude(s) for loop diagram number 3
          CALL R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,1))
          CALL R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,2))
          CALL R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,3))
          CALL R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,4))
          CALL FFV1P0_3(W(1,4),W(1,2),GC_5,ZERO,ZERO,W(1,8))
C         Counter-term amplitude(s) for loop diagram number 4
          CALL R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,5))
          CALL R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,6))
          CALL R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,7))
          CALL R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,8))
C         Counter-term amplitude(s) for loop diagram number 5
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQQ_1EPS,AMPL(2,9))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQQ_1EPS,AMPL(2,10))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQQ_1EPS,AMPL(2,11))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQQ_1EPS,AMPL(2,12))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB,AMPL(1,13))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQQ_1EPS,AMPL(2,14))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQT,AMPL(1,15))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQQ_1EPS,AMPL(2,16))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQG_1EPS,AMPL(2,17))
          CALL FFV1_0(W(1,4),W(1,3),W(1,5),R2_GQQ,AMPL(1,18))
C         Counter-term amplitude(s) for loop diagram number 7
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQQ_1EPS,AMPL(2,19))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQQ_1EPS,AMPL(2,20))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQQ_1EPS,AMPL(2,21))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQQ_1EPS,AMPL(2,22))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB,AMPL(1,23))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQQ_1EPS,AMPL(2,24))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQT,AMPL(1,25))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQQ_1EPS,AMPL(2,26))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQG_1EPS,AMPL(2,27))
          CALL FFV1_0(W(1,4),W(1,2),W(1,6),R2_GQQ,AMPL(1,28))
C         Counter-term amplitude(s) for loop diagram number 9
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQQ_1EPS,AMPL(2,29))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQQ_1EPS,AMPL(2,30))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQQ_1EPS,AMPL(2,31))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQQ_1EPS,AMPL(2,32))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB,AMPL(1,33))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQQ_1EPS,AMPL(2,34))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQT,AMPL(1,35))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQQ_1EPS,AMPL(2,36))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQG_1EPS,AMPL(2,37))
          CALL FFV1_0(W(1,1),W(1,3),W(1,8),R2_GQQ,AMPL(1,38))
C         Counter-term amplitude(s) for loop diagram number 13
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQQ_1EPS,AMPL(2,39))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQQ_1EPS,AMPL(2,40))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQQ_1EPS,AMPL(2,41))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQQ_1EPS,AMPL(2,42))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB,AMPL(1,43))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQQ_1EPS,AMPL(2,44))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQT,AMPL(1,45))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQQ_1EPS,AMPL(2,46))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQG_1EPS,AMPL(2,47))
          CALL FFV1_0(W(1,1),W(1,2),W(1,7),R2_GQQ,AMPL(1,48))
C         Counter-term amplitude(s) for loop diagram number 17
          CALL R2_GG_1_R2_GG_3_0(W(1,5),W(1,7),R2_GGQ,R2_GGB,AMPL(1
     $     ,49))
C         Counter-term amplitude(s) for loop diagram number 18
          CALL R2_GG_1_R2_GG_3_0(W(1,6),W(1,8),R2_GGQ,R2_GGB,AMPL(1
     $     ,50))
C         Counter-term amplitude(s) for loop diagram number 19
          CALL R2_GG_1_R2_GG_3_0(W(1,5),W(1,7),R2_GGQ,R2_GGT,AMPL(1
     $     ,51))
C         Counter-term amplitude(s) for loop diagram number 20
          CALL R2_GG_1_R2_GG_3_0(W(1,6),W(1,8),R2_GGQ,R2_GGT,AMPL(1
     $     ,52))
C         Counter-term amplitude(s) for loop diagram number 21
          CALL R2_GG_1_R2_GG_2_0(W(1,5),W(1,7),R2_GGG_1,R2_GGG_2
     $     ,AMPL(1,53))
C         Counter-term amplitude(s) for loop diagram number 22
          CALL R2_GG_1_R2_GG_2_0(W(1,6),W(1,8),R2_GGG_1,R2_GGG_2
     $     ,AMPL(1,54))
C         At this point, all CT amps needed for (QCD=6), i.e. of split
C          order ID=1, are computed.
          IF(FILTER_SO.AND.SQSO_TARGET.EQ.1) GOTO 2000
 2000     CONTINUE
          CT_REQ_SO_DONE=.TRUE.

C         Helas calls for the counterterm of type 'UVtree' in the UFO.
C         These are generated irrespectively of the produced loops.
C         In general, only wavefunction renormalization counterterms
C         (if needed by the loop UFO model) are of this type.
C         Quite often and in principle for all loop UFO models from 
C         FeynRules, there are none of these type of counterterms.

 3000     CONTINUE
          UVCT_REQ_SO_DONE=.TRUE.

          DO I=1,NCTAMPS
            DO J=1,NBORNAMPS
              CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0)
              IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1
              ITEMP = MG5_1_ML5SQSOINDEX(MG5_1_ML5SOINDEX_FOR_LOOP_AMP(
     $         I),MG5_1_ML5SOINDEX_FOR_BORN_AMP(J))
              IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.ITEMP) THEN
                DO K=1,3
                  TEMP2 = 2.0D0*HEL_MULT*DBLE(CFTOT*AMPL(K,I)
     $             *DCONJG(AMP(J)))
                  ANS(K,ITEMP)=ANS(K,ITEMP)+TEMP2
                  ANS(K,0)=ANS(K,0)+TEMP2
                ENDDO
              ENDIF
            ENDDO
          ENDDO

C         Coefficient construction for loop diagram with ID 3
          CALL FFV1L2_1(PL(0,0),W(1,5),GC_5,ZERO,ZERO,PL(0,1),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1,1))
          CALL FFV1L2_1(PL(0,1),W(1,7),GC_5,ZERO,ZERO,PL(0,2),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,1),4,COEFS,4,4,WL(1,0,1,2))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,2),2,4,1,1,4,55,H)
C         Coefficient construction for loop diagram with ID 4
          CALL FFV1L2_1(PL(0,0),W(1,6),GC_5,ZERO,ZERO,PL(0,3),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1,3))
          CALL FFV1L2_1(PL(0,3),W(1,8),GC_5,ZERO,ZERO,PL(0,4),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,3),4,COEFS,4,4,WL(1,0,1,4))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,4),2,4,2,1,4,56,H)
C         Coefficient construction for loop diagram with ID 5
          CALL FFV1L1P0_3(PL(0,0),W(1,3),GC_5,ZERO,ZERO,PL(0,5),COEFS)
          CALL MG5_1_UPDATE_WL_0_0(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1,5))
          CALL FFV1L3_2(PL(0,5),W(1,4),GC_5,ZERO,ZERO,PL(0,6),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,5),4,COEFS,4,4,WL(1,0,1,6))
          CALL FFV1L1_2(PL(0,6),W(1,5),GC_5,ZERO,ZERO,PL(0,7),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,6),4,COEFS,4,4,WL(1,0,1,7))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,7),2,4,3,1,1,57,H)
C         Coefficient construction for loop diagram with ID 6
          CALL FFV1L3_1(PL(0,0),W(1,3),GC_5,ZERO,ZERO,PL(0,8),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1,8))
          CALL FFV1L2P0_3(PL(0,8),W(1,4),GC_5,ZERO,ZERO,PL(0,9),COEFS)
          CALL MG5_1_UPDATE_WL_1_0(WL(1,0,1,8),4,COEFS,4,4,WL(1,0,1,9))
          CALL VVV1L2P0_1(PL(0,9),W(1,5),GC_4,ZERO,ZERO,PL(0,10),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,9),4,COEFS,4,4,WL(1,0,1
     $     ,10))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,10),2,4,3,1,1,58,H)
C         Coefficient construction for loop diagram with ID 7
          CALL FFV1L1P0_3(PL(0,0),W(1,2),GC_5,ZERO,ZERO,PL(0,11),COEFS)
          CALL MG5_1_UPDATE_WL_0_0(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,11))
          CALL FFV1L3_2(PL(0,11),W(1,4),GC_5,ZERO,ZERO,PL(0,12),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,11),4,COEFS,4,4,WL(1,0,1
     $     ,12))
          CALL FFV1L1_2(PL(0,12),W(1,6),GC_5,ZERO,ZERO,PL(0,13),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,12),4,COEFS,4,4,WL(1,0,1
     $     ,13))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,13),2,4,4,1,1,59,H)
C         Coefficient construction for loop diagram with ID 8
          CALL FFV1L3_1(PL(0,0),W(1,2),GC_5,ZERO,ZERO,PL(0,14),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,14))
          CALL FFV1L2P0_3(PL(0,14),W(1,4),GC_5,ZERO,ZERO,PL(0,15)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_0(WL(1,0,1,14),4,COEFS,4,4,WL(1,0,1
     $     ,15))
          CALL VVV1L2P0_1(PL(0,15),W(1,6),GC_4,ZERO,ZERO,PL(0,16)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,15),4,COEFS,4,4,WL(1,0,1
     $     ,16))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,16),2,4,4,1,1,60,H)
C         Coefficient construction for loop diagram with ID 9
          CALL FFV1L2P0_3(PL(0,0),W(1,1),GC_5,ZERO,ZERO,PL(0,17),COEFS)
          CALL MG5_1_UPDATE_WL_0_0(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,17))
          CALL FFV1L3_1(PL(0,17),W(1,3),GC_5,ZERO,ZERO,PL(0,18),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,17),4,COEFS,4,4,WL(1,0,1
     $     ,18))
          CALL FFV1L2_1(PL(0,18),W(1,8),GC_5,ZERO,ZERO,PL(0,19),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,18),4,COEFS,4,4,WL(1,0,1
     $     ,19))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,19),2,4,5,1,1,61,H)
C         Coefficient construction for loop diagram with ID 10
          CALL FFV1L3_2(PL(0,0),W(1,1),GC_5,ZERO,ZERO,PL(0,20),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,20))
          CALL FFV1L1P0_3(PL(0,20),W(1,3),GC_5,ZERO,ZERO,PL(0,21)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_0(WL(1,0,1,20),4,COEFS,4,4,WL(1,0,1
     $     ,21))
          CALL VVV1L2P0_1(PL(0,21),W(1,8),GC_4,ZERO,ZERO,PL(0,22)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,21),4,COEFS,4,4,WL(1,0,1
     $     ,22))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,22),2,4,5,1,1,62,H)
C         Coefficient construction for loop diagram with ID 11
          CALL FFV1L1P0_3(PL(0,20),W(1,2),GC_5,ZERO,ZERO,PL(0,23)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_0(WL(1,0,1,20),4,COEFS,4,4,WL(1,0,1
     $     ,23))
          CALL FFV1L3_2(PL(0,23),W(1,4),GC_5,ZERO,ZERO,PL(0,24),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,23),4,COEFS,4,4,WL(1,0,1
     $     ,24))
          CALL FFV1L1P0_3(PL(0,24),W(1,3),GC_5,ZERO,ZERO,PL(0,25)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_2_0(WL(1,0,1,24),4,COEFS,4,4,WL(1,0,1
     $     ,25))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,25),2,4,6,1,1,63,H)
C         Coefficient construction for loop diagram with ID 12
          CALL FFV1L3_1(PL(0,23),W(1,3),GC_5,ZERO,ZERO,PL(0,26),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,23),4,COEFS,4,4,WL(1,0,1
     $     ,26))
          CALL FFV1L2P0_3(PL(0,26),W(1,4),GC_5,ZERO,ZERO,PL(0,27)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_2_0(WL(1,0,1,26),4,COEFS,4,4,WL(1,0,1
     $     ,27))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,27),2,4,7,1,1,64,H)
C         Coefficient construction for loop diagram with ID 13
          CALL VVV1L2P0_1(PL(0,23),W(1,7),GC_4,ZERO,ZERO,PL(0,28)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,23),4,COEFS,4,4,WL(1,0,1
     $     ,28))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,28),2,4,8,1,1,65,H)
C         Coefficient construction for loop diagram with ID 14
          CALL FFV1L3_1(PL(0,17),W(1,2),GC_5,ZERO,ZERO,PL(0,29),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,17),4,COEFS,4,4,WL(1,0,1
     $     ,29))
          CALL FFV1L2_1(PL(0,29),W(1,7),GC_5,ZERO,ZERO,PL(0,30),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,29),4,COEFS,4,4,WL(1,0,1
     $     ,30))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,30),2,4,8,1,1,66,H)
C         Coefficient construction for loop diagram with ID 15
          CALL FFV1L2P0_3(PL(0,29),W(1,4),GC_5,ZERO,ZERO,PL(0,31)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_0(WL(1,0,1,29),4,COEFS,4,4,WL(1,0,1
     $     ,31))
          CALL FFV1L3_1(PL(0,31),W(1,3),GC_5,ZERO,ZERO,PL(0,32),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,31),4,COEFS,4,4,WL(1,0,1
     $     ,32))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,32),2,4,6,1,1,67,H)
C         Coefficient construction for loop diagram with ID 16
          CALL FFV1L3_1(PL(0,21),W(1,2),GC_5,ZERO,ZERO,PL(0,33),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,21),4,COEFS,4,4,WL(1,0,1
     $     ,33))
          CALL FFV1L2P0_3(PL(0,33),W(1,4),GC_5,ZERO,ZERO,PL(0,34)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_2_0(WL(1,0,1,33),4,COEFS,4,4,WL(1,0,1
     $     ,34))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,34),2,4,9,1,1,68,H)
C         Coefficient construction for loop diagram with ID 17
          CALL FFV1L2_1(PL(0,0),W(1,5),GC_5,MDL_MB,ZERO,PL(0,35),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,35))
          CALL FFV1L2_1(PL(0,35),W(1,7),GC_5,MDL_MB,ZERO,PL(0,36)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,35),4,COEFS,4,4,WL(1,0,1
     $     ,36))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,36),2,4,10,1,1,69,H)
C         Coefficient construction for loop diagram with ID 18
          CALL FFV1L2_1(PL(0,0),W(1,6),GC_5,MDL_MB,ZERO,PL(0,37),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,37))
          CALL FFV1L2_1(PL(0,37),W(1,8),GC_5,MDL_MB,ZERO,PL(0,38)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,37),4,COEFS,4,4,WL(1,0,1
     $     ,38))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,38),2,4,11,1,1,70,H)
C         Coefficient construction for loop diagram with ID 19
          CALL FFV1L2_1(PL(0,0),W(1,5),GC_5,MDL_MT,MDL_WT,PL(0,39)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,39))
          CALL FFV1L2_1(PL(0,39),W(1,7),GC_5,MDL_MT,MDL_WT,PL(0,40)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,39),4,COEFS,4,4,WL(1,0,1
     $     ,40))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,40),2,4,12,1,1,71,H)
C         Coefficient construction for loop diagram with ID 20
          CALL FFV1L2_1(PL(0,0),W(1,6),GC_5,MDL_MT,MDL_WT,PL(0,41)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,41))
          CALL FFV1L2_1(PL(0,41),W(1,8),GC_5,MDL_MT,MDL_WT,PL(0,42)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,41),4,COEFS,4,4,WL(1,0,1
     $     ,42))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,42),2,4,13,1,1,72,H)
C         Coefficient construction for loop diagram with ID 21
          CALL VVV1L2P0_1(PL(0,0),W(1,5),GC_4,ZERO,ZERO,PL(0,43),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,43))
          CALL VVV1L2P0_1(PL(0,43),W(1,7),GC_4,ZERO,ZERO,PL(0,44)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,43),4,COEFS,4,4,WL(1,0,1
     $     ,44))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,44),2,4,1,2,1,73,H)
C         Coefficient construction for loop diagram with ID 22
          CALL VVV1L2P0_1(PL(0,0),W(1,6),GC_4,ZERO,ZERO,PL(0,45),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),4,COEFS,4,4,WL(1,0,1
     $     ,45))
          CALL VVV1L2P0_1(PL(0,45),W(1,8),GC_4,ZERO,ZERO,PL(0,46)
     $     ,COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,45),4,COEFS,4,4,WL(1,0,1
     $     ,46))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,46),2,4,2,2,1,74,H)
C         Coefficient construction for loop diagram with ID 23
          CALL GHGHGL2_1(PL(0,0),W(1,5),GC_4,ZERO,ZERO,PL(0,47),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),1,COEFS,1,1,WL(1,0,1
     $     ,47))
          CALL GHGHGL2_1(PL(0,47),W(1,7),GC_4,ZERO,ZERO,PL(0,48),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,47),1,COEFS,1,1,WL(1,0,1
     $     ,48))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,48),2,1,1,1,1,75,H)
C         Coefficient construction for loop diagram with ID 24
          CALL GHGHGL2_1(PL(0,0),W(1,6),GC_4,ZERO,ZERO,PL(0,49),COEFS)
          CALL MG5_1_UPDATE_WL_0_1(WL(1,0,1,0),1,COEFS,1,1,WL(1,0,1
     $     ,49))
          CALL GHGHGL2_1(PL(0,49),W(1,8),GC_4,ZERO,ZERO,PL(0,50),COEFS)
          CALL MG5_1_UPDATE_WL_1_1(WL(1,0,1,49),1,COEFS,1,1,WL(1,0,1
     $     ,50))
          CALL MG5_1_CREATE_LOOP_COEFS(WL(1,0,1,50),2,1,2,1,1,76,H)
C         At this point, all loop coefficients needed for (QCD=6),
C          i.e. of split order ID=1, are computed.
          IF(FILTER_SO.AND.SQSO_TARGET.EQ.1) GOTO 4000
 4000     CONTINUE
          LOOP_REQ_SO_DONE=.TRUE.



        ENDIF
      ENDDO

C     Grouping of loop diagrams now done directly when creating the
C      LOOPCOEFS.

      DO I=0,NSQUAREDSO
        DO J=1,3
          BUFFR_BIS(J,I)=ANS(J,I)
        ENDDO
      ENDDO


C     MadLoop jumps to this label during stability checks when it
C      recomputes the same PS point with a different CTMode
 300  CONTINUE

C     Free cache when using IREGI
      IF(IREGIRECY.AND.MLREDUCTIONLIB(I_LIB).EQ.3) THEN
        CALL IREGI_FREE_PS()
      ENDIF



      DO I=0,NSQUAREDSO
        DO J=1,3
          ANS(J,I)=BUFFR_BIS(J,I)
        ENDDO
      ENDDO


      IF(SKIPLOOPEVAL) THEN
        GOTO 1226
      ENDIF

      DO I_SO=1,NSQUAREDSO
        DO J=1,NLOOPGROUPS
          S(I_SO,J)=.TRUE.
        ENDDO
        IF (FILTER_SO.AND.SQSO_TARGET.NE.I_SO) GOTO 5001
C       CutTools call for loop numbers 1,19,21
        CALL MG5_1_LOOP_2(5,7,DCMPLX(ZERO),DCMPLX(ZERO),2,I_SO,1)
C       CutTools call for loop numbers 2,20,22
        CALL MG5_1_LOOP_2(6,8,DCMPLX(ZERO),DCMPLX(ZERO),2,I_SO,2)
C       CutTools call for loop numbers 3,4
        CALL MG5_1_LOOP_3(3,4,5,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO)
     $   ,2,I_SO,3)
C       CutTools call for loop numbers 5,6
        CALL MG5_1_LOOP_3(2,4,6,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO)
     $   ,2,I_SO,4)
C       CutTools call for loop numbers 7,8
        CALL MG5_1_LOOP_3(1,3,8,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO)
     $   ,2,I_SO,5)
C       CutTools call for loop numbers 9,13
        CALL MG5_1_LOOP_4(1,2,4,3,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO
     $   ),DCMPLX(ZERO),2,I_SO,6)
C       CutTools call for loop numbers 10
        CALL MG5_1_LOOP_4(1,2,3,4,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO
     $   ),DCMPLX(ZERO),2,I_SO,7)
C       CutTools call for loop numbers 11,12
        CALL MG5_1_LOOP_3(1,2,7,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO)
     $   ,2,I_SO,8)
C       CutTools call for loop numbers 14
        CALL MG5_1_LOOP_4(1,3,2,4,DCMPLX(ZERO),DCMPLX(ZERO),DCMPLX(ZERO
     $   ),DCMPLX(ZERO),2,I_SO,9)
C       CutTools call for loop numbers 15
        CALL MG5_1_LOOP_2(5,7,DCMPLX(MDL_MB),DCMPLX(MDL_MB),2,I_SO,10)
C       CutTools call for loop numbers 16
        CALL MG5_1_LOOP_2(6,8,DCMPLX(MDL_MB),DCMPLX(MDL_MB),2,I_SO,11)
C       CutTools call for loop numbers 17
        CALL MG5_1_LOOP_2(5,7,DCMPLX(MDL_MT),DCMPLX(MDL_MT),2,I_SO,12)
C       CutTools call for loop numbers 18
        CALL MG5_1_LOOP_2(6,8,DCMPLX(MDL_MT),DCMPLX(MDL_MT),2,I_SO,13)
C       At this point, all reductions needed for (QCD=6), i.e. of
C        split order ID=1, are computed.
        IF(FILTER_SO.AND.SQSO_TARGET.EQ.1) GOTO 5000
        GOTO 5001
 5000   CONTINUE
        CTCALL_REQ_SO_DONE=.TRUE.
 5001   CONTINUE
      ENDDO

      DO I=1,NLOOPGROUPS
        LTEMP=.TRUE.
        DO K=1,NSQUAREDSO
          IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.K) THEN
            IF (.NOT.S(K,I)) LTEMP=.FALSE.
            DO J=1,3
              ANS(J,K)=ANS(J,K)+LOOPRES(J,K,I)
              ANS(J,0)=ANS(J,0)+LOOPRES(J,K,I)
            ENDDO
          ENDIF
        ENDDO
        IF((CTMODERUN.NE.-1).AND..NOT.CHECKPHASE.AND.(.NOT.LTEMP)) THEN
          WRITE(*,*) '##W03 WARNING Contribution ',I,' is unstable.'
        ENDIF
      ENDDO

 1226 CONTINUE

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
        IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN
C         TO KEEP TRACK OF THE FINAL ANSWER TO BE RETURNED DURING
C          CHECK PHASE
          DO I=0,NSQUAREDSO
            DO K=1,3
              BUFFR(K,I)=BUFFR(K,I)+ANS(K,I)
            ENDDO
          ENDDO
        ENDIF
C       SAVE RESULT OF EACH INDEPENDENT HELICITY FOR COMPARISON DURING
C        THE HELICITY FILTER SETUP
        HELSAVED(1,HELPICKED)=0.0D0
        HELSAVED(2,HELPICKED)=0.0D0
        HELSAVED(3,HELPICKED)=0.0D0
        DO I=1,NSQUAREDSO
          IF (CHOSEN_SO_CONFIGS(I)) THEN
            HELSAVED(1,HELPICKED)=HELSAVED(1,HELPICKED)+ANS(1,I)
            HELSAVED(2,HELPICKED)=HELSAVED(2,HELPICKED)+ANS(2,I)
            HELSAVED(3,HELPICKED)=HELSAVED(3,HELPICKED)+ANS(3,I)
          ENDIF
        ENDDO

        IF (CHECKPHASE.AND.NTRY.NE.0) THEN
C         SET THE HELICITY FILTER
          IF(.NOT.FOUNDHELFILTER) THEN
            HEL_INCONSISTENT=.FALSE.
            IF(MG5_1_ISZERO(DABS(HELSAVED(1,HELPICKED))+DABS(HELSAVED(2
     $       ,HELPICKED))+DABS(HELSAVED(3,HELPICKED)),REF/DBLE(NCOMB),
     $       -1,-1)) THEN
              IF(NTRY.EQ.1) THEN
                GOODHEL(HELPICKED)=-HELOFFSET
              ELSEIF(GOODHEL(HELPICKED).NE.-HELOFFSET) THEN
                WRITE(*,*) '##W02A WARNING Inconsistent zero helicit'
     $           //'y ',HELPICKED
                IF(HELINITSTARTOVER) THEN
                  WRITE(*,*) '##I01 INFO Initialization starting ove'
     $             //'r because of inconsistency in the helicit'
     $             //'y filter setup.'
                  NTRY=0
                ELSE
                  HEL_INCONSISTENT=.TRUE.
                ENDIF
              ENDIF
            ELSEIF(HELICITYFILTERLEVEL.GT.1) THEN
              DO H=1,HELPICKED-1
                IF(GOODHEL(H).GT.-HELOFFSET) THEN
C                 Be looser for helicity check, bring a factor 100
                  DUMMY=MG5_1_ISSAME(HELSAVED(1,HELPICKED),HELSAVED(1
     $             ,H),REF,.FALSE.)
                  IF(DUMMY.NE.0) THEN
                    IF(NTRY.EQ.1) THEN
C                     Set the matching helicity to be contributing
C                      once more
                      GOODHEL(H)=GOODHEL(H)+DUMMY
C                     Use an offset to clearly show it is linked to an
C                      other one and to avoid overlap
                      GOODHEL(HELPICKED)=-H-HELOFFSET
C                     Make sure we have paired this hel config to the
C                      same one last PS point
                    ELSEIF(GOODHEL(HELPICKED).NE.(-H-HELOFFSET)) THEN
                      WRITE(*,*) '##W02B WARNING Inconsistent matchin'
     $                 //'g helicity ',HELPICKED
                      IF(HELINITSTARTOVER) THEN
                        WRITE(*,*) '##I01 INFO Initialization startin'
     $                   //'g over because of inconsistency in th'
     $                   //'e helicity filter setup.'
                        NTRY=0
                      ELSE
                        HEL_INCONSISTENT=.TRUE.
                      ENDIF
                    ENDIF
                  ENDIF
                ENDIF
              ENDDO
            ENDIF
            IF(HEL_INCONSISTENT) THEN
C             This helicity has unstable filter so we will always
C              compute it by itself.
C             We therefore also need to remove it from the multiplicati
C             ve factor of the corresponding helicity.
              IF(GOODHEL(HELPICKED).LT.-HELOFFSET) THEN
                GOODHEL(-GOODHEL(HELPICKED)-HELOFFSET)=GOODHEL(
     $           -GOODHEL(HELPICKED)-HELOFFSET)-1
              ENDIF
C             If several helicities were matched to that one, we need
C              to chose another one as reference and redirect the
C              others to this new one
C             Of course if it is one, then we do not need to do
C              anything (because with HELINITSTARTOVER=.FALSE. we only
C              support exactly identical Hels.)
              IF(GOODHEL(HELPICKED).GT.-HELOFFSET.AND.GOODHEL(HELPICKED
     $         ).NE.1) THEN
                NEWHELREF=-1
                DO H=1,NCOMB
                  IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN
                    IF (NEWHELREF.EQ.-1) THEN
                      NEWHELREF=H
                      GOODHEL(H)=GOODHEL(HELPICKED)-1
                    ELSE
                      GOODHEL(H)=-NEWHELREF-HELOFFSET
                    ENDIF
                  ENDIF
                ENDDO
              ENDIF
C             In all cases, from now on this helicity will be computed
C              independantly of the others.
C             In particular, it is the only thing to do if the
C              helicity was flagged not contributing.
              GOODHEL(HELPICKED)=1
            ENDIF
          ENDIF

C         SET THE LOOP FILTER
          IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN
            DO I=1,NLOOPGROUPS
              DO J=1,NSQUAREDSO
                IF(.NOT.MG5_1_ISZERO(ABS(LOOPRES(1,J,I))+ABS(LOOPRES(2
     $           ,J,I))+ABS(LOOPRES(3,J,I)),(REF*1.0D-4),I,J)) THEN
                  IF(NTRY.EQ.1) THEN
                    GOODAMP(J,I)=.TRUE.
                    LOOPFILTERBUFF(J,I)=.TRUE.
                  ELSEIF(.NOT.LOOPFILTERBUFF(J,I)) THEN
                    WRITE(*,*) '##W02 WARNING Inconsistent loop amp '
     $               ,I,'.'
                    IF(LOOPINITSTARTOVER) THEN
                      WRITE(*,*) '##I01 INFO Initialization startin'
     $                 //'g over because of inconsistency in the loo'
     $                 //'p filter setup.'
                      NTRY=0
                    ELSE
                      GOODAMP(J,I)=.TRUE.
                    ENDIF
                  ENDIF
                ENDIF
              ENDDO
            ENDDO
          ENDIF
        ELSEIF (.NOT.HELDOUBLECHECKED.AND.NTRY.NE.0)THEN
C         DOUBLE CHECK THE HELICITY FILTER
          IF (GOODHEL(HELPICKED).EQ.-HELOFFSET) THEN
            IF (.NOT.MG5_1_ISZERO(DABS(HELSAVED(1,HELPICKED))
     $       +DABS(HELSAVED(2,HELPICKED))+DABS(HELSAVED(2,HELPICKED))
     $       ,REF/DBLE(NCOMB),-1,-1)) THEN
              WRITE(*,*) '##W15 Helicity filter could not be successfu'
     $         //'lly double checked.'
              WRITE(*,*) '##One reason for this is that you might hav'
     $         //'e changed sensible parameters which affected wha'
     $         //'t are the zero helicity configurations.'
              WRITE(*,*) '##MadLoop will try to reset the Helicit'
     $         //'y filter with the next PS points it receives.'
              NTRY=0
              OPEN(29,FILE=HELFILTERFN,ERR=348)
 348          CONTINUE
              CLOSE(29,STATUS='delete')
            ENDIF
          ENDIF
          IF (GOODHEL(HELPICKED).LT.-HELOFFSET.AND.NTRY.NE.0) THEN
            IF(MG5_1_ISSAME(HELSAVED(1,HELPICKED),HELSAVED(1,ABS(GOODHE
     $       L(HELPICKED)+HELOFFSET)),REF,.TRUE.).EQ.0) THEN
              WRITE(*,*) '##W15 Helicity filter could not be successfu'
     $         //'lly double checked.'
              WRITE(*,*) '##One reason for this is that you might hav'
     $         //'e changed sensible parameters which affected th'
     $         //'e helicity dependance relations.'
              WRITE(*,*) '##MadLoop will try to reset the Helicit'
     $         //'y filter with the next PS points it receives.'
              NTRY=0
              OPEN(30,FILE=HELFILTERFN,ERR=349)
 349          CONTINUE
              CLOSE(30,STATUS='delete')
            ENDIF
          ENDIF
C         SET HELDOUBLECHECKED TO .TRUE. WHEN DONE
C         even if it failed we do not want to redo the check afterwards
C          if HELINITSTARTOVER=.FALSE.
          IF (HELPICKED.EQ.NCOMB.AND.(NTRY.NE.0.OR..NOT.HELINITSTARTOVE
     $     R)) THEN
            DONEHELDOUBLECHECK=.TRUE.
          ENDIF
        ENDIF

C       GOTO NEXT HELICITY OR FINISH
        IF(HELPICKED.NE.NCOMB) THEN
          HELPICKED=HELPICKED+1
          MP_DONE=.FALSE.
          GOTO 205
        ELSE
C         Useful printout
C         do I=1,NCOMB
C         write(*,*) 'HELSAVED(1,',I,')=',HELSAVED(1,I)
C         write(*,*) 'HELSAVED(2,',I,')=',HELSAVED(2,I)
C         write(*,*) 'HELSAVED(3,',I,')=',HELSAVED(3,I)
C         write(*,*) '   GOODHEL(',I,')=',GOODHEL(I)
C         ENDDO
          DO I=0,NSQUAREDSO
            DO K=1,3
              ANS(K,I)=BUFFR(K,I)
            ENDDO
          ENDDO
          IF(NTRY.EQ.0) THEN
            NATTEMPTS=NATTEMPTS+1
            IF(NATTEMPTS.EQ.MAXATTEMPTS) THEN
              WRITE(*,*) '##E01 ERROR Could not initialize the filter'
     $         //'s in ',MAXATTEMPTS,' trials'
              STOP 1
            ENDIF
          ENDIF
        ENDIF
      ENDIF

      DO K=1,3
        DO I=0,NSQUAREDSO
          ANS(K,I)=ANS(K,I)/DBLE(IDEN)
          IF (USERHEL.NE.-1) THEN
            ANS(K,I)=ANS(K,I)*HELAVGFACTOR
          ENDIF
        ENDDO
      ENDDO


      IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.EQ.
     $ -1)) THEN
        STAB_INDEX=STAB_INDEX+1
        IF(DOING_QP_EVALS.AND.MLREDUCTIONLIB(I_LIB).EQ.1) THEN
C         NOW,ONLY CUTTOOLS PROVIDES QP
          DO I=0,NSQUAREDSO
            DO K=1,3
              QP_RES(K,I,STAB_INDEX)=ANS(K,I)
            ENDDO
          ENDDO
        ELSE
          DO I=0,NSQUAREDSO
            DO K=1,3
              DP_RES(K,I,STAB_INDEX)=ANS(K,I)
            ENDDO
          ENDDO
        ENDIF

        IF(DOING_QP_EVALS.AND.MLREDUCTIONLIB(I_LIB).EQ.1) THEN
          BASIC_CT_MODE=4
        ELSE
          BASIC_CT_MODE=1
        ENDIF

C       BEGINNING OF THE DEFINITIONS OF THE DIFFERENT EVALUATION
C        METHODS

        IF(.NOT.EVAL_DONE(2)) THEN
          EVAL_DONE(2)=.TRUE.
          IF(LOOPLIBS_DIRECTEST(MLREDUCTIONLIB(I_LIB)))THEN
            CTMODE=BASIC_CT_MODE+1
            GOTO 300
          ELSE
C           If some TIR library would not support the loop direction
C            test (they all do for now), then we would just copy the
C            answer from mode 1 and carry on.
            STAB_INDEX=STAB_INDEX+1
            IF(DOING_QP_EVALS)THEN
              DO I=0,NSQUAREDSO
                DO K=1,3
                  QP_RES(K,I,STAB_INDEX)=ANS(K,I)
                ENDDO
              ENDDO
            ELSE
              DO I=0,NSQUAREDSO
                DO K=1,3
                  DP_RES(K,I,STAB_INDEX)=ANS(K,I)
                ENDDO
              ENDDO
            ENDIF
          ENDIF
        ENDIF

        CTMODE=BASIC_CT_MODE

        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
     $   .1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
          EVAL_DONE(3)=.TRUE.
          CALL MG5_1_ROTATE_PS(PS,P,1)
          IF (DOING_QP_EVALS) CALL MG5_1_MP_ROTATE_PS(MP_PS,MP_P,1)
          GOTO 200
        ENDIF

        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
     $   .2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
          EVAL_DONE(4)=.TRUE.
          CALL MG5_1_ROTATE_PS(PS,P,2)
          IF (DOING_QP_EVALS) CALL MG5_1_MP_ROTATE_PS(MP_PS,MP_P,2)
          GOTO 200
        ENDIF

        CALL MG5_1_ROTATE_PS(PS,P,0)
        IF (DOING_QP_EVALS) CALL MG5_1_MP_ROTATE_PS(MP_PS,MP_P,0)

C       END OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS

        IF(DOING_QP_EVALS.AND.MLREDUCTIONLIB(I_LIB).EQ.1) THEN
          CALL MG5_1_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS)
          DO I=0,NSQUAREDSO
            ACCURACY(I)=ACC(I)
          ENDDO
          RET_CODE_H=3
          RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.TRUE.
     $     ,.TRUE.)
          IF(MAXVAL(ACC).GE.MLSTABTHRES) THEN
            I_QP_LIB=I_QP_LIB+1
            IF(I_QP_LIB.GT.QP_NLOOPLIB.OR.INDEX_QP_TOOLS(I_QP_LIB
     $       ).EQ.0)THEN
              RET_CODE_H=4
              RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.TRUE.
     $         ,.FALSE.)
              NEPS=NEPS+1
              CALL MG5_1_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP)
              IF(NEPS.LE.10) THEN
                WRITE(*,*) '##W03 WARNING An unstable PS point was'
     $           ,       ' detected.'
                IF (NSQUAREDSO.NE.1) THEN
                  WRITE(*,*) '##Accuracies for each split orde'
     $             //'r, starting with the summed case'
                  WRITE(*,*) '##DP accuracies (for each split orde'
     $             //'r): ',(TEMP1(I),I=0,NSQUAREDSO)
                  WRITE(*,*) '##QP accuracies (for each split orde'
     $             //'r): ',(ACC(I),I=0,NSQUAREDSO)
                ELSE
                  WRITE(*,*) '##DP accuracy: ',TEMP1(1)
                  WRITE(*,*) '##QP accuracy: ',ACC(1)
                ENDIF
                DO J=0,NSQUAREDSO
                  IF (NSQUAREDSO.NE.1.OR.J.NE.0) THEN
                    IF (J.EQ.0) THEN
                      WRITE(*,*) 'Details for all split orders summe'
     $                 //'d :'
                    ELSE
                      WRITE(*,*) 'Details for split order index : ',J
                    ENDIF
                    WRITE(*,*) 'Best estimate (fin,1eps,2eps):',(ANS(I
     $               ,J),I=1,3)
                    WRITE(*,*) 'Finite double precision evaluations :'
     $               ,(DP_RES(1,J,I),I=1,N_DP_EVAL)
                    WRITE(*,*) 'Finite quad   precision evaluations :'
     $               ,(QP_RES(1,J,I),I=1,N_QP_EVAL)
                  ENDIF
                ENDDO
                WRITE(*,*) 'PS point specification :'
                WRITE(*,*) 'Renormalization scale MU_R=',MU_R
                DO I=1,NEXTERNAL
                  WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I)
     $             ,P(3,I)
                ENDDO
              ENDIF
              IF(NEPS.EQ.10) THEN
                WRITE(*,*) 'Further output of the details of thes'
     $           //'e unstable PS points will now be suppressed.'
              ENDIF
            ELSE
              I_LIB=INDEX_QP_TOOLS(I_QP_LIB)
              EVAL_DONE(1)=.TRUE.
              DO I=2,MAXSTABILITYLENGTH
                EVAL_DONE(I)=.FALSE.
              ENDDO
              STAB_INDEX=0
              IF(NROTATIONS_QP.GE.1)THEN
                GOTO 200
              ELSE
                GOTO 300
              ENDIF
            ENDIF
          ENDIF
        ELSEIF(.NOT.DOING_QP_EVALS)THEN
          CALL MG5_1_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS)
          IF(MAXVAL(ACC).GE.MLSTABTHRES) THEN
            I_LIB=I_LIB+1
            IF((I_LIB.GT.NLOOPLIB.OR.MLREDUCTIONLIB(I_LIB).EQ.0
     $       ).AND.QP_TOOLS_AVAILABLE)THEN
              I_LIB=INDEX_QP_TOOLS(1)
              I_QP_LIB=1
              DOING_QP_EVALS=.TRUE.
              EVAL_DONE(1)=.TRUE.
              DO I=2,MAXSTABILITYLENGTH
                EVAL_DONE(I)=.FALSE.
              ENDDO
              STAB_INDEX=0
              CTMODE=4
              GOTO 200
            ELSEIF(I_LIB.LE.NLOOPLIB.AND.MLREDUCTIONLIB(I_LIB).GT.0
     $       )THEN
              EVAL_DONE(1)=.TRUE.
              DO I=2,MAXSTABILITYLENGTH
                EVAL_DONE(I)=.FALSE.
              ENDDO
              STAB_INDEX=0
              IF(NROTATIONS_DP.GE.1)THEN
                GOTO 200
              ELSE
                GOTO 300
              ENDIF
            ELSE
              DO I=0,NSQUAREDSO
                ACCURACY(I)=ACC(I)
              ENDDO
              RET_CODE_H=4
              RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $         ,.FALSE.)
              NEPS=NEPS+1
              IF(NEPS.LE.10) THEN
                WRITE(*,*) '##W03 WARNING An unstable PS point was'
     $           ,       ' detected.'
                WRITE(*,*) '##W03 WARNING No quadruple precision wil'
     $           //'l be used.'
                IF (NSQUAREDSO.NE.1) THEN
                  WRITE(*,*) 'Accuracies for each split orde'
     $             //'r, starting with the summed case'
                  WRITE(*,*) 'DP accuracies (for each split order): '
     $             ,(ACC(I),I=0,NSQUAREDSO)
                ELSE
                  WRITE(*,*) 'DP accuracy: ',ACC(1)
                ENDIF
                DO J=0,NSQUAREDSO
                  IF (NSQUAREDSO.NE.1.OR.J.NE.0) THEN
                    IF (J.EQ.0) THEN
                      WRITE(*,*) 'Details for all split orders summe'
     $                 //'d :'
                    ELSE
                      WRITE(*,*) 'Details for split order index : ',J
                    ENDIF
                    WRITE(*,*) 'Best estimate (fin,1eps,2eps):',(ANS(I
     $               ,J),I=1,3)
                    WRITE(*,*) 'Finite double precision evaluations :'
     $               ,(DP_RES(1,J,I),I=1,N_DP_EVAL)
                  ENDIF
                ENDDO
                WRITE(*,*) 'PS point specification :'
                WRITE(*,*) 'Renormalization scale MU_R=',MU_R
                DO I=1,NEXTERNAL
                  WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I)
     $             ,P(3,I)
                ENDDO
              ENDIF
              IF(NEPS.EQ.10) THEN
                WRITE(*,*) 'Further output of the details of thes'
     $           //'e unstable PS points will now be suppressed.'
              ENDIF
            ENDIF
          ELSE
            DO I=0,NSQUAREDSO
              ACCURACY(I)=ACC(I)
            ENDDO
            RET_CODE_H=2
            RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $       ,.TRUE.)
          ENDIF
        ENDIF
      ELSE
        RET_CODE_H=1
        DO I=0,NSQUAREDSO
          ACCURACY(I)=-1.0D0
        ENDDO
        RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $   ,.FALSE.)
      ENDIF

 9999 CONTINUE

C     Finalize the return code
      IF (MP_DONE_ONCE) THEN
        RET_CODE_T=2
      ELSE
        RET_CODE_T=1
      ENDIF
      IF(CHECKPHASE.OR..NOT.HELDOUBLECHECKED) THEN
        RET_CODE_H=1
        RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $   ,.FALSE.)
        RET_CODE_T=RET_CODE_T+2
        DO I=0,NSQUAREDSO
          ACCURACY(I)=-1.0D0
        ENDDO
      ENDIF

C     Finally for the summed result in ANS(1:3,0), make sure to only
C     consider the squared order asked for by the user.
C     Notice that this filtering using CHOSEN_SO_CONFIGS happens
C     here only while everywhere else one always considers the sum.
      DO J=1,3
        ANS(J,0)=0.0D0
      ENDDO
      DO I=1,NSQUAREDSO
        IF (CHOSEN_SO_CONFIGS(I)) THEN
          DO J=1,3
            ANS(J,0)=ANS(J,0)+ANS(J,I)
          ENDDO
        ENDIF
      ENDDO

C     Reinitialize the default threshold if it was specified by the
C      user
      IF (USER_STAB_PREC.GT.0.0D0) THEN
        MLSTABTHRES=MLSTABTHRES_BU
        CTMODEINIT=CTMODEINIT_BU
      ENDIF

C     Reinitialize the check phase logicals and the filters if check
C      bypassed
      IF (BYPASS_CHECK) THEN
        CHECKPHASE = OLD_CHECKPHASE
        HELDOUBLECHECKED = OLD_HELDOUBLECHECKED
        DO I=1,NCOMB
          GOODHEL(I)=OLD_GOODHEL(I)
        ENDDO
        DO I=1,NSQUAREDSO
          DO J=1,NLOOPGROUPS
            GOODAMP(I,J)=OLD_GOODAMP(I,J)
          ENDDO
        ENDDO
      ENDIF
      END

C     --=========================================--
C     General Helper functions and subroutine
C     for the main sloopmatrix subroutine
C     --=========================================--

      LOGICAL FUNCTION MG5_1_ISZERO(TOTEST, REFERENCE_VALUE, LOOP
     $ , SOINDEX)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=13)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)
C     
C     ARGUMENTS 
C     
      REAL*8 TOTEST, REFERENCE_VALUE
      INTEGER LOOP, SOINDEX
C     
C     GLOBAL 
C     
      INCLUDE 'MadLoopParams.inc'
      COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
      LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
      COMMON/MG5_1_LOOPRES/LOOPRES,S
C     ----------
C     BEGIN CODE
C     ----------
      IF(ABS(REFERENCE_VALUE).EQ.0.0D0) THEN
        MG5_1_ISZERO=.FALSE.
        WRITE(*,*) '##E02 ERRROR Reference value for comparison i'
     $   //'s zero.'
        STOP 1
      ELSE
        MG5_1_ISZERO=((ABS(TOTEST)/ABS(REFERENCE_VALUE)).LT.ZEROTHRES)
      ENDIF

      IF(LOOP.NE.-1) THEN
        IF((.NOT.MG5_1_ISZERO).AND.(.NOT.S(SOINDEX,LOOP))) THEN
          WRITE(*,*) '##W01 WARNING Contribution ',LOOP,' of spli'
     $     //'t order ',SOINDEX,' is detected as contributing wit'
     $     //'h CR=',(ABS(TOTEST)/ABS(REFERENCE_VALUE)),' but i'
     $     //'s unstable.'
        ENDIF
      ENDIF

      END

      INTEGER FUNCTION MG5_1_ISSAME(RESA,RESB,REF,USEMAX)
      IMPLICIT NONE
C     This function compares the result from two different helicity
C      configuration A and B
C     It returns 0 if they are not related and (+/-wgt) if A=(+/-wgt)*B
C     .
C     For now, the only wgt implemented is the integer 1 or -1.
C     If useMax is .TRUE., it uses all implemented weights no matter
C      what is HELINITSTARTOVER
C     
C     CONSTANTS
C     
      INTEGER MAX_WGT_TO_TRY
      PARAMETER (MAX_WGT_TO_TRY=2)
C     
C     ARGUMENTS 
C     
      REAL*8 RESA(3), RESB(3)
      REAL*8 REF
      LOGICAL USEMAX
C     
C     LOCAL VARIABLES
C     
      LOGICAL MG5_1_ISZERO
      INTEGER I,J
      INTEGER N_WGT_TO_TRY
      INTEGER WGT_TO_TRY(MAX_WGT_TO_TRY)
      DATA WGT_TO_TRY/1,-1/
C     
C     INCLUDES
C     
      INCLUDE 'MadLoopParams.inc'
C     ----------
C     BEGIN CODE
C     ----------
      MG5_1_ISSAME=0

C     If the helicity can be constructed progressively while allowing
C      inconsistency, then we only allow for weight one comparisons.
      IF (.NOT.HELINITSTARTOVER.AND..NOT.USEMAX) THEN
        N_WGT_TO_TRY=1
      ELSE
        N_WGT_TO_TRY=MAX_WGT_TO_TRY
      ENDIF

      DO I=1,N_WGT_TO_TRY
        DO J=1,3
          IF (MG5_1_ISZERO(ABS(RESB(J)),REF,-1,-1)) THEN
            IF(.NOT.MG5_1_ISZERO(ABS(RESB(J))+ABS(RESA(J)),REF,-1,
     $       -1)) THEN
              GOTO 1231
            ENDIF
C           Be looser for helicity comparison, so bring a factor 100
          ELSEIF(.NOT.MG5_1_ISZERO(ABS((RESA(J)/RESB(J))-DBLE(WGT_TO_TR
     $     Y(I))),1.0D0,-1,-1)) THEN
            GOTO 1231
          ENDIF
        ENDDO
        MG5_1_ISSAME = WGT_TO_TRY(I)
        RETURN
 1231   CONTINUE
      ENDDO
      END

      SUBROUTINE MG5_1_COMPUTE_ACCURACY(FULLLIST, LENGTH, ACC
     $ , ESTIMATE)
      IMPLICIT NONE
C     
C     PARAMETERS 
C     
      INTEGER MAXSTABILITYLENGTH
      COMMON/MG5_1_STABILITY_TESTS/MAXSTABILITYLENGTH
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)
C     
C     ARGUMENTS 
C     
      REAL*8 FULLLIST(3,0:NSQUAREDSO,MAXSTABILITYLENGTH)
      INTEGER LENGTH
      REAL*8 ACC(0:NSQUAREDSO), ESTIMATE(0:3,0:NSQUAREDSO)
C     
C     LOCAL VARIABLES 
C     
      LOGICAL MASK(MAXSTABILITYLENGTH)
      LOGICAL MASK3(3)
      DATA MASK3/.TRUE.,.TRUE.,.TRUE./
      INTEGER I,J,K
      REAL*8 AVG
      REAL*8 DIFF
      REAL*8 ACCURACIES(3)
      REAL*8 LIST(MAXSTABILITYLENGTH)

C     ----------
C     BEGIN CODE
C     ----------
      DO I=1,LENGTH
        MASK(I)=.TRUE.
      ENDDO
      DO I=LENGTH+1,MAXSTABILITYLENGTH
        MASK(I)=.FALSE.
C       For some architectures, it is necessary to initialize all the
C        elements of fulllist(i,j)
C       Beware that if the length provided is incorrect, then this can
C        corrup the fulllist given in argument.
        DO J=0,NSQUAREDSO
          DO K=1,3
            FULLLIST(K,J,I)=0.0D0
          ENDDO
        ENDDO
      ENDDO

      DO K=0,NSQUAREDSO

        DO I=1,3
          DO J=1,MAXSTABILITYLENGTH
            LIST(J)=FULLLIST(I,K,J)
          ENDDO
          DIFF=MAXVAL(LIST,1,MASK)-MINVAL(LIST,1,MASK)
          AVG=(MAXVAL(LIST,1,MASK)+MINVAL(LIST,1,MASK))/2.0D0
          ESTIMATE(I,K)=AVG
          IF (AVG.EQ.0.0D0) THEN
            ACCURACIES(I)=DIFF
          ELSE
            ACCURACIES(I)=DIFF/ABS(AVG)
          ENDIF
        ENDDO

C       The technique below is too sensitive, typically to
C       unstablities in very small poles
C       acc(k)=MAXVAL(ACCURACIES,1,MASK3)
C       The following is used instead
        ACC(K) = 0.0D0
        AVG = 0.0D0
        DO I=1,3
          ACC(K) = ACC(K) + ACCURACIES(I)*ABS(ESTIMATE(I,K))
          AVG = AVG + ESTIMATE(I,K)
        ENDDO
        IF (AVG.NE.0.0D0) THEN
          ACC(K)  = ACC(K) / ( ABS(AVG) / 3.0D0)
        ENDIF

      ENDDO

      END

      SUBROUTINE MG5_1_SET_N_EVALS(N_DP_EVALS,N_QP_EVALS)

      IMPLICIT NONE
      INTEGER N_DP_EVALS, N_QP_EVALS

      INCLUDE 'MadLoopParams.inc'

      IF(CTMODERUN.LE.-1) THEN
        N_DP_EVALS=2+NROTATIONS_DP
        N_QP_EVALS=2+NROTATIONS_QP
      ELSE
        N_DP_EVALS=1
        N_QP_EVALS=1
      ENDIF

      IF(N_DP_EVALS.GT.20.OR.N_QP_EVALS.GT.20) THEN
        WRITE(*,*) 'ERROR:: Increase hardcoded maxstabilitylength.'
        STOP 1
      ENDIF

      END

C     THIS SUBROUTINE SIMPLY SET THE GLOBAL PS CONFIGURATION GLOBAL
C      VARIABLES FROM A GIVEN VARIABLE IN DOUBLE PRECISION
      SUBROUTINE MG5_1_SET_MP_PS(P)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
      COMMON/MG5_1_MP_PSPOINT/MP_PS,MP_P
      REAL*8 P(0:3,NEXTERNAL)

      DO I=1,NEXTERNAL
        DO J=0,3
          MP_PS(J,I)=P(J,I)
        ENDDO
      ENDDO
      CALL MG5_1_MP_IMPROVE_PS_POINT_PRECISION(MP_PS)
      DO I=1,NEXTERNAL
        DO J=0,3
          MP_P(J,I)=MP_PS(J,I)
        ENDDO
      ENDDO

      END

C     --=========================================--
C     Functions for dealing with the ordering
C     and indexing of split order contributions
C     --=========================================--

      SUBROUTINE MG5_1_GET_NSQSO_LOOP(NSQSO)
C     
C     Simple subroutine returning the number of squared split order
C     contributions returned in ANS when calling sloopmatrix 
C     
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)

      INTEGER NSQSO

      NSQSO=NSQUAREDSO

      END

      SUBROUTINE MG5_1_GET_ANSWER_DIMENSION(ANS_DIM)
C     
C     MadLoop subroutines return an array of dimension ANS(0:3,0:ANS_DI
C     M)
C     In order for the user program to be able to correctly declare
C      this
C     array when calling MadLoop, this subroutine returns its dimension
C     
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)
      INTEGER ANS_DIM

      INCLUDE 'nsqso_born.inc'

      ANS_DIM=MAX(NSQSO_BORN,NSQUAREDSO)

      END

      INTEGER FUNCTION MG5_1_ML5SOINDEX_FOR_SQUARED_ORDERS(ORDERS)
C     
C     This functions returns the integer index identifying the split
C      orders list passed in argument which correspond to the values
C      of the following list of couplings (and in this order):
C     ['QCD']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQSO
      PARAMETER (NSO=1, NSQSO=1)
C     
C     ARGUMENTS
C     
      INTEGER ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I,J
      INTEGER SQPLITORDERS(NSQSO,NSO)
      DATA (SQPLITORDERS(  1,I),I=  1,  1) /    6/
      COMMON/MG5_1_ML5SQPLITORDERS/SQPLITORDERS
C     
C     BEGIN CODE
C     
      DO I=1,NSQSO
        DO J=1,NSO
          IF (ORDERS(J).NE.SQPLITORDERS(I,J)) GOTO 1009
        ENDDO
        MG5_1_ML5SOINDEX_FOR_SQUARED_ORDERS = I
        RETURN
 1009   CONTINUE
      ENDDO

      WRITE(*,*) 'ERROR:: Stopping function MG5_1_ML5SOINDEX_FOR_SQUAR'
     $ //'ED_ORDERS'
      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
      STOP

      END

      INTEGER FUNCTION MG5_1_ML5SOINDEX_FOR_BORN_AMP(AMPID)
C     
C     For a given born amplitude number, it returns the ID of the
C      split orders it has
C     
C     CONSTANTS
C     
      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=2)
C     
C     ARGUMENTS
C     
      INTEGER AMPID
C     
C     LOCAL VARIABLES
C     
      INTEGER BORNAMPORDERS(NBORNAMPS)
      DATA (BORNAMPORDERS(I),I=  1,  2) /    1,    1/
C     -----------
C     BEGIN CODE
C     -----------
      IF (AMPID.GT.NBORNAMPS) THEN
        WRITE(*,*) 'ERROR:: Born amplitude ID ',AMPID,' above th'
     $   //'e maximum ',NBORNAMPS
      ENDIF
      MG5_1_ML5SOINDEX_FOR_BORN_AMP = BORNAMPORDERS(AMPID)

      END

      INTEGER FUNCTION MG5_1_ML5SOINDEX_FOR_LOOP_AMP(AMPID)
C     
C     For a given loop amplitude number, it returns the ID of the
C      split orders it has
C     
C     CONSTANTS
C     
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=76)
C     
C     ARGUMENTS
C     
      INTEGER AMPID
C     
C     LOCAL VARIABLES
C     
      INTEGER LOOPAMPORDERS(NLOOPAMPS)
      DATA (LOOPAMPORDERS(I),I=  1,  5) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I=  6, 10) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 11, 15) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 16, 20) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 21, 25) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 26, 30) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 31, 35) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 36, 40) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 41, 45) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 46, 50) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 51, 55) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 56, 60) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 61, 65) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 66, 70) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 71, 75) /    2,    2,    2,    2,    2/
      DATA (LOOPAMPORDERS(I),I= 76, 76) /    2/
C     -----------
C     BEGIN CODE
C     -----------
      IF (AMPID.GT.NLOOPAMPS) THEN
        WRITE(*,*) 'ERROR:: Loop amplitude ID ',AMPID,' above th'
     $   //'e maximum ',NLOOPAMPS
      ENDIF
      MG5_1_ML5SOINDEX_FOR_LOOP_AMP = LOOPAMPORDERS(AMPID)

      END


      INTEGER FUNCTION MG5_1_ML5SQSOINDEX(ORDERINDEXA, ORDERINDEXB)
C     
C     This functions plays the role of the interference matrix. It can
C      be hardcoded or 
C     made more elegant using hashtables if its execution speed ever
C      becomes a relevant
C     factor. From two split order indices, it return the corresponding
C      index in the squared 
C     order canonical ordering.
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=2)
C     
C     ARGUMENTS
C     
      INTEGER ORDERINDEXA, ORDERINDEXB
C     
C     LOCAL VARIABLES
C     
      INTEGER I, SQORDERS(NSO)
      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
      DATA (AMPSPLITORDERS(  1,I),I=  1,  1) /    2/
      DATA (AMPSPLITORDERS(  2,I),I=  1,  1) /    4/
      COMMON/MG5_1_ML5AMPSPLITORDERS/AMPSPLITORDERS
C     
C     FUNCTION
C     
      INTEGER MG5_1_ML5SOINDEX_FOR_SQUARED_ORDERS
C     
C     BEGIN CODE
C     
      DO I=1,NSO
        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
     $   NDEXB,I)
      ENDDO
      MG5_1_ML5SQSOINDEX=MG5_1_ML5SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
      END

C     This is the inverse subroutine of ML5SOINDEX_FOR_SQUARED_ORDERS.
C      Not directly useful, but provided nonetheless.
      SUBROUTINE MG5_1_ML5GET_SQUARED_ORDERS_FOR_SOINDEX(SOINDEX
     $ ,ORDERS)
C     
C     This functions returns the orders identified by the squared
C      split order index in argument. Order values correspond to
C      following list of couplings (and in this order):
C     ['QCD']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQSO
      PARAMETER (NSO=1, NSQSO=1)
C     
C     ARGUMENTS
C     
      INTEGER SOINDEX, ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I
      INTEGER SQPLITORDERS(NSQSO,NSO)
      COMMON/MG5_1_ML5SQPLITORDERS/SQPLITORDERS
C     
C     BEGIN CODE
C     
      IF (SOINDEX.GT.0.AND.SOINDEX.LE.NSQSO) THEN
        DO I=1,NSO
          ORDERS(I) =  SQPLITORDERS(SOINDEX,I)
        ENDDO
        RETURN
      ENDIF

      WRITE(*,*) 'ERROR:: Stopping function MG5_1_ML5GET_SQUARED_ORDER'
     $ //'S_FOR_SOINDEX'
      WRITE(*,*) 'Could not find squared orders index ',SOINDEX
      STOP

      END SUBROUTINE

C     This is the inverse subroutine of getting amplitude SO orders.
C      Not directly useful, but provided nonetheless.
      SUBROUTINE MG5_1_ML5GET_ORDERS_FOR_AMPSOINDEX(SOINDEX,ORDERS)
C     
C     This functions returns the orders identified by the split order
C      index in argument. Order values correspond to following list of
C      couplings (and in this order):
C     ['QCD']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NAMPSO
      PARAMETER (NSO=1, NAMPSO=2)
C     
C     ARGUMENTS
C     
      INTEGER SOINDEX, ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I
      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
      COMMON/MG5_1_ML5AMPSPLITORDERS/AMPSPLITORDERS
C     
C     BEGIN CODE
C     
      IF (SOINDEX.GT.0.AND.SOINDEX.LE.NAMPSO) THEN
        DO I=1,NSO
          ORDERS(I) =  AMPSPLITORDERS(SOINDEX,I)
        ENDDO
        RETURN
      ENDIF

      WRITE(*,*) 'ERROR:: Stopping function MG5_1_ML5GET_ORDERS_FOR_AM'
     $ //'PSOINDEX'
      WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX
      STOP

      END SUBROUTINE


C     This function is not directly useful, but included for completene
C     ss
      INTEGER FUNCTION MG5_1_ML5SOINDEX_FOR_AMPORDERS(ORDERS)
C     
C     This functions returns the integer index identifying the
C      amplitude split orders passed in argument which correspond to
C      the values of the following list of couplings (and in this
C      order):
C     ['QCD']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NAMPSO
      PARAMETER (NSO=1, NAMPSO=2)
C     
C     ARGUMENTS
C     
      INTEGER ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I,J
      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
      COMMON/MG5_1_ML5AMPSPLITORDERS/AMPSPLITORDERS
C     
C     BEGIN CODE
C     
      DO I=1,NAMPSO
        DO J=1,NSO
          IF (ORDERS(J).NE.AMPSPLITORDERS(I,J)) GOTO 1009
        ENDDO
        MG5_1_ML5SOINDEX_FOR_AMPORDERS = I
        RETURN
 1009   CONTINUE
      ENDDO

      WRITE(*,*) 'ERROR:: Stopping function MG5_1_ML5SOINDEX_FOR_AMPOR'
     $ //'DERS'
      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
      STOP

      END

C     --=========================================--
C     Definition of additional access routines
C     --=========================================--

      SUBROUTINE MG5_1_FORCE_STABILITY_CHECK(ONOFF)
C     
C     This function can be called by the MadLoop user so as to always
C      have stability
C     checked, even during initialisation, when calling the *_thres
C      routines.
C     
      LOGICAL ONOFF

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      DATA BYPASS_CHECK, ALWAYS_TEST_STABILITY /.FALSE.,.FALSE./
      COMMON/MG5_1_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

      ALWAYS_TEST_STABILITY = ONOFF

      END SUBROUTINE

      SUBROUTINE MG5_1_SET_AUTOMATIC_TIR_CACHE_CLEARING(ONOFF)
C     
C     This function can be called by the MadLoop user so as to
C      manually chose when
C     to reset the TIR cache.
C     
      IMPLICIT NONE

      INCLUDE 'MadLoopParams.inc'

      LOGICAL ONOFF

      LOGICAL AUTOMATIC_TIR_CACHE_CLEARING
      DATA AUTOMATIC_TIR_CACHE_CLEARING/.TRUE./
      COMMON/MG5_1_RUNTIME_OPTIONS/AUTOMATIC_TIR_CACHE_CLEARING

      INTEGER N_DP_EVAL, N_QP_EVAL
      COMMON/MG5_1_N_EVALS/N_DP_EVAL,N_QP_EVAL

      WRITE(*,*) 'Warning: No TIR caching implemented. Call t'
     $ //'o SET_AUTOMATIC_TIR_CACHE_CLEARING did nothing.'
      END SUBROUTINE

      SUBROUTINE MG5_1_SET_COUPLINGORDERS_TARGET(SOTARGET)
      IMPLICIT NONE
C     
C     This routine can be accessed by an external user to set the
C      squared split order target.
C     If set to a value different than -1, the code will try to avoid
C      computing anything which
C     does not contribute to contributions of squared split orders
C      SQSO_TARGET and below.
C     This can considerably speed up the code. However, keep in mind
C      that any contribution of
C     'squared order index' larger than SQSO_TARGET cannot be trust.
C     
C     ARGUMENTS
C     
      INTEGER SOTARGET
C     
C     GLOBAL
C     
      INTEGER SQSO_TARGET
      COMMON/MG5_1_SOCHOICE/SQSO_TARGET
C     ----------
C     BEGIN CODE
C     ----------
      SQSO_TARGET = SOTARGET
      END

      SUBROUTINE MG5_1_SLOOPMATRIXHEL(P,HEL,ANS)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)
C     
C     ARGUMENTS 
C     
      REAL*8 P(0:3,NEXTERNAL)
      REAL*8 ANS(0:3,0:NSQUAREDSO)
      INTEGER HEL, USERHEL
      COMMON/MG5_1_USERCHOICE/USERHEL
C     ----------
C     BEGIN CODE
C     ----------
      USERHEL=HEL
      CALL MG5_1_SLOOPMATRIX(P,ANS)
      END

      SUBROUTINE MG5_1_SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED
     $ ,PREC_FOUND,RET_CODE)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)
C     
C     ARGUMENTS 
C     
      REAL*8 P(0:3,NEXTERNAL)
      REAL*8 ANS(0:3,0:NSQUAREDSO)
      INTEGER HEL, RET_CODE
      REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I
C     
C     GLOBAL VARIABLES
C     
      REAL*8 USER_STAB_PREC
      COMMON/MG5_1_USER_STAB_PREC/USER_STAB_PREC

      INTEGER H,T,U
      REAL*8 ACCURACY(0:NSQUAREDSO)
      COMMON/MG5_1_ACC/ACCURACY,H,T,U

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/MG5_1_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

C     ----------
C     BEGIN CODE
C     ----------
      USER_STAB_PREC = PREC_ASKED

      CALL MG5_1_SLOOPMATRIXHEL(P,HEL,ANS)
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0)
     $ ) THEN
        BYPASS_CHECK = .TRUE.
        CALL MG5_1_SLOOPMATRIXHEL(P,HEL,ANS)
        BYPASS_CHECK = .FALSE.
C       Make sure we correctly return an initialization-type T code
        IF (T.EQ.2) T=4
        IF (T.EQ.1) T=3
      ENDIF

C     Reset it to default value not to affect next runs
      USER_STAB_PREC = -1.0D0

      DO I=0,NSQUAREDSO
        PREC_FOUND(I)=ACCURACY(I)
      ENDDO
      RET_CODE=100*H+10*T+U

      END

      SUBROUTINE MG5_1_SLOOPMATRIX_THRES(P,ANS,PREC_ASKED,PREC_FOUND
     $ ,RET_CODE)
C     
C     Inputs are:
C     P(0:3, Nexternal)  double  :: Kinematic configuration (E,px,py,pz
C     )
C     PEC_ASKED          double  :: Target relative accuracy, -1 for
C      default
C     
C     Outputs are:
C     ANS(3)             double  :: Result (finite, single pole,
C      double pole) 
C     PREC_FOUND         double  :: Relative accuracy estimated for
C      the result
C     Returns -1 if no stab test could be performed.
C     RET_CODE			 integer :: Return code. See below for details
C     
C     Return code conventions: RET_CODE = H*100 + T*10 + U
C     
C     H == 1
C     Stability unknown.
C     H == 2
C     Stable PS (SPS) point.
C     No stability rescue was necessary.
C     H == 3
C     Unstable PS (UPS) point.
C     Stability rescue necessary, and successful.
C     H == 4
C     Exceptional PS (EPS) point.
C     Stability rescue attempted, but unsuccessful.
C     
C     T == 1
C     Default computation (double prec.) was performed.
C     T == 2
C     Quadruple precision was used for this PS point.
C     T == 3
C     MadLoop in initialization phase. Only double precision used.
C     T == 4
C     MadLoop in initialization phase. Quadruple precision used.
C     
C     U == 0
C     Not stable.
C     U == 1
C     Stable with CutTools in double precision.
C     U == 2
C     Stable with PJFry++.
C     U == 3
C     Stable with IREGI.
C     U == 4
C     Stable with Golem95
C     U == 9
C     Stable with CutTools in quadruple precision.         
C     
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=1)
C     
C     ARGUMENTS 
C     
      REAL*8 P(0:3,NEXTERNAL)
      REAL*8 ANS(0:3,0:NSQUAREDSO)
      REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
      INTEGER RET_CODE
C     
C     LOCAL VARIABLES
C     
      INTEGER I
C     
C     GLOBAL VARIABLES
C     
      REAL*8 USER_STAB_PREC
      COMMON/MG5_1_USER_STAB_PREC/USER_STAB_PREC

      INTEGER H,T,U
      REAL*8 ACCURACY(0:NSQUAREDSO)
      COMMON/MG5_1_ACC/ACCURACY,H,T,U

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/MG5_1_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

C     ----------
C     BEGIN CODE
C     ----------
      USER_STAB_PREC = PREC_ASKED
      CALL MG5_1_SLOOPMATRIX(P,ANS)
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0)
     $ ) THEN
        BYPASS_CHECK = .TRUE.
        CALL MG5_1_SLOOPMATRIX(P,ANS)
        BYPASS_CHECK = .FALSE.
C       Make sure we correctly return an initialization-type T code
        IF (T.EQ.2) T=4
        IF (T.EQ.1) T=3
      ENDIF

C     Reset it to default value not to affect next runs
      USER_STAB_PREC = -1.0D0
      DO I=0,NSQUAREDSO
        PREC_FOUND(I)=ACCURACY(I)
      ENDDO
      RET_CODE=100*H+10*T+U

      END

C     The subroutine below perform clean-up duties for MadLoop like
C      de-allocating
C     arrays
      SUBROUTINE MG5_1_EXIT_MADLOOP()
      CONTINUE
      END


