      SUBROUTINE ESC96(R,TYPE,EFAC,FI0,FI1,FI2,VPOT)
************************************************************************
**    Nijmegen Extended Soft-Core (ESC) NN potential in coordinate space
**    including pi-meson BW, TMO, and pair potentials
**    Version 1.1: October 1999
**    Email: info@nn-online.org
**    WWW: http://nn-online.org
**    References: Rijken & Stoks, Phys. Rev. C 54 (1996), 2851, 2869.
**
**    INPUT:  R       in fermi
**    -----   TYPE    'PP', 'NP', 'PN', or 'NN' (character*2)
**            EFAC    M/E factor for energy dependent one-pion exchange
**
**    OUTPUT: VPOT    2x2 potential matrix in MeV on LSJ-basis
**    ------  FI0,FI1,FI2 2x2 nonlocal contribution and its derivatives
**                        with dimensions [0], [MeV], [MeV**2],
**                        including the tensor nonlocal contribution.
**
**    COMMON/EMANHP/PHNAME has to be filled beforehand.
**           PHNAME is character*3 and contains the name of the phase
**           shift in the spectral notation.
**           - singlets:           1S0  1P1  1D2  1F3  1G4 ...
**           - triplets uncoupled: 3P0  3P1  3D2  3F3  3G4 ...
**           - triplets coupled:        3C1  3C2  3C3  3C4 ...
**             where 3C1 denotes 3S1 - epsilon1 - 3D1 channel
**                   3C2 denotes 3P2 - epsilon2 - 3F2 channel ...
**    DATA statements:
**         IPAIR = 1: pair potentials included
**         IPIMES= 1: pi-meson potentials included
**         IAXIAL= 1: axial OBE potentials included
**         NLOC  = 1: standard nonlocalities
**                 2: nonlocalities in spin-spin and tensor of OPE
************************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER PHNAME*3,PHNAM0*3, TYPE*2,TYP0*2
      LOGICAL FIRST, CALC, NWGRID, niso0,niso1
      INTEGER SPIN
      REAL*8 VPOT(2,2), FI0(2,2),FI1(2,2),FI2(2,2)
      PARAMETER(IPTMAX=2000)
      COMMON/EMANHP/PHNAME
      COMMON/POT2/  UC,USS,UT,ULS,ULSA,UQ12,FIC0,FIC1,FIC2
      COMMON/VOPE/  VPIS,VPIT,VPIQ,FPIS0,FPIS1,FPIS2,FPIT0,FPIT1,FPIT2
      COMMON/VENER/ PLMEV2,VEPIS,VEPIT
      COMMON/VPIMES/VCPI,VSPI,VTPI,VLPI, VCPM,VSPM,VTPM,VLPM
      COMMON/VPAIR/ V1C,V1S,V1T,V1L, V2C,V2S,V2T,V2L
      COMMON/POTFSV/ICAL,IPT,POTMAT(IPTMAX,10),OPEMAT(IPTMAX,2)
      DATA NLOC/1/, NADIA/0/, IPAIR/1/, IPIMES/1/, IAXIAL/0/
      DATA NOWRIT/1/
      DATA R0/-1D0/, ISO0/-1/, PHNAM0/'***'/, TYP0/'XX'/
      data niso0/.true./, niso1/.true./, rw0/1d-3/, eps/1d-4/
      SAVE FIRST,CALC,NWGRID, NCHAN,SPIN,L,J,ISO, TJ,TJM,TJP,SQJJ

      DO 50 I=1,2*2
        VPOT(I,1)=0D0
        FI0(I,1)=0D0
        FI1(I,1)=0D0
        FI2(I,1)=0D0
   50 CONTINUE

      IF(PHNAME.NE.PHNAM0) THEN
        IF(NOWRIT.EQ.1) PHNAM0=PHNAME
        NCHAN=1
        IF(PHNAME(2:2).EQ.'C') NCHAN=2
        IF(PHNAME(1:1).EQ.'1') SPIN=0
        IF(PHNAME(1:1).EQ.'3') SPIN=1
        READ(PHNAME,'(2X,I1)') J
        L=J
        IF(PHNAME.EQ.'3P0') L=1
        IF(NCHAN.EQ.2) L=J-1
        ISO=MOD(SPIN+L+1,2)
        TJ=2D0*J+1D0
        TJM=(J-1D0)
        TJP=(J+2D0)
        SQJJ=DSQRT(J*(J+1D0))
      ENDIF

      IF(NOWRIT.EQ.1) THEN
***     Potential is not stored in POTMAT
        CALL NYMISO(R,TYPE,ISO,NLOC,NADIA)
        VC  =UC
        VSS =USS
        VT  =UT
        VLS =ULS
        VLSA=ULSA
        VQ12=UQ12
        IF(IPAIR.EQ.1) THEN
          CALL PAIR(R,TYPE,ISO)
          VC =VC + (V1C+V2C)
          VSS=VSS+ (V1S+V2S)
          VT =VT + (V1T+V2T)
          VLS=VLS+ (V1L+V2L)
        ENDIF
        IF(IPIMES.EQ.1) THEN
          CALL PIMES(R,TYPE,ISO)
          VC =VC + (VCPI+VCPM)
          VSS=VSS+ (VSPI+VSPM)
          VT =VT + (VTPI+VTPM)
          VLS=VLS+ (VLPI+VLPM)
        ENDIF
        IF(IAXIAL.EQ.1) THEN
          CALL AXIAL(R,ISO,SPIN,VAXSS,VAXT,VAXLS,FI0AX,FI1AX,FI2AX)
          VSS =VSS +VAXSS
          VT  =VT  +VAXT
          VLS =VLS +VAXLS
          FIC0=FIC0+FI0AX
          FIC1=FIC1+FI1AX
          FIC2=FIC2+FI2AX
        ENDIF
        GOTO 13
      ENDIF

      NWGRID=.FALSE.
      IF(PHNAME.NE.PHNAM0 .AND.
     .  (PHNAME.EQ.'1S0' .OR. PHNAM0.EQ.'1S0')) NWGRID=.TRUE.
      IF(ISO.NE.ISO0 .OR. TYPE.NE.TYP0 .OR. NWGRID) FIRST=.TRUE.
      IF(ICAL.EQ.0) FIRST=.TRUE.
      IF(FIRST) THEN
        CALC=.TRUE.
        IPT=1
      ENDIF

      IF(.NOT.FIRST .AND. R.LT.R0) THEN
        CALC=.FALSE.
        IPT=1
      ENDIF
      ICAL=1
      FIRST=.FALSE.
      ISO0=ISO
      TYP0=TYPE
      PHNAM0=PHNAME
      R0=R
      IF(.NOT.CALC) THEN
   11   R0=POTMAT(IPT,1)
        IF(R0.LT.R) THEN
          IPT=IPT+1
          IF(IPT.GT.IPTMAX) GOTO 12
          GOTO 11
        ELSEIF(DABS(R0-R).LT.1D-6) THEN
          VC  =POTMAT(IPT,2)
          VSS =POTMAT(IPT,3)
          VT  =POTMAT(IPT,4)
          VLS =POTMAT(IPT,5)
          VLSA=POTMAT(IPT,6)
          VQ12=POTMAT(IPT,7)
          FIC0=POTMAT(IPT,8)
          FIC1=POTMAT(IPT,9)
          FIC2=POTMAT(IPT,10)
c         FIS0=POTMAT(IPT,11)
c         FIS1=POTMAT(IPT,12)
c         FIS2=POTMAT(IPT,13)
c         FIT0=POTMAT(IPT,14)
c         FIT1=POTMAT(IPT,15)
c         FIT2=POTMAT(IPT,16)
c         IF(NADIA.EQ.1) VEPIS=POTMAT(IPT,17)
c         IF(NADIA.EQ.1) VEPIT=POTMAT(IPT,18)
          VPIS=OPEMAT(IPT,1)
          VPIT=OPEMAT(IPT,2)
          GOTO 13
        ENDIF
   12   STOP '*** NIJMESC96: stored grid file differs with input'
      ELSE
        CALL NYMISO(R,TYPE,ISO,NLOC,NADIA)
        VC  =UC
        VSS =USS
        VT  =UT
        VLS =ULS
        VLSA=ULSA
        VQ12=UQ12
        IF(IPAIR.EQ.1) THEN
          CALL PAIR(R,TYPE,ISO)
          VC =VC + (V1C+V2C)
          VSS=VSS+ (V1S+V2S)
          VT =VT + (V1T+V2T)
          VLS=VLS+ (V1L+V2L)
        ENDIF
        IF(IPIMES.EQ.1) THEN
          CALL PIMES(R,TYPE,ISO)
          VC =VC + (VCPI+VCPM)
          VSS=VSS+ (VSPI+VSPM)
          VT =VT + (VTPI+VTPM)
          VLS=VLS+ (VLPI+VLPM)
        ENDIF
        IF(IAXIAL.EQ.1) THEN
          CALL AXIAL(R,ISO,SPIN,VAXSS,VAXT,VAXLS,FI0AX,FI1AX,FI2AX)
          VSS =VSS +VAXSS
          VT  =VT  +VAXT
          VLS =VLS +VAXLS
          FIC0=FIC0+FI0AX
          FIC1=FIC1+FI1AX
          FIC2=FIC2+FI2AX
        ENDIF
        POTMAT(IPT,1) =R
        POTMAT(IPT,2) =VC
        POTMAT(IPT,3) =VSS
        POTMAT(IPT,4) =VT
        POTMAT(IPT,5) =VLS
        POTMAT(IPT,6) =VLSA
        POTMAT(IPT,7) =VQ12
        POTMAT(IPT,8) =FIC0
        POTMAT(IPT,9) =FIC1
        POTMAT(IPT,10)=FIC2
c       POTMAT(IPT,11)=FIS0
c       POTMAT(IPT,12)=FIS1
c       POTMAT(IPT,13)=FIS2
c       POTMAT(IPT,14)=FIT0
c       POTMAT(IPT,15)=FIT1
c       POTMAT(IPT,16)=FIT2
c       IF(NADIA.EQ.1) POTMAT(IPT,17)=VEPIS
c       IF(NADIA.EQ.1) POTMAT(IPT,18)=VEPIT
        OPEMAT(IPT,1)=VPIS
        OPEMAT(IPT,2)=VPIT
        IPT=IPT+1
        IF(IPT.GT.IPTMAX) STOP '*** NIJMESC96: POTMAT is too small'
      ENDIF

   13 CONTINUE

      IF(NADIA.EQ.1) THEN
***     Non-adiabatic contribution to OPE
        VSS=VSS+VPIS + PLMEV2*VEPIS
        VT =VT +VPIT + PLMEV2*VEPIT
      ELSEIF(NLOC.EQ.1) THEN
***     Energy-dependent (M/E-factor) one-poin exchange
c       CALL PIOME(R,TYPE,EFAC,VPIS,VPIT)
        VSS=VSS+VPIS*EFAC
        VT =VT +VPIT*EFAC
      ENDIF
      IF(NCHAN.EQ.1) THEN
        IF(SPIN.EQ.0) THEN
          VPOT(1,1) = VC - 3D0*VSS - J*(J+1)*VQ12
        ELSEIF(L.EQ.J) THEN
          VPOT(1,1) = VC + VSS + 2D0*VT - VLS + (1D0-J*(J+1))*VQ12
        ELSEIF(PHNAME.EQ.'3P0') THEN
          VPOT(1,1) = VC + VSS - 4D0*VT - 2D0*VLS + 4D0*VQ12
        ENDIF
      ELSE
        VPOT(1,1) = VC + VSS - 2D0*(J-1)/(2*J+1)*VT +
     .              (J-1)*VLS + (J-1)*(J-1)*VQ12
        VPOT(2,2) = VC + VSS - 2D0*(J+2)/(2*J+1)*VT -
     .              (J+2)*VLS + (J+2)*(J+2)*VQ12
        VPOT(1,2) = 6D0*DSQRT(J*(J+1)*1D0)/(2*J+1)*VT
        VPOT(2,1) = VPOT(1,2)
      ENDIF

      IF(NLOC.EQ.1 .AND. NADIA.EQ.0) THEN
        FI0(1,1) = FIC0
        FI1(1,1) = FIC1
        FI2(1,1) = FIC2
        IF(NCHAN.EQ.2) THEN
          FI0(2,2) = FIC0
          FI1(2,2) = FIC1
          FI2(2,2) = FIC2
        ENDIF
c     ELSEIF(NLOC.EQ.2 .OR. NADIA.EQ.1) THEN
c       IF(NCHAN.EQ.1) THEN
c         IF(SPIN.EQ.0) THEN
c           FI0(1,1) = FIC0 - 3D0*FIS0
c           FI1(1,1) = FIC1 - 3D0*FIS1
c           FI2(1,1) = FIC2 - 3D0*FIS2
c         ELSEIF(L.EQ.J) THEN
c           FI0(1,1) = FIC0 + FIS0 + 2D0*FIT0
c           FI1(1,1) = FIC1 + FIS1 + 2D0*FIT1
c           FI2(1,1) = FIC2 + FIS2 + 2D0*FIT2
c         ELSEIF(PHNAME.EQ.'3P0') THEN
c           FI0(1,1) = FIC0 + FIS0 - 4D0*FIT0
c           FI1(1,1) = FIC1 + FIS1 - 4D0*FIT1
c           FI2(1,1) = FIC2 + FIS2 - 4D0*FIT2
c         ENDIF
c       ELSEIF(NCHAN.EQ.2) THEN
c         FI0(1,1) = FIC0 + FIS0 - 2D0*TJM/TJ*FIT0
c         FI1(1,1) = FIC1 + FIS1 - 2D0*TJM/TJ*FIT1
c         FI2(1,1) = FIC2 + FIS2 - 2D0*TJM/TJ*FIT2
c         FI0(1,2) = 6D0*SQJJ/TJ*FIT0
c         FI1(1,2) = 6D0*SQJJ/TJ*FIT1
c         FI2(1,2) = 6D0*SQJJ/TJ*FIT2
c         FI0(2,1) = 6D0*SQJJ/TJ*FIT0
c         FI1(2,1) = 6D0*SQJJ/TJ*FIT1
c         FI2(2,1) = 6D0*SQJJ/TJ*FIT2
c         FI0(2,2) = FIC0 + FIS0 - 2D0*TJP/TJ*FIT0
c         FI1(2,2) = FIC1 + FIS1 - 2D0*TJP/TJ*FIT1
c         FI2(2,2) = FIC2 + FIS2 - 2D0*TJP/TJ*FIT2
c       ENDIF
      ENDIF

      if( ((niso0.and.iso.eq.0) .or. (niso1.and.iso.eq.1))
     .    .and. dabs(r-rw0).lt.eps) then
        eptf=1d0+2d0*fi0(1,1)
        v11=vc-3d0*vss-j*(j+1)*vq12
        v33=vc+vss+2d0*vt-vls+(1d0-j*(j+1))*vq12
        write(*,'(a6,2x,d15.8,2(a14,f12.5))') 'r(fm)=',r,'eff-sing=',
     .    v11/eptf-(fi1(1,1)/eptf)**2/939d0,'eff-trip=',
     .    v33/eptf-(fi1(1,1)/eptf)**2/939d0
        write(*,'(a4,i1,5x,5f13.5)') 'iso=',iso,
     .           vc,vss,vt,vls,vq12
        write(*,'(6x,a3,1x,5f13.5)') 'OBE',
     .           uc,uss+vpis*efac,ut+vpit*efac,uls,uq12+vpiq*efac
        write(*,99) ' pi-pi',vcpi,vspi,vtpi,vlpi
        write(*,99) 'pi-mes',vcpm,vspm,vtpm,vlpm
        write(*,99) '  pair',v1c+v2c,v1s+v2s,v1t+v2t,v1l+v2l
        write(*,99) '1-pair',v1c,v1s,v1t,v1l
        write(*,99) '2-pair',v2c,v2s,v2t,v2l
        if(iso.eq.0) niso0=.false.
        if(iso.eq.1) niso1=.false.
   99   format(3x,a6,1x,4f13.5)
      endif

      RETURN
      END
************************************************************************
**          NIJMEGEN  NUCLEON - NUCLEON  POTENTIAL  PROGRAM           **
**          -----------------------------------------------           **
************************************************************************
      SUBROUTINE NYMISO(R,TYPE,ISO,NLOC,NADIA)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 NEUTM
      CHARACTER PHNAME*3,PHNAM0*3, TYPE*2,TYP0*2
      LOGICAL FIRST, BRDEP,BRDRO,BRDROC
      common/hoiho/ scal
      COMMON/EMANHP/PHNAME
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/PARAMS/PAR(8,5)
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/KAPPAV/AKRO,AKOM,AKFI
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/BRDRHO/GAMRO,THRRO,NRO,BRDRO,GAMROC,THRROC,NROC,BRDROC
      COMMON/BRDEPS/GAMEP,THREP0,THREPC,NEP,BRDEP
      COMMON/SCALIN/AMT,AMPV,AMT24
      COMMON/AMSCAL/PIOMS2,PIOMS
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/CUTTEN/ALAMT1,ALAVT1,ALAMT0,ALAVT0
      COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4,
     .              GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4
      COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2,
     .              AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2,
     .              AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2
      COMMON/GPAIR/ GPIETA,GPIETP,GSISI, GPISIG,GPIPOM,
     .              GPIPI0,GPIPI1,FPIPI1,GPIRO0,GPIRO1,GPIOM
      COMMON/HBARC/ HBC
      DATA FPPPI0/0.0745D0/, FIRST/.TRUE./
      DATA TYP0/'XX'/, PHNAM0/'***'/, RGT/25D0/
      DATA PI/3.14159265358979323846D0/, SRPI/1.7724538509055160273D0/
      SAVE ISIGN,CONV,PROTM,NEUTM

      IF(.NOT.FIRST) GOTO 10
      FIRST = .FALSE.
      CONV  = PI/180D0
**    Nucleon and meson masses (Particle Data Group 1992)
**    Only include meson-mass difference for pions
      HBC   = 197.327053D0
      PROTM = 938.27231D0
      NEUTM = 939.56563D0
      AMT   = 938.27231D0
      AMT24 = 4D0*AMT*AMT
      AMPV  = 139.5679D0
      AMPIC = 139.5679D0
      AMPI  = 134.9743D0
      PIM = (AMPI+2D0*AMPIC)/3D0
      AMETA = 547.45D0
      AMETP = 957.75D0
      ROM   = 768.1D0
      AMROC = ROM
      AMRO  = ROM
      AMOM  = 781.95D0
      AMFI  =1019.413D0
      A0M   = 982.7D0
      AMA0C = A0M
      AMA0  = A0M
      AMEP  = 760D0
      AMF0  = 974.1D0
      AVSC  = 1D0/AMROC**2
      FAC   = 1D0
      AWPIC = AMPIC
      AWROC = AMROC
      AWA0C = AMA0C
      PIOMS = AMPV
      PIOMS2= AMPV*AMPV
**    Broad rho meson: spectral density to two effective Yukawa's
**    Yukawa's fitted to PHI0C 0.001 - 2 fm (steps 0.005, ALAM1=850 MeV)
      BRDRO = .FALSE.
      GAMRO = 151.5D0
      THRRO = PIM+PIM
      NRO   = 1
      ARO   = 0.35765D0
      AMR1  = 667.877D0
      BRO   = 0.47396D0
      AMR2  = 914.234D0
      BRDROC= BRDRO
      GAMROC= GAMRO
      THRROC= THRRO
      NROC  = NRO
      AROC  = ARO
      AMRC1 = AMR1
      AWRC1 = AMRC1
      BROC  = BRO
      AMRC2 = AMR2
      AWRC2 = AMRC2
      AVSC1 = 1D0/AMRC1**2
      AVSC2 = 1D0/AMRC2**2
**    Broad epsilon meson: spectral density to two effective Yukawa's
**    Yukawa's fitted to PHI0C 0.001 - 2 fm (steps 0.005, ALAM0=800 MeV)
      BRDEP = .FALSE.
      GAMEP = 800D0
      THREP0= AMPI+AMPI
      THREPC= AMPIC+AMPIC
      NEP   = 0
      AEPS  = 0.16914D0
      AME1  = 487.913D0
      BEPS  = 0.61302D0
      AME2  =1021.402D0

      CALL ESCPAR
C              pseudovec  vector    tensor   scalar    pomeron
C              FPI        GRO       AKRO     GA0       GA2D
C              PV1        GV1       FV1      GEP       GPOM
C              ALPV       ALVE      ALVM     GF0       AMPOM
C              THPV       THV       GA0C%    THGS      AMA2D
C     Cut-off  ALAMP1     ALAMV1    ALAM1    ALAMS1    SCAL
C              ALAMP0     ALAMV0    ALAM0    ALAMS0    GPIOM
C     Pairs    GPIETA     GPIETP    GSISI    GPISIG    GPIPOM
C              GPIPI0     GPIPI1    FPIPI1   GPIRO0    GPIRO1

**    (joined) cutoffs for pseudoscalar, vector, scalar
      ALAMP1 = PAR(5,1)
      ALAMV1 = PAR(5,2)
      ALAMS1 = PAR(5,4)
      ALAM1  = PAR(5,3)
      ALAMP0 = PAR(6,1)
      ALAMV0 = PAR(6,2)
      ALAMS0 = PAR(6,4)
      ALAM0  = PAR(6,3)
      IF(ALAM1.NE.0D0) THEN
        ALAMP1=ALAM1
        ALAMV1=ALAM1
        ALAMS1=ALAM1
        IF(ALAM0.EQ.0D0) ALAM0=ALAM1
      ENDIF
      IF(ALAM0.NE.0D0) THEN
        ALAMP0=ALAM0
        ALAMV0=ALAM0
        ALAMS0=ALAM0
      ENDIF
      alamt1=alamv1
      alavt1=alamv1
      alamt0=alamv0
      alavt0=alamv0
