text:
C----------------------------------------------------------------------
      SUBROUTINE RCLOS()
C     DUMMY IF HBOOK IS USED
C----------------------------------------------------------------------
      END

C----------------------------------------------------------------------
      SUBROUTINE HWABEG
C     USER'S ROUTINE FOR INITIALIZATION
C----------------------------------------------------------------------
      INCLUDE 'HERWIG65.INC'
      include 'reweight0.inc'
      integer nwgt,max_weight,nwgt_analysis
      common/cnwgt/nwgt
      common/c_analysis/nwgt_analysis
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
      character*15 weights_info(max_weight)
      common/cwgtsinfo/weights_info
      REAL*8 pi
c      parameter (pi=3.14160E0)
      parameter (pi=3.141592653589793d0)
      integer ivlep1,ivlep2
      common/vvlin/ivlep1,ivlep2
      character*5 cc(3)
      data cc/'     ',' cuts',' CUTS'/
      integer maxset,maxplot,STEP,Nscales,Npdfpairs,
     # j2,k2,isc,jsc,kk1,i2,nps,nng,kk
      common/cnumplots/maxset,maxplot,STEP,Nscales,Npdfpairs

c different setups (no cuts, cuts, cuts+isolation)
c kinematic distributions we look at
      maxset=3
      maxplot=5+5+1
      Nscales=3
      Npdfpairs=20
      STEP=maxplot*(Nscales*Nscales+2*Npdfpairs)
c
      if((nwgt-1).ne.(Nscales*Nscales+2*Npdfpairs))then
        write(*,*)'Error #2',nwgt,Nscales,Npdfpairs
        CALL HWWARN('HWABEG',501)
      endif
      nwgt_analysis=nwgt
c
      call inihist

c loop over each setup
      do j2=1,maxset
      k2=(j2-1)*STEP
