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