**    pseudovector couplings
c     FPI   = PAR(1,1)
      FPI   = DSQRT(FPPPI0*DEXP(-(AMPI/ALAMP1)**2))
      PAR(1,1)=FPI
      PV1   = PAR(2,1)
      ALPV  = PAR(3,1)
      THPV  = PAR(4,1)*CONV
**    vector couplings
      GRO   = PAR(1,2)
      GV1   = PAR(2,2)
      ALVE  = PAR(3,2)
      THV   = PAR(4,2)*CONV
**    tensor couplings
      AKRO  = PAR(1,3)
      FV1   = PAR(2,3)
      ALVM  = PAR(3,3)
**    scalar couplings
      GA0   = PAR(1,4)
      GEP   = PAR(2,4)
      GF0   = PAR(3,4)
      THGS  = PAR(4,4)*CONV
**    diffractive contribution
      GA2D  = PAR(1,5)
      GPOM  = PAR(2,5)
      AMPOM = PAR(3,5)
      AMA2D = PAR(4,5)
      IF(AMA2D.EQ.0D0) AMA2D=AMPOM
      AMPOM2= AMPOM*AMPOM
      AMPOM4= AMPOM2*AMPOM2
      AMA2D2= AMA2D*AMA2D
      AMA2D4= AMA2D2*AMA2D2
**    useful parameter
      scal  = PAR(5,5)
**    Pair couplings
      GPIETA = PAR(7,1)
      GPIETP = PAR(7,2)
      GSISI  = PAR(7,3)
      GPISIG = PAR(7,4)
      GPIPOM = PAR(7,5)
      GPIPI0 = PAR(8,1)
      GPIPI1 = PAR(8,2)
      FPIPI1 = PAR(8,3)
      GPIRO0 = PAR(8,4)
      GPIRO1 = PAR(8,5)
      GPIOM  = PAR(6,5)
**    define the physical coupling constants
      CALL NYMCPFF
      CALL CHIRAL

   10 CONTINUE
      IF(TYPE.EQ.TYP0) GOTO 20
      TYP0=TYPE
      IF(TYPE.EQ.'PP') THEN
        AMY = PROTM
        AMN = PROTM
        I3Y = 1
        I3N = 1
      ELSEIF(TYPE.EQ.'NP') THEN
        AMY = NEUTM
        AMN = PROTM
        I3Y =-1
        I3N = 1
      ELSEIF(TYPE.EQ.'PN') THEN
        AMY = PROTM
        AMN = NEUTM
        I3Y = 1
        I3N =-1
      ELSEIF(TYPE.EQ.'NN') THEN
        AMY = NEUTM
        AMN = NEUTM
        I3Y =-1
        I3N =-1
      ENDIF
      CALL AMFACNN
      ISIGN= I3Y*I3N
      FPI2 = DABS(FPI2)*ISIGN
      GRO2 = DABS(GRO2)*ISIGN
      FRO2 = DABS(FRO2)*ISIGN
      GA02 = DABS(GA02)*ISIGN
      GFPRO= DABS(GFPRO)*ISIGN
      GFMRO= DABS(GFMRO)*ISIGN

**    isospin dependence
   20 FACISO= 4D0*ISO-2D0
      GA0C2 = FACISO*GA0*GA0
      FPIC2 = FACISO*FPI*FPI/(AMPV*AMPV)
      GROC2 = FACISO*GRO*GRO
      FROC2 = FACISO*FRO*FRO/AMT24
      GFPROC= FACISO*(GRO*FRO+FRO*GRO)/(2D0*AMT)
      GFMROC= FACISO*(GRO*FRO-FRO*GRO)/(2D0*AMT)
**    Phenomenological phase-shift difference in PP and NP 1S0
      IF(PHNAME.EQ.'1S0' .AND. (TYPE.EQ.'NP'.OR.TYPE.EQ.'PN'))
     .   GA0C2=GA0C2*(1D0+PAR(4,3))

      FPOM2 = (GPOM*AMPOM/AMT)**2 * 4D0*AMPOM/SRPI
      FA2D2 = (GA2D*AMA2D/AMT)**2 * 4D0*AMA2D/SRPI

      RH = R/HBC
      RGT= R
      CALL RSPACFF(RH,TYPE,NADIA)
      CALL POTFNN(RH,TYPE,ISIGN,FACISO,NLOC,NADIA)

      RETURN
      END
************************************************************************
      SUBROUTINE ESCPAR
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/PARAMS/PAR(8,5)
      DIMENSION PARTR(5,8)
      DATA PARTR/
     . .2695464D+00,.6198737D+00,.3705890D+01,.1868811D+01,.0000000D+00
     .,.1632696D+00,.3419249D+01,-.2432022D+00,.5810781D+00,.1954388D+01
     .,.3550000D+00,.1000000D+01,.400000D+00,-.2457693D+01,.3091017D+03
     .,-.230000D+02,.3750000D+02,.2076657D-01,.3790000D+02,.0000000D+00
     .,.8524012D+03,.8733411D+03,.0000000D+00,.6953030D+03,.0000000D+00
     .,.9275149D+03,.1036299D+04,.0000000D+00,.7734198D+03,.0000000D+00
     .,.0000000D+00,.0000000D+00,-.4820407D+00,.5269871D-01,.000000D+00
     .,-.4120379D+00,.5831442D-01,.2161068D+00,.000000D+00,.5962126D+00/
c      total with 1s0 np : 6177.2837686797  4761 f^{2}_{pi}=  0.074500
c      total without     : 5990.9917387491  4761 f^{2}_{pi}=  0.074500

      DO 1 I=1,8
      DO 1 J=1,5
        PAR(I,J)=PARTR(J,I)
    1 CONTINUE
      WRITE(6,*) ' Nijmegen extended soft-core potential parameters'
      DO 4 I=1,8
  4     WRITE(6,5) (PAR(I,J),J=1,5)
  5   FORMAT(5(1X,D15.7))
      RETURN
      END
************************************************************************
      SUBROUTINE NYMCPFF
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/KAPPAV/AKRO,AKOM,AKFI
      COMMON/SCALIN/AMT,AMPV,AMT24
      DATA SR3/1.7320508075688772935D0/
      data iprint/0/

**    pseudovector coupling constants
      PV8  = FPI * (4D0*ALPV-1D0)/SR3
      COST = DCOS(THPV)
      SINT = DSIN(THPV)
      FETA = COST*PV8 - SINT*PV1
      FETP = SINT*PV8 + COST*PV1
      FPI2 = FPI*FPI/(AMPV*AMPV)
      FETA2= FETA*FETA/(AMPV*AMPV)
      FETP2= FETP*FETP/(AMPV*AMPV)

**    vector coupling constants
      GV8  = GRO * (4D0*ALVE-1D0)/SR3
      COST = DCOS(THV)
      SINT = DSIN(THV)
      GFI  = COST*GV8 - SINT*GV1
      GOM  = SINT*GV8 + COST*GV1
      GRO2 = GRO*GRO
      GOM2 = GOM*GOM
      GFI2 = GFI*GFI

**    tensor coupling constants
      FRO  = AKRO*GRO
      FV8  = (GRO+FRO) * (4D0*ALVM-1D0)/SR3 - GV8
      COST = DCOS(THV)
      SINT = DSIN(THV)
      FFI  = COST*FV8 - SINT*FV1
      FOM  = SINT*FV8 + COST*FV1
      AKOM = FOM/GOM
      AKFI = FFI/GFI
      FRO2 = FRO*FRO/AMT24
      FOM2 = FOM*FOM/AMT24
      FFI2 = FFI*FFI/AMT24

      GFPRO = (GRO*FRO+FRO*GRO)/(2D0*AMT)
      GFPOM = (GOM*FOM+FOM*GOM)/(2D0*AMT)
      GFPFI = (GFI*FFI+FFI*GFI)/(2D0*AMT)
      GFMRO = (GRO*FRO-FRO*GRO)/(2D0*AMT)
      GFMOM = (GOM*FOM-FOM*GOM)/(2D0*AMT)
      GFMFI = (GFI*FFI-FFI*GFI)/(2D0*AMT)

**    scalar coupling constants
      COST = DCOS(THGS)
      SINT = DSIN(THGS)
      GS8  = COST*GF0 + SINT*GEP
      GS1  =-SINT*GF0 + COST*GEP
      ALGS = (SR3*GS8+GA0)/(4D0*GA0)
      GA02 = GA0*GA0
      GEP2 = GEP*GEP
      GF02 = GF0*GF0
      if(iprint.eq.0) then
        iprint=1
        write(*,'(a,5f13.6)') ' ga0,gep,gf0:',ga0,gep,gf0,gs8,gs1
        write(*,'(a,3f13.6)') ' fpi,fet,fx0:',fpi,feta,fetp
        write(*,'(a,3f13.6)') ' gro,gom,gfi:',gro,gom,gfi
        write(*,'(a,3f13.6)') ' fro,fom,ffi:',fro,fom,ffi
        write(*,'(a,3f13.6)') ' kro,kom,kfi:',akro,akom,akfi
      endif

      RETURN
      END
************************************************************************
      SUBROUTINE CHIRAL
C*    Imposes theoretical constraints on the pair coupling constants
C*    based on chiral Langrangians (see notes V. Stoks),
C*    using the linear sigma model.
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/SCALIN/AMT,AMPV,AMT24
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/GPAIR/ GPIETA,GPIETP,GSISI, GPISIG,GPIPOM,
     .              GPIPI0,GPIPI1,FPIPI1,GPIRO0,GPIRO1,GPIOM
      COMMON/PARAMS/PAR(8,5)
      DATA FPIDEC/92.4D0/, GAXI/1.2573D0/
      DATA PI4/12.5663706143592D0/, SQ4PI/3.54490770181103D0/
      DATA ICAL/0/, SQ2/1.4142135623731D0/, SQ6/2.4494897427832D0/
      DATA IC/-1/, iw/0/
      SAVE AMPV2,ROM2,AMA1,A1M2,ZM1,CHIRPAR,GROCHIR

      IF(ICAL.EQ.0) THEN
        ICAL=1
        AMPV2=AMPV*AMPV
        ROM2=ROM*ROM
        AMA1=1230D0
        A1M2=AMA1*AMA1
        ZM1 =ROM2/A1M2
        CHIRPAR = ZM1*(GAXI-IC*ZM1)/(2D0*FPIDEC**2)
        GROCHIR = DSQRT(1D0-ZM1)*ROM/(2D0*FPIDEC)
      ENDIF
c     if(iw.eq.0) write(*,*) ' axial:',gaxi/(2d0*fpidec*sq4pi),fpi/ampv
c     if(iw.eq.0) write(*,*) 'vector:',grochir/sq4pi,gro

      GPIPI1 =-IC*AMPV2*CHIRPAR/PI4
      FPIPI1 = GPIPI1 * FRO/GRO

      GPISIG = (A1M2/ROM2-2D0)*AMA1/ROM*AMPV2*CHIRPAR/PI4

      GPIRO1 = GRO*AMPV/FPIDEC*(GAXI-IC)/SQ4PI

      GPIPI0 =-(AMT/FPIDEC/PI4-GEP/SQ4PI*AMA1/ROM)*AMPV/(2D0*FPIDEC)
      GSISI  = 3D0*ZM1*GPIPI0
      GPIETA = 2D0*ROM/AMA1*GPIPI0
      gpieta = 0d0

      if(iw.eq.0) then
        iw=1
        write(*,'(a,2f13.6)') ' gpipi,fpipi:',gpipi1,fpipi1
        write(*,'(a,2f13.6)') ' gpisi,gpiro:',gpisig,gpiro1
        write(*,'(a,3f13.6)') ' p02,si2,pet:',gpipi0,gsisi,gpieta
      endif

C*    Store the pair couplings back into the parameter matrix
      PAR(7,1) = GPIETA
      PAR(7,2) = GPIETP
      PAR(7,3) = GSISI
      PAR(7,4) = GPISIG
      PAR(7,5) = GPIPOM
      PAR(8,1) = GPIPI0
      PAR(8,2) = GPIPI1
      PAR(8,3) = FPIPI1
      PAR(8,4) = GPIRO0
      PAR(8,5) = GPIRO1
      PAR(6,5) = GPIOM

      RETURN
      END
************************************************************************
      SUBROUTINE AMFACNN
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2,
     .              AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2,
     .              AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2

      AMY2  = AMY*AMY
      AMYI  = 1D0/AMY
      AMYI2 = AMYI*AMYI
      AMN2  = AMN*AMN
      AMNI  = 1D0/AMN
      AMNI2 = AMNI*AMNI
      AMYPN = AMY+AMN
      AMYMN = AMY-AMN
      AMYPN2= AMYPN*AMYPN
      AMYMN2= AMYMN*AMYMN
      AY2PN2= AMY2+AMN2
      AY2MN2= AMY2-AMN2
      AMYN  = AMY*AMN
      AMYNI = AMYI*AMNI
      AMYNI2= AMYNI*AMYNI
      AYPNI2= AMYI2+AMNI2
      AYMNI2= AMYI2-AMNI2
      REDM = AMYN/AMYPN
      ALF  = 4D0*REDM/AMYPN
      IF(AMY.EQ.AMN) ALF = 1D0
      RETURN
      END
************************************************************************
      SUBROUTINE RSPACFF(RH,TYPE,NADIA)
**      pseudoscalar functions P1-4: pi+, pi0, eta, eta'
**      scalar functions       S1-4: a0+, a00, epsilon, f0
**      vector functions      VV1-4: g*g in rho+, rho0, omega, phi
**                            VT1-4: g*f in rho+, rho0, omega, phi
**                            TT1-4: f*f in rho+, rho0, omega, phi
**                              VV5: g*g scalar part of rho+
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER TYPE*2
      LOGICAL BRDRO,BRDROC,BRDEP
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/BRDRHO/GAMRO,THRRO,NRO,BRDRO,GAMROC,THRROC,NROC,BRDROC
      COMMON/BRDEPS/GAMEP,THREP0,THREPC,NEP,BRDEP
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/CUTTEN/ALAMT1,ALAVT1,ALAMT0,ALAVT0
      COMMON/SPCTRL/RHB,AMTHR2,AMMES,GAMMES,ALABDA,N
      COMMON/RVMAT/ P1(6),P2(6),P3(6),P4(6),S1(6),S2(6),S3(6),S4(6),
     .              VV1(6),VV2(6),VV3(6),VV4(6),VT1(6),VT2(6),
     .              VT3(6),VT4(6),TT1(6),TT2(6),TT3(6),TT4(6),VV5(6)
      COMMON/RBMAT/ RES(6,2),REV(6,6)
      COMMON/RPINA/ F3PIS(3,2),F3PIT(3,2)
      RHB = RH

C**   Pseudoscalar mesons pi, eta, eta'
      CALL REGGEFF(RH,AMPI,ALAMP1,P2)
      IF(NADIA.EQ.1) CALL RNADIA(RH,TYPE,P2,F3PIS(1,2),F3PIT(1,2))
      CALL REGGEFF(RH,AMETA,ALAMP0,P3)
      CALL REGGEFF(RH,AMETP,ALAMP0,P4)

c**   Scalar mesons a0(980), epsilon, f0(975)
      CALL REGGEFF(RH,AMA0,ALAMS1,S2)
      IF(BRDEP) THEN
        AMMES = AMEP
        GAMMES= GAMEP
        ALABDA= ALAMS0
        N = NEP
        AMTHR2=THREP0*THREP0
        CALL BROAD(THREP0,RES(1,1),6)
        AMTHR2=THREPC*THREPC
        CALL BROAD(THREPC,RES(1,2),6)
        DO 10 I=1,6
          S3(I) = (RES(I,1) + 2D0*RES(I,2))/3D0
   10   CONTINUE
      ELSE
        CALL REGGEFF(RH,AME1,ALAMS0,RES(1,1))
        CALL REGGEFF(RH,AME2,ALAMS0,RES(1,2))
        DO 11 I=1,6
          S3(I) = AEPS*RES(I,1) + BEPS*RES(I,2)
   11   CONTINUE
      ENDIF
      CALL REGGEFF(RH,AMF0,ALAMS0,S4)

C**   Vector mesons rho, omega, phi
      IF(BRDRO) THEN
        AMTHR2= THRRO*THRRO
        AMMES = AMRO
        GAMMES= GAMRO
        N = NRO
        ALABDA= ALAMV1
        CALL BROAD(THRRO,VV2,6)
        IF(ALAMT1.NE.ALAMV1) THEN
          ALABDA= ALAVT1
          CALL BROAD(THRRO,VT2,6)
          ALABDA= ALAMT1
          CALL BROAD(THRRO,TT2,6)
        ELSE
          DO 21 I=1,6
            VT2(I) = VV2(I)
            TT2(I) = VV2(I)
   21     CONTINUE
        ENDIF
      ELSE
        CALL REGGEFF(RH,AMR1,ALAMV1,REV(1,1))
        CALL REGGEFF(RH,AMR2,ALAMV1,REV(1,2))
        IF(ALAMT1.NE.ALAMV1) THEN
          CALL REGGEFF(RH,AMR1,ALAVT1,REV(1,3))
          CALL REGGEFF(RH,AMR2,ALAVT1,REV(1,4))
          CALL REGGEFF(RH,AMR1,ALAMT1,REV(1,5))
          CALL REGGEFF(RH,AMR2,ALAMT1,REV(1,6))
        ENDIF
        DO 22 I=1,6
          VV2(I) = ARO*REV(I,1) + BRO*REV(I,2)
          IF(ALAMT1.NE.ALAMV1) THEN
            VT2(I) = ARO*REV(I,3) + BRO*REV(I,4)
            TT2(I) = ARO*REV(I,5) + BRO*REV(I,6)
          ELSE
            VT2(I) = ARO*REV(I,1) + BRO*REV(I,2)
            TT2(I) = ARO*REV(I,1) + BRO*REV(I,2)
          ENDIF
   22   CONTINUE
      ENDIF
      CALL REGGEFF(RH,AMOM,ALAMV0,VV3)
      CALL REGGEFF(RH,AMFI,ALAMV0,VV4)
      IF(ALAMT0.NE.ALAMV0) THEN
        CALL REGGEFF(RH,AMOM,ALAVT0,VT3)
        CALL REGGEFF(RH,AMOM,ALAMT0,TT3)
        CALL REGGEFF(RH,AMFI,ALAVT0,VT4)
        CALL REGGEFF(RH,AMFI,ALAMT0,TT4)
      ELSE
        DO 23 I=1,6
          VT3(I) = VV3(I)
          TT3(I) = VV3(I)
          VT4(I) = VV4(I)
          TT4(I) = VV4(I)
   23   CONTINUE
      ENDIF

C**   Isospin-1 mesons (pi+, rho+, a0+)
      IF(TYPE.EQ.'NP' .OR. TYPE.EQ.'PN') THEN
        CALL REGGEFF(RH,AWPIC,ALAMP1,P1)
        IF(NADIA.EQ.1) CALL RNADIA(RH,TYPE,P1,F3PIS(1,1),F3PIT(1,1))
        IF(BRDROC) THEN
          AMTHR2= THRROC*THRROC
          AMMES = AMROC
          GAMMES= GAMROC
          N = NROC
          ALABDA= ALAMV1
          CALL BROAD(THRROC,VV1,6)
          IF(ALAMT1.NE.ALAMV1) THEN
            ALABDA= ALAVT1
            CALL BROAD(THRROC,VT1,6)
            ALABDA= ALAMT1
            CALL BROAD(THRROC,TT1,6)
          ELSE
            DO 31 I=1,6
              VT1(I) = VV1(I)
              TT1(I) = VV1(I)
   31       CONTINUE
          ENDIF
          DO 32 I=1,6
            VV5(I) = AVSC*VV1(I)
   32     CONTINUE
        ELSE
          DO 33 I=1,6
            VV1(I) = AROC*REV(I,1) + BROC*REV(I,2)
            IF(ALAMT1.NE.ALAMV1) THEN
              VT1(I) = AROC*REV(I,3) + BROC*REV(I,4)
              TT1(I) = AROC*REV(I,5) + BROC*REV(I,6)
            ELSE
              VT1(I) = AROC*REV(I,1) + BROC*REV(I,2)
              TT1(I) = AROC*REV(I,1) + BROC*REV(I,2)
            ENDIF
            VV5(I) = AVSC1*AROC*REV(I,1) + AVSC2*BROC*REV(I,2)
   33     CONTINUE
        ENDIF
        CALL REGGEFF(RH,AWA0C,ALAMS1,S1)
      ENDIF

      RETURN
      END
