C  Code diff.f
C
C  MJS, 3.1.97
C
C  Reads data from file1, created using accfit.f
C
C  Obtains values of:
C     ZET=factor in expression for diffusion coefficient;
C     ROSS=Rosseland-mean opaity in cgs;
C     GRAD=radiative acceleration in cgs.
C  For any specified values of:
C     n=depth-point in stellar model;
C     CHI=abundance multiplier.
C  Output  written to file2.
C
C  Contents  of file1:
C    NS=number of depth-points;
C    MC=number of input values of CHIL=LOG10(CHI);
C    CN=smallest input values  of  CHIL;
C    DC=interval in input CHIL;
C    For depth-points n=1 to NS:
C       FLT(n)=LOG10(T), T  in K;
C       FLR(n)=LOG100(RHO), RHO  in  cgs;
C       For CHIL  points, m=1 to MC:
C           ZET(n,m)=factor ZETA
C           RS(n,m)=Rosseland mean (cgs);
C           GR(n,m)=radiative acceleratios (gs).
C
C  NOTE. T and RHO are not used by this code, but are read and written.
C
C  Parameters:
C      IPS=maximum value of NS;
C      IPC=maximum value of MC.
      PARAMETER(ips=2000,IPC=20)
C
      DIMENSION FLT(IPS),FLR(IPS),ZET(IPS,IPC),RS(IPS,IPC),GR(IPS,IPC)
C
      COMMON/CFIT/AR(IPS,IPC),BR(IPS,IPC),CR(IPS,IPC),DR(IPS,IPC),
     +            AG(IPS,IPC),BG(IPS,IPC),CG(IPS,IPC),DG(IPS,IPC),
     +            AZ(IPS,IPC),BZ(IPS,IPC),CZ(IPS,IPC),DZ(IPS,IPC),
     +            CN,DC,CL,CH,MC
C
C  Open  files
      CHARACTER*64 FILE1, FILE2
      WRITE(6,6000)
      READ(5,5000)FILE1
      OPEN(1,FILE=FILE1,STATUS='OLD',ERR=100)
      WRITE(6,6001)
      READ(5,5000)FILE2
      OPEN(2,FILE=FILE2,STATUS='NEW',ERR=200)
      WRITE(2,2000)
C
C  Read FILE1
      READ(1,*)NS,CN,DC,MC
      CALL CHECK(NS,MC)
      READ(1,*)(NN,FLT(N),FLR(N),
     +(ZET(N,M),M=1,MC),(RS(N,M),M=1,MC),(GR(N,M),M=1,MC),N=1,NS)
      CLOSE(1)
C
C  Get parameters for cubic-spline fits
      CALL FIT(NS,ZET,RS,GR)
C
C  Set limits on CHIL=LOG10(CHI).
C    Input data give CHIL in range CN to  CN+DC*(MC-1).
C    Allow  interpolations for intervals of .5*DC outside of that range.
C  Lower limit
      CL=CN-.5*DC
C  Higher limit
      CH=CN+DC*(MC-.5)
C
C  Make interpolations.
C
C  Read depth-point
   1  WRITE(6,6002)
      READ(5,*)N
      IF(N.LE.0)THEN
	 CLOSE(2)
         STOP
      ENDIF
      IF(N.GT.NS)THEN
         WRITE(6,6003) 
         GOTO 1
      ENDIF
      WRITE(6,6004)N,FLT(N),FLR(N)
C
C  Read required value of CHI
   2  WRITE(6,6005)
      READ(5,*)CHI
	PRINT*,' ****** CHI=',CHI
      IF(CHI.LE.0.)GOTO 1
      CHIL=LOG10(CHI)
      IF(CHIL.LT.CL.OR.CHIL.GT.CH)THEN
         WRITE(6,6006)CL,CH
         GOTO  2
      ENDIF
