c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c cc fitcry routine ccc c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine fitcry( aa , z , gamma , tem , rhom , qcrst ) implicit real*8 (a-h,o-z) real*8 mue,i,j,k,l,ip,jp,kp,lp dimension a (0:4),b (4),e(0:4),f(4),i (0:4),j (4),p(0:4),q(4) & ,ap(0:4),bp(4) ,ip(0:4),jp(4) 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 'cry-he.par' elseif ( idint(aa) .eq. 12 ) & then include 'cry-c.par' elseif ( idint(aa) .eq. 16 ) & then include 'cry-o.par' elseif ( idint(aa) .eq. 20 ) & then include 'cry-ne.par' elseif ( idint(aa) .eq. 24 ) & then include 'cry-mg.par' elseif ( idint(aa) .eq. 28 ) & then include 'cry-si.par' elseif ( idint(aa) .eq. 32 ) & then include 'cry-s.par' elseif ( idint(aa) .eq. 40 ) & then include 'cry-ca.par' elseif ( idint(aa) .eq. 56 ) & then include 'cry-fe.par' else call errcall( 3 ) return endif c u = 2.d0*pi*( log10(rho) - 3.d0 )/9.d0 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,171)------- fb = a( 0 )/2.d0 + c*u +d do 100 m = 1, 4 fm = float( m) fb = fb + a( m )*cos( fm*u ) + b( m )*sin( fm*u ) 100 continue c c -------F(u,5000)------- ft = e( 0 )/2.d0 + g*u +h do 110 m = 1, 4 fm = float( m) ft = ft + e( m )*cos( fm*u ) + f( m )*sin( fm*u ) 110 continue c c -------G(u,171)------- gb = i( 0 )/2.d0 + k*u +l do 120 m = 1, 4 fm = float( m) gb = gb + i( m )*cos( fm*u ) + j( m )*sin( fm*u ) 120 continue c c -------G(u,5000)------- gt = p( 0 )/2.d0 + r*u +s do 130 m = 1, 4 fm = float( m) gt = gt + p( m )*cos( fm*u ) + q( m )*sin( fm*u ) 130 continue c c vp= alphap0 + alphap1*gamma**(-1.d0/3.d0) & + alphap2*gamma**(-2.d0/3.d0) & + alphap3*gamma**(-1.d0) c wp= betap0 + betap1 *gamma**(-1.d0/3.d0) & + betap2 *gamma**(-2.d0/3.d0) & + betap3 *gamma**(-1.d0) c c -------F'(u,171)------- fp = ap( 0 )/2.d0 + cp*u +dp do 140 m = 1, 4 fm = float( m) fp = fp + ap( m )*cos( fm*u ) + bp( m )*sin( fm*u ) 140 continue c c -------G'(u,171)------- gp = ip( 0 )/2.d0 + kp*u +lp do 150 m = 1, 4 fm = float( m) gp = gp + ip( m )*cos( fm*u ) + jp( m )*sin( fm*u ) 150 continue c c -------- F_lat(u,Gamma) & G_lat(u,Gamma) ----------- flat = ( 1.d0 - v )*fb + v*ft glat = ( 1.d0 - w )*gb + w*gt c -------- F_phn(u,Gamma) & G_phn(u,Gamma) ----------- fphn = vp*fp gphn = wp*gp 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 vband = 3.56d-2*z*(rho6/aa)**(1.d0/3.d0)/t8 c -------------------------------------------------- if ( vband .ge. 150 ) & then fband = 0.d0 else fband = exp(-2.d0*vband) c write(*,*) 'vband & fband ',vband,fband endif c ------------------------------------------------- qlat = 0.5738d0*z**2.d0/aa*t8**6.d0*rho*fband & *(0.5d0*(( cv **2.d0 + ca **2.d0 ) & + dn*( cvp**2.d0 + cap**2.d0 ))*flat & - 0.5d0*(( cv **2.d0 - ca **2.d0 ) & + dn*( cvp**2.d0 - cap**2.d0 ))*glat) c c qphn = 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 ))*fphn & - 0.5d0*(( cv **2.d0 - ca **2.d0 ) & + dn*( cvp**2.d0 - cap**2.d0 ))*gphn) c qcrst = qlat + qphn return end