************************************************************************
      SUBROUTINE REGGEFF(RH,AM,ALAM,FIFUNS)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FIFUNS(6)
      DATA ZMATCH/1D-2/, SQPI/1.7724538509055160273D0/

      RATM=AM/ALAM
      RATM2=RATM*RATM
      RATM3=RATM*RATM2
      EXPM=FDEXP(RATM2)
      AM3=AM*AM*AM
      AM5=AM*AM*AM3
      Z=RH*ALAM

      IF(Z.GT.ZMATCH) THEN
        XLAM=0.5D0*RH*ALAM
        XLAM2=XLAM*XLAM
        ELAM=FDEXP(-XLAM2)
        X=AM*RH
        ERFCM=FDERFC(RATM-XLAM)
        ERFCP=FDERFC(RATM+XLAM)
        EPX=0D0
        IF(ERFCP.NE.0) EPX=FDEXP(X)
        EMX=FDEXP(-X)

        PHI0C = EXPM*(EMX*ERFCM-EPX*ERFCP)/(2D0*X)
        PHI0L =(EXPM*(EMX*ERFCM*(1D0+X)-EPX*ERFCP*(1D0-X))-
     .          4D0/SQPI*XLAM*ELAM ) / (2D0*X*X*X)
        PHI0T =(EXPM*(EMX*ERFCM*(1D0+X+X*X/3D0)
     .               -EPX*ERFCP*(1D0-X+X*X/3D0))
     .         -8D0/SQPI*XLAM*(0.5D0+XLAM2/3D0)*ELAM) / (2D0*X*X*X)
        DEL1C =-ELAM/(RATM3*2D0*SQPI)
        DEL1L =-ELAM/(RATM2*RATM3*4D0*SQPI)
        DEL1T =-ELAM*XLAM2/(RATM2*RATM3*6D0*SQPI)
        DEL2C = 3D0*(DEL1T-DEL1L)
      ELSE
        A=RATM
        B=-FDERFC(A)*A*EXPM*SQPI
        A2=A*A
        A3=A2*A
        A4=A2*A2
        A5=A3*A2
        A6=A4*A2
        A8=A6*A2
        A10=A8*A2
        Z2=Z*Z
        Z4=Z2*Z2
        Z6=Z4*Z2
        Z8=Z6*Z2
        PHI0C = (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z8/5806080D0+
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z6/40320D0+
     .          (4*A4*B+4*A4-2*A2+3)*Z4/480D0+
     .          (2*A2*B+2*A2-1)*Z2/12D0+B+1D0
        PHI0C = PHI0C/(SQPI*A)
        PHI0T = (32*A10*B+32*A10-16*A8+24*A6-60*A4+210*A2-945)
     .           *Z8/47900160D0+
     .          (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z6/362880D0+
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z4/5040D0+
     .          (4*A4*B+4*A4-2*A2+3)*Z2/180D0
        PHI0T = PHI0T/(SQPI*A3)
        PHI0L =-(32*A10*B+32*A10-16*A8+24*A6-60*A4+210*A2-945)
     .           *Z8/127733760D0-
     .          (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z6/725760D0-
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z4/6720D0-
     .          (4*A4*B+4*A4-2*A2+3)*Z2/120D0-(2*A2*B+2*A2-1D0)/6D0
        PHI0L = PHI0L/(SQPI*A3)
        DEL1C =-Z8/12288D0+Z6/768D0-Z4/64D0+Z2/8D0-0.5D0
        DEL1C = DEL1C/(SQPI*A3)
        DEL1T = Z8/9216D0-Z6/768D0+Z4/96D0-Z2/24D0
        DEL1T = DEL1T/(SQPI*A5)
        DEL1L =-Z8/24576D0+Z6/1536D0-Z4/128D0+Z2/16D0-0.25D0
        DEL1L = DEL1L/(SQPI*A5)
        DEL2C = 3D0*(DEL1T-DEL1L)
      ENDIF

      FIFUNS(1) = AM*PHI0C
      FIFUNS(2) = AM3*(PHI0C+DEL1C)
      FIFUNS(3) = AM5*(PHI0C+DEL1C+DEL2C)
      FIFUNS(4) = AM3*PHI0L
      FIFUNS(5) = AM5*(PHI0L+DEL1L)
      FIFUNS(6) = AM3*PHI0T
c     FIFUNS(7) = AM5*(PHI0T+DEL1T)

      RETURN
      END
************************************************************************
      SUBROUTINE BROAD(AMTHR,SUM,NDIM)
**      Calculates integrals for broad mesons using Gauss integration
**      on (mthr,infinity) by mapping (mthr,3*mthr) on (-1,1) with 16
**      points, and (3*mthr,20*mthr) on (-1,1) with 8 points.
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SUM(NDIM),FUNCW(6),YB(3),NP(2),X(20,2),W(20,2)
      DATA ICAL/0/, YB/1D0,3D0,20D0/, NP/16,8/, NI/2/
      SAVE X,W

      IF(ICAL.EQ.0) THEN
        CALL GAUSS(NP(1),X(1,1),W(1,1),1)
        CALL GAUSS(NP(2),X(1,2),W(1,2),1)
        ICAL=1
      ENDIF

      DO 1 I=1,NDIM
        SUM(I)=0D0
    1 CONTINUE
      DO 2 IN=1,NI
        BMA=0.5D0*(YB(IN+1)-YB(IN))*AMTHR
        BPA=0.5D0*(YB(IN+1)+YB(IN))*AMTHR
        DO 3 IY=1,NP(IN)
          AMY=BMA*X(IY,IN)+BPA
          CALL BRDFUN(FUNCW,AMY)
          DO 4 K=1,NDIM
            SUM(K)=SUM(K)+FUNCW(K)*W(IY,IN)*BMA
    4     CONTINUE
    3   CONTINUE
    2 CONTINUE

      RETURN
      END
************************************************************************
      SUBROUTINE BRDFUN(FINTGR,AMP)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FI(6),FINTGR(6)
      COMMON/SPCTRL/RB,AMTHR2,AMMES,GAMMES,ALAM,N
      DATA PI/3.14159265358979323846D0/

      AMP2=AMP*AMP
      AMMES2=AMMES*AMMES
      FACG=AMMES*GAMMES/(AMMES2-AMTHR2)**(N+0.5D0)
      SPEC=FACG*(AMP2-AMTHR2)**(N+0.5D0)/
     .     ( (AMP2-AMMES2)**2+FACG*FACG*(AMP2/AMMES2)**(2D0*N)*
     .       (AMP2-AMTHR2)**(2D0*N+1D0) ) / PI
      CALL REGGEFF(RB,AMP,ALAM,FI)
      DO 1 I=1,6
        FINTGR(I)=2D0*AMP*SPEC*FI(I)
    1 CONTINUE

      RETURN
      END
************************************************************************
      SUBROUTINE RNADIA(RH,TYPE,F2PI,F3PIS,F3PIT)
C*     Calculates the I_{3} functions with first and second derivatives
C*     of spin-spin and tensor OPE functions.
C*     These derivatives can be written in terms of
C*           PHI^{1}_{C}, PHI^{2}_{C}, PHI^{1}_{SO}, PHI^{0}_{T}
C-----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 F2PI(6),FILAM(6), X(20),W(20), SUM(4), F3PIS(3),F3PIT(3)
      CHARACTER TYPE*2
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      DATA NP/12/, ICAL/0/, TWONPI/0.636619772367581343D0/

      IF(ICAL.EQ.0) THEN
        CALL GAUSS(NP,X,W,1)
        DO 1 I=1,NP
          W(I)=2D0*W(I)/(1D0-X(I))**2
          X(I)=(1D0+X(I))/(1D0-X(I))
    1   CONTINUE
        ICAL=1
      ENDIF

      IF(TYPE.EQ.'PP' .OR. TYPE.EQ.'NN') THEN
        AMES=AMPI
      ELSEIF(TYPE.EQ.'NP' .OR. TYPE.EQ.'PN') THEN
        AMES=AWPIC
      ENDIF

      DO 10 K=1,4
   10   SUM(K)=0D0
      DO 11 IY=1,NP
        PLAM=X(IY)*PIM
        PLAM2=PLAM*PLAM
        CALL PHIFUN(RH,AMES,ALAMP1,PLAM,FILAM)
        SUM(1)=SUM(1)+(F2PI(2)-FILAM(2))/PLAM2*W(IY)
        SUM(2)=SUM(2)+(F2PI(3)-FILAM(3))/PLAM2*W(IY)
        SUM(3)=SUM(3)+(F2PI(5)-FILAM(5))/PLAM2*W(IY)
        SUM(4)=SUM(4)+(F2PI(6)-FILAM(6))/PLAM2*W(IY)
   11 CONTINUE
**    spin (0th, 1st, 2nd) and tensor (0th, 1st, 2nd)
      F3PIS(1) = PIM*TWONPI * SUM(1)
      F3PIS(2) =-PIM*TWONPI * SUM(3)*RH
      F3PIS(3) = PIM*TWONPI *(SUM(2)+2D0*SUM(3))
      F3PIT(1) = PIM*TWONPI * SUM(4)
      F3PIT(2) =-PIM*TWONPI *(3D0/RH*SUM(4)+RH/3D0*SUM(3))
      F3PIT(3) = PIM*TWONPI *(12D0/(RH*RH)*SUM(4)
     .                       +(5D0*SUM(3)+SUM(2))/3D0)

      RETURN
      END
***********************************************************************
      SUBROUTINE PHIFUN(RH,AM,ALAM,PLAM,PHI)
**     Basic function phi^{n}_{X} stored as
**     PHI = (m phi^{0}_{C}, m^{3} phi^{1}_{C}, m^{5} phi^{2}_{C},
**            m^{3} phi^{0}_{SO}, m^{5} phi^{1}_{SO},
**            m^{3} phi^{0}_{T}, m^{5} phi^{1}_{T})
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION PHI(6)
      DATA SRPI/1.7724538509055160273D0/

      AMP  = DSQRT(AM*AM+PLAM*PLAM)
      AMP2 = AMP*AMP
      AMP3 = AMP*AMP2
      RATM = AMP/ALAM
      RATM2= RATM*RATM
      RATM3= RATM*RATM2
      EXPM = FDEXP(AM*AM/(ALAM*ALAM))
      EXPL = FDEXP(-PLAM*PLAM/(ALAM*ALAM))
      XLAM = 0.5D0*RH*ALAM
      XLAM2= XLAM*XLAM
      ELAM = FDEXP(-XLAM2)
      X = AMP*RH
      X2= X*X
      X3= X*X2
      ERFCM = FDERFC(RATM-XLAM)
      ERFCP = FDERFC(RATM+XLAM)
      EPX = 0D0
      IF(ERFCP.NE.0) EPX = FDEXP(X)
      EMX = FDEXP(-X)

      PHI0C = EXPM*(EMX*ERFCM-EPX*ERFCP)/(2D0*X)
      PHI0L = (EXPM*(EMX*ERFCM*(1D0+X)-EPX*ERFCP*(1D0-X))
     .      - 4D0/SRPI*XLAM*ELAM*EXPL ) / (2D0*X3)
      PHI0T =(EXPM*(EMX*ERFCM*(1D0+X+X2/3D0)-EPX*ERFCP*(1D0-X+X2/3D0))
     .      - 8D0/SRPI*XLAM*(0.5D0+XLAM2/3D0)*ELAM*EXPL ) / (2D0*X3)
      DEL1C =-ELAM/(RATM3*2D0*SRPI) * EXPL
      PHI1C = PHI0C+DEL1C
      DEL1L =-ELAM/(RATM2*RATM3*4D0*SRPI) * EXPL
      PHI1L = PHI0L+DEL1L
      DEL1T =-ELAM*XLAM2/(RATM2*RATM3*6D0*SRPI) * EXPL
      PHI1T = PHI0T+DEL1T
      PHI2C = PHI1C+3D0*(DEL1T-DEL1L)
      PHI(1) = AMP *PHI0C
      PHI(2) = AMP3*PHI1C
      PHI(3) = AMP3*PHI2C*AMP2
      PHI(4) = AMP3*PHI0L
      PHI(5) = AMP3*PHI1L*AMP2
      PHI(6) = AMP3*PHI0T
c     PHI(7) = AMP3*PHI1T*AMP2

      RETURN
      END
************************************************************************
      SUBROUTINE POTFNN(RH,TYPE,ISIGN,FACISO,NLOC,NADIA)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER TYPE*2
      DIMENSION DUM1(9),DUM2(9)
      COMMON/POT2/  VC,VSS,VT,VLS,VLSA,VQ12,FIC0,FIC1,FIC2
      COMMON/VOPE/  VPIS,VPIT,VPIQ,FPIS0,FPIS1,FPIS2,FPIT0,FPIT1,FPIT2
      COMMON/VENER/ PLMEV2,VEPIS,VEPIT
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/SCALIN/AMT,AMPV,AMT24
      COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4,
     .              GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4
      COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2,
     .              AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2,
     .              AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2
      COMMON/RVMAT/ P1(6),P2(6),P3(6),P4(6),S1(6),S2(6),S3(6),S4(6),
     .              VV1(6),VV2(6),VV3(6),VV4(6),VT1(6),VT2(6),
     .              VT3(6),VT4(6),TT1(6),TT2(6),TT3(6),TT4(6),VV5(6)
      COMMON/RPINA/ F3PIS(3,2),F3PIT(3,2)
      EQUIVALENCE (VC,DUM1(1)), (VPIS,DUM2(1))

      DO 1 ID=1,9
        DUM1(ID)=0D0
        DUM2(ID)=0D0
    1 CONTINUE

      RH2=RH*RH

C*    Pseudoscalar-meson potential
      F2PION=FPI2
      IF(NADIA.EQ.1) THEN
**      Nonadiabatic correction to OPE, energy multiplication outside
        FISP0 =-0.5D0*FPI2*F3PIS(1,2)/3D0
        FISP1 =-0.5D0*FPI2*F3PIS(2,2)/3D0
        FISP2 =-0.5D0*FPI2*F3PIS(3,2)/3D0
        FITP0 =-0.5D0*FPI2*F3PIT(1,2)
        FITP1 =-0.5D0*FPI2*F3PIT(2,2)
        FITP2 =-0.5D0*FPI2*F3PIT(3,2)
        VEPIS = 0.5D0*FPI2*F3PIS(1,2)/(3D0*REDM)
        VEPIT = 0.5D0*FPI2*F3PIT(1,2)/REDM
      ELSEIF(NLOC.EQ.1) THEN
**      Energy-dependent OPE, added separately
        CONTINUE
      ELSEIF(NLOC.EQ.2) THEN
**      Nonlocal spin-spin and tensor
        FFNL = REDM*FPI2*AMYNI/2D0
        FPIS0=-FFNL*P2(2)/3D0
        FPIS1= FFNL*P2(5)*RH/3D0
        FPIS2=-FFNL*(P2(3)+2D0*P2(5))/3D0
        FPIT0=-FFNL*P2(6)
        FPIT1= FFNL*(3D0*P2(6)/RH+P2(5)*RH/3D0)
        FPIT2=-FFNL*(12D0*P2(6)/RH2+(5D0*P2(5)+P2(3))/3D0)
      ENDIF
      VPIS = FPI2*P2(2)/3D0
      VPIT = FPI2*P2(6)
      VPVS = (FETA2*P3(2)+FETP2*P4(2))/3D0
      VPVT = (FETA2*P3(6)+FETP2*P4(6))
      FPI2=F2PION

C*    Scalar-meson potential
      VSC  =-(GA02*S2(1)+GEP2*S3(1)+GF02*S4(1)) +
     .       (GA02*S2(2)+GEP2*S3(2)+GF02*S4(2)) * AYPNI2/8D0
      VSLS =-(GA02*S2(4)+GEP2*S3(4)+GF02*S4(4)) * AYPNI2/4D0
      VSLA =-(GA02*S2(4)+GEP2*S3(4)+GF02*S4(4)) * AYMNI2/4D0
      VSQ12=-(GA02*S2(6)+GEP2*S3(6)+GF02*S4(6)) * AMYNI2/16D0
     .       *3D0/RH2
      FICS0= REDM*(GA02*S2(1)+GEP2*S3(1)+GF02*S4(1))*AMYNI/2D0
      FICS1=-REDM*(GA02*S2(4)+GEP2*S3(4)+GF02*S4(4))*AMYNI/2D0*RH
      FICS2= REDM*(GA02*S2(2)+GEP2*S3(2)+GF02*S4(2))*AMYNI/2D0
     .     + REDM*(GA02*S2(4)+GEP2*S3(4)+GF02*S4(4))*AMYNI

C*    Diffractive contribution
      XI0 = AMPOM*RH
      XI02 = XI0*XI0
      EXI0 = 0D0
      IF(XI02.LT.170.) EXI0 = FDEXP(-XI02)
      XI1 = AMA2D*RH
      XI12 = XI1*XI1
      EXI1 = 0D0
      IF(XI12.LT.170.) EXI1 = FDEXP(-XI12)
      VDC  = FPOM2 * (1D0+AMPOM2*(1.5D0-XI02)*AYPNI2/2D0) * EXI0 +
     .       FA2D2 * (1D0+AMA2D2*(1.5D0-XI12)*AYPNI2/2D0) * EXI1*ISIGN
      VDLS = (FPOM2 * AMPOM2 * EXI0 +
     .        FA2D2 * AMA2D2 * EXI1*ISIGN) * AYPNI2/2D0
      VDLA = (FPOM2 * AMPOM2 * EXI0 +
     .        FA2D2 * AMA2D2 * EXI1*ISIGN) * AYMNI2/2D0
      VDQ12= (FPOM2 * AMPOM4 * EXI0 +
     .        FA2D2 * AMA2D4 * EXI1*ISIGN) * AMYNI2/4D0
      FICD0=-REDM*(FPOM2*EXI0 + ISIGN*FA2D2*EXI1) * AMYNI/2D0
      FICD1= REDM*(FPOM2*EXI0*XI02 + ISIGN*FA2D2*EXI1*XI12) * AMYNI/RH
      FICD2= REDM*(FPOM2*EXI0*(1D0-2D0*XI02)*AMPOM2 +
     .       ISIGN*FA2D2*EXI1*(1D0-2D0*XI12)*AMA2D2) * AMYNI

C*    Vector-meson potential
      VVC  = (GRO2*VV2(1)+GOM2*VV3(1)+GFI2*VV4(1))
     .      +(GRO2*VV2(2)+GOM2*VV3(2)+GFI2*VV4(2)) * AMYNI/(2D0*ALF)
     .      +(gro2*vv2(3)+gom2*vv3(3)+gfi2*vv4(3)) * amyni2/64d0
     .      +(GFPRO*VT2(2)+GFPOM*VT3(2)+GFPFI*VT4(2)) / (4D0*REDM)
     .      +(gfpro*vt2(3)+gfpom*vt3(3)+gfpfi*vt4(3))*amyni/(16d0*redm)
     .      +(GFMRO*VT2(2)+GFMOM*VT3(2)+GFMFI*VT4(2)) * AMYMN*AMYNI/4D0
     .      +(FRO2*TT2(3)+FOM2*TT3(3)+FFI2*TT4(3))*AMYNI/4D0
      VVS  = (GRO2*VV2(2)+GOM2*VV3(2)+GFI2*VV4(2)) * AMYNI/6D0
     .      +(GFPRO*VT2(2)+GFPOM*VT3(2)+GFPFI*VT4(2))*AMYPN*AMYNI/6D0
     .      -(GFMRO*VT2(2)+GFMOM*VT3(2)+GFMFI*VT4(2))*AMYMN*AMYNI/6D0
     .      +(FRO2*TT2(2)+FOM2*TT3(2)+FFI2*TT4(2))*2D0/3D0
C    .      +(FRO2*TT2(3)+FOM2*TT2(3)+FFI2*TT2(3))*AMYNI/12D0
      VVT  =-(GRO2*VV2(6)+GOM2*VV3(6)+GFI2*VV4(6)) * AMYNI/4D0
     .      -(GFPRO*VT2(6)+GFPOM*VT3(6)+GFPFI*VT4(6))*AMYPN*AMYNI/4D0
     .      +(GFMRO*VT2(6)+GFMOM*VT3(6)+GFMFI*VT4(6))*AMYMN*AMYNI/4D0
     .      -(FRO2*TT2(6)+FOM2*TT3(6)+FFI2*TT4(6))
C    .      -(FRO2*TT2(7)+FOM2*TT2(7)+FFI2*TT2(7))*AMYNI/8D0
      VVLS =-(GRO2*VV2(4)+GOM2*VV3(4)+GFI2*VV4(4))*(.5D0+1D0/ALF)*AMYNI
     .      -(gro2*vv2(5)+gom2*vv3(5)+gfi2*vv4(5))*amyni2/16d0
     .      -(GFPRO*VT2(4)+GFPOM*VT3(4)+GFPFI*VT4(4))/REDM
     .      -(gfpro*vt2(5)+gfpom*vt3(5)+gfpfi*vt4(5))*amyni/(4d0*redm)
     .      -(FRO2*TT2(5)+FOM2*TT3(5)+FFI2*TT4(5))*(.5D0+1D0/ALF)*AMYNI
      VVLA =-(GRO2*VV2(4)+GOM2*VV3(4)+GFI2*VV4(4))*AYMNI2/4D0
     .      +(GFMRO*VT2(4)+GFMOM*VT3(4)+GFMFI*VT4(4))/REDM
     .      +(FRO2*TT2(5)+FOM2*TT3(5)+FFI2*TT4(5))*AYMNI2/4D0
      VVQ12=((GRO2*VV2(6)+GOM2*VV3(6)+GFI2*VV4(6))
     .      +(GFPRO*VT2(6)+GFPOM*VT3(6)+GFPFI*VT4(6))*4D0*AMYPN
     .      +(FRO2*TT2(6)+FOM2*TT3(6)+FFI2*TT4(6))*8D0*AMYPN2)
     .      *AMYNI2/16D0 * 3D0/RH2
      FICV0= REDM*(GRO2*VV2(1)+GOM2*VV3(1)+GFI2*VV4(1)) *
     .            (AYPNI2+AMYNI)/2D0
      FICV1=-REDM*(GRO2*VV2(4)+GOM2*VV3(4)+GFI2*VV4(4)) *
     .            (AYPNI2+AMYNI)/2D0 * RH
      FICV2= REDM*(GRO2*VV2(2)+GOM2*VV3(2)+GFI2*VV4(2)) *
     .            (AYPNI2+AMYNI)/2D0
     .     + REDM*(GRO2*VV2(4)+GOM2*VV3(4)+GFI2*VV4(4)) *(AYPNI2+AMYNI)

      IF(TYPE.EQ.'NP' .OR. TYPE.EQ.'PN') THEN
**      Pseudoscalar: pi+
        F2PION=FPIC2
          IF(NADIA.EQ.1) THEN
**        Nonadiabatic correction to OPE, energy multiplication outside
            FPIS0=FPIS0 - 0.5D0*FPIC2*F3PIS(1,1)/3D0
            FPIS1=FPIS1 - 0.5D0*FPIC2*F3PIS(2,1)/3D0
            FPIS2=FPIS2 - 0.5D0*FPIC2*F3PIS(3,1)/3D0
            FPIT0=FPIT0 - 0.5D0*FPIC2*F3PIT(1,1)
            FPIT1=FPIT1 - 0.5D0*FPIC2*F3PIT(2,1)
            FPIT2=FPIT2 - 0.5D0*FPIC2*F3PIT(3,1)
            VEPIS=VEPIS + 0.5D0*FPIC2*F3PIS(1,1)/(3D0*REDM)
            VEPIT=VEPIT + 0.5D0*FPIC2*F3PIT(1,1)/REDM
          ELSEIF(NLOC.EQ.1) THEN
**          Energy-dependent OPE, energy multiplication outside
            CONTINUE
          ELSEIF(NLOC.EQ.2) THEN
            FFNL = REDM*FPIC2*(AMYNI-AYPNI2/4D0)
            FPIS0=FPIS0 - FFNL*P1(2)/3D0
            FPIS1=FPIS1 + FFNL*P1(5)*RH/3D0
            FPIS2=FPIS2 - FFNL*(P1(3)+2D0*P1(5))/3D0
            FPIT0=FPIT0 - FFNL*P1(6)
            FPIT1=FPIT1 + FFNL*(3D0*P1(6)/RH+P1(5)*RH/3D0)
            FPIT2=FPIT2 - FFNL*(12D0*P1(6)/RH2+(5D0*P1(5)+P1(3))/3D0)
          ENDIF
          VPIS = VPIS + FPIC2*P1(2)/3D0
          VPIT = VPIT + FPIC2*P1(6)
c       ENDIF
        FPIC2=F2PION

**      Scalar: a0+
        VSC  = VSC  - (GA0C2*S1(1) - GA0C2*S1(2)*AMYNI/4D0)
        VSLS = VSLS - GA0C2*S1(4)*AMYNI/2D0
        VSQ12= VSQ12- GA0C2*S1(6)*AMYNI2/16D0 * 3D0/RH2
        FICS0= FICS0- REDM*GA0C2*S1(1) * (AYPNI2/4D0-AMYNI)
        FICS1= FICS1+ REDM*GA0C2*S1(4) * (AYPNI2/4D0-AMYNI)*RH
        FICS2= FICS2- REDM*GA0C2*(AYPNI2/4D0-AMYNI)*(S1(2)+2D0*S1(4))

**      Diffractive: pomeron (A2)
        VDC  = VDC  + FACISO*FA2D2*
     .                (1D0+AMA2D2*(1.5D0-XI12)*AMYNI) * EXI1
        VDLS = VDLS + FACISO*FA2D2 * AMA2D2*AMYNI * EXI1
        VDQ12= VDQ12+ FACISO*FA2D2 * AMA2D4*AMYNI2/4D0 * EXI1
        FICD0= FICD0+ REDM * FACISO*FA2D2*(AYPNI2/4D0-AMYNI)*EXI1
        FICD1= FICD1- REDM * FACISO*FA2D2*(AYPNI2/4D0-AMYNI)*EXI1
     .                     * 2D0*XI12/RH
        FICD2= FICD2- REDM * FACISO*FA2D2*(AYPNI2/4D0-AMYNI)*EXI1 *
     .                AMA2D2*(2D0-4D0*XI12)

**      Vector: rho+
        VVC  = VVC + GROC2*(VV1(1)+AMYNI*VV1(2)/(2D0*ALF))
     .             + groc2*amyni2/64d0*vv1(3)
     .       + GFPROC*AMYPN*AYPNI2/8D0*VT1(2)
     .       + gfproc*amypn*amyni2*vt1(3)/16d0
     .       + FROC2*AMYMN2*AYPNI2/8D0*TT1(2)
     .       + froc2*aypni2/8d0*tt1(3)
        VVS  = VVS   + (1D0/ALF-1D0/3D0)*AMYNI/4D0 *
     .         (GROC2*VV1(2)+AMYPN*GFPROC*VT1(2)+AMYPN2*FROC2*TT1(2))
C    .       + FROC2*AYPNI2/24D0*TT1(3)
        VVT  = VVT   - AMYNI/4D0 *
     .         (GROC2*VV1(6)+AMYPN*GFPROC*VT1(6)+AMYPN2*FROC2*TT1(6))
C    .       - FROC2*AYPNI2/16D0*TT1(7)
        VVLS = VVLS  - GROC2/2D0*(AYPNI2+AMYNI)*VV1(4)
     .               - groc2*amyni2/16d0*vv1(5)
     .               - GFPROC/2D0*AMYPN*AYPNI2*VT1(4)
     .               - gfproc/4d0*amypn*amyni2*vt1(5)
     .               - FROC2/2D0*AY2MN2**2*AMYNI2*TT1(4)
     .               - froc2/2d0*tt1(5)*(aypni2+amyni)
        VVLA = VVLA  + GFMROC/REDM*VT1(4)
        VVQ12= VVQ12 + AMYNI2/16D0 * 3D0/RH2 * (GROC2*VV1(6)
     .               + 4D0*AMYPN*GFPROC*VT1(6)+8D0*AMYPN2*FROC2*TT1(6))
        FICV0= FICV0 + REDM*GROC2*(AYPNI2+AMYNI)/2D0*VV1(1)
        FICV1= FICV1 - REDM*GROC2*(AYPNI2+AMYNI)/2D0*VV1(4)*RH
        FICV2= FICV2 + REDM*GROC2*(AYPNI2+AMYNI)/2D0
     .               * (VV1(2)+2D0*VV1(4))
C       Second part of vector-vector potential (same as scalar)
        COUPL= GROC2*AMYMN2
        VVC  = VVC  - COUPL*(VV5(1) - VV5(2)*AMYNI/4D0)
        VVLS = VVLS - COUPL*VV5(4)*AMYNI/2D0
        VVQ12= VVQ12- COUPL*VV5(6)*AMYNI2/16D0 * 3D0/RH2
        FICV0= FICV0- REDM*COUPL*(AYPNI2/4D0-AMYNI)*VV5(1)
        FICV1= FICV1+ REDM*COUPL*(AYPNI2/4D0-AMYNI)*VV5(4)*RH
        FICV2= FICV2- REDM*COUPL*(AYPNI2/4D0-AMYNI)
     .              * (VV5(2)+2D0*VV5(4))
      ENDIF

      VC   = VSC         + VVC   + VDC
      VSS  =        VPVS + VVS
      VT   =        VPVT + VVT
      VLS  = VSLS        + VVLS  + VDLS
      VLSA = VSLA        + VVLA  + VDLA
      VQ12 = VSQ12       + VVQ12 + VDQ12

      FIC0 = FICS0        + FICV0 + FICD0
      FIC1 = FICS1        + FICV1 + FICD1
      FIC2 = FICS2        + FICV2 + FICD2

      RETURN
      END
***********************************************************************
      FUNCTION FDEXP(X)
      IMPLICIT REAL*8(A-Z)
      IF(X.LE.-100D0) THEN
        FDEXP=0D0
      ELSE
        FDEXP=DEXP(X)
      ENDIF
      RETURN
      END
***********************************************************************
      FUNCTION FDERFC(X)
      IMPLICIT REAL*8(A-Z)
      IF(X.GE. 10D0) THEN
        FDERFC=0D0
      ELSEIF(X.LE.-10D0) THEN
        FDERFC=2D0
      ELSE
        FDERFC=DERFC(X)
      ENDIF
      RETURN
      END
***********************************************************************
      SUBROUTINE DFI2DR(ND,RH,AM,ALAM,PLAM,D0FI,DFI)
**     Basic function  I_{2} = m phi^{0}_{C}  and its ND derivatives.
**     Each function is scaled by the charged pion mass squared (PIOMS2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION DFI(*)
      COMMON/AMSCAL/PIOMS2
      DATA ZLIM/1D-2/, RATLIM/1D0/, SQPI/1.7724538509055160273D0/

      AMP  =DSQRT(AM*AM+PLAM*PLAM)
      AMP2 =AMP*AMP
      AMP3 =AMP*AMP2
      RATM =AMP/ALAM
      RATM2=RATM*RATM
      RATM3=RATM*RATM2
      RATM5=RATM2*RATM3
      EXPM =FDEXP(AM*AM/(ALAM*ALAM))
      EXPL =FDEXP(-PLAM*PLAM/(ALAM*ALAM))
      Z=RH*ALAM
      IF(Z.GT.ZLIM .OR. RATM.GT.RATLIM) THEN
        XLAM = 0.5D0*RH*ALAM
        XLAM2= XLAM*XLAM
        ELAM = FDEXP(-XLAM2)
        X = AMP*RH
        X2= X*X
        X3= X*X2
        ERFCM = FDERFC(RATM-XLAM)
        ERFCP = FDERFC(RATM+XLAM)
        EPX = 0D0
        IF(ERFCP.NE.0) EPX = FDEXP(X)
        EMX = FDEXP(-X)
        PHI0C=EXPM*(EMX*ERFCM-EPX*ERFCP)/(2D0*X)
        PHI0L=(EXPM*(EMX*ERFCM*(1D0+X)-EPX*ERFCP*(1D0-X))
     .       - 4D0/SQPI*XLAM*ELAM*EXPL ) / (2D0*X3)
        PHI0T=(EXPM*(EMX*ERFCM*(1D0+X+X2/3D0)-EPX*ERFCP*(1D0-X+X2/3D0))
     .       - 8D0/SQPI*XLAM*(0.5D0+XLAM2/3D0)*ELAM*EXPL ) / (2D0*X3)
        DEL1C=-ELAM/(RATM3*2D0*SQPI) * EXPL
        DEL1L=-ELAM/(RATM5*4D0*SQPI) * EXPL
        DEL1T=-ELAM*XLAM2/(RATM5*6D0*SQPI) * EXPL
      ELSE
        A=RATM
        B=-FDERFC(A)*A*EXPM*SQPI
        A2=A*A
        A3=A2*A
        A4=A2*A2
        A5=A3*A2
        A6=A4*A2
        A8=A6*A2
        A10=A8*A2
        Z2=Z*Z
        Z4=Z2*Z2
        Z6=Z4*Z2
        Z8=Z6*Z2
        PHI0C = (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z8/5806080D0+
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z6/40320D0+
     .          (4*A4*B+4*A4-2*A2+3)*Z4/480D0+
     .          (2*A2*B+2*A2-1)*Z2/12D0+B+1D0
        PHI0C = PHI0C/(SQPI*A) * EXPL
        PHI0T = (32*A10*B+32*A10-16*A8+24*A6-60*A4+210*A2-945)
     .           *Z8/47900160D0+
     .          (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z6/362880D0+
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z4/5040D0+
     .          (4*A4*B+4*A4-2*A2+3)*Z2/180D0
        PHI0T = PHI0T/(SQPI*A3) * EXPL
        PHI0L =-(32*A10*B+32*A10-16*A8+24*A6-60*A4+210*A2-945)
     .           *Z8/127733760D0-
     .          (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z6/725760D0-
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z4/6720D0-
     .          (4*A4*B+4*A4-2*A2+3)*Z2/120D0-(2*A2*B+2*A2-1D0)/6D0
        PHI0L = PHI0L/(SQPI*A3) * EXPL
        DEL1C =-Z8/12288D0+Z6/768D0-Z4/64D0+Z2/8D0-0.5D0
        DEL1C = DEL1C/(SQPI*A3) * EXPL
        DEL1T = Z8/9216D0-Z6/768D0+Z4/96D0-Z2/24D0
        DEL1T = DEL1T/(SQPI*A5) * EXPL
        DEL1L =-Z8/24576D0+Z6/1536D0-Z4/128D0+Z2/16D0-0.25D0
        DEL1L = DEL1L/(SQPI*A5) * EXPL
      ENDIF

      D0FI  = AMP*PHI0C / PIOMS2
      DFI(1)=-AMP3*PHI0L*RH / PIOMS2
      IF(ND.EQ.1) RETURN

      PHI1C = PHI0C+DEL1C
      DFI(2)= AMP3*(PHI1C+2D0*PHI0L) / PIOMS2
      IF(ND.EQ.2) RETURN

      PHI1L = PHI0L+DEL1L
      DFI(3)=-AMP3*(AMP2*PHI1L*RH+6D0*PHI0T/RH) / PIOMS2
      IF(ND.EQ.3) RETURN

      DFI(4)= AMP3*AMP2*(PHI1C+3D0*(DEL1T-DEL1L))/PIOMS2 - 4D0/RH*DFI(3)
      IF(ND.EQ.4) RETURN

      END
***********************************************************************
      SUBROUTINE DFI4DR(ND,RH,AM,ALAM,D0FI,DFI)
**     Basic function  I_{4} = [-d/dm^{2}+1/Lambda^{2}] I_{2}(m,r)
**     and its ND derivatives.
**     Each function is scaled by the charged pion mass squared (PIOMS2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION DFI(*)
      COMMON/AMSCAL/PIOMS2
      DATA SQPI/1.7724538509055160273D0/

      RATM = AM/ALAM
      RATM2= RATM*RATM
      EXPM = FDEXP(RATM2)
      XLAM = 0.5D0*RH*ALAM
      X = AM*RH
      ERFCM = FDERFC(RATM-XLAM)
      ERFCP = FDERFC(RATM+XLAM)
      EPX = 0D0
      IF(ERFCP.NE.0) EPX = FDEXP(X)
      EMX = FDEXP(-X)

      D0FI  = 0.25D0*EXPM*(EMX*ERFCM+EPX*ERFCP)/AM / PIOMS2
      DFI(1)=-0.25D0*EXPM*(EMX*ERFCM-EPX*ERFCP) / PIOMS2
      IF(ND.EQ.1) RETURN

      AM2  = AM*AM
      XLAM2= XLAM*XLAM
      ELAM = FDEXP(-XLAM2)
      FORM= AM*ELAM/(2D0*SQPI*RATM * PIOMS2)
      DFI(2) = D0FI*AM2 - FORM
      IF(ND.EQ.2) RETURN

      DFI(3) = DFI(1)*AM2 + FORM*ALAM*XLAM
      IF(ND.EQ.3) RETURN

      END
***********************************************************************
      SUBROUTINE I2NYM(RH,AM,ALAM,PLAM,FI)
**     Basic function  I_{2} = m phi^{0}_{C}  of the Nijmegen model.
**     The function is scaled by the charged pion mass squared (PIOMS2)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/AMSCAL/PIOMS2
      DATA ZLIM/1D-2/, RATLIM/1D0/, SQPI/1.7724538509055160273D0/

      AMP =DSQRT(AM*AM+PLAM*PLAM)
      RATM=AMP/ALAM
      EXPM=FDEXP(AM*AM/(ALAM*ALAM))
      EXPL=FDEXP(-PLAM*PLAM/(ALAM*ALAM))
      Z=RH*ALAM
      IF(Z.GT.ZLIM .OR. RATM.GT.RATLIM) THEN
        XLAM = 0.5D0*RH*ALAM
        XLAM2= XLAM*XLAM
        ELAM = FDEXP(-XLAM2)
        X = AMP*RH
        ERFCM = FDERFC(RATM-XLAM)
        ERFCP = FDERFC(RATM+XLAM)
        EPX = 0D0
        IF(ERFCP.NE.0) EPX = FDEXP(X)
        EMX = FDEXP(-X)
        PHI0C=EXPM*(EMX*ERFCM-EPX*ERFCP)/(2D0*X)
      ELSE
        A=RATM
        B=-FDERFC(A)*A*EXPM*SQPI
        A2=A*A
        A4=A2*A2
        A6=A4*A2
        A8=A6*A2
        Z2=Z*Z
        Z4=Z2*Z2
        Z6=Z4*Z2
        Z8=Z6*Z2
        PHI0C = (16*A8*B+16*A8-8*A6+12*A4-30*A2+105)*Z8/5806080D0+
     .          (8*A6*B+8*A6-4*A4+6*A2-15)*Z6/40320D0+
     .          (4*A4*B+4*A4-2*A2+3)*Z4/480D0+
     .          (2*A2*B+2*A2-1)*Z2/12D0+B+1D0
        PHI0C = PHI0C/(SQPI*A) * EXPL
      ENDIF
      FI = AMP*PHI0C / PIOMS2

      RETURN
      END
***********************************************************************
      SUBROUTINE I2PHOT(RH,ALAM,PLAM,FI)
**     Basic function  I_{2} = lim (m-->0) m phi^{0}_{C}
**     (for massless photon), scaled by the charged pion mass squared
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/AMSCAL/PIOMS2

      AMP  = PLAM
      RATM = AMP/ALAM
      XLAM = 0.5D0*RH*ALAM
      X = AMP*RH
      ERFCM = FDERFC(RATM-XLAM)
      ERFCP = FDERFC(RATM+XLAM)
      EPX = 0D0
      IF(ERFCP.NE.0) EPX = FDEXP(X)
      EMX = FDEXP(-X)

      PHI0C = (EMX*ERFCM-EPX*ERFCP)/(2D0*RH)
      FI  = PHI0C / PIOMS2

      RETURN
      END
***********************************************************************
      SUBROUTINE PAIR(R,TYPE,ISO)
C*    Two-meson-exchange potential form pair terms (seagull vertices).
C*    Reference : Rijken & Stoks, Phys.Rev.C54 (1996) 2869.
C*    Mean pion mass everywhere
C*
C*    INPUT:  R       in fermi
C*    -----   TYPE    'PP', 'NN', 'NP', 'PN'
C*            ISO     total isospin
C*
C*    OUTPUT: COMMON/VPAIR/V1C,V1S,V1T,V1L, V2C,V2S,V2T,V2L
C*    ------
C*    Pair interactions considered:
C*        J(PC)=0(++):   (\pi\pi)_{0}    (\sigma\sigma)     (\pi\eta)
C*        J(PC)=1(--):   (\pi\pi)_{1} [g and f couplings]
C*        J(PC)=1(++):   (\pi\rho)_{1}   (\pi\sigma)    (\pi Pomeron)
C*
C*       IPV = 1: 1/M corrections from vertex functions
C*       INA = 1: k1.k2/M corrections from nonadiabatic expansion
C*           = 2: q.(k1-k2)/M corrections from nonadiabatic expansion
C*
C*    Note that the total r-dependent function is divided by PIOMS4 !!
C-----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER TYPE*2
      REAL*8 FINT(8),F2PI(3),I2POM,F2ET(2),F2RO(1),F2ROT(1),
     .       F2SI1(2),F2SI2(2),F2SI(2),F2RO1(1),F2RO2(1)
      COMMON/VPAIR/ V1C,V1S,V1T,V1L, V2C,V2S,V2T,V2L
      COMMON/GPAIR/ GPIETA,GPIETP,GSISI, GPISIG,GPIPOM,
     .              GPIPI0,GPIPI1,FPIPI1,GPIRO0,GPIRO1,GPIOM
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/CUTTEN/ALAMT1,ALAVT1,ALAMT0,ALAVT0
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/KAPPAV/AKRO,AKOM,AKFI
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4,
     .              GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4
      COMMON/POMFUN/I2POM,DI2POM,DDI2POM
      COMMON/I2NADI/CSSNA,CPPNA,SOPPNA,CPENA,SOPENA,
     .              SPSNA,TPSNA,SPONA,TPONA
      COMMON/SCALIN/AMT,AMPV,AMT24
      COMMON/AMSCAL/PIOMS2,PIOMS
      COMMON/HBARC/ HBC
      DATA SQPI/1.7724538509055160273D0/, TWONPI/0.636619772367581343D0/
      DATA NP/12/, IPV/1/, INA/2/

      FISO=4D0*ISO-3D0
      RH = R/HBC
      RH2= RH*RH

C*    General functions:
      CALL DFI2DR(3,RH,PIM,ALAMP1,0D0,F0PI,F2PI)
      CALL DFI2DR(2,RH,AME1,ALAMS0,0D0,F0SI1,F2SI1)
      CALL DFI2DR(2,RH,AME2,ALAMS0,0D0,F0SI2,F2SI2)
      F0SI   =AEPS*F0SI1   +BEPS*F0SI2
      F2SI(1)=AEPS*F2SI1(1)+BEPS*F2SI2(1)
      F2SI(2)=AEPS*F2SI1(2)+BEPS*F2SI2(2)
ce    CALL DFI2DR(2,RH,AMETA,ALAMP0,0D0,F0ETA,F2ET)
      CALL DFI2DR(1,RH,AMR1,ALAMV1,0D0,F0RO1,F2RO1)
      CALL DFI2DR(1,RH,AMR2,ALAMV1,0D0,F0RO2,F2RO2)
      F0RO   =ARO*F0RO1   +BRO*F0RO2
      F2RO(1)=ARO*F2RO1(1)+BRO*F2RO2(1)
      CALL DFI2DR(1,RH,AMR1,ALAVT1,0D0,F0RO1,F2RO1)
      CALL DFI2DR(1,RH,AMR2,ALAVT1,0D0,F0RO2,F2RO2)
      F2ROT(1)=ARO*F2RO1(1)+BRO*F2RO2(1)
cp    XPOM=AMPOM*RH
cp    XP2=XPOM*XPOM
cp    I2POM=0D0
cp    IF(XP2.LT.170.) I2POM = 4D0*AMPOM/SQPI * (AMPOM/AMT)**2
cp   .                      * FDEXP(-XP2)/PIOMS2
cp    DI2POM =-2D0*AMPOM2*RH*I2POM
cp    DDI2POM=-2D0*AMPOM2*(1D0-2D0*XP2)*I2POM
      ALAM2= ALAMP1*ALAMP1
      PPI0 = (ALAM2/PIOMS2)*ALAMP1/(2D0*SQPI)*FDEXP(-ALAM2*RH2/4D0)

      CSSNA=F2SI(1)**2
      CPPNA=F2PI(2)**2+2D0/RH2*F2PI(1)**2
      SOPPNA=2D0/RH2*F2PI(1)**2
ce    CPENA=F2PI(2)*F2ET(2)+2D0/RH2*F2PI(1)*F2ET(1)
ce    SOPENA=2D0/RH2*F2PI(1)*F2ET(1)
      SPSNA=(F2PI(3)+2D0*(F2PI(2)-F2PI(1)/RH)/RH)*F2SI(1)
     .     -(F2PI(2)*F2SI(2)+2D0*F2PI(1)*F2SI(1)/RH2)
      TPSNA=(F2PI(3)-(F2PI(2)-F2PI(1)/RH)/RH)*F2SI(1)
     .     -(F2PI(2)*F2SI(2)-F2PI(1)*F2SI(1)/RH2)
cp    SPONA=(F2PI(3)+2D0*(F2PI(2)-2D0*F2PI(1)/RH)/RH)*DI2POM
cp   .      -F2PI(2)*DDI2POM
cp    TPONA=(F2PI(3)-(F2PI(2)-2D0*F2PI(1)/RH)/RH)*DI2POM
cp   .      -F2PI(2)*DDI2POM

C*    Lambda integrals:
      CALL PAIRINT(RH,1,NP,FINT)
      B11PI1=PIM*FINT(1)*TWONPI
      B11PI0=PIM*FINT(2)*TWONPI
      FINTPI=PIM*FINT(3)*TWONPI
      B00PI1=PIM*FINT(4)*TWONPI
      SISINA=PIM*FINT(5)*TWONPI
      B11SI =PIM*FINT(6)*TWONPI
      B11RO =PIM*FINT(7)*TWONPI
ce    B11ETA=PIM*FINT(8)*TWONPI
      IF(IPV.NE.0 .OR. INA.NE.0) THEN
        CALL PAIRINT(RH,2,NP,FINT)
        PI0PVC =PIM*FINT(1)*TWONPI
        PI0PVL =PIM*FINT(2)*TWONPI
        PI0NAC =PIM*FINT(3)*TWONPI
        PI0NAL =PIM*FINT(4)*TWONPI
ce      PIETPVC=PIM*FINT(5)*TWONPI
ce      PIETPVL=PIM*FINT(6)*TWONPI
ce      PIETNAC=PIM*FINT(7)*TWONPI
ce      PIETNAL=PIM*FINT(8)*TWONPI
        CALL PAIRINT(RH,3,NP,FINT)
        PISIPVS=PIM*FINT(1)*TWONPI
        PISIPVT=PIM*FINT(2)*TWONPI
        PISINAS=PIM*FINT(3)*TWONPI
        PISINAT=PIM*FINT(4)*TWONPI
cp      PIPOPVS=PIM*FINT(5)*TWONPI
cp      PIPOPVT=PIM*FINT(6)*TWONPI
cp      PIPONAS=PIM*FINT(7)*TWONPI
cp      PIPONAT=PIM*FINT(8)*TWONPI
      ENDIF
      CALL PAIRINT(RH,4,NP,FINT)
      B11SIS=PIM*FINT(1)*TWONPI
      B11SIT=PIM*FINT(2)*TWONPI
cp    B11POS=PIM*FINT(3)*TWONPI
cp    B11POT=PIM*FINT(4)*TWONPI

C***  One-pair potentials:
C***  --------------------

C*    Scalar 0++
      FF1=GPIPI0*PIOMS*FPI*FPI
      FF2=GSISI*PIOMS*GEP*GEP*PIOMS2
      FF3=FISO*GPIETA*PIOMS*FPI*FETA
      VC0PP = 6D0*FF1*F2PI(1)*F2PI(1) + 2D0*FF2*F0SI*F0SI
ce   .      + 2D0*FF3*F2PI(1)*F2ET(1)
      VL0PP = 0D0
      IF(IPV.EQ.1) THEN
ce      VC0PP = VC0PP + 6D0*FF1/AMT*PI0PVC + FF3/AMT*PIETPVC
ce      VL0PP = 6D0*FF1/AMT*PI0PVL + FF3/AMT*PIETPVL
        VC0PP = VC0PP + 6D0*FF1/AMT*PI0PVC
        VL0PP = 6D0*FF1/AMT*PI0PVL
      ENDIF
      IF(INA.GE.1) THEN
ce      VC0PP = VC0PP - (3D0*FF1*PI0NAC+FF3*PIETNAC-FF2*SISINA)/AMT
ce      IF(INA.EQ.2) VL0PP = VL0PP - (3D0*FF1*PI0NAL+FF3*PIETNAL)/AMT
        VC0PP = VC0PP - (3D0*FF1*PI0NAC-FF2*SISINA)/AMT
        IF(INA.EQ.2) VL0PP = VL0PP - 3D0*FF1*PI0NAL/AMT
      ENDIF

C*    Vector 1--
      FF1=FISO*GPIPI1*FPI*FPI
      FF2=FISO*4D0/3D0*(GPIPI1+FPIPI1)/AMT*FPI*FPI
      VC1MM = FF1*4D0*B11PI1
      VS1MM =-FF2*F2PI(1)/RH*(F2PI(1)/RH+2D0*F2PI(2))
      VT1MM =-FF2*F2PI(1)/RH*(F2PI(1)/RH-F2PI(2))
      VL1MM =-FF1*4D0/AMT*F2PI(1)**2/RH2
      IF(IPV.EQ.1) THEN
        VC1MM = VC1MM - FF1/AMT*2D0*(F2PI(2)+2D0*F2PI(1)/RH)*PPI0
        VL1MM = VL1MM - FF1/AMT*4D0/RH*F2PI(1)*PPI0
      ENDIF
      IF(INA.GE.1) THEN
        VC1MM = VC1MM - FF1*2D0/AMT*(F2PI(2)**2+2D0/RH2*F2PI(1)**2)
        IF(INA.EQ.2) VL1MM = VL1MM - FF1*4D0/AMT*F2PI(1)**2/RH2
      ENDIF

C*    Axial 1++
      FF1=FISO/3D0*GPISIG*GEP*PIOMS*FPI
      FF2=FISO/3D0*GPIPOM*GPOM*PIOMS*FPI
      VS1PP = FF1*2D0*((F2PI(2)+2D0*F2PI(1)/RH)*F0SI
     .                 -F2PI(1)*F2SI(1))
cp   .      - FF2*2D0*((F2PI(2)+2D0*F2PI(1)/RH)*I2POM-F2PI(1)*DI2POM)
      VT1PP = FF1*2D0*((F2PI(2)-F2PI(1)/RH)*F0SI
     .                 -F2PI(1)*F2SI(1))
cp   .      - FF2*2D0*((F2PI(2)-F2PI(1)/RH)*I2POM-F2PI(1)*DI2POM)
      IF(IPV.EQ.1) THEN
cp      VS1PP = VS1PP + (FF1*PISIPVS-FF2*PIPOPVS)/AMT
cp      VT1PP = VT1PP + (FF1*PISIPVT-FF2*PIPOPVT)/AMT
        VS1PP = VS1PP + FF1*PISIPVS/AMT
        VT1PP = VT1PP + FF1*PISIPVT/AMT
      ENDIF
      IF(INA.GE.1) THEN
cp      VS1PP = VS1PP - (FF1*PISINAS-FF2*PIPONAS)/AMT
cp      VT1PP = VT1PP - (FF1*PISINAT-FF2*PIPONAT)/AMT
        VS1PP = VS1PP - FF1*PISINAS/AMT
        VT1PP = VT1PP - FF1*PISINAT/AMT
      ENDIF
      FF3=FISO*2D0/3D0*GPIRO1*PIOMS2*FPI*GRO/AMT
      VS1PP = VS1PP - FF3*( (F2PI(2)+2D0*F2PI(1)/RH)*F0RO
     .                     -2D0*F2PI(1)*(F2RO(1)+AKRO*F2ROT(1)) )
      VT1PP = VT1PP - FF3*( (F2PI(2)-F2PI(1)/RH)*F0RO
     .                     +F2PI(1)*(F2RO(1)+AKRO*F2ROT(1)) )

      V1C = VC0PP+VC1MM
      V1S =       VS1MM+VS1PP
      V1T =       VT1MM+VT1PP
      V1L = VL0PP+VL1MM

C***  Two-pair potentials:
C***  --------------------

C*    Scalar 0++
      WC0PP = -3D0*GPIPI0*GPIPI0*PIOMS2*B11PI0
     .        -GSISI*GSISI*PIOMS2*B11SI
ce   .        -0.5D0*FISO*GPIETA*GPIETA*PIOMS2*B11ETA

C*    Vector 1--
      WC1MM =-FISO*GPIPI1*GPIPI1*(PPI0*FINTPI-2D0*B00PI1)

C*    Axial 1++
cp    WS1PP = FISO/6D0*(GPISIG*GPISIG*B11SIS - GPIPOM*GPIPOM*B11POS)
cp    WT1PP = FISO/6D0*(GPISIG*GPISIG*B11SIT - GPIPOM*GPIPOM*B11POT)
      WS1PP = FISO/6D0*GPISIG*GPISIG*B11SIS
      WT1PP = FISO/6D0*GPISIG*GPISIG*B11SIT
      WS1PP = WS1PP - FISO*GPIRO1*GPIRO1*PIOMS2*B11RO

      V2C = WC0PP+WC1MM
      V2S = WS1PP
      V2T = WT1PP
      V2L = 0D0

      RETURN
      END
************************************************************************
      SUBROUTINE PAIRINT(RH,NC,NP,SUM)
**      Calculates integrals for pair terms using Gauss integration
**      on (0,infinity) to (-1,1) for hyperbolic mapping with NP points
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION SUM(8),FUNCW(8),X(20),W(20),NDIM(4)
cpce  DATA ICAL/0/, NDIM/8,8,8,4/
      DATA ICAL/0/, NDIM/7,4,4,2/
      SAVE X,W

      IF(ICAL.EQ.0) THEN
        CALL GAUSS(NP,X,W,1)
        DO 1 I=1,NP
          W(I)=2D0*W(I)/(1D0-X(I))**2
          X(I)=(1D0+X(I))/(1D0-X(I))
    1   CONTINUE
        ICAL=1
      ENDIF

      KDIM=NDIM(NC)
      DO 10 K=1,KDIM
   10   SUM(K)=0D0
      DO 11 IY=1,NP
        CALL FUNCP(FUNCW,X(IY),RH,NC)
        DO 12 K=1,KDIM
   12     SUM(K)=SUM(K)+FUNCW(K)*W(IY)
   11 CONTINUE

      RETURN
      END
************************************************************************
      SUBROUTINE FUNCP(FINT,Y,RH,NC)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 FINT(8),F2PI(2),F3PI(3),I2POM,F2SI(2),F2ET(2),
     .       F1SI1(1),F1SI2(1),F2SI1(2),F2SI2(2)
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4,
     .              GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4
      COMMON/POMFUN/I2POM,DI2POM,DDI2POM
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/AMSCAL/PIOMS2,PIOMS
      COMMON/I2NADI/CSSNA,CPPNA,SOPPNA,CPENA,SOPENA,
     .              SPSNA,TPSNA,SPONA,TPONA

      RH2=RH*RH
      PLAM=Y*PIM
      PLAM2=PLAM*PLAM
      IF(NC.EQ.1) THEN
        CALL DFI2DR(2,RH,PIM,ALAMP1,PLAM,F0PI,F2PI)
        CALL DFI2DR(1,RH,AME1,ALAMS0,PLAM,F0SI1,F1SI1)
        CALL DFI2DR(1,RH,AME2,ALAMS0,PLAM,F0SI2,F1SI2)
        F0SI   =AEPS*F0SI1   +BEPS*F0SI2
        F2SI(1)=AEPS*F1SI1(1)+BEPS*F1SI2(1)
        CALL I2NYM(RH,AMR1,ALAMV1,PLAM,F0RO1)
        CALL I2NYM(RH,AMR2,ALAMV1,PLAM,F0RO2)
        F0RO=ARO*F0RO1+BRO*F0RO2
ce      CALL I2NYM(RH,AMETA,ALAMP0,PLAM,F0ETA)
**        -(k1.k2)B_{1,1}(pi,pi)_{1}      [one pair]
        FINT(1)=F2PI(1)*F2PI(1)
**        B_{1,1}(pi,pi)_{0}              [two pair]
        FINT(2)=F0PI*F0PI
**        F(pi)_{1}, B_{0,0}(pi,pi)_{1}   [two pair]
        FINT(3)=F0PI
        FINT(4)=PLAM2*F0PI*F0PI
**        (k_{1}.k_{2}) [ I I - F F ](sigma,sigma)  [one pair]
        FINT(5)=-(CSSNA-F2SI(1)**2)/PLAM2
**        B_{1,1}(sigma,sigma)            [two pair]
        FINT(6)=F0SI*F0SI
**        B_{1,1}(pi,rho)                 [two pair]
        FINT(7)=F0PI*F0RO
**        B_{1,1}(pi,eta)                 [two pair]
ce      FINT(8)=F0PI*F0ETA
      ELSEIF(NC.EQ.2) THEN
**      Pseudovector and nonadiabatic for pi-pi_{0}, pi-eta:
        CALL DFI2DR(2,RH,PIM,ALAMP1,PLAM,F0PI,F2PI)
ce      CALL DFI2DR(2,RH,AMETA,ALAMP0,PLAM,F0ET,F2ET)
**        1/2(k_{1}^{2}+k_{2}^{2}) B_{1,1}(pi,pi)
        FINT(1)=-(F2PI(2)+2D0*F2PI(1)/RH)*F0PI
**        -i/2(sigma1+sigma2).(k_{1}+k_{2})xq B_{1,1}(pi,pi)
        FINT(2)=-2D0/RH*F2PI(1)*F0PI
**        (k_{1}.k_{2})^{2} [ I I - F F ](pi,pi)
        FINT(3)=(CPPNA-F2PI(2)**2-2D0*(F2PI(1)/RH)**2)/PLAM2
**        i/2 q.(k1-k2)(sigma1+sigma2).(k1xk2) [ I I - F F ](pi,pi)
        FINT(4)=(SOPPNA-2D0*(F2PI(1)/RH)**2)/PLAM2
**        (k_{1}^{2}+k_{2}^{2}) B_{1,1}(pi,eta)
ce      FINT(5)=-(F2PI(2)+2D0*F2PI(1)/RH)*F0ET
ce   .          -F0PI*(F2ET(2)+2D0*F2ET(1)/RH)
**        -i(sigma1+sigma2).(k_{1}+k_{2})xq B_{1,1}(pi,eta)
ce      FINT(6)=-2D0/RH*(F2PI(1)*F0ET+F0PI*F2ET(1))
**        (k_{1}.k_{2})^{2} [ I I - F F ](pi,eta)
ce      FINT(7)=(CPENA-F2PI(2)*F2ET(2)-2D0/RH2*F2PI(1)*F2ET(1))/PLAM2
**        i/2 q.(k1-k2)(sigma1+sigma2).(k1xk2) [ I I - F F ](pi,eta)
ce      FINT(8)=(SOPENA-2D0/RH2*F2PI(1)*F2ET(1))/PLAM2
      ELSEIF(NC.EQ.3) THEN
**      Pseudovector and nonadiabatic for pi-sigma, pi-pomeron:
        CALL DFI2DR(3,RH,PIM,ALAMP1,PLAM,F0PI,F3PI)
        CALL DFI2DR(2,RH,AME1,ALAMS0,PLAM,F0SI1,F2SI1)
        CALL DFI2DR(2,RH,AME2,ALAMS0,PLAM,F0SI2,F2SI2)
        F2SI(1)=AEPS*F2SI1(1)+BEPS*F2SI2(1)
        F2SI(2)=AEPS*F2SI1(2)+BEPS*F2SI2(2)
**        3/2[si1.k1si2.k2+si1.k2si2.k1-2si1.k2si2.k2]B_{1,1}(pi,sig)
        FINT(1)=F0PI*(F2SI(2)+2D0*F2SI(1)/RH)-F3PI(1)*F2SI(1)
        FINT(2)=F0PI*(F2SI(2)-F2SI(1)/RH)-F3PI(1)*F2SI(1)
**        3/2k1.k2[2si1.k1si2.k1-si1.k1si2.k2-si1.k2si2.k1](II-FF)(pi,sig)
        SPSF=(F3PI(3)+2D0*(F3PI(2)-F3PI(1)/RH)/RH)*F2SI(1)
     .      -(F3PI(2)*F2SI(2)+2D0*F3PI(1)*F2SI(1)/RH2)
        TPSF=(F3PI(3)-(F3PI(2)-F3PI(1)/RH)/RH)*F2SI(1)
     .      -(F3PI(2)*F2SI(2)-F3PI(1)*F2SI(1)/RH2)
        FINT(3)=(SPSNA-SPSF) / PLAM2
        FINT(4)=(TPSNA-TPSF) / PLAM2
cp      FACPOM=FDEXP(-PLAM2/(4D0*AMPOM2))
**        3/2[si1.k1si2.k2+si1.k2si2.k1-2si1.k2si2.k2]B_{1,1}(pi,Pom)
cp      FINT(5)=FACPOM * (F0PI*(DDI2POM+2D0*DI2POM/RH)-F3PI(1)*DI2POM)
cp      FINT(6)=FACPOM * (F0PI*(DDI2POM-DI2POM/RH)-F3PI(1)*DI2POM)
**        3/2k1.k2[2si1.k1si2.k1-si1.k1si2.k2-si1.k2si2.k1](II-FF)(pi,Pom)
cp      FINT(7)=(SPONA-((F3PI(3)+2D0*(F3PI(2)-2D0*F3PI(1)/RH)/RH)
cp   .                  *DI2POM-F3PI(2)*DDI2POM)*FACPOM) / PLAM2
cp      FINT(8)=(TPONA-((F3PI(3)-(F3PI(2)-2D0*F3PI(1)/RH)/RH)
cp   .                  *DI2POM-F3PI(2)*DDI2POM)*FACPOM) / PLAM2
      ELSEIF(NC.EQ.4) THEN
        CALL DFI2DR(2,RH,PIM,ALAMP1,PLAM,F0PI,F2PI)
        CALL DFI2DR(2,RH,AME1,ALAMS0,PLAM,F0SI1,F2SI1)
        CALL DFI2DR(2,RH,AME2,ALAMS0,PLAM,F0SI2,F2SI2)
        F0SI   =AEPS*F0SI1   +BEPS*F0SI2
        F2SI(1)=AEPS*F2SI1(1)+BEPS*F2SI2(1)
        F2SI(2)=AEPS*F2SI1(2)+BEPS*F2SI2(2)
**        3 sigma1.(k1-k2)sigma2.(k1-k2)B_{1,1}(pi,sigma)
        FINT(1)=(F2PI(2)+2D0*F2PI(1)/RH)*F0SI
     .         - 2D0*F2PI(1)*F2SI(1)
     .         + F0PI*(F2SI(2)+2D0*F2SI(1)/RH)
        FINT(2)=(F2PI(2)-F2PI(1)/RH)*F0SI
     .         - 2D0*F2PI(1)*F2SI(1)
     .         + F0PI*(F2SI(2)-F2SI(1)/RH)
**        3 sigma1.(k1-k2)sigma2.(k1-k2)B_{1,1}(pi,Pomeron)
cp      FACPOM=FDEXP(-PLAM2/(4D0*AMPOM2))
cp      FINT(3)=( (F2PI(2)+2D0*F2PI(1)/RH)*I2POM - 2D0*F2PI(1)*DI2POM
cp   .           + F0PI*(DDI2POM+2D0*DI2POM/RH) ) * FACPOM
cp      FINT(4)=( (F2PI(2)-F2PI(1)/RH)*I2POM - 2D0*F2PI(1)*DI2POM
cp   .           + F0PI*(DDI2POM-DI2POM/RH) ) * FACPOM
      ENDIF

      RETURN
      END
************************************************************************
      SUBROUTINE PIMES(R,TYPE,ISO)
C*    Pion-meson-exchange potentials.
C*    Reference : Rijken & Stoks, Phys.Rev.C54 (1996) 2851.
C*    Mean pion mass everywhere
C*
C*    INPUT:  R       in fermi
C*    -----   TYPE    'PP', 'NN', 'NP', 'PN'
C*            ISO     total isospin
C*
C*    OUTPUT: COMMON/VPIMES/VCPI,VSPI,VTPI,VLPI, VCPM,VSPM,VTPM,VLPM
C*    ------
C*    Note that the total r-dependent function is divided by PIOMS4 !!
C*    Standard: leading order Brueckner-Watson (parallel and crossed)
C*              and Taketani-Machida-Ohnuma (parallel)
C*    IPV = 1 : 1/M contributions from pseudovector correction
C*        =-1 : pseudoscalar coupling ==> sign change for INA=3
C*              this has to be set also in OM2VK and OM2SCA
C*    INA = 1 : k1.k2/M contributions from nonadiabatic expansion
C*        = 2 : q.(k1-k2)/M contributions from nonadiabatic expansion
C*        = 3 : 1/M off-shell contribution in TMO diagrams
C*              this has to be set also in OM2VK and OM2SCA
C*    IM2 = 1 : 1/M^2 contributions from (1+f/g) in vector mesons
C*        = 2 : 1/M^2 contributions remaining in vector mesons
C*        = 3 : 1/M^2 contributions in scalar and diffractive mesons
C-----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER TYPE*2
      REAL*8 FINT(11)
      REAL*8 I2PI,I2ET,I2EP, I2RO,I2OM,I2FI, I2A0,I2SI,I2F0,
     .       I4PI,I4ET,I4EP, I4RO,I4OM,I4FI, I4A0,I4SI,I4F0,
     .       I2POM,I2A2D,I4POM,I4A2D
      REAL*8 D2PI(4),D2ET(3),D2EP(3), D4PI(3),D4ET(3),D4EP(3),
     .       D2RO(2),D2OM(2),D2FI(2), D4RO(1),D4OM(1),D4FI(1),
     .       D2A0(2),D2SI(1),D2F0(1), D4A0(1),D4SI(1),D4F0(1),
     .       F2RO1(2),F2RO2(2),F2SI1(1),F2SI2(1),
     .       F4RO1(1),F4RO2(1),F4SI1(1),F4SI2(1)
      REAL*8 DROK(1),DRO2K(2),DROK2(2),DOMK(1),DFIK(1),DOMK2(2),DFIK2(2)
      common/hoiho/scal
      COMMON/VPIMES/VCPI,VSPI,VTPI,VLPI, VCPM,VSPM,VTPM,VLPM
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/KAPPAV/AKRO,AKOM,AKFI
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4,
     .              GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4
      COMMON/SCALIN/AMT,AMPV,AMT24
      COMMON/AMSCAL/PIOMS2,PIOMS
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/CUTTEN/ALAMT1,ALAVT1,ALAMT0,ALAVT0
      COMMON/IITRMS/PPIIC,PPIIS,PPIIT,
     .              ETIIS,ETIIT, EPIIS,EPIIT,
     .              A0IIS,A0IIT, A2IIS,A2IIT,I2A2D,DI2A2D,DDI2A2D,
     .              RIIEES,RIIEET,RIIEMT,RIIEML,RIIMMC,RIIMMS,RIIMMT,
     .              OIIEMT,OIIMMS,OIIMMT, FIIEMT,FIIMMS,FIIMMT,
     .              ROIIM2S,ROIIM2T,ROIIM2L, ROIIS,ROIIT,ROIIS1,ROIIT1,
     .              OMIIS1,OMIIT1, FIIIS1,FIIIT1,
     .              A0IIM2S,A0IIM2T,A0IIM2L, A2IIM2S,A2IIM2T,A2IIM2L
      COMMON/HBARC/ HBC
      DATA SQPI/1.7724538509055160273D0/, TWONPI/0.636619772367581343D0/
      DATA SQ2/1.414213562373095D0/
      DATA NP/12/, IPV/1/, INA/3/, IM2/3/

**      isospin factors
      FISO=4D0*ISO-3D0
      FPA=3D0-2D0*FISO
      FCR=3D0+2D0*FISO

      RH = R/HBC
**      I2(m,r) functions
      CALL DFI2DR(4,RH,PIM,  ALAMP1,0D0,I2PI, D2PI)
      CALL DFI2DR(3,RH,AMETA,ALAMP0,0D0,I2ET, D2ET)
      CALL DFI2DR(3,RH,AMETP,ALAMP0,0D0,I2EP, D2EP)
      CALL DFI2DR(2,RH,A0M,  ALAMS1,0D0,I2A0, D2A0)
      CALL DFI2DR(1,RH,AME1, ALAMS0,0D0,F0SI1,F2SI1)
      CALL DFI2DR(1,RH,AME2, ALAMS0,0D0,F0SI2,F2SI2)
      CALL DFI2DR(1,RH,AMF0, ALAMS0,0D0,I2F0, D2F0)
      CALL DFI2DR(2,RH,AMR1, ALAMV1,0D0,F0RO1,F2RO1)
      CALL DFI2DR(2,RH,AMR2, ALAMV1,0D0,F0RO2,F2RO2)
      CALL DFI2DR(2,RH,AMOM, ALAMV0,0D0,I2OM, D2OM)
      CALL DFI2DR(2,RH,AMFI, ALAMV0,0D0,I2FI, D2FI)
      I2SI   =AEPS*F0SI1   +BEPS*F0SI2
      D2SI(1)=AEPS*F2SI1(1)+BEPS*F2SI2(1)
      I2RO   =ARO*F0RO1   +BRO*F0RO2
      D2RO(1)=ARO*F2RO1(1)+BRO*F2RO2(1)
      D2RO(2)=ARO*F2RO1(2)+BRO*F2RO2(2)
      XPOM=AMPOM*RH
      XP2=XPOM*XPOM
      I2POM=0D0
      IF(XP2.LT.170.) I2POM=FDEXP(-XP2)
      DI2POM =-2D0*AMPOM2*RH*I2POM
      DDI2POM=-2D0*AMPOM2*(1D0-2D0*XP2)*I2POM
cp    XA2=AMA2D*RH
cp    XP2=XA2*XA2
cp    I2A2D=0D0
cp    IF(XP2.LT.170.) I2A2D=FDEXP(-XP2)
cp    DI2A2D =-2D0*AMA2D2*RH*I2A2D
cp    DDI2A2D=-2D0*AMA2D2*(1D0-2D0*XP2)*I2A2D

      IF(IPV.NE.0 .OR. INA.NE.0) THEN
        PIM2 =PIM*PIM
        AMET2=AMETA*AMETA
        AMEP2=AMETP*AMETP
**      I0'(m,r) functions for 1/M corrections in pseudoscalars
        XLAM =0.5D0*ALAMP1*RH
        XLAM2=XLAM*XLAM
        ALAM2=ALAMP1*ALAMP1
        DI0PI=-(ALAM2/PIOMS2)*ALAM2*XLAM/(2D0*SQPI)*FDEXP(-XLAM2)
        XLAM =0.5D0*ALAMP0*RH
        XLAM2=XLAM*XLAM
        ALAM2=ALAMP0*ALAMP0
        DI0PS=-(ALAM2/PIOMS2)*ALAM2*XLAM/(2D0*SQPI)*FDEXP(-XLAM2)
      ENDIF

      IF(INA.NE.0) THEN
**      I4(m,r) functions for nonadiabatic corrections
        CALL DFI4DR(3,RH,PIM,  ALAMP1,I4PI, D4PI)
        CALL DFI4DR(3,RH,AMETA,ALAMP0,I4ET, D4ET)
        CALL DFI4DR(3,RH,AMETP,ALAMP0,I4EP, D4EP)
        CALL DFI4DR(1,RH,A0M,  ALAMS1,I4A0, D4A0)
        CALL DFI4DR(1,RH,AME1, ALAMS0,F0SI1,F4SI1)
        CALL DFI4DR(1,RH,AME2, ALAMS0,F0SI2,F4SI2)
        CALL DFI4DR(1,RH,AMF0, ALAMS0,I4F0, D4F0)
        CALL DFI4DR(1,RH,AMR1, ALAMV1,F0RO1,F4RO1)
        CALL DFI4DR(1,RH,AMR2, ALAMV1,F0RO2,F4RO2)
        CALL DFI4DR(1,RH,AMOM, ALAMV0,I4OM, D4OM)
        CALL DFI4DR(1,RH,AMFI, ALAMV0,I4FI, D4FI)
        I4SI   =AEPS*F0SI1   +BEPS*F0SI2
        D4SI(1)=AEPS*F4SI1(1)+BEPS*F4SI2(1)
        I4RO   =ARO*F0RO1   +BRO*F0RO2
        D4RO(1)=ARO*F4RO1(1)+BRO*F4RO2(1)
        I4POM =I2POM /(4D0*AMPOM2)
        DI4POM=DI2POM/(4D0*AMPOM2)
cp      I4A2D =I2A2D /(4D0*AMA2D2)
cp      DI4A2D=DI2A2D/(4D0*AMA2D2)
        IF(INA.GE.3) THEN
**        Common combination of pion function
          PIS=D2PI(2)+2D0*D2PI(1)/RH
          PIT=D2PI(2)-D2PI(1)/RH
        ENDIF
      ENDIF

      CALL OPSEU(RH,D2PI(1),D2PI(2),D2PI(1),D2PI(2),PPIIC,PPIIS,PPIIT)
      CALL TMEINT(RH,1,NP,FINT)
      BPIPIC = PIM*FINT(1)*TWONPI
      BPIPIS = PIM*FINT(2)*TWONPI
      BPIPIT = PIM*FINT(3)*TWONPI
      CALL OPSEU(RH,D2PI(1),D2PI(2),D2ET(1),D2ET(2),ETIIC,ETIIS,ETIIT)
      CALL OPSEU(RH,D2PI(1),D2PI(2),D2EP(1),D2EP(2),EPIIC,EPIIS,EPIIT)
      CALL OSCL(RH,D2PI(1),D2PI(2),I2A0,A0IIS,A0IIT)
      CALL OVC1(RH,D2PI(1),D2PI(2),I2RO,RIIEES,RIIEET)
cp    CALL OSCL(RH,D2PI(1),D2PI(2),I2A2D,A2IIS,A2IIT)
      CALL TMEINT(RH,2,NP,FINT)
      BPIETS = PIM*FINT(1) *TWONPI
      BPIETT = PIM*FINT(2) *TWONPI
      BPIEPS = PIM*FINT(3) *TWONPI
      BPIEPT = PIM*FINT(4) *TWONPI
      BPIA0S = PIM*FINT(5) *TWONPI
      BPIA0T = PIM*FINT(6) *TWONPI
      BROEES = PIM*FINT(7) *TWONPI
      BROEET = PIM*FINT(8) *TWONPI
cp    BPIA2S = PIM*FINT(9) *TWONPI
cp    BPIA2T = PIM*FINT(10)*TWONPI

      IF(IM2.GE.1) THEN
**      Vector-meson functions for 1+f/g and (1+f/g)**2
        CALL SPECV(1,RH,0D0,D2RO,DROK,DROK2,AKRO)
        CALL SPECV(2,RH,0D0,D2OM,DOMK,DOMK2,AKOM)
        CALL SPECV(3,RH,0D0,D2FI,DFIK,DFIK2,AKFI)
**      Operators for 1/M^2 corrections for (1+f/g) and (1+f/g)**2
        CALL OVC2(RH,D2PI(1),D2PI(2),DROK(1),DROK2(1),DROK2(2),
     .            RIIEMT,RIIMMS,RIIMMT)
        CALL OVC3(RH,D2PI(1),D2PI(2),DROK(1),DROK2(1),DROK2(2),
     .            RIIEML,RIIMMC)
        CALL OVC2(RH,D2PI(1),D2PI(2),DOMK(1),DOMK2(1),DOMK2(2),
     .            OIIEMT,OIIMMS,OIIMMT)
        CALL OVC2(RH,D2PI(1),D2PI(2),DFIK(1),DFIK2(1),DFIK2(2),
     .            FIIEMT,FIIMMS,FIIMMT)
        CALL TMEINT(RH,3,NP,FINT)
        BROEMT = PIM*FINT(1) *TWONPI/AMT24
        BROEML = PIM*FINT(2) *TWONPI/AMT24
        BROMMC = PIM*FINT(3) *TWONPI/AMT24
        BROMMS = PIM*FINT(4) *TWONPI/AMT24
        BROMMT = PIM*FINT(5) *TWONPI/AMT24
        BOMEMT = PIM*FINT(6) *TWONPI/AMT24
        BOMMMS = PIM*FINT(7) *TWONPI/AMT24
        BOMMMT = PIM*FINT(8) *TWONPI/AMT24
        BFIEMT = PIM*FINT(9) *TWONPI/AMT24
        BFIMMS = PIM*FINT(10)*TWONPI/AMT24
        BFIMMT = PIM*FINT(11)*TWONPI/AMT24
      ENDIF
      IF(IM2.GE.2) THEN
**      Vector-meson functions for 1+2f/g
        CALL SPECVR(RH,0D0,D2RO,DRO2K,AKRO)
**      Operators for remaining 1/M^2 corrections for vector mesons
        CALL OM2VEC(RH,D2PI(1),D2PI(2),DRO2K(1),DRO2K(2),
     .              ROIIM2S,ROIIM2T,ROIIM2L)
        CALL OM2VK(RH,D2PI(1),D2PI(2),D2PI(3),D2PI(4),
     .              I2RO,D2RO(1),D2RO(2), ROIIS,ROIIT)
        CALL OM2VK1(RH,D2PI(1),D2PI(2),D2PI(3),D2PI(4),I2RO,
     .              ROIIS1,ROIIT1)
        CALL OM2VK1(RH,D2PI(1),D2PI(2),D2PI(3),D2PI(4),I2OM,
     .              OMIIS1,OMIIT1)
        CALL OM2VK1(RH,D2PI(1),D2PI(2),D2PI(3),D2PI(4),I2FI,
     .              FIIIS1,FIIIT1)
        CALL TMEINT(RH,4,NP,FINT)
        BROM2S = PIM*FINT(1) *TWONPI/AMT24
        BROM2T = PIM*FINT(2) *TWONPI/AMT24
        BROM2L = PIM*FINT(3) *TWONPI/AMT24
        BROSK  = PIM*FINT(4) *TWONPI/AMT24
        BROTK  = PIM*FINT(5) *TWONPI/AMT24
        BROSK1 = PIM*FINT(6) *TWONPI/AMT24
        BROTK1 = PIM*FINT(7) *TWONPI/AMT24
        BOMSK1 = PIM*FINT(8) *TWONPI/AMT24
        BOMTK1 = PIM*FINT(9) *TWONPI/AMT24
        BFISK1 = PIM*FINT(10)*TWONPI/AMT24
        BFITK1 = PIM*FINT(11)*TWONPI/AMT24
      ENDIF
      IF(IM2.EQ.3) THEN
**      Operators for 1/M^2 corrections for scalar mesons
        CALL OM2SCA(RH,D2PI(1),D2PI(2),D2PI(3),D2A0(1),D2A0(2),
     .              A0IIM2S,A0IIM2T,A0IIM2L)
cp      CALL OM2SCA(RH,D2PI(1),D2PI(2),D2PI(3),DI2A2D,DDI2A2D,
cp   .              A2IIM2S,A2IIM2T,A2IIM2L)
        CALL TMEINT(RH,5,NP,FINT)
        BA0M2S = PIM*FINT(1)*TWONPI/AMT24
        BA0M2T = PIM*FINT(2)*TWONPI/AMT24
        BA0M2L = PIM*FINT(3)*TWONPI/AMT24
cp      BA2M2S = PIM*FINT(4)*TWONPI/AMT24
cp      BA2M2T = PIM*FINT(5)*TWONPI/AMT24
cp      BA2M2L = PIM*FINT(6)*TWONPI/AMT24
      ENDIF

C*    Two-pion exchange:
      FPI4=FPI*FPI*FPI*FPI
      VCPI = FPI4*0.5D0*(FPA-FCR)*BPIPIC
      VSPI =-FPI4*0.5D0*(FPA+FCR)*BPIPIS
      VTPI =-FPI4*0.5D0*(FPA+FCR)*BPIPIT
      VLPI = 0D0
      IF(IPV.EQ.1) THEN
        VCPI = VCPI + FCR*FPI4*(D2PI(1)*(PIM2*D2PI(1)-DI0PI))/AMT
        VLPI = VLPI + FCR*FPI4*(D2PI(2)+D2PI(1)/RH)*D2PI(1)*2D0/RH/AMT
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAPS(RH,D2PI(1),D2PI(2),D2PI(3),D4PI(1),D4PI(2),D4PI(3),
     .             CPIPI,SPIPI,TPIPI,OPIPI)
        VCPI = VCPI + FPI4/AMT * ( 0.5D0*FPA-FCR)*CPIPI
        VSPI = VSPI + FPI4/AMT * (-0.5D0*FPA-FCR)*SPIPI
        VTPI = VTPI + FPI4/AMT * (-0.5D0*FPA-FCR)*TPIPI
        IF(INA.GE.2) VLPI = VLPI + FPI4/AMT * 0.5D0*FPA*OPIPI
        IF(INA.GE.3) THEN
**        extra part form iterated OPE diagrams (i.e., TMO diagrams)
          VCPI = VCPI + IPV*FPA*FPI4
     .                  *(D2PI(1)*(PIM2*D2PI(1)-DI0PI))/2D0/AMT
          VLPI = VLPI + IPV*FPA*FPI4*(D2PI(2)-D2PI(1)/RH)*D2PI(1)/RH/AMT
        ENDIF
      ENDIF

      VCPM = 0D0
      VSPM = 0D0
      VTPM = 0D0
      VLPM = 0D0

C*    Pseudoscalar (eta):
      FF=FPI*FPI*FETA*FETA * FISO
      VSPM = VSPM - FF*2D0*BPIETS
      VTPM = VTPM - FF*2D0*BPIETT
      IF(IPV.EQ.1) THEN
        VCPM = VCPM + FF*( (PIM2+AMET2)*D2PI(1)*D2ET(1)
     .                     -DI0PI*D2ET(1)-D2PI(1)*DI0PS )/AMT
        VLPM = VLPM + FF*( (D2PI(2)+D2PI(1)/RH)*D2ET(1)
     .                    +(D2ET(2)+D2ET(1)/RH)*D2PI(1) )*2D0/RH/AMT
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAPS(RH,D2PI(1),D2PI(2),D2PI(3),D4ET(1),D4ET(2),D4ET(3),
     .             CPIET24,SPIET24,TPIET24,OPIET24)
        CALL ONAPS(RH,D4PI(1),D4PI(2),D4PI(3),D2ET(1),D2ET(2),D2ET(3),
     .             CPIET42,SPIET42,TPIET42,OPIET42)
        VCPM = VCPM - FF/AMT * (CPIET24+CPIET42)*0.5D0
        VSPM = VSPM - FF/AMT * (SPIET24+SPIET42)*1.5D0
        VTPM = VTPM - FF/AMT * (TPIET24+TPIET42)*1.5D0
        IF(INA.GE.2) VLPM = VLPM + FF/AMT * (OPIET24+OPIET42)*0.5D0
        IF(INA.GE.3) THEN
**        extra part form iterated diagrams (i.e., TMO diagrams)
          VCPM = VCPM + IPV*FF*( (PIM2+AMET2)*D2PI(1)*D2ET(1)
     .                          -DI0PI*D2ET(1)-D2PI(1)*DI0PS )/2D0/AMT
          VLPM = VLPM + IPV*FF*( (D2PI(2)-D2PI(1)/RH)*D2ET(1)
     .                          +(D2ET(2)-D2ET(1)/RH)*D2PI(1) )/RH/AMT
        ENDIF
      ENDIF
C*    Pseudoscalar (eta'):
      FF=FPI*FPI*FETP*FETP * FISO
      VSPM = VSPM - FF*2D0*BPIEPS
      VTPM = VTPM - FF*2D0*BPIEPT
      IF(IPV.EQ.1) THEN
        VCPM = VCPM + FF*( (PIM2+AMEP2)*D2PI(1)*D2EP(1)
     .                     -DI0PI*D2EP(1)-D2PI(1)*DI0PS )/AMT
        VLPM = VLPM + FF*( (D2PI(2)+D2PI(1)/RH)*D2EP(1)
     .                    +(D2EP(2)+D2EP(1)/RH)*D2PI(1) )*2D0/RH/AMT
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAPS(RH,D2PI(1),D2PI(2),D2PI(3),D4EP(1),D4EP(2),D4EP(3),
     .             CPIEP24,SPIEP24,TPIEP24,OPIEP24)
        CALL ONAPS(RH,D4PI(1),D4PI(2),D4PI(3),D2EP(1),D2EP(2),D2EP(3),
     .             CPIEP42,SPIEP42,TPIEP42,OPIEP42)
        VCPM = VCPM - FF/AMT * (CPIEP24+CPIEP42)*0.5D0
        VSPM = VSPM - FF/AMT * (SPIEP24+SPIEP42)*1.5D0
        VTPM = VTPM - FF/AMT * (TPIEP24+TPIEP42)*1.5D0
        IF(INA.GE.2) VLPM = VLPM + FF/AMT * (OPIEP24+OPIEP42)*0.5D0
        IF(INA.GE.3) THEN
**        extra part form iterated diagrams (i.e., TMO diagrams)
          VCPM = VCPM + IPV*FF*( (PIM2+AMEP2)*D2PI(1)*D2EP(1)
     .                          -DI0PI*D2EP(1)-D2PI(1)*DI0PS )/2D0/AMT
          VLPM = VLPM + IPV*FF*( (D2PI(2)-D2PI(1)/RH)*D2EP(1)
     .                          +(D2EP(2)-D2EP(1)/RH)*D2PI(1) )/RH/AMT
        ENDIF
      ENDIF

C*    Scalar (a0):
      FF=FPI*FPI*GA0*GA0*PIOMS2
      VSPM = VSPM + FF*(FPA-FCR)*BPIA0S
      VTPM = VTPM + FF*(FPA-FCR)*BPIA0T
      IF(IPV.EQ.1) THEN
        VSPM = VSPM - FF*FCR*D2PI(1)*D2A0(1) / (3D0*AMT)
        VTPM = VTPM - FF*FCR*D2PI(1)*D2A0(1) / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),D4A0(1),SPIA024,TPIA024)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),D2A0(1),SPIA042,TPIA042)
        VSPM = VSPM + FF/AMT * (0.5D0*FPA-FCR)*(SPIA024+SPIA042)
        VTPM = VTPM + FF/AMT * (0.5D0*FPA-FCR)*(TPIA024+TPIA042)
        IF(INA.GE.3) THEN
          VSPM = VSPM + FF*FPA*(PIS*I2A0-IPV*D2PI(1)*D2A0(1))/(6D0*AMT)
          VTPM = VTPM + FF*FPA*(PIT*I2A0-IPV*D2PI(1)*D2A0(1))/(6D0*AMT)
        ENDIF
      ENDIF
      IF(IM2.GE.3) THEN
        VSPM = VSPM + FF*(FPA-FCR)*BA0M2S
        VTPM = VTPM + FF*(FPA-FCR)*BA0M2T
        VLPM = VLPM + FF*(FPA-FCR)*BA0M2L
      ENDIF
C*    Scalar (epsilon):
      FF=FPI*FPI*GEP*GEP*PIOMS2 * FISO
      IF(IPV.EQ.1) THEN
        VSPM = VSPM - FF*D2PI(1)*D2SI(1) / (3D0*AMT)
        VTPM = VTPM - FF*D2PI(1)*D2SI(1) / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),D4SI(1),SPISI24,TPISI24)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),D2SI(1),SPISI42,TPISI42)
        VSPM = VSPM - FF/AMT * (SPISI24+SPISI42)/2D0
        VTPM = VTPM - FF/AMT * (TPISI24+TPISI42)/2D0
        IF(INA.GE.3) THEN
          VSPM = VSPM + FF*(PIS*I2SI-IPV*D2PI(1)*D2SI(1))/(6D0*AMT)
          VTPM = VTPM + FF*(PIT*I2SI-IPV*D2PI(1)*D2SI(1))/(6D0*AMT)
        ENDIF
      ENDIF
C*    Scalar (f0):
      FF=FPI*FPI*GF0*GF0*PIOMS2 * FISO
      IF(IPV.EQ.1) THEN
        VSPM = VSPM - FF*D2PI(1)*D2F0(1) / (3D0*AMT)
        VTPM = VTPM - FF*D2PI(1)*D2F0(1) / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),D4F0(1),SPIF024,TPIF024)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),D2F0(1),SPIF042,TPIF042)
        VSPM = VSPM - FF/AMT * (SPIF024+SPIF042)/2D0
        VTPM = VTPM - FF/AMT * (TPIF024+TPIF042)/2D0
        IF(INA.GE.3) THEN
          VSPM = VSPM + FF*(PIS*I2F0-IPV*D2PI(1)*D2F0(1))/(6D0*AMT)
          VTPM = VTPM + FF*(PIT*I2F0-IPV*D2PI(1)*D2F0(1))/(6D0*AMT)
        ENDIF
      ENDIF

