! list of potnr  implemented here 
!      potnr=3     =>   Nijm93
!      potnr=7     =>   Nijm I
!      potnr=8     =>   Nijm II
!      potnr=11    =>   AV18  (needs input file!)
!      potnr=12    =>   separable (needs input file!) (for testing) 
!      potnr=31    =>   CDB 2000
!      potnr=63    =>   Idaho N3LO
!      potnr=67    =>   Idaho N3LO 600

#define POTNR 31
      
      MODULE potNN
       USE constants
       
#if(POTNR==3)      
       USE nijm
       INTEGER,PARAMETER :: idpars=0
#endif

#if(POTNR==7)      
       USE nijm
       INTEGER,PARAMETER :: idpars=1
#endif

#if(POTNR==8)      
       USE nijm
       INTEGER,PARAMETER :: idpars=2
#endif

#if(POTNR==11)      
       USE argonne
#endif
       
#if(POTNR==12)      
       USE yamaguchi
#endif
       
#if(POTNR==31)      
       USE cdbonnpot
#endif

#if(POTNR==63)      
       USE idahon3lopot
#endif

#if(POTNR==67)      
       USE idahon3lo600pot
#endif
        PRIVATE 

        PUBLIC potNNlsj
                 

       CONTAINS 
       
!    calculates matrix element in pw projection 
!    of <p' (l's)j | V | p (ls)j >    in fm^2 
!    normalization for momentum eigen states is 
!       <p' (l's')j' | p (ls)j >  =delta_jj' delta_ll' delta_ss' delta(p-p')/(pp')
!    input parameters 
!       w1=p' in fm-1 
!       w2=p in fm-1
!       lfix is the lowest l in the coupled channel 
!       sfix = s
!       jfix = j 
!       iforce = 1,2,3   ! pp,np,nn 
!     output: 
!      w6,w7,w8,w9   in fm^2 
!      for an uncoupled channel w6 = V (the single matrix element) 
!      for a coupled channel 
!       w6=V00=<l'=j-1 | V | l=j-1 >   
!       w7=V20=<l'=j+1 | V | l=j-1 >
!       w8=V02=<l'=j-1 | V | l=j+1 >
!       w9=V22=<l'=j+1 | V | l=j+1 >  
              

      SUBROUTINE potNNlsj(w1,w2,lfix,sfix,jfix,V00,V02,V20,V22,iforcein)
       IMPLICIT NONE
       REAL(8) w1,w2,w6,w7,w8,w9,mfaktor,fact,V00,V02,V20,V22
       INTEGER lfix,sfix,jfix,iforce,tfix,iforcein,iforcemin,iforcemax
       LOGICAL koup
       CHARACTER(LEN=2) force
       REAL(8) VEM(2,2)
       LOGICAL,SAVE :: first=.true.
       
       IF(first) THEN
        WRITE(*,*) 'POTNR =',POTNR
        first=.false.        
       END IF

       iforce=iforcein
       V00=0.0_dpreal
       V02=0.0_dpreal
       V20=0.0_dpreal
       V22=0.0_dpreal
              
       IF(jfix.NE.0 .AND. jfix.NE.lfix) THEN
          koup=.true.
       ELSE
          koup=.false.
       END IF

       SELECT CASE(iforce)
        CASE(1) 
           mfaktor=mprot
           force='PP'
        CASE(2) 
           mfaktor=mnuc
           force='NP'
        CASE(3)
           mfaktor=mneu
           force='NN'
        CASE DEFAULT
           WRITE(*,*) 'iforce in potNNlsj'
           CALL abort
       END SELECT
       
       w6=0.0
       w7=0.0
       w8=0.0
       w9=0.0

C     check whether isospin and iforce are compatible       
       tfix=(1+(-1)**(lfix+sfix))/2
       IF(tfix.EQ.0.AND.((iforce.EQ.1).OR.(iforce.EQ.3))) THEN
        WRITE(*,*) 'check ISOSPIN BOOKKEEPING!!'
        RETURN                  ! potential is zero in this case 
       END IF
       
#if(POTNR==3 || POTNR==7 || POTNR==8 )       
       CALL potnijmlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $                 force,idpars)