C
C  Obtain Z=ZETA, R=ROSS, G=GRAD
      Z=ZETA(N,CHI)
      R=ROSS(N,CHI)
      G=GRAD(N,CHI)
C
C  Write results
      WRITE(6,6007)CHI,Z,R,G
      WRITE(2,2001)N,FLT(N),FLR(N),CHI,Z,R,G
C
C  Go to next  CHIL  case
      GOTO 2
C
C  Messages for errors in opening files
  100 WRITE(6,6008)FILE1
      STOP
  200 WRITE(6,6009)FILE2
      STOP
C
 2000 FORMAT(4X,'N',4X,'LOG10(T)',2X,'LOG10(RHO)',9X,'CHI',
     + 8X,'ZETA',8X,'ROSS',8X,'GRAD'/)
 2001 FORMAT(I5,2F12.3,1P,4E12.3)
 5000 FORMAT(A64)
 6000 FORMAT('  Read name of input file, up to 64 characters')
 6001 FORMAT('  Read name of OUTput file, up to 64 characters')
 6002 FORMAT('  Read depth-point N (read N.LE.0 to stop)')
 6003 FORMAT('  Should have  N.LE.NSTAR, NSTAR=',I5)
 6004 FORMAT('  N=',I5,', LOG10(T)=',F6.3,',  LOG10(RHO)=',
     + F9.3/)
 6005 FORMAT('  Read CHIL=CHI (read CHI=0 for',
     + ' next depth-point)')
 6006 FORMAT('  Require LOG10(CHI) in range ',F6.3,' to ',F6.3)
 6007 FORMAT(5x,'CHI=',F7.3,', ZETA=',1P,E10.3,', ROSS=',E10.3,
     + ', GRAD=',E10.3)
 6008 FORMAT('  Error in openning input file'/
     +  5X,A64)
 6009 FORMAT('  Error in openning output file'/
     +  5X,A64)
C
      END
C
C************************
      SUBROUTINE CHECK(NS,MC)
C
      PARAMETER(ips=2000,IPC=20)
C
C  Check dimensions
      IF(NS.GT.IPS)THEN
        WRITE(6,6000)NS
        STOP
      ELSEIF(MC.GT.IPC)THEN
        WRITE(6,6001)MC
        STOP
      ENDIF
C
C  Check MC .GE. 5
C  (Condition MC.GE.5 is required for interpolation routines)
      IF(MC.LT.5)THEN
	 WRITE(6,6002)MC
	 STOP
      ENDIF
C
      RETURN
C
 6000 FORMAT('  Require parameter IPS of at least',I5)
 6001 format('  Require parameter IPC of at least',i5)
 6002 FORMAT('  MC=',I2/',  Interpolation routines require MC .GE. 5')
C
      END
C
C********************************************
	SUBROUTINE FIT(NS,ZET,RS,GR)
C
C  Given
C    RS(N,M)=Rosseland mean
C    GR(N,M)=radiative acceleration
C  at the tabular points (N,M).
C
C  Obtain coefficients in cubic-spline fits required for the calculation
C  of ROSS and GRAD as continuous functions of CHIL=LOG10(CHI).
C
C  The procedure is to fit LOG(RS) and LOG(GR) to cubics in CHIL
C  within each inerval between input tabular points.
C  For all CHIL, the fitted functions are continuous and have 
C  continuous first and second derivatives.
C
      PARAMETER(ips=2000,IPC=20)
C
      DIMENSION ZET(IPS,IPC),RS(IPS,IPC),GR(IPS,IPC),X(IPS),P(IPS)
C
      COMMON/CFIT/AR(IPS,IPC),BR(IPS,IPC),CR(IPS,IPC),DR(IPS,IPC),
     +            AG(IPS,IPC),BG(IPS,IPC),CG(IPS,IPC),DG(IPS,IPC),
     +            AZ(IPS,IPC),BZ(IPS,IPC),CZ(IPS,IPC),DZ(IPS,IPC),
     +            CN,DC,CL,CH,MC