C*    Vector (rho):
      FF=FPI*FPI*GRO*GRO*PIOMS2
      VSPM = VSPM + FF * (FPA-FCR)*BROEES
      VTPM = VTPM + FF * (FPA-FCR)*BROEET
      IF(IPV.EQ.1) THEN
        VSPM = VSPM + FF*FCR*D2PI(1)*D2RO(1) / (3D0*AMT)
        VTPM = VTPM + FF*FCR*D2PI(1)*D2RO(1) / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),D4RO(1),SPIRO24,TPIRO24)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),D2RO(1),SPIRO42,TPIRO42)
        VSPM = VSPM - FF/AMT * (0.5D0*FPA-FCR)*(SPIRO24+SPIRO42)
        VTPM = VTPM - FF/AMT * (0.5D0*FPA-FCR)*(TPIRO24+TPIRO42)
        IF(INA.GE.3) THEN
          VSPM = VSPM - FF*FPA*(PIS*I2RO-IPV*D2PI(1)*D2RO(1))/(6D0*AMT)
          VTPM = VTPM - FF*FPA*(PIT*I2RO-IPV*D2PI(1)*D2RO(1))/(6D0*AMT)
        ENDIF
      ENDIF
      IF(IM2.GE.1) THEN
        VCPM = VCPM + FF * (FPA-FCR)*BROMMC
        VSPM = VSPM + FF * (FPA+FCR)*BROMMS
        VTPM = VTPM + FF * (FPA+FCR)*(BROEMT+BROMMT)
        VLPM = VLPM + FF * (FPA-FCR)*BROEML
        IF(IM2.GE.2) THEN
          VSPM = VSPM + FF*((FPA-FCR)*(BROSK+BROM2S)-(FPA+FCR)*BROSK1)
          VTPM = VTPM + FF*((FPA-FCR)*(BROTK+BROM2T)-(FPA+FCR)*BROTK1)
          VLPM = VLPM + FF*(FPA-FCR)*BROM2L
        ENDIF
      ENDIF
