SUBROUTINE COLIS(ID,T,COL,CLOC) C =============================== C C Driving procedure for evaluation of the collisional rates C C Input: T - temperature C Output: COL - array of quantities proportional to the C collisional rates in all transitions C for a given temperature (ie. at a given depth) C Precisely, COL(IT)*Nlow(IT) is the UPWARD C collisional rate of the transition IT C Output: CLOC - array of quantities proportional to the C collisional rates in all transitions C for a given temperature (ie. at a given depth) C Precisely, CLOC(IT)*Nupp(IT) is the DOWNWARD C collisional rate of the transition IT C C Procedure COLIS calls procedures COLH and COLHE for evaluating C the collisional rates in hydrogen and helium, C and itself evaluates collisional rates for other species C C Evaluation is controlled by input parameter ICOL(ITR). C Meaning of ICOL for all species, excluding hydrogen and helium: C C a) for ionization C ICOL = 0 - Seatons formula ; here the value of the photo- C ionization cross section at the threshold is C transmitted in array OSC0 C = 1 - Allen's formula; again, OSC0 has the meaning of C the necessary multiplicative parameter C = 2 - the so-called SIMPLE1 mode - see below C = 3 - the so-called SIMPLE2 mode - see below C b) for excitation C ICOL = 0 - Van Regemorter formula, with standard g(bar)=0.25 C = 1 - Van Regemorter formula, with "exact" g(bar) C = 2 - the so-called SIMPLE1 mode - see below C = 3 - the so-called SIMPLE2 mode - see below C = 4 - Eissner-Seaton formula - see below C C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ODFPAR.FOR' PARAMETER (EXPIA1=-0.57721566,EXPIA2=0.99999193, * EXPIA3=-0.24991055,EXPIA4=0.05519968, * EXPIA5=-0.00976004,EXPIA6=0.00107857, * EXPIB1=0.2677734343,EXPIB2=8.6347608925, * EXPIB3=18.059016973,EXPIB4=8.5733287401, * EXPIC1=3.9584969228,EXPIC2=21.0996530827, * EXPIC3=25.6329561486,EXPIC4=9.5733223454) DIMENSION COL(MTRANS),TYPEARR(MXTCOL),CCRATE(MCFIT), * CCTEMP(MCFIT),CLOC(MTRANS) C COMMON/CTRTEMP/ te C CREGER(X,U,A,GG)=19.7363*X*EXP(-U)/U*GG*A CSEATN(X,U,A)=1.55D13*X/ABS(U)*EXP(-U)*A CALLEN(X,U,A)=X*A*EXP(-U)/U/U CSMPL1(X,U,A)=5.465D-11*X*EXP(-U)*A CSMPL2(X,U,A)=5.465D-11*X*EXP(-U)*A*(1.+U) CUPSX(X,U,A)=8.631D-6/X*EXP(-U)*A C DO I=1,NTRANS CLOC(I)=0. COL(I)=0. END DO C C calculate collider's populations: e, p, H(1s) C ANE : Electron C ANP : Proton C ANHM: H- C ANH : H(1s) ANE=ELEC(ID) c IF (IELH.gt.0) THEN ! if H is an explicit atom ANP=POPUL(NNEXT(IELH),ID) ! Protons ANH=POPUL(NFIRST(IELH),ID) ! H(1s) else anh=ahtot-anp ENDIF IF (IELHM.gt.0) THEN ANHM=POPUL(NFIRST(IELHM),ID) ! if H- is an explicit atom else anhm=1.0353e-16/t/sqrt(t)*exp(8762.9/t)*anh*ane ENDIF C HKT=HK/T SRT=SQRT(T) T32=UN/T/SRT TK=HKT/H CSTD=0.25 T0=TEFF if(t0.gt.0.) then TT0=UN-T/T0 SRT0=SQRT(T0)/SRT end if C C Call procedures COLH and COLHE for hydrogen and helium C IF(IATH.NE.0) CALL COLH(ID,T,COL) IF(IATHE.NE.0) CALL COLHE(T,COL) IF(IATH.NE.0.OR.IATHE.NE.0) THEN DO I=1, NTRANS COL(I)=COL(I) * ANE IF (LINE(I)) THEN CLOC(I)=COL(I) * EXP(FR0(I)*HKT)*G(ILOW(I))/G(IUP(I)) ELSE CORR=UN NKE=NNEXT(IEL(ILOW(I))) IF(NKE.NE.IUP(I)) CORR=G(NKE)/G(IUP(I))* * EXP((ENION(NKE)-ENION(IUP(I)))*TK) CLOC(I)=COL(I) * ANE * SBF(ILOW(I))*CORR ENDIF ENDDO ENDIF C C Loop over all explicit species, excluding hydrogen and helium C DO 100 IAT=1,NATOM IF(IAT.EQ.IATH.OR.IAT.EQ.IATHE) GO TO 100 N0I=N0A(IAT) NKI=NKA(IAT) DO 50 I=N0I,NKI-1 IE=IEL(I) DO 40 J=I+1,NKI IT=ITRA(I,J) IF(IT.EQ.0) GO TO 40 IC=ICOL(IT) COL(IT)=0.0 CLOC(IT)=0.0 C1=OSC0(IT) C2=CPAR(IT) U0=FR0(IT)*HKT U0HM=U0-8752.072/T ! including H-minus potential ! U0P =U0-157821.5/T ! including H-proton potential! DO K=1, MXTCOL TYPEARR(K)=0 ENDDO IF(LINE(IT)) GO TO 30 C C the detailed balancing factor for an inverse process C CORR=UN NKE=NNEXT(IEL(I)) IF(NKE.NE.J) CORR=G(NKE)/G(J)* * EXP((ENION(NKE)-ENION(J))*TK) CINV=ANE*SBF(I)*CORR C C *********** tabulated data *************** C For collisional processes that change the total C charge of the target atom, there are three C processes considered here: C - TYPE 1 Electron Collisional ionization C - TYPE 2 Charge exchange with protons C - TYPE 3 Charge exchange with hydrogen C C There are several 'calculated' options in the code for C electron collisional excitation that are C neglected if TYPE 1 tabulated data are present. C C There is also an option for TYPE 2 that is calculated C here (ICOL ge 10) that is also overriden if TYPE 2 is C present. IF (IC.GE.1000) THEN IORICE=1 IF (IC.LT.0) IORICE=-1 IC=ABS(IC) ITYPE=IC/1000 IC=IORICE*(MOD(IC,1000)-1) !ICOL RECOVERED DO K=1, MXTCOL TYPEARR(K)=MOD(ITYPE/(2**(K-1)),2) ENDDO C C ****** START 'GENCOL' FOR IONIZATION ********** C C ****** ELECTRON COLLISIONAL IONIZATION ***** IF (TYPEARR(1).EQ.1) THEN NX=0 DO K=1,MCFIT IF (CTEMP(1,K,IT).NE.0) THEN CCRATE(K)=log(CRATE(1,K,IT)) CCTEMP(K)=CTEMP(1,k,IT) NX=NX+1 ELSE ! CLEAN CC**** ARRAYS CCRATE(K)=0. CCTEMP(K)=0. ENDIF ENDDO cs=ylintp(t,cctemp,ccrate,nx,mcfit) CS=ANE*exp(CS) COL(IT)=COL(IT) + CS !UPPWARD CLOC(IT)=CLOC(IT) + CS*CINV !DOWNWARD ENDIF C ****** CHARGE EXCHANGE WITH PROTONS ******* IF (TYPEARR(2).EQ.1) THEN NX=0 DO K=1,MCFIT IF (CTEMP(2,K,IT).NE.0) THEN CCRATE(K)=log(CRATE(2,K,IT)) CCTEMP(K)=CTEMP(2,K,IT) NX=NX+1 ELSE ! CLEAN CC**** ARRAYS CCRATE(K)=0. CCTEMP(K)=0. ENDIF ENDDO cs=ylintp(t,cctemp,ccrate,nx,mcfit) cs=exp(cs)*ANP COL(IT)=COL(IT) + CS ! UPPWARD CINH=G(I)/G(J)*0.5*EXP(U0P) CLOC(IT)=CLOC(IT) + CS*CINH*ANH !DOWNWARD ENDIF C ******* CHARGE EXCHANGE WITH HYDROGEN ***** IF (TYPEARR(3).EQ.1) THEN NX=0 DO K=1,MCFIT IF (CTEMP(3,K,IT).NE.0) THEN CCRATE(K)=log(CRATE(3,K,IT)) CCTEMP(K)=CTEMP(3,K,IT) NX=NX+1 ELSE ! CLEAN CC**** ARRAYS CCRATE(K)=0. CCTEMP(K)=0. ENDIF ENDDO cs=ylintp(t,cctemp,ccrate,nx,mcfit) cs=exp(cs)*ANH COL(IT)=COL(IT) + CS ! UPPWARD CINH=G(I)/G(J)*TWO*EXP(U0HM) CLOC(IT)=CLOC(IT) + CS*CINH*ANHM !DOWNWARD ENDIF C ************** END GENCOL ****************** IF (IC.EQ.-1) GO TO 40 ENDIF C C ********* Charge transfer with protons reactions- an old scheme C IF(IC.GE.10) THEN IF(TYPEARR(2).NE.1) THEN C radiative charge transfer ionization of neutrals in the C grround state with protons te=T CS=HCTIon(1,NUMAT(IAT)) CS=CS*ANP COL(IT)=COL(IT) + CS !UPPWARD CS=CS*0.5*G(I)/G(J) * EXP(U0P) CLOC(IT)=CLOC(IT) + CS*ANH !DOWNWARD ENDIF IC=IC-10 IF(TYPEARR(1).eq.1) GO TO 40 ENDIF C C ********* Electron collisional ionization C IF(IC.EQ.0) THEN CS=CSEATN(UN/SRT,U0,C1)*ANE COL(IT)=COL(IT)+CS ELSE IF(IC.EQ.1) THEN CS=CALLEN(T32,U0,C1)*ANE COL(IT)=COL(IT)+CS ELSE IF(IC.EQ.2) THEN CS=CSMPL1(SRT,U0,C1)*ANE COL(IT)=COL(IT)+CS ELSE IF(IC.EQ.3) THEN CS=CSMPL2(SRT,U0,C1)*ANE COL(IT)=COL(IT)+CS ELSE IF(IC.EQ.4) THEN ia=numat(iatm(i)) CS=cion(ia,iz(iel(i)),enion(i)*6.24298e11,t)*ANE COL(IT)=COL(IT)+CS ELSE IF(IC.EQ.5) THEN ia=numat(iatm(i)) izc=iz(ie) rno=16. ii=i-nfirst(ie)+1 call irc(ii,t,izc,rno,cs) CS=CS*ANE col(it)=COL(IT)+CS ELSE IF(IC.LT.0) THEN CALL CSPEC(I,J,IC,C1,C2,U0,T,CS) CS=CS*ANE COL(IT)=COL(IT)+CS END IF CLOC(IT)=CLOC(IT)+CS*CINV !DOWNWARD if(ic.eq.4) go to 40 C C collisional excitations from level I to higher, non-explicit C levels are lumped into the collisional ionization rate C (the so-called modified collision ionization rate) C N0Q=NQUANT(NLAST(IE))+1 N1Q=ICUP(IE) IF(N1Q.EQ.0) GO TO 40 IQ=NQUANT(I) REL=G(I)/2./IQ/IQ DO 20 JQ=N0Q,N1Q XJ=JQ U0=(ENION(I)-EH/XJ/XJ)*TK IF(JQ.LE.20) CC1=OSH(IQ,JQ)*REL IF(JQ.GT.20) CC1=OSH(IQ,20)*(20./XJ)**3*REL GG=CSTD if(u0.gt.35.) go to 20 IF(U0.LE.UN) THEN EXPIU0=-LOG(U0)+EXPIA1+U0*(EXPIA2+ * U0*(EXPIA3+U0*(EXPIA4+ * U0*(EXPIA5+U0*EXPIA6)))) ELSE EXPIU0=EXP(-U0)*((EXPIB1+U0*(EXPIB2+ * U0*(EXPIB3+ * U0*(EXPIB4+U0))))/(EXPIC1+U0*(EXPIC2+ * U0*(EXPIC3+U0*(EXPIC4+U0)))))/U0 END IF GG0=0.276*EXP(U0)*EXPIU0 IF(GG0.GT.GG) GG=GG0 CS=CREGER(T32,U0,CC1,GG)*ANE COL(IT)=COL(IT)+CS !UPPWARD CLOC(IT)=CLOC(IT)+CS*ANE*SBF(I)*WOP(I,ID)*CORR !DOWNWARD 20 CONTINUE GO TO 40 C C ********* Collisional excitation C 30 CONTINUE C C the detailed balancing factor for an inverse process C CINV=EXP(U0)*G(I)/G(J) C c ********** Tabulated Data IF (IC.GE.1000) THEN IORICE=1 IF (IC.LT.0) IORICE=-1 IC=ABS(IC) ITYPE=IC/1000 IC=IORICE*(MOD(IC,1000)-1) !ICOL RECOVERED DO K=1, MXTCOL TYPEARR(K)=MOD(ITYPE/(2**(K-1)),2) END DO C ****** START 'GENCOL' FOR EXCITATION ********** C C ****** ELECTRON COLLISIONAL EXCITATION ***** IF (TYPEARR(1).EQ.1) THEN NX=0 DO K=1,MCFIT IF (CTEMP(1,K,IT).NE.0) THEN CCRATE(K)=log(CRATE(1,K,IT)) CCTEMP(K)=CTEMP(1,K,IT) NX=NX+1 ELSE !CLEAN CC**** ARRAYS CCRATE(K)=0. CCTEMP(K)=0. ENDIF END DO cs=ylintp(t,cctemp,ccrate,nx,mcfit) CS=exp(CS)*ANE COL(IT)=COL(IT) + CS ! UPPWARD CLOC(IT)=CLOC(IT) + CS*CINV ! DOWNWARD END IF C ****** PROTON COLLISIONAL EXCITATION ******* IF (TYPEARR(2).EQ.1) THEN NX=0 DO K=1,MCFIT IF (CTEMP(2,K,IT).NE.0) THEN CCRATE(K)=log(CRATE(2,K,IT)) CCTEMP(K)=CTEMP(2,K,IT) NX=NX+1 ELSE !CLEAN CC**** ARRAYS CCRATE(K)=0. CCTEMP(K)=0. ENDIF END DO cs=ylintp(t,cctemp,ccrate,nx,mcfit) CS=exp(CS)*ANP COL(IT)=COL(IT) + CS ! UPPWARD CLOC(IT)=CLOC(IT) + CS*CINV ! DOWNWARD END IF C ******* HYDROGEN COLLISIONAL EXCITATION ***** IF (TYPEARR(3).EQ.1) THEN NX=0 DO K=1,MCFIT IF (CTEMP(3,K,IT).NE.0) THEN CCRATE(K)=log(CRATE(3,K,IT)) CCTEMP(K)=CTEMP(3,K,IT) NX=NX+1 ELSE !CLEAN CC**** ARRAYS CCRATE(K)=0. CCTEMP(K)=0. ENDIF END DO cs=ylintp(t,cctemp,ccrate,nx,mcfit) CS=exp(CS)*ANH COL(IT)=COL(IT) + CS ! UPPWARD CLOC(IT)=CLOC(IT) + CS*CINV ! DOWNWARD END IF C ************** END GENCOL ****************** IF (IC.EQ.-1) GO TO 40 END IF IF(IC.LE.1.AND.IC.GT.0) THEN GG=CSTD IF(IC.EQ.1) GG=C2 IF(U0.LE.UN) THEN EXPIU0=-LOG(U0)+EXPIA1+ * U0*(EXPIA2+U0*(EXPIA3+U0*(EXPIA4+ * U0*(EXPIA5+U0*EXPIA6)))) ELSE EXPIU0=EXP(-U0)*((EXPIB1+U0* * (EXPIB2+U0*(EXPIB3+ * U0*(EXPIB4+U0))))/(EXPIC1+U0*(EXPIC2+ * U0*(EXPIC3+U0*(EXPIC4+U0)))))/U0 END IF GG0=0.276*EXP(U0)*EXPIU0 IF(GG0.GT.GG) GG=GG0 CS=CREGER(T32,U0,C1,GG)*ANE COL(IT)=COL(IT) + CS ! UPPWARD ELSE IF(IC.EQ.2) THEN CS=CSMPL1(SRT,U0,C1*C2)*ANE COL(IT)=COL(IT) + CS ! UPPWARD ELSE IF(IC.EQ.3) THEN CS=CSMPL2(SRT,U0,C2)*ANE COL(IT)=COL(IT) + CS ! UPPWARD ELSE IF(IC.EQ.4) THEN CS=CUPSX(SRT,U0,C2/G(I))*ANE COL(IT)=COL(IT) + CS ! UPPWARD ELSE IF(IC.EQ.9) THEN CS=OMECOL(I,J)*SRT0*EXP(-U0*TT0) * ANE COL(IT)=COL(IT) + CS ! UPPWARD ELSE IF(IC.LT.0) THEN CALL CSPEC(I,J,IC,C1,C2,U0,T,CS) CS=CS*ANE COL(IT)=COL(IT) + CS ! UPPWARD END IF CLOC(IT)=CLOC(IT) + CS*CINV ! DOWNWARD 40 CONTINUE 50 CONTINUE 100 CONTINUE RETURN END