SUBROUTINE smdem2(itype, ear, ne, abun, dens, z, ninputt, inputt, & dem, ifl, qintherm, invelocity, photar, status) IMPLICIT NONE INTEGER itype, ne, ninputt, ifl, status REAL ear(0:*), abun(*), inputt(*), dem(*), photar(*) REAL dens, z, invelocity LOGICAL qintherm c Subroutine to calculate the summed emission from plasma with a c DEM given by the input array. This routine handles the case of c the model being interpolated from the new APEC format input files c with separate lines and continua. Checks for the APECTHERMAL and c APECVELOCITY keywords - note that the results of these could conflict c with values set in higher level routines so need to write out a warning c if this occurs. c Arguments : c itype I i: type of plasma emission file c 1 = R-S c 4 = APEC c ear R i: model energy ranges c ne I i: number of model energies c abun R i: abundances c dens R i: density (cm^-3) c z R i: redshift c ninputt I i: number of temperatures c inputt R i: temperatures c dem R i: emission measures for input temperatures c ifl I i: dataset number (unused at present) c qintherm L i: if true apply thermal broadening (set by higher c level routine) c invelocityR i: gaussian velocity broadening to apply (set by c higher level routine) c photar R r: spectrum c status I r: 0==OK INCLUDE '../../inc/xspec.inc' C Pointers to local arrays C itval Tabulated temperatures C icont Tabulated continua C icenerg Energies for tabulated continua C iconEl Atomic numbers for tabulated elements C iconpt Offsets in the cont array for each Z and kT C iconsiz Sizes in the cont array for each Z and kT C ipseudo Tabulated pseudo-continua C ipenerg Energies for tabulated pseudo-continua C ipsupt Offsets in the pseudo array for each Z and kT C ipsusiz Sizes in the pseudo array for each Z and kT C ilineE Tabulated line energies C ilineEm Tabulated line emissivities C ilineEl Tabulated line elements C ilinpt Offsets in the line arrays for each kT C ilinsiz Sizes in the line arrays for each kT INTEGER itval, icont, iconEl, icenerg, iconpt, iconsiz, ipseudo INTEGER ipenerg, ipsupt, ipsusiz, ilineE, ilineEm, ilineEl INTEGER ilinpt, ilinsiz c Local variables REAL velocity INTEGER nmtval, nmelt, itsave, lenp, ie, ierr LOGICAL qinit, qtherm CHARACTER*72 contxt CHARACTER*255 wrtstr CHARACTER*128 pname, pvalue, ovalue INTEGER lenact CHARACTER*128 fgmstr EXTERNAL fgmstr, lenact SAVE qinit, ovalue, itsave, nmtval, nmelt SAVE itval, icont, iconEl, icenerg, iconpt, iconsiz, ipseudo SAVE ipenerg, ipsupt, ipsusiz, ilineE, ilineEm, ilineEl SAVE ilinpt, ilinsiz DATA qinit, ovalue, itsave /.TRUE., ' ', 0/ C First check for values of APECTHERMAL and APECVELOCITY xset parameters c Check the APECTHERMAL variable - if the variable has not been set then c use the value input to the routine pname = 'APECTHERMAL' pvalue = fgmstr(pname) lenp = lenact(pvalue) CALL upc(pvalue) qtherm = qintherm IF ( lenp .GT. 0 ) THEN IF ( (pvalue(:lenp) .EQ. 'YES' .OR. & pvalue(:lenp) .EQ. 'ON') ) THEN qtherm = .TRUE. ELSEIF ((pvalue(:lenp) .EQ. 'NO' .OR. & pvalue(:lenp) .EQ. 'OFF') ) THEN qtherm = .FALSE. ENDIF ENDIF c and the APECVELOCITY variable - if it doesn't exist then use the value c input. If a non-zero value was input and the variable has been set then c write a warning pname = 'APECVELOCITY' pvalue = fgmstr(pname) lenp = lenact(pvalue) velocity = invelocity IF ( lenp .GT. 0 ) THEN READ(pvalue(:lenp), *, iostat=ierr) velocity IF ( ierr .NE. 0 ) THEN WRITE(wrtstr,*) & 'Failed to read value from APECVELOCITY - assuming ', & invelocity CALL xwrite(wrtstr, 10) velocity = invelocity ELSE IF ( invelocity .GT. 0. .AND. & ABS(velocity-invelocity) .GT. 1.e-5 ) THEN CALL xwrite( &'Warning: inconsistency between velocity broadening set by model', & 10) CALL xwrite( &' and APECVELOCITY variable - using that from model', 10) velocity = invelocity ENDIF ENDIF ENDIF WRITE(wrtstr, *) 'Thermal broadening : ', qtherm CALL xwrite(wrtstr, 20) WRITE(wrtstr, *) 'Velocity broadening : ', velocity CALL xwrite(wrtstr, 20) C Algorithm is as follows. The APED data format has two files. The continuum C file contains an extension for each tabulated temperature. Each extension C has a row for each element with vectors of continuum energies, continuum C emissivities, pseudo-continuum energies, and pseudo-continuum emissivities. C Emissivities are in ph/cm^3/s/keV and energies are not uniformly distributed. C The line file also contains an extension for each tabulated temperature. C Each extension has a row for each line with wavelengths, emissivities C (ph/cm^3/s), and the element from which the line comes. C C The continuum energies and emissivities are read into a single array for C all Z, and kT with pointer and size arrays tracking where the data for C a given (Z,kT) pair is found. The same system is used for the C pseudo-continuum data. The line energies, emissivities, and elements are C also read into a single array, this time for all kT. Pointer and size C arrays track where the lines for a given kT can be found. C C This routine loops over the temperatures requested, for each temperature C interpolating between the nearest temperatures tabulated. For the continuum C and pseudo-continuum it also loops over elements, weighting by the C abundances. The lines are added in with the line emissivity multiplied C by the abundance. C Initialize the output array DO ie = 1, ne photar(ie) = 0.0 ENDDO C Check whether we need to reinitialize. If so this is either the first C time through or we have switched between plasma emission models or we C have switched APEC model files. IF ( itype .NE. itsave ) THEN qinit = .TRUE. itsave = itype ENDIF pname = 'APECROOT' pvalue = fgmstr(pname) IF ( pvalue .NE. ovalue ) THEN qinit = .TRUE. ovalue = pvalue ENDIF C If the first time through then load the model file. If this file has C already been loaded then LDAPED just returns the pointers to the C appropriate memory. IF ( qinit ) THEN CALL ldaped(itype, itval, nmtval, nmelt, icont, icenerg, & iconEl, iconpt, iconsiz, ipseudo, ipenerg, ipsupt, & ipsusiz, ilineE, ilineEm, ilineEl, ilinpt, & ilinsiz, status) contxt = 'Failure in LDAPED' IF ( status .NE. 0 ) GOTO 999 IF ( itype .EQ. 1 .AND. nmelt .NE. 13 ) THEN CALL xwrite('R-S file read does not have 12 elements', 5) RETURN ENDIF IF ( itype .EQ. 4 .AND. nmelt .NE. 14 ) THEN CALL xwrite('APEC file read does not have 13 elements', 5) RETURN ENDIF ENDIF c Now run the routine which actually does all the work CALL sumape(ear, ne, abun, z, ninputt, inputt, dem, nmtval, & MEMR(itval), nmelt, MEMR(icont), MEMR(icenerg), & MEMI(iconEl), MEMI(iconpt), MEMI(iconsiz), & MEMR(ipseudo), MEMR(ipenerg), MEMI(ipsupt), & MEMI(ipsusiz), MEMR(ilineE), MEMR(ilineEm), & MEMI(ilineEl), MEMI(ilinpt), MEMI(ilinsiz), qtherm, & velocity, photar, status) contxt = 'Failure in SUMAPE' IF ( status .NE. 0 ) GOTO 999 qinit = .FALSE. 999 CONTINUE IF ( status .NE. 0 ) THEN WRITE(wrtstr,'(''Error '', i3, '' in SMDEM2'')') status CALL xwrite(wrtstr, 10) CALL xwrite(contxt, 10) ENDIF RETURN END