C*    Vector (omega):
      FF=FPI*FPI*GOM*GOM*PIOMS2 * FISO
      IF(IPV.EQ.1) THEN
        VSPM = VSPM + FF*D2PI(1)*D2OM(1) / (3D0*AMT)
        VTPM = VTPM + FF*D2PI(1)*D2OM(1) / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),D4OM(1),SPIOM24,TPIOM24)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),D2OM(1),SPIOM42,TPIOM42)
        VSPM = VSPM + FF/AMT * (SPIOM24+SPIOM42)/2D0
        VTPM = VTPM + FF/AMT * (TPIOM24+TPIOM42)/2D0
        IF(INA.GE.3) THEN
          VSPM = VSPM - FF*(PIS*I2OM-IPV*D2PI(1)*D2OM(1))/(6D0*AMT)
          VTPM = VTPM - FF*(PIT*I2OM-IPV*D2PI(1)*D2OM(1))/(6D0*AMT)
        ENDIF
      ENDIF
      IF(IM2.GE.1) THEN
        VSPM = VSPM + FF * 2D0*BOMMMS
        VTPM = VTPM + FF * 2D0*(BOMEMT+BOMMMT)
        IF(IM2.GE.2) THEN
          VSPM = VSPM - FF * 2D0*BOMSK1
          VTPM = VTPM - FF * 2D0*BOMTK1
        ENDIF
      ENDIF
