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