C
        DC2=DC**(-2)
	DC3=DC**(-3)
C
	DO N=1,NS
C
C          PARAMETERS FOR ZETA
	   DO M=1,MC
	      X(M)=LOG(ZET(N,M))
	   ENDDO
	   CALL PRIME(X,P,DC,MC)
	   DO M=1,MC-1
	      AZ(N,M)=X(M)
	      BZ(N,M)=P(M)
	      CZ(N,M)=(3.*(X(M+1)-X(M))-DC*(2.*P(M)+P(M+1)))*DC2
	      DZ(N,M)=(-2.*(X(M+1)-X(M))+DC*(P(M)+P(M+1)))*DC3
	   ENDDO
C
C          PARAMETERS FOR ROSS
	   DO M=1,MC
	      X(M)=LOG(RS(N,M))
	   ENDDO
	   CALL PRIME(X,P,DC,MC)
	   DO M=1,MC-1
	      AR(N,M)=X(M)
	      BR(N,M)=P(M)
	      CR(N,M)=(3.*(X(M+1)-X(M))-DC*(2.*P(M)+P(M+1)))*DC2
	      DR(N,M)=(-2.*(X(M+1)-X(M))+DC*(P(M)+P(M+1)))*DC3
	   ENDDO
C
C          PARAMETERS FOR GRAD
	   DO M=1,MC
	      X(M)=LOG(GR(N,M))
	   ENDDO
	   CALL PRIME(X,P,DC,MC)
	   DO M=1,MC-1
	      AG(N,M)=X(M)
	      BG(N,M)=P(M)
	      CG(N,M)=(3.*(X(M+1)-X(M))-DC*(2.*P(M)+P(M+1)))*DC2
	      DG(N,M)=(-2.*(X(M+1)-X(M))+DC*(P(M)+P(M+1)))*DC3
	   ENDDO
C
	ENDDO
C
	WRITE(6,6000)
C
	RETURN
C
 6000   FORMAT('  Done FIT')
C
	END
C
C***********************************************************************
	SUBROUTINE PRIME(X,P,DEL,N)
C
C  The function X is tabulated at values X(I) corresponding to equal
C  intervals DEL in an independent variable.
C
C  PRIME gives P(I) to be the first derivative of X to be used in 
C  cubic-spline fits.
C
	PARAMETER(IPC=20)
C
	DIMENSION X(IPC),P(IPC),B(IPC),D(IPC),E(IPC)
C
	P(1)=(-1.5*X(1)+2.*X(2)-0.5*X(3))/DEL
	P(N)=(+1.5*X(N)-2.*X(N-1)+0.5*X(N-2))/DEL
C
	DO I=2,N-1
	   B(I)=(X(I+1)-X(I-1))*3./DEL
	ENDDO
C
	D(2)=-.25
	E(2)=.25*(B(2)-P(1))
	DO I=3,N-2
	   D(I)=-1./(D(I-1)+4.)
	   E(I)=-D(I)*(B(I)-E(I-1))
	ENDDO
C
	P(N-1)=(B(N-1)-P(N)-E(N-2))/(D(N-2)+4.)
	DO I=N-2,2,-1
	   P(I)=D(I)*P(I+1)+E(I)
	ENDDO
C
	RETURN
	END
C
C*************************************************************
      FUNCTION ZETA(N,CHI)
C
C  Computes ZETA, given N and CHI
C
      PARAMETER(ips=2000,IPC=20)
C
      COMMON/CFIT/AR(IPS,IPC),BR(IPS,IPC),CR(IPS,IPC),DR(IPS,IPC),
     +            AG(IPS,IPC),BG(IPS,IPC),CG(IPS,IPC),DG(IPS,IPC),
     +            AZ(IPS,IPC),BZ(IPS,IPC),CZ(IPS,IPC),DZ(IPS,IPC),
     +            CN,DC,CL,CH,MC