C*    Vector (phi):
      FF=FPI*FPI*GFI*GFI*PIOMS2 * FISO
      IF(IPV.EQ.1) THEN
        VSPM = VSPM + FF*D2PI(1)*D2FI(1) / (3D0*AMT)
        VTPM = VTPM + FF*D2PI(1)*D2FI(1) / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),D4FI(1),SPIFI24,TPIFI24)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),D2FI(1),SPIFI42,TPIFI42)
        VSPM = VSPM + FF/AMT * (SPIFI24+SPIFI42)/2D0
        VTPM = VTPM + FF/AMT * (TPIFI24+TPIFI42)/2D0
        IF(INA.GE.3) THEN
          VSPM = VSPM - FF*(PIS*I2FI-IPV*D2PI(1)*D2FI(1))/(6D0*AMT)
          VTPM = VTPM - FF*(PIT*I2FI-IPV*D2PI(1)*D2FI(1))/(6D0*AMT)
        ENDIF
      ENDIF
      IF(IM2.GE.1) THEN
        VSPM = VSPM + FF * 2D0*BFIMMS
        VTPM = VTPM + FF * 2D0*(BFIEMT+BFIMMT)
        IF(IM2.GE.2) THEN
          VSPM = VSPM - FF * 2D0*BFISK1
          VTPM = VTPM - FF * 2D0*BFITK1
        ENDIF
      ENDIF

