c********************************************************************* SUBROUTINE collision(nb,zb,mij,xi,ro,drox,t,omega11,zij,zpij,zsij, 1 domega11,dzij,dzpij,dzsij) c routine private du module mod_evol c calcul des intégrales de collision et des coefficients de c résistance à l'aide des coefficients de Paquette c et al. ApJS 61, 177, 1986 c on a extrait des tables de Paquette et al. les valeurs c aux points de table, puis on effectue une interpolation c par sbsp1dn pour toutes les tables en distinguant potentiels c attractifs et répulsifs, on a étendu l'intervalle jusqu'à psi=4 c les tables ont été scannées par P. Somlyo (OCA) c Auteur: P.Morel, Département J.D. Cassini, O.C.A. c CESAM2k c entrées c nb: nombre d'éléments+électrons c zb : charges c mij : masses réduites, (i,j indices des particules) c xi : abondances c t : température c ro : densité c drox : 1/ro dro/dx_1 dérivée de la densité / H par mole à ro près c sorties c om11(i,j) : omega11_ij c zij(i,j),zpij(i,j),zsij(i,j) : coefficients de résistance c d...(i,j,k) : derivees / xk c----------------------------------------------------------------------- USE mod_donnees, ONLY : amu, echarg, kbol, pi USE mod_kind USE mod_numerique, ONLY : bsp1dn, no_croiss IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION (:,:) :: mij REAL (kind=dp), INTENT(in), DIMENSION (:) :: zb, xi REAL (kind=dp), INTENT(in) :: ro, t, drox INTEGER, INTENT(in) :: nb REAL (kind=dp), INTENT(out), DIMENSION (:,:,:) :: domega11, dzij, 1 dzpij, dzsij REAL (kind=dp), INTENT(out), DIMENSION (:,:) :: omega11, zij, 1 zpij, zsij REAL (kind=dp), DIMENSION (:,:), ALLOCATABLE :: dfxdx, gam0 REAL (kind=dp), DIMENSION (:), ALLOCATABLE :: dld, dgam INTEGER, PARAMETER :: m=4 REAL (kind=dp), SAVE, DIMENSION(4,60) :: pa, pr REAL (kind=dp), SAVE, DIMENSION(60+m) :: psit REAL (kind=dp), SAVE, DIMENSION(60) :: psi REAL (kind=dp), DIMENSION(4) :: f, dfdx REAL (kind=dp), SAVE :: e4kpi, ld0, a00, xi00 REAL (kind=dp) :: bid, bid1, bid2, bid3, ld1, a01, 1 ld, a0, lambda, gam, eps INTEGER, SAVE :: knot INTEGER :: i, j, k, l LOGICAL, SAVE :: init=.TRUE. c------------------------------------------------------------------------ c les psi entre -7 et 4 DATA(psi(i), i=1,60) 1/-7.00000D+00,-6.81356D+00,-6.62712D+00,-6.44068D+00,-6.25424D+00, 2 -6.06780D+00,-5.88136D+00,-5.69492D+00,-5.50847D+00,-5.32203D+00, 3 -5.13559D+00,-4.94915D+00,-4.76271D+00,-4.57627D+00,-4.38983D+00, 4 -4.20339D+00,-4.01695D+00,-3.83051D+00,-3.64407D+00,-3.45763D+00, 5 -3.27119D+00,-3.08475D+00,-2.89831D+00,-2.71186D+00,-2.52542D+00, 6 -2.33898D+00,-2.15254D+00,-1.96610D+00,-1.77966D+00,-1.59322D+00, 7 -1.40678D+00,-1.22034D+00,-1.03390D+00,-8.47458D-01,-6.61017D-01, 8 -4.74576D-01,-2.88136D-01,-1.01695D-01, 8.47458D-02, 2.71186D-01, 9 4.57627D-01, 6.44068D-01, 8.30508D-01, 1.01695D+00, 1.20339D+00, 1 1.38983D+00, 1.57627D+00, 1.76271D+00, 1.94915D+00, 2.13559D+00, 2 2.32203D+00, 2.50847D+00, 2.69492D+00, 2.88136D+00, 3.06780D+00, 3 3.25424D+00, 3.44068D+00, 3.62712D+00, 3.81356D+00, 4.00000D+00/ c potentiels repulsifs DATA((pr(i,j),i=1,4),j=1,10) 1/-6.05084D+00,-5.11871D+00,-3.86436D+00,-5.10215D+00,-5.90778D+00, 2 -8.10346D-01,-3.72631D+00,-4.95710D+00,-5.93827D+00,-1.09027D+00, 3 -3.58933D+00,-4.81272D+00,-5.62370D+00,-4.70052D+00,-3.45322D+00, 4 -4.66904D+00,-5.48273D+00,-4.56270D+00,-3.31797D+00,-4.52606D+00, 5 -5.34248D+00,-4.42570D+00,-3.18362D+00,-4.38381D+00,-5.20298D+00, 6 -4.28953D+00,-3.05019D+00,-4.24230D+00,-5.06424D+00,-4.15423D+00, 7 -2.91771D+00,-4.10155D+00,-4.92627D+00,-4.01980D+00,-2.78618D+00, 8 -3.96157D+00,-3.56874D+00,-3.88626D+00,-2.65563D+00,-3.82239D+00/ DATA((pr(i,j),i=1,4),j=11,20) 1/-4.65269D+00,-3.75364D+00,-2.52608D+00,-3.68402D+00,-4.51712D+00, 2 -3.62193D+00,-2.39754D+00,-3.54645D+00,-4.38235D+00,-3.49116D+00, 3 -2.27004D+00,-3.40972D+00,-4.24842D+00,-3.36134D+00,-2.14359D+00, 4 -3.27382D+00,-4.11532D+00,-3.23248D+00,-2.01819D+00,-3.13878D+00, 5 -3.98304D+00,-3.10456D+00,-1.89388D+00,-3.00458D+00,-3.85159D+00, 6 -2.97764D+00,-1.77063D+00,-2.87122D+00,-3.72095D+00,-2.85167D+00, 7 -1.64847D+00,-2.73872D+00,-3.59114D+00,-2.72665D+00,-1.52738D+00, 8 -2.60707D+00,-3.46211D+00,-2.60258D+00,-1.40737D+00,-2.47623D+00/ DATA((pr(i,j),i=1,4),j=21,30) 1/-3.33384D+00,-2.47943D+00,-1.28841D+00,-2.34619D+00,-3.20629D+00, 2 -2.35718D+00,-1.17048D+00,-2.21692D+00,-3.07942D+00,-2.23579D+00, 3 -1.05356D+00,-2.08836D+00,-2.95315D+00,-2.11519D+00,-9.37581D-01, 4 -1.96046D+00,-2.82741D+00,-1.99533D+00,-8.22494D-01,-1.83314D+00, 5 -2.70209D+00,-1.87610D+00,-7.08219D-01,-1.70630D+00,-2.57705D+00, 6 -1.75740D+00,-5.94656D-01,-1.57983D+00,-2.45215D+00,-1.63909D+00, 7 -4.81680D-01,-1.45356D+00,-2.32718D+00,-1.52100D+00,-3.69138D-01, 8 -1.32730D+00,-2.20190D+00,-1.40291D+00,-2.56838D-01,-1.20083D+00/ DATA((pr(i,j),i=1,4),j=31,40) 1/-2.07603D+00,-1.28456D+00,-1.44549D-01,-1.07386D+00,-1.94919D+00, 2 -1.16565D+00,-3.19871D-02,-9.46030D-01,-1.82095D+00,-1.04578D+00, 3 8.11908D-02,-8.16930D-01,-1.69076D+00,-9.24498D-01, 1.95401D-01, 4 -6.86035D-01,-1.55797D+00,-8.01225D-01, 3.11145D-01,-5.52715D-01, 5 -1.42179D+00,-6.75283D-01, 4.29027D-01,-4.16206D-01,-1.28125D+00, 6 -5.45845D-01, 5.49770D-01,-2.75580D-01,-1.13515D+00,-4.11925D-01, 7 6.74237D-01,-1.29721D-01,-9.82097D-01,-2.72344D-01, 8.03439D-01, 8 2.27054D-02,-8.20375D-01,-1.25723D-01, 9.38549D-01, 1.83273D-01/ DATA((pr(i,j),i=1,4),j=41,50) 1/-6.48000D-01, 2.95099D-02, 1.08088D+00, 3.53792D-01,-4.62708D-01, 2 1.95076D-01, 1.23184D+00, 5.36265D-01,-2.62055D-01, 3.72719D-01, 3 1.39279D+00, 7.32752D-01,-4.36401D-02, 5.63959D-01, 1.56486D+00, 4 9.45095D-01, 1.94499D-01, 7.69677D-01, 1.74440D+00, 1.17445D+00, 5 4.53275D-01, 9.89546D-01, 1.94330D+00, 1.42058D+00, 7.31826D-01, 6 1.22142D+00, 2.14712D+00, 1.68111D+00, 1.02683D+00, 1.46107D+00, 7 2.35646D+00, 1.95101D+00, 1.33203D+00, 1.70268D+00, 2.56682D+00, 8 2.22311D+00, 1.63774D+00, 1.94032D+00, 2.77417D+00, 2.48981D+00/ DATA((pr(i,j),i=1,4),j=51,60) 1/ 1.93183D+00, 2.16959D+00, 2.97619D+00, 2.74523D+00, 2.20471D+00, 2 2.38878D+00, 3.17294D+00, 2.98668D+00, 2.45554D+00, 2.59881D+00, 3 3.36483D+00, 3.21495D+00, 2.68704D+00, 2.80145D+00, 3.55480D+00, 4 3.43179D+00, 2.90923D+00, 3.00093D+00, 3.74490D+00, 3.64304D+00, 5 3.12475D+00, 3.19827D+00, 3.93393D+00, 3.84962D+00, 3.33470D+00, 6 3.39367D+00, 4.12251D+00, 4.05247D+00, 3.54024D+00, 3.58748D+00, 7 4.31073D+00, 4.25233D+00, 3.74226D+00, 3.78000D+00, 4.49864D+00, 8 4.44977D+00, 3.94145D+00, 3.97146D+00, 4.68629D+00, 4.64525D+00/ c potentiels attractifs DATA((pa(i,j),i=1,4),j=1,10) 1/-4.93441D+00,-4.20096D+00,-3.09773D+00,-4.62717D+00,-4.94282D+00, 2 -4.16990D+00,-3.00014D+00,-4.69405D+00,-4.76088D+00,-4.00403D+00, 3 -2.84717D+00,-4.67806D+00,-4.70791D+00,-4.00179D+00,-2.86022D+00, 4 -4.61155D+00,-4.52790D+00,-3.83158D+00,-2.69255D+00,-4.42109D+00, 5 -4.35635D+00,-3.61785D+00,-2.46151D+00,-4.30253D+00,-4.27179D+00, 6 -3.51722D+00,-2.32512D+00,-4.09100D+00,-4.20588D+00,-3.39848D+00, 7 -2.16726D+00,-3.76524D+00,-4.11357D+00,-3.24841D+00,-1.98526D+00, 8 -3.53614D+00,-3.99670D+00,-3.12369D+00,-1.83742D+00,-3.44443D+00/ DATA((pa(i,j),i=1,4),j=11,20) 1/-3.87962D+00,-2.99692D+00,-1.69452D+00,-3.43331D+00,-3.75390D+00, 2 -2.85481D+00,-1.54092D+00,-3.37121D+00,-3.66677D+00,-2.71368D+00, 3 -1.38226D+00,-3.19712D+00,-3.57333D+00,-2.56342D+00,-1.21986D+00, 4 -3.02235D+00,-3.36031D+00,-2.36359D+00,-1.04488D+00,-2.80215D+00, 5 -3.13978D+00,-2.16774D+00,-8.75973D-01,-2.69128D+00,-2.75079D-01, 6 -2.01595D+00,-7.27202D-01,-2.57463D+00,-2.82613D+00,-1.85926D+00, 7 -5.83100D-01,-2.43375D+00,-2.66775D+00,-1.70861D+00,-4.46177D-01, 8 -2.33538D+00,-2.54727D+00,-1.57696D+00,-3.20621D-01,-2.26061D+00/ DATA((pa(i,j),i=1,4),j=21,30) 1/-2.41621D+00,-1.45077D+00,-2.02005D-01,-2.19280D+00,-2.26783D+00, 2 -1.32260D+00,-8.70453D-02,-2.08509D+00,-2.11486D+00,-1.19361D+00, 3 2.46472D-02,-1.92532D+00,-1.96140D+00,-1.06555D+00, 1.32932D-01, 4 -1.75235D+00,-1.81119D+00,-9.40306D-01, 2.37547D-01,-1.59151D+00, 5 -1.66614D+00,-8.19000D-01, 3.38400D-01,-1.44695D+00,-1.52612D+00, 6 -7.02071D-01, 4.35573D-01,-1.31362D+00,-1.38980D+00,-5.89455D-01, 7 5.29286D-01,-1.18559D+00,-1.26058D+00,-4.81316D-01, 6.19801D-01, 8 -1.05824D+00,-1.14156D+00,-3.77677D-01, 7.07441D-01,-9.28641D-01/ DATA((pa(i,j),i=1,4),j=31,40) 1/-1.02607D+00,-2.77467D-01, 7.92633D-01,-7.96285D-01,-9.08964D-01, 2 -1.79687D-01, 8.75828D-01,-6.59336D-01,-7.92023D-01,-8.40258D-02, 3 9.57448D-01,-5.20172D-01,-6.86472D-01, 8.77490D-03, 1.03785D+00, 4 -3.76711D-01,-5.78202D-01, 1.00806D-01, 1.11766D+00,-2.29466D-01, 5 -4.72205D-01, 1.92075D-01, 1.19737D+00,-7.99781D-02,-3.68228D-01, 6 2.83239D-01, 1.27759D+00, 7.16105D-02,-2.57953D-01, 3.75953D-01, 7 1.35909D+00, 2.27488D-01,-1.39680D-01, 4.71236D-01, 1.44268D+00, 8 3.88037D-01,-1.35509D-02, 5.69997D-01, 1.52925D+00, 5.51614D-01/ DATA((pa(i,j),i=1,4),j=41,50) 1/ 1.20033D-01, 6.73231D-01, 1.61985D+00, 7.18516D-01, 2.61232D-01, 2 7.82111D-01, 1.71567D+00, 8.90290D-01, 4.10756D-01, 8.97990D-01, 3 1.81803D+00, 1.06860D+00, 5.69701D-01, 1.02234D+00, 1.92836D+00, 4 1.25489D+00, 7.39285D-01, 1.15661D+00, 2.04806D+00, 1.45031D+00, 5 9.20587D-01, 1.30208D+00, 2.17839D+00, 1.65531D+00, 1.11427D+00, 6 1.45958D+00, 2.31708D+00, 1.86928D+00, 1.32032D+00, 1.62926D+00, 7 2.47365D+00, 2.09059D+00, 1.53774D+00, 1.81034D+00, 2.63814D+00, 8 2.31684D+00, 1.76441D+00, 2.00105D+00, 2.81221D+00, 2.54529D+00/ DATA((pa(i,j),i=1,4),j=51,60) 1/ 1.99704D+00, 2.19870D+00, 2.99346D+00, 2.77324D+00, 2.23174D+00, 2 2.40071D+00, 3.18019D+00, 2.99898D+00, 2.46425D+00, 2.60309D+00, 3 3.36799D+00, 3.21993D+00, 2.69385D+00, 2.81083D+00, 3.56740D+00, 4 3.44225D+00, 2.90799D+00, 3.02078D+00, 3.77385D+00, 3.66470D+00, 5 3.12370D+00, 3.22283D+00, 3.97120D+00, 3.87615D+00, 3.33379D+00, 6 3.42204D+00, 4.16662D+00, 4.08289D+00, 3.53944D+00, 3.61896D+00, 7 4.36044D+00, 4.28587D+00, 3.74155D+00, 3.81401D+00, 4.55297D+00, 8 4.48582D+00, 3.94081D+00, 4.00754D+00, 4.74443D+00, 4.68335D+00/ 2000 FORMAT(10es10.3) 2001 FORMAT(2i3,5es12.5) c initialisations IF(init)THEN init=.FALSE. CALL bsp1dn(4,pa,psi,psit,60,m,knot,.FALSE.,psi(1),l,f,dfdx) IF(no_croiss)THEN PRINT*,'Arrêt 1 dans collision' ; STOP ENDIF CALL bsp1dn(4,pr,psi,psit,60,m,knot,.FALSE.,psi(1),l,f,dfdx) IF(no_croiss)THEN PRINT*,'Arrêt 2 dans collision' ; STOP ENDIF c parties constantes e4kpi=pi*(echarg**2/2.d0/kbol)**2*sqrt(kbol/2.d0/pi) !Paqu (66) ld0=kbol/4.d0/pi*amu/echarg**2 !Cox et al. p. 1192 a00=3.d0/4.d0/pi*amu xi00=4.d0*kbol/echarg**2 !Paquette (67) ENDIF c allocations ALLOCATE(dfxdx(4,nb),dld(nb),dgam(nb),gam0(nb,nb)) c calcul des gam0, les charges zb varient suivant l'ionisation DO i=1,nb DO j=1,nb gam0(i,j)=xi00/zb(i)/zb(j) ENDDO ENDDO c calcul du lambda ld=0.d0 a0=0.d0 dld(1:nb)=zb(1:nb)**2 ld=sum(xi(1:nb)*dld(1:nb)) !longueur de Debye a0=sum(xi(1:nb-1)) ld1=sqrt(ld0/ld/ro*t) dld(1:nb)=-ld1/ld*dld(1:nb)/2.d0 dld(1)=dld(1)-ld1*drox/2.d0 a01=(a00/ro/a0)**(1.d0/3.d0) IF(ld1 > a01)THEN lambda=ld1 ELSE lambda=a01 dld(1:nb-1)=-a01/a0/3.d0 dld(nb)=0.d0 !electrons dld(1)=dld(1)-a01*drox/3.d0 ENDIF DO i=1,nb DO j=i,nb gam=gam0(i,j)*t*lambda !Paquette 67 dgam(1:nb)=gam/lambda*dld(1:nb) bid1=2.d0*gam ; bid2=1.d0+gam**2 bid3=log(bid2) !Paquette 68 gam=log(bid3) !Paquette 68 bid=bid1/bid2/bid3 dgam(1:nb)=bid*dgam(1:nb) bid=exp(gam) IF(gam > 4.d0)THEN f(1)=log(1.00141d0*bid-3.18209d0) !Paquette psi > 4 f(2)=log(0.99559d0*bid-1.29553d0) f(3)=log(1.99814d0*bid-0.64413d0) f(4)=log(1.99016d0*bid-4.56958d0) dfxdx(1,:)=1.00141d0*bid/f(1)*dgam(:) dfxdx(2,:)=0.99559d0*bid/f(2)*dgam(:) dfxdx(3,:)=1.99814d0*bid/f(3)*dgam(:) dfxdx(4,:)=1.99016d0*bid/f(4)*dgam(:) ELSE IF(zb(i)*zb(j) > 0.d0)THEN !répulsif CALL bsp1dn(4,pr,psi,psit,60,m,knot,.TRUE.,gam,l,f,dfdx) ELSE !attractif CALL bsp1dn(4,pa,psi,psit,60,m,knot,.TRUE.,gam,l,f,dfdx) ENDIF DO k=1,nb dfxdx(:,k)=dfdx*dgam(k) ENDDO ENDIF f=exp(f) DO k=1,nb dfxdx(:,k)=f*dfxdx(:,k) ENDDO c WRITE(*,2001)i,j,gam,f(1),f(2),f(3),f(4) eps=e4kpi*(zb(i)*zb(j)/t)**2*sqrt(t/mij(i,j)) !Paquette 66 omega11(i,j)=eps*f(1) !Paquette 65 zij(i,j)=1.d0-2.d0/5.d0*f(2)/f(1) !Paquette 23 zpij(i,j)=2.5d0-2.d0/5.d0*(5.d0*f(2)-f(3))/f(1)!Paquette 24 zsij(i,j)=f(4)/f(1) !Paquette 25 domega11(i,j,1:nb)=eps*dfxdx(1,1:nb) dzij(i,j,1:nb)=-2.d0/5.d0/f(1)*(dfxdx(2,1:nb) 1 -dfxdx(1,1:nb)*f(2)/f(1)) dzpij(i,j,1:nb)=-2.d0/5.d0/f(1)*(5.d0*dfxdx(2,1:nb) 1 -dfxdx(3,1:nb)-(5.d0*f(2)-f(3))*dfxdx(1,1:nb)/f(1)) dzsij(i,j,1:nb)=zsij(i,j)*(dfxdx(4,1:nb)/f(4) 1 -dfxdx(1,1:nb)/f(1)) IF(j /= i)THEN !symétrisation omega11(j,i)=omega11(i,j) zij(j,i)=zij(i,j) ; zpij(j,i)=zpij(i,j) ; zsij(j,i)=zsij(i,j) domega11(j,i,1:nb)=domega11(i,j,1:nb) dzij(j,i,1:nb)=dzij(i,j,1:nb) ; dzpij(j,i,1:nb)=dzpij(i,j,1:nb) dzsij(j,i,1:nb)=dzsij(i,j,1:nb) ENDIF ENDDO ENDDO DEALLOCATE(dfxdx,dld,dgam,gam0) RETURN END SUBROUTINE collision