C
      CHIL=LOG10(CHI)
      IF(CHIL.LT.CL.or.CHIL.GT.CH)THEN
	 WRITE(6,6000)N,CHIL,CL,CH
	 STOP
      ENDIF
C
      M=1+(CHIL-CN)/DC
      M=MAX(1,M)
      M=MIN(MC-1,M)
      U=CHIL-(CN+(M-1)*DC)
      ZETA=EXP(AZ(N,M)+U*(BZ(N,M)+U*(CZ(N,M)+U*DZ(N,M))))
C
      RETURN
C
 6000 FORMAT('  SUBROUTINE ZETA: N=',I5,', CHIL=',1P,E10.3/
     + '  CL=',E10.3,', CH=',E10.3/
     + '  ALLOWED RANGE FOR CHIL:'/
     + '     CL .LE. CHIL .LE. CH')
C
      END
C
C**********************************
      FUNCTION ROSS(N,CHI)
C
C  Computes ROSS, given N and CHI
C
      PARAMETER(ips=2000,IPC=20)
C
      COMMON/CFIT/AR(IPS,IPC),BR(IPS,IPC),CR(IPS,IPC),DR(IPS,IPC),
     +            AG(IPS,IPC),BG(IPS,IPC),CG(IPS,IPC),DG(IPS,IPC),
     +            AZ(IPS,IPC),BZ(IPS,IPC),CZ(IPS,IPC),DZ(IPS,IPC),
     +            CN,DC,CL,CH,MC
C
      CHIL=LOG10(CHI)
      IF(CHIL.LT.CL.or.CHIL.GT.CH)THEN
	 WRITE(6,6000)N,CHIL,CL,CH
	 STOP
      ENDIF
C
      M=1+(CHIL-CN)/DC
      M=MAX(1,M)
      M=MIN(MC-1,M)
      U=CHIL-(CN+(M-1)*DC)
      ROSS=EXP(AR(N,M)+U*(BR(N,M)+U*(CR(N,M)+U*DR(N,M))))
C
      RETURN
C
 6000 FORMAT('  SUBROUTINE ROSS: N=',I5,', CHIL=',1P,E10.3/
     + '  CL=',E10.3,', CH=',E10.3/
     + '  ALLOWED RANGE FOR CHIL:'/
     + '     CL .LE. CHIL .LE. CH')
C
      END
C
C**********************************
      FUNCTION GRAD(N,CHI)
C
C  Computes GRAD, given N and CHI
C
      PARAMETER(ips=2000,IPC=20)
C
      COMMON/CFIT/AR(IPS,IPC),BR(IPS,IPC),CR(IPS,IPC),DR(IPS,IPC),
     +            AG(IPS,IPC),BG(IPS,IPC),CG(IPS,IPC),DG(IPS,IPC),
     +            AZ(IPS,IPC),BZ(IPS,IPC),CZ(IPS,IPC),DZ(IPS,IPC),
     +            CN,DC,CL,CH,MC
C
      CHIL=LOG10(CHI)
      IF(CHIL.LT.CL.or.CHIL.GT.CH)THEN
	 WRITE(6,6000)N,CHIL,CL,CH
	 STOP
      ENDIF
C
      M=1+(CHIL-CN)/DC
      M=MAX(1,M)
      M=MIN(MC-1,M)
      U=CHIL-(CN+(M-1)*DC)
      GRAD=EXP(AG(N,M)+U*(BG(N,M)+U*(CG(N,M)+U*DG(N,M))))
C
      RETURN
C
 6000 FORMAT('  SUBROUTINE GRAD: N=',I5,', CHIL=',1P,E10.3/
     + '  CL=',E10.3,', CH=',E10.3/
     + '  ALLOWED RANGE FOR CHIL:'/
     + '     CL .LE. CHIL .LE. CH')
C
      END
C
C**********************************
