C		ACCFIT.FOR
C
C       M.J. SEATON.  1997  January  3
C
C  OPERATION.
C
C    FILES acc.zz GIVE BASIC ATOMIC DATA FOR THE CALCULATION OF 
C    RADIATIVE FORCES, ON FIXED GRIDS OF TEMPERATURE (T), ELECTRON
C    DENSITY (N_e) AND ABUNDANCE MULTIPLIER (CHI).
C
C    THIS CODE PROVIDES FOR INTERPOLATIONS TO THE TEMPERATURE AND 
C    MASS-DENSITY POINTS OF A STELLAR MODELS, AT A SPECIFIED GRID
C    OF POINTS FOR CHI
C
C
C  PARAMETERS
C
C    IPI     FOR TEMPERATURE POINTS ON ARCHIVED FILE
C    IPK     FOR ELECTRON-DENSITY POINTS
C    IPM     FOR INPUT MESH OF CHI VALUES, CHI_IN
C    IPOUT   FOR OUTPUT  MESH OF CHI VALUES, CHIOUT
C    IPSTAR  FOR TEMPERATURE-DENSITY POINTS OF STELLAR MODEL
C
	PARAMETER (IPI=100,IPK=60,IPM=10,ipstar=2000,IPOUT=50)
C
C
	COMMON/COP1/KS(IPI),KE(IPI),K0,K2,A(IPI,IPK)
	COMMON/COP2/IS(IPK),IE(IPK),JS(IPI),JE(IPI),I0,I2,J0,CT,CN,
     +  F(IPI,IPK),FX(IPI,IPK),FY(IPI,IPK),FXY(IPI,IPK)
C
	COMMON/CINT/N1,T,S1,S2
C
	COMMON/CIN/
     +  EP0(IPI,IPK),EP1(IPI,IPK),ZET(IPI,IPK),FMU0,FMU1,AMACC,
     +  NCHI,CHI_IN(IPM),
     +  ROSS(IPI,IPK,IPM),ROSSP(IPI,IPK,IPM),
     +  GAM(IPI,IPK,IPM),GAMP(IPI,IPK,IPM)
	COMMON/CABUND/FANACC
C
C
	DIMENSION FLT(IPSTAR),FLR(IPSTAR),RI(IPSTAR),CHIOUT(IPOUT),
     +  RSSOUT(IPSTAR,IPOUT),GMMOUT(IPSTAR,IPOUT),ZETOUT(IPSTAR,IPOUT)   
C
	CHARACTER*128 FILE1,FILE2,FILE3,FILE4
C
C  FILES NAMES:-
C
C       FILE          CONTENTS
C
C       FILE1         ACC FILE, BASIC ATOMIC DATA
C       FILE2         REQUIRED VALUES OF LOG10(T), LOG10(RHO) AND RI
C       FILE3         REQUIRED VALUES OF LOG10(CHI)
C       FILE4         OUTPUT FOR ZETA, ROSS AND GRAD
C
C  READ FILE NAMES
C
	READ(5,'(A)')FILE1
	READ(5,'(A)')FILE2
	READ(5,'(A)')FILE3
	READ(5,'(A)')FILE4
C
C  OPEN FILES
	OPEN(1,FILE=FILE1,STATUS='OLD',ERR=1001)
	OPEN(2,FILE=FILE2,STATUS='OLD',ERR=1002)
	OPEN(3,FILE=FILE3,STATUS='OLD',ERR=1003)
	OPEN(4,FILE=FILE4,STATUS='NEW',ERR=1004)
C
C  READ FILE1
	CALL READ1
	CLOSE(1) 
C
C  READ FILE2
	CALL READ2(I0,I2,CT,FLT,FLR,RI,TEFF,NSTAR)
	CLOSE(2)
C
C  READ FILE3
	CALL READ3(CHIOUT,MCHI,CHIL0,DCHIL)
	CLOSE(3)
C
	CT3=3.*CT
	C=-5.77975+CN*K0-CT3*I0
C
C       START LOOP ON OUTPUT VALUES OF CHI
	WRITE(6,*)
	DO 100 M=1,MCHI
C
	   CHI=CHIOUT(M)
	   FMU=FMU0+CHI*FMU1
C
C          STORE LOG10(R_OPAL) IN A(I,K)
	   DO I=1,I2
	      DO K=KS(I),KE(I)
	         EPA=EP0(I,K)+CHI*EP1(I,K)
	         A(I,K)=C+LOG10(FMU/EPA)+CN*K-CT3*I
	      ENDDO
	   ENDDO
C
C          GET PARAMETERS FOR INTEERPOLATING FROM CHI_IN TO CHIOUT
	   CALL CHNT(CHI,NCHI,CHI_IN)
C
C          DO INTERPOLATIONS FOR ROSS
C          GET F(I,K)=LOG10(ROSS) INTERPLATED TO LOOP-VALUE OF CHI
	   DO I=1,I2
	      DO K=KS(I),KE(I)
	if(ross(i,k,n1+1).le.0)then
	   print*,'  i,k,n1=',i,k,n1
	   print*,'  ross(i,k,n1+1)=',ross(i,k,n1+1)
	endif
	         CALL INTRP1(ROSS(I,K,N1),ROSS(I,K,N1+1),
     +                      ROSSP(I,K,N1),ROSSP(I,K,N1+1),F(I,K),0) !!!
              ENDDO
           ENDDO
C          GET RSSOUT(N,M)=LOG10(ROSS) ON STAR MESH
	   CALL SMOOTH
	   CALL CHNTP
	   CALL DERIV
	   DO N=1,NSTAR
	      CALL INTRP2(FLT(N),FLR(N),RSSOUT(N,M),IERR)
           ENDDO
C
C          DO INTERPOLATIONS FOR GAMMA
C          GET F(I,K)=LOG10(GAM) INTERPOLATED  TTO LOOP-VALUE OF CHI
	   DO I=1,I2
	      DO  K=KS(I),KE(I)
	if(gam(i,k,n1+1).le.0)then
	   print*,'  i,k,n1=',i,k,n1
	   print*,'  gam(i,k,n1+1)=',gam(i,k,n1+1)
	endif
		 CALL INTRP1(GAM(I,K,N1),GAM(I,K,N1+1),
     +           GAMP(I,K,N1),GAMP(I,K,N1+1),F(I,K),1)
	      ENDDO
	   ENDDO
C          GET GMMOUUT(N,M)=LOG10(GAMMA) ON STAR  MESH
	   CALL SMOOTH
	   CALL CHNTP
	   CALL DERIV
	   DO N=1,NSTAR
	      CALL INTRP2(FLT(N),FLR(N),GMMOUT(N,M),IERR)
           ENDDO
C
C       DO   INTERPOLATIONS FOR ZETA
C       GET F(I,K)=LOG10(ZETA)
        DO I=1,I2
	   DO K=KS(I),KE(I)
	      F(I,K)=LOG10(ZET(I,K))
	   ENDDO
	ENDDO
C       GET ZETOUT(N,M)=LOG10(ZETA) ON STAR MESH
	CALL  SMOOTH
	CALL CHNTP
	CALL DERIV
	DO N=1,NSTAR
	   CALL INTRP2(FLT(N),FLR(N),ZETOUT(N,M),IERR)
	ENDDO
C
	WRITE(6,602)LOG10(CHI)
C
  100   CONTINUE
