CCCCC/  Version 122611.  Fit to proton F2 and FL (or F1) data using target mass    /CCCCC 
CCCCC/  (TM) formalism described in I. Schienbein et al., J.Phys.G35:053101,2008.  /CCCCC 
CCCCC/  Subroutine takes input x, Q^2 and returns the full structure functions     /CCCCC
CCCCC/  which describe the data (F2,F1), as well as the TM corrected massless      /CCCCC
CCCCC/  limit structure functions (F2tmc, F1tmc).                                  /CCCCC
CCCCC/                                                                             /CCCCC
CCCCC/  passing dpar = parameter # returns results with parameter changed          /CCCCC
CCCCC/  by 1-sigma.  The errors vectors below are utilized for this                /CCCCC


      SUBROUTINE SFTM_SUB(x,q2,dpar,f2,f2tmc,f1,f1tmc,fl,fltmc)
      IMPLICIT none


      real*8 mp2,x,q2,t,f1,f1tmc,f2,f2tmc,fl,fltmc,r,rtmc
      real*8 parm1(5),parm2(5),parm2t(15),parm1t(15)
      real*8 eparm2t(15),eparm1t(15),parm2t2(15),parm1t2(15)
      integer i,j,k,dpar
      logical writeparm/.true./

      data parm2t /
     & -.21296E+01,0.25084E+01,0.10317E+00,0.65104E+01,0.21250E+00,
     & 0.88383E-01,0.38138E+01,-.40206E-01,-.49528E+00,0.55488E+01,
     & 0.27815E+00,-.49298E+01,0.28415E+00,-.51740E+00,0.99406E+00 /

      data parm1t /
     & -.20651E+01,0.10000E+01,0.00000E+00,0.00000E+00,0.78427E-01,
     & 0.00000E+00,0.00000E+00,-.29497E-01,0.00000E+00,0.00000E+00,
     & 0.60402E+00,-.48846E+01,0.34171E+00,-.34821E+00,0.00000E+00 /


CCC   Parameter error vectors   CCC


      data eparm2t /
     & 0.10982E-01,0.69624E-01,0.27532,0.51715,0.13371E-01,
     & 0.92819E-02,0.66797,0.26615E-02,0.50006E-01,0.71585,
     & 0.20970E-01,0.46174,0.31295E-01,0.75336E-01,0.18015 /

     
      data eparm1t /
     & 0.50737E-01,0.0,0.0,0.0,0.14325E-01,0.0,0.0,0.15097E-01,0.0,0.0,
     & 0.32949,2.7565,0.19065,0.10693,0.0 /

CCCCCCCCCCCCCCCCC

      parm2t2 = parm2t     !!!  Make copy for safe keeping
      parm1t2 = parm1t

      if(dpar.GT.0) then   !!!  change parameter dpar by 1-sigma
        parm2t(dpar) = parm2t(dpar) + eparm2t(dpar)
      elseif(dpar.LT.0) then
        parm1t(abs(dpar)) = parm1t(abs(dpar)) + eparm1t(abs(dpar))
      endif

c     write(24,*) parm2t(1)

c      write(6,*) dpar

      mp2 = 0.9382727*0.9382727

      t = 1.+4.*x*x*mp2/q2 


      parm1(1) =  parm1t(1)+parm1t(6)*exp(-1*q2/parm1t(7))
     &            +parm1t(8)*log(0.3**2+q2)
      parm1(2) =  parm1t(2)+parm1t(9)*exp(-1*q2/parm1t(10))
     &            +parm1t(11)*log(0.3**2+q2)
      parm1(3) =  parm1t(3)+parm1t(12)*exp(-1*q2/parm1t(13))
     &            +parm1t(14)*log(0.3**2+q2)
      parm1(4) =  parm1t(4)+parm1t(15)*log(0.3**2+q2)
      parm1(5) =  parm1t(5)


      parm2(1) =  parm2t(1)+parm2t(6)*exp(-1*q2/parm2t(7))
     &            +parm2t(8)*log(0.3**2+q2)
      parm2(2) =  parm2t(2)+parm2t(9)*exp(-1*q2/parm2t(10))
     &            +parm2t(11)*log(0.3**2+q2)
      parm2(3) =  parm2t(3)+parm2t(12)*exp(-1*q2/parm2t(13))
     &            +parm2t(14)*log(0.3**2+q2)
      parm2(4) =  parm2t(4)+parm2t(15)*log(0.3**2+q2)
      parm2(5) =  parm2t(5)

  
      parm2t = parm2t2     !!!  Restore original values for next call
      parm1t = parm1t2
  
      if(writeparm) then
        write(26,*) "$A_n^L$ ", " & ",parm1t(5)," & ","0.0"," & ","0.0",
     &            " & "," 0.0 "," \\"
        write(26,*) "$B_n^L$ ", " & ", parm1t(1)," & ","0.0"," & ",
     &             "0.0"," & ",parm1t(12)," \\"
        write(26,*) "$C_n^L$ ", " & ", parm1t(2)," & ",parm1t(6)," & ",
     &            parm1t(7), " & ", parm1t(10)," \\" 
        write(26,*) "$D_n^L$ ", " & ", parm1t(3)," & ",parm1t(8)," & ",
     &            parm1t(9), " & ",parm1t(11)," \\"
        write(26,*) "$E_n^L$ ", " & ",parm1t(4)," & ","0.0"," & ","0.0",
     &            " & "," 0.0 "," \\"
 
        write(26,*) "$A_n^2$ ", " & ",parm2t(5)," & ","0.0"," & ","0.0",
     &            " & "," 0.0 "," \\"
        write(26,*) "$B_n^2$ ", " & ", parm2t(1)," & ",parm2t(6)," & ",
     &            parm2t(7)," & ",parm2t(8)," \\"
        write(26,*) "$C_n^2$ ", " & ", parm2t(2)," & ",parm2t(9)," & ",
     &            parm2t(10), " & ", parm2t(11)," \\"
        write(26,*) "$D_n^2$ ", " & ", parm2t(3)," & ",parm2t(12)," & ",
     &            parm2t(13), " & ",parm2t(14)," \\"
        write(26,*) "$E_n^2$ ", " & ",parm2t(4)," & ","0.0"," & ","0.0",
     &            " & "," 0.0 "," \\"
      endif

      call sffit2(1,x,q2,1,parm1,parm2,fltmc)
      call sffit2(1,x,q2,0,parm1,parm2,fl)
      call sffit2(2,x,q2,1,parm2,parm2,f2tmc)
      call sffit2(2,x,q2,0,parm2,parm2,f2)

      f1tmc = (f2tmc-fltmc)/2./x
      f1 = (t*f2-fl)/2./x
      r = fl/2./x/f1
      rtmc = fltmc/2./x/f1tmc
      

 2001 format(6f8.5)

      end