#endif

#if(POTNR==9 )       
       CALL potreidlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $                 force)
#endif

#if(POTNR==10)          
       mfaktor=hbarc/41.47       
       call potlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup)
#endif

#if(POTNR==11)          
       call potav18lsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif

#if(POTNR==12)  
       mfaktor=hbarc/41.47
       call potyglsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9)
#endif

#if(POTNR==13)          
       mfaktor=hbarc/41.47       
       call potmtlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9)
#endif
#if(POTNR>24 && POTNR<31)
!     Aufruf fuer Chiral Potential NLO
         mfaktor=hbarc/41.47105
         CALL potchirallsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $         0,POTNR-24,0)
#endif
#if(POTNR==31)          
       call potcdblsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif      
#if(POTNR>33 && POTNR<40)
!     Aufruf fuer Chiral Potential NNLO
         mfaktor=hbarc/41.47105
         CALL potchirallsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $         1,3,0)        !        POTNR-33) nur ein cut-off
#endif
#if(POTNR>39 && POTNR<46)
!     Aufruf fuer Chiral Potential NNLO-DELTA
         mfaktor=hbarc/41.47105
         CALL potchirallsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $         2,3,0)         !       POTNR-39) nur ein cut-off
#endif
#if(POTNR>45 && POTNR<52)
!     Aufruf fuer Chiral Potential NNLO*1
         mfaktor=hbarc/41.47105
         CALL potchirallsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $         3,POTNR-45,0)
#endif
#if(POTNR>51 && POTNR<58)
!     Aufruf fuer Chiral Potential NNLO*2
         mfaktor=hbarc/41.47105
         CALL potchirallsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,
     $         4,POTNR-51,0)
#endif       
       

#if(POTNR==60)     
       mfaktor=hbarc/41.47      
       call potpestlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup)
#endif

#if(POTNR==61)     
       mfaktor=hbarc/41.47105       
       call potnnlolsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup)
#endif

#if(POTNR==62)          
       call potvlowklsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,
     $                  koup,iforce)
#endif

#if(POTNR==63)          
       call potn3lolsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif

#if(POTNR==66)           
       CALL potchn3lolsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,
     $                        koup,ostat,cutnum,iforce)
#endif

#if(POTNR==67)          
       call potn3lo600lsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif

#if(POTNR==68)     
       mfaktor=hbarc/41.47105       
       call potnnlodeltalsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup)
#endif

#if(POTNR==69)          
       call potav18lsjnoempot(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,
     $                        koup,iforce)
#endif

#if(POTNR==70)      
       CALL potnijmesclsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce) 
#endif

#if(POTNR==71)      
       mfaktor=197.33/41.467       
       CALL potc3dlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup) 
#endif

#if(POTNR==72)      
       CALL potvsrglsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,iforce) 
#endif

#if(POTNR==73)          
       call potn3lowlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif

#if(POTNR==74)          
       call potdelshlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,iforce)
#endif

#if(POTNR==75)   
       CALL potchn2lolocallsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,
     $                       koup,ostat,cutnum,cutnumsfr,iforce)
#endif  
#if(POTNR==76)          
       call potsogNNlsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif 
       
#if(POTNR==77)   
       CALL potichloclsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,
     $                       koup,ostat,cutnum,iforce)
#endif  
             
#if(POTNR==78)          
       call potsogNNTPElsj(w1,w2,lfix,sfix,jfix,w6,w7,w8,w9,koup,iforce)
#endif              
       
C     EMPOT switches to full EM in nn & pp system 
C     no EMPOT in np system
C     if EMPOT undefined and COULOMB defined 
C     then point Coulomb in pp system
C     otherwise no EM interaction 
C     make sure that AV18 does not get any EM interaction 

