c      ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      cc     Calculate f_photo                                    ccc
c      cc                   cfp.f                                  ccc
c      ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine cfp  ( fp, temp, lamb, xi )  
c
        implicit real*8 (a-h,o-z)
        real*8 lamb
        dimension temp( 1 ), lamb( 1 ), xi( 1 ), fp( 1 )
     &          , a(0:2), b(3), cc(0:2,0:6), d(0:2,1:5)
        common /param/ ndat, nn
        pi = 3.141592654d0
c
        do 100 n = 1 , ndat
          if     ( temp( n ) .lt. 1.d7 )
     &       then
               write(0,*) ' Temperature is less than 1.D7 (K)'
          elseif ( temp( n ) .lt. 1.d8 )
     &       then
               include 'c7-8.par'
               include 'd7-8.par'
               c   = 0.5654d0 + dlog10( temp( n )/1.d7 )
               tau = dlog10( temp( n )/1.d7 )
          elseif ( temp( n ) .lt. 1.d9 )
     &       then
               include 'c8-9.par'
               include 'd8-9.par'
               c   =  1.5654d0
               tau = dlog10( temp( n )/1.d8 )
             else
               include 'c9.par'
               include 'd9.par'
               c   =  1.5654d0
               tau = dlog10( temp( n )/1.d9 )
          endif
c
         b( 1 ) = 6.290d-3
         b( 2 ) = 7.483d-3
         b( 3 ) = 3.061d-4
c
         do  200 i = 0, 2
           a( i )  = 0.5d0*cc( i, 0 )
     &             + 0.5d0*cc( i, 6 )*dcos(1.d1*pi*tau)
          do 300 j = 1, 5
               pjt = pi*float( j )*tau
           a( i )  = a( i ) + cc( i, j )*dcos( 5.d0/3.d0 *pjt )
     &                      +  d( i, j )*dsin( 5.d0/3.d0 *pjt )
 300      continue
 200     continue
c
c        ---------------------------------------------------------
           if ( c*xi( n ) .gt. 1.d2 ) then 
              fp( n ) = 0.d0
c             write(0,*)'cfp.f: Warning',c*xi(n)
              goto 100
           endif
c        ---------------------------------------------------------
c
         fp( n ) = ( a( 0 ) + a( 1 )*xi( n ) + a( 2 )*xi( n )**2.d0)
     &           * dexp( -c*xi( n ))
     1           / ( xi( n )**3.d0 + b( 1 )/lamb( n )
     &                             + b( 2 )/lamb( n )**2.d0
     2                             + b( 3 )/lamb( n )**3.d0 ) 
 100    continue
        return
        end