C
C  CHANGE TO:
C     ZETOUT=FACTOR ZETA;
C     RSSOUT=ROSSELAND  MEAN IN  cgs;
C     GRAD=RADIATIVE ACCELERATION IN cgs.  
	CG0=1.89148E-15*TEFF**4/AMACC
	DO N=1,NSTAR
	   CG=CG0*RI(N)**(-2)
	   DO M=1,MCHI
	      FMU=FMU0+FMU1*CHIOUT(M)
	      RSSOUT(N,M)=RSSOUT(N,M)+7.226956-LOG10(FMU)
	      GMMOUT(N,M)=CG*FMU*10**(GMMOUT(N,M)+RSSOUT(N,M))
	      RSSOUT(N,M)=10**RSSOUT(N,M)
	      ZETOUT(N,M)=10**ZETOUT(N,M)
	   ENDDO
	ENDDO
C
C  WRITE OUTPUT FILE
	WRITE(4,4000)NSTAR,CHIL0,DCHIL,MCHI
	DO N=1,NSTAR
	   WRITE(4,4001)N,FLT(N),FLR(N)  
	   WRITE(4,4002)(ZETOUT(N,M),M=1,MCHI) 
	   WRITE(4,4002)(RSSOUT(N,M),M=1,MCHI)
	   WRITE(4,4002)(GMMOUT(N,M),M=1,MCHI)
        ENDDO
	CLOSE(4,STATUS='KEEP')
C
	WRITE(6,603)FILE4
C
	STOP
C
C  MESSAGES FOR ERRORS IN  OPENING FILES
 1001   WRITE(6,6001)FILE1
	STOP
 1002   WRITE(6,6002)FILE2
	STOP
 1003   WRITE(6,6003)FILE3
	STOP
 1004   WRITE(6,6004)FILE4
	STOP
C
C  CONTENTS OF OUTPUT FILE, FILE4
C     NSTAR=MUMBER OF STAR POINTS
C     CHIL0=SMALLEST VALUE OF LOG10(CHI)
C     DCHIL=INTERVAL IN LOG10(CHI)
C     MCHI=number of CHI values
C     Start loop n=1,NSTAR
C         N, FLT(N)=LOG10(T), FLR(N)=LOG10(RHO)
C         (ZETA(N,M),M=1,MCHI)
C         (ROSS(N,M),M=1,MCHI) 
C         (GRAD(N,M),M=1,MCHI)
C     WHERE:
C          ZETA=FACTOR FOR DIFFUSION COEFFICIENT;
C          ROSS=ROSSELAND-MEAN OPACITY IN cgs;
C          GRAD=RADIATIVE ACCELERRATION IN cgs.
C
C  NOTE ON ERROR CODE IERR
C  EXECUTION OF
C     CALL INTRP2(FLT(N),FLR(N),G,IERR)
C  RETURNS AN ERROR CODE IERR=1 IF FLT(N) OR FLR(N) ARE OUT OF
C  THE RANGE INCLUDED IN INPUT DATA FROM FILE1.
C  IF IERR=1:
C     A MESSAGE IS PRINTED BY INTRP2;
C     INTRP2 SETS G=9.999.
C
  602   FORMAT(10X,'Processed log(CHI)=',1p,e11.3)
  603   FORMAT(/2X,'Output file written:-'/2x,A128)
 6001   FORMAT('  Error opening FILE1'/1X,A128)
 6002   FORMAT('  Error opening FILE2'/1X,A128)
 6003   FORMAT('  Error opening FILE3'/1X,A128)
 6004   FORMAT('  Error opening FILE4'/1X,A128)
 4000   FORMAT(I5,1P,2E11.3,0P,I5)
 4001   FORMAT(I5,1P,2E11.3)
 4002   FORMAT(1P,7E11.3)
C
	END
C*****************************************************************
	SUBROUTINE READ1
C
C   READS AND PROCESSES BASIC ATOMIC DATA FROM FILE1
C
	parameter (ipi=100,ipk=60,IPM=10)
c
	CHARACTER*80 HEAD
C
	COMMON/COP1/KS(IPI),KE(IPI),K0,K2,A(IPI,IPK)
	COMMON/COP2/IS(IPK),IE(IPK),JS(IPI),JE(IPI),I0,I2,J0,CT,CN,
     +  F(IPI,IPK),FX(IPI,IPK),FY(IPI,IPK),FXY(IPI,IPK)
C
	COMMON/CIN/
     +  EP0(IPI,IPK),EP1(IPI,IPK),ZET(IPI,IPK),FMU0,FMU1,AMACC,
     +  NCHI,CHI_IN(IPM),
     +  ROSS(IPI,IPK,IPM),ROSSP(IPI,IPK,IPM),
     +  GAM(IPI,IPK,IPM),GAMP(IPI,IPK,IPM)
C
	COMMON/CABUND/FANACC
C
C
C   READ:
C     IZACC=NUCLEAR CHARGE OF SELECTED ELEMENT
C     NTOT=NUMBER OF FREQUENCY POINTS
C     NCHI=NUMBER OF  CHI  VALUES
C     FMU0, FMU1 SUCH THAT FMU=FMU0+CHI*FMU1
C     AMACC=ATOMIC MASS FOR SELECTED ELEMENT
C     FANACC=ABUNDANCE OF SELECTED ELEMENT FOR CHI=1
	READ(1,*)IZACC,NTOT,NCHI,FMU0,FMU1
	READ(1,*)AMACC,FANACC
C
C       VALUES OF CHI
C       (NOTE THAT CHI_IN(N)=INPUT VALUE OF CHI)
	READ(1,*)(CHI_IN(N),N=1,NCHI)
C
C       GET SCALE FACTORS CT, I0, I2
	READ(1,*)IT1,IT2,IT3
	ct=0.025*it3
	i0=it1/it3-1
	i2=it2/it3-i0
C
C       GET SCALE PARAMETERS K0 AND CN
	read(1,*)it,jn1,jn2,jn3
	cn=0.25*jn3
	k0=jn1/jn3-1
	backspace(1)
C
C       START LOOP ON TEMPERATURES
	DO 30 i=1,I2
C          READ TEMPERATURE INDEX IT AND ELECTRON-DENSITY INDICES JN1,JN2,JN3
	   READ(1,*)IT,JN1,JN2,JN3
	   CALL CHCKI(I,IT,IT1,IT3)
	   ks(i)=jn1/jn3-k0
	   ke(i)=jn2/jn3-k0
C
C          START LOOP ON DENSITIES          
	   DO 20 k=ks(i),ke(i)
C             READ DATA ON NUMBER OF ELECTRONS PER ATOM
C             ELECTRON PER ATOM = EP0+CHI*EP1
 	      READ(1,*)JN,EP0(i,k),EP1(i,k)
	      CALL  CHCKJ(K,JN,JN1,I,KS(I),JN3)
	      READ(1,*)ZET(I,K)
C             READ DATA FOR ROSS, GAM AND DERIVATIVES  
	      READ(1,*)(ROSS(I,K,N) ,N=1,NCHI)
	      READ(1,*)(ROSSP(I,K,N),N=1,NCHI)
	      READ(1,*)(GAM(I,K,N),  N=1,NCHI)
	      READ(1,*)(GAMP(I,K,N), N=1,NCHI)
C
   20      CONTINUE
