c subroutine to calculate the trivial strong-coupling corretion
c for B phase (according to the table in J.W. Serene and D. Rainer, 
c Phys. Rep. 101, 221 (1983).
c uses linear extrapolation for large dcpcn
c input
c   gap = bare energy gap of the B phase (unit=k_BT_c)
c   dcpcn=specific heat jump in B phase =deltaC/Cnormal at T=T_c
c   temp = temperature/superfluid transition temperature
c output:
c   gap= corrected gap 
      subroutine wcp(temp,dcpcn,gap)
      implicit none
      double precision temp,dcpcn,gap
      double precision c(5),x(0:10,5)
      double precision wt1,wt2,delta,wc1,wc2,corr
      integer j,it,itp,ic,icp
      data(c(   j),j=1,5)/1.43,1.6,1.8,2.0,2.2/
      data(x(10,j),j=1,5)/1.,1.056,1.115,1.171,1.221/
      data(x( 9,j),j=1,5)/1.,1.048,1.097,1.141,1.180/
      data(x( 8,j),j=1,5)/1.,1.041,1.083,1.119,1.150/
      data(x( 7,j),j=1,5)/1.,1.036,1.072,1.102,1.128/
      data(x( 6,j),j=1,5)/1.,1.032,1.063,1.089,1.112/
      data(x( 5,j),j=1,5)/1.,1.028,1.056,1.079,1.099/
      data(x( 4,j),j=1,5)/1.,1.026,1.051,1.073,1.091/
      data(x( 3,j),j=1,5)/1.,1.024,1.049,1.069,1.086/
      data(x( 2,j),j=1,5)/1.,1.024,1.048,1.068,1.085/
      data(x( 1,j),j=1,5)/1.,1.024,1.048,1.068,1.085/
      data(x( 0,j),j=1,5)/1.,1.024,1.048,1.068,1.085/
      it=int(temp*10.-0.00001)
      if(it.ge.10)return
      itp=it+1
      wt1=(0.1*float(itp)-temp)/0.1
      wt2=1.-wt1
      do 100 ic=1,4
      icp=ic+1
100   if(dcpcn.lt.c(icp)) goto 110
      ic=4
      icp=5
110   delta=c(icp)-c(ic)
      wc1=(c(icp)-dcpcn)/delta
      wc2=1.-wc1
      corr=wt1*(wc1*x(it ,ic)+wc2*x(it ,icp))
     x    +wt2*(wc1*x(itp,ic)+wc2*x(itp,icp))
      gap=gap*corr
      return
      end