c      cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      cc          fitliq routine                                 ccc
c      cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine fitliq( aa , z , gamma , tem , rhom , qliq ) 
       implicit real*8 (a-h,o-z)
       real*8 mue,i,j,k,l
       dimension a(0:5),b(5),e(0:5),f(5),i(0:5),j(5),p(0:5),q(5)
       common /wa/ n
c
       mue  = aa/z
       rho  = rhom*mue
       rho6 = rho/1.d6
       t8   = tem/1.d8
       pi   = 3.1415927d0
c
c      -- read parameters ----
       if     ( idint(aa) .eq. 4  ) 
     &    then 
                                   include 'liq-he.par'
       elseif ( idint(aa) .eq. 12 ) 
     &    then 
                                   include 'liq-c.par'
       elseif ( idint(aa) .eq. 16 ) 
     &    then 
                                   include 'liq-o.par'
       elseif ( idint(aa) .eq. 20 ) 
     &    then 
                                   include 'liq-ne.par'
       elseif ( idint(aa) .eq. 24 ) 
     &    then
                                   include 'liq-mg.par'
       elseif ( idint(aa) .eq. 28 ) 
     &    then
                                   include 'liq-si.par'
       elseif ( idint(aa) .eq. 32 ) 
     &    then
                                   include 'liq-s.par'
       elseif ( idint(aa) .eq. 40 ) 
     &    then
                                   include 'liq-ca.par'
       elseif ( idint(aa) .eq. 56 ) 
     &    then
                                   include 'liq-fe.par'
       else
           call errcall( 3 )
           return
       endif
c      -----------------------------------------------------
c
       u     = 2.d0*pi*( log10(rho) - 3.d0 )/1.d1
c
           v =  alpha0 + alpha1 *gamma**(-1.d0/3.d0)
     &                 + alpha2 *gamma**(-2.d0/3.d0) 
     &                 + alpha3 *gamma**(-1.d0)
c
           w =  beta0  + beta1 *gamma**(-1.d0/3.d0)
     &                 + beta2 *gamma**(-2.d0/3.d0) 
     &                 + beta3 *gamma**(-1.d0)
c
c     -------F(u,1)-------
       fb    =  a( 0 )/2.d0  + c*u +d   
      do 100 m = 1, 5
         fm  = float( m)
         fb  =  fb + a( m )*cos( fm*u ) + b( m )*sin( fm*u )
 100  continue
c
c     -------F(u,160)-------
       ft    =  e( 0 )/2.d0  + g*u +h   
      do 110 m = 1, 5
         fm  = float( m)
         ft  =  ft + e( m )*cos( fm*u ) + f( m )*sin( fm*u )
 110  continue
c
c     -------G(u,1)-------
       gb    =  i( 0 )/2.d0  + k*u +l   
      do 120 m = 1, 5
         fm  = float( m)
         gb  =  gb + i( m )*cos( fm*u ) + j( m )*sin( fm*u )
 120  continue
c
c     -------G(u,160)-------
       gt    =  p( 0 )/2.d0  + r*u +s   
      do 130 m = 1, 5
         fm  = float( m)
         gt  =  gt + p( m )*cos( fm*u ) + q( m )*sin( fm*u )
 130  continue
c
c     -------- F_liq(u,Gamma) & G_liq(u,Gamma) -----------
      fliq   = v*fb + ( 1.d0 - v )*ft
      gliq   = w*gb + ( 1.d0 - w )*gt
c      
            dn = float( n )
        snthtw = 0.2325d0 
            cv = 0.5d0 + 2.d0*snthtw
            ca = 0.5d0
           cvp = 1.d0  - cv
           cap = 1.d0  - ca
c
       qliq    = 0.5738d0*z**2.d0/aa*t8**6.d0*rho
     &          *(0.5d0*(( cv **2.d0 + ca **2.d0 )
     &               + dn*( cvp**2.d0 + cap**2.d0 ))*fliq
     &          - 0.5d0*(( cv **2.d0 - ca **2.d0 )
     &               + dn*( cvp**2.d0 - cap**2.d0 ))*gliq)
c
       return
       end