C
   30   CONTINUE
C
C
C     ADJUST INDICES KS(I), KE(I) FOR RANGE OF K
      DO 90 I=1,I2
         DO 80 K=KS(I),KE(I)
	    IF(EP0(I,K).LE.0)THEN
               KE(I)=K-1
               GOTO 90
            ENDIF
   80    CONTINUE
   90 CONTINUE
C
      K2=KE(I2)
C
      DO 100 I=I2-1,1,-1
         IF(KS(I).GT.KS(I+1))KS(I)=KS(I+1)
         IF(KE(I).GT.KE(I+1))KE(I)=KE(I+1)
  100 CONTINUE
C
C     GET INDICES IS(K), IE(K) FOR RANGE OF I
      KB=0
      DO 120 I=1,I2
         IF(KE(I).GT.KB)THEN
            KA=KB+1
            KB=KE(I)
            DO 110 K=KA,KB
               IS(K)=I
  110       CONTINUE
         ENDIF
  120 CONTINUE
      KA=K2+1
      DO 140 I=I2,1,-1
         IF(KS(I).LT.KA)THEN
            KB=KA-1
            KA=KS(I)
            DO 130 K=KA,KB
  130       IE(K)=I
         ENDIF
  140 CONTINUE
C
C
C  CHECK FOR GAM .GE. 0.
	NN=NCHI+1
	DO I=1,I2
	   DO K=KS(I),KE(I)
	      DO N=1,NCHI
		 IF(GAM(I,K,N).LE.0)NN=MIN(N,NN)
	      ENDDO
	   ENDDO
	ENDDO
	IF(NN.LE.NCHI)THEN
	   WRITE(6,601)NCHI,NN-1
	   NCHI=NN-1
	ENDIF
C
	WRITE(6,600)
C
	RETURN
C
  600 FORMAT(/2X,'Completed READ1')
  601 FORMAT(/2X,'Changed number of CHI values from',i3,
     + ' to',i3/2x,'in order to avoid negative values of GAM')  
C
C DATA OBTAINED IN READ1:-
C    AMACC=mass of accelerated ion (all mases in amu)
C    Mean  atom  mass=FMU0+CHI*FMU1
C    Input CHI, (CHI_IN(N), N=1,NCHI)
C    Temperature index, I=1 to I2 in steps of 1
C    Density index, K=1 to K2  in  steps of 1
C    For given I, K=KS(I) to KE(I)
C    For given K, I=IS(K) to IE(K)
C    LOG10(T)=CT*(I+I0)
C    LOG10(N_e)=CN*(K+K0)
C    FOR EACH I,K
C        Electrons per atom=EP0(I,K)+CHI*EP1(I,K)
C        Factor in diffusion coefficient=ZET(I,K)
C    For each I, K, N:
C        ROSS=Rosseland-mean cross-section  (in a.u.)
C        ROSSP=derivative w.r.t. CHI
C        GAM=gamma
C        GAMP=derivative w.r.t. CHI
C
	END
C************************************************
	SUBROUTINE READ2(I0,I2,CT,FLT,FLR,RI,TEFF,NSTAR)
C
C  READS DATA FOR STELLAR MODEL
C
	PARAMETER(ipstar=2000)
	DIMENSION FLT(IPSTAR),FLR(IPSTAR),RI(IPSTAR)
C
C  ALLOWED RANGE OF LOG(T)
	FLT1=CT*(1+I0)
	FLT2=CT*(I2+I0)
C
	N=0
	N0=0
	N1=0
	NSTAR=0
C
	READ(2,*)TEFF
    1   READ(2,*,END=2)GLT,GLR,GRI
	IF(GLT.LT.FLT1)THEN
	   N0=N0+1
	ELSEIF(GLT.GT.FLT2)THEN
	   N1=N1+1
	ELSEIF(NSTAR.LT.IPSTAR)THEN
	   NSTAR=NSTAR+1
	   N=N+1
	   FLT(N)=GLT
	   FLR(N)=GLR
	   RI(N)=GRI
	ELSE
	   N=N+1
	ENDIF
	GOTO 1
C
    2   WRITE(6,6000)NSTAR,FLT(1),FLR(1),FLT(NSTAR),FLR(NSTAR)
	IF(N0.GT.0)WRITE(6,6001)N0,FLT1
	IF(N1.GT.0)WRITE(6,6002)N1,FLT2
	IF(N.GT.NSTAR)WRITE(6,6003)N
C
C  DATA  OBTAINED IN READ2
C    TEFF, EFFECTIVE  TEMPERATURE
C    NSTAR, NUMBER  OF STAR POINTS
C    FOR N=1,NSTAR:
C      FLT(N)=LOG10(T)
C      FLR(N)=LOG10(RHO)
C      RI(N)=(r/R_*)
C
 6000   FORMAT(/2X,'Completed READ2'/
     +  /2X,'NSTAR=',I6/1P,
     +  2X,'FIRST, LOG10(T)=',E10.3,', LOG10(RHO)=',E10.3/
     +  2X,'LAST , LOG10(T)=',E10.3,', LOG10(RHO)=',E10.3/)
 6001   FORMAT(/2X,'OMITTED ',I5,'  POINTS WITH LOG10(T) LESS THAN',
     +  1P,E10.3)
 6002   FORMAT(/2X,'OMITTED ',I5,' POINTS WITH  LOG10(T) GREATER THAN',
     +  1P,E10.3)
 6003   FORMAT(/2X,'NEED PARAMETER IPSTAR=',I5,' TO INCLUDE ALL',
     +  ' STAR  POINTS'/)
C
	END
C******************************************************************
	SUBROUTINE READ3(CHIOUT,MCHI,CHIL0,DCHIL)
C
C  READS DATA FOR OUTPUT VALUES OF CHI
C
	PARAMETER (IPI=100,IPK=60,IPM=10,ipstar=2000,IPOUT=50)
	DIMENSION CHIOUT(IPOUT)
C
	COMMON/CIN/
     +  EP0(IPI,IPK),EP1(IPI,IPK),ZET(IPI,IPK),FMU0,FMU1,AMACC,
     +  NCHI,CHI_IN(IPM),
     +  ROSS(IPI,IPK,IPM),ROSSP(IPI,IPK,IPM),
     +  GAM(IPI,IPK,IPM),GAMP(IPI,IPK,IPM)
C
C  MCHI=NUMBER OF CHI  VALUES
	READ(3,*)CHIL0,DCHIL,MCHI
C       CHECK DIMENSIONS
	IF(MCHI.GT.IPOUT)THEN
	   WRITE(6,600)MCHI,IPOUT
	   STOP
	ENDIF
C
C  CHECK RANGE FOR CHI
C
C  SMALLEST LOG100(CHI) FROM FILE1
	CHLMIN=LOG10(CHI_IN(1))
C  AND HIGHEST
	CHLMAX=LOG10(CHI_IN(NCHI))
C  SMALLEST LOG10(CHI) FROM FILE3=CHIL00
C  AND HIGHEST
	CHILM=CHIL0+(MCHI-1)*DCHIL