#if(defined(EMPOT) && (POTNR!=11))  
       VEM(1:2,1:2)=0.0                
       IF(iforce.eq.1.or.iforce.EQ.3) THEN
        IF(jfix.eq.0) THEN
         CALL potempotlsj(w1,w2,lfix,sfix,jfix,
     X      VEM(1,1),VEM(1,2),VEM(2,1),VEM(2,2),
     X      koup,iforce)
        ELSE IF(jfix.ne.lfix) THEN
         IF((-1)**(lfix+sfix+1).le.0.0) THEN
          CALL potempotlsj(w1,w2,lfix,sfix,jfix,
     X       VEM(1,1),VEM(1,2),VEM(2,1),VEM(2,2),
     X       koup,iforce)
         ELSE
          VEM(1:2,1:2)=0.0
         END IF
        ELSE
         VEM(1:2,1:2)=0.0
         CALL potempotlsj(w1,w2,lfix,sfix,jfix,
     X      VEM(1,1),VEM(1,2),VEM(2,1),VEM(2,2),
     X      .false.,iforce)
         VEM(2,1)=0.0
         VEM(1,2)=0.0
         VEM(2,2)=0.0
        END IF                  ! coupled/j=0 or uncoupled 
        w6=w6+VEM(1,1)
        w7=w7+VEM(1,2)
        w8=w8+VEM(2,1)
        w9=w9+VEM(2,2)
       END IF                   ! pp or nn        
#else
#if(defined(COULOMB) && (POTNR!=11))
C     Coulombpotential aufrufen, falls pp Kanal and if not AV18
       IF(iforce.EQ.1) THEN
        VEM(1:2,1:2)=0.0
        IF(jfix.eq.0 .or. jfix.ne.lfix) THEN
         CALL potcoullsj(w1,w2,lfix,sfix,jfix,
     $      VEM(1,1),VEM(1,2),VEM(2,1),VEM(2,2),
     $      koup)
        ELSE
         CALL potcoullsj(w1,w2,lfix,sfix,jfix,
     $      VEM(1,1),VEM(1,2),VEM(2,1),VEM(2,2),
     $      .false.)
        END IF
        w6=w6+VEM(1,1)
        w7=w7+VEM(1,2)
        w8=w8+VEM(2,1)
        w9=w9+VEM(2,2)        
       END IF
#endif
#endif

       fact=1.0/(mfaktor*w1*w2)      ! output is fm-1 need fm^2

       
#ifdef NNAVEPOT
C  isospin to decide over how to average       
       tfix=(1+(-1)**(lfix+sfix))/2
       IF(tfix.eq.1) THEN
        SELECT CASE(iforce)
         CASE(1)
          V00=V00+w6*fact*PPAVEFAKT
          V02=V02+w7*fact*PPAVEFAKT
          V20=V20+w8*fact*PPAVEFAKT
          V22=V22+w9*fact*PPAVEFAKT
         CASE(2)
          V00=V00+w6*fact*NPAVEFAKT
          V02=V02+w7*fact*NPAVEFAKT
          V20=V20+w8*fact*NPAVEFAKT
          V22=V22+w9*fact*NPAVEFAKT
         CASE(3)
          V00=V00+w6*fact*NNAVEFAKT
          V02=V02+w7*fact*NNAVEFAKT
          V20=V20+w8*fact*NNAVEFAKT
          V22=V22+w9*fact*NNAVEFAKT
         CASE DEFAULT
          STOP 'problem with iforce in potNN averaging'
        END SELECT 
       ELSE
         V00=w6*fact
         V02=w7*fact
         V20=w8*fact
         V22=w9*fact
       END IF
       
       END DO  ! iforce
#else
       V00=w6*fact
       V02=w7*fact
       V20=w8*fact
       V22=w9*fact
#endif

      END SUBROUTINE

      END MODULE potNN