C*    Pi-Pomeron (isosinglet):
      FF=FPI*FPI*FPOM2 * FISO
      IF(IPV.EQ.1) THEN
        VSPM = VSPM + FF*D2PI(1)*DI2POM / (3D0*AMT)
        VTPM = VTPM + FF*D2PI(1)*DI2POM / (3D0*AMT)
      ENDIF
      IF(INA.NE.0) THEN
        CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),DI4POM,SPIPO24,TPIPO24)
        CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),DI2POM,SPIPO42,TPIPO42)
        VSPM = VSPM + FF/AMT * (SPIPO24+SPIPO42)/2D0
        VTPM = VTPM + FF/AMT * (TPIPO24+TPIPO42)/2D0
        IF(INA.GE.3) THEN
          VSPM = VSPM - FF*(PIS*I2POM-IPV*D2PI(1)*DI2POM)/(6D0*AMT)
          VTPM = VTPM - FF*(PIT*I2POM-IPV*D2PI(1)*DI2POM)/(6D0*AMT)
        ENDIF
      ENDIF
C*    Pi-diffractive (a2):
cp    FF=FPI*FPI*FA2D2
cp    VSPM = VSPM - FF*(FPA-FCR)*BPIA2S
cp    VTPM = VTPM - FF*(FPA-FCR)*BPIA2T
cp    IF(IPV.EQ.1) THEN
cp      VSPM = VSPM + FF*FCR*D2PI(1)*DI2A2D / (3D0*AMT)
cp      VTPM = VTPM + FF*FCR*D2PI(1)*DI2A2D / (3D0*AMT)
cp    ENDIF
cp    IF(INA.NE.0) THEN
cp      CALL ONAVS(RH,D2PI(1),D2PI(2),D2PI(3),DI4A2D,SPIA224,TPIA224)
cp      CALL ONAVS(RH,D4PI(1),D4PI(2),D4PI(3),DI2A2D,SPIA242,TPIA242)
cp      VSPM = VSPM - FF/AMT * (0.5D0*FPA-FCR)*(SPIA224+SPIA242)
cp      VTPM = VTPM - FF/AMT * (0.5D0*FPA-FCR)*(TPIA224+TPIA242)
cp      IF(INA.GE.3) THEN
cp        VSPM = VSPM - FF*FPA*(PIS*I2A2D-IPV*D2PI(1)*DI2A2D)/(6D0*AMT)
cp        VTPM = VTPM - FF*FPA*(PIT*I2A2D-IPV*D2PI(1)*DI2A2D)/(6D0*AMT)
cp      ENDIF
cp    ENDIF
cp    IF(IM2.GE.3) THEN
cp      VSPM = VSPM - FF*(FPA-FCR)*BA2M2S
cp      VTPM = VTPM - FF*(FPA-FCR)*BA2M2T
cp      VLPM = VLPM - FF*(FPA-FCR)*BA2M2L
cp    ENDIF

C*    Scalar-scalar exchange (epsilon, pomeron, and diffractive):
***   Use: 2/pi int dl/l**2 (1-exp[-l**2/2ma**2])=sqrt(2/pi)/ma
cp    FF=0.5D0*FA2D2*FA2D2*(FPA-FCR)*(SQ2/SQPI)/AMA2D
cp    VCPM = VCPM + FF*I2A2D*I2A2D
cp    IF(IM2.GE.3) THEN
cp      VCPM = VCPM - FF*2D0*(DDI2A2D+2D0/RH*DI2A2D)*I2A2D / AMT24
cp      VLPM = VLPM - FF*4D0*I2A2D*DI2A2D/RH / AMT24
cp    ENDIF
cp    IF(INA.NE.0) THEN
cp      FFSS=(GEP*GEP*PIOMS2)**2
cp      FFSP=-GEP*GEP*PIOMS2*FPOM2
cp      FFSA=-GEP*GEP*PIOMS2*FA2D2*FISO
cp      FFPP= FPOM2*FPOM2
cp      FFPA= FPOM2*FA2D2*FISO
cp      FFAA= FA2D2*FA2D2*(FCR-0.5D0*FPA)
cp      VCPM = VCPM + ( FFSS*D2SI(1)*D4SI(1)
cp   .                 +FFSP*(D2SI(1)*DI4POM+D4SI(1)*DI2POM)
cp   .                 +FFSA*(D2SI(1)*DI4A2D+D4SI(1)*DI2A2D)
cp   .                 +FFPP*DI2POM*DI4POM
cp   .                 +FFPA*(DI2POM*DI4A2D+DI4POM*DI2A2D)
cp   .                 +FFAA*DI2A2D*DI4A2D*2D0 ) / (2D0*AMT)
cp    ENDIF

      RETURN
      END
************************************************************************
      SUBROUTINE TMEINT(RH,NC,NP,SUM)
**      Calculates integrals for pi-meson terms using Gauss integration
**      on (0,infinity) to (-1,1) for hyperbolic mapping with NP points
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION SUM(11),FUNCW(11),X(20),W(20),NDIM(5)
cp    DATA ICAL/0/, NDIM/3,10,11,11,6/
      DATA ICAL/0/, NDIM/3,8,11,11,3/
      SAVE X,W

      IF(ICAL.EQ.0) THEN
        CALL GAUSS(NP,X,W,1)
        DO 1 I=1,NP
          W(I)=2D0*W(I)/(1D0-X(I))**2
          X(I)=(1D0+X(I))/(1D0-X(I))
    1   CONTINUE
        ICAL=1
      ENDIF

      KDIM=NDIM(NC)
      DO 10 K=1,KDIM
   10   SUM(K)=0D0
      DO 11 IY=1,NP
        CALL FUNCT(FUNCW,X(IY),RH,NC)
        DO 12 K=1,KDIM
   12     SUM(K)=SUM(K)+FUNCW(K)*W(IY)
   11 CONTINUE

      RETURN
      END
************************************************************************
      SUBROUTINE FUNCT(FINT,Y,RH,NC)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 FINT(11), F2PI(2),F2PI3(3),F2PI4(4),
     .       F2ET(2),F2EP(2),F2A0(2), F2RO(2),F2OM(2),F2FI(2),
     .       F2RO1(2),F2RO2(2), I2A2D
      REAL*8 FROK(1),FRO2K(2),FROK2(2),FOMK(1),FFIK(1),FOMK2(2),FFIK2(2)
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/KAPPAV/AKRO,AKOM,AKFI
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4,
     .              GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4
      COMMON/IITRMS/PPIIC,PPIIS,PPIIT,
     .              ETIIS,ETIIT, EPIIS,EPIIT,
     .              A0IIS,A0IIT, A2IIS,A2IIT,I2A2D,DI2A2D,DDI2A2D,
     .              RIIEES,RIIEET,RIIEMT,RIIEML,RIIMMC,RIIMMS,RIIMMT,
     .              OIIEMT,OIIMMS,OIIMMT, FIIEMT,FIIMMS,FIIMMT,
     .              ROIIM2S,ROIIM2T,ROIIM2L, ROIIS,ROIIT,ROIIS1,ROIIT1,
     .              OMIIS1,OMIIT1, FIIIS1,FIIIT1,
     .              A0IIM2S,A0IIM2T,A0IIM2L, A2IIM2S,A2IIM2T,A2IIM2L

      PLAM=Y*PIM
      PLAM2=PLAM*PLAM
      IF(NC.EQ.1) THEN
***        FINT = dlam/lam**2 [(I2 operator I2) - (F operator F)](pi,pi)
        CALL DFI2DR(2,RH,PIM, ALAMP1,PLAM,F0PI, F2PI)
        CALL OPSEU(RH,F2PI(1),F2PI(2),F2PI(1),F2PI(2),PIFFC,PIFFS,PIFFT)
        FINT(1) = (PPIIC-PIFFC) / PLAM2
        FINT(2) = (PPIIS-PIFFS) / PLAM2
        FINT(3) = (PPIIT-PIFFT) / PLAM2
      ELSEIF(NC.EQ.2) THEN
***        FINT = dlam/lam**2 [(I2 operator I2) - (F operator F)]
***        Leading-order terms
        CALL DFI2DR(2,RH,PIM, ALAMP1,PLAM,F0PI, F2PI)
        CALL DFI2DR(2,RH,AMETA,ALAMP0,PLAM,F0ET,F2ET)
        CALL DFI2DR(2,RH,AMETP,ALAMP0,PLAM,F0EP,F2EP)
        CALL I2NYM(RH,AMR1,ALAMV1,PLAM,F0RO1)
        CALL I2NYM(RH,AMR2,ALAMV1,PLAM,F0RO2)
        F0RO=ARO*F0RO1+BRO*F0RO2
        CALL I2NYM(RH,A0M,ALAMS1,PLAM,F0A0)
        CALL OPSEU(RH,F2PI(1),F2PI(2),F2ET(1),F2ET(2),ETFFC,ETFFS,ETFFT)
        CALL OPSEU(RH,F2PI(1),F2PI(2),F2EP(1),F2EP(2),EPFFC,EPFFS,EPFFT)
        CALL OSCL(RH,F2PI(1),F2PI(2),F0A0,A0FFS,A0FFT)
        CALL OVC1(RH,F2PI(1),F2PI(2),F0RO,RFFEES,RFFEET)
cp      F2A2D=FDEXP(-PLAM2/(4D0*AMA2D2))*I2A2D
cp      CALL OSCL(RH,F2PI(1),F2PI(2),F2A2D,A2FFS,A2FFT)
        FINT(1) = (ETIIS-ETFFS) / PLAM2
        FINT(2) = (ETIIT-ETFFT) / PLAM2
        FINT(3) = (EPIIS-EPFFS) / PLAM2
        FINT(4) = (EPIIT-EPFFT) / PLAM2
        FINT(5) = (A0IIS-A0FFS) / PLAM2
        FINT(6) = (A0IIT-A0FFT) / PLAM2
        FINT(7) = (RIIEES-RFFEES) / PLAM2
        FINT(8) = (RIIEET-RFFEET) / PLAM2
cp      FINT(9) = (A2IIS-A2FFS) / PLAM2
cp      FINT(10)= (A2IIT-A2FFT) / PLAM2
      ELSEIF(NC.EQ.3) THEN
***        FINT = dlam/lam**2 [(I2 operator I2) - (F operator F)]
***        1/M**2 for (1+\kappa) and (1+\kappa)**2 in vector
        CALL DFI2DR(2,RH,PIM, ALAMP1,PLAM,F0PI, F2PI)
        CALL DFI2DR(2,RH,AMR1,ALAMV1,PLAM,F0RO1,F2RO1)
        CALL DFI2DR(2,RH,AMR2,ALAMV1,PLAM,F0RO2,F2RO2)
        F2RO(1)=ARO*F2RO1(1)+BRO*F2RO2(1)
        F2RO(2)=ARO*F2RO1(2)+BRO*F2RO2(2)
        CALL DFI2DR(2,RH,AMOM,ALAMV0,PLAM,F0OM,F2OM)
        CALL DFI2DR(2,RH,AMFI,ALAMV0,PLAM,F0FI,F2FI)
**      Vector-meson functions for 1+f/g and (1+f/g)**2
        CALL SPECV(1,RH,PLAM,F2RO,FROK,FROK2,AKRO)
        CALL SPECV(2,RH,PLAM,F2OM,FOMK,FOMK2,AKOM)
        CALL SPECV(3,RH,PLAM,F2FI,FFIK,FFIK2,AKFI)
**      Operators for 1/M^2 corrections for (1+f/g) and (1+f/g)**2
        CALL OVC2(RH,F2PI(1),F2PI(2),FROK(1),FROK2(1),FROK2(2),
     .            RFFEMT,RFFMMS,RFFMMT)
        CALL OVC3(RH,F2PI(1),F2PI(2),FROK(1),FROK2(1),FROK2(2),
     .            RFFEML,RFFMMC)
        CALL OVC2(RH,F2PI(1),F2PI(2),FOMK(1),FOMK2(1),FOMK2(2),
     .            OFFEMT,OFFMMS,OFFMMT)
        CALL OVC2(RH,F2PI(1),F2PI(2),FFIK(1),FFIK2(1),FFIK2(2),
     .            FFFEMT,FFFMMS,FFFMMT)
        FINT(1) = (RIIEMT-RFFEMT) / PLAM2
        FINT(2) = (RIIEML-RFFEML) / PLAM2
        FINT(3) = (RIIMMC-RFFMMC) / PLAM2
        FINT(4) = (RIIMMS-RFFMMS) / PLAM2
        FINT(5) = (RIIMMT-RFFMMT) / PLAM2
        FINT(6) = (OIIEMT-OFFEMT) / PLAM2
        FINT(7) = (OIIMMS-OFFMMS) / PLAM2
        FINT(8) = (OIIMMT-OFFMMT) / PLAM2
        FINT(9) = (FIIEMT-FFFEMT) / PLAM2
        FINT(10)= (FIIMMS-FFFMMS) / PLAM2
        FINT(11)= (FIIMMT-FFFMMT) / PLAM2
      ELSEIF(NC.EQ.4) THEN
***        FINT = dlam/lam**2 [(I2 operator I2) - (F operator F)]
***        Remaining 1/M**2 terms in vector
        CALL DFI2DR(4,RH,PIM, ALAMP1,PLAM,F0PI,F2PI4)
        CALL DFI2DR(2,RH,AMR1,ALAMV1,PLAM,F0RO1,F2RO1)
        CALL DFI2DR(2,RH,AMR2,ALAMV1,PLAM,F0RO2,F2RO2)
        F0RO=ARO*F0RO1+BRO*F0RO2
        F2RO(1)=ARO*F2RO1(1)+BRO*F2RO2(1)
        F2RO(2)=ARO*F2RO1(2)+BRO*F2RO2(2)
        CALL I2NYM(RH,AMOM,ALAMV0,PLAM,F0OM)
        CALL I2NYM(RH,AMFI,ALAMV0,PLAM,F0FI)
**      Vector-meson function for 1+2f/g
        CALL SPECVR(RH,PLAM,F2RO,FRO2K,AKRO)
        CALL OM2VEC(RH,F2PI4(1),F2PI4(2),FRO2K(1),FRO2K(2),
     .              ROFFM2S,ROFFM2T,ROFFM2L)
        CALL OM2VK(RH,F2PI4(1),F2PI4(2),F2PI4(3),F2PI4(4),
     .             F0RO,F2RO(1),F2RO(2), ROFFS,ROFFT)
        CALL OM2VK1(RH,F2PI4(1),F2PI4(2),F2PI4(3),F2PI4(4),F0RO,
     .              ROFFS1,ROFFT1)
        CALL OM2VK1(RH,F2PI4(1),F2PI4(2),F2PI4(3),F2PI4(4),F0OM,
     .              OMFFS1,OMFFT1)
        CALL OM2VK1(RH,F2PI4(1),F2PI4(2),F2PI4(3),F2PI4(4),F0FI,
     .              FIFFS1,FIFFT1)
        FINT(1) = (ROIIM2S-ROFFM2S) / PLAM2
        FINT(2) = (ROIIM2T-ROFFM2T) / PLAM2
        FINT(3) = (ROIIM2L-ROFFM2L) / PLAM2
        FINT(4) = (ROIIS -ROFFS ) / PLAM2
        FINT(5) = (ROIIT -ROFFT ) / PLAM2
        FINT(6) = (ROIIS1-ROFFS1) / PLAM2
        FINT(7) = (ROIIT1-ROFFT1) / PLAM2
        FINT(8) = (OMIIS1-OMFFS1) / PLAM2
        FINT(9) = (OMIIT1-OMFFT1) / PLAM2
        FINT(10)= (FIIIS1-FIFFS1) / PLAM2
        FINT(11)= (FIIIT1-FIFFT1) / PLAM2
      ELSEIF(NC.EQ.5) THEN