C
C  CHECK THAT CHIL0.GE.CHLMIN
	IF(CHIL0.LT.CHLMIN)THEN
	   DO M=2,MCHI
	     IF(CHIL0+(M-1)*DCHIL.GE.CHLMIN)THEN
		C0=CHIL0
		M0=MCHI
		CHIL0=CHIL0+(M-1)*DCHIL
		MCHI=MCHI-M+1
		WRITE(6,6000)C0,CHLMIN,CHIL0,M0,MCHI
		GOTO 1
	     ENDIF
	   ENDDO
	   WRITE(6,6001)CHILM,CHLMIN
	   STOP
	ENDIF
C
C  CHECK THAT CHILM.LE.CHLMAX
    1	IF(CHILM.GT.CHLMAX)THEN
	   DO M=MCHI-1,1,-1
	      IF(CHIL0+DCHIL*(M-1).LE.CHLMAX)THEN
		WRITE(6,6002)CHILM,CHLMAX,MCHI,M
		MCHI=M
	        GOTO 2
	      ENDIF
	    ENDDO
	    WRITE(6,6003)CHIL0,CHLMAX
	    STOP
	ENDIF
    2   CONTINUE
C
C  STORE CHI IN ARRAY CHIOUT(M)
	DO  M=1,MCHI
	   CHIOUT(M)=10**(CHIL0+(M-1)*DCHIL)
        ENDDO
C
C  CHECK THAT MCHI.GE.5 (A CONDITI0N REQUIRED FOR INTERPOLATION ROUTINES)
	IF(MCHI.LT.5)THEN
	    WRITE(6,6005)MCHI
	    STOP
	ENDIF
C
	WRITE(6,6004)CHIL0,DCHIL,MCHI
C
	RETURN
C
C  DATA OBTAINED IN READ3
C  CHI DATA FOR OUTPUT FILE
C    CHIL0=LOWEST  LOG10(CHI)
C    DCHIL=ITERVAL IN LOG10(CHI)
C    MCHI=NUMBER OF CHI VALUES
C    (CHIOUT(M),M=1,MCHI)=VALUES OF CHI
C
  600 FORMAT(/5X,'MCHI=',I5,', IPOUT=',I5/
     + 5X,'Require parameter IPOUT .GE. MCHI')
 6000 FORMAT('  Lowest LOG10(CHI) from FILE3=CHILO=',1P,E10.3/
     +       '  Lowest LOG10(CHI) from FILE1=',E10.3/
     +       '  Changed to CHIL0=',E10.3/
     +       '  and changed number of CHI values from ',0P,I5,
     +       ' to ',I5)
 6001 FORMAT('  Highest LOG10(CHI) from FILE3=',1P,E10.3/
     +       '  Lowest  LOG10(CHI) from FILE1=',E10.3/
     +       '  STOP')
 6002 FORMAT('  Highest LOG10(CHI) from FILE3=',1P,E10.3/
     +       '  Highest LOG10(CHI) from FILE1=',E10.3/
     +       '  Reduced MCHI from ',0P,I5,' to ',I5)
 6003 FORMAT('  Lowest LOG10(CHI) from FILE3=',1P,E10.3/
     +       '  Highest LOG10(CHI) from FILE1=',E10.3/
     +       '  STOP')
 6004 FORMAT(/'  LOG10(CHI) from READ3:'/
     + 5X,'lowest=',1P,E10.3/
     + 5X,'interval=',E10.3/
     + 5X,'number of values=',0P,I5)      
 6005 FORMAT('  Have MCHI=',I2,', require MCHI.GE.5'/
     +       '  STOP')
C
	END
C*******************************************************************
	SUBROUTINE CHNT(CHI,NCHI,CHI_IN)
C
C  SETS UP VARIABLES N1,T,S1 AND S2 FOR INTERPOLATIONS IN CHI.
C  THOSE VARIABLES KEPT IN /CINT/
C
	PARAMETER(IPM=10)
	DIMENSION CHI_IN(IPM)
	COMMON/CINT/N1,T,S1,S2
C
	IF(CHI.LE.CHI_IN(1))THEN
	   N1=1
        ELSEIF(CHI.GE.CHI_IN(NCHI))THEN
	   N1=NCHI-1
	ELSE
	   DO 10 N=NCHI,1,-1
	      IF(CHI.GE.CHI_IN(N))THEN
		 N1=N
		 GOTO 11
	      ENDIF
   10      CONTINUE
   11      CONTINUE
	ENDIF
C
	S=LOG10(CHI_IN(N1+1)/CHI_IN(N1))
	T=LOG10(CHI/CHI_IN(N1))/S
	S1=S*CHI_IN(N1)
	S2=S*CHI_IN(N1+1)
C
	RETURN
C
	END
C
C***********************************************************	
	SUBROUTINE INTRP1(Y1,Y2,Y1P,Y2P,Y,irg) !!!
C
C  DOES INTERPOLATIONS IN CHI, USING CUBIC FITS.
C
	COMMON/CINT/N1,T,S1,S2
C
	if(y1.le.0)then
	print*,' intrp1, y1=',y1
	if(irg.eq.0)then
	     print*,' ROSS'
	else
	     print*,'  GAM'
	endif
	print*,' stop'
	stop
	endif
	if(y2.le.0)then
	print*,'  intrp1, y2=',y2
	if(irg.eq.0)then
	     print*,' ROSS'
	else
	     print*,'  GAM'
	endif
	print*,' stop'
	stop
	endif
	A=LOG10(Y1)
	B=LOG10(Y2)
	C=S1*(Y1P/Y1)
	D=S2*(Y2P/Y2)
C
	Y=A+T*(C+T*(3.*(B-A)-2.*C-D+T*(2.*(A-B)+C+D)))
C
	RETURN
C
	END
C
C****************************************************************************
      SUBROUTINE SMOOTH
C
C  SMOOTHING ROUTINE AS USED IN OPFIT.FOR
C
C  SMOOTHING OF OP DATA FOR LOG10(T).GE.4.1
C  (BETTER NOT SMOOTHED FOR SMALLER LOG10(T))
C
      PARAMETER (IPI=100,IPK=60,IPN=30,IPJ=60)
      COMMON/COP1/KS(IPI),KE(IPI),K0,K2,A(IPI,IPK)
      COMMON/COP2/IS(IPK),IE(IPK),JS(IPI),JE(IPI),I0,I2,J0,CT,CN,
     + F(IPI,IPK),FX(IPI,IPK),FY(IPI,IPK),G(IPI,IPK)
C
      DIMENSION ALP(11)
      DATA ALP/
     + -0.0844897959,-0.0048979592,+0.0073469388,+0.0012244898,
     +  0.3379591837,+0.0195918367,-0.0293877551,+0.4787755102,
     +  0.0277551020,-0.0416326531,-0.0069387755/
      DIMENSION BET(11)
      DATA BET/
     + -0.0048979592,-0.0661224490,-0.0293877551,+0.0195918367,
     +  0.2644897959,+0.1175510204,-0.0783673469,+0.0277551020,
     +  0.3746938776,+0.1665306122,-0.1110204082/
      DIMENSION GAM(6)
      DATA GAM/+0.0073469388,-0.0293877551,-0.0416326531,
     +         +0.1175510204,+0.1665306122,+0.2359183673/