c loop over scales and PDFs; skip the first which is equal to (1,1)
       do kk=2,nwgt_analysis
         kk1=k2+maxplot*(kk-2)
         call mbook(kk1+ 1 ,'l+ pt      '//cc(j2)//weights_info(kk),
     #              4.d0,0.d0,400.d0)
         call mbook(kk1+ 2 ,'ll pt      '//cc(j2)//weights_info(kk),
     #              4.d0,0.d0,400.d0)
         call mbook(kk1+ 3 ,'ll m       '//cc(j2)//weights_info(kk),
     #              4.d0,0.d0,400.d0)
         call mbook(kk1+ 4 ,'ll EE      '//cc(j2)//weights_info(kk),
     #              12.d0,0.d0,1200.d0)
         call mbook(kk1+ 5 ,'pTl1+pTl2  '//cc(j2)//weights_info(kk),
     #              10.d0,0.d0,1000.d0)
c The centers of the bins correspond to the Mellin moments (eg, the bin
c centered at 1 is the first Mellin moment). The bin at zero should be
c equal to one (ie, zeroth Mellin moment), as a way of cross check
         call mbook(kk1+ 6 ,'l+ pt      '//cc(j2)//weights_info(kk),
     #              1.d0,-0.5d0,4.5d0)
         call mbook(kk1+ 7 ,'ll pt      '//cc(j2)//weights_info(kk),
     #              1.d0,-0.5d0,4.5d0)
         call mbook(kk1+ 8 ,'ll m       '//cc(j2)//weights_info(kk),
     #              1.d0,-0.5d0,4.5d0)
         call mbook(kk1+ 9 ,'ll EE      '//cc(j2)//weights_info(kk),
     #              1.d0,-0.5d0,4.5d0)
         call mbook(kk1+ 10,'pTl1+pTl2  '//cc(j2)//weights_info(kk),
     #              1.d0,-0.5d0,4.5d0)
         call mbook(kk1+ 11,'total      '//cc(j2)//weights_info(kk),
     #              1.d0,-0.5d0,4.5d0)
       enddo
      enddo
cccccccccccccccccccccccccccccccccccccccc
c Here is how scales are being scanned:
c       yfactR(1)=1.d0
c       yfactR(2)=rw_Rscale_down
c       yfactR(3)=rw_Rscale_up
c       yfactF(1)=1.d0
c       yfactF(2)=rw_Fscale_down
c       yfactF(3)=rw_Fscale_up
ccccccccccccccccccccccccccccccccccccccc
      return
      END

C----------------------------------------------------------------------
      SUBROUTINE HWAEND
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
C----------------------------------------------------------------------
      INCLUDE 'HERWIG65.INC'
      REAL*8 SONE
      PARAMETER (SONE=1.d0)
      REAL*8 XNORM
      INTEGER I,IVLEP1,IVLEP2
      COMMON/VVLIN/IVLEP1,IVLEP2
      integer maxset,maxplot,STEP,Nscales,Npdfpairs,
     # j2,k2,isc,jsc,kk1,i2,nps,nng
      common/cnumplots/maxset,maxplot,STEP,Nscales,Npdfpairs
      INTEGER KK,nwgt_analysis
      common/c_analysis/nwgt_analysis

c$$$      OPEN(UNIT=99,NAME='HERQQ.TOP',STATUS='UNKNOWN')
      OPEN(UNIT=99,FILE='HERQQ.TOP',STATUS='UNKNOWN')

C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
      XNORM=1.D3/DFLOAT(NEVHEP)
      DO I=1,5000
 	CALL MFINAL3(I)
        CALL MCOPY(I,I+5000)
        CALL MOPERA(I+5000,'F',I+5000,I+5000,XNORM,0.d0)
      ENDDO
c Normalize Mellin moments (i.e. divide plots '6'-'10' by the total x-section plot '11')
c NOTE: below '5' == (maxplot-1)/2 ; '6'='5+1' !!!
      do j2=1,maxset
       k2=(j2-1)*STEP
c loop over scales and PDFs; skip the first which is equal to (1,1)
       do kk=2,nwgt_analysis
         kk1=k2+maxplot*(kk-2)+5
         do i=1,5
           CALL MOPERA(kk1+i+5000,'/',kk1+6+5000,kk1+i+5000,SONE,SONE)
         enddo
       enddo
      enddo
c
      DO I=1,5000
 	CALL MFINAL3(I+5000)
      ENDDO
c
c loop over each setup
      do j2=1,maxset
      k2=(j2-1)*STEP
c loop over scales and PDFs; skip the first which is equal to (1,1)
       do kk=2,nwgt_analysis
         kk1=k2+maxplot*(kk-2)
         call multitop(5000+kk1+ 1 ,4999,2,3,'l+ pt      ',' ','LOG')
         call multitop(5000+kk1+ 2 ,4999,2,3,'ll pt      ',' ','LOG')
         call multitop(5000+kk1+ 3 ,4999,2,3,'ll m       ',' ','LOG')
         call multitop(5000+kk1+ 4 ,4999,2,3,'ll EE      ',' ','LOG')
         call multitop(5000+kk1+ 5 ,4999,2,3,'pTl1+pTl2  ',' ','LOG')
c Mellin
         call multitop(5000+kk1+ 6 ,4999,2,3,'l+ pt      ',' ','LOG')
         call multitop(5000+kk1+ 7 ,4999,2,3,'ll pt      ',' ','LOG')
         call multitop(5000+kk1+ 8 ,4999,2,3,'ll m       ',' ','LOG')
         call multitop(5000+kk1+ 9 ,4999,2,3,'ll EE      ',' ','LOG')
         call multitop(5000+kk1+ 10,4999,2,3,'pTl1+pTl2  ',' ','LOG')
         call multitop(5000+kk1+ 11,4999,2,3,'total      ',' ','LOG')
       enddo
      enddo
c
      CLOSE(99)
      END

C----------------------------------------------------------------------
      SUBROUTINE HWANAL
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
C----------------------------------------------------------------------
      INCLUDE 'HERWIG65.INC'
      include 'reweight0.inc'
      integer nwgt_analysis,max_weight
      common/c_analysis/nwgt_analysis
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
      double precision ww(max_weight),www(max_weight)
      common/cww/ww
      DOUBLE PRECISION HWVDOT,PSUM(4)
      INTEGER ICHSUM,ICHINI,IHEP
      LOGICAL DIDSOF
      LOGICAL BHADRN,BMESON,BBARYON
      INTEGER NL,NN,IST,ID,I,J,IBHAD,mjet,mjet30,mbjet,mbjet30,
     # ILP,ILN,kk,ie,IL(30),jbjmax(30)
      DOUBLE PRECISION ptj1,ptj2,ptjb1,ptjb2,etaj1,etaj2,
     # etajb1,etajb2,dr1,dr2,ptpairlep,detalep,dphilep,dRlep,xmlep,
     # getdr,getptv4,getpseudorapv4,getdelphiv4,getdrv4,getinvmv4,
     # sumEE,sumPT,var,wgt0,PLEP(4,30),ptlep(30),etalep(30),
     # pairlep(4),obs(30)
      integer maxset,maxplot,STEP,Nscales,Npdfpairs,
     # j2,k2,isc,jsc,k1,iobs,ibin,kk1,i2,nps,nng
      common/cnumplots/maxset,maxplot,STEP,Nscales,Npdfpairs

      LOGICAL acceptance,isolation
c jet stuff
      INTEGER NMAX
      PARAMETER (NMAX=4000)
      INTEGER NERANGE
      PARAMETER (NERANGE=1)
      LOGICAL ISAB(NMAX),ISABJET(NMAX),KEEPJ(NMAX,NERANGE)
      INTEGER NJET,JET(NMAX),IPOS(NMAX),NUMBINJET(NMAX),
     # jjmax(30,NERANGE)
      DOUBLE PRECISION PALG,RFJ,SYCUT,PP(4,NMAX),PJET(4,NMAX),
     # PTJET(NMAX),ETAJET(NMAX)
      REAL*8 WWW0
      INTEGER IVLEP1,IVLEP2
      COMMON/VVLIN/IVLEP1,IVLEP2
c pseudorapidity ranges defined as follows:
c   etarangeL(ie)<|eta|<etarangeU(ie),    ie=1,NERANGE
      double precision etarangeL(1),etarangeU(1)
      data etarangeL/0.0d0/
      data etarangeU/2.4d0/
c
      IF (IERROR.NE.0) RETURN
      DO I=1,nwgt_analysis
         WWW(I)=EVWGT*ww(i)/ww(1)
      ENDDO
c
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
C EFFECT, SO THROW THE EVENT AWAY
      IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN
        CALL HWWARN('HWANAL',111)
        GOTO 999
      ENDIF
      WWW0=EVWGT
      CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PSUM)
      CALL HWVSCA(4,-1D0,PSUM,PSUM)
      ICHSUM=0
      ICHINI=ICHRG(IDHW(1))+ICHRG(IDHW(2))
      DIDSOF=.FALSE.
      NL=0
      NN=0
      IBHAD=0
      DO 100 IHEP=1,NHEP
        IF (IDHW(IHEP).EQ.16) DIDSOF=.TRUE.
        IF (ISTHEP(IHEP).EQ.1) THEN
          CALL HWVSUM(4,PHEP(1,IHEP),PSUM,PSUM)
          ICHSUM=ICHSUM+ICHRG(IDHW(IHEP))
        ENDIF
        IST=ISTHEP(IHEP)
        ID=IDHEP(IHEP)
        CALL BHAD(ID,BHADRN,BMESON,BBARYON)
c Assume the relevant B-hadrons have been set stable in the driver
        IF((ABS(ID).EQ.11.OR.ABS(ID).EQ.13).AND.IST.EQ.1)THEN
          NL=NL+1
          IL(NL)=IHEP
          DO I=1,4
             PLEP(I,NL)=PHEP(I,IHEP)
          ENDDO
        ELSEIF(ABS(ID).GT.100.AND.IST.EQ.1) THEN
          NN=NN+1
          IF (NN.GT.NMAX) STOP 'Too many particles [hadrons]!'
          IPOS(NN)=IHEP
          DO I=1,4
             PP(I,NN)=PHEP(I,IHEP)
          ENDDO
        ENDIF
        IF(BHADRN.AND.IST.EQ.1)THEN
          IBHAD=IBHAD+1
          ISAB(IHEP)=.TRUE.
        ELSE
          ISAB(IHEP)=.FALSE.
        ENDIF
  100 CONTINUE
C CHECK MOMENTUM AND CHARGE CONSERVATION
      IF (HWVDOT(3,PSUM,PSUM).GT.1.d-4*PHEP(4,1)**2) THEN
         CALL HWUEPR
         CALL HWWARN('HWANAL',112)
         GOTO 999
      ENDIF
      IF (ICHSUM.NE.ICHINI) THEN
         CALL HWUEPR
         CALL HWWARN('HWANAL',113)
         GOTO 999
      ENDIF
c
      IF( (NL.LT.2.OR.IBHAD.LT.2).AND.
     #    IERROR.EQ.0) THEN
         CALL HWUEPR
         WRITE(*,*)NL,IBHAD
         CALL HWWARN('HWANAL',503)
      ENDIF
C---CLUSTER THE EVENT
      palg=-1.d0
      rfj=0.5d0
      sycut=0.d0
      call fastjetppgenkt(pp,nn,rfj,sycut,palg,pjet,njet,jet)
c Check order in pt, and keep only jet with |eta|<2.4
      mjet=0
      mjet30=0
      do i=1,njet
        ptjet(i)=getptv4(pjet(1,i))
        if(i.gt.1)then
          if (ptjet(i).gt.ptjet(i-1)) then
            write (*,*) "Error in fastjet: jets should be ordered in pt"
            CALL HWWARN('HWANAL',501)
          endif
        endif
        numbinjet(i)=0
        etajet(i)=getpseudorapv4(pjet(1,i))
        do ie=1,NERANGE
          keepj(i,ie)=.false.
          if( abs(etajet(i)).ge.etarangeL(ie) .and.
     #        abs(etajet(i)).lt.etarangeU(ie) )then
            keepj(i,ie)=.true.
            if(ie.eq.1)then
              mjet=mjet+1
              if(ptjet(i).ge.30.d0)mjet30=mjet30+1
            endif
          endif
        enddo
      enddo
c Check whether a jet is a light or a b jet.
c Do so only for jets within pseudorapidity acceptance
      do i=1,nn
        if(jet(i).ne.0)then
          if(keepj(jet(i),1))then
            ID=IDHEP(IPOS(i))
            CALL BHAD(ID,BHADRN,BMESON,BBARYON)
            if(BHADRN)then
              numbinjet(jet(i))=numbinjet(jet(i))+1
              if(.not.ISAB(IPOS(i)))then
                CALL HWWARN('HWANAL',502)
              endif
            endif
          endif
        elseif(jet(i).gt.njet)then
          CALL HWWARN('HWANAL',504)
        endif
      enddo
c
      mbjet=0
      mbjet30=0
c The two hardest b jets within acceptance. Keep in mind the jets are ordered
c in pt, so these are just the first two encountered
      jbjmax(1)=0
      jbjmax(2)=0
      do i=1,njet
        if(keepj(i,1))then
          if(numbinjet(i).eq.0)then
            ISABJET(i)=.FALSE.
          else
            ISABJET(i)=.TRUE.
            mbjet=mbjet+1
            if(jbjmax(2).eq.0.and.jbjmax(1).ne.0)jbjmax(2)=i
            if(jbjmax(1).eq.0)jbjmax(1)=i
            if(ptjet(i).ge.30.d0)mbjet30=mbjet30+1
          endif
        endif
      enddo
c Sanity checks
      if( (jbjmax(1).eq.jbjmax(2).and.jbjmax(1).ne.0) .or.
     #    (jbjmax(1).eq.0.and.jbjmax(2).ne.0) .or.
     #    (mbjet.ge.2.and.jbjmax(2).eq.0) .or.
     #    (mbjet.eq.1.and.jbjmax(2).ne.0) .or.
     #    (mbjet.eq.1.and.jbjmax(1).eq.0) )then
        write(*,*)njet,mjet,mbjet
        write(*,*)jbjmax(1),jbjmax(2)
        CALL HWWARN('HWANAL',504)
      endif
c Tag the first two jets different from the hardest b jets -- these
c can be light jets, but also softer b jets. Also check isolation
c from leptons of all jets harder than 30 GeV
      isolation=.true.
      do ie=1,NERANGE
        jjmax(1,ie)=0
        jjmax(2,ie)=0
        do i=1,njet
          if(keepj(i,ie).and.jbjmax(1).ne.i.and.jbjmax(2).ne.i)then
            if(jjmax(2,ie).eq.0.and.jjmax(1,ie).ne.0)jjmax(2,ie)=i
            if(jjmax(1,ie).eq.0)jjmax(1,ie)=i
            if(ptjet(i).ge.30.d0)then
              if(ie.eq.1)then
                dr1=getdr(PLEP(4,1),PLEP(1,1),PLEP(2,1),PLEP(3,1),
     #                    pjet(4,i),pjet(1,i),pjet(2,i),pjet(3,i))
                dr2=getdr(PLEP(4,2),PLEP(1,2),PLEP(2,2),PLEP(3,2),
     #                    pjet(4,i),pjet(1,i),pjet(2,i),pjet(3,i))
                isolation=isolation .and.
     #                    (dr1.gt.0.4d0).and.(dr2.gt.0.4d0)
              endif
            endif
          endif
        enddo
      enddo
c
      if(jbjmax(1).ne.0)then
        ptjb1=ptjet(jbjmax(1))
        etajb1=etajet(jbjmax(1))
      else
        ptjb1=-1.d8
        etajb1=-1.d8
      endif
      if(jbjmax(2).ne.0)then
        ptjb2=ptjet(jbjmax(2))
        etajb2=etajet(jbjmax(2))
      else
        ptjb2=-1.d8
        etajb2=-1.d8
      endif
      if(jjmax(1,1).ne.0)then
        ptj1=ptjet(jjmax(1,1))
        etaj1=etajet(jjmax(1,1))
      else
        ptj1=-1.d8
        etaj1=-1.d8
      endif
      if(jjmax(2,1).ne.0)then
        ptj2=ptjet(jjmax(2,1))
        etaj2=etajet(jjmax(2,1))
      else
        ptj2=-1.d8
        etaj2=-1.d8
      endif
c Use the first two leptons in the event record
      sumEE=0.d0
      sumPT=0.d0
      do i=1,2
        ptlep(i)=getptv4(plep(1,i))
        etalep(i)=getpseudorapv4(plep(1,i))
        sumEE=sumEE+plep(4,i)
        sumPT=sumPT+ptlep(i)
      enddo

      if(IDHEP(IL(1)).LT.0)then
        ILP=1
        ILN=2
      else
        ILP=2
        ILN=1
      endif
      do j=1,4
        pairlep(j)=plep(j,1)+plep(j,2)
      enddo
      ptpairlep=getptv4(pairlep)
      dphilep=getdelphiv4(plep(1,1),plep(1,2))
      xmlep=getinvmv4(pairlep)
      detalep=etalep(ILP)-etalep(ILN)
      dRlep=getdrv4(plep(1,1),plep(1,2))
c
      acceptance=
     #  ptlep(1).ge.20.d0.and.abs(etalep(1)).le.2.4d0 .and.
     #  ptlep(2).ge.20.d0.and.abs(etalep(2)).le.2.4d0 .and.
     #  ptjb1.ge.30.d0.and.abs(etajb1).le.2.4d0 .and.
     #  ptjb2.ge.30.d0.and.abs(etajb2).le.2.4d0
c
      obs(1)=ptlep(ILP)
      obs(2)=ptpairlep
      obs(3)=xmlep
      obs(4)=sumEE
      obs(5)=sumPT
c loop over each setup
      do j2=1,maxset
      k2=(j2-1)*STEP
c loop over scales and PDFs; skip the first which is equal to (1,1)
       do kk=2,nwgt_analysis
        kk1=k2+maxplot*(kk-2)
c NOTE: '3' == maxset !
        if ((j2.eq.1).or.
     #     ((j2.eq.2).and.acceptance).or.
     #     ((j2.eq.3).and.acceptance.and.isolation)) then
         WWW0=WWW(kk)
c
         do ibin=0,4
          k1=kk1+11
          call mfill(k1,dfloat(ibin),WWW0)
          do iobs=1,5
            if(ibin.eq.0)then
              var=1.d0
            else
              var=obs(iobs)**ibin
            endif
            wgt0=WWW0*var
            k1=kk1+iobs+5
            call mfill(k1,dfloat(ibin),wgt0)
            if(ibin.eq.0)then
c Fill observables just once
              k1=kk1+iobs
              call mfill(k1,obs(iobs),WWW0)
            endif
          enddo
         enddo
c
        endif
       enddo
      enddo
C
 999  RETURN
      END

C----------------------------------------------------------------------
      SUBROUTINE BHAD(IDPDG,BHADRN,BMESON,BBARYON)
C     TEST FOR A B-FLAVOURED HADRON
C----------------------------------------------------------------------
      INTEGER IDPDG,ID,IM,IB
      LOGICAL BHADRN,BMESON,BBARYON
C
      ID=ABS(IDPDG)
      IM=MOD(ID/100,100)
      IB=MOD(ID/1000,100)
      BMESON=IM.EQ.5
      BBARYON=IB.EQ.5
      IF(BMESON.AND.BBARYON)CALL HWWARN('BHAD  ',500)
      BHADRN=BMESON.OR.BBARYON
 999  END

      FUNCTION IHADR(ID)
c Returns the PDG code of the heavier quark in the hadron of PDG code ID
      IMPLICIT NONE
      INTEGER IHADR,ID,ID1
C
      IF(ID.NE.0)THEN
        ID1=ABS(ID)
        IF(ID1.GT.10000)ID1=ID1-1000*INT(ID1/1000)
        IHADR=ID1/(10**INT(LOG10(DFLOAT(ID1))))
      ELSE
        IHADR=0
      ENDIF
      RETURN
      END

      function getrapidity(en,pl)
      implicit none
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
      parameter (tiny=1.d-8)
c
      xplus=en+pl
      xminus=en-pl
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
        if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
          y=0.5d0*log( xplus/xminus )
        else
          y=sign(1.d0,pl)*1.d8
        endif
      else
        y=sign(1.d0,pl)*1.d8
      endif
      getrapidity=y
      return
      end

      function getpseudorap(en,ptx,pty,pl)
      implicit none
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
      parameter (tiny=1.d-5)
c
      pt=sqrt(ptx**2+pty**2)
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
        eta=sign(1.d0,pl)*1.d8
      else
        th=atan2(pt,pl)
        eta=-log(tan(th/2.d0))
      endif
      getpseudorap=eta
      return
      end

      function getinvm(en,ptx,pty,pl)
      implicit none
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
      parameter (tiny=1.d-5)
c
      tmp=en**2-ptx**2-pty**2-pl**2
      if(tmp.gt.0.d0)then
        tmp=sqrt(tmp)
      elseif(tmp.gt.-tiny)then
        tmp=0.d0
      else
        write(*,*)'Attempt to compute a negative mass'
        stop
      endif
      getinvm=tmp
      return
      end

      function getdelphi(ptx1,pty1,ptx2,pty2)
      implicit none
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
      parameter (tiny=1.d-5)
c
      pt1=sqrt(ptx1**2+pty1**2)
      pt2=sqrt(ptx2**2+pty2**2)
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
        tmp=ptx1*ptx2+pty1*pty2
        tmp=tmp/(pt1*pt2)
        if(abs(tmp).gt.1.d0+tiny)then
          write(*,*)'Cosine larger than 1'
          stop
        elseif(abs(tmp).ge.1.d0)then
          tmp=sign(1.d0,tmp)
        endif
        tmp=acos(tmp)
      else
        tmp=1.d8
      endif
      getdelphi=tmp
      return
      end

      function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
      implicit none
      real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
     # getpseudorap,getdelphi
c
      deta=getpseudorap(en1,ptx1,pty1,pl1)-
     #     getpseudorap(en2,ptx2,pty2,pl2)
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
      getdr=sqrt(dphi**2+deta**2)
      return
      end

      function getdry(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
      implicit none
      real*8 getdry,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
     # getrapidity,getdelphi
c
      deta=getrapidity(en1,pl1)-
     #     getrapidity(en2,pl2)
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
      getdry=sqrt(dphi**2+deta**2)
      return
      end

      function getptv(p)
      implicit none
      real*8 getptv,p(5)
c
      getptv=sqrt(p(1)**2+p(2)**2)
      return
      end

      function getpseudorapv(p)
      implicit none
      real*8 getpseudorapv,p(5)
      real*8 getpseudorap
c
      getpseudorapv=getpseudorap(p(4),p(1),p(2),p(3))
      return
      end

      function getrapidityv(p)
      implicit none
      real*8 getrapidityv,p(5)
      real*8 getrapidity
c
      getrapidityv=getrapidity(p(4),p(3))
      return
      end

      function getinvmv(p)
      implicit none
      real*8 getinvmv,p(5)
      real*8 getinvm
c
      getinvmv=getinvm(p(4),p(1),p(2),p(3))
      return
      end

      function getptv4(p)
      implicit none
      real*8 getptv4,p(4)
c
      getptv4=sqrt(p(1)**2+p(2)**2)
      return
      end

      function getpseudorapv4(p)
      implicit none
      real*8 getpseudorapv4,p(4)
      real*8 getpseudorap
c
      getpseudorapv4=getpseudorap(p(4),p(1),p(2),p(3))
      return
      end

      function getrapidityv4(p)
      implicit none
      real*8 getrapidityv4,p(4)
      real*8 getrapidity
c
      getrapidityv4=getrapidity(p(4),p(3))
      return
      end

      function getdrv4(p1,p2)
      implicit none
      real*8 getdrv4,p1(4),p2(4)
      real*8 getdr
c
      getdrv4=getdr(p1(4),p1(1),p1(2),p1(3),
     #              p2(4),p2(1),p2(2),p2(3))
      return
      end

      function getinvmv4(p)
      implicit none
      real*8 getinvmv4,p(4)
      real*8 getinvm
c
      getinvmv4=getinvm(p(4),p(1),p(2),p(3))
      return
      end

      function getdelphiv4(p1,p2)
      implicit none
      real*8 getdelphiv4,p1(4),p2(4)
      real*8 getdelphi
c
      getdelphiv4=getdelphi(p1(1),p1(2),
     #                      p2(1),p2(2))
      return
      end

      function getcosv4(q1,q2)
      implicit none
      real*8 getcosv4,q1(4),q2(4)
      real*8 xnorm1,xnorm2,tmp
c
      if(q1(4).lt.0.d0.or.q2(4).lt.0.d0)then
        getcosv4=-1.d10
        return
      endif
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
        tmp=-1.d10
      else
        tmp=q1(1)*q2(1)+q1(2)*q2(2)+q1(3)*q2(3)
        tmp=tmp/(xnorm1*xnorm2)
        if(abs(tmp).gt.1.d0.and.abs(tmp).le.1.001d0)then
          tmp=sign(1.d0,tmp)
        elseif(abs(tmp).gt.1.001d0)then
          write(*,*)'Error in getcosv4',tmp
          stop
        endif
      endif
      getcosv4=tmp
      return
      end

      subroutine getperpenv4(q1,q2,qperp)
c Normal to the plane defined by \vec{q1},\vec{q2}
      implicit none
      real*8 q1(4),q2(4),qperp(4)
      real*8 xnorm1,xnorm2
      integer i
c
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
        do i=1,4
          qperp(i)=-1.d10
        enddo
      else
        qperp(1)=q1(2)*q2(3)-q1(3)*q2(2)
        qperp(2)=q1(3)*q2(1)-q1(1)*q2(3)
        qperp(3)=q1(1)*q2(2)-q1(2)*q2(1)
        do i=1,3
          qperp(i)=qperp(i)/(xnorm1*xnorm2)
        enddo
        qperp(4)=1.d0
      endif
      return
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
module comoving_module

    use amrex_fort_module, only : rt => amrex_real

    integer,  save :: comoving_type
    real(rt), save :: comoving_OmM, comoving_OmB, comoving_OmN, comoving_h

end module comoving_module

max_stars_repo_name: burlen/Nyx

max_issues_repo_name: burlen/Nyx

max_forks_repo_name: burlen/Nyx


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.
C*******************************************************************************
C
      integer function SLEN(STRING)
C !DEC$ ATTRIBUTES DLLEXPORT::SLEN
C*******************************************************************************
C
C    $Id: slen.f,v 1.1 2002/05/02 18:52:20 andreasw Exp $
C
C-------------------------------------------------------------------------------
C                                 Title
C-------------------------------------------------------------------------------
C
CT    Simple function that determines the length of a FORTRAN string
C
C-------------------------------------------------------------------------------
C                             Author, date
C-------------------------------------------------------------------------------
C
CA    Andreas Waechter      05/01/02  Release as version IPOPT 2.0
C
C-------------------------------------------------------------------------------
C                             Parameter list
C-------------------------------------------------------------------------------
C
C     Name     I/O   Type   Meaning
CP    STRING    I    C*(*)  String whos length is to be determined
CP    SLEN      O    I      Position of last non-space character in STRING
C
C*******************************************************************************
C
C                              Declarations
C
C*******************************************************************************
C
      IMPLICIT NONE
C
C-------------------------------------------------------------------------------
C                             Parameter list
C-------------------------------------------------------------------------------
C
      character*(*) STRING
C
C-------------------------------------------------------------------------------
C                            Local variables
C-------------------------------------------------------------------------------
C
      integer i
C
C*******************************************************************************
C
C                           Executable Statements
C
C*******************************************************************************
C
      if( STRING.eq.' ' ) then
         SLEN = 0
      else
         do SLEN = len(STRING), 1, -1
            if( STRING(slen:slen).ne.' ' ) goto 150
         enddo
      endif
 150  continue
      return
      end

ext: f

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
PROGRAM ApplicationDriver

  USE KindModule, ONLY: &
    DP, Half, One, Two, &
    Pi, TwoPi, FourPi
  USE ProgramHeaderModule, ONLY: &
    iX_B0, iX_B1, iX_E0, iX_E1, &
    nDimsX, nDOFX
  USE ProgramInitializationModule, ONLY: &
    InitializeProgram, &
    FinalizeProgram
  USE ReferenceElementModuleX, ONLY: &
    InitializeReferenceElementX, &
    FinalizeReferenceElementX
  USE ReferenceElementModuleX_Lagrange, ONLY: &
    InitializeReferenceElementX_Lagrange, &
    FinalizeReferenceElementX_Lagrange
  USE InputOutputModuleHDF, ONLY: &
    WriteFieldsHDF
  USE GeometryComputationModule, ONLY: &
    ComputeGeometryX
  USE GravitySolutionModule_Newtonian_PointMass, ONLY: &
    ComputeGravitationalPotential
  USE FluidFieldsModule, ONLY: &
    nCF, nPF, nAF, &
    uCF, uPF, uAF, &
    uDF
  USE TimeSteppingModule_SSPRK, ONLY: &
    InitializeFluid_SSPRK, &
    FinalizeFluid_SSPRK, &
    UpdateFluid_SSPRK
  USE Euler_SlopeLimiterModule_NonRelativistic_IDEAL, ONLY: &
    InitializeSlopeLimiter_Euler_NonRelativistic_IDEAL, &
    FinalizeSlopeLimiter_Euler_NonRelativistic_IDEAL, &
    ApplySlopeLimiter_Euler_NonRelativistic_IDEAL
  USE Euler_PositivityLimiterModule_NonRelativistic_IDEAL, ONLY: &
    InitializePositivityLimiter_Euler_NonRelativistic_IDEAL, &
    FinalizePositivityLimiter_Euler_NonRelativistic_IDEAL, &
    ApplyPositivityLimiter_Euler_NonRelativistic_IDEAL
  USE GeometryFieldsModule, ONLY: &
    nGF, uGF
  USE EquationOfStateModule, ONLY: &
    InitializeEquationOfState, &
    FinalizeEquationOfState
  USE InitializationModule, ONLY: &
    InitializeFields, ApplyPerturbations
  USE Euler_UtilitiesModule, ONLY: &
    ComputeFromConserved_Euler, &
    ComputeTimeStep_Euler
  USE Euler_dgDiscretizationModule, ONLY: &
    ComputeIncrement_Euler_DG_Explicit
  USE Euler_TallyModule_NonRelativistic, ONLY: &
    InitializeTally_Euler_NonRelativistic, &
    FinalizeTally_Euler_NonRelativistic, &
    ComputeTally_Euler_NonRelativistic

  IMPLICIT NONE

  INCLUDE 'mpif.h'

  REAL(DP), PARAMETER :: mDot   = FourPi
  REAL(DP), PARAMETER :: Mass   = Half
  REAL(DP), PARAMETER :: rShock = One
  REAL(DP), PARAMETER :: Gamma  = 4.0_DP / 3.0_DP
  REAL(DP), PARAMETER :: Mach   = 1.0d2

  LOGICAL  :: wrt
  INTEGER  :: iCycle, iCycleD
  INTEGER  :: nX(3), nNodes
  REAL(DP) :: t, dt, t_end, dt_wrt, t_wrt, wTime
  INTEGER  :: iErr

  nX     = [ 256, 1, 1 ]
  nNodes = 2

  CALL InitializeProgram &
         ( ProgramName_Option &
             = 'StandingAccretionShock', &
           nX_Option &
             = nX, &
           swX_Option &
             = [ 1, 0, 0 ], &
           bcX_Option &
             = [ 11, 3, 1 ], &
           xL_Option &
             = [ 0.2d0, 0.0d0, 0.0d0 ], &
           xR_Option &
             = [ 2.0d0, Pi,    TwoPi ], &
           zoomX_Option &
             = [ 1.0025_DP, 1.0_DP, 1.0_DP ], &
           nNodes_Option &
             = nNodes, &
           CoordinateSystem_Option &
             = 'SPHERICAL', &
           BasicInitialization_Option &
             = .TRUE. )

  t_end   = 5.0d+1
  dt_wrt  = 2.5d-1
  iCycleD = 100

  CALL InitializeReferenceElementX

  CALL InitializeReferenceElementX_Lagrange

  CALL ComputeGeometryX &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF )

  CALL ComputeGravitationalPotential &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, Mass )

  CALL InitializeFluid_SSPRK( nStages = 2 )

  CALL InitializeEquationOfState &
         ( EquationOfState_Option = 'IDEAL', &
           Gamma_IDEAL_Option = Gamma )

  CALL InitializeSlopeLimiter_Euler_NonRelativistic_IDEAL &
         ( BetaTVD_Option                   = 1.15_DP, &
           BetaTVB_Option                   = 0.00_DP, &
           SlopeTolerance_Option            = 1.00d-6, &
           UseSlopeLimiter_Option           = .TRUE., &
           UseCharacteristicLimiting_Option = .TRUE., &
           UseTroubledCellIndicator_Option  = .TRUE., &
           UseConservativeCorrection_Option = .TRUE., &
           LimiterThresholdParameter_Option = 0.015_DP )

  CALL InitializePositivityLimiter_Euler_NonRelativistic_IDEAL &
         ( UsePositivityLimiter_Option = .TRUE., &
           Verbose_Option = .TRUE., &
           Min_1_Option = 1.0d-12, &
           Min_2_Option = 1.0d-12 )

  CALL InitializeFields &
         ( mDot, Mass, rShock, Gamma, Mach )

  CALL ApplyPerturbations( 1.4_DP, 1.6_DP, 0, 2.0_DP )

  CALL ApplySlopeLimiter_Euler_NonRelativistic_IDEAL &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, uDF )

  CALL ApplyPositivityLimiter_Euler_NonRelativistic_IDEAL &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF )

  CALL ComputeFromConserved_Euler &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, uPF, uAF )

  CALL WriteFieldsHDF &
         ( 0.0_DP, WriteGF_Option = .TRUE., WriteFF_Option = .TRUE. )

  ! --- Evolve ---

  wTime = MPI_WTIME( )

  t     = 0.0_DP
  t_wrt = dt_wrt
  wrt   = .FALSE.

  CALL InitializeTally_Euler_NonRelativistic &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF )

  CALL ComputeTally_Euler_NonRelativistic &
       ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, Time = t, &
         SetInitialValues_Option = .TRUE., Verbose_Option = .TRUE. )

  iCycle = 0
  DO WHILE ( t < t_end )

    iCycle = iCycle + 1

    CALL ComputeTimeStep_Euler &
           ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, &
             CFL = 0.3_DP / ( nDimsX * DBLE( nNodes - 1 ) + One ), &
             TimeStep = dt )

    IF( t + dt > t_end )THEN

      dt = t_end - t

    END IF

    IF( t + dt > t_wrt )THEN

      dt    = t_wrt - t
      t_wrt = t_wrt + dt_wrt
      wrt   = .TRUE.

    END IF

    IF( MOD( iCycle, iCycleD ) == 0 )THEN

      WRITE(*,'(A8,A8,I8.8,A2,A4,ES13.6E3,A1,A5,ES13.6E3)') &
          '', 'Cycle = ', iCycle, '', 't = ',  t, '', 'dt = ', dt

    END IF

    CALL UpdateFluid_SSPRK &
           ( t, dt, uGF, uCF, uDF, ComputeIncrement_Euler_DG_Explicit )

    t = t + dt

    IF( wrt )THEN

      CALL ComputeFromConserved_Euler &
             ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, uPF, uAF )

      CALL WriteFieldsHDF &
             ( t, WriteGF_Option = .TRUE., WriteFF_Option = .TRUE. )

      CALL ComputeTally_Euler_NonRelativistic &
             ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, Time = t, &
               Verbose_Option = .TRUE. )

      wrt = .FALSE.

    END IF

  END DO

  CALL ComputeFromConserved_Euler &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, uPF, uAF )

  CALL WriteFieldsHDF &
         ( t, WriteGF_Option = .TRUE., WriteFF_Option = .TRUE. )

  CALL ComputeTally_Euler_NonRelativistic &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, Time = t, &
           Verbose_Option = .TRUE. )

  CALL FinalizeTally_Euler_NonRelativistic

  wTime = MPI_WTIME( ) - wTime

  WRITE(*,*)
  WRITE(*,'(A6,A,I8.8,A,ES12.6E2,A)') &
    '', 'Finished ', iCycle, ' Cycles in ', wTime, ' s'
  WRITE(*,*)

  ! --- Finalize ---

  CALL FinalizePositivityLimiter_Euler_NonRelativistic_IDEAL

  CALL FinalizeSlopeLimiter_Euler_NonRelativistic_IDEAL

  CALL FinalizeEquationOfState

  CALL FinalizeFluid_SSPRK

  CALL FinalizeReferenceElementX_Lagrange

  CALL FinalizeReferenceElementX

  CALL FinalizeProgram

  WRITE(*,*)
  WRITE(*,'(2x,A)') 'git info'
  WRITE(*,'(2x,A)') '--------'
  WRITE(*,*)
  WRITE(*,'(2x,A)') 'git branch:'
  CALL EXECUTE_COMMAND_LINE( 'git branch' )
  WRITE(*,*)
  WRITE(*,'(2x,A)') 'git describe --tags:'
  CALL EXECUTE_COMMAND_LINE( 'git describe --tags' )
  WRITE(*,*)
  WRITE(*,'(2x,A)') 'git rev-parse HEAD:'
  CALL EXECUTE_COMMAND_LINE( 'git rev-parse HEAD' )
  WRITE(*,*)
  WRITE(*,'(2x,A)') 'date:'
  CALL EXECUTE_COMMAND_LINE( 'date' )
  WRITE(*,*)

END PROGRAM ApplicationDriver

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


! { dg-do run }
!
! PR 45521: [F08] GENERIC resolution with ALLOCATABLE/POINTER and PROCEDURE
!
! Contributed by <wangmianzhi1@linuxmail.org>

  interface test
    procedure testAlloc
    procedure testPtr
  end interface

  integer, allocatable :: a1
  integer, pointer :: a2

  if (.not.test(a1)) STOP 1
  if (test(a2)) STOP 2

contains

  logical function testAlloc(obj)
    integer, allocatable :: obj
    testAlloc = .true.
  end function

  logical function testPtr(obj)
    integer, pointer :: obj
    testPtr = .false.
  end function

end


text:
*DECK ZUNI2
      SUBROUTINE ZUNI2 (ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
     +   TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  ZUNI2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to ZBESI and ZBESK
C***LIBRARY   SLATEC
C***TYPE      ALL (CUNI2-A, ZUNI2-A)
C***AUTHOR  Amos, D. E., (SNL)
C***DESCRIPTION
C
C     ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
C     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
C     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
C
C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
C     Y(I)=CZERO FOR I=NLAST+1,N
C
C***SEE ALSO  ZBESI, ZBESK
C***ROUTINES CALLED  D1MACH, ZABS, ZAIRY, ZUCHK, ZUNHJ, ZUOIK
C***REVISION HISTORY  (YYMMDD)
C   830501  DATE WRITTEN
C   910415  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  ZUNI2
C     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
C    *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
      DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
     * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR,
     * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII,
     * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI,
     * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI,
     * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR,
     * CYI, D1MACH, ZABS, CAR, SAR
      INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
     * NN, NUF, NW, NZ, IDUM
      DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3),
     * CSRR(3), CYR(2), CYI(2)
      EXTERNAL ZABS
      DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
      DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
     * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/
      DATA HPI, AIC  /
     1      1.57079632679489662D+00,     1.265512123484645396D+00/
C***FIRST EXECUTABLE STATEMENT  ZUNI2
      NZ = 0
      ND = N
      NLAST = 0
C-----------------------------------------------------------------------
C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
C     EXP(ALIM)=EXP(ELIM)*TOL
C-----------------------------------------------------------------------
      CSCL = 1.0D0/TOL
      CRSC = TOL
      CSSR(1) = CSCL
      CSSR(2) = CONER
      CSSR(3) = CRSC
      CSRR(1) = CRSC
      CSRR(2) = CONER
      CSRR(3) = CSCL
      BRY(1) = 1.0D+3*D1MACH(1)/TOL
C-----------------------------------------------------------------------
C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
C-----------------------------------------------------------------------
      ZNR = ZI
      ZNI = -ZR
      ZBR = ZR
      ZBI = ZI
      CIDI = -CONER
      INU = FNU
      ANG = HPI*(FNU-INU)
      C2R = COS(ANG)
      C2I = SIN(ANG)
      CAR = C2R
      SAR = C2I
      IN = INU + N - 1
      IN = MOD(IN,4) + 1
      STR = C2R*CIPR(IN) - C2I*CIPI(IN)
      C2I = C2R*CIPI(IN) + C2I*CIPR(IN)
      C2R = STR
      IF (ZI.GT.0.0D0) GO TO 10
      ZNR = -ZNR
      ZBI = -ZBI
      CIDI = -CIDI
      C2I = -C2I
   10 CONTINUE
C-----------------------------------------------------------------------
C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
C-----------------------------------------------------------------------
      FN = MAX(FNU,1.0D0)
      CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
     * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
      IF (KODE.EQ.1) GO TO 20
      STR = ZBR + ZETA2R
      STI = ZBI + ZETA2I
      RAST = FN/ZABS(STR,STI)
      STR = STR*RAST*RAST
      STI = -STI*RAST*RAST
      S1R = -ZETA1R + STR
      S1I = -ZETA1I + STI
      GO TO 30
   20 CONTINUE
      S1R = -ZETA1R + ZETA2R
      S1I = -ZETA1I + ZETA2I
   30 CONTINUE
      RS1 = S1R
      IF (ABS(RS1).GT.ELIM) GO TO 150
   40 CONTINUE
      NN = MIN(2,ND)
      DO 90 I=1,NN
        FN = FNU + (ND-I)
        CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI,
     *   ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
        IF (KODE.EQ.1) GO TO 50
        STR = ZBR + ZETA2R
        STI = ZBI + ZETA2I
        RAST = FN/ZABS(STR,STI)
        STR = STR*RAST*RAST
        STI = -STI*RAST*RAST
        S1R = -ZETA1R + STR
        S1I = -ZETA1I + STI + ABS(ZI)
        GO TO 60
   50   CONTINUE
        S1R = -ZETA1R + ZETA2R
        S1I = -ZETA1I + ZETA2I
   60   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        RS1 = S1R
        IF (ABS(RS1).GT.ELIM) GO TO 120
        IF (I.EQ.1) IFLAG = 2
        IF (ABS(RS1).LT.ALIM) GO TO 70
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
        APHI = ZABS(PHIR,PHII)
        AARG = ZABS(ARGR,ARGI)
        RS1 = RS1 + LOG(APHI) - 0.25D0*LOG(AARG) - AIC
        IF (ABS(RS1).GT.ELIM) GO TO 120
        IF (I.EQ.1) IFLAG = 1
        IF (RS1.LT.0.0D0) GO TO 70
        IF (I.EQ.1) IFLAG = 3
   70   CONTINUE
C-----------------------------------------------------------------------
C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
C     EXPONENT EXTREMES
C-----------------------------------------------------------------------
        CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM)
        CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM)
        STR = DAIR*BSUMR - DAII*BSUMI
        STI = DAIR*BSUMI + DAII*BSUMR
        STR = STR + (AIR*ASUMR-AII*ASUMI)
        STI = STI + (AIR*ASUMI+AII*ASUMR)
        S2R = PHIR*STR - PHII*STI
        S2I = PHIR*STI + PHII*STR
        STR = EXP(S1R)*CSSR(IFLAG)
        S1R = STR*COS(S1I)
        S1I = STR*SIN(S1I)
        STR = S2R*S1R - S2I*S1I
        S2I = S2R*S1I + S2I*S1R
        S2R = STR
        IF (IFLAG.NE.1) GO TO 80
        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
        IF (NW.NE.0) GO TO 120
   80   CONTINUE
        IF (ZI.LE.0.0D0) S2I = -S2I
        STR = S2R*C2R - S2I*C2I
        S2I = S2R*C2I + S2I*C2R
        S2R = STR
        CYR(I) = S2R
        CYI(I) = S2I
        J = ND - I + 1
        YR(J) = S2R*CSRR(IFLAG)
        YI(J) = S2I*CSRR(IFLAG)
        STR = -C2I*CIDI
        C2I = C2R*CIDI
        C2R = STR
   90 CONTINUE
      IF (ND.LE.2) GO TO 110
      RAZ = 1.0D0/ZABS(ZR,ZI)
      STR = ZR*RAZ
      STI = -ZI*RAZ
      RZR = (STR+STR)*RAZ
      RZI = (STI+STI)*RAZ
      BRY(2) = 1.0D0/BRY(1)
      BRY(3) = D1MACH(2)
      S1R = CYR(1)
      S1I = CYI(1)
      S2R = CYR(2)
      S2I = CYI(2)
      C1R = CSRR(IFLAG)
      ASCLE = BRY(IFLAG)
      K = ND - 2
      FN = K
      DO 100 I=3,ND
        C2R = S2R
        C2I = S2I
        S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
        S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
        S1R = C2R
        S1I = C2I
        C2R = S2R*C1R
        C2I = S2I*C1R
        YR(K) = C2R
        YI(K) = C2I
        K = K - 1
        FN = FN - 1.0D0
        IF (IFLAG.GE.3) GO TO 100
        STR = ABS(C2R)
        STI = ABS(C2I)
        C2M = MAX(STR,STI)
        IF (C2M.LE.ASCLE) GO TO 100
        IFLAG = IFLAG + 1
        ASCLE = BRY(IFLAG)
        S1R = S1R*C1R
        S1I = S1I*C1R
        S2R = C2R
        S2I = C2I
        S1R = S1R*CSSR(IFLAG)
        S1I = S1I*CSSR(IFLAG)
        S2R = S2R*CSSR(IFLAG)
        S2I = S2I*CSSR(IFLAG)
        C1R = CSRR(IFLAG)
  100 CONTINUE
  110 CONTINUE
      RETURN
  120 CONTINUE
      IF (RS1.GT.0.0D0) GO TO 140
C-----------------------------------------------------------------------
C     SET UNDERFLOW AND UPDATE PARAMETERS
C-----------------------------------------------------------------------
      YR(ND) = ZEROR
      YI(ND) = ZEROI
      NZ = NZ + 1
      ND = ND - 1
      IF (ND.EQ.0) GO TO 110
      CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
      IF (NUF.LT.0) GO TO 140
      ND = ND - NUF
      NZ = NZ + NUF
      IF (ND.EQ.0) GO TO 110
      FN = FNU + (ND-1)
      IF (FN.LT.FNUL) GO TO 130
C      FN = CIDI
C      J = NUF + 1
C      K = MOD(J,4) + 1
C      S1R = CIPR(K)
C      S1I = CIPI(K)
C      IF (FN.LT.0.0D0) S1I = -S1I
C      STR = C2R*S1R - C2I*S1I
C      C2I = C2R*S1I + C2I*S1R
C      C2R = STR
      IN = INU + ND - 1
      IN = MOD(IN,4) + 1
      C2R = CAR*CIPR(IN) - SAR*CIPI(IN)
      C2I = CAR*CIPI(IN) + SAR*CIPR(IN)
      IF (ZI.LE.0.0D0) C2I = -C2I
      GO TO 40
  130 CONTINUE
      NLAST = ND
      RETURN
  140 CONTINUE
      NZ = -1
      RETURN
  150 CONTINUE
      IF (RS1.GT.0.0D0) GO TO 140
      NZ = N
      DO 160 I=1,N
        YR(I) = ZEROR
        YI(I) = ZEROI
  160 CONTINUE
      RETURN
      END


module intaodmod
!$$$ module documentation block
!                .      .    .                                       .
! module:    intaodmod    module for intaod and its tangent linear intaod_tl
!   prgmmr: hclin            org: ncar/mmm         date: 2010-10-20
!
! abstract: module for intaod
!
! program history log:
!   2010-10-20  hclin   - modified from intrad for total aod
!   2016-05-18  guo     - replaced ob_type with polymorphic obsNode through type casting
!
! subroutines included:
!   sub intaod_
!
! variable definitions:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use m_obsNode , only: obsNode
  use m_aeroNode , only: aeroNode
  use m_aeroNode , only: aeroNode_typecast
  use m_aeroNode , only: aeroNode_nextcast
  implicit none

  private
  public :: intaod

  interface intaod; module procedure &
            intaod_
  end interface

contains

  subroutine intaod_(aerohead,rval,sval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intaod
!   prgmmr: hclin          org: ncar/mmm                date: 2010-10-20
!
! abstract: apply aod operator and adjoint
!           nonlinear qc is not implemented yet.
!
! program history log:
!   2010-10-20  hclin   - modified from intrad for total aod
!   2014-12-03  derber  - modify so that use of obsdiags can be turned off
!
!   input argument list:
!     aerohead  - obs type pointer to obs structure
!     rval
!     sval
!
!   output argument list:
!     rval
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind,r_quad
    use aeroinfo, only: aerojacnames,nsigaerojac
    use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
    use jfunc, only: jiter
    use gridmod, only: latlon11,nsig
    use constants, only: zero,half,one,r3600
    use gsi_bundlemod, only: gsi_bundle
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use gsi_bundlemod, only: gsi_bundleputvar
    use gsi_chemguess_mod, only: gsi_chemguess_get
    use mpeu_util, only: getindex
    use m_obsdiagNode, only: obsdiagNode_set
    implicit none

! Declare passed variables
    class(obsNode),pointer,intent(in) :: aerohead
    type(gsi_bundle), intent(in   ) :: sval
    type(gsi_bundle), intent(inout) :: rval

! Declare local variables
    integer(i_kind) j1,j2,j3,j4,i1,i2,i3,i4,n,k,ic,nn
    integer(i_kind) istatus,naero
    integer(i_kind),dimension(nsig) :: i1n,i2n,i3n,i4n
    real(r_kind) val
    real(r_kind) w1,w2,w3,w4
!   real(r_kind) cg_aero,p0,wnotgross,wgross
    type(aeroNode), pointer :: aeroptr

    real(r_kind),pointer,dimension(:) :: sv_chem
    real(r_kind),pointer,dimension(:) :: rv_chem
    real(r_kind),dimension(nsigaerojac) :: tval,tdir

  ! Inquire about chemistry
    naero = size(aerojacnames)
    if ( naero <= 0 ) return

    !aeroptr => aerohead
    aeroptr => aeroNode_typecast(aerohead)
    do while (associated(aeroptr))
       j1=aeroptr%ij(1)
       j2=aeroptr%ij(2)
       j3=aeroptr%ij(3)
       j4=aeroptr%ij(4)
       w1=aeroptr%wij(1)
       w2=aeroptr%wij(2)
       w3=aeroptr%wij(3)
       w4=aeroptr%wij(4)

       do k=1,nsigaerojac
          tval(k)=zero
       end do

!  Begin Forward model
!  calculate temperature, q, ozone, sst vector at observation location
       i1n(1) = j1
       i2n(1) = j2
       i3n(1) = j3
       i4n(1) = j4
       do k=2,nsig
          i1n(k) = i1n(k-1)+latlon11
          i2n(k) = i2n(k-1)+latlon11
          i3n(k) = i3n(k-1)+latlon11
          i4n(k) = i4n(k-1)+latlon11
       enddo
       do ic = 1, naero
          call gsi_bundlegetpointer (sval,trim(aerojacnames(ic)),sv_chem,istatus)
          do k=1,nsig
             i1 = i1n(k)
             i2 = i2n(k)
             i3 = i3n(k)
             i4 = i4n(k)
             tdir(k+nsig*(ic-1))=  w1* sv_chem(i1)+w2* sv_chem(i2)+ &
                                   w3* sv_chem(i3)+w4* sv_chem(i4)
          end do
          nullify(sv_chem)
       end do

!  begin channel specific calculations
       do nn=1,aeroptr%nlaero
          ic=aeroptr%icx(nn)

!       include observation increment
          val=zero

!       Include contributions from atmospheric jacobian
          do k=1,nsigaerojac
             val=val+tdir(k)*aeroptr%daod_dvar(k,nn)
          end do

          if(luse_obsdiag)then
             if (lsaveobsens) then
                val = val*aeroptr%err2(nn)*aeroptr%raterr2(nn)
                !-- aeroptr%diags(nn)%ptr%obssen(jiter) = val
                call obsdiagNode_set(aeroptr%diags(nn)%ptr,jiter=jiter,obssen=val)
             else
                !-- if (aeroptr%luse) aeroptr%diags(nn)%ptr%tldepart(jiter) = val
                if (aeroptr%luse) call obsdiagNode_set(aeroptr%diags(nn)%ptr,jiter=jiter,tldepart=val)
             endif
          endif

          if (l_do_adjoint) then
             if ( .not. lsaveobsens) then

                val=val-aeroptr%res(nn)

!             Multiply by variance.
              !if (nlnqc_iter .and. pg_aero(ic) > tiny_r_kind .and. &
              !                     b_aero(ic)  > tiny_r_kind) then
              !   cg_aero=cg_term/b_aero(ic)
              !   wnotgross= one-pg_aero(ic)*varqc_iter
              !   wgross = varqc_iter*pg_aero(ic)*cg_aero/wnotgross
              !   p0   = wgross/(wgross+exp(-half*aeroptr%err2(nn)*val*val))
              !   val = val*(one-p0)
              !endif

                val = val*aeroptr%err2(nn)*aeroptr%raterr2(nn)
             endif

!          Begin adjoint

!          Extract contributions from atmospheric jacobian
             do k=1,nsigaerojac
                tval(k)=tval(k)+aeroptr%daod_dvar(k,nn)*val
             end do

          endif
       end do

       if (l_do_adjoint) then
!    Distribute adjoint contributions over surrounding grid points

          do ic = 1, naero
             call gsi_bundlegetpointer (rval,trim(aerojacnames(ic)),rv_chem,istatus)
             if ( istatus /= 0 ) then
                write(6,*) 'error gsi_bundlegetpointer in intaod for ', aerojacnames(ic)
             end if
             do k=1,nsig
                n = k + nsig*(ic-1)
                i1 = i1n(k)
                i2 = i2n(k)
                i3 = i3n(k)
                i4 = i4n(k)
                rv_chem(i1) = rv_chem(i1) + w1*tval(n)
                rv_chem(i2) = rv_chem(i2) + w2*tval(n)
                rv_chem(i3) = rv_chem(i3) + w3*tval(n)
                rv_chem(i4) = rv_chem(i4) + w4*tval(n)
             end do
             nullify(rv_chem)
          end do
       endif ! < l_do_adjoint >

       !aeroptr => aeroptr%llpoint
       aeroptr => aeroNode_nextcast(aeroptr)
!       call stop2(999)
    end do

    return
  end subroutine intaod_

end module intaodmod


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYALEM
C...Calculates the running alpha_electromagnetic.

      DOUBLE PRECISION FUNCTION PYALEM(Q2)

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)

C...Commonblocks.
      include 'inc/pydat1'

C...Calculate real part of photon vacuum polarization.
C...For leptons simplify by using asymptotic (Q^2 >> m^2) expressions.
C...For hadrons use parametrization of H. Burkhardt et al.
C...See R. Kleiss et al, CERN 89-08, vol. 3, pp. 129-131.
      AEMPI=PARU(101)/(3D0*PARU(1))
      IF(MSTU(101).LE.0.OR.Q2.LT.2D-6) THEN
        RPIGG=0D0
      ELSEIF(MSTU(101).EQ.2.AND.Q2.LT.PARU(104)) THEN
        RPIGG=0D0
      ELSEIF(MSTU(101).EQ.2) THEN
        RPIGG=1D0-PARU(101)/PARU(103)
      ELSEIF(Q2.LT.0.09D0) THEN
        RPIGG=AEMPI*(13.4916D0+LOG(Q2))+0.00835D0*LOG(1D0+Q2)
      ELSEIF(Q2.LT.9D0) THEN
        RPIGG=AEMPI*(16.3200D0+2D0*LOG(Q2))+
     &  0.00238D0*LOG(1D0+3.927D0*Q2)
      ELSEIF(Q2.LT.1D4) THEN
        RPIGG=AEMPI*(13.4955D0+3D0*LOG(Q2))+0.00165D0+
     &  0.00299D0*LOG(1D0+Q2)
      ELSE
        RPIGG=AEMPI*(13.4955D0+3D0*LOG(Q2))+0.00221D0+
     &  0.00293D0*LOG(1D0+Q2)
      ENDIF

C...Calculate running alpha_em.
      PYALEM=PARU(101)/(1D0-RPIGG)
      PARU(108)=PYALEM

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      subroutine flush(un)
      integer un
      end

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


  subroutine exchng1(a, nx, s, e, comm1d, nbrbottom, nbrtop)
  use mpi
  integer nx, s, e, comm1d, nbrbottom, nbrtop
  double precision a(0:nx+1,s-1:e+1)
  integer ierr, req(4)
!
  call MPI_IRECV (&
          a(1,s-1), nx, MPI_DOUBLE_PRECISION, nbrbottom, 0, &
          comm1d, req(1), ierr)
  call MPI_IRECV (&
          a(1,e+1), nx, MPI_DOUBLE_PRECISION, nbrtop, 1, &
          comm1d, req(2), ierr)
  call MPI_ISEND (&
          a(1,e), nx, MPI_DOUBLE_PRECISION, nbrtop, 0, &
          comm1d, req(3), ierr)
  call MPI_ISEND (&
          a(1,s), nx, MPI_DOUBLE_PRECISION, nbrbottom, 1, &
          comm1d, req(4), ierr)
!
  call MPI_WAITALL (4, req, MPI_STATUSES_IGNORE, ierr)
  return
  end


! { dg-do run }
! { dg-options "-fdump-tree-original" }
! Test that all string length calculations are
! optimized away.
program main
  character (len=999) :: c
  character (len=5) :: unit
  unit = '    '
  read (unit=unit,fmt='(I5)') i   ! Hide from optimizers
  j = 7
  c = '123456789'
  if (len(c(         3        : 5                )) /= 3) call abort ! Case  1
  if (len(c(     i*(i+1)      :     (i+1)*i + 2  )) /= 3) call abort ! Case  2
  if (len(c(     i*(i+1)      : 2 + (i+1)*i      )) /= 3) call abort ! Case  3
  if (len(c(     i*(i+1) + 2  :     (i+1)*i + 3  )) /= 2) call abort ! Case  4
  if (len(c( 2 + i*(i+1)      :     (i+1)*i + 3  )) /= 2) call abort ! Case  5
  if (len(c(     i*(i+1) + 2  : 3 + (i+1)*i      )) /= 2) call abort ! Case  6
  if (len(c( 2 + i*(i+1)      : 3 + (i+1)*i      )) /= 2) call abort ! Case  7
  if (len(c(     i*(i+1) - 1  :     (i+1)*i + 1  )) /= 3) call abort ! Case  8
  if (len(c(     i*(i+1) - 1  : 1 + (i+1)*i      )) /= 3) call abort ! Case  9
  if (len(c(     i*(i+1)      :     (i+1)*i -(-1))) /= 2) call abort ! Case 10
  if (len(c(     i*(i+1) +(-2):     (i+1)*i - 1  )) /= 2) call abort ! Case 11
  if (len(c(     i*(i+1) + 2  :     (i+1)*i -(-4))) /= 3) call abort ! Case 12
  if (len(c(     i*(i+1) - 3  :     (i+1)*i - 1  )) /= 3) call abort ! Case 13
  if (len(c(13 - i*(i+1)      :15 - (i+1)*i      )) /= 3) call abort ! Case 14
  if (len(c(     i*(i+1) +(-1):     (i+1)*i      )) /= 2) call abort ! Case 15
  if (len(c(-1 + i*(i+1)      :     (i+1)*i      )) /= 2) call abort ! Case 16
  if (len(c(     i*(i+1) - 2  :     (i+1)*i      )) /= 3) call abort ! Case 17
  if (len(c(  (i-2)*(i-3)     :   (i-3)*(i-2)    )) /= 1)       call abort ! Case 18
end program main
! { dg-final { scan-tree-dump-times "_abort" 0 "original" } }


text:
      Program DRIVER
c**************************************************************************
c     This is the driver for the whole calculation
c**************************************************************************
      use extra_weights
      use mint_module
      use FKSParams
      implicit none
C
C     CONSTANTS
C
      double precision zero
      parameter       (ZERO = 0d0)
      include 'nexternal.inc'
      include 'genps.inc'
C
C     LOCAL
C
      integer i,j,k,l,l1,l2,kchan
      character*130 buf
c
c     Global
c
cc
      include 'run.inc'
      include 'coupl.inc'

c     Vegas stuff
      integer         nndim
      common/tosigint/nndim

      real*8 sigint
      external sigint

      integer irestart
      logical savegrid

      logical            flat_grid
      common/to_readgrid/flat_grid                !Tells if grid read from file

      integer i_momcmp_count
      double precision xratmax
      common/ccheckcnt/i_momcmp_count,xratmax

      character*4      abrv
      common /to_abrv/ abrv
      integer n_mp, n_disc
      integer ixi_i,iphi_i,iy_ij,vn
      logical putonshell
      logical unwgt
      double precision evtsgn
      common /c_unwgt/evtsgn,unwgt

      logical SHsep
      logical Hevents
      common/SHevents/Hevents
      character*10 dum
c statistics for MadLoop
      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1(0:9)
      common/ups_stats/ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1

      double precision virtual_over_born
      common/c_vob/virtual_over_born
      include 'orders.inc'

c timing statistics
      include "timing_variables.inc"
      real*4 tOther, tTot

c PineAPPL
      logical pineappl
      common /for_pineappl/ pineappl
c stats for granny_is_res
      double precision deravg,derstd,dermax,xi_i_fks_ev_der_max
     &     ,y_ij_fks_ev_der_max
      integer ntot_granny,derntot,ncase(0:6)
      common /c_granny_counters/ ntot_granny,ncase,derntot,deravg,derstd
     &     ,dermax,xi_i_fks_ev_der_max,y_ij_fks_ev_der_max

      logical useitmax
      common/cuseitmax/useitmax

C-----
C  BEGIN CODE
C-----
c Write the process PID in the log.txt files (i.e., to the screen)
      write (*,*) getpid()

      useitmax=.false. ! to be overwritten in open_output_files.f if need be
c
c     Setup the timing variable
c
      call cpu_time(tBefore)
      fixed_order=.true.
      nlo_ps=.false.

c     Read general MadFKS parameters
c
      call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
      min_virt_fraction_mint=min_virt_fraction
      do kchan=1,maxchannels
         do i=0,n_ave_virt
            average_virtual(i,kchan)=0d0
         enddo
         virtual_fraction(kchan)=max(virt_fraction,min_virt_fraction)
      enddo
      n_ord_virt=amp_split_size

c
c     Read process number
c
      ntot_granny=0
      derntot=0
      do i=0,6
         ncase(i)=0
      enddo
      ntot=0
      nsun=0
      nsps=0
      nups=0
      neps=0
      n100=0
      nddp=0
      nqdp=0
      nini=0
      n10=0
      do i=0,9
        n1(i)=0
      enddo

      call setrun                !Sets up run parameters
      call setpara('param_card.dat')   !Sets up couplings and masses
      call setcuts               !Sets up cuts and particle masses
      call printout              !Prints out a summary of paramaters
      call run_printout          !Prints out a summary of the run settings
      call initcluster
      call check_amp_split
c
c     Get user input
c
      write(*,*) "getting user params"
      call get_user_params(ncalls0,itmax,imode)
      if(imode.eq.0)then
        flat_grid=.true.
      else
        flat_grid=.false.
      endif
      ndim = 3*(nexternal-nincoming)-4
      if (abs(lpp(1)) .ge. 1) ndim=ndim+1
      if (abs(lpp(2)) .ge. 1) ndim=ndim+1
      nndim=ndim
c Don't proceed if muF1#muF2 (we need to work out the relevant formulae
c at the NLO)
      if( ( fixed_fac_scale .and.
     #       (muF1_over_ref*muF1_ref_fixed) .ne.
     #       (muF2_over_ref*muF2_ref_fixed) ) .or.
     #    ( (.not.fixed_fac_scale) .and.
     #      muF1_over_ref.ne.muF2_over_ref ) )then
        write(*,*)'NLO computations require muF1=muF2'
        stop
      endif
      write(*,*) "about to integrate ", ndim,ncalls0,itmax
c PineAPPL
      if (imode.eq.0) pineappl=.False. ! overwrite when starting completely fresh
      if(pineappl) then
         write(6,*) "Initializing PineAPPL ..."
c     Set flavor map, starting from all possible
c     parton lumi configurations defined in initial_states_map.dat
         call setup_flavourmap
c     Fill the number of combined matrix elements for given initial state luminosity
         call find_iproc_map
         write(6,*) "   ... done."
      endif
      if (abrv(1:4).eq.'virt') then
         only_virt=.true.
      else
         only_virt=.false.
      endif
c     Prepare the MINT folding
      do j=1,ndimmax
         if (j.le.ndim) then
            ifold(j)=1
         else
            ifold(j)=0
         endif
      enddo
      ifold_energy=ndim-2
      ifold_yij=ndim-1
      ifold_phi=ndim
c
      i_momcmp_count=0
      xratmax=0.d0
      unwgt=.false.
      call addfil(dum)
      if (imode.eq.-1.or.imode.eq.0) then
         if(imode.eq.0)then
c Don't safe the reweight information when just setting up the grids.
            doreweight=.false.
            do_rwgt_scale=.false.
            do_rwgt_pdf=.false.
         else
            doreweight=do_rwgt_scale.or.do_rwgt_pdf.or.store_rwgt_info
         endif
c
         write (*,*) 'imode is ',imode

         if (ickkw.eq.-1) then
            min_virt_fraction=1d0
            do kchan=1,nchans
               virtual_fraction(kchan)=1d0
            enddo
         endif
         call mint(sigint)
         call topout
         call deallocate_weight_lines
c
      else
         write (*,*) 'Unknown imode',imode
         stop
      endif

      if (ntot.ne.0) then
         write(*,*) "Satistics from MadLoop:"
         write(*,*)
     &        "  Total points tried:                              ",ntot
         write(*,*)
     &        "  Stability unknown:                               ",nsun
         write(*,*)
     &        "  Stable PS point:                                 ",nsps
         write(*,*)
     &        "  Unstable PS point (and rescued):                 ",nups
         write(*,*)
     &        "  Exceptional PS point (unstable and not rescued): ",neps
         write(*,*)
     &        "  Double precision used:                           ",nddp
         write(*,*)
     &        "  Quadruple precision used:                        ",nqdp
         write(*,*)
     &        "  Initialization phase-space points:               ",nini
         write(*,*)
     &        "  Unknown return code (100):                       ",n100
         write(*,*)
     &        "  Unknown return code (10):                        ",n10
         write(*,*)
     &        "  Unit return code distribution (1):               "
         do j=0,9
           if (n1(j).ne.0) then
              write(*,*) "#Unit ",j," = ",n1(j)
           endif
         enddo
      endif

      write (*,*) 'counters for the granny resonances'
      write (*,*) 'ntot     ',ntot_granny
      if (ntot_granny.gt.0) then
         do i=0,6
            write (*,*) '% icase ',i,' : ',ncase(i)/dble(ntot_granny)
         enddo
         write (*,*) 'average,std dev. and max of derivative:',deravg
     &        ,sqrt(abs(derstd-deravg**2)),dermax
         write (*,*)
     &        'and xi_i_fks and y_ij_fks corresponding to max of der.',
     &        xi_i_fks_ev_der_max,y_ij_fks_ev_der_max
      endif
      call cpu_time(tAfter)
      tTot = tAfter-tBefore
      tOther = tTot - (tBorn+tGenPS+tReal+tCount+tIS+tFxFx+tf_nb+tf_all
     &     +t_as+tr_s+tr_pdf+t_plot+t_cuts+t_MC_subt+t_isum+t_p_unw
     $     +t_write)
      write(*,*) 'Time spent in Born : ',tBorn
      write(*,*) 'Time spent in PS_Generation : ',tGenPS
      write(*,*) 'Time spent in Reals_evaluation: ',tReal
      write(*,*) 'Time spent in MCsubtraction : ',t_MC_subt
      write(*,*) 'Time spent in Counter_terms : ',tCount
      write(*,*) 'Time spent in Integrated_CT : ',tIS-tOLP
      write(*,*) 'Time spent in Virtuals : ',tOLP
      write(*,*) 'Time spent in FxFx_cluster : ',tFxFx
      write(*,*) 'Time spent in Nbody_prefactor : ',tf_nb
      write(*,*) 'Time spent in N1body_prefactor : ',tf_all
      write(*,*) 'Time spent in Adding_alphas_pdf : ',t_as
      write(*,*) 'Time spent in Reweight_scale : ',tr_s
      write(*,*) 'Time spent in Reweight_pdf : ',tr_pdf
      write(*,*) 'Time spent in Filling_plots : ',t_plot
      write(*,*) 'Time spent in Applying_cuts : ',t_cuts
      write(*,*) 'Time spent in Sum_ident_contr : ',t_isum
      write(*,*) 'Time spent in Pick_unwgt : ',t_p_unw
      write(*,*) 'Time spent in Write_events : ',t_write
      write(*,*) 'Time spent in Other_tasks : ',tOther
      write(*,*) 'Time spent in Total : ',tTot

      open (unit=12, file='res.dat',status='unknown')
      do kchan=0,nchans
         write (12,*)ans(1,kchan),unc(1,kchan),ans(2,kchan),unc(2,kchan)
     $        ,itmax,ncalls0,tTot
      enddo
      close(12)

      if(i_momcmp_count.ne.0)then
        write(*,*)'     '
        write(*,*)'WARNING: genps_fks code 555555'
        write(*,*)i_momcmp_count,xratmax
      endif

      end

      block data timing
c timing statistics
      include "timing_variables.inc"
      data tOLP/0.0/
      data tGenPS/0.0/
      data tBorn/0.0/
      data tIS/0.0/
      data tReal/0.0/
      data tCount/0.0/
      data tFxFx/0.0/
      data tf_nb/0.0/
      data tf_all/0.0/
      data t_as/0.0/
      data tr_s/0.0/
      data tr_pdf/0.0/
      data t_plot/0.0/
      data t_cuts/0.0/
      data t_MC_subt/0.0/
      data t_isum/0.0/
      data t_p_unw/0.0/
      data t_write/0.0/
      end

      double precision function sigint(xx,vegas_wgt,ifl,f)
      use weight_lines
      use extra_weights
      use mint_module
      implicit none
      include 'nexternal.inc'
      include 'nFKSconfigs.inc'
      include 'run.inc'
      include 'orders.inc'
      include 'fks_info.inc'
      double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
     $     ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
      integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min,iamp
     $     ,nFKS_max,izero,ione,itwo,mohdr,i,iran_picked
      parameter (izero=0,ione=1,itwo=2,mohdr=-100)
      logical passcuts,passcuts_nbody,passcuts_n1body,sum,firsttime
      data firsttime/.true./
      external passcuts
      integer             ini_fin_fks(maxchannels)
      common/fks_channels/ini_fin_fks
      data sum /.false./
      integer         nndim
      common/tosigint/nndim
      logical       nbody
      common/cnbody/nbody
      double precision p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
     $     ,pswgt_cnt(-2:2),jac_cnt(-2:2)
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
      double precision p_born(0:3,nexternal-1)
      common /pborn/   p_born
      double precision virtual_over_born
      common/c_vob/virtual_over_born
      logical                calculatedBorn
      common/ccalculatedBorn/calculatedBorn
      character*4      abrv
      common /to_abrv/ abrv
c PineAPPL
      logical pineappl
      common /for_pineappl/ pineappl
      double precision       wgt_ME_born,wgt_ME_real
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
      integer ini_fin_fks_map(0:2,0:fks_configs)
      save ini_fin_fks_map
      if (new_point .and. ifl.ne.2) then
         pass_cuts_check=.false.
      endif
      if (firsttime) then
         firsttime=.false.
         call setup_ini_fin_fks_map(ini_fin_fks_map)
         write (*,*) 'initial-final FKS maps:'
         write (*,*) 0 ,':',ini_fin_fks_map(0,:)
         write (*,*) 1 ,':',ini_fin_fks_map(1,:)
         write (*,*) 2 ,':',ini_fin_fks_map(2,:)
      endif
      if (ifl.ne.0) then
         write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
         stop 1
      endif
      if (pineappl .and. sum) then
         write (*,*) 'WARNING: PineAPPL only possible '/
     &        /'with MC over FKS directories',pineappl,sum
         write (*,*) 'Switching to MC over FKS directories'
         sum=.false.
      endif
      sigint=0d0
      icontr=0
      do iamp=0,amp_split_size
         virt_wgt_mint(iamp)=0d0
         born_wgt_mint(iamp)=0d0
      enddo
      virtual_over_born=0d0
      wgt_me_born=0d0
      wgt_me_real=0d0
      if (ickkw.eq.-1) H1_factor_virt=0d0
      if (ickkw.eq.3) call set_FxFx_scale(0,p)
      call update_vegas_x(xx,x)
      call get_MC_integer(max(ini_fin_fks(ichan),1)
     $     ,ini_fin_fks_map(ini_fin_fks(ichan),0),iran_picked,vol)
      nFKS_picked=ini_fin_fks_map(ini_fin_fks(ichan),iran_picked)

c The nbody contributions
      if (abrv.eq.'real') goto 11
      nbody=.true.
      calculatedBorn=.false.
      call get_born_nFKSprocess(nFKS_picked,nFKS_born)
      call update_fks_dir(nFKS_born)
      if (ini_fin_fks(ichan).eq.0) then
         jac=1d0
      else
         jac=0.5d0
      endif
      call generate_momenta(nndim,iconfig,jac,x,p)
      if (p_born(0,1).lt.0d0) goto 12
      call compute_prefactors_nbody(vegas_wgt)
      call set_cms_stuff(izero)
      passcuts_nbody=passcuts(p1_cnt(0,1,0),rwgt)
      if (passcuts_nbody) then
         pass_cuts_check=.true.
         if (ickkw.eq.3) call set_FxFx_scale(1,p1_cnt(0,1,0))
         call set_alphaS(p1_cnt(0,1,0))
         call include_multichannel_enhance(1)
         if (abrv(1:2).ne.'vi') then
            call compute_born
         endif
         if (abrv.ne.'born') then
            call compute_nbody_noborn
         endif
      endif

 11   continue
c The n+1-body contributions (including counter terms)
      if ( abrv(1:4).eq.'born' .or.
     $     abrv(1:4).eq.'bovi' .or.
     $     abrv(1:2).eq.'vi' ) goto 12
      nbody=.false.
      if (sum) then
         nFKS_min=1
         nFKS_max=ini_fin_fks_map(ini_fin_fks(ichan),0)
         MC_int_wgt=1d0
      else
         nFKS_min=iran_picked
         nFKS_max=iran_picked
         MC_int_wgt=1d0/vol
      endif
      do i=nFKS_min,nFKS_max
         iFKS=ini_fin_fks_map(ini_fin_fks(ichan),i)
         calculatedBorn=.false.
         ! MZ this is a temporary fix for processes without
         ! soft singularities associated to the initial state
         ! DO NOT extend this fix to event generation
         wgt_me_born=0d0
         wgt_me_real=0d0
         jac=MC_int_wgt
         call update_fks_dir(iFKS)
         call generate_momenta(nndim,iconfig,jac,x,p)
         if (p_born(0,1).lt.0d0) cycle
         call compute_prefactors_n1body(vegas_wgt,jac)
         call set_cms_stuff(izero)
         passcuts_nbody =passcuts(p1_cnt(0,1,0),rwgt)
         call set_cms_stuff(mohdr)
         passcuts_n1body=passcuts(p,rwgt)
         if (passcuts_nbody .and. abrv.ne.'real') then
            pass_cuts_check=.true.
            call set_cms_stuff(izero)
            if (ickkw.eq.3) call set_FxFx_scale(2,p1_cnt(0,1,0))
            call set_alphaS(p1_cnt(0,1,0))
            call include_multichannel_enhance(3)
            call compute_soft_counter_term(0d0)
            call set_cms_stuff(ione)
            call compute_collinear_counter_term(0d0)
            call set_cms_stuff(itwo)
            call compute_soft_collinear_counter_term(0d0)
         endif
         if (passcuts_n1body) then
            pass_cuts_check=.true.
            call set_cms_stuff(mohdr)
            if (ickkw.eq.3) call set_FxFx_scale(3,p)
            call set_alphaS(p)
            call include_multichannel_enhance(2)
            call compute_real_emission(p,1d0)
         endif
      enddo

 12   continue
c Include PDFs and alpha_S and reweight to include the uncertainties
      if (ickkw.eq.-1) call include_veto_multiplier
      call include_PDF_and_alphas
      if (doreweight) then
         if (do_rwgt_scale .and. ickkw.ne.-1) call reweight_scale
         if (do_rwgt_scale .and. ickkw.eq.-1) call reweight_scale_NNLL
         if (do_rwgt_pdf) call reweight_pdf
      endif

      if (pineappl) then
         if (sum) then
            write (*,*) 'ERROR: PineAPPL only possible '/
     &           /'with MC over FKS directories',pineappl,sum
            stop 1
         endif
         call fill_pineappl_weights(vegas_wgt)
      endif

c Importance sampling for FKS configurations
      if (sum) then
         call get_wgt_nbody(sig)
         call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
     $        ,abs(sig))
      else
         call get_wgt_no_nbody(sig)
         call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
     $        ,abs(sig)*vol)
      endif

c Finalize PS point
      call fill_plots
      call fill_mint_function(f)
      return
      end

      subroutine update_fks_dir(nFKS)
      implicit none
      include 'run.inc'
      integer nFKS
      integer              nFKSprocess
      common/c_nFKSprocess/nFKSprocess
      nFKSprocess=nFKS
      call fks_inc_chooser()
      call leshouche_inc_chooser()
      call setcuts
      call setfksfactor(.false.)
      if (ickkw.eq.3) call configs_and_props_inc_chooser()
      return
      end

      subroutine setup_ini_fin_FKS_map(ini_fin_FKS_map)
      implicit none
      include 'nexternal.inc'
      include 'nFKSconfigs.inc'
      include 'fks_info.inc'
      integer ini_fin_FKS_map(0:2,0:fks_configs),iFKS
      ini_fin_FKS_map(0,0)=0
      ini_fin_FKS_map(1,0)=0
      ini_fin_FKS_map(2,0)=0
      do iFKS=1,fks_configs
         ini_fin_FKS_map(0,0)=ini_fin_FKS_map(0,0)+1
         ini_fin_FKS_map(0,ini_fin_FKS_map(0,0))=iFKS
         if (fks_j_d(iFKS).le.nincoming .and.
     $       fks_j_d(iFKS).gt.0) then
            ini_fin_FKS_map(2,0)=ini_fin_FKS_map(2,0)+1
            ini_fin_FKS_map(2,ini_fin_FKS_map(2,0))=iFKS
         elseif (fks_j_d(iFKS).gt.nincoming .and.
     $           fks_j_d(iFKS).le.nexternal) then
            ini_fin_FKS_map(1,0)=ini_fin_FKS_map(1,0)+1
            ini_fin_FKS_map(1,ini_fin_FKS_map(1,0))=iFKS
         else
            write (*,*) 'ERROR in setup_ini_fin_FKS_map',fks_j_d(iFKS)
     $           ,nincoming,iFKS
            stop 1
         endif
      enddo
      return
      end

      subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)
      implicit none
      include 'nexternal.inc'
      include 'nFKSconfigs.inc'
      include 'fks_info.inc'
      integer nFKS_in,nFKS_out,iFKS,iiFKS,nFKSprocessBorn(fks_configs)
      logical firsttime
      data firsttime /.true./
      save nFKSprocessBorn
c
      if (firsttime) then
         firsttime=.false.
         do iFKS=1,fks_configs
            nFKSprocessBorn(iFKS)=0
            if ( need_color_links_D(iFKS) .or.
     &           need_charge_links_D(iFKS) )then
               nFKSprocessBorn(iFKS)=iFKS
            endif
            if (nFKSprocessBorn(iFKS).eq.0) then
c     try to find the process that has the same j_fks but with i_fks a
c     gluon
               do iiFKS=1,fks_configs
                  if ( (need_color_links_D(iiFKS) .or.
     &                  need_charge_links_D(iiFKS)) .and.
     &                 fks_j_D(iFKS).eq.fks_j_D(iiFKS) ) then
                     nFKSprocessBorn(iFKS)=iiFKS
                     exit
                  endif
               enddo
            endif
c     try to find the process that has the j_fks initial state if
c     current j_fks is initial state (and similar for final state j_fks)
            if (nFKSprocessBorn(iFKS).eq.0) then
               do iiFKS=1,fks_configs
                  if ( need_color_links_D(iiFKS) .or.
     &                 need_charge_links_D(iiFKS) ) then
                     if ( fks_j_D(iiFKS).le.nincoming .and.
     &                    fks_j_D(iFKS).le.nincoming ) then
                        nFKSprocessBorn(iFKS)=iiFKS
                        exit
                     elseif ( fks_j_D(iiFKS).gt.nincoming .and.
     &                        fks_j_D(iFKS).gt.nincoming ) then
                        nFKSprocessBorn(iFKS)=iiFKS
                        exit
                     endif
                  endif
               enddo
            endif
c     If still not found, just pick any one that has a soft singularity
            if (nFKSprocessBorn(iFKS).eq.0) then
               do iiFKS=1,fks_configs
                  if ( need_color_links_D(iiFKS) .or.
     &                 need_charge_links_D(iiFKS) ) then
                     nFKSprocessBorn(iFKS)=iiFKS
                  endif
               enddo
            endif
c     if there are no soft singularities at all, just do something trivial
            if (nFKSprocessBorn(iFKS).eq.0) then
               nFKSprocessBorn(iFKS)=iFKS
            endif
         enddo
         write (*,*) 'Total number of FKS directories is', fks_configs
         write (*,*) 'For the Born we use nFKSprocesses:'
         write (*,*)  nFKSprocessBorn
      endif
      if (nFKSprocessBorn(nFKS_in).eq.0) then
         write(*,*) 'Could not find the correct map to Born '/
     &        /'FKS configuration for the NLO FKS '/
     &        /'configuration', nFKS_in
         stop 1
      else
         nFKS_out=nFKSprocessBorn(nFKS_in)
      endif
      return
      end

      subroutine update_vegas_x(xx,x)
      use mint_module
      implicit none
      integer i
      double precision xx(ndimmax),x(99),ran2
      external ran2
      integer         nndim
      common/tosigint/nndim
      character*4      abrv
      common /to_abrv/ abrv
      do i=1,99
         if (abrv.eq.'born'.or.abrv(1:2).eq.'vi') then
            if(i.le.nndim-3)then
               x(i)=xx(i)
            elseif(i.le.nndim) then
               x(i)=ran2()      ! Choose them flat when not including real-emision
            else
               x(i)=0.d0
            endif
         else
            if(i.le.nndim)then
               x(i)=xx(i)
            else
               x(i)=0.d0
            endif
         endif
      enddo
      return
      end

c
      subroutine get_user_params(ncall,nitmax,irestart)
c**********************************************************************
c     Routine to get user specified parameters for run
c**********************************************************************
      use mint_module
      implicit none
c
c     Constants
c
      include 'genps.inc'
      include 'nexternal.inc'
      include 'nFKSconfigs.inc'
      include 'fks_info.inc'
      include 'run.inc'
      include 'orders.inc'
c
c     Arguments
c
      integer ncall,nitmax
c
c     Local
c
      integer i, j, kchan
      double precision dconfig(maxchannels)
c
c     Global
c
      integer ini_fin_fks(maxchannels)
      common /fks_channels/ini_fin_fks
      integer           isum_hel
      logical                   multi_channel
      common/to_matrix/isum_hel, multi_channel
      integer           use_cut
      common /to_weight/use_cut

      integer        lbw(0:nexternal)  !Use of B.W.
      common /to_BW/ lbw

      character*5 abrvinput
      character*4 abrv
      common /to_abrv/ abrv

      logical nbody
      common/cnbody/nbody

c
c To convert diagram number to configuration
c
      include 'born_conf.inc'
c
c Vegas stuff
c
      integer irestart,itmp
      character * 70 idstring
      logical savegrid

      logical unwgt
      double precision evtsgn
      common /c_unwgt/evtsgn,unwgt
      logical fillh
      integer mc_hel,ihel
      double precision volh
      common/mc_int2/volh,mc_hel,ihel,fillh
      integer random_offset_split
      common /c_random_offset_split/ random_offset_split
      logical done
      character*140 buffer
c-----
c  Begin Code
c-----
      unwgt=.false.
      open (unit=83,file='input_app.txt',status='old')
      done=.false.
      nchans=0
      do while (.not. done)
         read(83,'(a)',err=222,end=222) buffer
         if (buffer(1:7).eq.'NPOINTS') then
            buffer=buffer(10:100)
            read(buffer,*) ncall
            write (*,*) 'Number of phase-space points per iteration:',ncall
         elseif(buffer(1:11).eq.'NITERATIONS') then
            read(buffer(14:),*) nitmax
            write (*,*) 'Maximum number of iterations is:',nitmax
         elseif(buffer(1:8).eq.'ACCURACY') then
            read(buffer(11:),*) accuracy
            write (*,*) 'Desired accuracy is:',accuracy
         elseif(buffer(1:10).eq.'ADAPT_GRID') then
            read(buffer(13:),*) use_cut
            write (*,*) 'Using adaptive grids:',use_cut
         elseif(buffer(1:12).eq.'MULTICHANNEL') then
            read(buffer(15:),*) i
            if (i.eq.1) then
               multi_channel=.true.
               write (*,*) 'Using Multi-channel integration'
            else
               multi_channel=.false.
               write (*,*) 'Not using Multi-channel integration'
            endif
         elseif(buffer(1:12).eq.'SUM_HELICITY') then
            read(buffer(15:),*) i
            if (nincoming.eq.1) then
               write (*,*) 'Sum over helicities in the virtuals'/
     $              /' for decay process'
               mc_hel=0
            elseif (i.eq.0) then
               mc_hel=0
               write (*,*) 'Explicitly summing over helicities'/
     $              /' for the virtuals'
            else
               mc_hel=1
               write(*,*) 'Do MC over helicities for the virtuals'
            endif
            isum_hel=0
         elseif(buffer(1:6).eq.'NCHANS') then
            read(buffer(9:),*) nchans
            write (*,*) 'Number of channels to integrate together:'
     $           ,nchans
            if (nchans.gt.maxchannels) then
               write (*,*) 'Too many integration channels to be '/
     $              /'integrated together. Increase maxchannels',nchans
     $              ,maxchannels
               stop 1
            endif
         elseif(buffer(1:7).eq.'CHANNEL') then
            if (nchans.le.0) then
               write (*,*) '"NCHANS" missing in input files'/
     $              /' (still zero)',nchans
               stop
            endif
            read(buffer(10:),*) (dconfig(kchan),kchan=1,nchans)
            do kchan=1,nchans
               iconfigs(kchan) = int(dconfig(kchan))
               if ( nint(dconfig(kchan)*10)-iconfigs(kchan)*10.eq.0 )
     $              then
                  ini_fin_fks(kchan)=0
               elseif ( nint(dconfig(kchan)*10)-iconfigs(kchan)*10.eq.1
     $                 ) then
                  ini_fin_fks(kchan)=1
               elseif ( nint(dconfig(kchan)*10)-iconfigs(kchan)*10.eq.2
     $                 ) then
                  ini_fin_fks(kchan)=2
               else
                  write (*,*) 'ERROR: invalid configuration number',dconfig
                  stop 1
               endif
               do i=1,mapconfig(0)
                  if (iconfigs(kchan).eq.mapconfig(i)) then
                     iconfigs(kchan)=i
                     exit
                  endif
               enddo
            enddo
            write(*,*) 'Running Configuration Number(s): '
     $           ,(iconfigs(kchan),kchan=1,nchans)
            write(*,*) 'initial-or-final',(ini_fin_fks(kchan),kchan=1,nchans)
         elseif(buffer(1:5).eq.'SPLIT') then
            read(buffer(8:),*) random_offset_split
            write (*,*) 'Splitting channel:',random_offset_split
         elseif(buffer(1:8).eq.'WGT_MULT') then
            read(buffer(11:),*) wgt_mult
            write (*,*) 'Weight multiplier:',wgt_mult
         elseif(buffer(1:8).eq.'RUN_MODE') then
            read(buffer(11:),*) abrvinput
            if(abrvinput(5:5).eq.'0')then
               nbody=.true.
            else
               nbody=.false.
            endif
            abrv=abrvinput(1:4)
            write (*,*) "doing the ",abrv," of this channel"
            if(nbody)then
               write (*,*) "integration Born/virtual with Sfunction=1"
            else
               write (*,*) "Normal integration (Sfunction != 1)"
            endif
         elseif(buffer(1:7).eq.'RESTART') then
            read(buffer(10:),*) irestart
            if (irestart.eq.0) then
               write (*,*) 'RESTART: Fresh run'
            elseif(irestart.eq.-1) then
               write (*,*) 'RESTART: Use old grids, but refil plots'
            elseif(irestart.eq.1) then
               write (*,*) 'RESTART: continue with existing run'
            else
               write (*,*) 'RESTART:',irestart
            endif
         endif
         cycle
 222     done=.true.
      enddo
      close(83)

      if (fks_configs.eq.1) then
         if (pdg_type_d(1,fks_i_d(1)).eq.-21.and.abrv.ne.'born') then
C Two cases can occur
C   1) the process has been generated with the LOonly flav
C   2) the process has only virtual corrections, e.g. z > v v [QED]
C the two cases can be distinguished by looking at the values
C  of AMP_SPLIT_SIZE, AMP_SPLIT_SIZE_BORN (if they are ==, it is 1))
           if (amp_split_size.eq.amp_split_size_born) then
             write (*,*) 'Process generated with [LOonly=QCD]. '/
     $           /'Setting abrv to "born".'
             abrv='born'
             if (ickkw.eq.3) then
               write (*,*) 'FxFx merging not possible with'/
     $              /' [LOonly=QCD] processes'
               stop 1
             endif
           else
             write (*,*) 'Process only with virtual corrections'/
     $           /'Setting abrv to "bovi".'
             abrv='bovi'
           endif
         endif
      endif
