C============================================================================
C                nCTEQ Parton Distribution Functions:
C     
C     SetNPDF calls the PDS file by name, and uses ReadPDS
C	from CTEQ 6 distribution
C
C     FIO: 3 June 2011
C===========================================================================
      Subroutine SetNPDF (filename,ierr)
      Implicit Double Precision (A-H,O-Z)
      Character filename*(*)
      Character Tablefile*132
      Logical fmtpds  !*** New/Old Format of PDS file
      Common /Setchange/ Isetch
      save
      Isetch=1  !*** Tell Ctq6Pdf there is a new set:
      ierr=0
C ------------------------------------------------------------------------
      CALL TRMSTR(filename,Ifilename)  !*** find length of directory name
      Tablefile=filename(1:Ifilename)
      IU= NextUn()
      Open(IU, File=Tablefile, Status='OLD', Err=100)
      fmtpds=.TRUE.   !*** USE NEW PDS FORMAT   
      Call Readpds (IU,fmtpds)
      Close (IU)
      Return

 100  Print *, ' Data file ', Tablefile, ' cannot be opened '
     >  //'in SetNPDF!!'
c      Stop
      ierr=1
      return
      End
C============================================================================
C                nCTEQ Parton Distribution Functions:
C     
C     Set2NPDF allows to call two PDS files by their names, 
C     and uses ReadPDS from CTEQ 6 distribution.
C
C     KK: 14 June 2011
C===========================================================================
      Subroutine Set2NPDF (filename1,filename2,ierr)
      Implicit Double Precision (A-H,O-Z)
      Character filename1*(*)
      Character filename2*(*)
      Character Tablefile*132
      Logical fmtpds  !*** New/Old Format of PDS file
      Common /Setchange/ Isetch
	Integer i
C ====
      PARAMETER (MXX = 201, MXQ = 25, MXF = 6, MaxVal=4)
      PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
      Common
     > / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX)
     > / CtqPar2 / Nx, Nt, NfMx, MxVal
     > / XQrange / Qini, Qmax, Xmin
     > / QCDtable /  Alambda, Nfl, Iorder
     > / Masstbl / Amass(6)
C ====
      Common
     > / CtqPar1f / F_Al(2), F_XV(0:MXX,2), F_TV(0:MXQ,2), F_UPD(MXPQX,2)
     > / CtqPar2f / F_Nx(2), F_Nt(2), F_NfMx(2), F_MxVal(2)
     > / XQrangef / F_Qini(2), F_Qmax(2), F_Xmin(2)
     > / QCDtablef /  F_Alambda(2), F_Nfl(2), F_Iorder(2)
     > / Masstblf / F_Amass(6,2)

      save
      Isetch=1  !*** Tell Ctq6Pdf there is a new set:
      ierr=0
C ------------------------------------------------------------------------
      CALL TRMSTR(filename1,Ifilename)  !*** find length of directory name
      Tablefile=filename1(1:Ifilename)
      IU= NextUn()
      Open(IU, File=Tablefile, Status='OLD', Err=100)
      fmtpds=.TRUE.   !*** USE NEW PDS FORMAT   
      Call Readpds (IU,fmtpds)
      Close (IU)
	
	F_Al(1) = Al

	do i=0,MXX
	F_XV(i,1)= XV(i)
	end do

	do i=0,MXQ
	F_TV(i,1)= TV(i)
	end do

	do i=0,MXPQX
	F_UPD(i,1)= UPD(i)
	end do

	F_Nx(1) = Nx
	F_Nt(1) = Nt
	F_NfMx(1) = NfMx
	F_MxVal(1) = MxVal
			
	F_Qini(1) = Qini
	F_Qmax(1) = Qmax
	F_Xmin(1) = Xmin

	F_Alambda(1) = Alambda
	F_Nfl(1) = Nfl
	F_Iorder(1)  = Iorder

C ------------------------------------------------------------------------
      CALL TRMSTR(filename2,Ifilename)  !*** find length of directory name
      Tablefile=filename2(1:Ifilename)
      IU= NextUn()
      Open(IU, File=Tablefile, Status='OLD', Err=100)
      fmtpds=.TRUE.   !*** USE NEW PDS FORMAT   
      Call Readpds (IU,fmtpds)
      Close (IU)

	F_Al(2) = Al
	
	do i=0,MXX
	F_XV(i,2)= XV(i)
	end do
	
	do i=0,MXQ
	F_TV(i,2)= TV(i)
	end do
	
	do i=0,MXPQX
	F_UPD(i,2)= UPD(i)
	end do
	
	F_Nx(2) = Nx
	F_Nt(2) = Nt
	F_NfMx(2) = NfMx
	F_MxVal(2) = MxVal
			
	F_Qini(2) = Qini
	F_Qmax(2) = Qmax
	F_Xmin(2) = Xmin
	
	F_Alambda(2) = Alambda
	F_Nfl(2) = Nfl
	F_Iorder(2)  = Iorder

      Return

 100  Print *, ' Data file ', Tablefile, ' cannot be opened '
     >  //'in SetNPDF!!'
c      Stop
      ierr=1
      return
      End