C
         H0(I,J)=F(I,J)
         H1(I,J)=
     +    ALP(1)*( F(I-2,J  )+F(I+2,J  ) )
     +   +ALP(2)*( F(I-2,J+1)+F(I+2,J+1)+F(I-2,J+3)+F(I+2,J+3)
     +          +F(I-1,J+4)+F(I+1,J+4) )
     +   +ALP(3)*( F(I-2,J+2)+F(I+2,J+2) )
     +   +ALP(4)*( F(I-2,J+4)+F(I+2,J+4) )
     +   +ALP(5)*( F(I-1,J  )+F(I+1,J  ) )
     +   +ALP(6)*( F(I-1,J+1)+F(I+1,J+1)+F(I-1,J+3)+F(I+1,J+3) )
     +   +ALP(7)*( F(I-1,J+2)+F(I+1,J+2) )
     +   +ALP(8)*  F(I  ,J  )
     +   +ALP(9)*( F(I  ,J+1)+F(I  ,J+3) )
     +   +ALP(10)* F(I  ,J+2) +ALP(11)*F(I  ,J+4)
         H2(I,J)=
     +    BET(1)*( F(I-2,J-1)+F(I+2,J-1)+F(I-2,J+3)+F(I+2,J+3) )
     +   +BET(2)*( F(I-2,J  )+F(I+2,J  ) )
     +   +BET(3)*( F(I-2,J+1)+F(I+2,J+1) )
     +   +BET(4)*( F(I-2,J+2)+F(I+2,J+2)+F(I-1,J-1)+F(I+1,J-1)
     +            +F(I-1,J+3)+F(I+1,J+3) )
     +   +BET(5)*( F(I-1,J  )+F(I+1,J  ) )
     +   +BET(6)*( F(I-1,J+1)+F(I+1,J+1) )
     +   +BET(7)*( F(I-1,J+2)+F(I+1,J+2) )
     +   +BET(8)*( F(I  ,J-1)+F(I  ,J+3) )
     +   +BET(9)*F(I  ,J  ) +BET(10)*F(I  ,J+1) +BET(11)*F(I  ,J+2)
C
         H3(I,J)=
     +    GAM(1)*( F(I-2,J-2)+F(I-2,J+2)+F(I+2,J-2)+F(I+2,J+2) )
     +   +GAM(2)*( F(I-2,J+1)+F(I-2,J-1)+F(I-1,J-2)+F(I-1,J+2)
     +            +F(I+1,J-2)+F(I+1,J+2)+F(I+2,J-1)+F(I+2,J+1) )
     +   +GAM(3)*( F(I-2,J  )+F(I+2,J  )+F(I  ,J-2)+F(I  ,J+2) )
     +   +GAM(4)*( F(I-1,J-1)+F(I-1,J+1)+F(I+1,J-1)+F(I+1,J+1) )
     +   +GAM(5)*( F(I-1,J  )+F(I  ,J-1)+F(I  ,J+1)+F(I+1,J  ) )
     +   +GAM(6)*  F(I  ,J  )
         H4(I,J)=
     +    BET(1)*( F(I-2,J+1)+F(I+2,J+1)+F(I-2,J-3)+F(I+2,J-3) )
     +   +BET(2)*( F(I-2,J  )+F(I+2,J  ) )
     +   +BET(3)*( F(I-2,J-1)+F(I+2,J-1) )
     +   +BET(4)*( F(I-2,J-2)+F(I+2,J-2)+F(I-1,J+1)+F(I+1,J+1)
     +            +F(I-1,J-3)+F(I+1,J-3) )
     +   +BET(5)*( F(I-1,J  )+F(I+1,J  ) )
     +   +BET(6)*( F(I-1,J-1)+F(I+1,J-1) )
     +   +BET(7)*( F(I-1,J-2)+F(I+1,J-2) )
     +   +BET(8)*( F(I  ,J+1)+F(I  ,J-3) )
     +   +BET(9)*F(I  ,J  ) +BET(10)*F(I  ,J-1) +BET(11)*F(I  ,J-2)
         H5(I,J)=
     +    ALP(1)*( F(I-2,J  )+F(I+2,J  ) )
     +   +ALP(2)*( F(I-2,J-1)+F(I+2,J-1)+F(I-2,J-3)+F(I+2,J-3)
     +            +F(I-1,J-4)+F(I+1,J-4) )
     +   +ALP(3)*( F(I-2,J-2)+F(I+2,J-2) )
     +   +ALP(4)*( F(I-2,J-4)+F(I+2,J-4) )
     +   +ALP(5)*( F(I-1,J  )+F(I+1,J  ) )
     +   +ALP(6)*( F(I-1,J-1)+F(I+1,J-1)+F(I-1,J-3)+F(I+1,J-3) )
     +   +ALP(7)*( F(I-1,J-2)+F(I+1,J-2) )
     +   +ALP(8)*  F(I  ,J  )
     +   +ALP(9)*( F(I  ,J-1)+F(I  ,J-3) )
     +   +ALP(10)* F(I  ,J-2) +ALP(11)*F(I  ,J-4)
         H6(I,J)=
     +    ALP(1)*( F(I  ,J-2)+F(I  ,J+2) )
     +   +ALP(2)*( F(I+1,J-2)+F(I+1,J+2)+F(I+3,J-2)+F(I+3,J+2)
     +          +F(I+4,J-1)+F(I+4,J+1) )
     +   +ALP(3)*( F(I+2,J-2)+F(I+2,J+2) )
     +   +ALP(4)*( F(I+4,J-2)+F(I+4,J+2) )
     +   +ALP(5)*( F(I  ,J-1)+F(I  ,J+1) )
     +   +ALP(6)*( F(I+1,J-1)+F(I+1,J+1)+F(I+3,J-1)+F(I+3,J+1) )
     +   +ALP(7)*( F(I+2,J-1)+F(I+2,J+1) )
     +   +ALP(8)*  F(I  ,J  )
     +   +ALP(9)*( F(I+1,J  )+F(I+3,J  ) )
     +   +ALP(10)* F(I+2,  J) +ALP(11)*F(I+4,J  )
         H7(I,J)=
     +    BET(1)*( F(I-1,J-2)+F(I-1,J+2)+F(I+3,J-2)+F(I+3,J+2) )
     +   +BET(2)*( F(I  ,J-2)+F(I  ,J+2) )
     +   +BET(3)*( F(I+1,J-2)+F(I+1,J+2) )
     +   +BET(4)*( F(I+2,J-2)+F(I+2,J+2)+F(I-1,J-1)+F(I-1,J+1)
     +            +F(I+3,J-1)+F(I+3,J+1) )
     +   +BET(5)*( F(I  ,J-1)+F(I  ,J+1) )
     +   +BET(6)*( F(I+1,J-1)+F(I+1,J+1) )
     +   +BET(7)*( F(I+2,J-1)+F(I+2,J+1) )
     +   +BET(8)*( F(I-1,J  )+F(I+3,J  ) )
     +   +BET(9)*F(I  ,J  ) +BET(10)*F(I+1,J  ) +BET(11)*F(I+2,J  )
         H8(I,J)=
     +    BET(1)*( F(I+1,J-2)+F(I+1,J+2)+F(I-3,J-2)+F(I-3,J+2) )
     +   +BET(2)*( F(I  ,J-2)+F(I  ,J+2) )
     +   +BET(3)*( F(I-1,J-2)+F(I-1,J+2) )
     +   +BET(4)*( F(I-2,J-2)+F(I-2,J+2)+F(I+1,J-1)+F(I+1,J+1)
     +            +F(I-3,J-1)+F(I-3,J+1) )
     +   +BET(5)*( F(I  ,J-1)+F(I  ,J+1) )
     +   +BET(6)*( F(I-1,J-1)+F(I-1,J+1) )
     +   +BET(7)*( F(I-2,J-1)+F(I-2,J+1) )
     +   +BET(8)*( F(I+1,J  )+F(I-3,J  ) )
     +   +BET(9)*F(I  ,J  ) +BET(10)*F(I-1,J  ) +BET(11)*F(I-2,J  )
         H9(I,J)=
     +    ALP(1)*( F(I  ,J-2)+F(I  ,J+2) )
     +   +ALP(2)*( F(I-1,J-2)+F(I-1,J+2)+F(I-3,J-2)+F(I-3,J+2)
     +            +F(I-4,J-1)+F(I-4,J+1) )
     +   +ALP(3)*( F(I-2,J-2)+F(I-2,J+2) )
     +   +ALP(4)*( F(I-4,J-2)+F(I-4,J+2) )
     +   +ALP(5)*( F(I  ,J-1)+F(I  ,J+1) )
     +   +ALP(6)*( F(I-1,J-1)+F(I-1,J+1)+F(I-3,J-1)+F(I-3,J+1) )
     +   +ALP(7)*( F(I-2,J-1)+F(I-2,J+1) )
     +   +ALP(8)*  F(I  ,J  )
     +   +ALP(9)*( F(I-1,J  )+F(I-3,J  ) )
     +   +ALP(10)* F(I-2,J  ) +ALP(11)*F(I-4,J  )
