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