C============================================================================
C                nCTEQ Parton Distribution Functions:
C     
C     CtqNpdf make the (A,Z) isospin correction and calls Ctq6Pdf
C     If (A,Z)=1, the "nuclear proton" is returned
C
C	CtqNpdf(iparton,x,Q,iA,iZ)
C
C	Arguments : iparton - ID of a parton as in CTEQ 6 PDFs
C			x,Q     - Bjorken x & energy scale
C			iA, iZ  - atomic number & number of protons used for 
C				    the isospin correction (Note: to use PDFs 
C				    of partons of bound nucleons in a nucleus 
C      			    with A,Z, one has to read in the file 
C				    fitname_A_Z.pds with 
C				    SetNPDF('fitname_A_Z.pds',ierr) and set 
C				    here the correct values of A,Z)
C
C     FIO: 3 June 2011
C===========================================================================
      function  CtqNpdf(iparton,x,Q,iA,iZ)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      A=Float(iA)
      Z=Float(iZ)

      if(iparton.eq.-2) then
       result=Ctq6Pdf(-2,x,Q)*Z/A + Ctq6Pdf(-1,x,Q)*(A-Z)/A
      elseif(iparton.eq.-1) then
       result=Ctq6Pdf(-1,x,Q)*Z/A + Ctq6Pdf(-2,x,Q)*(A-Z)/A
      elseif(iparton.eq.1)  then
       result=Ctq6Pdf( 1,x,Q)*Z/A + Ctq6Pdf( 2,x,Q)*(A-Z)/A
      elseif(iparton.eq.2)  then
       result=Ctq6Pdf( 2,x,Q)*Z/A + Ctq6Pdf( 1,x,Q)*(A-Z)/A
      else 
        result=Ctq6Pdf (iparton,x,Q)
      endif

      CtqNpdf=result
      return
      end
C============================================================================
C                nCTEQ Parton Distribution Functions:
C     
C     Ctq2Npdf makes the same isospin correction as CtqNpdf but allows
C	to switch between 2 pds files initialized by Set2NPDF.
C
C	Ctq2Npdf(iset,iparton,x,Q,iA,iZ)
C
C	Arguments : iset    - value 1 or 2 - ID of the pds file from 
C				    Set2NPDF(1,2,ierr)
C	            iparton - ID of a parton as in CTEQ 6 PDFs
C			x,Q     - Bjorken x & energy scale
C			iA, iZ  - atomic number & number of protons used for 
C				    the isospin correction (Note: to use PDFs 
C				    of partons of bound nucleons in a nucleus 
C      			    with A,Z, one has to read in the file 
C				    fitname_A_Z.pds with 
C				    SetNPDF('fitname_A_Z.pds',ierr) and set 
C				    here the correct values of A,Z)
C
C     KK: 14 June 2011
C===========================================================================

      function  Ctq2Npdf(iset,iparton,x,Q,iA,iZ)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
	Integer i
C ====
      PARAMETER (MXX = 201, MXQ = 25, MXF = 6, MaxVal=4)
      PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
      Common
     > / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX)
     > / CtqPar2 / Nx, Nt, NfMx, MxVal
     > / XQrange / Qini, Qmax, Xmin
     > / QCDtable /  Alambda, Nfl, Iorder
     > / Masstbl / Amass(6)
C ====
      Common
     > / CtqPar1f / F_Al(2), F_XV(0:MXX,2), F_TV(0:MXQ,2), F_UPD(MXPQX,2)
     > / CtqPar2f / F_Nx(2), F_Nt(2), F_NfMx(2), F_MxVal(2)
     > / XQrangef / F_Qini(2), F_Qmax(2), F_Xmin(2)
     > / QCDtablef /  F_Alambda(2), F_Nfl(2), F_Iorder(2)
     > / Masstblf / F_Amass(6,2)
C ====
      Common /Setchange/ Isetch  !*** set non-zero when set is changed
      data isetOLD /-1/
      save
c -------------------------------------------------------------------
      if((Isetch.ne.0).or.(iset.ne.isetOLD)) then 

	Al = F_Al(iset)
	
	do i=0,MXX
	XV(i) = F_XV(i,iset)
	end do
	
	do i=0,MXQ
	TV(i) = F_TV(i,iset)
	end do
	
	do i=0,MXPQX
	UPD(i) = F_UPD(i,iset)
	end do
	
	Nx = F_Nx(iset)
	Nt = F_Nt(iset)
	NfMx = F_NfMx(iset)
	MxVal = F_MxVal(iset)
			
	Qini = F_Qini(iset)
	Qmax = F_Qmax(iset)
	Xmin = F_Xmin(iset)
	
	Alambda = F_Alambda(iset)
	Nfl = F_Nfl(iset)
	Iorder = F_Iorder(iset)

        endif !*** END OF INITIALIZATION 
c -------------------------------------------------------------------
	
      A=Float(iA)
      Z=Float(iZ)

      if(iparton.eq.-2) then
       result=Ctq6Pdf(-2,x,Q)*Z/A + Ctq6Pdf(-1,x,Q)*(A-Z)/A
      elseif(iparton.eq.-1) then
       result=Ctq6Pdf(-1,x,Q)*Z/A + Ctq6Pdf(-2,x,Q)*(A-Z)/A
      elseif(iparton.eq.1)  then
       result=Ctq6Pdf( 1,x,Q)*Z/A + Ctq6Pdf( 2,x,Q)*(A-Z)/A
      elseif(iparton.eq.2)  then
       result=Ctq6Pdf( 2,x,Q)*Z/A + Ctq6Pdf( 1,x,Q)*(A-Z)/A
      else 
        result=Ctq6Pdf (iparton,x,Q)
      endif

      Ctq2Npdf=result

      isetOLD=iset
      return
      end
c =======================================================