C
C
C  SET LOWER LIMIT ON LOG10(T)
      IF(CT*(1+I0).GE.4.101)THEN
	 I1=1
      ELSE
         I1=4.101/CT-I0
      ENDIF
C
      DO 20 I=I1,I2
         DO 10 K=KS(I),KE(I)
         IF    (K.GE.3   .AND.I.LE.IE(K-2)-2.AND.   ! N=2
     +          K.LE.K2-2.AND.I.GE.IS(K+2)+2.AND.   ! M=2
     +          I.GE.3   .AND.K.LE.KE(I-2)-2.AND.   ! P=2
     +          I.LE.I2-2.AND.K.GE.KS(I+2)+2)THEN   ! Q=2
                   G(I,K)=H3(I,K)
         ELSEIF(K.GE.2   .AND.I.LE.IE(K-1)-2.AND.   ! N=2
     +          K.LE.K2-3.AND.I.GE.IS(K+3)+2.AND.   ! M=2
     +          I.GE.3   .AND.K.LE.KE(I-2)-3.AND.   ! P=1
     +          I.LE.I2-2.AND.K.GE.KS(I+2)+1)THEN   ! Q=3
                   G(I,K)=H2(I,K)
         ELSEIF(K.GE.4   .AND.I.LE.IE(K-3)-2.AND.   ! N=2
     +          K.LE.K2-1.AND.I.GE.IS(K+1)+2.AND.   ! M=2 
     +          I.GE.3.  .AND.K.LE.KE(I-2)-1.AND.   ! P=3
     +          I.LE.I2-2.AND.K.GE.KS(I+2)+3)THEN   ! Q=1
                   G(I,K)=H4(I,K)  
         ELSEIF(K.GE.3   .AND.I.LE.IE(K-2)-3.AND.   ! N=1
     +          K.LE.K2-2.AND.I.GE.IS(K+2)+1.AND.   ! M=3
     +          I.GE.2.  .AND.K.LE.KE(I-1)-2.AND.   ! P=2
     +          I.LE.I2-3.AND.K.GE.KS(I+3)+2)THEN   ! Q=2
                   G(I,K)=H7(I,K)
         ELSEIF(K.GE.3.  .AND.I.LE.IE(K-2)-1.AND.   ! N=3
     +          K.LE.K2-2.AND.I.GE.IS(K+2)+3.AND.   ! M=1
     +          I.GE.4.  .AND.K.LE.KE(I-3)-2.AND.   ! P=2
     +          I.LE.I2-1.AND.K.GE.KS(I+1)+2)THEN   ! Q=2
                   G(I,K)=H8(I,K)
         ELSEIF(K.GE.1   .AND.I.LE.IE(K  )-2.AND.   ! N=2
     +          K.LE.K2-4.AND.I.GE.IS(K+4)+2.AND.   ! M=2
     +          I.GE.3   .AND.K.LE.KE(I-2)-4.AND.   ! P=0
     +          I.LE.I2-2.AND.K.GE.KS(I+2)  )THEN   ! Q=4
                   G(I,K)=H1(I,K)  
         ELSEIF(K.GE.5.  .AND.I.LE.IE(K-4)-2.AND.   ! N=2
     +          K.LE.K2  .AND.I.GE.IS(K  )+2.AND.   ! M=2
     +          I.GE.3.  .AND.K.LE.KE(I-2)  .AND.   ! P=4
     +          I.LE.I2-2.AND.K.GE.KS(I+2)+4)THEN   ! Q=0
                   G(I,K)=H5(I,K)
         ELSEIF(K.GE.3   .AND.I.LE.IE(K-2)-4.AND.   ! N=0
     +          K.LE.K2-2.AND.I.GE.IS(K+2)  .AND.   ! M=4
     +          I.GE.1.  .AND.K.LE.KE(I  )-2.AND.   ! P=2
     +          I.LE.I2-4.AND.K.GE.KS(I+4)+2)THEN   ! Q=2
                   G(I,K)=H6(I,K) 
         ELSEIF(K.GE.3   .AND.I.LE.IE(K-2)  .AND.   ! N=4
     +          K.LE.K2-2.AND.I.GE.IS(K+2)+4.AND.   ! M=0
     +          I.GE.5.  .AND.K.LE.KE(I-4)-2.AND.   ! P=2
     +          I.LE.I2  .AND.K.GE.KS(I  )+2)THEN   ! Q=2
                   G(I,K)=H9(I,K)
         ELSE
                   G(I,K)=H0(I,K)
         ENDIF
   10    CONTINUE
   20 CONTINUE
C
C
      DO 60 I=I1,I2
         DO 50 K=KS(I),KE(I)
   50    F(I,K)=G(I,K)
   60 CONTINUE
C
      RETURN
C
      END
C
C*******************************************************************
      SUBROUTINE CHNTP
C
C  INTERPOLATES TO MESH WITH CONSTANT INTERVALS IN LOG10(R)
C  WHERE R IS THE OPAL VARIABLE RHO/T6**3.
C
C  ADAPTED FROM OPFIT.FOR
C
C     NEW DENSITY INDEX IS J. 
C     FOR EACH I, J=JS(I) TO JE(I) IN STEPS OF 1.
C     LOG10(R)=0.5*(J+J0)
C     NEW F(I,J)=LOG10(OPACITY) INTERPOLATED TO (I,J) MESH.
C
      PARAMETER(IPI=100,IPK=60,IPJ=60)
      COMMON/COP1/KS(IPI),KE(IPI),K0,K2,A(IPI,IPK)
      COMMON/COP2/IS(IPK),IE(IPK),JS(IPI),JE(IPI),I0,I2,J0,CT,CN,
     + F(IPI,IPK),FX(IPI,IPK),FY(IPI,IPK),FXY(IPI,IPK)
C
      DIMENSION AA(IPJ),B(IPJ),C(IPJ)
C
      L(I,K)=LINT(2.*A(I,K)+0.5)
