C CODE ADD.FOR C C M.J. SEATON, 28.10.96 C C READS FILES STAGES.zz C CREATES FILES ACC.zz C C THE FILES STAGES.zz HAVE DATA FOR INDIVIDUAL IONISATTION STGAES. C THE FILES ACC.zz HAVE DATA SUMMED OVER IONISATION STAGES. C THE WEIGHTS USED ARE SET IN SUBROUTINE WEIGHTS. C THE DEPENDENCE OF DIFFUSION COEFFICIENTS ON CHARGE IS C SET IN SUBROUTINE ZETA. C C PARAMETERS: C IPN FOR NUMBER OF CHI VALUES C IPZ FOR NUCEAR CHARGE PARAMETER(IPN=10,IPZ=28) C DIMENSION CHI(IPN),GAM(0:IPZ,IPN),GAMP(0:IPZ,IPN), + ROSS(IPN),ROSSP(IPN),PHI(0:IPZ) DIMENSION W(0:IPZ),ZET(0:IPZ) DIMENSION GAMB(IPN),GAMPB(IPN) C CHARACTER*128 STGFL,ACCFL c !!! logical neg !!! neg=.false. !!! C C OPEN FILES C PRINT*,' READ NAME OF INPUT FILE (FILE STAGES.zz)' READ(5,'(A)')STGFL PRINT*,' READ NAME OF OUTPUT FILE (FILE ACC.zz)' READ(5,'(A)')ACCFL C OPEN(1,FILE=STGFL,STATUS='OLD') OPEN(2,FILE=ACCFL,STATUS='NEW') C C IZ=NUCLEAR CHARGE OF SELECTED ELEMENT, k. C NTOT=NUMBER OF FREQUENCY POINTS FOR u IN RANGE 0 TO 20. C AMACC=ATOM MASS FOR k IN amu. C FACC=STANDARD ABUNDANCE OF ELEMENT, k. C CHI IS THE FACTOR BY WHICH ABUNDANCE OF k IS MODIFED. C NCHI=NUMBER OF TABULAR VALUES OF CHI. C THE MEAN ATOM MASS IS FMU=FMU0+CHI*FMU1 IN amu. C TEMPERATURE INDEX, I=I1,I2,I3 (LOG10(T)=0.025*I) C READ(1,*)IZ,NTOT,NCHI,FMU0,FMU1 CALL CHECKD(IPZ,IZ,IPN,NCHI) READ(1,*)AMACC,FACC READ(1,*)(CHI(N),N=1,NCHI) READ(1,*)I1,I2,I3 C WRITE(2,2000)IZ,NTOT,NCHI,FMU0,FMU1 WRITE(2,2001)AMACC,FACC WRITE(2,2002)(CHI(N),N=1,NCHI) WRITE(2,2003)I1,I2,I3 C C START LOOP ON TEMPERATURES C DO II=I1,I2,I3 C C FOR EACH I, DENSITY INDEX IS J=J1,J2,J3 (LOG10(N_e)=0.25*J). C READ(1,*,END=100)I,J1,J2,J3 CALL CHECKI(II,I) C C T=TEMPERATURE T=10**(0.025*I) c WRITE(2,2004)I,J1,J2,J3 C C START LOOP ON DENSITIES. C DO JJ=J1,J2,J3 C C THE NUMNER OF ELECTRONS PER ATOM FOR THE MIXTURE IS EPA=EP0+CHI*EP1. C THE IONISATION STAGES ARE JZ=JZ1 TO JZ2 (JZ=0 FOR NEUTRAL ATOM). C THE IONISATION FRACTIONS ARE PHI(JZ). C ROSS(N)=ROSSELAND-MEAN CROSS SECTION (IN ATOMIC UNITS) FOR C ABUNDANCE MULTIPLIER CHI(N). C ROSSP(N)=DERIVATIVE OF ROSS WITH RESPECT TO CHI. C READ(1,*,END=101)J,EP0,EP1,JZ1,JZ2 CALL CHECKJ(I,JJ,J) READ(1,*)(PHI(JZ),JZ=JZ1,JZ2) READ(1,*)(ROSS(N),N=1,NCHI) READ(1,*)(ROSSP(N),N=1,NCHI) C C FN=ELECTRON DENSITY FN=10**(0.25*J) C C START LOOP ON ION STAGES. C DO JJZ=JZ1,JZ2 C C GAM(JZ,N)= gamma PARAMETER FOR ION STAGE JZ AND ABUNDANCE N. C GAMP(JZ,N)=DERIVATIVE WITH RESSPECT TO CHI. C READ(1,*)JZ,(GAM(JZ,N),N=1,NCHI) CALL CHCKJZ(I,J,JJZ,JZ) READ(1,*)JZ,(GAMP(JZ,N),N=1,NCHI) C ENDDO C CALL ZETA(T,FN,JZ1,JZ2,ZET) CALL ZETBAR(JZ1,JZ2,PHI,ZET,ZETB) CALL WGHTS(JZ1,JZ2,W) CALL GAMBAR(JZ1,JZ2,NCHI,GAM,GAMP,W,GAMB,GAMPB) C WRITE(2,2005)J,EP0,EP1 WRITE(2,2006)ZETB WRITE(2,2006)(ROSS(N),N=1,NCHI) WRITE(2,2006)(ROSSP(N),N=1,NCHI) WRITE(2,2006)(GAMB(N),N=1,NCHI) WRITE(2,2006)(GAMPB(N),N=1,NCHI) c !!!!!! do n=1,nchi if(gamb(n).le.0)then if(.not.neg)then write(6,6003) neg=.true. endif write(6,6004)i,j,n,gamb(n) endif enddo c !!!!!! C ENDDO ENDDO C GOTO 200 C C TERMINATION BEFORE END OF LOOP ON I. 100 WRITE(6,6001)I GOTO 200 C C TERMINATION BEFORE END OF LOOP ON J. 101 WRITE(6,6002)I,J C C TERMINATION C 200 CLOSE(1) CLOSE(2) WRITE(6,6000)STGFL,ACCFL STOP C 2000 FORMAT(I5,I7,I5,1P,2E14.5) 2001 FORMAT(1P,2E12.5) 2002 FORMAT(1P,8E9.1) 2003 FORMAT(3I5) 2004 FORMAT(I5,5X,3I5) 2005 FORMAT(I5,5X,1P,2E12.4) 2006 FORMAT(1P,7E11.3) 6000 FORMAT(/5X,'READ FILE:'/5X,A128/5X,'WRITTEN FILE;'/ + 5X,A128/) 6001 FORMAT(/5X,'DID NOT COMPLETE LOOP ON I, DONE TO', + ' I=',I3/) 6002 FORMAT(/5X,'FOR I=',I3,' DID NOT COMPLETE LOOP ON J', + ', DONE TO J=',I3/) 6003 format(/10X,' NEGATIVE GAM:'// + 4x,'I',4X,'J',4X,'N',9X,'GAM'/) 6004 FORMAT(3I5,1P,E12.3) C END C******************************************************** C SUBROUTINE CHECKD(IPZ,IZ,IPN,NCHI) C IF(IZ.GT.IPZ)THEN WRITE(6,6000)IZ,IPZ STOP ENDIF C IF(NCHI.GT.IPN)THEN WRITE(6,6001)NCHI,IPN STOP ENDIF C 6000 FORMAT(/5X,' IZ=',I3,', IPZ=',I3/ + 5X,'REQUIRE IZ.LE.IPZ. STOPPING'/) 6001 FORMAT(/5X,' NCHI=',I3,', IPN=',I3/ + 5X,'REQUIRE NCHI.LE.IPN. STOPPING'/) C RETURN END C C************************************************** C SUBROUTINE CHECKI(II,I) C IF(II.NE.I)THEN WRITE(6,6000)II,I STOP ENDIF C 6000 FORMAT(/5X,'In loop DO II=I1,I2,I3:'/ + 5X,'For II=',I3,' reads I=',I3/ + 5x,'Data error. STOP') C RETURN END C************************************************* SUBROUTINE CHECKJ(I,JJ,J) C IF(JJ.NE.J)THEN WRITE(6,6000)I,JJ,J STOP ENDIF C 6000 FORMAT(/5X,'In loop DO JJ=J1,J2,J3 for I=',I3/ + 5X,'For JJ=',i3,' reads J=',I3/ + 5x,'Data error. STOP') C RETURN END C*************************************************** SUBROUTINE CHCKJZ(I,J,JJZ,JZ) C IF(JJZ.NE.JZ)THEN WRITE(6,6000)I,J,JJZ,JZ STOP ENDIF C RETURN C 6000 FORMAT(/5X,'In loop D0 JJZ=JZ1, JZ2 for I=',I3, + ' and J=',I3/5x,'For JJZ=',I3,' reads JZ=',I3/ + 5X,'Data error. STOP'/) C END C***************************************************** C SUBROUTINE ZETA(T,FN,JZ1,JZ2,ZET) C C OBTAINS FACTORS zeta GIVING DEPENDENCE OF DIFFUSION COEFFICIENT C ON IONNISATION STAGE. C C FOR POSITIVE IONS (JZ>0) USES FORMULAE FROM ALLER AND CHAPMAN, C ApJ, 132, 461, 1960. C FOR JZ=0 USES APPROXIMATION RECOMMENDED BY GONZALEZ, LEBLANC, C ARTRU AND MICHAUD, A&A, 297, 223, 1995. C C NOTE. THIS SUBROUTINE CAN BE CHANGED IF THE USER PREFERS SOME C OTHER EXPRESSIONS FOR THE COEFFICIENTS. C PARAMETER(IPZ=28) DIMENSION ZET(0:IPZ) C B=2.7285E8*T**3/FN IF(JZ1.EQ.0)ZET(0)=100./MAX(1.,LOG(1.+B)) DO JZ=MAX(1,JZ1),JZ2 ZM2=1./FLOAT(JZ*JZ) ZET(JZ)=ZM2/MAX(1.,LOG(1.+B*ZM2)) ENDDO C RETURN END C C*********************************************** C SUBROUTINE WGHTS(JZ1,JZ2,W) C C CALCULATES THE WEIGHTS BY WHICH THE RADIATIVE ACCELERATIONS FOR C EACH HIONISATION STAGE ARE TO BE MULTIPLIED. C C THIS ROUTINE USES WEIGHTS W(JZ)=1. C IT CAN BE CHANGED IF THE USER PREFERS SOME OTHER WEIGHTS. C PARAMETER(IPZ=28) DIMENSION W(0:IPZ) C DO JZ=JZ1,JZ2 W(JZ)=1. ENDDO C RETURN END C C********************************************** C SUBROUTINE ZETBAR(JZ1,JZ2,PHI,ZET,ZETB) C C CALCULATES ZETB=MEAN VALUE OF zeta, USING COEFFICIENTS ZET C CALCULATED IN SUBROUTINE ZETA. C PARAMETER(IPZ=28) DIMENSION PHI(0:IPZ),ZET(0:IPZ) C ZETB=0. DO JZ=JZ1,JZ2 ZETB=ZETB+PHI(JZ)*ZET(JZ) ENDDO C RETURN END C C*********************************************** C SUBROUTINE GAMBAR(JZ1,JZ2,NCHI,GAM,GAMP,W,GAMB,GAMPB) C C CALCULATES GAMB=MEAN VALUE OF gamma USING WEIGHTS CALCULATED C IN SUBROUTINE WGHTS. C PARAMETER (IPZ=28,IPN=10) DIMENSION GAM(0:IPZ,IPN),GAMP(0:IPZ,IPN),W(0:IPZ), + GAMB(IPN),GAMPB(IPN) C DO N=1,NCHI GAMB(N)=0. GAMPB(N)=0. DO JZ=JZ1,JZ2 GAMB(N)=GAMB(N)+W(JZ)*GAM(JZ,N) GAMPB(N)=GAMPB(N)+W(JZ)*GAMP(JZ,N) ENDDO ENDDO C RETURN END C C******************************************************