C N. R. BADNELL v3.21 08 DEC 2022 C C PROGRAM ADASEX/J C C----------------------------------------------------------------------- C C THIS CODE IS FOR USE WITH LSRM/ICFT/BPRM/DARC & AS-DW/PWB RUNS. C C IT READS AN R-MATRIX OMEGA/U FILE OF COLLISION STRENGTHS, C CONVOLUTES/AVERAGES THEM AND OUTPUTS TO AN ADAS ADF04 FILE. C DEFAULT IS A MAXWELL CONVOLUTION FOR A TYPE-3 UPSILON FILE. C OPTIONAL KAPPA/DRUYVESTEYN CONVOLUTION FOR A TYPE-4 DOWN/UPSILON FILE. C OPTIONAL INTERVAL-AVERAGE FOR A BINNED TYPE-5 ADF04 FILE. C OPTIONAL READ OF TYPE-5 ADF04 FILE (RM-BINNED OR AS-DW/PWB) AND C CONVOLUTION TO TYPE-3/4. (NO R-MATRIX OMEGA/U FILE USED THEN.) C C THE CONVOLUTION USES A 2-POINT QUADRATURE FORMULA THRU THE C TABULATED ENERGY RANGE (SO THERE IS NO MESH RESTRICTION). C IF INFINITE ENERGY POINTS ARE PRESENT IT USES A 2-POINT INTERPOLATION C IN BURGESS-TULLY SPACE TO EXTEND THE ENERGY MESH, C FOR ALLOWED TRANSITIONS: DIPOLE AND (IF PRESENT) BORN C AND EXTRAPOLATION FOR FORBIDDEN TRANSITIONS. C C THE USER *MUST* FLAG IF BORN LIMITS ARE *NOT* PRESENT, BY SETTING C IBORN = 0 WHEN LS ARE NOT GOOD QUANTUM NUMBERS (IC/BP) OR -2 (LS/JK). C IF NO DIPOLE PRESENT, ONLY THE EXPLICITLY TABULATED COLLISION STRENGTH C MAY BE INTEGRATED. THERE IS NO OTHER ATTEMPT AT FITTING/EXTRAPOLATION. C C----------------------------------------------------------------------- C PROGRAM MAIN C IMPLICIT REAL*8(A-H,O-Z) C REAL*4 OMEGIN REAL*4 TARRY(2),TIME C REAL*8 MAXE,KAPPA C INTEGER*8 MXD1,MXD2,MXD3,MXDM C !MAX DIMENSIONS PARAMETER (NMPTS=250000) !ENERGY POINTS PARAMETER (NMPIN=3000) !IN BATCHES OF AT LEAST 500 PARAMETER (NMBIN=501) !TOTAL NO. BIN ENERGIES PARAMETER (NSTMX=5000) !STATES/LEVELS PARAMETER (NTRMX=NSTMX*(NSTMX+1)/2) !TRANSITIONS PARAMETER (NMTMP=20) !OUTPUT TEMPERATURES PARAMETER (NLNMX=500) !NO. OF COMMENT LINES C PARAMETER(MXD1=NMPIN/NTRMX,MXD2=NTRMX/NMPIN,MXD3=MXD1+MXD2, X MXDM=NMPIN*MXD1/MXD3+NTRMX*MXD2/MXD3+1) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DTEN=10.0D0) PARAMETER (D0PT1=0.1D0) PARAMETER (D0PT2=0.2D0) PARAMETER (D1PT5=1.5D0) PARAMETER (D2PT5=2.5D0) PARAMETER (D1P10=1.0D10) PARAMETER (D1M4=1.0D-4) PARAMETER (D1M30=1.0D-30) !ADAS ZERO C PARAMETER (DMIN=1.0D-99) !ADF04 MIN PARAMETER (DMAX=9.99D99) !ADF04 MAX C PARAMETER (DFSC=DONE/137.03599976D0) PARAMETER (DALF=DFSC*DFSC) C PARAMETER (DKMIN=1.5D0+1.D-10) !MIN KAPPA PARAMETER (DKMAX=1.0D10) !MAX KAPPA PARAMETER (DXMIN=0.1D0) !MIN XDRY PARAMETER (DXMAX=1.0D3) !MAX XDRY C PARAMETER (DELIN0=0.025D0) !*TAZ = DEFAULT LIN ENERGY STEP PARAMETER (DELOG0=1.08D0) !DEFAULT LOG ENERGY STEP FACTOR C PARAMETER (ABVTHR=1.0D-2) PARAMETER (ETHRSH0=1.D-10) PARAMETER (TOLE=1.D+4) !COULD ALLOW ADASEXJ ENERG PARAMETER (TOLH0=1.D-10) C PARAMETER (CONRCM=1.097373D+05) !RYDBERGS TO CM^-1 PARAMETER (CONRYK=1.578885D5) !RYDBERGS TO KELVIN PARAMETER (CONRYEV=13.6058) !RYDBERGS TO EV PARAMETER (CONRAD=2.008123d9) !ALPHA**3/4*HBAR C LOGICAL BBORN,BEQUAL,BEXP,BFORM,BINF,BLAST,BSHORT,BREAD,BOMTHR X ,BREV,BRMX,BDW,BPWB,btop,b37 C CHARACTER IONTRM*2,IEL*2 CHARACTER CONF(NSTMX)*15,TRM*15,XCONF(NSTMX)*31 CHARACTER CMANT(NMBIN+2)*5,CEXP(NMBIN+2)*3 CHARACTER NAME*30,DATE*30,DATE8*8,COMMENT(NLNMX)*80 CHARACTER LBCD(10)*1,SBCD(10)*1,ISBCD*1,LTBCD*1 C CHARACTER(LEN=1) CAR1 CHARACTER(LEN=8) CANT(21),CBLNK8 CHARACTER(LEN=19) C300 CHARACTER(LEN=41) F542 CHARACTER(LEN=41) F543 CHARACTER(LEN=33) F570 CHARACTER(LEN=49) F580 CHARACTER(LEN=49) F581 CHARACTER(LEN=29) F600 CHARACTER(LEN=43) F761 CHARACTER(LEN=35) F762 CHARACTER(LEN=200) CARD C ALLOCATABLE AR(:,:),UPS(:,:,:),OMEGIN(:,:) C ALLOCATABLE :: TEMPE37(:),E37(:,:),F37(:,:) C DIMENSION EMSH(0:NMPTS),EJ(0:NMPIN),NOMRD(NMPTS),OMEGA(0:MXDM) C DIMENSION ENAT(NSTMX),ENERG(NSTMX),FJT(NSTMX),IZERO(NSTMX) X ,JXTWO(NSTMX),LAT(NSTMX),ISAT(NSTMX),INDEX(NSTMX) X ,ITOT(NSTMX),LIN(NTRMX),LFN(NTRMX) C DIMENSION TKADAS(NMTMP),FREL(NMTMP) C DIMENSION EBIN(0:NMBIN),XMANT(0:NMBIN+1),IEXP(0:NMBIN+1) X ,XB(0:NMBIN) C C HISTORIC NAMELIST NAMED ADASEX TO MIRROR ADASEX CODE, AND UNCHANGED C WHEN ADASEXJ DEVELOPED: C NAMELIST/ADASEX/ BEQUAL,BEXP,BOMTHR,BREV,BDW,BPWB,btop X,EASYM,EFITF,FIPOT,eshft,XBTEST X,IBORN,IEL,IEINCR,IONE,IONTRM,IRMPS,ITCC,IREL !,ITYPE X,MXEIN,MXTMP X,NDELBT,NSTATES,NUMTMP X,NLOWMN,NLOWMX,NUPMN,NUPMX X,NBIN,EMIN,EMAX,MAXE,DELIN,NLIN !DEFINE BINNED ENERGY MESH X,ICON,KAPPA,XDRY,IDUP1 !DEFINE NON-MAXWELLIAN C X,EFITH,EFITL,GFMIN,IELAS,ILOW,IPTOMG,IRATES !NO LONGER IN USE X,NLEVS,NTERM,NTERMS !EQUIVALENT TO NSTATES X,IRDTMP !NO LONGER NEEDED, SEE NUMTMP C C ALLOW USE OF NAMELIST ADASEXJ TO MIRROR ADASEXJ CODE, DROP HISTORIC C (NO LONGER USED/NEEDED) OPTIONS: C NAMELIST/ADASEXJ/BEQUAL,BEXP,BOMTHR,BREV,BDW,BPWB,btop X,EASYM,EFITF,FIPOT,eshft,XBTEST X,IBORN,IEL,IEINCR,IONE,IONTRM,IRMPS,ITCC,IREL !,ITYPE X,MXEIN,MXTMP X,NDELBT,NSTATES,NUMTMP X,NLOWMN,NLOWMX,NUPMN,NUPMX X,NBIN,EMIN,EMAX,MAXE,DELIN,NLIN !DEFINE BINNED ENERGY MESH X,ICON,KAPPA,XDRY,IDUP1 !DEFINE NON-MAXWELLIAN C DATA EINF/1.0D6/ !NORMAL Z-SCALED INF E DATA NTADAS/13/ !NUMBER ADAS TEMPS DATA (TKADAS(I),I=1,13)/2.0D+02,5.0D+02,1.0D+03,2.0D+03,5.0D+03, X1.0D+04,2.0D+04,5.0D+04,1.0D+05,2.0D+05,5.0D+05,1.0D+06,2.0D+06/ DATA (LBCD(I),I=1,10)/"S","P","D","F","G","H","I","K","L","M"/ DATA (SBCD(I),I=1,7)/'1','2','3','4','5','6','7'/ DATA IFLAGB/0/ DATA CBLNK8/' '/ C EQUIVALENCE(NTERM,NTERMS,NLEVS,NSTATES) !ALLOW OLD NAMES FOR NOW C PI=ACOS(-DONE) C F600='(1PE14.8,6E11.3/(14X,6E11.3))' !SUPPRESS IFORT REMARK C C SOME HIGH LEVEL LOGICAL SWITCHES ("INTERNAL" CANNOT BE SET BY USER). C BEXP=.FALSE. !TRUE=1.0E+00, FALSE=1.0+00 BINF=.FALSE. !INTERNAL: INF E FLAG BRMX=.FALSE. !INTERNAL: R-MX ELSE DW/PWB BDW=.FALSE. !USER MUST SET .T. IF THRESHOLD MISSING E.G.NEUTRAL BPWB=.FALSE. !SET TRUE ONLY IF NON TYPE 1 BOMTHR=.TRUE. !R-MX OMEGAS START "AT" THRESHOLDS BREV=.FALSE. !REVERSE INDEX ORDER ON R-MX OUTPUT adf04 btop=.false.!if(bdw)compare allowed omegas with Born limit & reset XBTEST=5 !APPLY btop ONLY FOR X.GT.XBTEST lrglam=30 !default AS-RMX top-up L/J, for logging only C C UNIT NOS IN USE: C C 5 - adasex/j.in STANDARD INPUT C 6 - adasex/j.out STANDARD OUTPUT C C 7 - omega/u/OMEGA/U COLLISION STRENGTH INPUT C OR C 8 - adf04_om RM-BINNED/AS-DW/PWB COLLISION STRENGTH INPUT C C 9 - adf04_om RM-BINNED COLLISION STRENGTH OUTPUT C AND/OR C 9 - adf04_ups EFFECTIVE COLLISION STRENGTH OUTPUT C C 11 - adf37 NUMERICAL ENERGY DISTRIBUTION FUNCTION (ICON=3) C C 98 - adf04xxx SCRATCH I/O OPEN(98,STATUS='SCRATCH',FORM='FORMATTED') C 99 - omega SCRATCH I/O C INQUIRE (FILE='adasexj.in',EXIST=BFORM) IF(BFORM)THEN OPEN(5,FILE='adasexj.in') OPEN(6,FILE='adasexj.out') ITCC=1 ELSE INQUIRE (FILE='adasex.in',EXIST=BFORM) IF(BFORM)THEN OPEN(5,FILE='adasex.in') OPEN(6,FILE='adasex.out') ITCC=-1 ELSE STOP '*** ERROR: NEITHER adasexj.in NOR adasex.in FOUND!!!' ENDIF ENDIF C C QUICK LOOK-UP TO SEE WHICH FILES TO OPEN (NSTATES ONLY) C NBIN=0 NSTATES=0 !COMPULSORY C NAMLIN=1 READ(5,ADASEX,END=1) GO TO 2 1 REWIND(5) NAMLIN=2 READ(5,ADASEXJ) 2 CONTINUE C REWIND(5) !GET REST OF INPUT LATER C IF(NSTATES.EQ.0)THEN WRITE(6,*) '*** ERROR: NSTATES MUST BE SPECIFIED IN adasexj.in' WRITE(6,*) ' SET NSTATES .GT. 0 TO AVERAGE FIRST NSTATES OF' X ,' AN R-MATRIX OMEGA FILE' WRITE(6,*) ' SET NSTATES .LT. 0 (ANY VALUE) TO CONVOLUTE A' X ,' TYPE-5 ADF04 FILE adf04_om' WRITE(0,*) '*** ERROR: NSTATES MUST BE SPECIFIED IN adasexj.in' X ,' - SEE adasexj.out FOR DETAILS' STOP ENDIF C BREAD=NSTATES.LT.0 C IF(BREAD)THEN !READ AN EXISTING TYPE 5 adf04 FILE C INQUIRE (FILE='adf04_om',EXIST=BFORM) !INFILE IF(BFORM)THEN OPEN(8,FILE='adf04_om',FORM='FORMATTED') ELSE WRITE(6,*)'***ERROR: NSTATES.LT.0 BUT adf04_om FILE NOT FOUND' STOP '***ERROR: NSTATES.LT.0 BUT adf04_om FILE NOT FOUND' ENDIF C IREAD=8 !GET LEVEL INFO FROM adf04_om ***IGNORE ANY ON UNIT5*** C READ(IREAD,571)IEL,IQ,NZED,IQ1,FIPOT,IONTRM !GET HEADER C !SINCE BYPASS omega SET-UP NELC=NZED-IQ NOMWRT=0 NOMT=0 C OPEN(9,FILE='adf04_ups',FORM='FORMATTED') !OUTFILE C ELSE !READ AN R-MATRIX "OMEGA" FILE C INQUIRE (FILE='omega',EXIST=BFORM) !INFILE IF(BFORM)THEN WRITE(6,"(/'*** USING FORMATTED omega FILE')") open(7,file='omega') OPEN(99,STATUS='SCRATCH',FORM='FORMATTED') ELSE INQUIRE (FILE='omegau',EXIST=BFORM) IF(BFORM)THEN WRITE(6,"(/'*** USING UNFORMATTED omegau FILE')") open(7,file='omegau',FORM='UNFORMATTED') BFORM=.FALSE. ELSE INQUIRE (FILE='OMEGA',EXIST=BFORM) IF(BFORM)THEN WRITE(6,"(/'*** USING FORMATTED OMEGA FILE')") open(7,file='OMEGA') OPEN(99,STATUS='SCRATCH',FORM='FORMATTED') ELSE INQUIRE (FILE='OMEGAU',EXIST=BFORM) IF(BFORM)THEN WRITE(6,"(/'*** USING UNFORMATTED OMEGAU FILE')") open(7,file='OMEGAU',FORM='UNFORMATTED') BFORM=.FALSE. ELSE WRITE(6,"(/'*** ERROR: NO SUITABLE OMEGA FILE FOUND')") STOP '*** ERROR: NO SUITABLE OMEGA FILE FOUND' ENDIF ENDIF ENDIF ENDIF C BRMX=.TRUE. C IF(BFORM)READ(7,*)NZED,NELC IF(.NOT.BFORM)READ(7)NZED,NELC C IQ=NZED-NELC !omega SCALED CHARGE IQ1=IQ+1 !adf04 SCALED CHARGE IREAD=5 !USER SUPPLIES CONFIG/LEVEL INFO ON adasex/j.in C IF(NBIN.EQ.0)THEN !OUTFILE(S) OPEN(9,FILE='adf04_ups') !SOLE ELSE OPEN(9,FILE='adf04_om') !.LT.0 WILL DO _ups AFTER ENDIF C ENDIF C TAZ=IQ*IQ IF(TAZ.EQ.DZERO)TAZ=DONE FIQ2=IQ1*IQ1 IONE=IQ IF(IONE.GT.0)IONE=1 C C----------------------------------------------------------------------- C C READ INPUT FROM NAMELIST DATA: C C (ANY ENERGIES INPUT ARE ASSUMED TO BE IN UNSCALED RYDBERGS IF POSITIVE C AND STGF Z-SCALED IF NEGATIVE (& THEN RE-SET POSITIVE!) RECALL THAT C STGF Z-SCALING IS ONE LESS THAN ADAS adf04 (TEMPERATURE) Z-SCALING, C NEUTRALS EXCEPTED. ANY TEMPERATURES INPUT HERE ARE UNSCALED KELVIN.) C C C THE FIRST IS *COMPULSORY* C C NSTATES -- IF .GT. 0 THEN OF THE NAST TERMS/LEVELS INCLUDED IN THE C omega FILE THE FIRST NSTATES OF THEM WILL HAVE AVERAGED C TRANSITION DATA CALCULATED AND TERM/LEVEL INFORMATION C WRITTEN TO THE adf04_om/ups FILE. C CORRESPONDINGL, THERE MUST BE *EXACTLY* NTERMS/NLEVS C OF CONFIG/TERM/LEVEL INFO SPECIFIED IN adasex/j.in. C (HISTORIC NLEVS (JK/IC) OR NTERMS (LS) ARE RECOGNIZED.) C -- IF .LT. 0 THEN AN EXISTING TYPE 5 adf04_om FILE IS READ C AND CONVERTED TO TYPE 3/4 adf04_ups. CURRENTLY, THE C VALUE OF NSTATES DOES NOT MATTER, BUT FOR FUTURE SAFETY C SET ITS MAGNITUDE .GT. THE NO. OF STATES IN adf04_om: C THE CONFIG/TERM/LEVEL INFO IS TAKEN FROM adf04_om, ANY C USER DATA IN adasex/j.in IS CURRENTLY IGNORED (STC). C *NO DEFAULT* (INITIALIZED TO ZERO STOPS) C C C THE NEXT 3 SHOULD BE SPECIFIED IF READING AN omega FILE, C FOR AN adf04 FILE TO BE READILY USABLE BY ADAS. C IF READING AN adf04_om FILE THEN ANY SPECIFIED HERE WILL *OVERWRITE*. C C IEL -- SYMBOL FOR ELEMENT - TWO CHARACTERS C C FIPOT -- IONIZATION POTENTIAL IN CM-1 C C IONTRM -- TERM SPECIFICATION FOR RESULTANT ION - TWO CHARACTERS C C C --- OPTIONAL --- C C THE ELECTRON DISTRIBUTION IS MAXWELLIAN BY DEFAULT (ICON=0), BUT C C KAPPA -- IF .GT. 1.5 THEN USES A KAPPA-DISTRIBUTION (ICON=1). C C XDRY -- IF .GE. 1.0 USES A DRUYVESTEYN DISTRIBUTION (ICON=2). C C ICON -- IF .EQ. 3 THEN READ AN adf37 NUMERICAL DISTRIBUTION. C DEFAULT: -1 IS RESET TO MAXWELLIAN, KAPPA OR XDRY, C I.E. USER SHOULD ONLY SET FOR .EQ. 3. C OTHER VALUES ARE RESERVED FOR POSSIBLE FUTURE USE. C C IREL -- IF .NE. 0 THEN APPLY RELATIVISTIC (JUTTNER) CORRECTION C TO EFFECTIVE COLLISION STRENGTHS. C DEFAULT: 0 (WE LEAVE THIS TO ADAS TO APPLY) C C ITCC -- IF .GT. 0 THEN ASSUMES INPUT omega FROM AN IC CALCULATN C E.G. ICFT, BREIT-PAULI OR DARC. ENERGIES ARE TAKEN FROM C THE omega FILE. C -- IF .EQ. 0 THEN ASSUMES INPUT FROM A PURE JK-COUPLING C CALCULATION WITH STGICF. LEVEL ENERGIES ARE TAKEN FROM C adasexj.in SINCE omega ENERGIES ARE DEGENERATE WITHIN A C TERM. C -- IF .LT. 0 THEN ASSUMES INPUT omega FROM AN LS CALCULATN C AND ENERGIES ARE TAKEN FROM THE omega FILE. C adasex.in IS READ FOR CONFIGURATION INFO ONLY. C DEFAULT: 1 (IC adasexj.in) OR -1 (LS adasex.in) C C NLOWMN, NLOWMX, NUPMN, NUPMX -- CAN BE USED TO RESTRICT OUTPUT C TO A SUBSET OF TRANSITIONS FROM LOWER STATES BETWEEN C NLOWMN & NLOWMX TO UPPER STATES BETWEEN NUPMN & NUPMX. C DEFAULT: ALL PRESENT IN omega/adf04_om. C C IONE -- .EQ. 0 IF omega/adf04_om CONTAIN ELASTIC TRANSITIONS C .EQ. 1 IF THEY DON'T CONTAIN ELASTIC TRANSITIONS C DEFAULT: 0 (NEUTRALS) OR 1 (IONS). C C NUMTMP -- NUMBER OF TEMPERATURES AT WHICH EFFECTIVE COLLISION C STRENGTHS SHOULD BE CALCULATED. C IF .GT. 0 THEY ARE TAKEN FROM THE DEFAULT ADAS VALUES. C IF .LT. 0 THEN READ-IN -NUMTMP TEMPERATURES IN C DEGREES KELVIN DIRECTLY FOLLOWING THE NAMELIST LINE(S). C DEFAULT: 10 C C MXTMP -- IF EQUAL TO ZERO, THEN DO NOT LET THE TEMPERATURES IN C THE RUN EXCEED 1/2 MAXIMUM ENERGY IN THE OMEGA FILE; IF C NOT EQUAL TO ZERO, NO SUCH LIMIT. BUT A WARNING WILL BE C PRINTED IF THE TEMPERATURE EXCEEDS THE MAXIMUM ENERGY. C DEFAULT: 1 C C IEINCR -- IS THE ENERGY INDEX INCREMENT IN THE QUADRATURE FOR THE C EFFECTIVE COLLISION STRENGTHS. THUS, IEINCR .GT.1 SKIPS C IEINCR-1 MESH POINTS BETWEEN EACH MESH POINT IN THE C FINITE ENERGY REGION. THIS OPTION CAN BE USED TO C GENERATE adf04 FILES WITH EFFECTIVE COLLISION STRENGTHS C CALCULATED USING DIFFERENT NUMBERS OF MESH POINTS. C BY COMPARING SUCH FILES, ONE CAN INVESTIGATE C CONVERGENCE WITH THE DENSITY OF MESH POINTS, IF NOT C ALREADY DONE SO BY BUILDING AN "INTERLACED" omega FILE. C DEFAULT: 1 (ALL ENERGIES USED) C C NDELBT -- USE THE BURGESS-TULLY METHOD FOR THE FIRST NDELBT C QUADRATURE STEPS. THEREAFTER, USE THE TRAPEZOIDAL RULE. C N.B. BOTH TREAT THE COLLISION STRENGTH AS LINEAR STILL. C .LT. 0 IS RESET TO ALL. C DEFAULT: 20 IF RAW COLLISION STRENGTH (I.E. NON-BINNED) C -1 (I.E. ALL) IF CONVOLUTING BINNED DATA, C EXCEPT 2 IF DRUYVESTEYN DISTRIB. C N.B. BINNED OMEGA IS CONSTANT, INSTEAD OF LINEAR (RAW) C AND SO TRAPEZOIDAL GIVES NO SPEED-UP IF CAN INTEGRATE C DISTRIBUTION ANALYTICALLY (I.E. MAXWELLIAN AND KAPPA). C C OPTIONS RELATING TO HANDLING HIGH ENERGIES (CONVOLUTED): C C EASYM -- THE 2-POINT HIGH ENERGY INTERPOLATION IN BURGESS-TULLY C SPACE USES THE INFINITE ENERGY POINT AND THE LAST C FINITE ENERGY POINT. IF FOR SOME REASON THE LAST FINITE C ENERGY POINTS ARE UNREPRESENTATIVE, SET EASYM TO BE THE C MAXIMUM "SAFE" ENERGY. C ALL POINTS ARE STILL INCLUDED IN THE FINITE ENERGY C INTEGRATION THOUGH. SO, THEY SHOULD BE EXCLUDED PRIOR C TO RUNNING ADASEXJ IF THEY ARE THAT "BAD". C DEFAULT: USES THE LAST FINITE ENERGY POINT. C C IBORN -- IF .GT. 0 THEN BORN LIMITS ARE PRESENT AND SO ZERO C LIMITS CORRESPOND TO FORBIDDEN TRANSITIONS, WHICH MAYBE C EXTRAPOLATED AS 1/E**|IBORN| - *** SEE ALSO EFITF ***. C -- IF .LE. 0 THEN BORN LIMITS ARE NOT PRESENT. THEN C IF ITCC .GT. 0 (IC)THEN TRANSITIONS WITH ZERO LIMIT ARE C TAKEN TO BE ALLOWED AND EXTRAPOLATED AS 1/E**|IBORN| C (RECOMMEND IBORN=0 THEN!) C IF ITCC .LE. 0 (LS/JK)THEN TRANSITIONS WITH ZERO LIMITS C ARE TAKEN TO BE ALLOWED IF SPIN-ALLOWED AND SO ARE C EXTRAPOLATED AS A CONSTANT, AND ARE TAKEN TO BE C FORBIDDEN IF SPIN-CHANGING AND SO ARE EXTRAPOLATED AS C 1/E**|IBORN|. (RECOMMEND IBORN=-2 THEN.) C DEFAULT: 2 (BORN LIMITS PRESENT.) C C EFITF -- ONLY USE POINTS ABOVE EFITF (UNSCALED RYDBERGS) TO C DETERMINE THE HIGH ENERGY BEHAVIOUR (POWER LAW) OF C FORBIDDEN TRANSITIONS. (BOUNDED BY 1/E**1 AND 1/E**3.) C .EQ.0, INTERNALLY SET TO 2/3 MAX SCATTERING ENERGY. C N.B. EFITF=MAX(EFITF,ENAT(NAST)) IS ENFORCED. C 1/E**IBORN CAN BE FORCED BY EFIT LARGE. SEE ALSO IRMPS. C DEFAULT: 0 C C IRMPS -- IF .NE. 0 THEN THE INPUT omega FILE IS FROM AN RMPS C CALCULATION AND EFITF=MAX(EFITF,ENAT(NSTATES)) IS C ENFORCED INSTEAD SINCE ENAT(NAST) CAN BE VERY LARGE. C IF .LT. 0 THEN 1/E**IBORN BEHAVIOUR IS FORCED. THIS C CAN BE USED FOR RMPS OR NON-RMPS CALCULATIONS AS C IRMPS HAS NO OTHER EFFECT. C DEFAULT:0 C C MXEIN -- MAX NUMBER OF ENEGIES READ-IN AT A TIME FROM omega. C SUFFICIENT ENERGIES ABOVE ENAT(NAST) SHOULD BE TREATED C TOGETHER IN THE FINAL BATCH TO ENABLE EFITF USAGE. C DEFAULT: CHOSEN SO AS TO LEAVE AT LEAST 500 POINTS IN C THE LAST BATCH. C C C OPTIONS RELATED TO BINNED COLLISION STRENGTHS: C (BINNED ENERGY MESH IS ALSO USED AS THE INTERPOLATION MESH FOR DW/PWB) C C NBIN -- .GT. 0 NUMBER OF LOGARITHMIC SPACED BIN ENERGIES. C .LT. 0 LOOP BACK-UP AND CONVOLUTE BINNED (NBIN=-NBIN). C IF |NBIN| .LT. 10 THEN INTERNALLY RESET TO ~100. C (EXACT VALUE DEPENDS ON THE WIDTH OF RESONANCE REGION.) C DEFAULT: 0 (NO BINNED DATA) C C EMIN,MAX-- ENERGY RANGE OF NBIN. C DEFAULT: TKADAS(1)/10, ENAT(NSTATES)*TAZ (UNSCALED RYD) C C NLIN -- NUMBER OF LINEARLY SPACED BIN ENERGIES TO FOLLOW LOG. C DEFAULT: INITIALLY, ALL ALLOWED BY DIMENSION. THIS IS C REDUCED SUBSEQUENTLY WHEN THE LAST ENERGY IS C KNOWN FROM THE omega FILE. C C DELIN -- LINEAR STEP. C DEFAULT: 0.025*TAZ (UNSCALED RYD). C C MAXE -- FINAL BIN ENERGY. ONLY USED IF DELIN.EQ.0. C DEFAULT: 0 (NOT USED THEN) C C N.B. THE LINEAR MESH IS SET UPON UTILIZING TWO OUT OF NLIN,DELIN,MAXE. C C BEQUAL -- .TRUE. FORCES THE SAME MAX FINAL ENERGY FOR ALL *FINAL* C STATES. SINCE THE omega FILE DEFINES A COMMON MAXIMUM C *INITIAL* ENERGY, THIS DROPS PROGRESSIVELY MORE DATA C FOR LOWER AND LOWER FINAL STATES. THE NUMBER OF BINNED C COLLISION STRENGTHS IS THE SAME FOR ALL TRANSITIONS C AND IS EQUAL TO THE NUMBER OF NON-ZERO BIN ENERGIES. C .FALSE. KEEPS ALL omega DATA FOR ALL TRANSITIONS. C THE NUMBER OF BINNED COLLISION STRENGTHS DECREASES C WITH INCREASING FINAL STATE. THE NUMBER PRESENT FOR C THE HIGHEST UPPER STATE IS THE COMMON NUMBER (FINAL C ENERGY RANGE) WRITTEN FOR ALL TRANSITIONS BY .TRUE. C ***DEFAULT: .FALSE. (HIGH-T CONVOLUTED UPSILONS THEN MOST C CLOSELY MATCH THOSE FROM CONVOLUTING omega DIRECTLY.) C C----------------------------------------------------------------------- C C THE FOLLOWING ARE NO LONGER EXPLICITLY USED, THEY ARE PRESENT IN THE C NAMELIST FOR BACKWARDS COMPATIBILTY. THEY ARE EITHER IGNORED - THOSE C (LOW-LEVEL) VARIABLES COMMENTED-OUT - OR THE CODE WILL STOP IF USER C ATTEMPTS TO RE-SET HIGH-LEVEL ONES. C ILOW=1 !NO LEAST SQUARES FIT IRATES=0 !NO RATES IPTOMG=0 !NO DEBUG PRINT IELAS=1 !USE IONE C GFMIN=-1 !RE-INSTATE TO FLAG WEAK DIPOLE TO ALTERNATE C-VALUE? C EFITL=DZERO !USE EFITF, AS ONLY FORBIDDEN "FIT" C EFITH=D1P10 !SUBSUMED INTO EASYM C C----------------------------------------------------------------------- C C SET DEFAULTS FOR POSSIBLE USER INPUT: C C NSTATES=0 !COMPULSORY C IF(.NOT.BREAD)THEN !DO NOT OVERWRITE IEL=' ' !SHOULD BE SET FIPOT=DZERO !SHOULD BE SET IONTRM=' ' !SHOULD BE SET ENDIF C COLD ITYPE=3 !INTERNAL ONLY NOW IREL=0 ICON=-1 !FOR NUMERICAL DISTRIB (TBD) KAPPA=-DONE XDRY=-DONE IDUP1=1 !=2 SUPPRESS DOWNSILON C ITCC=-1 OR 1 !DEPENDING ON ADASEX/J NLOWMN=1 NLOWMX=NSTMX NUPMN=1 NUPMX=NSTMX NUMTMP=9999 !DEFAULT RE-SET TO 10 MXTMP=1 IRDTMP=0 !CAN USE NUMTMP.LT.0 TO FLAG READ IEINCR=1 IF(BREAD.OR.NBIN.NE.0)THEN !COARSE MESH NDELBT=2 ELSE !FINE MESH NDELBT=20 ENDIF EASYM=D1P10 IBORN=2 EFITF=DZERO IRMPS=0 MXEIN=-1 eshft=dzero C C BINNED ENERGY MESH C C NBIN=0 !NO. LOG ENERGIES EMIN=0 !FROM EMAX=0 !TO NLIN=0 !NO. LINEAR ENERGIES DELIN=0 !OR LINEAR STEP MAXE=0 !TO BEQUAL=.FALSE. !FORCE SAME NO BINS FOR ALL TRANSITIONS C C IF(NAMLIN.EQ.1)READ(5,ADASEX) IF(NAMLIN.EQ.2)READ(5,ADASEXJ) C C IF(IRATES.NE.0) X STOP 'RATES OUTPUT OPTION NOT PRESENT IN THIS VERSION OF ADASEXJ' C IF(ILOW.NE.1) X STOP 'LEAST SQUARES FIT NOT PRESENT IN THIS VERSION OF ADASEXJ' C IF(IELAS.LE.0) X STOP ' IELAS NOT USED, SET IONE = 1/0 TO FLAG IN/ELASTIC' C IF(IPTOMG.NE.0) XWRITE(0,*)'NO FIT/EXTRP OMEGA OUTPUT SINCE WE INTERP IN B-T SPACE' C IF(NSTATES.EQ.0) X STOP 'NSTATES MUST BE SET BY USER (AND BE .NE. ZERO!)' C c IF(BREAD.AND.NBIN.NE.0)STOP'CANNOT READ *AND* WRITE BINNED FILES' C C C SET INTERNAL ITYPE SWITCH, WHICH WILL ULTIMATELY LABEL C THE OUTPUT adf04 FILE. C C ITYPE -- THE TYPE OF ADF04 FILE TO BE WRITTEN: C .EQ. 3 FOR STANDARD MAXWELLIAN UPSILON FILE C .EQ. 4 FOR NON-MAXWELLIAN DISTRIBUTION C .EQ. 5 FOR INTERVAL AVERAGED COLLISION STRENGTHS C .EQ.-5 DITTO, AND CONVOLUTES BINNED OMEGAS WITH THE C ELECTRON DISTRIBUTION FUNCTION. C IF(NBIN.EQ.0)THEN ITYPE=3 ELSEIF(NBIN.GT.0)THEN ITYPE=5 ELSE ITYPE=-5 NBIN=-NBIN ENDIF C ITYPE0=ITYPE ITYPE=IABS(ITYPE) C ICON0=ICON IF(ICON0.LT.0)THEN IF(KAPPA.GE.DKMIN)THEN ICON=1 ELSEIF(XDRY.GE.DXMIN)THEN ICON=2 ELSE ICON=0 ENDIF ELSEIF(ICON0.EQ.1)THEN IF(KAPPA.LT.DKMIN) X STOP 'MUST SET A VALID KAPPA VALUE FOR THE K-DISTRIBUTION' ELSEIF(ICON0.EQ.2)THEN IF(XDRY.LT.DXMIN) STOP X 'MUST SET A VALID XDRY VALUE FOR THE DRUYVESTEYN DISTRIBUTION' ELSEIF(ICON0.EQ.3)THEN C INQUIRE (FILE='adf37',EXIST=B37) IF(.NOT.B37)STOP '*** USER ICON=3 BUT adf37 FILE NOT FOUND...' ICON=3 C OPEN(11,FILE='adf37') READ(11,*) !SKIP HEADER C READ(11,*)IFORM37 !=1 SUPERPOSITION =2 NUMERICAL IF(IFORM37.EQ.1)STOP 'NOT CODED FOR SUPERPOSITION' READ(11,*)IUNITS37 !=1 KELVIN, =2 EV IF(IUNITS37.EQ.1)THEN CONE=CONRYK ELSE CONE=CONRYEV ENDIF C READ(11,*)NENG37 !NO. ENERGIES FOR ADF37 F(E) C READ(11,*)JTEMP37 !NO. TEMPERATURES C IF(NUMTMP.GT.0)JTEMP37=MIN(JTEMP37,NUMTMP) !ALLOW USER RESTRICT IF(JTEMP37.GT.NMTMP)THEN WRITE(0,*)'*** FOR ADF37, INCREASE NMTMP TO',JTEMP37 STOP'*** INCREASE NMTMP FOR ADF37' ENDIF C ALLOCATE (TEMPE37(JTEMP37)) ALLOCATE (E37(NENG37,JTEMP37),F37(NENG37,JTEMP37)) C READ(11,*)IDUM !CUT-OFF THRESHOLD (UNUSED) READ(11,*)IDUM,DUM !HIGH ENERGY BEHAVIOUR (UNUSED) C DO J=1,JTEMP37 TEMPE37(J)=J !SIMPLE INDEX LABEL READ(11,*)(E37(N,J),N=1,NENG37) READ(11,*)(F37(N,J),N=1,NENG37) DO N=1,NENG37 !CONVERT TO RYD IF(E37(N,J).EQ.DZERO)THEN F37(N,J)=DZERO ELSE E37(N,J)=E37(N,J)/CONE F37(N,J)=F37(N,J)*CONE/SQRT(E37(N,J)) !SO FBAR NOW ENDIF ENDDO CALL ADF37(IPRINT,TEMPE37(J),ICON,NENG37,E37(1,J),F37(1,J)) ENDDO CLOSE(11) IRDTMP=-1 !SWITCH-OFF ANY USER READS IF(NUMTMP.LT.0)THEN !READ ANY USER TEMPS FIRST TO POSITION NUMTMP=-NUMTMP READ(5,*)(TDUM,N=1,NUMTMP) ENDIF ELSEIF(ICON0.GT.3)THEN STOP 'ICON.GT.3 IS RESERVED FOR POSSIBLE FUTURE USE' ENDIF C IF(XDRY.GT.DXMAX)STOP X 'MUST SET A SENSIBLE XDRY VALUE FOR THE DRUYVESTEYN DISTRIBUTION' IF(KAPPA.GT.DKMAX)THEN !RE-SET TO MAXWELLIAN KAPPA=-1 ICON=0 ENDIF C IF(ICON.GT.0)THEN IDUP=2 IF(ITYPE.EQ.3)ITYPE=4 ELSE IDUP=1 IF(ITYPE.EQ.4)ITYPE=3 ENDIF IF(IDUP1.LE.0)THEN !SUPPRESS UPSILON IDUP1=1 IDUP=1 ENDIF IDUP0=IDUP !HOLD C IF(ITYPE0.NE.5)THEN IF(ICON.EQ.0)WRITE(6,520) IF(ICON.EQ.1)WRITE(6,521)KAPPA IF(ICON.EQ.2)WRITE(6,522)XDRY IF(ICON.EQ.3)WRITE(6,523) ENDIF C IF(ICON*IREL.NE.0)THEN WRITE(6,776) WRITE(0,776) ENDIF C !NO LONGER USER SET, SO CANNOT HAPPEN... C IF(ITYPE.NE.1.AND.ITYPE.NE.3.AND.ITYPE.NE.4.AND.ITYPE.NE.5) C X STOP 'ILLEGAL ADF04 ITYPE' C IF(IONE.NE.0)IONE=IONE/IABS(IONE) C BBORN=.FALSE. !INFINITE ENERGY IF(IBORN.GT.0)THEN BBORN=.TRUE. ELSEIF(IBORN.EQ.0)THEN IF(ITCC.LT.0)THEN WRITE(*,531) WRITE(6,531) ENDIF ELSE IF(ITCC.GT.0)THEN WRITE(*,532)IBORN WRITE(6,532)IBORN ENDIF ENDIF IF(IBORN.NE.0.AND.IABS(IBORN).NE.2)THEN WRITE(*,533)IBORN WRITE(6,533)IBORN ENDIF C IF(XBTEST.LE.DZERO)XBTEST=5 !APPLY btop ONLY FOR X.GT.XBTEST C IF(NDELBT.LT.0)THEN NDELBT=NMPTS ELSEIF(NDELBT.EQ.0)THEN WRITE(*,534) WRITE(6,534) ENDIF C IF(NUMTMP.EQ.9999)NUMTMP=10 !SET DEFAULT IF(NMTMP.GT.NMBIN)THEN !NOT LIKELY UNLESS USER SWITCHES-OFF BINNED WRITE(6,535)NMBIN,NMTMP STOP 'DIMENSION NMBIN MUST BE AT LEAST NMTMP' ENDIF IF(NUMTMP.LT.0)THEN NUMTMP=-NUMTMP IRDTMP=1 ENDIF IF(NUMTMP.GT.NMTMP)THEN WRITE(6,536)NUMTMP,NMTMP NUMTMP=NMTMP ENDIF TMIN0=TKADAS(1) TMAX0=TKADAS(NTADAS) IF(ICON.EQ.3)THEN !NOW TRANSFER TEMP INFO NUMTMP=JTEMP37 DO I=1,NUMTMP TKADAS(I)=TEMPE37(I)*CONRYK !SINCE TKADAS IS K HERE ENDDO ELSEIF(IRDTMP.NE.0)THEN IRDTMP=1 READ(5,*)(TKADAS(I),I=1,NUMTMP) JTEMP37=NUMTMP NENG37=1 ELSE IF(NUMTMP.GT.13)THEN WRITE(6,538)NUMTMP NUMTMP=13 ENDIF IF(BREAD)THEN !SINCE BYPASS omega READS IRDTMP=-1 DO I=1,NUMTMP TKADAS(I)=TKADAS(I)*FIQ2 ENDDO ENDIF JTEMP37=NUMTMP NENG37=1 ENDIF IF(ICON0.NE.3)ALLOCATE (E37(NENG37,JTEMP37),F37(NENG37,JTEMP37)) IF(NUMTMP.GT.20)WRITE(6,537) C C READ (CONFIGURATION) TERM SPECIFICATIONS AND J-VALUES C FROM adasex/j.in *OR* FROM adf04_om C IF(BREAD)THEN NSTIN=NSTMX !MAX ALLOWED BY DIMEN SINCE WE DON'T KNOW ELSE NSTIN=NSTATES !USER HAS FLAGGED HOW MANY TO USE IF(NSTIN.GE.NSTMX)THEN !AS MAY ACCESS NAST+1 WRITE(6,551)NSTIN STOP 'INCREASE THE PARAMETER NSTMX' ENDIF ENDIF C C----------------------------------------------------------------------- C C CONF(I) -- CONFIG/TERM SPECIFICATION FOR THE ITH STATE - UP TO C 15 CHARACTERS. NOTE THE TOTAL VALUES OF 2S+1 AND C AND L ARE DETERMINED FROM THE LAST FOUR CHARACTERS C WITHIN TERM, AND THESE FOUR CHARACTERS MUST START C WITH (, AND END WITH ), AND MUST BE IN THE FORM: C (1S),OR (3P), OR (5D), ETC. C FJT(I) -- J-VALUE FOR THE ITH STATE C C CASE ITCC.EQ.0 (JK-COUPLING) READ ALSO ENERGIES, AND INDEX VALUES: C C ENERG(I) -- ENERGY OF THE ITH LEVEL IN CM-1. C INDEX(I) -- THIS IS THE LEVEL INDEX USED TO ASSOCIATE THE INPUT C LEVELS WITH THE LEVELS READ FROM THE omega FILE C GENERATED USING STGICF. THIS IS ESPECIALLY USEFUL C IN A PURE JK COUPLING RUN WHERE THE LEVELS ARE C ORDERED FIRST BY TERM AND THEN BY INCREASING J C VALUE. IF THE ORDER OF THE INPUT LEVELS AND THOSE C IN THE omega FILE ARE THE SAME THEN THE INDICES CAN C BE LEFT BLANK AND THE CODE WILL AUTOMATICALLY SET C THEM EQUAL TO THE LEVEL NUMBER READ IN. C C CASE ITCC.NE.0 (LS/ICFT/BPRM/DARC) C C ENERGY AND INDEX ARE IGNORED. IF THE CONFIGURATION/TERM/LEVEL ORDER C IN adasex/j.in DOES NOT MATCH THAT (IMPLIED) IN omega THEN C CUT-AND-PASTE THE CORRECT ORDER INTO adasex/j.in. C C----------------------------------------------------------------------- C C CHECK FOR OLD OR NEW STYLE TERM/CONFIG LIST C READ(IREAD,501)TRM,CONF(1),CONF(2) BACKSPACE(IREAD) C CARD=' -1' !DEFAULT TERMINATOR BSHORT=.TRUE. DO I=LEN(TRM),1,-1 IF(TRM(I:I).EQ.')')GO TO 5 ENDDO C C ASSUME NEW STYLE, AND NOT BLANK OLD-STYLE... READ(IREAD,502)TRM,CONF(1),CONF(2) BACKSPACE(IREAD) TRM=CONF(1) IF(TRM(6:6).EQ.'0')THEN !SHORT CA F542='(I5,1X,A15,4X,I1,1X,I1,1X,F10.1,1X,F15.1)' F580="(I5,1X,A15,2X,' (',I1,')',I1,'(',F10.1,')',F15.1)" c write(0,*)'short ca' GO TO 3 ENDIF TRM=CONF(2) IF(TRM(7:7).EQ.'0')THEN !LONG CA F543='(I5,1X,A31,4X,I1,1X,I1,1X,F10.1,1X,F15.1)' F581="(I5,1X,A31,2X,' (',I1,')',I1,'(',F10.1,')',F15.1)" c write(0,*)'long ca' GO TO 4 ENDIF DO I=LEN(TRM),1,-1 IF(TRM(I:I).EQ.')')THEN F543='(I5,1X,A31,4X,I1,1X,I1,1X,F4.1,1X,F21.4)' F581="(I5,1X,A31,2X,' (',I1,')',I1,'(',F4.1,')',F21.4)" c write(0,*)'long ls/ic' GO TO 4 ENDIF ENDDO F542='(I5,1X,A15,4X,I1,1X,I1,1X,F4.1,1X,F21.4)' F580="(I5,1X,A15,2X,' (',I1,')',I1,'(',F4.1,')',F21.4)" C SHORT 3 IF(BREAD)THEN DO L=1,NSTIN READ(IREAD,F542)IN BACKSPACE(IREAD) IF(IN.GT.0)THEN READ(IREAD,F542)INDEX(L),CONF(L),ISAT(L),LAT(L),FJT(L) X ,ENERG(L) C WRITE(6,F580)INDEX(L),CONF(L),ISAT(L),LAT(L),FJT(L) C X ,ENERG(L) ELSEIF(IN.LT.0)THEN READ(IREAD,590)CARD GO TO 15 ELSE STOP 'NON-STANDARD ADF04 FILE' ENDIF ENDDO WRITE(6,551)NSTIN*2 STOP 'INCREASE THE PARAMETER NSTMX' ELSE !NOT CA DO L=1,NSTIN READ(IREAD,542)INDEX(L),CONF(L),ISAT(L),LAT(L),FJT(L),ENERG(L) C WRITE(6,580)INDEX(L),CONF(L),ISAT(L),LAT(L),FJT(L),ENERG(L) ENDDO GO TO 15 ENDIF C LONG 4 BSHORT=.FALSE. IF(BREAD)THEN DO L=1,NSTIN READ(IREAD,F543)IN BACKSPACE(IREAD) IF(IN.GT.0)THEN READ(IREAD,F543)INDEX(L),XCONF(L),ISAT(L),LAT(L),FJT(L) X ,ENERG(L) C WRITE(6,F581)INDEX(L),XCONF(L),ISAT(L),LAT(L),FJT(L) C X ,ENERG(L) ELSEIF(IN.LT.0)THEN READ(IREAD,590)CARD GO TO 15 ELSE STOP 'NON-STANDARD ADF04 FILE' ENDIF ENDDO WRITE(6,551)NSTIN*2 STOP 'INCREASE THE PARAMETER NSTMX' ELSE !NOT CA DO L=1,NSTIN READ(IREAD,543)INDEX(L),XCONF(L),ISAT(L),LAT(L),FJT(L) X ,ENERG(L) C WRITE(6,581)INDEX(L),XCONF(L),ISAT(L),LAT(L),FJT(L) C X ,ENERG(L) ENDDO GO TO 15 ENDIF C C OLD STYLE (NOT CA) 5 IF(BREAD)STOP 'CANNOT PROCESS OLD STYLE INPUT HERE' LNSTR=I !SINCE RIGHT JUSTIFIED DO 10 L=1,NSTATES READ(IREAD,500)CONF(L),FJT(L),ENERG(L),INDEX(L) C WRITE(6,580)CONF(L),FJT(L),ENERG(L),INDEX(L) IF(ITCC.NE.0.OR.INDEX(L).EQ.0)INDEX(L)=L TRM=CONF(L) ISBCD=TRM(LNSTR-2:LNSTR-2) LTBCD=TRM(LNSTR-1:LNSTR-1) DO I=1,7 IF(SBCD(I).EQ.ISBCD)THEN ISAT(L)=I GO TO 8 ENDIF ENDDO STOP 'INPUT OF TERMS INVALID' 8 DO I=1,10 IF(LBCD(I).EQ.LTBCD)THEN LAT(L)=I-1 GO TO 10 ENDIF ENDDO STOP 'INPUT OF TERMS INVALID' 10 ENDDO C C----------------------------------------------------------------------- C 15 CONTINUE C C POSSIBLE EXIT TO READ adf04_om RATHER THAN omega C IF(BREAD)THEN C NSTATES=L-1 !WE NOW HAVE THE ACTUAL NUMBER OF LEVELS IF(NSTATES.EQ.NSTIN)THEN !LIKELY NEED TO INCREASE MSTMX WRITE(6,551)2*NSTIN STOP 'INCREASE THE PARAMETER NSTMX' ENDIF C DO L=1,NSTATES ITOT(L)=0 !BIN LENGTH, TO BE SET IWJ=NINT(2*FJT(L)+1) LX=INDEX(L) JXTWO(LX)=IWJ ENDDO C !SANITY CHECK READ(IREAD,*)DUM,ITYPIN,EDUM C IF(NINT(DUM).NE.IQ1)STOP 'IQ1 MIS-MATCH' !SHOULD NOT HAPPEN C IF(ITYPIN.NE.1.AND.ITYPIN.NE.5)THEN !CHECK OMEGA NOT UPSILON WRITE(6,541)ITYPIN STOP 'adf04_om IS NOT ITYPE=1 OR 5!' ENDIF C C IMPORTANT TO KNOW IF R-MX OR DW/PWB (NO PROBLEM DISTINGUISHING LATTER) C IF(.NOT.BPWB)BPWB=ITYPIN.EQ.1 !PWB SHOULD BE ITYPIN=1 C IF(.NOT.BDW)BDW=EDUM.EQ.0 !ALLOW OVERRIDE E.G. NEUTRALS C BRMX=.NOT.BDW.AND..NOT.BPWB .and.itypin.eq.5 C IF(NZED.EQ.NELC.AND.BRMX)THEN !CANNOT DETECT NEUTRAL AS-DW, SO WRITE(6,*)'*** ATTENTION: ASSUMING R-MATRIX adf04_om INPUT!' WRITE(6,*)'*** SET BDW=.TRUE. FOR NEUTRAL AS-DW INPUT' WRITE(*,*)'*** ATTENTION: ASSUMING R-MATRIX adf04_om INPUT!' WRITE(*,*)'*** SET BDW=.TRUE. FOR NEUTRAL AS-DW INPUT' ENDIF c btop=btop.and.bdw C C SET OUTPUT ITYPE C (ITYPE0 NOT USED BEFORE RE-SET, NOR ITYPE.EQ.4) C IF(BRMX)THEN !RE-SET ITYPE=3 !TO SKIP BIN SET-UP (WILL READ) ELSE !USES "BINS" TO INTERPOLATE ITYPE=5 BEQUAL=.TRUE. !SINCE (AS) USES COMMON FINAL MESH IF(BPWB.AND.IRMPS.EQ.0)IRMPS=-1 ENDIF C C GET (MAX) NUMBER OF COLLISION ENERGIES C KMAX=0 17 KMAX=KMAX+20 READ(IREAD,'(I5)')NTEST IF(NTEST.EQ.0)GO TO 17 IF(KMAX.GE.NMBIN)THEN !+1 FOR INF E PT WRITE(6,585)KMAX STOP 'INCREASE NMBIN' ENDIF BACKSPACE(IREAD) READ(IREAD,'(I4)')ITEST BACKSPACE(IREAD) C C TEST FOR OLD adf04 WITH I4.GT.999 LEVELS - NEW SWITCHES TO I5 FOR GAP C NEED TO WORK ON FORMAT I/O TESTS... C IF(ITEST.NE.0.AND.NSTATES.GE.1000)THEN !NEED TO RE-SYNC TEST I/O WRITE(6,*)'***INTERNAL ERROR, EXPECTING I5 FOR .GE.1000 LEVELS' STOP 'INTERNAL ERROR' ENDIF IF(ITEST.EQ.0.AND.NSTATES.LT.1000)THEN !NEED TO RE-SYNC TEST I/O WRITE(6,*)'***INTERNAL ERROR, EXPECTING I4 FOR .LT.1000 LEVELS' STOP 'INTERNAL ERROR' ENDIF C C TRY AND CATCH ANY IN/ELASTIC INCONSISTANCY C READ(IREAD,*)J,I BACKSPACE(IREAD) C IF(IONE.EQ.1.AND.J.EQ.I)THEN WRITE(6,*)'*** ELASTIC TRANSITIONS PRESENT, SET IONE=0' STOP '*** ELASTIC TRANSITIONS PRESENT, BUT NOT EXPECTED!' ENDIF IF(IONE.EQ.0.AND.J.NE.I)THEN WRITE(6,*)'*** ELASTIC TRANSITIONS NOT PRESENT, SET IONE=1' STOP '*** ELASTIC TRANSITIONS NOT PRESENT, BUT EXPECTED!' ENDIF C NTLINE=KMAX/20 DO N=1,NTLINE BACKSPACE(IREAD) ENDDO C C DETERMINE MAX NUMBER OF TRANSITIONS REQUIRED. C (THERE MAY BE LESS ON FILE.) C LIMN=MAX(1,NLOWMN) LIMX=MIN(NSTATES-IONE,NLOWMX) LFMX=MIN(NSTATES,NUPMX) C NTRAN=0 DO LI=LIMN,LIMX LFMN=MAX(LI+IONE,NUPMN) DO LF=LFMN,LFMX NTRAN=NTRAN+1 ENDDO ENDDO C NTRANS=NTRAN C ENDIF C C----------------------------------------------------------------------- C C READ INFORMATION DESCRIBING DATA C C NAME -- NAME OF PERSON GENERATING THE EXCITATION DATA - UP TO C 30 CHARACTERS C DATE -- DATE - UP TO 30 CHARACTERS C COMMENT -- DATA COMMENTARY - UP TO 72 CHARACTERS PER LINE C END COMMENTS WITH A PERIOD IN THE FIRST COLUMN OF C THE LAST LINE. LIMITED TO NLNMX LINES, INCLUDING THE C PERIOD ON THE LAST LINE. C C----------------------------------------------------------------------- C IF(IREAD.EQ.5)THEN !SO BYPASS IF BREAD READ(5,510)NAME READ(5,510)DATE I=1 20 COMMENT(I)(1:6)='C ' READ(5,530)COMMENT(I)(7:80) IF(COMMENT(I)(7:7).EQ.'.')GO TO 30 I=I+1 IF(I.GT.NLNMX)THEN WRITE(6,540)NLNMX STOP '*** INCREASE THE PARAMETER NLNMX (OR SHORTEN COMMENTS)' ENDIF GO TO 20 30 NMLNS=I-1 ENDIF C C WE HAVE ALL USER INPUT C C----------------------------------------------------------------------- C C SET-UP ANY BIN ENERGIES C EFITL=ENERG(NSTATES)/(CONRCM*TAZ) C NTOT=0 IF(ITYPE.EQ.5)THEN IF(EMAX.LT.DZERO)EMAX=-EMAX*TAZ !UNSCALE IF(EMIN.LT.DZERO)EMIN=-EMIN*TAZ !UNSCALE IF(DELIN.LT.DZERO)DELIN=-DELIN*TAZ !UNSCALE IF(MAXE.LT.DZERO)MAXE=-MAXE*TAZ !UNSCALE IF(EMAX.EQ.DZERO)THEN EMAX=EFITL*TAZ !MAX SPECTROSCOPIC IF(IONE.EQ.1)EMAX=EMAX-ENERG(2)/CONRCM !FINAL ENERGY ENDIF IF(EMIN.EQ.DZERO)THEN EMIN=FIQ2*TMIN0/(10*CONRYK) !FIQ2 AS ADAS TEMP if(icon.eq.1.and.kappa.lt.d2pt5) x emin=emin/(done-log(kappa-d1pt5)) if(icon.eq.2.and.xdry.lt.done)emin=emin*xdry**3 ENDIF C T0=LOG10(EMIN) IF(NBIN.LT.10)THEN COLD NBIN=101 ERES=LOG10(DELOG0) !DO NOT "Z-SCALE"! NBIN1=(LOG10(EMAX)-T0)/ERES NBIN=NBIN1+1 ELSE NBIN1=NBIN-1 !NO. OF BINS ERES=(LOG10(EMAX)-T0)/NBIN1 ENDIF IF(NBIN.GT.NMBIN)THEN WRITE(6,585)NBIN STOP 'INCREASE THE PARAMETER NMBIN' ENDIF C N=0 EBIN(0)=DZERO REWIND(98) DO N0=0,NBIN1 T=N0 T=T0+T*ERES T=10**T WRITE(98,730)T !ROUND TO 3 S.F. BACKSPACE(98) READ(98,750)T,IT T=T*DTEN**IT !UNSCALED TO MATCH EJ IF(ABS(T-EBIN(N)).GT.DZERO)THEN N=N+1 EBIN(N)=T ENDIF ENDDO IF(NBIN.GT.N)THEN write(*,*)nbin-n,' degenerate bin energy cases dropped' NBIN=N NBIN1=NBIN-1 ENDIF ERES0=10**ERES ERES1=EBIN(NBIN)-EBIN(NBIN1) IF(NLIN.LE.0)THEN NTOT=NMBIN NLIN=NTOT-NBIN NLIN=-NLIN ELSE NTOT=NBIN+NLIN IF(NTOT.GT.NMBIN)THEN WRITE(6,585)NTOT STOP 'INCREASE THE PARAMETER NMBIN' ENDIF ENDIF IF(DELIN.EQ.DZERO)THEN IF(MAXE.GT.DZERO)THEN DELIN=(MAXE-EMAX)/IABS(NLIN) IF(NLIN.LT.0.AND.DELIN.GT.D0PT1)THEN NLIN=(MAXE-EMAX)/D0PT1 NLIN=NLIN+1 NTOT=NBIN+NLIN IF(NTOT.GT.NMBIN)THEN WRITE(6,585)NTOT STOP 'INCREASE THE PARAMETER NMBIN' ENDIF ENDIF ELSE DELIN=DELIN0*TAZ ENDIF DELIN=MAX(DELIN,ERES1) ENDIF IF(DELIN.GT.ERES1)THEN IF(DELIN/TAZ.GT.D0PT1)WRITE(*,*)'*** DELIN TOO LARGE? ',DELIN ELSE IF(DELIN/TAZ.GT.D0PT2)WRITE(*,*)'*** DELIN=ERES1 TOO LARGE? ' X ,DELIN ENDIF NLIN=IABS(NLIN) REWIND(98) DO N=1,NLIN T=EBIN(NBIN)+DELIN*N WRITE(98,730)T BACKSPACE(98) READ(98,750)T,IT EBIN(NBIN+N)=T*DTEN**IT !UNSCALED TO MATCH EJ ENDDO IF(BREAD)THEN WRITE(6,586)NBIN,EMIN,ERES0,EMAX,ERES1,EBIN(NTOT),DELIN ELSE WRITE(6,587)NBIN,EMIN,ERES0,EMAX,ERES1,EBIN(NTOT),DELIN DO N=1,NSTATES !INITIAL MAX NO BINS FOR ALL STATES ITOT(N)=NTOT ENDDO ENDIF ENDIF C IF(BREAD)THEN C KXX=MAX(KMAX,NTOT,NUMTMP) KXX=KXX+1 !FOR INF E C ALLOCATE (AR(NSTATES,NSTATES) X ,UPS(NUMTMP,IDUP,NTRANS),OMEGIN(KXX,NTRANS)) C NBIN=NTOT !HOLD NMLNS=0 GO TO 170 !BYPASS OMEGA >>>>>>>>>>>>>>>>>>>> ENDIF C C*********************************************************************** C C START READ/PROCESSING OF R-MATRIX omega FILE C C*********************************************************************** C C READ BASIC DATA C C IF(BFORM)READ(7,*)NZED,NELC C IF(.NOT.BFORM)READ(7)NZED,NELC C IF(BFORM)READ(7,*)NAST,MXE,NOMWRT IF(.NOT.BFORM)READ(7)NAST,MXE,NOMWRT C IF(NOMWRT.NE.0)THEN NOMT=IABS(NOMWRT) ELSE NOMT=(NAST*(NAST+1-2*IONE))/2 ENDIF C IF(NAST.GT.NSTMX)THEN WRITE(6,550)NAST STOP 'INCREASE THE PARAMETER NSTMX' ENDIF IF(NSTATES.GT.NAST)THEN WRITE(6,555)NAST,NSTATES STOP 'NO. OF CC STATES EXCEEDS THOSE IN adasex/j.in...' ENDIF IF(MXE.GT.NMPTS)THEN WRITE(6,560)MXE STOP 'INCREASE THE PARAMETER NMPTS' ENDIF C IF(BFORM)READ(7,*)(IZERO(I),JXTWO(I),I=1,NAST) IF(.NOT.BFORM)READ(7)(IZERO(I),JXTWO(I),I=1,NAST) IF(BFORM)READ(7,*)(ENAT(I),I=1,NAST) IF(.NOT.BFORM)READ(7)(ENAT(I),I=1,NAST) C IF(ITCC.GE.0.AND.IZERO(1).NE.0) X STOP ' IC/JK RUN FLAGGED BUT omega FILE CONTAINS LS DATA' IF(ITCC.LT.0.AND.IZERO(1).EQ.0) X STOP ' LS RUN FLAGGED BUT omega FILE CONTAINS IC/JK DATA' C IF(IRMPS.EQ.0)THEN EFITL=ENAT(NAST) ELSE EFITL=ENAT(NSTATES) ENDIF c if(eshft.gt.dzero)eshft=eshft/taz eshft=abs(eshft) C C CARRY-OUT BASIC COMPATIBILTY CHECKS ON adasex/j.in AND omega FILES C DO L=1,NSTATES C IWJ=NINT(2*FJT(L)+1) LX=INDEX(L) JXTWO(LX)=IABS(JXTWO(LX)) IF(ITCC.LT.0)JXTWO(LX)=(2*JXTWO(LX)+1)*IABS(IZERO(LX)) !LS C IF(IWJ.NE.JXTWO(LX))THEN WRITE(6,575)L,IWJ,LX,IZERO(LX),JXTWO(LX) STOP 'MIS-MATCH OF TARGET STATES: adasex/j.in AND omega' ENDIF C IF(ITCC.NE.0)THEN !IC/LS EE=ENERG(L)/(TAZ*CONRCM) IF(ABS(ENAT(LX)-EE).GT.TOLE)THEN !ENERG CONTAINS OBSERVED? WRITE(6,576)EE,ENAT(LX) ELSE !IGNORE ENERG(L)=ENAT(LX)*TAZ*CONRCM ENDIF ENDIF C ENDDO C C DETERMINE NUMBER OF TRANSITONS TO BE READ AT ANY COLLISION ENERGY C IF(NOMWRT.GE.0)THEN DO I=1,MXE NOMRD(I)=NOMT ENDDO ELSE IF(BFORM)READ(7,*)E0 !CHECK FIRST E IF(.NOT.BFORM)READ(7)E0 !CHECK FIRST E BACKSPACE(7) I=1 IF(E0.GE.ENAT(NAST))GO TO 41 !ALL OPEN - EARLY EXIT IF(BFORM)THEN I0=1 NO=NOPEN(E0,ENAT,NAST,IONE,I0) NO=MIN(NO,NOMT) READ(7,*)E0,(DUM,N=1,NO) READ(7,*)E !=E0+EINCR EINCR=E-E0 E0=E0-EINCR EINCH=EINCR/2 NREC=2+(NO-1)/6 DO N=1,NREC BACKSPACE(7) ENDDO ENDIF I0=1 DO I=1,MXE IF(BFORM)THEN E=E0+I*EINCR WRITE(99,F600)E BACKSPACE(99) READ(99,F600)E BACKSPACE(99) ELSE READ(7)E ENDIF IF(E.GE.ENAT(NAST))GO TO 41 NOMRD(I)=NOPEN(E,ENAT,NAST,IONE,I0) NOMRD(I)=MIN(NOMRD(I),NOMT) ENDDO I=MXE+1 41 I0=I DO I=I0,MXE NOMRD(I)=NOMT ENDDO IF(.NOT.BFORM)THEN REWIND(7) NREC=4 DO N=1,NREC READ(7) ENDDO ENDIF ENDIF C C SET-UP BATCHES OF ENERGIES TO BE READ SO AS NOT TO EXCEED C THE NMPIN DIMENSION OF OMEGIN. C (MXE.LE.NMPTS STILL REQUIRED.) C IF(MXEIN.LE.0)THEN MXEIN=NMPIN ELSE IF(MXEIN.GT.NMPIN)THEN WRITE(6,*)'USER INPUT ERROR: REDUCE MXEIN TO ',NMPIN STOP 'REDUCE MXEIN' ENDIF ENDIF IF(MXE.LE.MXEIN)THEN NBATCH=1 MBATCH=MXE ELSE III=1 IF(BFORM)III=0 MBATCH=MXEIN-III IF(BBORN.AND.IRMPS.GE.0)THEN !FOR FORBIDDEN POWER LAW 55 NBATCH=1+(MXE-1-III)/MBATCH MNB=MBATCH*(NBATCH-1) MHIGH=MXE-MNB IF(MHIGH.GT.0.AND.MHIGH.LT.500)THEN MBATCH=MBATCH-MBATCH/5 IF(NBATCH.GT.100.OR.MBATCH.LT.500)THEN WRITE(6,595)MXEIN STOP ' BATCH ENERGY SET-UP OUT OF CONTROL (MXEIN)' ENDIF GO TO 55 ENDIF ENDIF ENDIF C ALLOCATE (OMEGIN(MBATCH,NOMT)) C C SET-UP TEMPERATURES C IF(IRDTMP.EQ.0)THEN !CONVERT INTERNAL Z-SCALED TO ACTUAL TEMPS DO IT=1,NUMTMP TKADAS(IT)=TKADAS(IT)*FIQ2 ENDDO ENDIF C C SET-UP TRANSITIONS AND INDEXING C C MAX VALUE OF NTRAN USED IN NON-ALLOCATABLE DIMENSION TEST: C NTRAN=(NSTATES*(NSTATES+1-2*IONE))/2 C NTRAN=0 C IF(.NOT.BREV)THEN !STANDARD ORDER C LIMN=MAX(1,NLOWMN) LIMX=MIN(NSTATES-IONE,NLOWMX) LFMX=MIN(NSTATES,NUPMX) C DO LI=LIMN,LIMX C LFMN=MAX(LI+IONE,NUPMN) C DO LF=LFMN,LFMX C NTRAN=NTRAN+1 LIN(NTRAN)=LI LFN(NTRAN)=LF C ENDDO C ENDDO C ELSE !FLAG REVERSE (AS) ORDER C LIMN=MAX(1,NLOWMN) LFMN=MAX(1+IONE,NUPMN) LFMX=MIN(NSTATES,NUPMX) C DO LF=LFMN,LFMX C LIMX=MIN(LF-IONE,NLOWMX) C DO LI=LIMN,LIMX C NTRAN=NTRAN+1 LIN(NTRAN)=LI LFN(NTRAN)=LF C ENDDO C ENDDO C ENDIF C NTRANS=NTRAN C C INITIALIZE UPSILON C KMAX=NUMTMP IF(ITYPE.EQ.5)THEN KMAX=NTOT IDUP=1 !BIN FIRST ELSEIF(ICON.GT.0)THEN !ITYPE=3->4 ITYPE=4 ENDIF KXX=MAX(KMAX,NUMTMP) C ALLOCATE (AR(NSTATES,NSTATES),UPS(KXX,IDUP,NTRANS)) C IF(BOMTHR)THEN !FILL-IN FROM THRESHOLD UPSMIN=DZERO ELSE !SKIP AS STARTING AT "HIGH-E" UPSMIN=D1M30 ENDIF C DO N=1,NTRANS DO I=1,IDUP DO K=1,KMAX UPS(K,I,N)=UPSMIN ENDDO ENDDO ENDDO C C RE-ENTRY POINT FOR NEW ENERGY BATCH C IE1=0 C 60 CONTINUE ! <<<<<<<<<<<<<<<<<<<< C IE0=IE1+1 IE1=IE1+MBATCH IE1=MIN(IE1,MXE) BLAST=IE1.EQ.MXE C C READ ELECTRON-IMPACT EXCITATION COLLISION STRENGTHS C IOM=0 IF(BFORM)THEN DO I=IE0,IE1 IOM=IOM+1 READ(7,F600)EMSH(I),(OMEGIN(IOM,N),N=1,NOMRD(I)) IF(NOMWRT.LT.0)THEN ET=E0+I*EINCR IF(ABS(EMSH(I)-ET).GT.EINCH)THEN IF(NOMRD(I).LT.NOMT)THEN WRITE(6,*)'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS' WRITE(6,*)I,NOMRD(I),ET,EMSH(I) STOP'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS' ENDIF ENDIF ENDIF ENDDO ELSE DO I=IE0,IE1 IOM=IOM+1 READ(7)EMSH(I),(OMEGA(N),N=1,NOMRD(I)) !OMEGA*8 DO N=1,NOMRD(I) OMEGIN(IOM,N)=OMEGA(N) !OMEGIN*4 ENDDO ENDDO ENDIF IF(ITCC.EQ.0.AND.NOMWRT.LT.0)THEN !ZERO IOM=0 DO I=IE0,IE1 IOM=IOM+1 DO N=NOMRD(I)+1,NOMT OMEGIN(IOM,N)=0 ENDDO ENDDO ENDIF C IFIN=IE1 !FINAL (GLOBAL) ENERGY IOMINF=IOM !FINAL (BATCH) ENERGY, OMEGIN SHOULD HOLD INF E LIMIT C C TEST ENERGY MESH FOR SEQUENTIAL POINTS C DO I=IE0+1,IE1 H=EMSH(I)-EMSH(I-1) IF(H.LT.-TOLH0)THEN WRITE(6,*)'*** MESH ERROR:',I-1,EMSH(I-1),I,EMSH(I) STOP '*** MESH ERROR, NEED INCREASING ENERGIES!' ENDIF ENDDO C C CHECK TO SEE IF WE START AT (FIRST) THRESHOLD C IF(IE0.EQ.1)THEN H=EMSH(2)-EMSH(1) T=EMSH(1)-ENAT(1+IONE) IF(T.GE.5*H.AND.BOMTHR)THEN IF(NZED.EQ.NELC.AND.IONE.EQ.0)THEN WRITE(6,604)T,H STOP X'*** OMEGA FILE DOES NOT START AT THRESHOLD?, CHECK IONE?' ELSE WRITE(6,605)T,H STOP X'*** OMEGA FILE DOES NOT START AT THRESHOLD, SET BOMTHR=.FALSE.' ENDIF ELSEIF(T.LT.5*H.AND..NOT.BOMTHR)THEN WRITE(6,606)T,H STOP '*** OMEGA FILE STARTS AT THRESHOLD, SET BOMTHR=.TRUE.' ENDIF ENDIF C C THAT WAS THE FINAL BATCH C IF(BLAST)THEN !HAVE ALL ENERGIES NOW C IF(EMSH(MXE).GT.EINF-TOLH0)THEN BINF=.TRUE. EINF=EMSH(MXE) !JUST IN CASE TINF=EINF*TAZ IFIN=MXE-1 !LAST FINITE ENERGY ELSEIF(EMSH(MXE).GE.3*ENAT(NAST)/2)THEN WRITE(6,610) STOP 'ERROR: INFINITE ENERGY INFO MISSING' ELSE WRITE(6,611) ENDIF C EFIN=EMSH(IFIN) C ENDIF C C*********************************************************************** C C END READ/PROCESSING OF R-MATRIX omega FILE C C*********************************************************************** C C C RE-ENTRY POINT FOR CONVOLUTION OF DW/PWB OMEGAS. C 65 CONTINUE ! <<<<<<<<<<<<<<<<<<<< C IF(BLAST)THEN C C REDUCE ANY BIN ENERGY RANGE C IF(ITYPE.EQ.5.OR.ITYPE0.EQ.-3)THEN NLAST=MIN(NSTATES,NUPMX) IF(BEQUAL)THEN !RESTRICT TO SHORTEST FINAL MESH KLAST=1 ELSE !ALLOW FULL (UNEQUAL) RANGE FOR ALL STATES KLAST=NLAST NLAST=1+IONE !CAN RESTRICT GROUND IF NO ELEASTIC ENDIF IF(BRMX)THEN NLAST0=NLAST ELSE NLAST0=1 ENDIF EFN=EFIN*TAZ-ENERG(NLAST0)/CONRCM IF(EFN.GT.EBIN(NTOT))THEN WRITE(6,612)EBIN(NTOT),EFN STOP '*** LINEAR ENERGY MESH TOO SHORT' ELSE DO K=1,KLAST DO N=NTOT,1,-1 IF(EFN.GT.EBIN(N))GO TO 67 ENDDO 67 NTOT=N !STOP A TAD SHORT C N=N+1 !EXTRAP TO FILL BIN ITOT(K)=N EFN=EFIN*TAZ-ENERG(K+1)/CONRCM ENDDO NTOT=ITOT(1) DO K=KLAST+1,NLAST ITOT(K)=NTOT ENDDO WRITE(6,613)NTOT,EBIN(NTOT) !,ifin,efin IF(ITYPE0.EQ.-5.AND.NTOT.GE.NMPIN)THEN WRITE(6,614)NTOT STOP '*** INCREASE NMPIN' ENDIF ENDIF KMAX=NTOT IF(ITYPE.EQ.5)THEN !GO DETERMINE BINNED COLLISIONS STRENGTHS GO TO 90 ELSE !INTERPOLATE & CONVOLUTE DW/PWB GO TO 70 ENDIF ENDIF C ENDIF C C RE-ENTRY POINT FOR CONVOLUTION OF BINNED/READ RMX OMEGAS. C 70 CONTINUE ! <<<<<<<<<<<<<<<<<<<< C IF(BLAST)THEN C C ADJUST SCALINGS C TOLH=TOLH0*TAZ IF(EASYM.GT.DZERO)EASYM=EASYM/TAZ !CONVERT TO Z-SCALED IF(EFITF.GT.DZERO)EFITF=EFITF/TAZ !CONVERT TO Z-SCALED EASYM=ABS(EASYM) !ELSE REMOVE SIGN EFITF=ABS(EFITF) !ELSE REMOVE SIGN C C SET GLOBAL ENERGY RANGE TO BE USED TO DETERMINE THE POWER LAW FOR C HIGH ENERGY FORBIDDEN TRANSITIONS C IF(EFITF.EQ.DZERO)EFITF=EFIN*2/3 EFITF=MAX(EFITF,EFITL) IF(EASYM.LT.EINF)WRITE(6,615)EFITF*TAZ,EASYM*TAZ C C CHECK TO SEE IF USER IS REDUCING THE MINIMUM TEMPERATURE BELOW C THE ADAS MINIMUM AND WARN IF NECESSARY C IF(IRDTMP.NE.0)THEN TMIN0=TMIN0*FIQ2 TMIN=TKADAS(1) DO IT=2,NUMTMP IF(TKADAS(IT).LT.TMIN)TMIN=TKADAS(IT) ENDDO IF(TMIN.LT.TMIN0/2)THEN WRITE(6,547)TMIN,TMIN0 IF(NDELBT.LE.10)THEN WRITE(*,548)NDELBT WRITE(6,548)NDELBT ENDIF ENDIF ENDIF C C CHECK TO SEE THAT THE MAXIMUM TEMPERATURE DOES NOT EXCEED C HALF THE MAXIMUM SCATTERING ENERGY. C IF IT DOES, THEN IF MXTMP.EQ.0, REDUCE THE MAXIMUM TEMPERATURE, C ELSE JUST WARN. C EMXTMP=EFIN*TAZ*CONRYK/2 C DO IT=1,NUMTMP IF(EMXTMP.LT.TKADAS(IT))GO TO 75 ENDDO GO TO 80 75 IF(MXTMP.EQ.0)THEN !REDUCE WRITE(6,616)TKADAS(NUMTMP),TKADAS(IT-1) NUMTMP=IT-1 ELSE !JUST WARN IF(EMXTMP.LT.TKADAS(NUMTMP)/2) X WRITE(6,618)TKADAS(NUMTMP),EMXTMP*2 ENDIF 80 KMAX=NUMTMP C ENDIF C C C*********************************************************************** C C *** BEGIN LOOP OVER ALL TRANSITIONS *** C TO BIN OR CONVOLUTE COLLISION STRENGTHS C C*********************************************************************** C 90 CONTINUE C C IF EBIN IS ALREADY EJ, JUST TRANSFER ONCE. C DW/PWB MUST START AT THRESHOLD (PWB MAY NEED TWEAKING) C THIS WILL BE AT POSITION 1, UNLIKE R-MX BINNED WHICH IS AT POSITION 0. C IF(ITYPE0.EQ.-3)THEN I00=1 C IE0=1 !ALREADY SET IF(.NOT.BRMX)THEN T=EMSH(IE0) IFIN0=IFIN IF(ITYPIN.EQ.1)THEN !TRANSFER THRESHOLD PARAM IF(BPWB.AND.NZED.NE.NELC)CALL REMAPX(EMSH,OMEGA,IFIN) DO I=IE0,IFIN XB(I)=EMSH(I)-DONE ENDDO T=ABS(XB(1)) ENDIF IF(T.GT.ETHRSH0)THEN IF(NZED.NE.NELC)THEN WRITE(6,620)EMSH(IE0) STOP '***ERROR: IONIC DW/PWB MUST START AT THRESHOLD!' ENDIF C AS OMITS THRESHOLD FOR NEUTRALS, SO I00=0 EMSH(0)=0 XMANT(0)=0 XB(0)=0 ENDIF IOMINFO=IOMINF IOMINF=ITOT(1)+1 ENDIF C NTOT=ITOT(1) DO I=IE0,NTOT EJ(I)=EBIN(I) ENDDO C EJ(0)=0 ENDIF C do 140 NTRAN=1,NTRANS C LI=LIN(NTRAN) LF=LFN(NTRAN) C LOI=INDEX(LI) LOF=INDEX(LF) C IF(LOI.LT.LOF)THEN ILI=LOI ILF=LOF ELSE !GET ORIGINAL ORDER ILI=LOF ILF=LOI ENDIF C NTR=0 IF(NOMWRT.LT.0)NTR=((ILF-1)*(ILF-2*IONE))/2+ILI IF(NOMWRT.GT.0)NTR=ILF+NAST*(ILI-1)-(ILI*(ILI-1+2*IONE))/2 IF(NTR.GT.NOMT)GO TO 140 C EF=ENERG(LF)/(CONRCM*TAZ) EI=ENERG(LI)/(CONRCM*TAZ) ETHRSH=EF-EI !Z-SCALED ETHRSH=MAX(ETHRSH,ETHRSH0) !FOR DEGENERATE LEVELS ETHRSH=ETHRSH*TAZ C IEASYM=0 IEFT=0 C C IF WE ALREADY HAVE EJ, THEN JUST TRANSFER OMEGIN TO OMEGA C (INTERPOLATE FIRST IF DW/PWB) C IF(ITYPE0.EQ.-3)THEN C NTR=NTRAN !AS JUST MIRRORED UPS NTOT=ITOT(LF) C IF(.NOT.BRMX)THEN !DW/PWB DO I=IE0,IFIN0 XMANT(I)=OMEGIN(I,NTR) ENDDO OMEGA(0)=XMANT(I00) C IF(ITYPIN.EQ.1)THEN DO I=IE0,IFIN EMSH(I)=XB(I)*ETHRSH ENDDO IF(NZED.NE.NELC)CALL REMAPX(EMSH,XMANT,IFIN) !REMAP OMEGA ENDIF C OMGINF=OMEGIN(IOMINFO,NTR) C CALL OMGINT(EMSH,XMANT,I00,IFIN,ETHRSH,TINF,OMGINF X ,EJ,OMEGA,IE0,NTOT) C OMEGIN(IOMINF,NTR)=OMGINF C DO I=IE0,NTOT OMEGIN(I,NTR)=(OMEGA(I)+OMEGA(I-1))/2 !AVERAGE ENDDO ENDIF C DO I=IE0,NTOT OMEGA(I)=OMEGIN(I,NTR) EMSHI=EJ(I)/TAZ+EF !REALLY SHOULD TEST FINAL E IF(EMSHI.LT.EFITF)IEFT=I IF(EMSHI.LT.EASYM)IEASYM=I ENDDO IE=NTOT-IE0+1 C ELSE C C SET-UP FINAL STATE MESH ENERGY, AND TRANSFER OMEGIN TO OMEGA C IE=0 IOM=0 DO I=IE0,IFIN !EXCLUDES ANY INF E IOM=IOM+1 IF(EMSH(I).GE.EF+eshft)THEN IE=IE+1 EJ(IE)=(EMSH(I)-EF)*TAZ OMEGA(IE)=OMEGIN(IOM,NTR) IF(EMSH(I).LT.EFITF)IEFT=IE IF(EMSH(I).LT.EASYM)IEASYM=IE ENDIF ENDDO C ENDIF C IF(IE.EQ.0)GO TO 140 !ALL CLOSED IN THIS BATCH MXETOT=IE C C PICK-UP CONTRIBUTION FROM THRESHOLD TO FIRST TABULATED POINT C (RAW OMEGA ONLY) C MXE0=1 IF(UPS(1,1,NTRAN).EQ.DZERO)THEN !FIRST TIME IF(EJ(1).LT.ABVTHR*TAZ)THEN MXE0=0 EJ(0)=0 IF(NZED.EQ.NELC)THEN !NEUTRAL IS ZERO OMEGA(0)=0 ELSE IF(IE.EQ.1)THEN OMEGA(0)=OMEGA(1) ELSE T=(EJ(0)-EJ(1))/(EJ(2)-EJ(1)) OMEGA(0)=(DONE-T)*OMEGA(1)+T*OMEGA(2) ENDIF ENDIF ELSEIF(LI.EQ.LIMN)THEN WRITE(6,640)LF,EJ(1) ENDIF ENDIF C C----------------------------------------------------------------------- C C BINNED COLLISION STRENGTHS C C----------------------------------------------------------------------- C IF(ITYPE.EQ.5)THEN C NTOT=ITOT(LF) C CALL OMGBIN(EJ,OMEGA,MXE0,MXETOT,IEINCR,EBIN,UPS(1,1,NTRAN),NTOT) C IF(BLAST)THEN AR(LI,LF)=DZERO AR(LF,LI)=DZERO IF(BINF)THEN OMGINF=OMEGIN(IOMINF,NTR) IF(OMGINF.NE.DZERO)THEN IF(BBORN.AND.OMGINF.LT.DZERO.OR..NOT.BBORN)THEN !DIPOLE CINF=ABS(OMGINF)/LOG(TINF) !CINF=4*SLIN/3 DELE=(EF-EI)*TAZ !AS ETHRSH MAYBE ADJUSTED AR(LI,LF)=CONRAD*(DELE**3)*CINF/JXTWO(LOF) AR(LF,LI)=-CINF ELSE !BORN IFLAGB=1 !FLAG LIMIT FOUND AR(LF,LI)=OMGINF ENDIF ENDIF ENDIF ENDIF C GO TO 140 !WE ARE DONE C ENDIF C C----------------------------------------------------------------------- C C SET BURGESS-TULLY FINITE E INTERP POINT (USER MAY RESTRICT LAST E) C IF(BINF)THEN IF(IEASYM.EQ.0)IEASYM=IE !USE LAST FINITE E OMGASYM=OMEGA(IEASYM) UASYME=EJ(IEASYM)/ETHRSH ENDIF C C SET BURGESS-TULLY INFINITE E INTERP POINT, AND SET TRANSITION NTYPE C IF(BLAST)THEN C NTYPE=0 !BURGESS TRANSITION TYPE CINF=DZERO AR(LF,LI)=DZERO IF(ITYPE0.NE.-3)AR(LI,LF)=DZERO C IF(BINF)THEN !SET-UP TO ADD-IN HIGH-E CONTRIB C OMGINF=OMEGIN(IOMINF,NTR) C IF(OMGINF.NE.DZERO)THEN IF(BBORN.AND.OMGINF.LT.DZERO.OR..NOT.BBORN)THEN !DIPOLE NTYPE=1 ISPN=-1 CINF=ABS(OMGINF)/LOG(TINF) !CINF=4*SLIN/3 AR(LF,LI)=-CINF CASYM=OMGASYM/LOG(UASYME+EXP(DONE)) IF(ITYPE0.NE.-3)THEN DELE=(EF-EI)*TAZ !AS ETHRSH MAYBE ADJUSTED AR(LI,LF)=CONRAD*(DELE**3)*CINF/JXTWO(LOF) ENDIF ELSE !BORN IFLAGB=1 !FLAG LIMIT FOUND NTYPE=2 ISPN=0 CINF=OMGINF AR(LF,LI)=CINF CASYM=OMGASYM ENDIF ELSE IF(ITCC.LE.0.AND..NOT.BBORN.AND. !LS/JK X IABS(ISAT(LF)).EQ.IABS(ISAT(LI)))THEN ISPN=0 ELSE ISPN=IABS(IBORN) ENDIF ALF=ISPN CASYM=OMGASYM IF(ISPN.EQ.0)THEN !CONSTANT, ONLY DIPOLE LIMITS NTYPE=2 ELSE !FORBIDDEN NTYPE=3 C IF(BDW)IEFT=IEASYM-3 !IF ORIG. NOT INTERPOLATED MXE2=IEASYM-IEFT IF(IRMPS.GE.0.AND.MXE2.GT.0)THEN IF(MXE2.LT.3)THEN IF(BRMX)THEN WRITE(6,670)MXE2,MXEIN,IBORN STOP '*** INSUFFICIENT HIGH ENERGY POINTS' ELSE IRMPS=-1 WRITE(6,671)IBORN ENDIF ELSE MXE1=IEASYM-MXE2/3 TJ=OMEGA(MXE1)*OMEGA(IEASYM) IF(TJ.GT.DZERO)THEN TJ=LOG(OMEGA(MXE1)/OMEGA(IEASYM)) X /LOG(EJ(IEASYM)/EJ(MXE1)) ISPN0=NINT(TJ) IF(ISPN0.NE.ISPN)THEN IF(ISPN0.LT.ISPN)THEN ISPN=ISPN-1 T=ISPN ALF=MAX(T,TJ) ENDIF IF(ISPN0.GT.ISPN)THEN ISPN=ISPN+1 T=ISPN ALF=MIN(T,TJ) ENDIF ELSE ALF=TJ ENDIF ENDIF ENDIF ENDIF CASYM=CASYM*UASYME**ALF ENDIF ENDIF C ELSE !LOW-E UPS1 CONTRIB ONLY ENDIF ENDIF C C----------------------------------------------------------------------- C C CALCULATE EFFECTIVE COLLISION STRENGTHS AS A FUNCTION OF TEMPERATURE C C----------------------------------------------------------------------- C C EXPLICIT LOW ENERGY CONTRIBUTION C DO IT=1,NUMTMP C TEMP=TKADAS(IT)/CONRYK !KELVIN TO RYD C IF(ICON.EQ.0)THEN CALL MAX1(EJ,OMEGA,MXE0,MXETOT,IEINCR,TEMP,TOLH,NDELBT X ,DOWN1) ELSEIF(ICON.EQ.1)THEN CALL KAP1(EJ,ETHRSH,OMEGA,MXE0,MXETOT,IEINCR,TEMP,TOLH,NDELBT X ,KAPPA,DOWN1,UP1) ELSEIF(ICON.EQ.2)THEN CALL DRY1(EJ,ETHRSH,OMEGA,MXE0,MXETOT,IEINCR,TEMP,TOLH,NDELBT X ,XDRY,DOWN1,UP1) ELSEIF(ICON.EQ.3)THEN CALL NUM1(EJ,ETHRSH,OMEGA,MXE0,MXETOT,IEINCR,TEMP,TOLH,NDELBT X ,NENG37,E37(1,IT),F37(1,IT),DOWN1,UP1) ENDIF C IF(IDUP1.EQ.1)UPS(IT,1,NTRAN)=UPS(IT,1,NTRAN)+DOWN1 IF(IDUP.EQ.2)UPS(IT,2,NTRAN)=UPS(IT,2,NTRAN)+UP1 C ENDDO C C INTERPOLATED HIGH ENERGY PART C IF(BINF)THEN C EJI=EJ(MXETOT) UASYM=UASYME+DONE c c write(6,*)lf,li,mxetot,eji,ethrsh,ieasym,uasym,omgasym,casym,cinf C DO IT=1,NUMTMP C TEMP=TKADAS(IT)/CONRYK !KELVIN TO RYD C IF(ICON.EQ.0)THEN CALL MAX2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMP X ,DOWN2) ELSEIF(ICON.EQ.1)THEN CALL KAP2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMP X ,KAPPA,DOWN2,UP2) ELSEIF(ICON.EQ.2)THEN CALL DRY2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMP X ,XDRY,DOWN2,UP2) ELSEIF(ICON.EQ.3)THEN CALL NUM2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMP X ,NENG37,E37(1,IT),F37(1,IT),DOWN2,UP2) ENDIF C IF(IDUP1.EQ.1)UPS(IT,1,NTRAN)=UPS(IT,1,NTRAN)+DOWN2 IF(IDUP.EQ.2)UPS(IT,2,NTRAN)=UPS(IT,2,NTRAN)+UP2 C ENDDO C ENDIF C C----------------------------------------------------------------------- C 140 continue C C*********************************************************************** C C *** END LOOP OVER ALL TRANSITIONS *** C TO BIN OR CONVOLUTE COLLISION STRENGTHS C C*********************************************************************** C C IF(.NOT.BLAST)THEN IF(.NOT.BFORM)THEN !MINOR REFINEMENT BACKSPACE(7) IE1=IE1-1 ENDIF GO TO 60 !GO AND GET NEXT BATCH OF ENERGIES ENDIF C C END OF ENERGIES, CHECK BORN CONSISTANCY. C IF(BINF.AND.BBORN.AND.IFLAGB.EQ.0)THEN WRITE(6,680) STOP 'BORN LIMITS EXPECTED, BUT NONE FOUND' ENDIF C C APPLY RELATIVISTIC (JUTTNER) CORRECTION TO DISTRIBUTION. C IF(IREL.NE.0.AND.ITYPE.NE.5)THEN C NU=2 MU=4*NU*NU C DO IT=1,NUMTMP C T=TKADAS(IT)/CONRYK !KELVIN TO RYD THETA=T*DALF/2 !KT(a.u.)/C**2 T8=THETA/8 C START-OFF WITH ABRAMOWITZ & STEGUN 9.7.2, AVOIDS SMALL THETA OVERFLOW KSUM=-5/LOG10(THETA) KSUM=KSUM+1 KSUM=MIN(KSUM,10) T=DONE SUM=DONE DO K=1,KSUM T=T*(MU-(2*K-1)**2)*T8/K SUM=SUM+T c write(*,*)k,t,sum,done/sum ENDDO C IF NOT CONVERGED, SAFE TO EVALUATE EXPLICITLY NOW AS THETA LARGE IF(ABS(T/SUM).GT.D1M4)then c write(*,*)tkadas(it),theta,sum,t,t/sum TT=DONE/THETA FX=SQRT(PI*THETA/2)/(BESSK(2,TT)*EXP(TT)) FREL(IT)=FX ELSE FREL(IT)=DONE/SUM ENDIF c write(*,*)it,ksum,theta,t/sum,done/sum,fx C ENDDO C DO NTRAN=1,NTRANS C DO IT=1,NUMTMP C IF(IDUP1.EQ.1)UPS(IT,1,NTRAN)=UPS(IT,1,NTRAN)*FREL(IT) IF(IDUP.EQ.2)UPS(IT,2,NTRAN)=UPS(IT,2,NTRAN)*FREL(IT) C ENDDO C ENDDO ENDIF C C CHECK FOR PROBLEMS (NEGATIVE UPSILONS) C IFLAG1=0 IFLAG2=0 C DO NTRAN=1,NTRANS C IF(IDUP1.EQ.1)THEN DO IT=1,NUMTMP IF(UPS(IT,1,NTRAN).LT.DZERO)THEN WRITE(6,672)LFN(NTRAN),LIN(NTRAN) X ,(UPS(II,1,NTRAN),II=1,NUMTMP) IFLAG1=IFLAG1+1 GO TO 145 ENDIF 145 ENDDO ENDIF C IF(IDUP.EQ.2)THEN DO IT=1,NUMTMP IF(UPS(IT,2,NTRAN).LT.DZERO)THEN WRITE(6,672)LFN(NTRAN),LIN(NTRAN) X ,(UPS(II,2,NTRAN),II=1,NUMTMP) IFLAG2=IFLAG2+1 GO TO 146 ENDIF 146 ENDDO ENDIF C ENDDO C IF(IFLAG1.GT.0)THEN WRITE(6,673)IFLAG1 ENDIF IF(IFLAG2.GT.0)THEN WRITE(6,673)IFLAG2 ENDIF IF(IFLAG1+IFLAG2.GT.0)THEN WRITE(0,*)' *** WARNING: THERE ARE NEGATIVE UPSILONS ***' WRITE(0,*)' *** SEE adasexj.out FOR DETAILS ***' ENDIF C C C*********************************************************************** C C BEGIN I/O OF ADAS adf04 FILE C C*********************************************************************** C 170 CONTINUE !ENTRY POINT IF READING adf04_om C C----------------------------------------------------------------------- C C SET-UP I/O FORMATS C IF(NSTATES.LT.1000)THEN IF(KMAX.LE.20)THEN IF(BEXP)THEN F761='(F5.2,I5,8X,20(1PE10.2))' F762='(2I4,22(1PE10.2))' ELSE IF(BREAD)THEN F761='(F5.2,I5,6X,20(F5.2,I3))' F762='(2I4,22(F5.2,I3))' ELSE F761='(F5.2,I5,6X,20( A5,A3))' F762='(2I4,22( A5,A3))' ENDIF ENDIF ELSE IF(BEXP)THEN F761='(F5.2,I5,8X,20(1PE10.2)/(18X,20(1PE10.2)))' F762='(2I4,21(1PE10.2)/(18X,20(1PE10.2)))' !WRAP LAST POINT ELSE IF(BREAD)THEN F761='(F5.2,I5,6X,20(F5.2,I3)/(16X,20(F5.2,I3)))' F762='(2I4,21(F5.2,I3)/(16X,20(F5.2,I3)))'!WRAP LAST POINT ELSE F761='(F5.2,I5,6X,20( A5,A3)/(16X,20( A5,A3)))' F762='(2I4,21( A5,A3)/(16X,20( A5,A3)))'!WRAP LAST POINT ENDIF ENDIF ENDIF IF(BEXP)THEN C300='(2I4,10X,21(2X,A8))' ELSE C300='(2I4, 8X,21( A8))' ENDIF ELSE IF(KMAX.LE.20)THEN IF(BEXP)THEN F761='(F5.2,I5,10X,20(1PE10.2))' F762='(2I5,22(1PE10.2))' ELSE IF(BREAD)THEN F761='(F5.2,I5, 8X,20(F5.2,I3))' F762='(2I5,22(F5.2,I3))' ELSE F761='(F5.2,I5, 8X,20( A5,A3))' F762='(2I5,22( A5,A3))' ENDIF ENDIF ELSE IF(BEXP)THEN F761='(F5.2,I5,10X,20(1PE10.2)/(20X,20(1PE10.2)))' F762='(2I5,21(1PE10.2)/(20X,20(1PE10.2)))' !WRAP LAST POINT ELSE IF(BREAD)THEN F761='(F5.2,I5, 8X,20(F5.2,I3)/(18X,20(F5.2,I3)))' F762='(2I5,21(F5.2,I3)/(18X,20(F5.2,I3)))'!WRAP LAST POINT ELSE F761='(F5.2,I5, 8X,20( A5,A3)/(18X,20( A5,A3)))' F762='(2I5,21( A5,A3)/(18X,20( A5,A3)))'!WRAP LAST POINT ENDIF ENDIF ENDIF IF(BEXP)THEN C300='(2I5,10X,21(2X,A8))' ELSE C300='(2I5, 8X,21( A8))' ENDIF ENDIF C C----------------------------------------------------------------------- C C I/O ENERGIES/LEVELS+TEMPERATURES C IF(BREAD)THEN !READ C IF(BEXP)THEN READ(IREAD,F761,err=180)DUM,IDUM,(EMSH(I),I=1,KMAX) ELSE READ(IREAD,F761,err=180)DUM,IDUM,(EMSH(I),IEXP(I),I=1,KMAX) ENDIF 180 CONTINUE C !RE-CONSTRUCT ENERGIES DO I=1,KMAX !MAX POSS (IF LINE FULL) IF(.NOT.BEXP)EMSH(I)=EMSH(I)*DTEN**IEXP(I) IF(EMSH(I).EQ.0.AND.I.GT.1)GO TO 190 !FOR DW ITYPE=5 ENDDO C 190 NTOT=I-1 !ACTUAL (MAX) NUMBER ITOT(1)=NTOT C IF(BEQUAL)THEN !E.G. DW/PWB DO N=2,NSTATES ITOT(N)=NTOT ENDDO ENDIF C IF(ITYPIN.EQ.1)THEN !EMSH=X TOO LARGE IN GENERAL EFIN=TMAX0/(2*CONRYK) !EXCEPT NEAR-DEGENERATE LEVELS ELSE EFIN=EMSH(NTOT)+(ENERG(1+IONE)-ENERG(1))/CONRCM EFIN=EFIN/TAZ ENDIF C BINF=.TRUE. !MAY ALL BE ZERO, BUT FIELD WILL EXIST MXE=NTOT+1 !ALLOW SPACE FOR INF E OMG C EMSH(MXE)=EINF !NOW USE EINF EXPLICITLY, NOT EMSH(MXE) TINF=EINF*TAZ IOMINF=MXE !SINCE ONLY ONE BATCH HERE BLAST=.TRUE. C ITYPE0=-5 !RE-SET C ELSE !WRITE C F570="(A2,'+',I2,2I10,F15.0,'(',A2,')')" IFIPOT=FIPOT IF(ABS(FIPOT-IFIPOT).GT.1.D-4) X F570(21:21)='4' WRITE(9,F570)IEL,IQ,NZED,IQ1,FIPOT,IONTRM !HEADER C DO L=1,NSTATES !LEVELS IF(BSHORT)THEN WRITE(9,F580)L,CONF(L),ISAT(L),LAT(L),FJT(L),ENERG(L) ELSE WRITE(9,F581)L,XCONF(L),ISAT(L),LAT(L),FJT(L),ENERG(L) ENDIF ENDDO C WRITE(9,590)CARD !TERMINATOR C FIQ1=IQ1 IF(BEXP)THEN IF(ITYPE.EQ.5)THEN WRITE(9,F761)FIQ1,ITYPE,(EBIN(N),N=1,KMAX) !DROP ZERO ELSE WRITE(9,F761)FIQ1,ITYPE,(TKADAS(IT),IT=1,KMAX) ENDIF ELSE REWIND(98) IF(ITYPE.EQ.5)THEN WRITE(98,730)(EBIN(N),N=1,KMAX) !DROP ZERO ELSE WRITE(98,730)(TKADAS(IT),IT=1,KMAX) ENDIF REWIND(98) READ(98,740)(CMANT(IT),CEXP(IT),IT=1,KMAX) WRITE(9,F761)FIQ1,ITYPE,(CMANT(IT),CEXP(IT),IT=1,KMAX) ENDIF C ENDIF C C LOOP-OVER ALL TRANSITIONS TO I/O (OMEGA AND/OR UPSILON) FROM/TO adf04 C IF(BREAD)THEN !GET TYPE 5 INPUT TO PROCESS C NTRAN=0 NTRANS=(NSTATES*(NSTATES+1-2*IONE))/2 !MAX POSS ON FILE C DO NTR=1,NTRANS C 192 READ(IREAD,C300)J0,I0 BACKSPACE(IREAD) C IF(J0.LT.0)GO TO 200 !WE ARE DONE C KF=MAX(J0,I0) !GET UPPER KMAX=ITOT(KF) C IF(KMAX.EQ.0)THEN !NEED TO DETERMINE LENGTH C NT=1 NRD=21 C 194 READ(IREAD,C300)JJ,II,(CANT(I),I=1,NRD) C IF(JJ*JJ.NE.JJ*J0.OR.II*II.NE.II*I0)THEN !NEXT TRANSITION I=2 !INF WAS IN FINAL POSITION READ IF(CANT(21).EQ.CBLNK8)I=1 !BUT NOT IN POS 21 GO TO 195 ENDIF C NRD=20 DO I=1,NRD IF(CANT(I).EQ.CBLNK8)GO TO 195 !FIRST BLANK TERMINATES ENDDO KMAX=KMAX+NRD NT=NT+1 I0=II J0=JJ GO TO 194 C !WE HAVE THE CORRECT LENGTH 195 KMAX=KMAX+I-2 !EXC. INFINITE POINT ITOT(KF)=KMAX C DO N=1,NT BACKSPACE(IREAD) ENDDO C ENDIF C IF(BEXP)THEN READ(IREAD,F762)JJ,II,ARIF,(OMEGA(K),K=1,KMAX),ARFI ELSE READ(IREAD,F762)JJ,II,(XMANT(K),IEXP(K),K=0,KMAX+1) ENDIF C IF(JJ.GE.II)THEN LI=II LF=JJ ELSE !YOU NEVER KNOW LI=JJ LF=II ENDIF C IF(LF.GT.NUPMX.OR.LF.LT.NUPMN.OR. X LI.GT.NLOWMX.OR.LI.LT.NLOWMN)GO TO 192 C NTRAN=NTRAN+1 C LIN(NTRAN)=LI LFN(NTRAN)=LF C DO I=1,IDUP DO K=1,NUMTMP UPS(K,I,NTRAN)=-D1M30 !SO THRESHOLD BIN NOT RE-ADJUSTED ENDDO ENDDO C IF(BEXP)THEN AR(LI,LF)=ARIF AR(LF,LI)=ARFI DO K=1,KMAX OMEGIN(K,NTRAN)=OMEGA(K) ENDDO ELSE AR(LI,LF)=XMANT(0)*DTEN**IEXP(0) AR(LF,LI)=XMANT(KMAX+1)*DTEN**IEXP(KMAX+1) DO K=1,KMAX OMEGIN(K,NTRAN)=XMANT(K)*DTEN**IEXP(K) ENDDO ENDIF IF(BINF)THEN CINF=AR(LF,LI) IF(CINF.LT.DZERO)CINF=CINF*LOG(TINF) OMEGIN(IOMINF,NTRAN)=CINF ENDIF c c fix top-up failure on high-e near/degenerate Born allowed transitions c if(btop.and.ar(lf,li).gt.dzero)then !Born de=(energ(lf)-energ(li))/conrcm de=max(de,d1m30) om=omegin(2,ntran) do k=3,kmax !skip threshold point rat=emsh(k)/de om1=om om=omegin(k,ntran) if(rat.gt.xbtest.and.om.lt.0.8*om1.and.om1.lt.ar(lf,li)) x then if(lrglam.gt.2*rat)then i7x=73 else i7x=74 endif write(i7x,762)jj,ii,k,emsh(k),de,rat, x (omegin(kk,ntran),kk=1,kmax),ar(lf,li) 762 format(2i5,i3,22(1pe10.2)) cdebug call flush(i7x) kmin=k-1 !max(2,k-1) if no threshold skip go to 196 endif enddo kmin=kmax+1 if(om1.gt.1.01*ar(lf,li))then !allow for 3 s.f. round-off if(om.lt.0.99*ar(lf,li))kmin=kmax elseif(om1.lt.0.99*ar(lf,li))then if(om.lt.0.99*om1)kmin=kmax endif if(kmin.eq.kmax)then de=(energ(lf)-energ(li))/conrcm de=max(de,d1m30) rat=emsh(kmax)/de write(75,762)jj,ii,kmax,emsh(kmax),de,rat, x (omegin(kk,ntran),kk=1,kmax),ar(lf,li) endif 196 do kk=kmin,kmax omegin(kk,ntran)=ar(lf,li) enddo endif C ENDDO C 200 NTRANS=NTRAN C 201 READ(IREAD,F762)IDUM !TERMINATOR IF(IDUM.NE.-1)GO TO 201 READ(IREAD,F762)IDUM,IDUM !TERMINATOR I=NMLNS+1 210 READ(IREAD,780,END=220)COMMENT(I) !COMMENTS CAR1=COMMENT(I)(1:1) IF(CAR1.NE.'C'.AND.CAR1.NE.'c'.AND.CAR1.NE.'!')GO TO 220 I=I+1 IF(I.GT.NLNMX)THEN WRITE(6,540)NLNMX STOP 'INCREASE THE PARAMETER NLNMX (OR SHORTEN COMMENTS)' ENDIF GO TO 210 220 NMLNS=I-1 C ELSE !OUTPUT TYPE 5 OR 3/4 C DO NTRAN=1,NTRANS C LI=LIN(NTRAN) LF=LFN(NTRAN) C AR(LI,LF)=MAX(AR(LI,LF),D1M30) C IF(ITYPE.EQ.5)KMAX=ITOT(LF) C DO I=IDUP1,IDUP DO K=KMAX,1,-1 IF(UPS(K,I,NTRAN).LT.DMIN)THEN DO J=K,1,-1 UPS(J,I,NTRAN)=DMIN ENDDO GO TO 230 ENDIF ENDDO 230 CONTINUE ENDDO C IF(IDUP.EQ.2)THEN DO K=KMAX,1,-1 IF(UPS(K,2,NTRAN).GT.DMAX)THEN DO J=K,1,-1 UPS(J,2,NTRAN)=DMAX ENDDO GO TO 240 ENDIF ENDDO 240 CONTINUE ENDIF C IF(BEXP)THEN IF(IDUP1.EQ.1) X WRITE(9,F762)LF,LI,AR(LI,LF),(UPS(K,1,NTRAN),K=1,KMAX) X ,AR(LF,LI) IF(IDUP.EQ.2) X WRITE(9,F762)LI,LF,D1M30,(UPS(K,2,NTRAN),K=1,KMAX) X ,AR(LF,LI) ELSE REWIND(98) IF(IDUP1.EQ.1) X WRITE(98,730)AR(LI,LF),(UPS(K,1,NTRAN),K=1,KMAX),AR(LF,LI) IF(IDUP.EQ.2) X WRITE(98,730)D1M30,(UPS(K,2,NTRAN),K=1,KMAX),AR(LF,LI) REWIND(98) IF(IDUP1.EQ.1)THEN READ(98,740)(CMANT(K),CEXP(K),K=1,KMAX+2) WRITE(9,F762)LF,LI,(CMANT(K),CEXP(K),K=1,KMAX+2) ENDIF IF(IDUP.EQ.2)THEN READ(98,740)(CMANT(K),CEXP(K),K=1,KMAX+2) WRITE(9,F762)LI,LF,(CMANT(K),CEXP(K),K=1,KMAX+2) ENDIF ENDIF C ENDDO C WRITE(9,F762)-1 !TERMINATOR WRITE(9,F762)-1,-1 !TERMINATOR C !COMMENTS WRITE(9,770) IF(BPWB.AND.NZED.NE.NELC)WRITE(9,771) IF(ITYPE.EQ.4)THEN IF(ICON.EQ.1)WRITE(9,772)KAPPA IF(ICON.EQ.2)WRITE(9,773)XDRY IF(ICON.EQ.3)WRITE(9,774) !NUMERICAL ENDIF IF(IREL.NE.0)WRITE(9,775) DO I=1,NMLNS WRITE(9,780)COMMENT(I) ENDDO WRITE(9,790)NAME,DATE C ENDIF C C*********************************************************************** C C END I/O OF ADAS adf04 FILE C C*********************************************************************** C C C NOW LOOP BACKUP AND CONVOLUTE OMEGAS WITH ELECTRON DISTRIB. FUNCTION C (UNLESS USER INPUT NBIN.GT.0 FOR JUST BIN OUTPUT.) C NEED TO RE-ADJUST (FOOL) VARIOUS VARIABLES. C IF(ITYPE0.EQ.-5)THEN C ITYPE=3 ITYPE0=-ITYPE IF(ICON.GT.0)ITYPE=4 C IF(ICON.EQ.0)THEN !NO TIME SAVING WITH TRAP C NDELBT=-NDELBT-1 !SO DISABLE NDELBT=-NMPTS-1 ELSEIF(ICON.EQ.1)THEN !NO TIME SAVING WITH TRAP C NDELBT=-NDELBT-1 !SO DISABLE NDELBT=-NMPTS-1 ELSEIF(ICON.EQ.2)THEN !ALLOW TRAP FOR SPEED NDELBT=-NDELBT-1 ELSEIF(ICON.GE.3)THEN !ALLOW TRAP FOR SPEED NDELBT=-NDELBT-1 ENDIF C IEINCR=1 !DON'T SKIP ENERGY BINS IE0=1 C EJ(IE0-1)=EBIN(0) !=0 AS CODED, BUT IN PRINCIPLE NTOT=ITOT(1) IFIN=NTOT !LAST FINITE BIN ENERGY INDEX C IF(BREAD)THEN !BINNED OMEGAS ALREADY IN PLACE C BREAD=.FALSE. IF(BRMX)THEN !JUGGLE MESH DO I=IE0,IFIN EBIN(I)=EMSH(I) ENDDO ELSE NTOT=NBIN !RESTORE ENDIF C ELSE !TRANSFER BACK C CLOSE(9) OPEN(9,FILE='adf04_ups') !NEW OUTFILE C MXE=IFIN IF(BINF)THEN !SET INFINITE ENERGY OMG AT LAST BIN+1 MXE=MXE+1 IOMINF=MXE !SINCE ONLY ONE BATCH HERE ENDIF C DEALLOCATE (OMEGIN) ALLOCATE (OMEGIN(MXE,NTRANS)) C DO NTRAN=1,NTRANS DO I=IE0,IFIN OMEGIN(I,NTRAN)=UPS(I,1,NTRAN) ENDDO IF(BINF)THEN LI=LIN(NTRAN) LF=LFN(NTRAN) CINF=AR(LF,LI) IF(CINF.LT.DZERO)CINF=CINF*LOG(TINF) OMEGIN(IOMINF,NTRAN)=CINF ENDIF ENDDO C IF(IDUP0.GT.IDUP)THEN IDUP=IDUP0 DEALLOCATE (UPS) ALLOCATE (UPS(NUMTMP,IDUP,NTRANS)) ENDIF C DO N=1,NTRANS DO I=1,IDUP DO IT=1,NUMTMP UPS(IT,I,N)=D1M30 !SO THRESHOLD BIN NOT RE-ADJUSTED ENDDO ENDDO ENDDO C ENDIF C IF(BRMX.AND.NMLNS+4.LT.NLNMX)THEN N=NMLNS+1 COMMENT(N)='C ' N=N+1 COMMENT(N)='C '// X '*** ATTENTION: THESE UPSILONS HAVE BEEN DETERMINED FROM THE' N=N+1 COMMENT(N)='C '// X '*** ENERGY-AVERAGED OMEGAS. MORE PRECISE VALUES MAY BE OBTAINED' N=N+1 COMMENT(N)='C '// X '*** BY CONVOLUTING THE ORIGINAL RAW COLLISION STRENGTHS.' NMLNS=N ENDIF NAME='NAME: ADASEXJ.F' CALL DATE_AND_TIME(DATE8) DATE='DATE: '//DATE8(7:7)//DATE8(8:8)//'/'//DATE8(5:5)// X DATE8(6:6)//'/'//DATE8(3:3)//DATE8(4:4) C C WE HAVE RMX/DW/PWB COLLISION STRENGTHS SO HEAD BACK UP AND CONVOLUTE C IF(NUMTMP.GT.0)THEN IF(BRMX)THEN GO TO 70 ! >>>>>>>>>>>>>>>>>> ELSE GO TO 65 ! >>>>>>>>>>>>>>>>>> ENDIF ENDIF C ENDIF C C DE-ALLOCATE C DEALLOCATE (AR,UPS,OMEGIN) C C TIME C DUM=DTIME(TARRY) TIME=TARRY(1) C TIME=TIME/60 WRITE(6,800)TIME C C----------------------------------------------------------------------- C C FORMATS C 500 FORMAT(A15,F5.1,F11.0,I5) 501 FORMAT(3A15) 502 FORMAT(5X,3A15) 510 FORMAT(A30) 520 FORMAT(/'*** CONVOLUTION ELECTRON DISTRIBUTION IS MAXWELLIAN') 521 FORMAT(/'*** CONVOLUTION ELECTRON DISTRIBUTION IS KAPPA, WITH K=' X,F7.2) 522 FORMAT(/'*** CONVOLUTION ELECTRON DISTRIBUTION IS DRUYVESTEYN,' X,' WITH X=',F6.1) 523 FORMAT(/' CONVOLUTION ELECTRON DISTRIBUTION IS NUMERICAL (adf37)') 530 FORMAT(A74) 531 FORMAT(/'*** ATTENTION: YOU HAVE FLAGGED LS-COUPLING (ITCC < 0) ' X,'BUT HAVE SET IBORN = 0'/'*** RECOMMEND SETTING IBORN = -2 FOR ' X,'EXTRAPOLATION OF FORBIDDEN TRANSITIONS'/'***NON-DIPOLE ALLOWED' X,' WILL STILL BE TAKEN TO BE CONSTANT AT HIGH ENERGIES') 532 FORMAT(/'*** ATTENTION: YOU HAVE FLAGGED IC-COUPLING (ITCC > 0) ', X'BUT HAVE SET IBORN =',I3,' (<0)'/'*** RECOMMEND SETTING IBORN= 0' X,' FOR CONSTANT EXTRAPOLATION OF NON-DIPOLE ALLOWED TRANSITIONS'/ X ' ***(CURRENTLY, THERE IS NO WAY TO DISTINGUISH FORBIDDEN' X,' TRANSITIONS.)') 533 FORMAT(/'*** ATTENTION: YOU HAVE SET IBORN =',I3 X/'***RECOMMEND SETTING |IBORN| = 0 FOR EXTRAPOLATION OF' X,' FORBBIDEN TRANSITIONS.') 534 FORMAT(/'*** STRONG WARNING: YOU HAVE CHOSEN TRAPEZOIDAL RULE' X ,' FROM E=0!') 535 FORMAT('*** ERROR: DIMENSION NMBIN MUST BE AT LEAST NMTMP:',2I5) 536 FORMAT(/5X,'*** ATTENTION, THE NUMBER OF TEMPERATURES HAS BEEN' X ,' REDUCED ON DIMENSION GROUNDS FROM',I4,' TO',I4/) 537 FORMAT(/5X,'*** ATTENTION, THE NUMBER OF TEMPERATURES EXCEEDS' X ,' THE MAX (20) THAT CAN BE HANDLED BY ADAS'/5X,'*** THE ADF04' X ,' SHOULD NOT BE PUT INTO THE ADAS DATABASE.'/) 538 FORMAT(/5X,'*** ATTENTION, THE NUMBER OF TEMPERATURES HAS BEEN' X ,' REDUCED TO THE ADAS DEFAULT I.E. FROM',I4,' TO 13'/) 540 FORMAT(/5X, 1 '********************** FATAL INPUT ERROR! *******************'/ 2 5X,'THE COMMENTARY IS LONGER THAN ',I2,' LINES. INCREASE THE '/ 3 5X,'VALUE OF NLNMX IN THE CODE OR SHORTEN THE COMMENTS'/5X, 4 '************************************************************'/) 541 FORMAT('*** ERROR: adf04_om IS NOT A COLLISION STRENGTH FILE:' X,' ITYPE=',I3) 542 FORMAT(I5,1X,A15,4X,I1,1X,I1,1X,F4.1,1X,F21.4) 543 FORMAT(I5,1X,A31,4X,I1,1X,I1,1X,F4.1,1X,F21.4) 547 FORMAT(/'*** WARNING, YOUR MINIMUM TEMPERATURE IS LESS THAN THE' X,' ADAS DEFAULT ONE: ',1P,2E10.2) 548 FORMAT(/'*** STRONG WARNING, YOU LIKELY NEED TO INCREASE NDELBT' X,' ABOVE ITS DEFAULT VALUE OF',I5) 550 FORMAT(/5X, 1 '****************** FATAL OMEGA INPUT ERROR! ****************'/ 2 5X,'THE NUMBER OF STATES IN THE CLOSE-COUPLING CALCULATION IS '/ 2 5X,'LARGER THAN THE CURRENT DIMENSION. INCREASE THE PARAMETER'/ 3 5X,' NSTMX TO ',I3/5X, 4 '************************************************************'/) 551 FORMAT(/5X, 1 '****************** FATAL INPUT ERROR! ****************'/ 2 5X,'THE NUMBER OF adf04 TERMS/LEVELS IS LARGER THAN THE CURRENT'/ 2 5X,' DIMENSION. INCREASE THE PARAMETER NSTMX TO ',5X,I3/5X, 4 '************************************************************'/) 555 FORMAT(/5X, 1 '********************** FATAL INPUT ERROR! *******************'/ 2 5X,'THE NUMBER OF STATES IN THE CLOSE-COUPLING CALCULATION IS', 3 ' EQUAL TO ',I5/5X,'BUT THE NUMBER OF LEVELS IN THE FILE', 4 ' adasexj.in IS EQUAL TO ',I5) 560 FORMAT(/5X, 1 '****************** FATAL OMEGA INPUT ERROR! ****************'/ 2 5X,'THE NUMBER OF ENERGY POINTS IN THE MESH IS LARGER THAN THE'/ 2 5X,' CURRENT DIMENSION. INCREASE THE PARAMETER NMPTS TO ',I6/5X, 4 '************************************************************'/) 571 FORMAT(A2,1X,I2,2I10,F15.0,1X,A2,1X) 575 FORMAT('MIS-MATCH OF TARGET STATES: adasex/j.in AND omega'/5I5) 576 FORMAT('MIS-MATCH OF TARGET ENERGIES: adasex/j.in AND omega'/ X 2F16.4) C 580 FORMAT(I5,1X,A15,2X,' (',I1,')',I1,'(',F4.1,')',F21.4) C 581 FORMAT(I5,1X,A31,2X,' (',I1,')',I1,'(',F4.1,')',F21.4) 585 FORMAT(/5X, 1 '****************** FATAL INPUT ERROR! ****************'/ 2 5X,'THE NUMBER OF BIN ENERGIES IS LARGER THAN THE CURRENT '/ 2 5X,'DIMENSION. INCREASE THE PARAMETER NMBIN TO ',I6/5X, 4 '************************************************************'/) 586 FORMAT(/' INTERNAL BIN ENERGY SET-UP: NBIN=',I6/ X' STARTING AT EMIN=',1PE12.4,3X,'WITH WIDTH FACTOR=',1PE12.4/ X' ENDING AT EMAX=',1PE12.4,3X,' WITH BIN WIDTH=',1PE12.4/ X'+LINEAR TO E(MAX)=',1PE12.4,3X,' WITH BIN WIDTH=',1PE12.4) 587 FORMAT(/' BINNED COLLISION STRENGTHS WRITTEN: NBIN=',I6/ X' STARTING AT EMIN=',1PE12.4,3X,'WITH WIDTH FACTOR=',1PE12.4/ X' ENDING AT EMAX=',1PE12.4,3X,' WITH BIN WIDTH=',1PE12.4/ X'+LINEAR TO E(MAX)=',1PE12.4,3X,' WITH BIN WIDTH=',1PE12.4) 590 FORMAT(A200) 595 FORMAT(/' BATCH ENERGY SET-UP OUT OF CONTROL, TRY OVERRIDING THE' X,' CURRENT VALUE OF MXEIN=',I5,' OR FORCE 1/E**2 FOR FORBIDDEN' X,' TRANSITIONS (SEE EFITF OR IRMPS)') C 600 FORMAT(1PE14.8,6(E11.3)/(14X,6(E11.3))) 604 FORMAT(/'*** OMEGA FILE DOES NOT START AT THRESHOLD?,',1P2E10.3/ X '*** CHECK IONE, OR RE-SET BOMTHR=.FALSE. TO AVOID EXTRAPOLATION' X,' TO THRESHOLD') 605 FORMAT(/'*** OMEGA FILE DOES NOT START AT THRESHOLD,',1P2E10.3/ X '*** RE-SET BOMTHR=.FALSE. TO AVOID EXTRAPOLATION TO THRESHOLD') 606 FORMAT(/'*** OMEGA FILE STARTS CLOSE TO THRESHOLD,',1P2E10.3/ X '*** RE-SET BOMTHR=.TRUE. TO ALLOW EXTRAPOLATION TO THRESHOLD') 610 FORMAT(/5X, 1 '****************** FATAL OMEGA INPUT ERROR! ****************'/ 2 5X,'THE omega FILE MUST CONTAIN OMEGAS FOR DIPOLE TRANSITIONS', 3 5X,' CALCULATED FROM A LONG RANGE DIPOLE APPROXIMATION',5X, 4 '************************************************************'/) 611 FORMAT(/'*** NO HIGH ENERGY POINTS, UPSILON INCOMPLETE!') 612 FORMAT(/'*** ERROR: LINEAR ENERGY MESH TOO SHORT:',2F12.3/ X'*** INCREASE NMBIN AND/OR SET MAXE FOR MORE DETAILS') 613 FORMAT(/' TOTAL NUMBER OF BIN ENERGIES=',I4/' LAST COMMON FINAL ', X'ENERGY=',1PE12.4) 614 FORMAT(/'*** INCREASE NMPIN TO ',I5,' HOLD ALL BINNED OMEGAS') 615 FORMAT(/' FORBIDDEN TRANSITION POWER LAW FROM THE ENERGY RANGE :', X 2F10.3,' RYD'/) 616 FORMAT(/5X,'****************************** WARNING!', 1 '******************************',/ 2 5X,'THE MAXIMUM TEMPERATURE WAS DECREASED FROM ',1PE9.2, 3 ' K TO ',E9.2,' K'/ 4 5X,'IN LIGHT OF THE MAXIMUM ENERGY CONTAINED IN THE omega', 5 ' FILE, A HIGHER',/ 6 5X,'TEMPERATURE COULD YIELD INACCURATE EFFECTIVE COLLISION' 7 ,' STRENGTHS',/ 9 5X,'*****************************************************', * '*****************'/) 618 FORMAT(/5X,'****************************** WARNING!', 1 '******************************',/ 2 5X,'THE MAXIMUM TEMPERATURE SHOULD BE DECREASED FROM ', 3 1PE9.2,' K TO ',E9.2,' K'/ 4 5X,'IN LIGHT OF THE MAXIMUM ENERGY CONTAINED IN THE omega', 5 ' FILE, A HIGHER',/ 6 5X,'TEMPERATURE COULD YIELD INACCURATE EFFECTIVE COLLISION' 7 ,' STRENGTHS',/ 9 5X,'*****************************************************', * '*****************'/) 620 FORMAT(/'***ERROR: IONIC DW/PWB MUST START AT THRESHOLD, EMSH(1)=' X,1PE9.2) 640 FORMAT(/'*** WARNING: OMEGA FILE INCOMPLETE? ***'/ X'*** THE FIRST SCATTERING ENERGY ABOVE LEVEL',I5/ X'*** IS AT',F10.4,' WHICH IS TOO LARGE FOR ADAS TEMPERATURES'/) 670 FORMAT(/'*** ENERGY MESH PROBLEM, INSUFFICIENT POINTS IN LAST', X ' BATCH, OR TOO FEW HIGH ENERGY POINTS: MXE2=',I4/ X '*** UNABLE TO DETERMINE HIGH ENERGY POWER FOR FORBIDDEN', X ' TRANSITIONS'/'*** EITHER RESET MXEIN, CURRENTLY =',I5/ X 'OR SET IRMPS.LT.0 TO FORCE A 1/E**IBORN BEHAVIOUR'/ X '*** CURRENTLY, IBORN=',I3,' (IBORN=2 RECOMMENDED)') 671 FORMAT(/'*** FORCING 1/E**IBORN=',I2,' BEHAVIOUR FOR FORBIDDEN' X,' TRANSITIONS, SINCE TOO FEW HIGH ENERGY POINTS TO DETERMINE IT') 672 FORMAT(2I4,1P,21(E10.2)/(8X,21(E10.2))) 673 FORMAT(/' *** WARNING: THERE ARE',I6,' TRANSITIONS WITH NEGATIVE' X,' UPSILONS ***') 680 FORMAT(/'*** ERROR: BORN LIMITS EXPECTED (IBORN.GT.0)', X ' BUT NONE HAVE BEEN FOUND - RE-RUN WITH IBORN=0') 730 FORMAT(22(1PE9.2)) 740 FORMAT(22(A5,1X,A3)) 750 FORMAT(F5.2,1X,I3) 770 FORMAT('C',79('-')/'C') 771 FORMAT('C PLANE WAVE BORN COLLISION STRENGTHS WERE REMAPPED' X/'C') 772 FORMAT('C *** UP/DOWNSILONS FOR A KAPPA DISTRIBUTION, WITH K=' X,F6.1/'C') 773 FORMAT('C *** UP/DOWNSILONS FOR A DRUYVESTEYN DISTRIBUTION,' X,' WITH X=',F6.1/'C') 774 FORMAT('C *** UP/DOWNSILONS FOR A NUMERICAL DISTRIBUTION...') 775 FORMAT('C JUTTNER RELATIVISTIC CORRECTION APPLIED TO THE', X ' DISTRIBUTION'/'C') 776 FORMAT(/' *** ATTENTION: YOU ARE APPLYING A JUTTNER RELATIVISTIC ' X ,'CORRECTION TO A NON-MAXWELLIAN DISTRIBUTION...') 780 FORMAT(A80) 790 FORMAT('C'/'C',5X,A30,4X,A30/'C',79('-')) 800 FORMAT(/'CPU TIME =',F9.3, ' MIN') C END C C*********************************************************************** C SUBROUTINE OMGBIN(E,OMG,IEI,IEF,IEINCR,EBIN,SBIN,NBIN) C C----------------------------------------------------------------------- C C BIN RAW COLLISION STRENGTHS READ FROM AN R-MATRIX omega FILE. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (D1M30=1.0D-30) !ADAS ZERO C DIMENSION E(0:*),OMG(0:*),EBIN(0:*),SBIN(*) C IF(SBIN(1).EQ.DZERO)SBIN(1)=D1M30 C IMN=IEI IMX=IEF-IEINCR C DO N=1,NBIN IF(E(IMN).LE.EBIN(N))GO TO 10 ENDDO c write(6,*)iei,e(iei),nbin,ebin(nbin),ief,e(ief) c STOP 'INITIAL BIN ENERGY OUT OF RANGE...' !ONLY IF SINGLE BATCH RETURN C 10 SUM=DZERO C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) OMGAVG=(OMG(I)+OMG(IPM))/2 !FOR OMEGA MESH .LLT. BIN IF(E(IPM).GT.EBIN(N))THEN !ELSE H0=EBIN(N)-E(I) !STILL .GT. 0 HERE T=H0/(2*H) SUM0=OMG(I)*(DONE-T)+OMG(IPM)*T SUM0=SUM0*H0 SUM=SUM+SUM0 SBIN(N)=SBIN(N)+SUM/(EBIN(N)-EBIN(N-1)) SUM=-SUM0 20 N=N+1 IF(N.GT.NBIN)RETURN IF(E(IPM).GT.EBIN(N))THEN !OMEGA MESH COARSER THAN BIN H0=EBIN(N)-EBIN(N-1) EAVG=(EBIN(N)+EBIN(N-1))/2 T=(EAVG-E(I))/H SUM0=OMG(I)*(DONE-T)+OMG(IPM)*T SUM=SUM-SUM0*H0 SBIN(N)=SBIN(N)+SUM0 GO TO 20 ENDIF ENDIF SUM=SUM+OMGAVG*H ENDDO C !COMPLETE ANY FINAL PARTIALLY FILLED BIN IF(N.EQ.NBIN)THEN H=EBIN(N)-EBIN(N-1) H0=EBIN(N)-E(IPM) T=H0/(2*H) SUM0=OMG(IPM)*(DONE+T)-OMG(IPM-IEINCR)*T SUM0=SUM0*H0 SUM=SUM+SUM0 ENDIF C SBIN(N)=SBIN(N)+SUM/(EBIN(N)-EBIN(N-1)) C RETURN C END C C*********************************************************************** C SUBROUTINE OMGINT(EMSH,OMEGIN,MXE0,MXE1,ETHRSH,TINF,OMGINF X ,EJ,OMEGA,IEI,IEF) C C----------------------------------------------------------------------- C C INTERPOLATE PWB/DW COLLISION STRENGTHS ONTO CONVOLUTION ENERGY MESH. C PWB (ONLY) WILL BE EXTRAPOLATED FOR NEAR DEGENERATE LEVELS. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) C PARAMETER (NLAG=2) !ORDER OF LAGRANGE C DIMENSION EMSH(0:*),OMEGIN(0:*),EJ(0:*),OMEGA(0:*) C ME=MXE0 MEX=MXE1-NLAG+1 NLAGH=(NLAG+1)/2+MXE0 c C DO IE=IEI,IEF C IF(EJ(IE).GT.EMSH(MXE1))GO TO 10 !FOR PWB C 5 IF(EJ(IE).GT.EMSH(NLAGH).AND.ME.LT.MEX)THEN ME=ME+1 NLAGH=NLAGH+1 go to 5 !case X-mesh de~0 ENDIF C CALL POLINT(EMSH(ME),OMEGIN(ME),NLAG,EJ(IE),OMEGA(IE),DELY) C ENDDO C RETURN C C ASSUME BORN/DIPOLE LIMITS EXIST FOR PWB, AND INTERPOLATE ACCORDINGLY. C SO NTYPE=3 NOT USED PRESENTLY. C 10 CONTINUE C UASYME=EMSH(MXE1)/ETHRSH UASYM=UASYME+DONE OMGASYM=OMEGIN(MXE1) C IF(OMGINF.LT.DZERO)THEN !DIPOLE NTYPE=1 EE=EXP(DONE) TX=LOG(UASYME+EE) XASYM=DONE-DONE/TX CASYM=OMGASYM/LOG(UASYME+EE) CINF=ABS(OMGINF)/LOG(TINF) !CINF=4*SLIN/3 ELSEIF(OMGINF.GT.DZERO)THEN !BORN NTYPE=2 TX=UASYM+DONE XASYM=UASYM/TX CASYM=OMGASYM CINF=OMGINF ELSE NTYPE=3 CASYM=OMGASYM*UASYME**2 ENDIF C DO IX=IE,IEF C UE=EJ(IX)/ETHRSH U=UE+DONE C IF(NTYPE.EQ.1)THEN !DIPOLE XRL=LOG(UE+EE) XR=DONE-DONE/XRL OMG=CASYM+CINF*(XR-XASYM)*XRL OMG=OMG*TX ELSEIF(NTYPE.EQ.2)THEN !BORN U1=U+DONE XR=U/U1 OMG=CASYM/U1+CINF*(XR-XASYM) OMG=OMG*TX ELSEIF(NTYPE.EQ.3)THEN !FORBIDDEN OMG=CASYM/UE**2 ENDIF C OMEGA(IX)=OMG C ENDDO C RETURN C END C C*********************************************************************** C SUBROUTINE POLINT(XA,YA,N,X,Y,DY) C C----------------------------------------------------------------------- C C NEVILLE'S ALGORITHM FOR LAGRANGE INTERPOLATION C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (NMAX=10) C PARAMETER(DZERO=0.0D0) C DIMENSION XA(N),YA(N) DIMENSION C(NMAX),D(NMAX) C IF(N.GT.NMAX)STOP 'SR.POLINT: TOO MANY INTERPOLATION POINTS' C NS=1 DIF=ABS(X-XA(1)) DO I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) ENDDO C Y=YA(NS) NS=NS-1 C DO M=1,N-1 DO I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.DZERO)STOP 'SR.POLINT: NUMERICAL FAILURE' DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN ENDDO IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY ENDDO C RETURN END C C*********************************************************************** C SUBROUTINE REMAPX(XMSH,OMEGA,MXX) C C----------------------------------------------------------------------- C C ON THE FIRST CALL, REMAPS THRESHOLD UNIT PARAMETER X'=X+A/(B+X), C WHERE A=3, B=1 CORRESPONDS TO COWAN'S REMAPPING. C ON SUBSEQUENT CALLS, REMAPS THE CORRESPONDING PWB COLLISION STRENGTHS. C C----------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) C PARAMETER (A=3.0D0) PARAMETER (B=1.0D0) PARAMETER (XTHRSH=1.01D0) C DIMENSION XMSH(0:*),OMEGA(0:*) C SAVE IW,W1,W2,ITHRSH,MDROP C DATA ICALL/0/ C C REMAP X-VALUES. C IF(ICALL.EQ.0)THEN C IW=-1 MXX0=MXX MXX=0 C DO M=1,MXX0 XP=XMSH(M) X=(XP-B)**2+4*(B*XP-A) IF(X.GE.DZERO)THEN X=(SQRT(X)+XP-B)/2 IF(X.GE.DONE)THEN IF(MXX.EQ.0)THEN IW=M-1 X1=XMSH(IW) X2=XMSH(IW+1) ENDIF MXX=MXX+1 XMSH(MXX)=X ENDIF ENDIF ENDDO C C DETERMINE INTERPOLATE WEIGHTS FOR X=1 (IF NECESSARY). C IF(XMSH(1).GT.XTHRSH)THEN ITHRSH=1 DO M=MXX,1,-1 XMSH(M+1)=XMSH(M) ENDDO XMSH(1)=DONE XP=DONE+A/(B+DONE) W1=(X2-XP)/(X2-X1) W2=(XP-X1)/(X2-X1) ELSE ITHRSH=0 ENDIF C MDROP=MXX0-MXX MXX=MXX+ITHRSH !ADD SPACE FOR THRESHOLD C c write(*,*)'reducing number of x-values from ',mxx0,' to ', mxx c write(*,*) (m,xmsh(m),m=1,mxx) C ICALL=1 RETURN C ENDIF C C REMAP OMEGAS, DISCARDING MDROP ORIGINAL X-VALUES AND C INSERTING THRESHOLD, IF NECESSARY C IF(ITHRSH.GT.0)OMEGA(1)=W1*OMEGA(IW)+W2*OMEGA(IW+1) c write(*,*)i,omega(iw),omega(1),omega(iw+1) C DO M=1,MXX-ITHRSH OMEGA(M+ITHRSH)=OMEGA(M+MDROP) ENDDO C RETURN END C C*********************************************************************** C INTEGER FUNCTION NOPEN(E,ENAT,NAST,IONE,I0) C C----------------------------------------------------------------------- C C DETERMINES THE NUMBER OF OPEN TARGET STATES AT A GIVEN INCIDENT ENERGY C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION ENAT(NAST) C DO I=I0,NAST IF(E.LT.ENAT(I))GO TO 1 ENDDO I=NAST+1 1 I0=I-1 NOPEN=(I0*(I0-2*IONE+1))/2 C RETURN END C C*********************************************************************** C SUBROUTINE MAX1(E,OMG,IEI,IEF,IEINCR,TEMP,TOLH,NDELBT0,DOWNM1) C C----------------------------------------------------------------------- C C *** MAXWELLIAN DISTRIBUTION. TEMP IS THE ELECTRON TEMPERATURE. C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTHS WHICH LIE BETWEEN ENERGIES E(IEI) AND E(IEF), C WHERE E IS THE SCATTERED ENERGY, I.E. E=0 AT THRESHOLD. C C IF THESE ARE BINNED COLLISION STRENGTHS, I.E. CONSTANT ACROSS THE C BIN WIDTH, THEN USE EXACT ANALYTIC DISTRIBUTION FUNCTION INTEGRAL. C (TRAPEZOIDAL IS ALSO PRESENT FOR TEST PURPOSES - ONLY THE INTEGRAND C IS AVERAGED THEN.) C C OTHERWISE, USE A 2-POINT QUADRATURE IN WHICH THE COLLISION STRENGTH C IS TAKEN TO BE PIECEWISE LINEAR. N.B. DO NOT FORGET THIS POINT! C C THE FIRST NDELBT STEPS FROM THRESHOLD USE THE BURGESS & TULLY METHOD C (ASTRON. ASTROPHYS. 254 436-452 (1992). THIS TREATS THE DISTRIBUTION C FUNCTION EXACTLY, BUT IT IS RELATIVELY SLOW. C C THE REMAINING STEPS (AT HIGHER ENERGIES) REVERT TO THE USUAL FASTER C TRAPEZOIDAL METHOD. THIS TREATS THE DISTRIBUTION FUNCTION AS C PIECEWISE LINEAR AS WELL. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) C PARAMETER (YPAR=20) !MAX E/KT C DIMENSION E(0:*),OMG(0:*) C NDELBT=NDELBT0 IF(NDELBT0.LT.0)NDELBT=NDELBT+1 !ADJUST FLAG C SUM=DZERO C IMN=IEI C C----------------------------------------------------------------------- C C QUADRATURE FOR BINNED DATA - COLLISION STRENGTH IS CONSTANT OVER BIN C WIDTH (IEINCR NORMALLY FORCED TO 1 HERE.) C C----------------------------------------------------------------------- C IF(NDELBT0.LT.0)THEN C C START-OFF WITH EXACT INTEGRATION OF DISTRIBUTION C IMX=IMN-(NDELBT+1)*IEINCR IMX=MIN(IMX,IEF) C Y2=E(IMN-IEINCR)/TEMP !NORMALLY ZERO Y2=EXP(-Y2) !NORMALLY ONE C DO I=IMN,IMX,IEINCR Y1=Y2 Y2=E(I)/TEMP IF(Y2.LT.YPAR)THEN Y2=EXP(-Y2) SUM=SUM+OMG(I)*(Y1-Y2) ELSE IMX=IEF GO TO 5 ENDIF ENDDO C 5 IF(IMX.EQ.IEF)THEN !NORMAL EXIT DOWNM1=SUM RETURN ENDIF C C COMPLETE WITH TRAPEZOIDAL RULE (TEST PURPOSES C.F. OTHER DISTRIBS) C SUM=SUM*(2*TEMP) C IMN=IMX IMX=IEF-IEINCR C Y2=E(IMN)/TEMP Y2=EXP(-Y2) C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) Y1=Y2 Y2=E(IPM)/TEMP IF(Y2.LT.YPAR)THEN Y2=EXP(-Y2) SUM=SUM+OMG(IPM)*(Y1+Y2)*H !AVERAGE DISTRIB ONLY ELSE GO TO 10 ENDIF ENDDO C 10 DOWNM1=SUM/(2*TEMP) RETURN C ENDIF C C----------------------------------------------------------------------- C C QUADRATURE FOR RAW COLLISION STRENGTH (OMEGA NOT BINNED) C C----------------------------------------------------------------------- C H=E(IMN+IEINCR)-E(IMN) IF(ABS(H).GT.TOLH)THEN J=E(IMN)/H IMX=IMN+(NDELBT-J-1)*IEINCR ELSE IMX=IMN+(NDELBT-1)*IEINCR ENDIF IMX=MIN(IMX,IEF-IEINCR) C C START-OFF WITH BURGESS-TULLY INTERPOLATION C Y=E(IMN)/TEMP O2=OMG(IMN)*EXP(-Y) C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) IF(ABS(H).GT.TOLH)THEN !CASE REPEATED ENERGIES DEL=H/TEMP EXPDEL=EXP(+DEL) T=(DONE-EXPDEL)/DEL Y=E(IPM)/TEMP IF(Y.LT.YPAR)THEN O1=O2 O2=OMG(IPM)*EXP(-Y) SUM=SUM+O1*(DONE+T/EXPDEL)-O2*(DONE+T) ELSE DOWNM1=SUM RETURN ENDIF ENDIF ENDDO C SUM=SUM*(2*TEMP) C C COMPLETE WITH TRAPEZOIDAL RULE (FASTER) C IMN=MAX(IMN,IMX+IEINCR) IMX=IEF-IEINCR C Y=E(IMN)/TEMP O2=OMG(IMN)*EXP(-Y) C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) Y=E(IPM)/TEMP IF(Y.LT.YPAR)THEN O1=O2 O2=OMG(IPM)*EXP(-Y) SUM=SUM+(O1+O2)*H !AVERAGE WHOLE INTEGRAND ELSE GO TO 20 ENDIF ENDDO C 20 DOWNM1=SUM/(2*TEMP) C RETURN END C C*********************************************************************** C SUBROUTINE MAX2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMP,DOWNM2) C C----------------------------------------------------------------------- C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTH WITH ENERGIES GREATER THAN (FINAL) EJI. C IF INFINITE ENERGY POINT EXISTS (DIPOLE MUST EXIST AND, IDEALLY, BORN) C THEN INTERPOLATE IN BURGESS-TULLY SPACE, ELSE EXTRAPOLATE (FORBIDDEN), C USING A 2-POINT FORMULA. THE SECOND (LOWER) ENERGY INTERPOLATION POINT C IS GIVEN BY UASYM WITH SCALED OMEGA BY CASYM. THIS (UASYM) NORMALLY C CORRESPONDS TO THE ENERGY EJI, BUT DOES NOT HAVE TO. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) C PARAMETER (YPAR0=15) !MAX E/KT PARAMETER (H0=0.2D0) !LOG STEP C EE=EXP(DONE) C YPAR=YPAR0 IF(NTYPE.EQ.1)YPAR=3*YPAR/2 IF(NTYPE.EQ.3)YPAR=2*YPAR/3 C DOWNM2=DZERO C EJF=YPAR*TEMP !MAX FINAL E C IF(EJF.LE.EJI)RETURN !QUICK RETURN C H=H0 !LOG(EJ) STEP N=LOG(EJF/EJI)/H !NO. STEPS NSIMP=2*(N/2)+1 !ODD C IF(NSIMP.LT.3)RETURN !QUICK RETURN C C SET BURGESS X-VALUE FOR INTERPOLATION (OTHER VALUE IS X=1) C IF(CINF.NE.DZERO)THEN IF(NTYPE.EQ.1)THEN !DIPOLE TX=LOG(UASYM-DONE+EE) XASYM=DONE-DONE/TX ELSEIF(NTYPE.EQ.2)THEN !INTERPOLATE BORN TX=UASYM+DONE XASYM=UASYM/TX ELSE STOP 'TRANSITION TYPE NOT FOUND FOR CINF .NE. ZERO' ENDIF C TX=DONE/(DONE-XASYM) ENDIF C C SIMPSON'S RULE QUADRATURE C X1=LOG(EJI) X=X1-H IFAC=1 C SUM=DZERO C DO N=1,NSIMP C X=X+H EJ=EXP(X) UE=EJ/ETHRSH U=UE+DONE C C INTERPOLATE C IF(NTYPE.EQ.1)THEN !DIPOLE XRL=LOG(UE+EE) XR=DONE-DONE/XRL OMG=CASYM+CINF*(XR-XASYM)*XRL OMG=OMG*TX ELSEIF(NTYPE.EQ.2)THEN IF(CINF.EQ.DZERO)THEN !CONSTANT (NO BORN) EXTRAP OMG=CASYM ELSE !INTERPOLATE BORN U1=U+DONE XR=U/U1 OMG=CASYM/U1+CINF*(XR-XASYM) OMG=OMG*TX ENDIF ELSEIF(NTYPE.EQ.3)THEN !FORBIDDEN (WE HAD BORN) EXTRAP OMG=CASYM/UE**ALF !ALF APPROX =2 ENDIF C Y=EJ/TEMP SUM=SUM+IFAC*EXP(-Y)*EJ*OMG !CONVERTING FROM LOG(EJ) TO EJ C IF(IFAC.LT.3)THEN IFAC=4 ELSE IFAC=2 ENDIF C ENDDO C SUM=SUM-EXP(-Y)*EJ*OMG !CORRECT LAST WEIGHT C DOWNM2=SUM*H/(3*TEMP) !FINALIZE C RETURN C END C C*********************************************************************** C SUBROUTINE KAP1(E,ETHRSH,OMG,IEI,IEF,IEINCR,TEMPE,TOLH,NDELBT0 X ,KAPPA,DOWNK1,UPK1) C C----------------------------------------------------------------------- C C *** KAPPA DISTRIBUTION VERSION OF MAX1. TEMPE IS THE EFFECTIVE TEMP. C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTHS WHICH LIE BETWEEN ENERGIES E(IEI) AND E(IEF), C WHERE E IS THE SCATTERED ENERGY, I.E. E=0 AT THRESHOLD. C C IF THESE ARE BINNED COLLISION STRENGTHS, I.E. CONSTANT ACROSS THE C BIN WIDTH, THEN USE EXACT ANALYTIC DISTRIBUTION FUNCTION INTEGRAL. C (TRAPEZOIDAL IS ALSO PRESENT FOR TEST PURPOSES - ONLY THE INTEGRAND C IS AVERAGED THEN.) C C OTHERWISE, USE A 2-POINT QUADRATURE IN WHICH THE COLLISION STRENGTH C IS TAKEN TO BE PIECEWISE LINEAR. N.B. DO NOT FORGET THIS POINT! C C THE FIRST NDELBT STEPS FROM THRESHOLD USE THE BURGESS & TULLY METHOD C (ASTRON. ASTROPHYS. 254 436-452 (1992). THIS TREATS THE DISTRIBUTION C FUNCTION EXACTLY, BUT IT IS RELATIVELY SLOW. C C THE REMAINING STEPS (AT HIGHER ENERGIES) REVERT TO THE USUAL FASTER C TRAPEZOIDAL METHOD. THIS TREATS THE DISTRIBUTION FUNCTION AS C LINEAR AS WELL. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C REAL*8 KAPPA,KAPPAP,KAPPAM C LOGICAL BLOG C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (D0PT5=0.5D0) PARAMETER (D1PT5=1.5D0) C PARAMETER (ARGFCT=57.0D0) PARAMETER (EXPFCT=650.0D0) PARAMETER (YPAR=1.D16) !MAX (1+E/KT)**KAPPA C DIMENSION E(0:*),OMG(0:*) C NDELBT=NDELBT0 IF(NDELBT0.LT.0)NDELBT=NDELBT+1 !ADJUST FLAG C TEMPK=TEMPE*(KAPPA-D1PT5) !=KAPPA*CHARACTERISTIC ENERGY KAPPAP=KAPPA+DONE KAPPAM=KAPPA-DONE IF(KAPPAP.LT.ARGFCT)THEN FRACK=GAMA(KAPPAP)/(GAMA(KAPPA-D0PT5)*SQRT(KAPPA-D1PT5)) ELSE FRACK=KAPPA*SQRT(KAPPA/(KAPPA-D1PT5)) ENDIF C TL=(E(IEF)+ETHRSH)/TEMPK !WORST CASE SCENARIO TL=KAPPAP*LOG(DONE+TL) T=ETHRSH/TEMPE BLOG=T.GT.EXPFCT.OR.TL.GT.EXPFCT C IF(BLOG)THEN EEE=DZERO !CANNOT REPRESENT EEE=T !OR PROCEED WITH LOGS ELSE EEE=EXP(T) ENDIF C SUMD=DZERO SUMU=DZERO C IMN=IEI C C----------------------------------------------------------------------- C C QUADRATURE FOR BINNED DATA - COLLISION STRENGTH IS CONSTANT OVER BIN C WIDTH (IEINCR NORMALLY FORCED TO 1 HERE.) C C----------------------------------------------------------------------- C IF(NDELBT0.LT.0)THEN C C START-OFF WITH EXACT INTEGRATION OF DISTRIBUTION C IMX=IMN-(NDELBT+1)*IEINCR IMX=MIN(IMX,IEF) C Y2D=E(IMN-IEINCR)/TEMPK !NORMALLY ZERO Y2D=(DONE+Y2D)**KAPPA IF(Y2D.LT.YPAR)THEN Y2D=DONE/(Y2D*KAPPA) IF(EEE.NE.DZERO)THEN Y2U=(E(IMN-IEINCR)+ETHRSH)/TEMPK IF(BLOG)THEN TL=KAPPA*LOG(DONE+Y2U) T=-TL+EEE IF(T.LT.EXPFCT)THEN Y2U=EXP(T)/KAPPA ELSE Y2U=DZERO ENDIF ELSE Y2U=(DONE+Y2U)**KAPPA Y2U=DONE/(Y2U*KAPPA) ENDIF ENDIF ENDIF C DO I=IMN,IMX,IEINCR Y1D=Y2D Y2D=E(I)/TEMPK Y2D=(DONE+Y2D)**KAPPA IF(Y2D.LT.YPAR)THEN Y2D=DONE/(Y2D*KAPPA) SUMD=SUMD+OMG(I)*(Y1D-Y2D) IF(EEE.NE.DZERO)THEN Y1U=Y2U Y2U=(E(I)+ETHRSH)/TEMPK IF(BLOG)THEN TL=KAPPA*LOG(DONE+Y2U) T=-TL+EEE IF(T.LT.EXPFCT)THEN Y2U=EXP(T)/KAPPA ELSE Y2U=DZERO Y1U=DZERO ENDIF ELSE Y2U=(DONE+Y2U)**KAPPA Y2U=DONE/(Y2U*KAPPA) ENDIF SUMU=SUMU+OMG(I)*(Y1U-Y2U) ENDIF ELSE !WE ARE DONE - CONVERGED IMX=IEF GO TO 5 ENDIF ENDDO C 5 IF(IMX.EQ.IEF)THEN !NORMAL EXIT DOWNK1=FRACK*SUMD IF(EEE.NE.DZERO)THEN UPK1=FRACK*SUMU IF(.NOT.BLOG)UPK1=UPK1*EEE ELSE UPK1=DZERO ENDIF RETURN ENDIF C C COMPLETE WITH TRAPEZOIDAL RULE (TEST PURPOSES C.F. OTHER DISTRIBS) C SUMD=SUMD*(2*TEMPK) SUMU=SUMU*(2*TEMPK) C IMN=IMX IMX=IEF-IEINCR C Y2D=E(IMN)/TEMPK Y2D=(DONE+Y2D)**KAPPAP IF(Y2D.LT.YPAR)THEN Y2D=DONE/Y2D IF(EEE.NE.DZERO)THEN Y2U=(E(IMN)+ETHRSH)/TEMPK IF(BLOG)THEN TL=KAPPAP*LOG(DONE+Y2U) T=-TL+EEE IF(T.LT.EXPFCT)THEN Y2U=EXP(T) ELSE Y2U=DZERO ENDIF ELSE Y2U=(DONE+Y2U)**KAPPAP Y2U=DONE/Y2U ENDIF ENDIF ENDIF C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) Y1D=Y2D Y2D=E(IPM)/TEMPK Y2D=(DONE+Y2D)**KAPPAP IF(Y2D.LT.YPAR)THEN Y2D=DONE/Y2D SUMD=SUMD+OMG(IPM)*(Y1D+Y2D)*H !AVERAGE DISTRIB ONLY IF(EEE.NE.DZERO)THEN Y1U=Y2U Y2U=(E(IPM)+ETHRSH)/TEMPK IF(BLOG)THEN TL=KAPPAP*LOG(DONE+Y2U) T=-TL+EEE IF(T.LT.EXPFCT)THEN Y2U=EXP(T) ELSE Y2U=DZERO Y1U=DZERO ENDIF ELSE Y2U=(DONE+Y2U)**KAPPAP Y2U=DONE/Y2U ENDIF SUMU=SUMU+OMG(IPM)*(Y1U+Y2U)*H !AVERAGE DISTRIB ONLY ENDIF ELSE !WE ARE DONE - CONVERGED GO TO 10 ENDIF ENDDO C 10 DOWNK1=FRACK*SUMD/(2*TEMPK) IF(EEE.NE.DZERO)THEN UPK1=FRACK*SUMU/(2*TEMPK) IF(.NOT.BLOG)UPK1=UPK1*EEE ELSE UPK1=DZERO ENDIF RETURN C ENDIF C C----------------------------------------------------------------------- C C QUADRATURE FOR RAW COLLISION STRENGTH (OMEGA NOT BINNED) C C----------------------------------------------------------------------- C H=E(IMN+IEINCR)-E(IMN) IF(ABS(H).GT.TOLH)THEN J=E(IMN)/H IMX=IMN+(NDELBT-J-1)*IEINCR ELSE IMX=IMN+(NDELBT-1)*IEINCR ENDIF IMX=MIN(IMX,IEF-IEINCR) C C START-OFF WITH BURGESS-TULLY INTERPOLATION C Y2D=E(IMN)/TEMPK U2D=(DONE+Y2D*KAPPA)/KAPPAM Y2D=(DONE+Y2D)**KAPPA IF(Y2D.LT.YPAR)THEN Y2D=DONE/(Y2D*KAPPA) IF(EEE.NE.DZERO)THEN Y2U=(E(IMN)+ETHRSH)/TEMPK U2U=(DONE+Y2U*KAPPA)/KAPPAM IF(BLOG)THEN TL=KAPPA*LOG(DONE+Y2U) T=-TL+EEE IF(T.LT.EXPFCT)THEN Y2U=EXP(T)/KAPPA ELSE Y2U=DZERO ENDIF ELSE Y2U=(DONE+Y2U)**KAPPA Y2U=DONE/(Y2U*KAPPA) ENDIF ENDIF ENDIF C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) IF(ABS(H).GT.TOLH)THEN !CASE REPEATED ENERGIES DEL=H/TEMPK U1D=U2D Y1D=Y2D Y2D=E(IPM)/TEMPK U2D=(DONE+Y2D*KAPPA)/KAPPAM Y2D=(DONE+Y2D)**KAPPA IF(Y2D.LT.YPAR)THEN Y2D=DONE/(Y2D*KAPPA) W1=(OMG(IPM)-OMG(I))/DEL W0=(OMG(I)*E(IPM)-OMG(IPM)*E(I))/H SUMD=SUMD+Y1D*(W0+W1*U1D)-Y2D*(W0+W1*U2D) IF(EEE.NE.DZERO)THEN U1U=U2U Y1U=Y2U Y2U=(E(IPM)+ETHRSH)/TEMPK U2U=(DONE+Y2U*KAPPA)/KAPPAM IF(BLOG)THEN TL=KAPPA*LOG(DONE+Y2U) T=-TL+EEE IF(T.LT.EXPFCT)THEN Y2U=EXP(T)/KAPPA ELSE Y2U=DZERO Y1U=DZERO ENDIF ELSE Y2U=(DONE+Y2U)**KAPPA Y2U=DONE/(Y2U*KAPPA) ENDIF W0=W0+ETHRSH*(OMG(I)-OMG(IPM))/H SUMU=SUMU+Y1U*(W0+W1*U1U)-Y2U*(W0+W1*U2U) ENDIF ELSE !WE ARE DONE - CONVERGED DOWNK1=FRACK*SUMD IF(EEE.NE.DZERO)THEN UPK1=FRACK*SUMU IF(.NOT.BLOG)UPK1=UPK1*EEE ELSE UPK1=DZERO ENDIF RETURN ENDIF ENDIF ENDDO C SUMD=SUMD*(2*TEMPK) SUMU=SUMU*(2*TEMPK) C C COMPLETE WITH TRAPEZOIDAL RULE (FASTER) C IMN=MAX(IEI,IMX+IEINCR) IMX=IEF-IEINCR C YD=E(IMN)/TEMPK YD=(DONE+YD)**KAPPAP IF(YD.LT.YPAR)THEN O2D=OMG(IMN)/YD IF(EEE.NE.DZERO)THEN YU=(E(IMN)+ETHRSH)/TEMPK IF(BLOG)THEN TL=KAPPAP*LOG(DONE+YU) T=-TL+EEE IF(T.LT.EXPFCT)THEN YU=EXP(T) ELSE YU=DZERO ENDIF ELSE YU=(DONE+YU)**KAPPAP YU=DONE/YU ENDIF O2U=OMG(IMN)*YU ENDIF ENDIF C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) YD=E(IPM)/TEMPK YD=(DONE+YD)**KAPPAP IF(YD.LT.YPAR)THEN O1D=O2D O2D=OMG(IPM)/YD SUMD=SUMD+(O1D+O2D)*H !AVERAGE WHOLE INTEGRAND IF(EEE.NE.DZERO)THEN YU=(E(IPM)+ETHRSH)/TEMPK IF(BLOG)THEN TL=KAPPAP*LOG(DONE+YU) T=-TL+EEE IF(T.LT.EXPFCT)THEN YU=EXP(T) ELSE YU=DZERO O2U=DZERO ENDIF ELSE YU=(DONE+YU)**KAPPAP YU=DONE/YU ENDIF O1U=O2U O2U=OMG(IPM)*YU SUMU=SUMU+(O1U+O2U)*H !AVERAGE WHOLE INTEGRAND ENDIF ELSE !WE ARE DONE - CONVERGED GO TO 20 ENDIF ENDDO C 20 DOWNK1=FRACK*SUMD/(2*TEMPK) IF(EEE.NE.DZERO)THEN UPK1=FRACK*SUMU/(2*TEMPK) IF(.NOT.BLOG)UPK1=UPK1*EEE ELSE UPK1=DZERO ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE KAP2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMPE X ,KAPPA,DOWNK2,UPK2) C C----------------------------------------------------------------------- C C *** KAPPA DISTRIBUTION VERSION OF MAX2. TEMPE IS THE EFFECTIVE TEMP. C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTH WITH ENERGIES GREATER THAN (FINAL) EJI. C IF INFINITE ENERGY POINT EXISTS (DIPOLE MUST EXIST AND, IDEALLY, BORN) C THEN INTERPOLATE IN BURGESS-TULLY SPACE, ELSE EXTRAPOLATE (FORBIDDEN), C USING A 2-POINT FORMULA. THE SECOND (LOWER) ENERGY INTERPOLATION POINT C IS GIVEN BY UASYM WITH SCALED OMEGA BY CASYM. THIS (UASYM) NORMALLY C CORRESPONDS TO THE ENERGY EJI, BUT DOES NOT HAVE TO. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C REAL*8 KAPPA,KAPPAP C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (D0PT5=0.5D0) PARAMETER (D1PT5=1.5D0) PARAMETER (D1P3=1.0D3) C PARAMETER (ARGFCT=57.0D0) PARAMETER (EXPFCT=650.0D0) C PARAMETER (YPAR0=1.D14) !MAX (1+E/KT)**KAPPA PARAMETER (H0=0.2D0) !LOG STEP C EE=EXP(DONE) C TEMPK=TEMPE*(KAPPA-D1PT5) !=KAPPA*CHARACTERISTIC ENERGY KAPPAP=KAPPA+DONE IF(KAPPAP.LT.ARGFCT)THEN FRACK=GAMA(KAPPAP)/(GAMA(KAPPA-D0PT5)*SQRT(KAPPA-D1PT5)) ELSE FRACK=KAPPA*SQRT(KAPPA/(KAPPA-D1PT5)) ENDIF C T=ETHRSH/TEMPE IF(T.LT.EXPFCT)THEN EEE=EXP(T) ELSE EEE=DZERO !CANNOT REPRESENT ENDIF C YPAR=YPAR0 IF(KAPPA.LT.1.501D0)YPAR=YPAR*100 IF(NTYPE.EQ.1)YPAR=100*YPAR IF(NTYPE.EQ.3)YPAR=YPAR/100 C DOWNK2=DZERO UPK2=DZERO C IF(KAPPA.LT.D1P3)THEN EJF=TEMPK*(YPAR**(DONE/KAPPAP)-DONE) !MAX FINAL E ELSE EJF=TEMPE*LOG(YPAR) ENDIF C IF(EJF.LE.EJI)RETURN !QUICK RETURN C H=H0 !LOG(EJ) STEP N=LOG(EJF/EJI)/H !NO. STEPS NSIMP=2*(N/2)+1 !ODD C IF(NSIMP.LT.3)RETURN !QUICK RETURN C C SET BURGESS X-VALUE FOR INTERPOLATION (OTHER VALUE IS X=1) C IF(CINF.NE.DZERO)THEN IF(NTYPE.EQ.1)THEN !DIPOLE TX=LOG(UASYM-DONE+EE) XASYM=DONE-DONE/TX ELSEIF(NTYPE.EQ.2)THEN !INTERPOLATE BORN TX=UASYM+DONE XASYM=UASYM/TX ELSE STOP 'TRANSITION TYPE NOT FOUND FOR CINF .NE. ZERO' ENDIF C TX=DONE/(DONE-XASYM) ENDIF C C SIMPSON'S RULE QUADRATURE C X1=LOG(EJI) X=X1-H IFAC=1 C SUMD=DZERO SUMU=DZERO C DO N=1,NSIMP C X=X+H EJ=EXP(X) UE=EJ/ETHRSH U=UE+DONE C C INTERPOLATE C IF(NTYPE.EQ.1)THEN !DIPOLE XRL=LOG(UE+EE) XR=DONE-DONE/XRL OMG=CASYM+CINF*(XR-XASYM)*XRL OMG=OMG*TX ELSEIF(NTYPE.EQ.2)THEN IF(CINF.EQ.DZERO)THEN !CONSTANT (NO BORN) EXTRAP OMG=CASYM ELSE !INTERPOLATE BORN U1=U+DONE XR=U/U1 OMG=CASYM/U1+CINF*(XR-XASYM) OMG=OMG*TX ENDIF ELSEIF(NTYPE.EQ.3)THEN !FORBIDDEN (WE HAD BORN) EXTRAP OMG=CASYM/UE**ALF !ALF APPROX =2 ENDIF C T=IFAC*EJ*OMG !CONVERTING FROM LOG(EJ) TO EJ YD=EJ/TEMPK YD=(DONE+YD)**KAPPAP SUMD=SUMD+T/YD C IF(EEE.NE.DZERO)THEN YU=(EJ+ETHRSH)/TEMPK YU=(DONE+YU)**KAPPAP SUMU=SUMU+T/YU ENDIF C IF(IFAC.LT.3)THEN IFAC=4 ELSE IFAC=2 ENDIF C ENDDO C SUMD=SUMD-EJ*OMG/YD !CORRECT LAST WEIGHT IF(EEE.NE.DZERO)SUMU=SUMU-EJ*OMG/YU !CORRECT LAST WEIGHT C DOWNK2=FRACK*SUMD*H/(3*TEMPK) !FINALIZE UPK2=FRACK*SUMU*EEE*H/(3*TEMPK) !FINALIZE C RETURN C END C C*********************************************************************** C SUBROUTINE DRY1(E,ETHRSH,OMG,IEI,IEF,IEINCR,TEMPE,TOLH,NDELBT0 X ,XDRY,DOWNX1,UPX1) C C----------------------------------------------------------------------- C C *** DRUYVESTEYN DISTRBN VERSION OF MAX1. TEMPE IS THE EFFECTIVE TEMP. C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTHS WHICH LIE BETWEEN ENERGIES E(IEI) AND E(IEF), C WHERE E IS THE SCATTERED ENERGY, I.E. E=0 AT THRESHOLD. C C IF THESE ARE BINNED COLLISION STRENGTHS, I.E. CONSTANT ACROSS THE C BIN WIDTH, THEN USE "EXACT" ANALYTIC DISTRIBUTION FUNCTION INTEGRAL. C THIS (INCOMPLETE GAMMA FUNCTION EVALUATION) IS SLOW, SO COMPLETE C WITH TRAPEZOIDAL RULE - ONLY THE INTEGRAND IS AVERAGED THEN. C C OTHERWISE, USE A 2-POINT QUADRATURE IN WHICH THE COLLISION STRENGTH C IS TAKEN TO BE PIECEWISE LINEAR. N.B. DO NOT FORGET THIS POINT! C C THE FIRST NDELBT STEPS FROM THRESHOLD USE THE BURGESS & TULLY METHOD C (ASTRON. ASTROPHYS. 254 436-452 (1992). THIS TREATS THE DISTRIBUTION C FUNCTION EXACTLY, BUT IT IS RELATIVELY SLOW. C C THE REMAINING STEPS (AT HIGHER ENERGIES) REVERT TO THE USUAL FASTER C TRAPEZOIDAL METHOD. THIS TREATS THE DISTRIBUTION FUNCTION AS C LINEAR AS WELL. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C LOGICAL BLOG C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DTWO=2.0D0) PARAMETER (DSIX=6.0D0) PARAMETER (D1PT5=1.5D0) PARAMETER (D2PT5=2.5D0) C PARAMETER (EXPFCT=650.0D0) PARAMETER (YPARL=80.0D0) C PARAMETER (YPAR=20) !MAX (E/KT)**X C DIMENSION E(0:*),OMG(0:*) C NDELBT=NDELBT0 IF(NDELBT0.LT.0)NDELBT=NDELBT+1 !ADJUST FLAG C PI=ACOS(-DONE) C=SQRT(PI/DSIX) T1X=DONE/XDRY T2X=DTWO/XDRY G1X=GAMA(T1X) G2X=GAMA(T2X) T5=D2PT5/XDRY T3=D1PT5/XDRY T5=GAMA(T5) T3=GAMA(T3) T53=T5/T3 C TEMPX=3*TEMPE/(2*T53) !=CHARACTERISTIC ENERGY/T53 FRACX=C*XDRY*SQRT(T53)/T3 C TL=(E(IEF)+ETHRSH)/TEMPX !WORST CASE SCENARIO TL=XDRY*LOG(TL) T=ETHRSH/TEMPE BLOG=T.GT.EXPFCT.OR.TL.GT.LOG(EXPFCT) C IF(BLOG)THEN C T=DZERO !CANNOT REPRESENT EEE=T !OR PROCEED WITH LOGS EEL=T ELSE EEE=EXP(T) EEL=DZERO ENDIF C SUMD=DZERO SUMU=DZERO C IMN=IEI C C----------------------------------------------------------------------- C C QUADRATURE FOR BINNED DATA - COLLISION STRENGTH IS CONSTANT OVER BIN C WIDTH (IEINCR NORMALLY FORCED TO 1 HERE.) C C----------------------------------------------------------------------- C IF(NDELBT0.LT.0)THEN !RESTRICT TO -NDELBT BINS AS GAMIN SLOW C C START-OFF WITH "EXACT" INTEGRATION OF DISTRIBUTION C IMX=IMN-(NDELBT+1)*IEINCR IMX=MIN(IMX,IEF) IF(IMX.LT.IMN)GO TO 10 C IF2D=1 IF2U=1 C Y2D=E(IMN-IEINCR)/TEMPX !NORMALLY ZERO IF(Y2D.NE.DZERO)THEN IF(XDRY*LOG(Y2D).LT.YPARL)THEN Y2D=Y2D**XDRY ELSE Y2D=2*YPAR ENDIF ENDIF IF(Y2D.LT.YPAR)THEN Y2D=GAMIND(IF2D,T1X,Y2D,DZERO) IF(EEE.NE.DZERO)THEN Y2U=(E(IMN-IEINCR)+ETHRSH)/TEMPX IF(BLOG)THEN TL=XDRY*LOG(Y2U) IF(TL.GT.YPARL)Y2U=DZERO ENDIF Y2U=Y2U**XDRY Y2U=GAMIND(IF2U,T1X,Y2U,EEL) ENDIF ELSE DOWNX1=DZERO UPX1=DZERO RETURN ENDIF C c write(6,*)'***',tempe,y2d,y2u DO I=IMN,IMX,IEINCR IF1D=IF2D Y1D=Y2D Y2D=E(I)/TEMPX IF(XDRY*LOG(Y2D).LT.YPARL)THEN Y2D=Y2D**XDRY ELSE Y2D=2*YPAR ENDIF IF(Y2D.LT.YPAR)THEN Y2D=GAMIND(IF2D,T1X,Y2D,DZERO) IF(IF1D*IF2D.LT.0)Y1D=Y1D-G1X SUMD=SUMD+OMG(I)*(Y2D-Y1D) IF(EEE.NE.DZERO)THEN IF1U=IF2U Y1U=Y2U Y2U=(E(I)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(Y2U).GT.YPARL)THEN Y1U=DZERO Y2U=DZERO GO TO 5 ENDIF ENDIF Y2U=Y2U**XDRY Y2U=GAMIND(IF2U,T1X,Y2U,EEL) IF(IF1U*IF2U.LT.0)Y1U=Y1U-G1X*EXP(EEL) 5 SUMU=SUMU+OMG(I)*(Y2U-Y1U) ENDIF c write(6,*)i,y2d,sumd,y2u,sumu ELSE !WE ARE DONE - CONVERGED IMX=IEF GO TO 10 ENDIF ENDDO C 10 IF(IMX.EQ.IEF)THEN !QUICK EXIT DOWNX1=FRACX*SUMD/XDRY IF(EEE.NE.DZERO)THEN UPX1=FRACX*SUMU/XDRY IF(.NOT.BLOG)UPX1=UPX1*EEE ELSE UPX1=DZERO ENDIF RETURN ENDIF C C COMPLETE WITH TRAPEZOIDAL (FASTER) C SUMD=SUMD*(2*TEMPX)/XDRY SUMU=SUMU*(2*TEMPX)/XDRY C IMN=IMX IMX=IEF-IEINCR C Y2D=E(IMN)/TEMPX IF(Y2D.NE.DZERO)THEN IF(XDRY*LOG(Y2D).LT.YPARL)THEN Y2D=Y2D**XDRY ELSE Y2D=2*YPAR ENDIF ENDIF IF(Y2D.LT.YPAR)THEN Y2D=EXP(-Y2D) IF(EEE.NE.DZERO)THEN Y2U=(E(IMN)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(Y2U).LT.YPARL)THEN Y2U=Y2U**XDRY T=-Y2U+EEL IF(T.LT.EXPFCT)THEN Y2U=EXP(T) ELSE Y2U=DZERO ENDIF ELSE Y2U=DZERO ENDIF ELSE Y2U=Y2U**XDRY Y2U=EXP(-Y2U) ENDIF ENDIF ELSE GO TO 20 ENDIF C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) Y1D=Y2D Y2D=E(IPM)/TEMPX IF(XDRY*LOG(Y2D).LT.YPARL)THEN Y2D=Y2D**XDRY ELSE Y2D=2*YPAR ENDIF IF(Y2D.LT.YPAR)THEN Y2D=EXP(-Y2D) SUMD=SUMD+OMG(IPM)*(Y1D+Y2D)*H !AVERAGE DISTRIB ONLY IF(EEE.NE.DZERO)THEN Y1U=Y2U Y2U=(E(IPM)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(Y2U).LT.YPARL)THEN Y2U=Y2U**XDRY T=-Y2U+EEL IF(T.LT.EXPFCT)THEN Y2U=EXP(T) ELSE Y2U=DZERO Y1U=DZERO ENDIF ELSE Y2U=DZERO Y1U=DZERO ENDIF ELSE Y2U=Y2U**XDRY Y2U=EXP(-Y2U) ENDIF SUMU=SUMU+OMG(IPM)*(Y1U+Y2U)*H !AVERAGE DISTRIB ONLY ENDIF ELSE !WE ARE DONE - CONVERGED GO TO 20 ENDIF ENDDO C 20 DOWNX1=FRACX*SUMD/(2*TEMPX) IF(EEE.NE.DZERO)THEN UPX1=FRACX*SUMU/(2*TEMPX) IF(.NOT.BLOG)UPX1=UPX1*EEE ELSE UPX1=DZERO ENDIF RETURN C ENDIF C C----------------------------------------------------------------------- C C QUADRATURE FOR RAW COLLISION STRENGTH (OMEGA NOT BINNED) C C----------------------------------------------------------------------- C H=E(IMN+IEINCR)-E(IMN) IF(ABS(H).GT.TOLH)THEN J=E(IMN)/H IMX=IMN+(NDELBT-J-1)*IEINCR ELSE IMX=IMN+(NDELBT-1)*IEINCR ENDIF IMX=MIN(IMX,IEF-IEINCR) C IF2D=1 IF2U=1 JF2D=1 JF2U=1 C C START-OFF WITH BURGESS-TULLY INTERPOLATION C Y2D=E(IMN)/TEMPX IF(Y2D.NE.DZERO)THEN IF(XDRY*LOG(Y2D).LT.YPARL)THEN Y2D=Y2D**XDRY ELSE Y2D=2*YPAR ENDIF ENDIF IF(Y2D.LT.YPAR)THEN U2D=GAMIND(JF2D,T2X,Y2D,DZERO) Y2D=GAMIND(IF2D,T1X,Y2D,DZERO) IF(EEE.NE.DZERO)THEN Y2U=(E(IMN)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(Y2U).GT.YPARL)THEN Y2U=DZERO U2U=DZERO ENDIF ENDIF Y2U=Y2U**XDRY U2U=GAMIND(JF2U,T2X,Y2U,EEL) Y2U=GAMIND(IF2U,T1X,Y2U,EEL) ENDIF ELSE DOWNX1=DZERO UPX1=DZERO RETURN ENDIF C c write(6,*)'***',tempe,y2d,u2d,y2u,u2u DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) IF(ABS(H).GT.TOLH)THEN !CASE REPEATED ENERGIES DEL=H/TEMPX IF1D=IF2D JF1D=JF2D U1D=U2D Y1D=Y2D Y2D=E(IPM)/TEMPX Y2D=E(IPM)/TEMPX IF(XDRY*LOG(Y2D).LT.YPARL)THEN Y2D=Y2D**XDRY ELSE Y2D=2*YPAR ENDIF IF(Y2D.LT.YPAR)THEN U2D=GAMIND(JF2D,T2X,Y2D,DZERO) IF(JF1D*JF2D.LT.0)U1D=U1D-G2X Y2D=GAMIND(IF2D,T1X,Y2D,DZERO) IF(IF1D*IF2D.LT.0)Y1D=Y1D-G1X W1=(OMG(IPM)-OMG(I))/DEL W0=(OMG(I)*E(IPM)-OMG(IPM)*E(I))/H SUMD=SUMD+W0*(Y2D-Y1D)+W1*(U2D-U1D) IF(EEE.NE.DZERO)THEN IF1U=IF2U JF1U=JF2U U1U=U2U Y1U=Y2U Y2U=(E(IPM)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(Y2U).GT.YPARL)THEN Y1U=DZERO Y2U=DZERO U1U=DZERO U2U=DZERO GO TO 25 ENDIF ENDIF Y2U=Y2U**XDRY U2U=GAMIND(JF2U,T2X,Y2U,EEL) IF(JF1U*JF2U.LT.0)U1U=U1U-G2X*EXP(EEL) Y2U=GAMIND(IF2U,T1X,Y2U,EEL) IF(IF1U*IF2U.LT.0)Y1U=Y1U-G1X*EXP(EEL) W0=W0+ETHRSH*(OMG(I)-OMG(IPM))/H 25 SUMU=SUMU+W0*(Y2U-Y1U)+W1*(U2U-U1U) ENDIF c write(6,*)y2d,u2d,sumd,y2u,u2u,sumu ELSE !WE ARE DONE - CONVERGED DOWNX1=FRACX*SUMD/XDRY IF(EEE.NE.DZERO)THEN UPX1=FRACX*SUMU/XDRY IF(.NOT.BLOG)UPX1=UPX1*EEE ELSE UPX1=DZERO ENDIF RETURN ENDIF ENDIF ENDDO C SUMD=SUMD*(2*TEMPX)/XDRY SUMU=SUMU*(2*TEMPX)/XDRY C C COMPLETE WITH TRAPEZOIDAL RULE (FASTER) C IMN=MAX(IMN,IMX+IEINCR) IMX=IEF-IEINCR C YD=E(IMN)/TEMPX IF(YD.NE.DZERO)THEN IF(XDRY*LOG(YD).LT.YPARL)THEN YD=YD**XDRY ELSE YD=2*YPAR ENDIF ENDIF IF(YD.LT.YPAR)THEN O2D=OMG(IMN)*EXP(-YD) IF(EEE.NE.DZERO)THEN YU=(E(IMN)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(YU).LT.YPARL)THEN YU=YU**XDRY T=-YU+EEL IF(T.LT.EXPFCT)THEN YU=EXP(T) ELSE YU=DZERO ENDIF ELSE YU=DZERO ENDIF ELSE YU=YU**XDRY YU=EXP(-YU) ENDIF O2U=OMG(IMN)*YU ENDIF ELSE GO TO 30 ENDIF C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) YD=E(IPM)/TEMPX IF(XDRY*LOG(YD).LT.YPARL)THEN YD=YD**XDRY ELSE YD=2*YPAR ENDIF IF(YD.LT.YPAR)THEN O1D=O2D O2D=OMG(IPM)*EXP(-YD) SUMD=SUMD+(O1D+O2D)*H !AVERAGE WHOLE INTEGRAND IF(EEE.NE.DZERO)THEN YU=(E(IPM)+ETHRSH)/TEMPX IF(BLOG)THEN IF(XDRY*LOG(YU).LT.YPARL)THEN YU=YU**XDRY T=-YU+EEL IF(T.LT.EXPFCT)THEN YU=EXP(T) ELSE YU=DZERO O2U=DZERO ENDIF ELSE YU=DZERO O2U=DZERO ENDIF ELSE YU=YU**XDRY YU=EXP(-YU) ENDIF O1U=O2U O2U=OMG(IPM)*YU SUMU=SUMU+(O1U+O2U)*H !AVERAGE WHOLE INTEGRAND ENDIF ELSE !WE ARE DONE - CONVERGED GO TO 30 ENDIF ENDDO C 30 DOWNX1=FRACX*SUMD/(2*TEMPX) IF(EEE.NE.DZERO)THEN UPX1=FRACX*SUMU/(2*TEMPX) IF(.NOT.BLOG)UPX1=UPX1*EEE ELSE UPX1=DZERO ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE DRY2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMPE X ,XDRY,DOWNX2,UPX2) C C----------------------------------------------------------------------- C C *** DRUYVESTEYN DISTRBN VERSION OF MAX2. TEMPE IS THE EFFECTIVE TEMP. C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTH WITH ENERGIES GREATER THAN (FINAL) EJI. C IF INFINITE ENERGY POINT EXISTS (DIPOLE MUST EXIST AND, IDEALLY, BORN) C THEN INTERPOLATE IN BURGESS-TULLY SPACE, ELSE EXTRAPOLATE (FORBIDDEN), C USING A 2-POINT FORMULA. THE SECOND (LOWER) ENERGY INTERPOLATION POINT C IS GIVEN BY UASYM WITH SCALED OMEGA BY CASYM. THIS (UASYM) NORMALLY C CORRESPONDS TO THE ENERGY EJI, BUT DOES NOT HAVE TO. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DSIX=6.0D0) PARAMETER (D1PT5=1.5D0) PARAMETER (D2PT5=2.5D0) C PARAMETER (EXPFCT=650.0D0) C PARAMETER (YPAR0=15) !MAX (E/KT)**X PARAMETER (H0=0.2D0) !LOG STEP C EE=EXP(DONE) C PI=ACOS(-DONE) C=SQRT(PI/DSIX) T5=D2PT5/XDRY T3=D1PT5/XDRY T5=GAMA(T5) T3=GAMA(T3) T53=T5/T3 C TEMPX=3*TEMPE/(2*T53) !=CHARACTERISTIC ENERGY/T53 FRACX=C*XDRY*SQRT(T53)/T3 C T=ETHRSH/TEMPE IF(T.LT.EXPFCT)THEN EEE=EXP(T) ELSE EEE=DZERO !CANNOT REPRESENT ENDIF C YPAR=YPAR0 IF(NTYPE.EQ.1)YPAR=3*YPAR/2 IF(NTYPE.EQ.3)YPAR=2*YPAR/3 IF(XDRY.LT.DONE)YPAR=YPAR*(DONE-LOG(XDRY)) C DOWNX2=DZERO UPX2=DZERO C EJF=TEMPX*YPAR**(DONE/XDRY) !MAX FINAL E C IF(EJF.LE.EJI)RETURN !QUICK RETURN C H=H0 !LOG(EJ) STEP IF(XDRY.GT.DONE)H=H/(DONE+LOG(XDRY)) N=LOG(EJF/EJI)/H !NO. STEPS NSIMP=2*(N/2)+1 !ODD C IF(NSIMP.LT.3)RETURN !QUICK RETURN C C SET BURGESS X-VALUE FOR INTERPOLATION (OTHER VALUE IS X=1) C IF(CINF.NE.DZERO)THEN IF(NTYPE.EQ.1)THEN !DIPOLE TX=LOG(UASYM-DONE+EE) XASYM=DONE-DONE/TX ELSEIF(NTYPE.EQ.2)THEN !INTERPOLATE BORN TX=UASYM+DONE XASYM=UASYM/TX ELSE STOP 'TRANSITION TYPE NOT FOUND FOR CINF .NE. ZERO' ENDIF C TX=DONE/(DONE-XASYM) ENDIF C C SIMPSON'S RULE QUADRATURE C X1=LOG(EJI) X=X1-H IFAC=1 C SUMD=DZERO SUMU=DZERO C DO N=1,NSIMP C X=X+H EJ=EXP(X) UE=EJ/ETHRSH U=UE+DONE C C INTERPOLATE C IF(NTYPE.EQ.1)THEN !DIPOLE XRL=LOG(UE+EE) XR=DONE-DONE/XRL OMG=CASYM+CINF*(XR-XASYM)*XRL OMG=OMG*TX ELSEIF(NTYPE.EQ.2)THEN IF(CINF.EQ.DZERO)THEN !CONSTANT (NO BORN) EXTRAP OMG=CASYM ELSE !INTERPOLATE BORN U1=U+DONE XR=U/U1 OMG=CASYM/U1+CINF*(XR-XASYM) OMG=OMG*TX ENDIF ELSEIF(NTYPE.EQ.3)THEN !FORBIDDEN (WE HAD BORN) EXTRAP OMG=CASYM/UE**ALF !ALF APPROX =2 ENDIF C T=IFAC*EJ*OMG !CONVERTING FROM LOG(EJ) TO EJ YD=EJ/TEMPX YD=YD**XDRY SUMD=SUMD+T*EXP(-YD) C IF(EEE.NE.DZERO)THEN YU=(EJ+ETHRSH)/TEMPX YU=YU**XDRY SUMU=SUMU+T*EXP(-YU) ENDIF C IF(IFAC.LT.3)THEN IFAC=4 ELSE IFAC=2 ENDIF C ENDDO C SUMD=SUMD-EJ*OMG*EXP(-YD) !CORRECT LAST WEIGHT IF(EEE.NE.DZERO)SUMU=SUMU-EJ*OMG*EXP(-YU) !CORRECT LAST WEIGHT C DOWNX2=FRACX*SUMD*H/(3*TEMPX) !FINALIZE UPX2=FRACX*SUMU*EEE*H/(3*TEMPX) !FINALIZE C RETURN C END C C*********************************************************************** C SUBROUTINE NUM1(E,ETHRSH,OMG,IEI,IEF,IEINCR,TEMPE,TOLH,NDELBT X ,NENG37,E37,F37,DOWNN1,UPN1) C C----------------------------------------------------------------------- C C *** NUMERICAL DISTRIBUTION (adf37). TEMPE IS THE ELECTRON TEMPERATURE. C C *** CURRENTLY DOES NOT INTEGRATE OUTSIDE OF THE NUMERICAL DISTRIBUTION C TABULATED ENERGY RANGE DEFINED BY E37(1) TO E37(NENG37)!!! C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTHS WHICH LIE BETWEEN ENERGIES E(IEI) AND E(IEF), C WHERE E IS THE SCATTERED ENERGY, I.E. E=0 AT THRESHOLD. C C CURRENTLY, TRAPEZOIDAL RULE IS USED FOR BOTH BINNED (NDELBT.LT.0) AND C NON-BINNED (NDELBT.GT.0) COLLISION STRENGTHS, AND ACROSS THE ENTIRE C ENERGY RANGE. (N.B. TOLH IS NOT USED.) C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) C PARAMETER (EXPFCT=650.0D0) PARAMETER (T2ORTPI=1.1283792D0) !2/sqrt(pi) C DIMENSION E(0:*),OMG(0:*),E37(*),F37(*) C DATA ICON/3/ !0 CAN TEST EXACT MAXELLIAN DATA XDUM/0/ !NOT USED C SUMD=DZERO SUMU=DZERO C T=ETHRSH/TEMPE IF(T.LT.EXPFCT)THEN EEE=EXP(T) ELSE EEE=DZERO !CANNOT REPRESENT UPN1=DZERO ENDIF C IMN=IEI IMX=IEF-IEINCR C C----------------------------------------------------------------------- C C QUADRATURE FOR BINNED DATA - COLLISION STRENGTH IS CONSTANT OVER BIN C WIDTH (IEINCR NORMALLY FORCED TO 1 HERE.) C C----------------------------------------------------------------------- C IF(NDELBT.LT.0)THEN C ED=E(IMN) Y2D=FBAR(ED,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) IF(EEE.NE.DZERO)THEN EU=ED+ETHRSH Y2U=FBAR(EU,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) ENDIF C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) Y1D=Y2D ED=E(IPM) Y2D=FBAR(ED,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) SUMD=SUMD+OMG(IPM)*(Y1D+Y2D)*H !AVERAGE DISTRIB ONLY IF(EEE.NE.DZERO)THEN Y1U=Y2U EU=ED+ETHRSH Y2U=FBAR(EU,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) SUMU=SUMU+OMG(IPM)*(Y1U+Y2U)*H !AVERAGE DISTRIB ONLY ENDIF ENDDO C DOWNN1=SUMD*SQRT(TEMPE)/(2*T2ORTPI) !AS FABR~T2ORTPI/TEMP**3/2 UPN1=DZERO IF(EEE.NE.DZERO) X UPN1=SUMU*EEE*SQRT(TEMPE)/(2*T2ORTPI) !AS FABR~T2ORTPI/TEMP**3/2 C RETURN C ENDIF C C----------------------------------------------------------------------- C C QUADRATURE FOR RAW COLLISION STRENGTH (OMEGA NOT BINNED) C C----------------------------------------------------------------------- C ED=E(IMN) O2D=OMG(IMN)*FBAR(ED,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) EU=ED+ETHRSH O2U=OMG(IMN)*FBAR(EU,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) C DO I=IMN,IMX,IEINCR IPM=I+IEINCR H=E(IPM)-E(I) O1D=O2D ED=E(IPM) O2D=OMG(IPM)*FBAR(ED,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) SUMD=SUMD+(O1D+O2D)*H !AVERAGE WHOLE INTEGRAND IF(EEE.NE.DZERO)THEN O1U=O2U EU=ED+ETHRSH O2U=OMG(IPM)*FBAR(EU,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37) SUMU=SUMU+(O1U+O2U)*H !AVERAGE WHOLE INTEGRAND ENDIF ENDDO C DOWNN1=SUMD*SQRT(TEMPE)/(2*T2ORTPI) !AS FABR~T2ORTPI/TEMP**3/2 IF(EEE.NE.DZERO) XUPN1=SUMU*EEE*SQRT(TEMPE)/(2*T2ORTPI) !AS FABR~T2ORTPI/TEMP**3/2 C RETURN END C C*********************************************************************** C SUBROUTINE NUM2(EJI,ETHRSH,NTYPE,UASYM,CASYM,CINF,ALF,TEMPE X ,NENG37,E37,F37,DOWNN2,UPN2) C C----------------------------------------------------------------------- C C *** NUMERICAL DISTRIBUTION (adf37). TEMPE IS THE ELECTRON TEMPERATURE. C C *** CURRENTLY DOES NOT INTEGRATE OUTSIDE OF THE NUMERICAL DISTRIBUTION C TABULATED ENERGY RANGE DEFINED BY E37(1) TO E37(NENG37)!!! C C CALCULATE THE CONTRIBUTION TO THE EFFECTIVE COLLISION STRENGTH FROM C THE COLLISION STRENGTH WITH ENERGIES GREATER THAN (FINAL) EJI. C IF INFINITE ENERGY POINT EXISTS (DIPOLE MUST EXIST AND, IDEALLY, BORN) C THEN INTERPOLATE IN BURGESS-TULLY SPACE, ELSE EXTRAPOLATE (FORBIDDEN), C USING A 2-POINT FORMULA. THE SECOND (LOWER) ENERGY INTERPOLATION POINT C IS GIVEN BY UASYM WITH SCALED OMEGA BY CASYM. THIS (UASYM) NORMALLY C CORRESPONDS TO THE ENERGY EJI, BUT DOES NOT HAVE TO. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) C PARAMETER (EXPFCT=650.0D0) PARAMETER (T2ORTPI=1.1283792D0) !2/sqrt(pi) C PARAMETER (H0=0.2D0) !LOG STEP C DIMENSION E37(*),F37(*) C DATA ICON/3/ !0 CAN TEST EXACT MAXELLIAN DATA XDUM/0/ !NOT USED C EE=EXP(DONE) C T=ETHRSH/TEMPE IF(T.LT.EXPFCT)THEN EEE=EXP(T) ELSE EEE=DZERO !CANNOT REPRESENT ENDIF C DOWNN2=DZERO UPN2=DZERO C C************************************************* C EJF=E37(NENG37) !<<<<<<<<<<<<<<<<3/2? GO TO 1 endif IF(EEE.LE.EOLD.OR.TEMPE.LT.TOLD)NOLD=1 EOLD=EEE DO N=NOLD,NENG37 IF(EEE.LT.E37(N))THEN N2=N+NLAGH-1 N1=N-NLAGH IF(N1.LE.0)THEN N2=NLAG N1=1 ELSEIF(N2.GT.NENG37)THEN N2=NENG37 N1=N2-NLAG+1 ENDIF NOLD=N GO TO 1 ENDIF ENDDO NOLD=NENG37 N2=NENG37 N1=N2-NLAG+1 1 CONTINUE DO N=N1,N2 DD0=DONE DO M=N1,N2 IF(N.NE.M)THEN DD0=DD0*(EEE-E37(M)) DD0=DD0/(E37(N)-E37(M)) ENDIF ENDDO c dd0=dd0*sqrt(e37(n)) !interp full f(e) e.g. SOL QUP=QUP+DD0*F37(N) ENDDO c if(eee.ne.dzero)qup=qup/sqrt(eee) !if interp full f QUP=QUP/T2ORTPI ENDIF C FBAR=T2ORTPI*QUP !*SQRT(EEE) FOR FULL F C RETURN END C C*********************************************************************** C SUBROUTINE ADF37(IPRINT,TEMPE37,ICON,NENG37,E37,F37) C C CHECK DISTRIBUTION NORMALIZATION AND TRY AND SET EFFECTIVE TEMPS C TO HELP LATER RR QUADRATURE (SEE SR.BODE). MORE COSMETIC FOR DR. C NORMALLY CALLED WITH ICON=3, AS ANALYTIC SHOULD BE KNOWN/EXACT. C BUT MAY HAVE ALTERNATE NUMERICAL ICON FLAGGED IN THE FUTURE. C C KNOWN ISSUES FOR NUMERICAL DISTRIBUTIONS: C C CONTRIBUTION FROM THE TAILS I.E.OUTSIDE THE RANGE OF ENERGY TABULATION C CURRENTLY, NO EXTRAPOLATION IS CARRIED-OUT - A MAXWELLIAN IS COMMENTED C OUT FOR THE ADF37 TEST CASE TO ILLUSTRATE ITS EFFECT/IMPORTANCE. C C INTERPOLATION OF THE DISTRIBUTION F(E) OR F(E)/SQRT(E) ETC (WHERE F(E) C IS NORMED TO UNITY). THE FORMER IS NECESSARY FOR THE ADF37 SOL EXAMPLE C (NORM IS WILDLY OFF OTHERWISE) BUT THE LATTER IS MUCH MORE ACCURATE C AT HIGH-E FOR THE MAXWELLIAN CASE. WOULD EXPECT THE LATTER TO BE C PREFERABLE. THE SOL EXAMPLE HAS VERY FEW POINTS AT LOW E, SO SUSPECT. C C SEE FN.FBAR FOR CURRENT EXTRAPOLATION & INTERPOLATION AND DEVELOPMENT. C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (NDIM1=1001) PARAMETER (NDIM37=19) C PARAMETER (IZERO=0) PARAMETER (DONE=1.0D0) C PARAMETER (TNORM=1.0D-3) !RENORM F(E) IF ERROR EXCEEDS TNORM PARAMETER (TMEAN=0.1D-0) !RECALC E*F(E) IF ERROR EXCEEDS TMEAN PARAMETER (ITMAX=5) !FOR ITERATION ON TEMPE=E*F(E) c c parameter (conryk=1.5789d5) !1.578885d5) !rydbergs to kelvin C DIMENSION E37(*),F37(*) C DATA XDUM/0/ !NOT USED C FMAX=-1 NFMAX=0 DO N=1,NENG37 F=F37(N)*SQRT(E37(N)) !AS F37 IS FBAR NOW IF(F.GT.FMAX)THEN FMAX=F NFMAX=N ENDIF ENDDO TEMPE37=2*E37(NFMAX) !CASE MAXWELLIAN C NSTEP=4001 C ITER=0 2 TEMPE=TEMPE37 ITER=ITER+1 C IF(IPRINT.GE.1)WRITE(0,"(/'T_EFF=',1P,D12.5)")TEMPE IF(IPRINT.GE.0)WRITE(6,"(/'T_EFF=',1P,D12.5)")TEMPE C EMAX=1000*TEMPE EMIN=TEMPE/1000 C EE=EMIN EMAX=LOG(EMAX) EMIN=LOG(EMIN) DE=(EMAX-EMIN)/(NSTEP-1) C SUMFF=0 SUMFM=0 SUMFFE=0 SUMFME=0 EE=0 !MAX(0,EE-EXP(DE)) FF=0 FM=0 E=EMIN-DE DO N=1,NSTEP E0=EE FF0=FF FM0=FM E=E+DE EE=EXP(E) TSQRTE=SQRT(EE) FF=FBAR(EE,TEMPE,ICON,XDUM,XDUM,NENG37,E37,F37)*TSQRTE FM=FBAR(EE,TEMPE,IZERO,XDUM,XDUM,NENG37,E37,F37)*TSQRTE IF(IPRINT.GE.2) X WRITE(77,'(I5,1P,3D12.3)')N,EE,FM/TSQRTE,FF/TSQRTE H=EE-E0 SUMFF=SUMFF+H*(FF+FF0) SUMFM=SUMFM+H*(FM+FM0) SUMFFE=SUMFFE+H*(FF*EE+FF0*E0) SUMFME=SUMFME+H*(FM*EE+FM0*E0) ENDDO C SUMFF=SUMFF/2 SUMFM=SUMFM/2 C IF(IPRINT.GE.1) XWRITE(0,"('MAX NORM=',1P,D12.5,3X,'DIST NORM=',D12.5)")SUMFM,SUMFF IF(IPRINT.GE.0) XWRITE(6,"('MAX NORM=',1P,D12.5,3X,'DIST NORM=',D12.5)")SUMFM,SUMFF C SUMFFE=2*SUMFFE/6 SUMFME=2*SUMFME/6 C IF(IPRINT.GE.1) XWRITE(0,"('MAXWELL TE =',1P,D12.5,3X,'DISTRIBN TE =',D12.5)") XSUMFME,SUMFFE !*conryk IF(IPRINT.GE.0) XWRITE(6,"('MAXWELL TE =',1P,D12.5,3X,'DISTRIBN TE =',D12.5)") XSUMFME,SUMFFE !*conryk C TEMPE37=SUMFFE C RAT=SUMFFE/TEMPE IF((ABS(RAT-DONE).GT.TMEAN.OR.ITER.EQ.1).AND.ITER.LT.ITMAX)GO TO 2 C IF(ABS(SUMFF-DONE).GT.TNORM)THEN IF(IPRINT.GE.0)WRITE(0,1000)J,SUMFF IF(IPRINT.GE.-1)WRITE(6,1000)J,SUMFF DO N=1,NENG37 F37(N)=F37(N)/SUMFF ENDDO ENDIF C RETURN C 1000 FORMAT(/'*** ATTENTION: RENORMALIZING ADF37 DISTRIBUTION FOR ' X ,'TEMPERATURE INDEX',I4,' BY FACTOR',F6.3) C END C C*********************************************************************** C REAL*8 FUNCTION GAMA(X) C C----------------------------------------------------------------------- C C FN.GAMA EVALUATES THE GAMMA FUNCTION OF ARGUMENT X, C USING LANCZOS' FORMULA. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (NMAX=6) PARAMETER (ARGFCT=120.0D0) !=D.P., 57. FOR S.P. C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DHALF=0.5D0) C DIMENSION COEF(6) C DATA COEF,STP,SER0 /76.18009172947146D0,-86.50532032941677D0, & 24.01409824083091D0, -1.231739572450155D0, & .1208650973866179D-2, -.5395239384953D-5, & 2.5066282746310005D0, 1.000000000190015D0/ C IF(X.GT.ARGFCT)THEN WRITE(6,*)'***ERROR IN FN.GAMA: ARGUMENT TOO LARGE (OVERFLOW) =' X ,X STOP '***ERROR IN FN.GAMA: ARGUMENT TOO LARGE (OVERFLOW)' ENDIF C IF(X.LE.DZERO)THEN WRITE(6,*)'***ERROR IN FN.GAMA: ARGUMENT IS NOT POSITIVE =',X STOP '***ERROR IN FN.GAMA: ARGUMENT IS NOT POSITIVE' ENDIF C G=NMAX-1 TE=X+DHALF TMP=G+TE TMP=EXP(-TMP)*TMP**TE C SER=SER0 Y=X C DO J=1,NMAX Y=Y+DONE SER=SER+COEF(J)/Y ENDDO C GAMA=TMP*STP*SER/X C END C C*********************************************************************** C REAL*8 FUNCTION GAMIND(IFLAG,A,X,E) C C----------------------------------------------------------------------- C C EVALUATES THE FOLLOWING UNREGULARIZED LOWER INCOMPLETE GAMMA FUNCTION: C C INT_0^X [EXP(-X+E)*X**(A-1)] DX C C I.E., WITH OPTIONAL ADDITIONAL CONSTANT FACTOR EXP(E) WHICH MAY HELP C KEEP THE EXPONENTIAL REPRESENTABLE. C C *** THIS VERSION IS FOR USE WHEN EVALUATING *DIFFERENCES* OF FUNCTIONS C AS SUCH, WHEN X .GE. A+1 IT SETS GAMIND=-GAMIND, ONLY. C IT DOES *NOT* ADD GAMMA(A) TO FORM THE FULL COMPLEMENT. C WHEN EVALUATING DIFFERENCES, THIS GAMMA(A) CANCELS EXACTLY C (IMPLICITLY) AND SO USE OF THIS ROUTINE AVOIDS CANCELLATION ERROR. C THE USER MUST CHECK THAT SUCCESSIVE VALUES (WHICH FORM THE C DIFFERENCE) HAVE THE SAME IFLAG SIGN. IF NOT, GAMMA(A) MUST BE C ADDED TO THE NEGATIVE IFLAG ONE OR SUBTRACTED FROM THE POSITIVE. C CANCELLATION SHOULD NOT BE AN ISSUE HERE AT X .EQ. A+1. C THE USE OF IFLAG AVOIDS UNDERFLOW ISSUES, FLOATING POINT PRECISION C AND ALLOWS THE INTERNAL SWITCH BETWEEN METHODS TO BE CHANGED C WITHOUT NEEDING TO BE SYNCH'ED WITH THE CALLING PROGRAM. C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DTWO=2.0D0) C PARAMETER (EXPFCT=650.0D0) PARAMETER (FPMIN=1.0D-30) PARAMETER (ITMAX=100) PARAMETER (TOL=3.0D-7) C IF(X.LT.DZERO)THEN WRITE(6,*)' *** ERROR IN FN.GAMIND: X=',X,' .LT. 0' STOP '*** ERROR IN FN.GAMIND: X .LT. 0' ENDIF C IF(A.LE.DZERO)THEN WRITE(6,*)' *** ERROR IN FN.GAMIND: A=',A,' .LE. 0' STOP '*** ERROR IN FN.GAMIND: A .LT. 0' ENDIF C IF(X.EQ.DZERO)THEN GAMIND=DZERO IFLAG=1 RETURN ENDIF C IF(X.LT.A+DONE)THEN C AP=A SUM=DONE/A DEL=SUM C DO I=1,ITMAX AP=AP+DONE DEL=DEL*X/AP SUM=SUM+DEL IF(ABS(DEL).LT.ABS(SUM)*TOL)THEN T=-X+A*LOG(X)+E IF(T.LT.EXPFCT)THEN GAMIND=SUM*EXP(T) ELSE GAMIND=DZERO ENDIF IFLAG=1 RETURN ENDIF ENDDO C ELSE C B=X+DONE-A C=DONE/FPMIN D=DONE/B H=D DO I=1,ITMAX AN=-I*(I-A) B=B+DTWO D=AN*D+B IF(ABS(D).LT.FPMIN)D=FPMIN C=B+AN/C IF(ABS(C).LT.FPMIN)C=FPMIN D=DONE/D DEL=D*C H=H*DEL IF(ABS(DEL-DONE).LT.TOL)THEN T=-X+A*LOG(X)+E IF(T.LT.EXPFCT)THEN GAMIND=-EXP(T)*H ELSE GAMIND=DZERO ENDIF IFLAG=-1 RETURN ENDIF ENDDO C ENDIF C WRITE(6,*)'*** FN.GAMIND: TOLERANCE IS TOO SMALL: INCREASE ITMAX' STOP '*** FN.GAMIND: TOLERANCE IS TOO SMALL: INCREASE ITMAX' C C IFLAG=0 C RETURN !UNCOMMENT IF JUST WARN C END C C*********************************************************************** C C----------------------------------------------------------------------- C C A COLLECTION OF INTEGER ORDER MODIFIED BESSEL FUNCTION ROUTINES. C C----------------------------------------------------------------------- C C*********************************************************************** C REAL*8 FUNCTION BESSK(N,X) IMPLICIT REAL*8(A-H,O-Z) C C----------------------------------------------------------------------- C C MODIFIED BESSEL FUNCTION K(N,X) GENERATED BY UP RECURRENCE RELATION: C C K(N+1,X)=(2N/X)*K(N,X)+K(N-1,X) (A&S: 9.6.26) C C----------------------------------------------------------------------- C IF(N.LT.2)STOP 'BAD ARGUMENT N IN BESSK' C TOX=2.0D0/X BKM=BESSK0(X) BK=BESSK1(X) DO J=1,N-1 BKP=BKM+J*TOX*BK BKM=BK BK=BKP ENDDO BESSK=BK C RETURN END C C*********************************************************************** C REAL*8 FUNCTION BESSK0(X) C C----------------------------------------------------------------------- C C MODIFIED BESSEL FUNCTION K(0,X) GENERATED FROM ABRAMOWITZ & STEGUN C POLYNOMIAL FITS 9.8.5 AND 9.8.6 C C----------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) C DATA P1,P2,P3,P4,P5,P6,P7/-0.57721566D0,0.42278420D0,0.23069756D0, X 0.3488590D-1,0.262698D-2,0.10750D-3,0.74D-5/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,-0.7832358D-1,0.2189568D-1, X -0.1062446D-1,0.587872D-2,-0.251540D-2,0.53208D-3/ C IF (X.LE.2.0D0) THEN Y=X*X/4.0D0 BESSK0=(-LOG(X/2.0D0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y* X (P6+Y*P7)))))) ELSE Y=(2.0D0/X) BESSK0=(EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y* X Q7)))))) ENDIF C RETURN END C C*********************************************************************** C REAL*8 FUNCTION BESSK1(X) IMPLICIT REAL*8(A-H,O-Z) C C----------------------------------------------------------------------- C C MODIFIED BESSEL FUNCTION K(1,X) GENERATED FROM ABRAMOWITZ & STEGUN C POLYNOMIAL FITS 9.8.7 AND 9.8.8 C C----------------------------------------------------------------------- C DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,0.15443144D0,-0.67278579D0, X -0.18156897D0,-0.1919402D-1,-0.110404D-2,-0.4686D-4/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,0.23498619D0,-0.3655620D-1, X 0.1504268D-1,-0.780353D-2,0.325614D-2,-0.68245D-3/ C IF (X.LE.2.0D0) THEN Y=X*X/4.0D0 BESSK1=(LOG(X/2.0D0)*BESSI1(X))+(1.0D0/X)*(P1+Y*(P2+Y*(P3+Y* X (P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=2.0D0/X BESSK1=(EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y* X Q7)))))) ENDIF C RETURN END C C*********************************************************************** C REAL*8 FUNCTION BESSI(N,X) IMPLICIT REAL*8(A-H,O-Z) C C----------------------------------------------------------------------- C C MODIFIED BESSEL FUNCTION I(N,X) GENERATED BY DOWN RECURRENCE RELATION: C C I(N+1,X)=-(2N/X)*I(N,X)+I(N-1,X) (A&S: 9.6.26) C C----------------------------------------------------------------------- C PARAMETER (BIGNO=1.0D10) PARAMETER (BIGNI=1.0D-10) C PARAMETER (IACC=40) C IF (N.LT.2)STOP 'BAD ARGUMENT N IN BESSI' C IF (X.EQ.0.0D0) THEN BESSI=0.0D0 ELSE TOX=2.0D0/ABS(X) BIP=0.0D0 BI=1.0D0 BESSI=0.0D0 M=2*((N+INT(SQRT(DFLOAT(IACC*N))))) DO J=M,1,-1 BIM=BIP+DFLOAT(J)*TOX*BI BIP=BI BI=BIM IF (ABS(BI).GT.BIGNO) THEN BESSI=BESSI*BIGNI BI=BI*BIGNI BIP=BIP*BIGNI ENDIF IF (J.EQ.N) BESSI=BIP ENDDO BESSI=BESSI*BESSI0(X)/BI IF (X.LT.0.0D0.AND.MOD(N,2).EQ.1) BESSI=-BESSI ENDIF C RETURN END C C*********************************************************************** C REAL*8 FUNCTION BESSI0(X) IMPLICIT REAL*8(A-H,O-Z) C C----------------------------------------------------------------------- C C MODIFIED BESSEL FUNCTION I(0,X) GENERATED FROM ABRAMOWITZ & STEGUN C POLYNOMIAL FITS 9.8.1 AND 9.8.2 C C----------------------------------------------------------------------- C DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0, X 1.2067492D0,0.2659732D0,0.360768D-1,0.45813D-2/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, X 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1, X 0.2635537D-1,-0.1647633D-1,0.392377D-2/ C IF (ABS(X).LT.3.75D0)THEN Y=(X/3.75D0)**2 BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))) ELSE AX=ABS(X) Y=3.75D0/AX BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y* X (Q7+Y*(Q8+Y*Q9)))))))) ENDIF C RETURN END C C*********************************************************************** C REAL*8 FUNCTION BESSI1(X) IMPLICIT REAL*8(A-H,O-Z) C C----------------------------------------------------------------------- C C MODIFIED BESSEL FUNCTION I(1,X) GENERATED FROM ABRAMOWITZ & STEGUN C POLYNOMIAL FITS 9.8.3 AND 9.8.4 C C----------------------------------------------------------------------- C DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0, X 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, X -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, X -0.2895312D-1,0.1787654D-1,-0.420059D-2/ C IF (ABS(X).LT.3.75D0) THEN Y=(X/3.75D0)**2 BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE AX=ABS(X) Y=3.75D0/AX BESSI1=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y* X (Q7+Y*(Q8+Y*Q9)))))))) IF(X.LT.0.0D0)BESSI1=-BESSI1 ENDIF C RETURN END