C
      J0=L(I2,KS(I2))-1
      JS(I2)=1
      DO 10 I=I2-1,1,-1
         JS(I)=MAX(JS(I+1),L(I,KS(I))-J0)
   10 CONTINUE
      JE(1)=L(1,KE(1))-J0
      DO 20 I=2,I2
         JE(I)=MIN(JE(I-1),L(I,KE(I))-J0)
   20 CONTINUE
C
      DO 50 I=1,I2
         N=0
         DO 30 K=KS(I),KE(I)
            N=N+1
            AA(N)=A(I,K)
            B(N)=F(I,K)
   30    CONTINUE
         CALL SPLINE(AA,B,N,C)
         DO 40 J=JS(I),JE(I)
            RL=.5*(J+J0)
            CALL SPLINT(AA,B,N,C,RL,F(I,J))
   40    CONTINUE
   50 CONTINUE
C
      IA=1
      JSA=JS(1)
      JEA=JE(1)
      DO 110 I=2,I2
         IF(JS(I).NE.JSA.OR.JE(I).NE.JEA)THEN
            IA=I
            JSA=JS(I)
            JEA=JE(I)
         ENDIF
  110 CONTINUE
C
      RETURN
C
C
      END
C
C*****************************************************************
C
      FUNCTION LINT(F)
C
C  CONVERSION OF FLOATING POINT TO INTEGER WITH DOWNWARD
C  TRUNCATION FOR F.LT.0
      LINT=F
      IF(F.LT.0.AND.F.NE.LINT)LINT=LINT-1
      RETURN
      END
C*******************************************************************
      SUBROUTINE SPLINE(X,Y,N,Y2)
C
C  CODE FOR SECOND DERIVATIVES USED IN SPLINE FIT.
C  ADAPTED FROM CODE GIVEN IN Recipes.
C
      PARAMETER (NMAX=60)
      DIMENSION X(N),Y(N),Y2(N),U(NMAX)
C
C     FIRST DERIVATIVES AT END POINTS USING QUADRATIC FIT
         YP1=((Y(3)-Y(1))*(X(2)-X(1))**2
     +   -(Y(2)-Y(1))*(X(3)-X(1))**2)/
     +   ((X(3)-X(1))*(X(2)-X(1))*(X(2)-X(3)))
         YPN=((Y(N-2)-Y(N))*(X(N-1)-X(N))**2
     +   -(Y(N-1)-Y(N))*(X(N-2)-X(N))**2)/
     +   ((X(N-2)-X(N))*(X(N-1)-X(N))*(X(N-1)-X(N-2)))
C
      Y2(1)=-0.5
      U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
      DO 11 I=2,N-1
        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        P=SIG*Y2(I-1)+2.
        Y2(I)=(SIG-1.)/P
        U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
     *      /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
11    CONTINUE
      QN=0.5
      UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.)
      DO 12 K=N-1,1,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
12    CONTINUE
      RETURN
      END
C*********************************************************************
      SUBROUTINE SPLINT(XA,YA,N,Y2A,X,Y)
C
C  CODE USING SPLINE FITS, ADAPTED FROM THAT GIVEN IN
C  Recipes.
C
      DIMENSION XA(N),YA(N),Y2A(N)
      KLO=1
      KHI=N
1     IF (KHI-KLO.GT.1) THEN
        K=(KHI+KLO)/2
        IF(XA(K).GT.X)THEN
          KHI=K
        ELSE
          KLO=K
        ENDIF
      GOTO 1
      ENDIF
      H=XA(KHI)-XA(KLO)
      IF (H.EQ.0.) PAUSE 'Bad XA input.'
      A=(XA(KHI)-X)/H
      B=(X-XA(KLO))/H
      Y=A*YA(KLO)+B*YA(KHI)+
     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.
      RETURN
      END
C*********************************************************************
C
      SUBROUTINE DERIV
C
                                                                        
      PARAMETER (IPI=100,IPK=60,IPN=30,IPJ=60)
      COMMON/COP2/IS(IPK),IE(IPK),JS(IPI),JE(IPI),I0,I2,J0,CT,CN,
     + F(IPI,IPK),FX(IPI,IPK),FY(IPI,IPK),FXY(IPI,IPK)
      DIMENSION C(IPI),IJS(IPK),IJE(IPK)
C
C  GET IJS
      JA=JE(1)+1
      DO 20 I=1,I2
         IF(JS(I).LT.JA)THEN
            JB=JA-1
            JA=JS(I)
            DO 10 J=JA,JB
   10       IJS(J)=I
         ENDIF
   20 CONTINUE
C  GET IJE
      JB=JS(I2)-1
      DO 40 I=I2,1,-1
         IF(JE(I).GT.JB)THEN
            JA=JB+1
            JB=JE(I)
            DO 30 J=JA,JB
   30       IJE(J)=I
         ENDIF
   40 CONTINUE
C
C  GET FX
      DO 70 J=JS(I2),JE(1)
         L=0
         DO 50 I=IJS(J),IJE(J)
            L=L+1
            C(L)=F(I,J)
   50    CONTINUE
         CALL GET(C,L)
         L=0
         DO 60 I=IJS(J),IJE(J)
            L=L+1
            FX(I,J)=C(L)
   60    CONTINUE
   70 CONTINUE
C
C  GET FY
      DO 100 I=1,I2
         L=0
         DO 80 J=JS(I),JE(I)
            L=L+1
            C(L)=F(I,J)
   80    CONTINUE
         CALL GET(C,L)
         L=0
         DO 90 J=JS(I),JE(I)
            L=L+1
            FY(I,J)=C(L)
   90    CONTINUE
  100 CONTINUE
C
C  GET FXY
      DO 130 I=1,I2
         L=0
         DO 110 J=JS(I),JE(I)
            L=L+1
            C(L)=FX(I,J)
  110    CONTINUE
         CALL GET(C,L)
         L=0
         DO 120 J=JS(I),JE(I)
            L=L+1
            FXY(I,J)=C(L)
  120    CONTINUE
  130 CONTINUE
C         
      DO 140 I=1,I2
         F  (I,JE(I)+1)=0.
         FX (I,JE(I)+1)=0.
         FY (I,JE(I)+1)=0.
         FXY(I,JE(I)+1)=0.
  140 CONTINUE
      DO 150 J=JS(I2),JE(1)
         F  (IJE(J)+1,J)=0.
         FX (IJE(J)+1,J)=0.
         FY (IJE(J)+1,J)=0.
         FXY(IJE(J)+1,J)=0.
  150 CONTINUE
C
      RETURN
C
      END
C
C******************************************************************
C
      SUBROUTINE GET(F,N)
C
C  SIMPLIFIED CODE FOR SPLINE COEFFICIENTS, FOR CASE OF INTERVALS
C  OF UNITY.
C  RETURNS DERIVATIVES OF ORIGINAL F IN LOCATION F
C
C
      PARAMETER (IPI=100)
      DIMENSION F(IPI),D(IPI),T(IPI)
C
      FP1=(-11.*F(1)+18.*F(2)-9.*F(3)+2.*F(4))/6.
      FPN=(11.*F(N)-18.*F(N-1)+9.*F(N-2)-2.*F(N-3))/6.
C
      D(1)=-.5
      T(1)=.5*(-F(1)+F(2)-FP1)
