C C PURPOSE: GENERAL ROUTINE FOR COMPUTING COLLISIONAL RATES C IN THE SYSNTHESIS CODE FOR M1 LINES, THIS SHOULD BE USED C ONLY FOR DETERMINING IONIZATION BALANCE C C INPUTS: C C OUTPUTS: C C COMMON: C CATOM CATMOS CATMO2 CCONST CINPUT CLU C COMMENTS: October 6, 1999, P. JUDGE C C C********************************************************************* C SUBROUTINE GENCOL C C GENERAL ROUTINE FOR COMPUTING COLLISIONAL RATES (VERSION MAR '89) C FORMAT IS EITHER GIVEN IN TERMS OF A COLLISION STRENGTH OMEGA C OR A CE(TE) RATE C C OMEGAS ARE USED FOR CHARGED IONS IN GENERAL, C SINCE THEY REMAIN ROUGHLY CONSTANT WITH TEMPERATURE. C RATE IS PROP. TO 1/SQRT(TE) * EXP(DELTA-E) C C CE(TE) VALUES ARE USED LARGELY FOR NEUTRALS FOR SAME REASON. C RATE IS PROP. TO SQRT(TE) * EXP(DELTA-E) C C CI(TE) VALUES ARE USED FOR IONIZATION RATES: THESE HAVE A C SQRT(TE) * EXP(DELTA-E) DEPENDENCE AS FOR CE(TE) C C ORIGINALLY CODED BY P.G. JUDGE, APRIL 2ND, 1987 C MODIFIED FOR FORTRAN-77 STANDARD BY P.G. JUDGE, MAY 21ST 1987 C: C: GENCOL 87-05-21 NEW ROUTINE: (PHILIP JUDGE) C: GENERAL COLLISONAL ROUTINE READING TABLES C: C: 89-03-24 MODFICATIONS: (PHILIP JUDGE) C: NEW INPUT 'CP' AND 'CH' ADDED FOR INELASTIC COLLISIONS C: WITH NEUTRAL AND IONIZED HYDROGEN. C: 'SEMI ' CHANGED SO THAT THE ABSORPTION OSCILLATOR STRENGTH C: IS REQUIRED FOR INPUT. C: C: 89-06-07 MODIFICATIONS: (MATS CARLSSON) C: DO LOOP VARIABLE IT CHANGED TO JT IN ORDER NOT TO CONFLICT C: WITH THE MAIN ITERATION VARIABLE IT IN ITER C: C: 90-04-10 MODIFICATIONS: (PHILIP JUDGE) C: LOW TEMPERATURE DIELECTRONIC RECOMBINATION OPTION ADDED C: FOLLOWING NUSSBAUMER AND STOREY (A+A ) C: NOTE THAT ONLY THE RATE DOWNWARD IS ADDED TO THE FIXED RATE C: C: 90-08-27 MODIFICATIONS: (MATS CARLSSON) C: DO LOOP 1 REPLACED BY A GOTO TO REMOVE RESTRICTION ON C: NUMBER OF LINES WITH COLLISIONAL DATA C: C: 92-08-10 MODIFICATIONS (MATS CARLSSON) C: INTEGER ARRAYS WITH ELEMENT NAME CHANGED TO CHARACTER ARRAYS C: C: 92-09-07 MODIFICATIONS (MATS CARLSSON) C: ADDED OPTIONS RECOMB AND ALPHA C: C: 94-01-24 MODIFICATIONS (PHILIP JUDGE) C: AR85 OPTIONS ADDED. ARNAUD AND ROTHENFLUG 1985 C: 94-02-22 MODIFICATIONS (PHILIP JUDGE) C: AR85-CEA FORMALISM ADDED C: 94-03-08 MODIFICATIONS (PHILIP JUDGE) C: AR85-CH, AR85-CH+, AR85-CHE, AR85-CHE+, FORMALISM ADDED C INCLUDE 'PREC' PARAMETER (MTGRD=50,MSHELL=5) INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLU' C C LOCAL VARIABLES: UP TO MTGRD POINTS IN TEMPERATURE GRID ALLOWED C THIS CAN BE INCREASED C DIMENSION TGRID(MTGRD),CGRID(MTGRD) DIMENSION CDI(5,MSHELL) CHARACTER *20 KEY CNST=0.0 C C OUTPUT FOR IWATOM .GT. 1 C IF(IWATOM .GT. 1)WRITE(LOUT,800) 800 FORMAT('1 COLLISION RATES FROM SUBROUTINE G E N C O L'// * 1X,'PARAMETERS ARE AS FOLLOWS:'// * 1X,'TEMP - TEMPERATURE GRID GIVEN IN FILE ATOM'/ * 1X,' (SPLINE INTEPOLATIONS ARE MADE)'// * 1X,'CE - NEUTRAL B-B COLLISION RATE WHERE'/ * 1X,' R(U-L) = CE * NE * G(L) * SQRT(TEMP) / G(U)'// * 1X,'OHM - ION B-B COLLISION RATE WHERE'/ * 1X,' R(U-L) = OHM * NE * 8.63E-6 / SQRT(TEMP) / G(U)'// * 1X,'SEMI - B-B COLLISIONAL RATE USING THE SEMI-EMPIRICAL') IF(IWATOM .GT. 1) WRITE(LOUT,801) 801 FORMAT(1X,' FORMULAE OF VAN REGEMORTER AND SHEVELKO'// * 1X,'CP - B-B COLLISIONAL RATE WITH PROTONS WHERE'/ * 1X,' R(U-L) = CP * NP '// * 1X,'CH - B-B COLLISIONAL RATE WITH NEUTRAL HYDROGEN WHERE'/ * 1X,' R(U-L) = CH * N(H I) '// * 1X,'CI - B-F IONIZATION RATE WHERE'/ * 1X,' R(L-U) = CI * NE * EXP(-DELTAE/KT) * SQRT(T)'// * 1X,'CALP - COLLISIONS WITH ALPHA PARTICLES'/ * 1X,' R(I-J) = CALP * N(NK)'/// * 1X,'***********************************************'/ * 1X,'*** ALL THE ABOVE HAVE REVERSE RATES IN LTE ***'/ * 1X,'*** THE FOLLOWING HAVE ZERO REVERSE RATES ***'/ * 1X,'*** UNLESS SPECIFIED EXPLICITLY ***'/ * 1X,'***********************************************'// * 1X,'CH0 - CHARGE TRANSFER RATE WITH NEUTRAL HYDROGEN'/ * 1X,' R(I-J) = CH0 * NH(LEVEL1)'// * 1X,'CH+ - CHARGE TRANSFER RATE WITH IONIZED HYDROGEN'/ * 1X,' R(I-J) = CH+ * N(H II)'/// * 1X,'RECO - RECOMBINATION RATE'/ * 1X,' R(J-I) = RECO * NE'/// * 1X,'CONS - CONSTANT RATE'/ * 1X,' R(J-I) = CONST '/// * 1X,'KEY ','FROM TO PARAMETERS... '/) C C ASSUME INITIALLY THAT DATA ARE GIVEN INDEPENDENT OF TEMPERATURE: C I.E. THAT NTEMP IS ORIGINALLY SET = 1 C NTEMP=1 C C READ THE KEYWORD 'KEY' AND ASSOCIATED PARAMETERS C 50 CONTINUE READ(LDUMS,100,END=301)KEY 100 FORMAT(A) CALL LJUST(KEY) IF(KEY(1:1) .EQ. ' ') GOTO 50 IF (KEY(1:3) .EQ. 'END')THEN GOTO 300 ELSE IF (KEY(1:4) .EQ. 'TEMP')THEN READ(LDUMS,*)NTEMP,(TGRID(JT),JT=1,MIN(NTEMP,MTGRD)) IF(NTEMP .GT. MTGRD)THEN CALL STOP(' GENCOL: WORK ARRAYS (TGRID) TOO SMALL') ENDIF IF (IWATOM .GT. 1)THEN WRITE(LOUT,803)KEY(1:4),(TGRID(JT),JT=1,NTEMP) 803 FORMAT(1X,A,8X,10(1P,1X,E10.3)) ENDIF GOTO 50 ELSE IF (KEY(1:4) .EQ. 'SEMI') THEN READ(LDUMS, *) IL,IH, FAB IF(IWATOM .GT. 1)WRITE(LOUT,8004)KEY(1:4),IL,IH 8004 FORMAT(1X,A,2(2X,I2)) ELSE IF (KEY(1:4) .EQ. 'LTDR') THEN READ(LDUMS, *) IL,IH, ALTDR,BLTDR,CLTDR,DLTDR,FLTDR IF(IWATOM .GT. 1)WRITE(LOUT,8005)KEY(1:4),IL,IH,ALTDR, * BLTDR,CLTDR,DLTDR,FLTDR 8005 FORMAT(1X,A,2(2X,I2),5(2X,F11.4)) C C ADDITION P. JUDGE 24-JAN-1994: BEGIN C ELSE IF (KEY(1:7) .EQ. 'AR85-RR') THEN READ(LDUMS, *) IL,IH,ARADAR,ETAAR ELSE IF (KEY(1:8) .EQ. 'AR85-CDI') THEN READ(LDUMS, *) IL,IH,NCDI IF(NCDI .GT. MSHELL) CALL STOP('GENCOL: NCDI .GT. MSHELL') READ(LDUMS,*) ((CDI(J,I),J=1,5),I=1,NCDI) C C ADDITION P. JUDGE 24-JAN-1994: END C C: C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS START: C: ELSE IF (KEY(1:8) .EQ. 'AR85-CEA') THEN READ(LDUMS, *) IL,IH,CEAFAK ILO=MIN( IL, IH ) IHI=MAX( IL, IH ) CALL AR85CEA(ILO,IHI,CEA) CEA=CEA*CEAFAK C: C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS END: C: C: C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS START: C: ALL FOUR CHARGE TRANSFER KEYWORDS HERE ELSE IF (KEY(1:7) .EQ. 'AR85-CH') THEN READ(LDUMS, *) IL,IH READ(LDUMS,*) AR85T1,AR85T2,AR85A,AR85B,AR85C,AR85D C: C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS END: C: ELSE IF (KEY(1:7) .EQ. 'SHULL82') THEN READ(LDUMS, *) IL,IH,ACOLSH,TCOLSH,ARADSH,XRADSH,ADISH, * BDISH,T0SH,T1SH IF(IWATOM .GT. 1)WRITE(LOUT,8006)KEY(1:6),IL,IH,ACOLSH, * TCOLSH,ARADSH,XRADSH,ADISH,BDISH,T0SH,T1SH 8006 FORMAT(1X,A,2(2X,I2),8(2X,1P,E11.4)) ELSE IF (KEY(1:6) .EQ. 'CORONA') THEN READ(LDUMS, *) IL,IH IF(IWATOM .GT. 1)WRITE(LOUT,8004)KEY(1:6),IL,IH ELSE IF (KEY(1:8) .EQ. 'CONSTANT') THEN READ(LDUMS, *) IL,IH, CNST IF(IWATOM .GT. 1)WRITE(LOUT,8014)KEY(1:4),IL,IH, CT 8014 FORMAT(1X,A,2(2X,I2),2X,1P,E11.4) ELSE IF (KEY(1:3) .EQ. 'OHM' * .OR. KEY(1:2) .EQ. 'CE' * .OR. KEY(1:2) .EQ. 'CP' * .OR. KEY(1:2) .EQ. 'CH' * .OR. KEY(1:4) .EQ. 'CALP' * .OR. KEY(1:4) .EQ. 'RECO' * .OR. KEY(1:3) .EQ. 'CH0' * .OR. KEY(1:3) .EQ. 'CH+' * .OR. KEY(1:2) .EQ. 'CI' )THEN READ(LDUMS,*)IL,IH,(CGRID(JT),JT=1,NTEMP) NINTEP=NTEMP C: C: 23-FEB-1994 P.G.JUDGE MODIFICATIONS START: C: SPECIAL CASE OF RAPID TE DEPENDENCE- INTERPOLATE LINEARLY IF(KEY(1:3) .EQ. 'CH+' .OR. KEY(1:3) .EQ. 'CH0') THEN NINTEP=2 c WRITE(LJOBLO,*)'LINEAR INTERPOLATION OF LOG10 ',KEY(1:3) c * , ' BETWEEN LEVELS', IL, IH DO 333 JT=1,NTEMP CGRID(JT)=LOG10(CGRID(JT)) 333 CONTINUE ENDIF C: C: 23-FEB-1994 P.G.JUDGE MODIFICATIONS END: C: IF(IWATOM .GT. 1)THEN WRITE(LOUT,802)KEY(1:4), IL, IH, (CGRID(JT),JT=1,NTEMP) 802 FORMAT(1X,A,2(2X,I2),10(1P,1X,E10.3)) ENDIF ELSE WRITE(LJOBLO,191)KEY(1:4) 191 FORMAT(' GENCOL: KEYWORD(1:4) IS (',A,')') CALL STOP(' GENCOL: UNKNOWN KEYWORD IN ATOMIC DATA') ENDIF C C IDENTIFY THE UPPER AND LOWER LEVELS: C ILO=MIN( IL, IH ) IHI=MAX( IL, IH ) C C CHECK TO SEE THAT THE LEVELS ARE BETWEEN 1 AND NK C IF(ILO .LT. 1 .OR. IHI .GT. NK) * CALL STOP(' GENCOL: LEVEL INDEX OUTSIDE RANGE (1,NK)') C C CHECK TO SEE IF THERE ARE LEVELS WITH ENERGIES NOT INCREASING WITH C THE LEVEL LABEL (THIS AFFECTS CI ETC.) C IF(EV(IHI) .LT. EV(ILO) )THEN WRITE(LJOBLO,1001)KEY(1:4),IHI,ILO 1001 FORMAT(' GENCOL: WARNING- KEYWORD ',A,2(2X,I2), * ' ENERGY OF UPPER LEVEL IS .LT. LOWER LEVEL'/ * ' WARNING- CHECK RESULTS CAREFULLY') ENDIF C C DEPTH DO-LOOP: C C C SEMI-EMPIRICAL COLLISIONS C IF (KEY(1:4) .EQ. 'SEMI')THEN EUPCM=EV(IHI)*EE/CC/HH ELOCM=EV(ILO)*EE/CC/HH FEM = FAB*G(ILO)/G(IHI) CT=SEMIC(ION(ILO),EUPCM,ELOCM,FEM,TEMP,III) CDN=CT*NE CUP = CDN * NSTAR(IHI) / NSTAR(ILO) C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CORONAL ION BALANCE DATA (FROM DATA READ IN OPACITY PACKAGE) C ELSE IF (KEY(1:6) .EQ. 'CORONA')THEN PE=BK*NE*TEMP CALL CORONR(ATOMID(1:3),ION(ILO),TEMP,PE,RECRAT,RIRAT) C(IHI,ILO) = C(IHI,ILO) + RECRAT*NE C(ILO,IHI) = C(ILO,IHI) + RIRAT*NE ELSE IF(TEMP .LT. TGRID(1)) THEN CT=CGRID(1) ELSE IF (TEMP .GT. TGRID(NTEMP)) THEN CT=CGRID(NTEMP) ELSE CT=SPLIN(TEMP,TGRID,CGRID,NTEMP,NINTEP) ENDIF ENDIF IF(KEY(1:3) .EQ. 'CH+' .OR. KEY(1:3) .EQ. 'CH0') CT=10.**CT IF(CT .LT. 0.) THEN WRITE(LJOBLO,*)'INTERPOLATED COLLISION PARAMETER < 0.', * ' AT TEMPERATURE ', TEMP WRITE(LJOBLO,*)'FOR TRANSITION ', KEY,' ',IL,' ', IH WRITE(LJOBLO,*)(TGRID(JT),JT=1,NTEMP) WRITE(LJOBLO,*)(CGRID(JT),JT=1,NTEMP) CALL STOP('NEGATIVE COLLISION RATES ENCOUNTERED') ENDIF C C OMEGAS ARE GIVEN (+VE IONS) C IF(KEY(1:3) .EQ. 'OHM') THEN CDN = 8.63E-06 * CT * NE / ( G(IHI)*SQRT(TEMP) ) CUP = CDN * NSTAR(IHI) / NSTAR(ILO) C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CE VALUES ARE GIVEN (NEUTRALS) C ELSE IF (KEY(1:2) .EQ. 'CE')THEN CDN = NE * CT * G(ILO) * SQRT(TEMP) / G(IHI) CUP= CDN * NSTAR(IHI) / NSTAR(ILO) C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CP VALUES ARE GIVEN (B-B COLLISIONS WITH PROTONS) C ELSE IF (KEY(1:2) .EQ. 'CP')THEN CDN = NH(6) * CT CUP= CDN * NSTAR(IHI) / NSTAR(ILO) C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CH VALUES ARE GIVEN (COLLISIONS WITH NEUTRAL HYDROGEN) C ELSE IF (KEY(1:3) .EQ. 'CH ')THEN CDN = NH(1) * CT CUP= CDN * NSTAR(IHI) / NSTAR(ILO) C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CI VALUES ARE GIVEN C ELSE IF (KEY(1:2) .EQ. 'CI')THEN DEKT= (EV(IHI)-EV(ILO)) * EK / TEMP CUP = NE * CT * EXP(-DEKT) * SQRT(TEMP) CDN = CUP * NSTAR(ILO) / NSTAR(IHI) C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CH OR CH+ VALUES GIVEN: CHARGE TRANSFER COLLISIONS WITH C HYDROGEN ATOMS/IONS NB: ORDERING IMPORTANT HERE C EXAMPLE INPUT: C C CH0 C 1 2 1.E-9 1.E-9 1.E-9 C C THIS ASSUMES THAT THE RATE 1.E-9*NH(1) WILL BE ADDED TO THE C COLLISION RATE FROM LEVEL 1 TO LEVEL 2 C ELSE IF (KEY(1:3) .EQ. 'CH0') THEN C(IL,IH)= NH(1) * CT + C(IL,IH) ELSE IF (KEY(1:3) .EQ. 'CH+') THEN C(IL,IH)= NH(6) * CT + C(IL,IH) C C LOW-TEMPERATURE DIELECTRONIC RECOMBINATION COEFFS. ARE GIVEN C FORMULA (9) OF NUSSBAUMER AND STOREY (A+A PAPER II) C ELSE IF(KEY(1:4) .EQ. 'LTDR') THEN T4LTDR=TEMP/1.E4 CDN = 1.E-12*NE* (ALTDR/T4LTDR + BLTDR + * CLTDR*T4LTDR +DLTDR*T4LTDR*T4LTDR) * * EXP(-FLTDR/T4LTDR)/T4LTDR/SQRT(T4LTDR) C(IHI,ILO) = CDN + C(IHI,ILO) C C ADDITION P. JUDGE 24-JAN-1994: BEGIN C C PARAMETRIC FITS TO DATA OF ARNAUD AND ROTHENFLUG 1985, C APJS 60, 425. C CDI IS 5 BY MSHELL ELEMENT ARRAY, CONTAINING PARAMETERS C I, A,B,C,D OF TABLE 1 OF ARNAUD AND ROTHENFLUG ELSE IF(KEY(1:7) .EQ. 'AR85-RR') THEN CDN= ARADAR*(TEMP/1.E4)**ETAAR C(IHI,ILO) = CDN + C(IHI,ILO) ELSE IF(KEY(1:8) .EQ. 'AR85-CDI') THEN CUP=0. DO 321 J=1,NCDI XJ=CDI(1,J)*EE/BK/TEMP FXJ=EXP(-XJ)*SQRT(XJ) *(CDI(2,J)+CDI(3,J)*(1.+XJ) + * (CDI(4,J)-XJ*(CDI(2,J)+CDI(3,J)*(2.+XJ)))*FONE(XJ) + * CDI(5,J)*XJ*FTWO(XJ)) CUP = CUP+ 6.69E-07 / CDI(1,J) / SQRT(CDI(1,J))*FXJ 321 CONTINUE IF(CUP .LT. 0.) CALL STOP('GENCOL: CDI-CUP .LT.0') C(ILO,IHI) = CUP*NE + C(ILO,IHI) C(IHI,ILO) = CUP*NE*NSTAR(ILO)/NSTAR(IHI) * + C(IHI,ILO) C C ADDITION P. JUDGE 24-JAN-1994: END C C: C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS START: C: ELSE IF(KEY(1:8) .EQ. 'AR85-CEA') THEN C(ILO,IHI) = CEA*NE + C(ILO,IHI) C C THIS IS INCORRECT SINCE POPULATION OF UPPER LEVEL IS THE UPPER LEVEL C OF THE DOUBLY EXCITED STATE, **, I.E. THIS SHOULD BE MULTIPLIED BY C G(**)/G(IHI) EXP(E(**) - E(IHI)). A BETTER APPROXIMATION WOULD SEEM TO C BE 0. TIMES THIS TO AVOID PROBLEMS. HENCE COMMENT THIS LINE OUT C C C(IHI,ILO) = CEA*NE*NSTAR(ILO)/NSTAR(IHI) C * + C(IHI,ILO) C: C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS END: C: C: C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS START: C: C: AR85-CH CHARGE TRANSFER RECOMBINATION WITH NEUTRAL HYDROGEN C: ELSE IF(KEY(1:7) .EQ. 'AR85-CH' .AND. KEY(1:8) .NE. * 'AR85-CH+' .AND. KEY(1:8) .NE. 'AR85-CHE') THEN IF(TEMP .GE. AR85T1 .AND. TEMP .LE. AR85T2) THEN T4=TEMP/1.E4 CUP = AR85A * 1.E-9 * T4**AR85B * (1. + AR85C * * EXP(AR85D*T4))*NH(1) C(IL,IH) = CUP + C(IL,IH) ENDIF C: C: AR85-CH+ CHARGE TRANSFER WITH IONIZED HYDROGEN C: ELSE IF(KEY(1:8) .EQ. 'AR85-CH+') THEN IF(TEMP .GE. AR85T1 .AND. TEMP .LE. AR85T2) THEN T4=TEMP/1.E4 CUP = AR85A * 1.E-9 * T4**AR85B * EXP(-AR85C*T4) * * EXP(-AR85D*EE/BK/TEMP)*NH(6) C(IL,IH) = CUP + C(IL,IH) ENDIF C: C: AR85-CHE CHARGE TRANSFER WITH NEUTRAL HELIUM C: C: C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS END: C: C C COMPUTATIONS FROM TABLES OF SHULL & VAN STEENBERG, C AP J SUPPL 48 P95 C ELSE IF(KEY(1:7) .EQ. 'SHULL82') THEN CDN=ARADSH*(TEMP/1.E4)**(-XRADSH) + * ADISH /TEMP/SQRT(TEMP) * EXP(-T0SH/TEMP)* * (1.E0+BDISH * (EXP(- T1SH/TEMP))) CDN=CDN*NE CUP=ACOLSH * SQRT(TEMP) * EXP( -TCOLSH / TEMP) * / (1.E0 + 0.1 * TEMP / TCOLSH) CUP=CUP*NE C* C* 3-BODY RECOMBINATION (HIGH DENSITY LIMIT) C* CDN=CDN+CUP*NSTAR(ILO)/NSTAR(IHI) C* C(IHI,ILO) = CDN + C(IHI,ILO) C(ILO,IHI) = CUP + C(ILO,IHI) C C CALP VALUES ARE GIVEN (COLLISIONS WITH ALPHA PARTICLES) C ATOM HAS TO BE HE C INITIA IS CALLED FOR FIRST RATE TO READ POPULATIONS C FROM RSTRT C C RECOMBINATION COEFFICIENTS ARE GIVEN C C(J-I)=NE*CT C ELSE IF (KEY(1:4) .EQ. 'RECO') THEN C(IHI,ILO)= NE*CT + C(IHI,ILO) C C SIMPLE ADDITION OF CONSTANT, FOR EXAMPLE FOR PHOTOIONIZATION C C(J-I)=CT C ELSE IF (KEY(1:8) .EQ. 'CONSTANT') THEN C(IL,IH)= CNST + C(IL,IH) C C MORE KEYWORDS CAN BE ADDED HERE C ENDIF GOTO 50 C 301 CALL STOP(' GENCOL: KEYWORD ''END'' NOT FOUND IN COLLISION FILE') 300 CONTINUE RETURN END C C*********************************************************************** C FUNCTION SEMIC(Z,EUP,ELO,FEM,TEMP,IFLAG) C C RETURNS COLLISION RATE DOWNWARDS USING SEMI-EMPIRICAL GF'S C C P.G. JUDGE, NOVEMBER 1987. C C INPUT: C Z - CHARGE OF ION+1 (E.G. FOR C II Z=2) C EUP - ENERGY OF UPPER LEVEL IN CM-1 C ELO - ENERGY OF LOWER LEVEL IN CM-1 C FEM - EMISSION OSCILLATOR STRENGTH ( = GF/G(UPPER)) C TEMP - TEMPERATURE IN DEGREES KELVIN C C OUTPUT: C IFLAG - = 1 FOR VAN REGEMORTER APPROXIMATION C IFLAG - = 2 FOR SHEVELKO APPROXIMATION C SEMIC- DOWNWARD COLLISION RATE IN CM3 /S C INCLUDE 'PREC' INTEGER Z DATA C1,RY /1.43882,109737.312/ C C INITIALIZE C DELTE=EUP-ELO BETA=DELTE*C1/TEMP SEMIC=0.0 C C IFLAG=1 => VAN REGEMORTER APPROXN. C IF(BETA .GE. 0.01)THEN IFLAG=1 SEMIC = 3.2E-7*FEM*(RY/DELTE)**1.5*SQRT(BETA)*P(Z,BETA) C C IFLAG=2 => SHEVELKO APPROX C ELSE IF(BETA .LT. 0.01)THEN IFLAG=2 SEMIC=1.74E-07*FEM*EXP(BETA)/SQRT(BETA*DELTE/RY) + *LOG( ELO / DELTE * SQRT(RY / C1 / TEMP)) ENDIF C RETURN END C C ********************************************************************** C FUNCTION P(Z,B) C C THERMAL P-FUNCTION FOR SEMIC EMPIRICAL COLLISION RATES C REFERENCE: SOBEL'MAN - 'ATOMIC PHYSICS II. ' C INCLUDE 'PREC' PARAMETER (N1=10) DIMENSION BETREF(N1),PNREF(N1),PCREF(N1) INTEGER Z DATA BETREF/-2.0,-1.699,-1.398,-1.0,-0.699,-0.398,0.0 + ,0.301,0.602,1.0/ DATA PNREF/1.160,0.956,0.758,0.493,0.331,0.209 + ,0.100,0.063,0.040,0.023/ DATA PCREF/1.160,0.977,0.788,0.554,0.403 + ,0.290,0.214,0.201,0.200,0.200/ C C LIMIT OF HIGH TEMPERATURE: KT >> E C IF (B .LT. 0.01)THEN P = 0.27566 * (0.577 + LOG(B)) RETURN C C INTERMEDIATE TEMPS (MOST IMPORTANT FOR EQM PLASMAS) C SPLINE INTERPOLATION ONTO LOGB GRID: C ELSE IF(B .GT. 0.01 .AND. B .LT. 10.0)THEN BB=LOG10(B) IF (Z .EQ. 1)P = SPLIN(BB,BETREF,PNREF,N1,N1) IF (Z .GT. 1)P = SPLIN(BB,BETREF,PCREF,N1,N1) RETURN C C LIMIT OF LOW TEMPERATURE: KT << E C ELSE IF(B .GT. 10.0)THEN IF(Z .EQ. 1)P = 0.066 / SQRT(B) IF(Z .GT. 1)P = 0.200 RETURN ENDIF END C C ********************************************************************** C SUBROUTINE CORONR(CEL,ION,TEMP,PRESS,RRATE,RIRATE) C C CORONAL ION BALANCE CALCULATION C NEW ROUTINE FOR USE WITH V2.0 OPACITY PACKAGE AND GENCOL C C INPUT: ELEMENT IDENTIFIER CEL (CHARACTER) C ION STAGE IDENTIFIER ION (INTEGER) C TEMP (ELECTRON TEMPERATURE IN K (RL) C PRESS (ELECTRON PRESSURE IN DYNE/CM2) (RL) C C OUTPUT: RRATE RECOMBINATION RATE (CM3/S) C RIRATE IONIZATION RATE (CM3/S) C C EXAMPLE: C CEL='C ', ION=1, WILL RETURN C THE RECOMBINATION RATE PER C II ATOM PER ELECTRON C AND THE IONIZATION RATE PER C I ATOM PER ELECTRON C C NOTES: C COMPUTATIONS FROM TABLES OF SHULL & VAN STEENBERG, C AP J SUPPL 48 P95 C NO PE-DEPENDENCE, THIS COULD BE ADDED WHEN COMPUTATIONS BECOME C AVAILABLE C IF THE ION FRACTION IS NOT FOUND THEN THE ROUTINE WILL STOP C C 89-03-28 NEW ROUTINE (PHILIP JUDGE) C C: CORONR 92-08-10 MODIFICATIONS (MATS CARLSSON) C: INTEGER ARRAYS WITH ELEMENT NAME CHANGED TO CHARACTER ARRAYS C: C: 92-09-09 MODIFICATIONS (MATS CARLSSON) C: INTEGER IEL IN WRITE STATEMENT CHANGED TO CHARACTER C: INCLUDE 'PREC' PARAMETER (MCOR=30) INTEGER NCOR,IONCOR(MCOR) DOUBLE PRECISION ARADC(MCOR),XRADC(MCOR),ADIC(MCOR),T0C(MCOR), * BDIC(MCOR),T1C(MCOR),ACOLC(MCOR),TCOLC(MCOR),TECOR,DECOR COMMON /CCORON/ TECOR,DECOR, * ARADC,XRADC,ADIC,T0C,BDIC,T1C,ACOLC,TCOLC,NCOR,IONCOR CHARACTER*3 CCOR COMMON /CCCOR/ CCOR(MCOR) COMMON /CLU/ LINPUT,LATOM,LATOM2,LATMOS,LDSCAL,LABUND,LOUT, * LTIME,LRSTRT,LDSCA2,LWMAT,LNIIT,LDUMS,LDUMI,LDUMC,LOPC, * LJNY,LINIT,LPHI,LJOBLO,LATHSE CHARACTER*(*) CEL C C IDENTIFY THE ELEMENT AND IONIZATION STAGE TO BE COMPUTED C DO 100 N=1,NCOR IF( CEL .EQ. CCOR(N) .AND. ION .EQ. IONCOR(N) )GOTO 99 100 CONTINUE WRITE(LJOBLO,1024) CEL,ION 1024 FORMAT(' CORONR: SPECIES',A3,I3,' NOT FOUND IN ABSDAT FILE') CALL STOP(' ') C C ELEMENT HAS BEEN COMPUTED- NOW RETURN THE RATES C 99 RRATE=ARADC(N)*(TEMP/1.E4)**(-XRADC(N)) + * ADIC(N) * (TEMP**(-1.5)) * EXP(-T0C(N)/TEMP)* * (1.E0+BDIC(N) * (EXP(- T1C(N)/TEMP))) RIRATE=ACOLC(N) * SQRT(TEMP) * EXP( -TCOLC(N) / TEMP) * / (1.E0 + 0.1 * TEMP / TCOLC(N)) RETURN END