c
c     Here I want to set up with B.W. we map and which we don't
c
      lbw(0)=0
      end
c


submodule(euler_interface_m) euler_prob_0058_m
    implicit none

contains

    module character(len=20) function euler0058()
        write (euler0058, "(i20)") ans()
    end function euler0058

    integer(int64) function ans()
        use euler_primes_m, only: is_prime
        implicit none

        integer(int64) :: n, k, i, arr(4)

        k = 0
        n = 2
        main_loop: do
            call diagonal_nums(n, arr)
            do i = 1, 3
                if (is_prime(arr(i))) k = k + 1
            end do

            associate (p_nums => real(k), tot_nums => real((n - 1)*4 + 1))
                if (p_nums/tot_nums < 0.1_sp) exit main_loop
            end associate
            n = n + 1
        end do main_loop

        ans = side_len(n)
    end function ans

    subroutine diagonal_nums(n, arr)
        integer(int64), intent(in) :: n
        integer(int64), intent(out) :: arr(:)
        integer(int64) :: bottom_left, diff

        bottom_left = (n*2 - 1)**2
        diff = (n - 1)*2
        arr(1) = bottom_left - 3*diff
        arr(2) = bottom_left - 2*diff
        arr(3) = bottom_left - diff
        arr(4) = bottom_left
    end subroutine diagonal_nums

    integer(int64) function side_len(n)
        integer(int64), intent(in) :: n

        side_len = n*2_int64 - 1_int64
    end function side_len

end submodule euler_prob_0058_m


text:
!===============================================================================
!  S H E B A - Shear-wave Birefringence Analysis
!===============================================================================
!  This software is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!===============================================================================
!
!  James Wookey, School of Earth Sciences, University of Bristol
!
!=======================================================================
      module array_sizes ! Utility module for F90/95 for SAC files
!=======================================================================
      implicit none

!        * parameters
         integer, parameter :: np1=181
         integer, parameter :: np2=41
         integer, parameter :: np2int=161
         integer, parameter :: np2XC=np2
         integer, parameter :: np2XCint=np2int
         integer, parameter :: npc=500
         integer, parameter :: np=200000
         integer, parameter :: nsurfmax = 1000 ; ! this one is for stack

!=======================================================================
      end module array_sizes
!=======================================================================
!     END OF sheba_config module
!=======================================================================

max_stars_repo_name: jwookey/sheba

max_issues_repo_name: jwookey/sheba

max_forks_repo_name: jwookey/sheba