C
      DO 10 J=2,N-1
         D(J)=-1./(4.+D(J-1))
         T(J)=-D(J)*(F(J-1)-2.*F(J)+F(J+1)-T(J-1))
   10 CONTINUE
C
      D(N)=(FPN+F(N-1)-F(N)-T(N-1))/(2.+D(N-1))
C
      DO 20 J=N-1,1,-1
         D(J)=D(J)*D(J+1)+T(J)
   20 CONTINUE
C
      DO 30 J=2,N-1
         F(J)=-F(J)+F(J+1)-2.*D(J)-D(J+1)
   30 CONTINUE
      F(1)=FP1
      F(N)=FPN
C
      RETURN
      END
C
C*********************************************************************
C
      SUBROUTINE INTRP2(FLT,FLRHO,G,IERR)
C      
      PARAMETER (IPI=100,IPK=60,IPN=30,IPJ=60)
      COMMON/COP2/IS(IPK),IE(IPK),JS(IPI),JE(IPI),I0,I2,J0,CT,CN,
     + F(IPI,IPK),FX(IPI,IPK),FY(IPI,IPK),FXY(IPI,IPK)
C
      DIMENSION B(16)
C
C  FUNCTION DEFINITIONS FOR CUBIC EXPANSION
C
      FF(S,T)=    B( 1)+T*(B( 2)+T*(B( 3)+T*B( 4)))
     +   +S*(     B( 5)+T*(B( 6)+T*(B( 7)+T*B( 8)))
     +   +S*(     B( 9)+T*(B(10)+T*(B(11)+T*B(12)))
     +   +S*(     B(13)+T*(B(14)+T*(B(15)+T*B(16))) )))
C
C
      FLR=FLRHO+18.-3.*FLT
      X=FLT/CT-I0
      Y=2.*FLR-J0
      I=X+1.E-5
      IF(ABS(X-I).LE.1.E-5)X=I
      J=Y+1.E-5
      IF(ABS(Y-J).LE.1.E-5)Y=J
C
      IF(X.GE.1.AND.X.LE.I2)THEN
         IF(Y.LT.JS(I).OR.Y.GT.JE(I))GOTO 2000
      ELSE
         GOTO 1000
      ENDIF
C
C     INTERPOLATE
C
C  GIVEN FUNCTIONS AND DERIVATIVES AT GRID POINTS, COMPUTE COEFFICIENTS.
      B(1)=F(I,J)
      B(2)=FY(I,J)
      B(3)=3*(-F(I,J)+F(I,J+1))-2*FY(I,J)-FY(I,J+1)
      B(4)=2*(F(I,J)-F(I,J+1))+FY(I,J)+FY(I,J+1)
C
      B(5)=FX(I,J)
      B(6)=FXY(I,J)
      B(7)=3*(-FX(I,J)+FX(I,J+1))-2*FXY(I,J)-FXY(I,J+1)
      B(8)=2*(FX(I,J)-FX(I,J+1))+FXY(I,J)+FXY(I,J+1)
C
      B(9)=3*(-F(I,J)+F(I+1,J))-2*FX(I,J)-FX(I+1,J)
      B(10)=3*(-FY(I,J)+FY(I+1,J))-2*FXY(I,J)-FXY(I+1,J)
      B(11)=9*(F(I,J)-F(I+1,J)+F(I+1,J+1)-F(I,J+1))
     + +6*(FX(I,J)-FX(I,J+1)+FY(I,J)-FY(I+1,J))
     + +4*FXY(I,J)
     + +3*(FX(I+1,J)-FX(I+1,J+1)-FY(I+1,J+1)+FY(I,J+1))
     + +2*(FXY(I,J+1)+FXY(I+1,J))
     + +FXY(I+1,J+1)
      B(12)=6*(-F(I,J)+F(I+1,J)-F(I+1,J+1)+F(I,J+1))
     + +4*(-FX(I,J)+FX(I,J+1))
     + +3*(-FY(I,J)+FY(I+1,J)+FY(I+1,J+1)-FY(I,J+1))
     + +2*(-FX(I+1,J)+FX(I+1,J+1)-FXY(I,J)-FXY(I,J+1))
     + -FXY(I+1,J)-FXY(I+1,J+1)
C
      B(13)=2*(F(I,J)-F(I+1,J))+FX(I,J)+FX(I+1,J)
      B(14)=2*(FY(I,J)-FY(I+1,J))+FXY(I,J)+FXY(I+1,J)
      B(15)=6*(-F(I,J)+F(I+1,J)-F(I+1,J+1)+F(I,J+1))
     + +4*(-FY(I,J)+FY(I+1,J))
     + +3*(-FX(I,J)-FX(I+1,J)+FX(I+1,J+1)+FX(I,J+1))
     + +2*(FY(I+1,J+1)-FY(I,J+1)-FXY(I,J)-FXY(I+1,J))
     + -FXY(I+1,J+1)-FXY(I,J+1)
      B(16)=4*(F(I,J)-F(I+1,J)+F(I+1,J+1)-F(I,J+1))
     + +2*(FX(I,J)+FX(I+1,J)-FX(I+1,J+1)-FX(I,J+1)
     +    +FY(I,J)-FY(I+1,J)-FY(I+1,J+1)+FY(I,J+1))
     + +FXY(I,J)+FXY(I+1,J)+FXY(I+1,J+1)+FXY(I,J+1)
C
C  GET G=LOG10(ROSS)
      U=X-I
      V=Y-J
      G=FF(U,V)
      RETURN
C
C     OUT OF T RANGE
 1000 WRITE(6,600)FLT
      GOTO 3000
C     RHO OUT OF RANGE
 2000 WRITE(6,601)FLT,FLRHO,FLR
 3000 IERR=1
      G=9.999
      RETURN
C
  600 FORMAT(2X,'LOG(T)=',1P,E9.2,' IS OUT OF RANGE')
  601 FORMAT(2X,'FOR LOG(T)=',1P,E9.2,', LOG(RHO)=',E9.2,
     + ' IS OUT OF RANGE (LOG(R)=',E9.2,')')
C
      END
C
C******************************************************************   
	SUBROUTINE  CHCKI(I,IT,I1,I3)
C
	IF(IT.NE.(I1+(I-1)*I3))THEN
	   WRITE(6,6000)I,IT,I1,I3
	   STOP
	ENDIF
C
 6000   FORMAT(/2X,'INDEXING ERROR IN LOOP 30 OF READ1'/
     +  2X,'I=',I3/2X,'IT=',I3/2X,'I1=',I3/2X,'I3=',I3/
     +  2X,'SHOULD HAVE IT=I1+(I-1)*I3'/)
C
	RETURN
	END
C
C****************************************************************
	SUBROUTINE CHCKJ(K,JN,JN1,I,KS,JN3)
C
	IF(JN.NE.(JN1+(K-KS)*JN3))THEN
	   WRITE(6,6000)I,K,JN,JN1,KS,JN3
	   STOP
	ENDIF
C
 6000   FORMAT(/2X,'INDEXING ERROR LOOP 20 OF  READ1, I=',I3/
     +  2X,'K=',I3/2X,'JN=',I3/2X,'KS(I)=',I3/2X,'JN3=',I3/
     +  2X,'SHOULD HAVE JN=JN1+(K-KS(I))*JN3'/)
C
	RETURN
	END
C
C******************************************
	
