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******************************************************