***        FINT = dlam/lam**2 [(I2 operator I2) - (F operator F)]
***        1/M**2 terms in scalar and diffractive
        CALL DFI2DR(3,RH,PIM,ALAMP1,PLAM,F0PI,F2PI3)
        CALL DFI2DR(2,RH,A0M,ALAMS1,PLAM,F0A0,F2A0)
        CALL OM2SCA(RH,F2PI3(1),F2PI3(2),F2PI3(3),F2A0(1),F2A0(2),
     .              A0FFM2S,A0FFM2T,A0FFM2L)
cp      CALL OM2SCA(RH,F2PI3(1),F2PI3(2),F2PI3(3),DI2A2D,DDI2A2D,
cp   .              A2FFM2S,A2FFM2T,A2FFM2L)
cp      FACA2D=FDEXP(-PLAM2/(4D0*AMA2D2))
        FINT(1) = (A0IIM2S-A0FFM2S) / PLAM2
        FINT(2) = (A0IIM2T-A0FFM2T) / PLAM2
        FINT(3) = (A0IIM2L-A0FFM2L) / PLAM2
cp      FINT(4) = (A2IIM2S-A2FFM2S*FACA2D) / PLAM2
cp      FINT(5) = (A2IIM2T-A2FFM2T*FACA2D) / PLAM2
cp      FINT(6) = (A2IIM2L-A2FFM2L*FACA2D) / PLAM2
      ENDIF

      RETURN
      END
************************************************************************
      SUBROUTINE OPSEU(R,F1,F2,G1,G2,CEN,SIG,TEN)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  (k1.k2)^2 F(r)G(r) = CEN
C*   ->  (sig1.k1xk2)(sig2.k1xk2) F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*
      CEN = F2*G2+2D0*F1*G1/(R*R)
      SIG = 2D0*(F1*G1/R+F1*G2+F2*G1)/(R*3D0)
      TEN = ((F1/R-F2)*G1+F1*(G1/R-G2))/(R*3D0)
      RETURN
      END
************************************************************************
      SUBROUTINE OSCL(R,F1,F2,G0,SIG,TEN)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  (sig1.k1)(sig2.k1) F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*
      SIG =-(F2+2D0*F1/R)*G0/3D0
      TEN =-(F2-F1/R)*G0/3D0
      RETURN
      END
************************************************************************
      SUBROUTINE OVC1(R,F1,F2,G0,SIG,TEN)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  -(sig1.k1)(sig2.k1) F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*
      SIG = (F2+2D0*F1/R)*G0/3D0
      TEN = (F2-F1/R)*G0/3D0
      RETURN
      END
************************************************************************
      SUBROUTINE OVC2(R,F1,F2,G1,H1,H2,TEM,SMM,TMM)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  [-k1^2(sig1.k1 sig2.k2 + sig1.k2 sig2.k1)
C*        +2 k1.k2(sig1.k1)(sig2.k1) ] F(r)G(r) =
C*                          = TEM S_{12}
C*   ->  [ (sig1.k1)(sig2.k1)k2^2 - k1^2 k2^2(sig1.sig2)
C*        +(sig1.k1xk2)(sig2.k1xk2) ] F(r)G(r) =
C*                          = SMM (sig1.sig2) + TMM S_{12}
C*
      TEM =-2D0*(F2-F1/R)*G1/R
      SMM =-((F2+F1/R)*(H2+H1/R)+2D0*F1*H1/(R*R))*2D0/3D0
      TMM = (2D0*(F2-F1/R)*H2-F2*(H2-H1/R))/3D0
      RETURN
      END
************************************************************************
      SUBROUTINE OVC3(R,F1,F2,G1,H1,H2,SOEM,CMM)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  -2i Q.k1xk2(sig1+sig2).k1 F(r)G(r) = SOEM L.S
C*   ->  [ k1^2 k2^2 - (k1.k2)^2 ] F(r)G(r) = CMM
C*
      SOEM= 4D0*F1*G1/(R*R)
      CMM = 2D0*(F1*H2+F2*H1+F1*H1/R)/R
      RETURN
      END
************************************************************************
      SUBROUTINE ONAPS(R,F1,F2,F3,G1,G2,G3,CEN,SIG,TEN,SO)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  (k1.k2)^3 F(r)G(r) = CEN
C*   ->  (k1.k2)(sig1.k1xk2)(sig2.k1xk2) F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*   ->  i(k1.k2)(sig1+sig2).(k1xk2) Q.(k1-k2) F(r)G(r) = SO L.S
C*
      R2=R*R
      CEN =-6D0*(F2-F1/R)*(G2-G1/R)/R2-F3*G3
      SIG = 2D0*((F1/R-F2+R*F3)*(G1/R-G2+R*G3)/R2-F3*G3)/3D0
      TEN =-((2D0*(F1/R-F2)+0.5D0*R*F3)*(2D0*(G1/R-G2)+0.5D0*R*G3)/R2
     .       -0.25D0*F3*G3)/3D0
      SO  =-4D0*(F2-F1/R)*(G2-G1/R)/R2
      RETURN
      END
************************************************************************
      SUBROUTINE ONAVS(R,F1,F2,F3,G1,SIG,TEN)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  (k1.k2)(sig1.k1)(sig2.k1) F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*
      SIG = (F3+2D0*(F2-F1/R)/R)*G1/3D0
      TEN = (F3-(F2-F1/R)/R)*G1/3D0
      RETURN
      END
************************************************************************
      SUBROUTINE OM2VEC(R,F1,F2,G1,G2,SIG,TEN,SO)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  [ -1/2 k1^2(sig1.k1 sig2.k2 + sig1.k2 sig2.k1)
C*         + (k2^2+k1.k2)(sig1.k1)(sig2.k1) ] F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*   ->  - i Q.k1xk2(sig1+sig2).k1 F(r)G(r) = SO L.S
C*
      SIG = (F2+2D0*F1/R)*(G2+2D0*G1/R) /3D0
      TEN = (F2-F1/R)*(G2-G1/R) /3D0
      SO  = 2D0*F1*G1/(R*R)
      RETURN
      END
************************************************************************
      SUBROUTINE OM2VK(R,F1,F2,F3,F4,G0,G1,G2,SIG,TEN)
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA IPV/1/,IPIMSQ/1/
C*    IPIMSQ=1 includes 1/M**2 terms in pseudovector vertices
C*
C*   ->  [( k1^2+3k1.k2+k2^2 )(sig1.k1)(sig2.k1)
C*        - 1/2 k1.k2(sig1.k1 sig2.k2 + sig1.k2 sig2.k1)
C*        + k1.k2(sig1.k1)(sig2.k1) ] F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*
      SIG = ( (F4+4D0*F3/R)*G0 + (F2+2D0*F1/R)*(G2+2D0*G1/R)
     .        +(3+IPIMSQ)*(F3+2D0*(F2-F1/R)/R)*G1
     .       - IPIMSQ*IPV*(F2*G2+2D0*F1*G1/(R*R)) )/3D0
      TEN = ( (F4+(F3-6D0*(F2-F1/R)/R)/R)*G0 + (F2-F1/R)*(G2+2D0*G1/R)
     .        +(3+IPIMSQ)*(F3-(F2-F1/R)/R)*G1
     .       - IPIMSQ*IPV*(F2*G2-F1*G1/(R*R)) )/3D0
      RETURN
      END
************************************************************************
      SUBROUTINE OM2VK1(R,F1,F2,F3,F4,G0,SIG,TEN)
      IMPLICIT REAL*8 (A-H,O-Z)
C*
C*   ->  k1^2(sig1.k1)(sig2.k1) F(r)G(r) = SIG (sig1.sig2) + TEN S_{12}
C*
      SIG = (F4+4D0*F3/R)*G0/3D0
      TEN = (F4+(F3-6D0*(F2-F1/R)/R)/R)*G0/3D0
      RETURN
      END
************************************************************************
      SUBROUTINE OM2SCA(R,F1,F2,F3,G1,G2,SIG,TEN,SO)
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA IPV/1/,IPIMSQ/1/
C*    IPIMSQ=1 includes 1/M**2 terms in pseudovector vertices
C*
C*   ->  [ (k2^2+2k1.k2)(sig1.k1)(sig2.k1)
C*        - 1/2 k1^2(sig1.k1 sig2.k2 + sig1.k2 sig2.k1)
C*        + 1/2 k1.k2(sig1.k1 sig2.k2 + sig1.k2 sig2.k1)
C*        - k1.k2(sig1.k1)(sig2.k1) ] F(r)G(r) =
C*                          = SIG (sig1.sig2) + TEN S_{12}
C*   ->  - i Q.k1xk2(sig1+sig2).k1 F(r)G(r) = SO L.S
C*
      SIG = ( (F2+2D0*F1/R)*(G2+2D0*G1/R)
     .       + (1-IPIMSQ)*(F3+2D0*(F2-F1/R)/R)*G1
     .       + IPIMSQ*IPV*(F2*G2+2D0*F1*G1/(R*R)) ) /3D0
      TEN = ( (F2-F1/R)*(G2-G1/R) + (1-IPIMSQ)*(F3-(F2-F1/R)/R)*G1
     .       + IPIMSQ*IPV*(F2*G2-F1*G1/(R*R)) ) /3D0
      SO  = 2D0*F1*G1/(R*R)
      RETURN
      END
************************************************************************
      SUBROUTINE SPECVR(RH,PLAM,D2RO,DRO2K,KAPPA)
**    (1+2f/g) function for rho exchange
**    Relevant when vector(g) and tensor(f) form factors are different
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 KAPPA,D2RO(2),DRO2K(2),FRO1(2),FRO2(2)
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/CUTTEN/ALAMT1,ALAVT1,ALAMT0,ALAVT0

      CALL DFI2DR(2,RH,AMR1,ALAVT1,PLAM,F0RO1,FRO1)
      CALL DFI2DR(2,RH,AMR2,ALAVT1,PLAM,F0RO2,FRO2)

      DRO2K(1)=D2RO(1)+2D0*KAPPA*(ARO*FRO1(1)+BRO*FRO2(1))
      DRO2K(2)=D2RO(2)+2D0*KAPPA*(ARO*FRO1(2)+BRO*FRO2(2))

      RETURN
      END
************************************************************************
      SUBROUTINE SPECV(N,RH,PLAM,D2VV,D2VT,D2TT,KAPPA)
**    (1+f/g) and (1+f/g)**2 functions for rho/omega/phi exchange
**    Relevant when vector(g) and tensor(f) form factors are different
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 KAPPA,D2VV(2),D2VT(1),D2TT(2), DUM1(2),DUM2(2),
     .       FRO1(2),FRO2(2)
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2,
     .              AVSC1,AVSC2, AEPS,AME1,BEPS,AME2
      COMMON/CUTTEN/ALAMT1,ALAVT1,ALAMT0,ALAVT0

      IF(N.EQ.1) THEN
**      rho exchange
        CALL DFI2DR(2,RH,AMR1,ALAVT1,PLAM,F0RO1,FRO1)
        CALL DFI2DR(2,RH,AMR2,ALAVT1,PLAM,F0RO2,FRO2)
        DUM1(1)=ARO*FRO1(1)+BRO*FRO2(1)
        DUM1(2)=ARO*FRO1(2)+BRO*FRO2(2)
        CALL DFI2DR(2,RH,AMR1,ALAMT1,PLAM,F0RO1,FRO1)
        CALL DFI2DR(2,RH,AMR2,ALAMT1,PLAM,F0RO2,FRO2)
        DUM2(1)=ARO*FRO1(1)+BRO*FRO2(1)
        DUM2(2)=ARO*FRO1(2)+BRO*FRO2(2)
      ELSEIF(N.EQ.2) THEN
**      omega exchange
        CALL DFI2DR(2,RH,AMOM,ALAVT0,PLAM,F0OM,DUM1)
        CALL DFI2DR(2,RH,AMOM,ALAMT0,PLAM,F0OM,DUM2)
      ELSEIF(N.EQ.3) THEN
**      phi exchange
        CALL DFI2DR(2,RH,AMFI,ALAVT0,PLAM,F0FI,DUM1)
        CALL DFI2DR(2,RH,AMFI,ALAMT0,PLAM,F0FI,DUM2)
      ELSE
**      STOP '*** SPECV in PIMES: incorrect call'
      ENDIF

      D2VT(1)=D2VV(1)+KAPPA*DUM1(1)
      D2TT(1)=D2VV(1)+KAPPA*(2D0*DUM1(1)+KAPPA*DUM2(1))
      D2TT(2)=D2VV(2)+KAPPA*(2D0*DUM1(2)+KAPPA*DUM2(2))

      RETURN
      END
***********************************************************************
      SUBROUTINE GAUSS(N,X,W,NT)
C*    This subroutine calculates the weight-factors and argument
C*    distributions for a Gaussian integration in (-1 , +1).
C*    N must be 4,8,12,16, or 20.
C*    otherwise the weight-factors and distributions will be those
C*    of a two-point Simpson integration.
C*    NT=1 : N-points Gauss
C*    NT=2 : N-points Gauss on (-1,0) and N-points Gauss on (0,+1)
C-----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION WGAUS1(30), WGAUS2(30), WGAUSS(60), X(50), W(50)
      EQUIVALENCE (WGAUS1(1),WGAUSS(1)), (WGAUS2(1),WGAUSS(31))

      DATA WGAUS1
     .    /.339981043584856D0, .652145154862546D0,
     .     .861136311594053D0, .347854845137454D0,
     .     .183434642495650D0, .362683783378362D0,
     .     .525532409916329D0, .313706645877887D0,
     .     .796666477413627D0, .222381034453374D0,
     .     .960289856497536D0, .101228536290376D0,
     .     .125233408511469D0, .249147045813403D0,
     .     .367831498998180D0, .233492536538355D0,
     .     .587317954286617D0, .203167426723066D0,
     .     .769902674194305D0, .160078328543346D0,
     .     .904117256370475D0, .106939325995318D0,
     .     .981560634246719D0, .047175336386512D0,
     .     .0950125098376374D0, .1894506104550685D0,
     .     .2816035507792589D0, .1826034150449236D0,
     .     .4580167776572274D0, .1691565193950025D0/
      DATA WGAUS2  /
     .     .6178762444026437D0, .1495959888165767D0,
     .     .7554044083550030D0, .1246289712555339D0,
     .     .8656312023878317D0, .0951585116824928D0,
     .     .9445750230732326D0, .0622535239386479D0,
     .     .9894009349916499D0, .0271524594117541D0,
     .     .0765265211334973D0, .1527533871307259D0,
     .     .2277858511416451D0, .1491729864726037D0,
     .     .3737060887154196D0, .1420961093183821D0,
     .     .5108670019508271D0, .1316886384491766D0,
     .     .6360536807265150D0, .1181945319615184D0,
     .     .7463319064601508D0, .1019301198172404D0,
     .     .8391169718222188D0, .0832767415767047D0,
     .     .9122344282513259D0, .0626720483341091D0,
     .     .9639719272779138D0, .0406014298003869D0,
     .     .9931285991850949D0, .0176140071391521D0/

      IF(MOD(N,4).NE.0) GOTO 2

      K=N/4-1
      INIT=2*K*(K+1)+1
      M=N/2
      DO 1 I=1,M
        J=INIT-2+I*2
        X(M+I)=WGAUSS(J)
        X(M+1-I)=-WGAUSS(J)
        W(M+I)=WGAUSS(J+1)
        W(M+1-I)=WGAUSS(J+1)
    1 CONTINUE
      IF(NT.EQ.1) RETURN
      DO 10 I=1,N
        X(N+I)=(X(I)+1D0)/2D0
        X(I)=(X(I)-1D0)/2D0
        W(N+I)=W(I)/2D0
        W(I)=W(N+I)
   10 CONTINUE
      RETURN

***   Simpson's rule: only odd N
    2 H=2D0/(N-1D0)
      X(1)=-1D0
      W(1)=H/3D0
      W23=(2D0*H)/3D0
      W43=2D0*W23
      NN=(N-1)/2
      DO 3 I=1,NN
        J=2*I
        X(J)=X(J-1)+H
        W(J)=W43
        X(J+1)=X(J)+H
        W(J+1)=W23
    3 CONTINUE
      X(N)=1D0
      W(N)=W(1)
      RETURN
      END
***********************************************************************
      SUBROUTINE AXIAL(R,ISO,SPIN,VSIG,VTEN,VLS,FI0AX,FI1AX,FI2AX)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FI(7,3)
      INTEGER SPIN
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/MEANM/ PIM,ROM,A0M
      COMMON/SCALIN/AMT,AMPV,AMT24
      COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2,
     .              AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2,
     .              AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2
      COMMON/HBARC/ HBC
      DATA SR3/1.7320508075688772935D0/, PI/3.14159265358979323846D0/
      DATA AX1/0D0/, THAX/0D0/, ALAX/0.4D0/
      DATA AMA1/1260D0/, AMF1/1285D0/, AME/1420D0/

      RH=R/HBC
c     FA1= FPI*(AMA1/AMPV) * DSQRT(AMA1**2-ROM**2)/ROM
      FA1= FPI*(AMA1/AMPV)
      AX8= FA1*(4D0*ALAX-1D0)/SR3
      COST=DCOS(THAX*PI/180D0)
      SINT=DSIN(THAX*PI/180D0)
      FF1 =COST*AX8-SINT*AX1
      FE  =SINT*AX8+COST*AX1

      FA12=(4D0*ISO-3D0)*FA1*FA1
      FF12=FF1*FF1
      FE2 =FE *FE
      FAC1=FA12/(AMA1*AMA1)
      FAC2=FF12/(AMF1*AMF1)
      FAC3=FE2 /(AME *AME )

      CALL REGGEFF(RH,AMA1,ALAMP1,FI(1,1))
      CALL REGGEFF(RH,AMF1,ALAMP0,FI(1,2))
      CALL REGGEFF(RH,AME ,ALAMP0,FI(1,3))

      VSIG=-FA12*(FI(1,1)+AMYNI/6D0*FI(2,1))
     .     -FF12*(FI(1,2)+AMYNI/6D0*FI(2,2))
     .     -FE2 *(FI(1,3)+AMYNI/6D0*FI(2,3))
     .     +(FAC1*FI(2,1)+FAC2*FI(2,2)+FAC3*FI(2,3))/3D0
      VTEN=+(FA12*FI(6,1)+FF12*FI(6,2)+FE2 *FI(6,3))*AMYNI/4D0
     .     +(FAC1*FI(6,1)+FAC2*FI(6,2)+FAC3*FI(6,3))
      VLS =-(FA12*FI(4,1)+FF12*FI(4,2)+FE2 *FI(4,3))*AMYNI/2D0

      FNONLOC=REDM*(4D0*SPIN-3D0)*5D0/6D0*AMYNI
      FI0AX=-FNONLOC*(FA12*FI(1,1)+FF12*FI(1,2)+FE2*FI(1,3))
      FI1AX= FNONLOC*(FA12*FI(4,1)+FF12*FI(4,2)+FE2*FI(4,3))*RH
      FI2AX=-FNONLOC*(FA12*FI(2,1)+FF12*FI(2,2)+FE2*FI(2,3))
     .      -FNONLOC*(FA12*FI(4,1)+FF12*FI(4,2)+FE2*FI(4,3))*2D0

      RETURN
      END
***********************************************************************
      SUBROUTINE PIOME(R,TYPE,EFAC,VPIS,VPIT)
C*    This subroutine computes the pion potential from the new
C*    Nijmegen model with the M/E factor.
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 PHI(7),PHIC(7)
      CHARACTER TYPE*2
      COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0,
     .              AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC
      COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2,
     .              ALVE,THV ,GV1, GRO,GOM,GFI,  GRO2,GOM2,GFI2,
     .              ALVM,     FV1, FRO,FOM,FFI,  FRO2,FOM2,FFI2,
     .              ALGS,THGS,GS1, GA0,GEP,GF0,  GA02,GEP2,GF02,
     .              GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI,
     .              FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2
      COMMON/CUTOFF/ALAM1,ALAMP1,ALAMV1,ALAMS1,
     .              ALAM0,ALAMP0,ALAMV0,ALAMS0
      COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2,
     .              AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2,
     .              AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2
      COMMON/HBARC/ HBC

      RH=R/HBC
      CALL REGGEFF(RH,AMPI,ALAMP1,PHI)
      FF=FPI2*EFAC
      VPIS = FF * PHI(2)/3D0
      VPIT = FF * PHI(6)

      IF(TYPE.EQ.'NP' .OR. TYPE.EQ.'PN') THEN
        CALL REGGEFF(RH,AMPIC,ALAMP1,PHIC)
        FF=FPIC2/ALF * EFAC
        VPIS = VPIS + FF*PHIC(2)/3D0
        VPIT = VPIT + FF*PHIC(6)
      ENDIF

      RETURN
      END
