C BLF IS A PROGRAM DESIGNED TO MODEL SEVERAL PETROLOGIC PROCESSES C AS THEY MIGHT OCCUR IN A DRY UPPER CRUSTAL MAGMA CHAMBER. THESE C PROCESSES INCLUDE IN SITU FRACTIONAL AND EQUILIBRIUM XLIZATION, C MAGMA CHAMBER RECHARGE, CRUSTAL ASSIMILATION, AND PERIODIC TAPPING C OF THE MAGMA CHAMBER INCLUDING SYSTEMS OPEN OR CLOSED TO OXYGEN. C THIS PROGRAM WAS WRITTEN BY ROGER L.NIELSEN, COLLEGE OF OCEANOGRAPHY, C OCEAN. ADMIN. 104, OREGON STATE UNIV, CORVALLIS, OREGON 97731-5503 C phone 503-737-3023. C THIS VERSION OF BLF INCORPORATES TEMPERATURE DEPENDENT, C COMPOSITIONALLY INDEPENDENT, SINGLE COMPONENT DISTRIBUTION COEFFICIENTS C FROM EXPERIMENTAL DATA. THE PROGRAM ALSO HAS THE OPTION C OF RANDOMIZING THE MIXING AND SAMPLING PARAMETERS ABOUT A GAUSIAN MEAN C VALUE. PERIODICITY IS INDEPENDENT OF ASSIMIATION C C C VERSION 6/6/91 C program BLF IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NCOMP=10) DIMENSION XL(10),UL(10),DULDT(10),DULDX(10),TP(4) DIMENSION F(60),CUM(60),G(14,50),X(14,50),W(15,50),WT(50),DK(50) DIMENSION T(60),SLOPE(9,60),YCEPT(9,60),FRAC(60),CPXTR(60) DIMENSION OLTR(60),PLTR(60),SPTR(60),OPXTR(60),ILMTR(60),PIGTR(60) DIMENSION XCAO(4),XFEO(4),XMGO(4),H(100),DZ(100) DIMENSION RDM(20),RD(20),DND(10001),APTR(60) CHARACTER LINE(10)*80,ELEM(60)*10 REAL LR,MNOL,NA,KSP,IFE,IMG,ICR,ILMTR,OEN,OFS,OWO,JD REAL OF,OM,OLNI,OT,MIX,MGPLAG C C INITIALIZE PYX SOLUTION PARAMETERS AND RELATIVE STD STATES C DO 1 I=1,100 H(I)=0. DZ(I)=0. 1 CONTINUE START=0. SIZEMAG=1. MIX=0.0 PCXL=0. PCXLSL=0. R=8.3143 DZ(1)=1.D0 DZ(2)=1.D-12 DZ(3)=2.D-04 DZ(4)=4.D-02 DZ(43)= 31216. DZ(44)=-.00610 DZ(45)=25484. DZ(46)=.08120 DZ(47)=20697. DZ(48)=-.00235 DZ(49)=16941. DZ(50)=.00592 DZ(16)=-.09557 DZ(17)=-.07903 DZ(18)=.08807 DZ(19)=.04559 DZ(21)=36.641121 DZ(22)=-.089 DZ(20)=24.028327 H(35)=20000. H(36)=80000. H(43)=-38140. H(45)=9900. H(47)=0. H(46)=-9450. H(39)=-23000. H(32)=17000. C THESE ARE RELATIVE TO OPX!! C C FOLLOWING ARE REVISED, NIELSEN FEB,86 H(1)=-150172.9 H(2)=-52039.2 H(3)=-68094.1 H(4)=-18466.06 H(5)=-153733.9 H(6)=-53882.2 H(7)=-62833.1 H(8)=15085.9 H(10)=-126394.0 H(11)=-105216.0 XMGO(1)=0.54/2. XMGO(2)=0.95/2. XMGO(3)=0.90/2. XCAO(1)=0.45/2. XCAO(2)=0.04/2. XCAO(3)=0.09/2. DZ(5)=2.49429 DZ(10)=83.30929 DZ(11)= 33.5897 DZ(12)=39.57607 DZ(6)=15.11971 DZ(23)= 85.21929 DZ(24)= 35.26577 DZ(25)= 40.10607 DZ(30)= 65.68297 DZ(29)= 57.4697 H(48)=-109810. C DO 3 I=1,3 TP(I)=1450 3 CONTINUE DO 10 I=1,49 F(I)=0. CUM(I)=0. T(I)=1350. DO 5 J=1,14 RDM(J)=0. RD(J)=0. G(J,I)=0. W(J,I)=0. X(J,I)=0. 5 CONTINUE 10 CONTINUE LR=0 SYTOT=1. C1=0. C=0. ST=0. Z=0. START=0. C INPUT FILE OF RANDOM # OPEN(UNIT=1,FILE='RANDOM.DAT',STATUS='OLD') DO 7 I=1,9999 READ(1,*) DND(I) 7 CONTINUE CLOSE (UNIT=1) C C OPEN FILE FOR INPUT DATA C C OPEN(UNIT=1,FILE='BLF.DAT',STATUS='OLD') OPEN(UNIT=2,FILE='BLFPLOT',STATUS='unknown') OPEN(UNIT=3,FILE='newds3.dat',STATUS='old') C C INPUT TITLE - LINE(1) C TOTAL # OF ELEMENTS - NT C STARTING COMPOSITION - G(1,I) C ELEMENT NAMES - ELEM(I) C MOLECULAR WEIGHTS - WT(I) C C READ(1,*)LINE(1) READ(1,*)LINE(2) READ(3,*)LINE(2) READ(1,*) NT READ(1,*)LINE(2) READ(3,*) NT READ(3,*)LINE(2) DO 20 I=1,13 READ(1,*) ELEM(I),G(1,I) READ(3,*) ELEM(I),WT(I) 20 CONTINUE C INPUT TRACE ELEMENTS AND THEIR ATOMIC WEIGHTS C G(1,I) I=14-NT WT(I) I=14-NT C C INPUT SLOPES AND Y INTERCEPTS FOR TRACE ELEMENT PARTITIONING C C SLOPE(J,I) YCEPT(J,I) C READ(1,*)LINE(2) READ(3,*)LINE(2) READ(3,*)LINE(2) DO 30 I=14,NT READ(1,*) ELEM(I),G(1,I) READ(3,*) ELEM(I),WT(I) DO 25 J=1,8 READ(3,*) SLOPE(J,I),YCEPT(J,I) 25 CONTINUE 30 CONTINUE C C INPUT ASSIMILANT COMPOSITION - G(3,I) C C READ(1,*)LINE(2) DO 35 I=1,NT READ(1,*) ELEM(I),G(3,I) 35 CONTINUE C C INPUT FRACTIONATION FACTORS FOR OLIVINE, OPX, CPX, PLAG, C SPINEL ILMENTE AND ILMENITE. A VALUE OF 1 REPRESENTS PERFECT FRAC- C TIONAL CRYSTALLIZATION. 0. MODELS EQUILIBRIUM CRYSTALLIZATION. C C C INPUT RECHARGE, ASSIMILATION,ERUPTION FACTORS C RECHG ASSIM ERUPT C C THESE FACTORS ARE REFERENCED TO THE CRYSTALLIZATION INCREMENT C A VALUE OF 1 REPRESENTS A RECHARGE ETC. OF AN AMOUNT EQUAL C TO THAT REMOVED FROM THE LIQUID SYSTEM BY CRYSTALLIZATION C C INPUT THE OXYGEN FUGACITY - THE FUGACITY IS SET TO PARALLEL C THE QFM CURVE. ENTER THE # OF NATURAL LOG UNITS ABOVE OR C BELOW THE BUFFER YOU WISH TO SET THE FO2 C C INPUT THE MIXING PROCESSES YOU WISH TO MODEL C THE CHOICES ARE EITHER PERIODIC OR CONTINUOUS C CONTINUOUS ASSUMES MIXING AND FRACTIONATION ARE OCCURING C SIMULTANEOUSLY. PERIODIC MIXING ASSUMES THAT EPISODES C OF MIXING ARE SEPARATED BY PERIODS OF FRACTIONATION. C C THE PERIODS BETWEEN MIXING EPISODES IS SET BY THE VARIABLE C (PER). THIS SHOULD BE SET TO 0. WHEN CONTINUOUS MIXING IS C BEING MODELLED. C READ(1,*) LINE(2) READ(1,*) (FRAC(I),I=1,8) READ(1,*) LINE(2),RECHGI,ASSIMI,ERUPTI READ(1,*) LINE(2),DFO2 READ(1,*) LINE(2),PERI RECHG=RECHGI ASSIM=ASSIMI ERUPT=ERUPTI PER=PERI C C INPUT HOW OFTEN YOU WISH PERIODICITY OF RECHARGE C C INPUT THE CRYSTALLIZATION INCREMENT -C2X C SUGGESTED VALUE IS 0.2 % - ENTER AS % C WHEN MODELLING COMPATIBLE ELEMENTS A SMALLER VALUE C E.G. 0.01-0.05% IS MORE APPROPRIATE. C C READ(1,*) LINE(2),UTI READ(1,*) LINE(2),C2X C C READ IN RANDOMIZATION PARAMETERS C DO 37 I=1,5 READ(1,*) LINE(2),RD(I) 37 CONTINUE UT=UTI READ(1,*) LINE(2),IC READ(1,*) LINE(2),SOLMAG CLOSE (UNIT=1) CLOSE(UNIT=3) C C OUTPUT STARTING CONDITIONS TO FILE C C OPEN(UNIT=3,FILE='BLF.out',STATUS='unknown') WRITE(3,40) LINE(1) 40 FORMAT(4X,A,/) WRITE(2,41) LINE(1) 41 FORMAT(A,/,I2) WRITE(3,45) 45 FORMAT(4X,'STARTING COMPOSITION ASSIMILANT',/) WRITE(2,*) 'PCTFR',CHAR(09),'TEMP',CHAR(09) WRITE(2,*) CHAR(09),CHAR(09),(ELEM(I),CHAR(09),I=1,NT) DO 52 I=1,NT WRITE(3,50)ELEM(I),G(1,I),G(3,I) 50 FORMAT(A,4X,F7.2,10X,F7.2) 52 CONTINUE WRITE(3,55) 55 FORMAT(' OLIVINE OPX CPX PLAG ', :' SPINEL ILMENITE PIG APATITE',/) WRITE(3,60) (FRAC(I),I=1,8) 60 FORMAT(8F9.2) WRITE(3,65) RECHG,ASSIM,ERUPT 65 FORMAT(/,'RECHARGE ',F5.1,/,'ASSIMILATION ',F5.1,/, :'ERUPTION ',F5.1,/) WRITE(3,70) DFO2 70 FORMAT(/,'FUGACITY OF O2 IS SET',F5.1,' LOG UNITS FROM QFM',/) UT=UT/100. PER=PER/100. C2X=C2X/100. UTI=UTI/100. PERI=PERI/100. WRITE(3,71)(RD(I),I=1,5) 71 FORMAT(' PERIOD BLF % XLS RECH ASSIM :. ERUPT ',/,2F7.2,' ',3F7.2) WRITE(3,73) UT,PER,C2X 73 FORMAT(' RECHARGE INTERVAL ',F10.4,' FRACTION ',/,' PERIODIC : MIX ',F10.4,/,' CRYSTALLIZATION INCREMENT ',F10.4,/) DO 1274 I=14,NT WRITE(3,1271) ELEM(I) 1271 FORMAT(A,' OLIVINE OPX CPX PLAG SPINEL ', :' ILMEN PIG APATITE') WRITE(3,1272)(SLOPE(J,I),J=1,8) 1272 FORMAT('SLOPE',8F9.2) WRITE(3,1273)(YCEPT(J,I),J=1,8) 1273 FORMAT('YCEPT',8F9.2,/) 1274 CONTINUE TG=0. TA=.000001 G(1,8)=G(1,8)/10000. G(3,8)=G(3,8)/10000. G(1,11)=G(1,11)/10000 G(3,11)=G(3,11)/10000 DO 75 I=14,NT G(1,I)=G(1,I)/10000 G(3,I)=G(3,I)/10000 75 CONTINUE UT=UTI PER=PERI IF(RD(1).EQ.1.) UT=UTI*DND(IC)*2 IF (RD(2).EQ.1.) PER=PERI*DND(IC+36)*2 WRITE(*,76) UT,PER 76 FORMAT(' periodicity of rech and in situ ',2F12.5,/) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CONVERSION TO MOLE FRACTION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 80 I=1,NT G(1,I)=G(1,I)/WT(I) G(3,I)=G(3,I)/WT(I) TG=TG+G(1,I) TA=TA+G(3,I) 80 CONTINUE DO 85 I=1,NT G(1,I)=G(1,I)/TG G(3,I)=G(3,I)/TA G(4,I)=G(1,I) G(5,I)=G(1,I) G(6,I)=G(1,I) 85 CONTINUE RTIO=DLOG((G(1,10)/G(1,9))/2.D0) DO 95 I=1,8 DO 90 K=14,NT IF (YCEPT(I,K).EQ.0.) YCEPT(I,K)=-15. 90 CONTINUE 95 CONTINUE C C BEGINNING OF MAIN ITERATION LOOP C C 100 IF (PER.EQ.0.) GO TO 103 PRINT *,' ' PRINT *,'MIX',MIX,'PER',PER PRINT *,' ' IF (MIX.GE.PER) GO TO 9000 IF (PER.GT.0.) GO TO 110 103 SUM=0. DO 105 I=1,NT SUM=SUM+G(4,I) 105 CONTINUE DO 107 I=1,NT X(10,I)=0. G(4,I)=G(4,I)-(ERUPT*C*(G(4,I)/SUM))+C*(RECHG*G(5,I)+ASSIM*G(3,I)) 107 CONTINUE 110 SYTOT=0. DO 120 I=1,NT SYTOT=SYTOT+G(4,I) X(10,I)=0. G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 120 CONTINUE XLTOT=0. DO 125 I=1,8 XLTOT=XLTOT+F(I) 125 CONTINUE C=C2X C1=C IF(DFO2.LE.-5.) GO TO 127 FETOT=G(1,10)+G(1,9) G(1,9)=FETOT/(2.*EXP(RTIO)+1.) G(1,10)=FETOT-G(1,9) 127 VL=51.-50.*(G(1,9)/(G(1,9)+G(1,10))) C PRINT *,VL,'VL',G(1,9),'FE+2',G(1,10),'FE+3' C PRINT *,RTIO,FETOT IF(VL.LT.5.) GO TO 130 VL=5. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CALCULATION OF MELT COMPONENTS C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 130 Z=0. DO 132 I=1,NT Z=Z+G(1,I) 132 CONTINUE ALUM=G(1,3)/Z TITAN=G(1,7)/Z SILICA=(G(1,4)*1.1)/Z SIO2=G(1,4)/Z PHOS=G(1,12)/Z SODA=G(1,1)/Z CLC=G(1,6)/Z POTA=G(1,5)/Z ZRLIQ=G(1,21)/Z BLIQ=G(1,22)/Z FLIQ=G(1,38)/Z TALIQ=G(1,39)/Z C WRITE(*,*) Z, G(1,21),G(1,22) C WRITE(*,*) ZRLIQ,BLIQ,FLIQ,TALIQ Y=G(1,1)+G(1,4)+G(1,5) S=G(1,2)+G(1,3)+G(1,6)+G(1,7)+G(1,8)+G(1,9)+G(1,10)+G(1,11)+ :G(1,12)+G(1,13)-G(1,1)-G(1,5) G(1,3)=(G(1,3)-G(1,1)-G(1,5))/S G(1,2)=G(1,2)/S G(1,6)=G(1,6)/S G(1,1)=G(1,1)/Y G(1,4)=G(1,4)/Y G(1,5)=G(1,5)/Y G(1,7)=G(1,7)/S G(1,8)=G(1,8)/S G(1,9)=G(1,9)/S G(1,10)=G(1,10)/S G(1,11)=G(1,11)/S G(1,12)=G(1,12)/S G(1,13)=G(1,13)/S DO 140 I=14,NT G(1,I)=G(1,I)/S 140 CONTINUE C PRINT *,'BEFORE OLIVINE' CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C OLIVINE LIQUIDUS TEMPERATURE CALCULATION SECTION C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC D2=1000. D=0. ITER=0 D2=D2/2. C PRINT *,'BEFORE 150','TEMPERATURE ',T(1) 150 OF=EXP(5655.5/T(1)-4.29)*G(1,9)*G(1,4)**0.5 CC PRINT *,'AFTER 150' OM=EXP(5578.0/T(1)-2.96)*G(1,2)*G(1,4)**0.5 C PRINT *,'AFTER OM',OM MNOL=EXP(8655./T(1)-6.5)*G(1,13) C PRINT *,'AFTER MNOL' OLNI=EXP(11230./T(1)-5.7)*G(1,11)*G(1,4)**0.5 OT=0.018*G(1,7) RC=EXP(9242./T(1)-5.47)*G(1,8)/5. CC=0.002 ITER=ITER+1 IF (ITER.GT.400) go to 9999 D=2./3.-OF-OM-OLNI-RC-OT-CC-MNOL T(1)=T(1)-D2*D IF(ABS(D).GT..0002) GO TO 150 DO 160 I=14,NT OLTR(I)=EXP(SLOPE(1,I)/T(1)+YCEPT(1,I))*G(1,I) 160 CONTINUE WRITE(*,165) OM*300./2.,T(1)-273. 165 FORMAT (/,' OLIVINE ',F7.2,5X,F7.1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C PYROXENE LIQUIDUS CALCULATION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC XL(1)=G(1,2) XL(3)=G(1,6) XL(2)=G(1,9) XL(4)=G(1,4) DO 168 INDEXA=1,3 H(30)=TP(INDEXA) H(22)=XMGO(INDEXA)*2. H(23)=XCAO(INDEXA)*2. CALL PXCALC(INDEXA,XL,UL,DULDT,DULDX,TP,H,DZ,R) TP(INDEXA)=H(30) XMGO(INDEXA)=H(22)/2. XCAO(INDEXA)=H(23)/2. XFEO(INDEXA)=0.5D0-XMGO(INDEXA)-XCAO(INDEXA) 168 CONTINUE C=C2X IF(XCAO(3).GT.0.10) TP(3)=TP(3)-1.5 IF(XCAO(1).LT.0.10) TP(1)=TP(1)-1.5 IF(TP(1).GT.TP(2).AND.TP(1).GT.TP(3)) GO TO 220 IF(TP(2).GT.TP(1).AND.TP(2).GT.TP(3)) GO TO 170 IF(TP(3).GT.TP(2).AND.TP(3).GT.TP(2)) GO TO 185 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ORTHOPYROXENE MINOR ELEM CALC C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 170 T(2)=TP(2) PA=EXP(30000./T(2)-23.65)*G(1,3) PC=EXP(10526./T(2)-6.3)*G(1,8)*5. PT=(4934./T(2)+105.6881*ALUM**2-23.7*ALUM+0.7897* :(2*CLC+SODA+POTA)/SIO2+0.3799*XCAO(2)-2.4404)*TITAN PMN=EXP(8600/T(2)-6.42)*G(1,13) PN=EXP(22727./T(2)-14.65)*G(1,11) WRITE(*,175) XMGO(2)*200.,XCAO(2)*200.,T(2)-273. 175 FORMAT(' ORTHOPYROX ',F5.1,' EN ',F5.1,' WO ',F6.1,' TEMP') DO 180 I=14,NT OPXTR(I)=EXP(SLOPE(2,I)/T(2)+YCEPT(2,I))*G(1,I) 180 CONTINUE DTI=PT/TITAN T(3)=1300. T(7)=1300. X(10,2)=1. X(10,4)=.5-PA/2. X(10,6)=XCAO(2) X(10,7)=PT X(10,8)=PC X(10,9)=XFEO(2) X(10,11)=PN X(10,13)=PMN DO 545 I=14,NT X(10,I)=OPXTR(I) 545 CONTINUE XMCA=XMGO(2)*XCAO(2) XCAMGL=G(1,6)*G(1,2)*G(1,4) GAF=G(1,3)*G(1,9) X(10,15)=YCEPT(2,15)*((XMCA*GAF*G(1,15))/(XFEO(2)*XCAMGL)) X(10,20)=YCEPT(2,20)*((XMCA*GAF*G(1,20))/(XFEO(2)*XCAMGL)) X(10,25)=YCEPT(2,25)*((XMCA*GAF*G(1,25))/(XFEO(2)*XCAMGL)) X(10,26)=YCEPT(2,26)*((XMCA*GAF*G(1,26))/(XFEO(2)*XCAMGL)) X(10,27)=YCEPT(2,27)*((XMCA*GAF*G(1,27))/(XFEO(2)*XCAMGL)) X(10,28)=YCEPT(2,28)*((XMCA*GAF*G(1,28))/(XFEO(2)*XCAMGL)) X(10,30)=YCEPT(2,30)*((XMCA*GAF*G(1,30))/(XFEO(2)*XCAMGL)) X(10,31)=YCEPT(2,31)*((XMCA*GAF*G(1,31))/(XFEO(2)*XCAMGL)) X(10,32)=YCEPT(2,32)*((XMCA*GAF*G(1,32))/(XFEO(2)*XCAMGL)) X(10,33)=YCEPT(2,33)*((XMCA*GAF*G(1,33))/(XFEO(2)*XCAMGL)) X(10,34)=YCEPT(2,34)*((XMCA*GAF*G(1,34))/(XFEO(2)*XCAMGL)) X(10,35)=YCEPT(2,35)*((XMCA*GAF*G(1,35))/(XFEO(2)*XCAMGL)) X(10,36)=YCEPT(2,36)*((XMCA*GAF*G(1,36))/(XFEO(2)*XCAMGL)) X(10,37)=YCEPT(2,37)*((XMCA*GAF*G(1,37))/(XFEO(2)*XCAMGL)) X(10,21)=(SLOPE(2,21)*DTI+YCEPT(2,21))*ZRLIQ X(10,22)=(SLOPE(2,22)*DTI+YCEPT(2,22))*BLIQ X(10,38)=(SLOPE(2,38)*DTI+YCEPT(2,38))*FLIQ X(10,39)=(SLOPE(2,39)*DTI+YCEPT(2,39))*TALIQ DO 547 I=3,13 X(10,2)=X(10,2)-X(10,I) 547 CONTINUE GO TO 259 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C PIGEONITE MINOR ELEMENT CALC C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 185 T(2)=TP(3) PA=EXP(30000./T(2)-23.65)*G(1,3) PC=EXP(10526./T(2)-6.3)*G(1,8)*5. PT=(4934./T(2)+105.6881*ALUM**2-23.7*ALUM+0.7897* :(2*CLC+SODA+POTA)/SIO2+0.3799*XCAO(3)-2.4404)*TITAN PMN=EXP(8600/T(2)-6.42)*G(1,13) PN=EXP(22727./T(2)-14.65)*G(1,11) JD=0. WRITE(*,205)XMGO(3)*200.,XCAO(3)*200.,T(2)-273. 205 FORMAT(' PIGEONITE ',F5.1,' EN ',F5.1,' WO ',F6.1,' TEMP') T(7)=T(2) T(2)=1300. T(3)=1300. DO 210 I=14,NT PIGTR(I)=EXP(SLOPE(7,I)/T(7)+YCEPT(7,I))*G(1,I) 210 CONTINUE DTI=PT/TITAN X(10,2)=1.-JD X(10,4)=.5-PA/2. X(10,6)=XCAO(3) X(10,7)=PT X(10,8)=PC/2. X(10,9)=XFEO(3) X(10,11)=PN X(10,13)=PMN X(10,1)=JD X(10,3)=PA DO 634 I=14,NT X(10,I)=PIGTR(I) 634 CONTINUE XMCA=XMGO(3)*XCAO(3) XCAMGL=G(1,6)*G(1,2)*G(1,4) GAF=G(1,3)*G(1,9) X(10,15)=YCEPT(7,15)*((XMCA*GAF*G(1,15))/(XFEO(3)*XCAMGL)) X(10,20)=YCEPT(7,20)*((XMCA*GAF*G(1,20))/(XFEO(3)*XCAMGL)) X(10,25)=YCEPT(7,25)*((XMCA*GAF*G(1,25))/(XFEO(3)*XCAMGL)) X(10,26)=YCEPT(7,26)*((XMCA*GAF*G(1,26))/(XFEO(3)*XCAMGL)) X(10,27)=YCEPT(7,27)*((XMCA*GAF*G(1,27))/(XFEO(3)*XCAMGL)) X(10,28)=YCEPT(7,28)*((XMCA*GAF*G(1,28))/(XFEO(3)*XCAMGL)) X(10,30)=YCEPT(7,30)*((XMCA*GAF*G(1,30))/(XFEO(3)*XCAMGL)) X(10,31)=YCEPT(7,31)*((XMCA*GAF*G(1,31))/(XFEO(3)*XCAMGL)) X(10,32)=YCEPT(7,32)*((XMCA*GAF*G(1,32))/(XFEO(3)*XCAMGL)) X(10,33)=YCEPT(7,33)*((XMCA*GAF*G(1,33))/(XFEO(3)*XCAMGL)) X(10,34)=YCEPT(7,34)*((XMCA*GAF*G(1,34))/(XFEO(3)*XCAMGL)) X(10,35)=YCEPT(7,35)*((XMCA*GAF*G(1,35))/(XFEO(3)*XCAMGL)) X(10,36)=YCEPT(7,36)*((XMCA*GAF*G(1,36))/(XFEO(3)*XCAMGL)) X(10,37)=YCEPT(7,37)*((XMCA*GAF*G(1,37))/(XFEO(3)*XCAMGL)) X(10,21)=(SLOPE(7,21)*DTI+YCEPT(7,21))*ZRLIQ X(10,22)=(SLOPE(7,22)*DTI+YCEPT(7,22))*BLIQ X(10,38)=(SLOPE(7,38)*DTI+YCEPT(7,38))*FLIQ X(10,39)=(SLOPE(7,39)*DTI+YCEPT(7,39))*TALIQ DO 635 I=3,13 X(10,2)=X(10,2)-X(10,I) 635 CONTINUE GO TO 259 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CLINOPYROXENE LIQUIDUS CALCULATION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 220 T(2)=TP(1) CR=EXP(36000./T(2)-24.7)*G(1,8)*5. CT=(4934./T(2)+105.6881*ALUM**2-23.7*ALUM+0.7897* :(2*CLC+SODA+POTA)/SIO2+0.3799*XCAO(1)-2.4404)*TITAN CMN=EXP(9600./T(2)-7.548)*G(1,13) CA=(502./T(2)+16.104*CT-0.0872*(2*CLC+SODA+POTA) :/SIO2+0.1888*XCAO(1)-0.1937)*ALUM*0.8 CNI=EXP(22727./T(2)-14.65)*G(1,11) JD=EXP(6800./T(2)-7.496)*G(1,1) C PRINT *,' CT, CA', CT, CA WRITE(*,250)XMGO(1)*200,XCAO(1)*200.,T(2)-273. 250 FORMAT(' HI CA PYROXENE',F5.1,' EN ',F5.1,' WO ',F6.1,' TEMP') T(3)=T(2) T(2)=1300. T(7)=1300. DO 255 I=14,NT CPXTR(I)=EXP(SLOPE(3,I)/T(3)+YCEPT(3,I))*G(1,I) 255 CONTINUE DTI=CT/TITAN X(10,2)=1.-JD X(10,4)=.5-2*CT-(CA-2*CT)/2. X(10,6)=XCAO(1) X(10,7)=CT X(10,8)=CR/2. X(10,9)=XFEO(1) X(10,11)=CNI X(10,13)=CMN X(10,1)=JD X(10,3)=CA DO 565 I=14,NT X(10,I)=CPXTR(I) 565 CONTINUE denom=(G(1,6)*G(1,2)*G(1,4)*XFEO(1)) XMCA=(XMGO(1)*XCAO(1)*G(1,3)*G(1,9))/denom C WRITE(*,*) ' XMGO(1) ', ' XCAO(1) ', ' XFEO(1) ' C WRITE(*,*) XMGO(1), XCAO(1), XFEO(1), T(3) X(10,15)=EXP(SLOPE(3,15)/T(3)+YCEPT(3,15))*XMCA*G(1,15) X(10,20)=EXP(SLOPE(3,20)/T(3)+YCEPT(3,20))*XMCA*G(1,20) X(10,25)=EXP(SLOPE(3,25)/T(3)+YCEPT(3,25))*XMCA*G(1,25) X(10,26)=EXP(SLOPE(3,26)/T(3)+YCEPT(3,26))*XMCA*G(1,26) X(10,27)=EXP(SLOPE(3,27)/T(3)+YCEPT(3,27))*XMCA*G(1,27) X(10,28)=EXP(SLOPE(3,28)/T(3)+YCEPT(3,28))*XMCA*G(1,28) X(10,30)=EXP(SLOPE(3,30)/T(3)+YCEPT(3,30))*XMCA*G(1,30) X(10,31)=EXP(SLOPE(3,31)/T(3)+YCEPT(3,31))*XMCA*G(1,31) X(10,32)=EXP(SLOPE(3,32)/T(3)+YCEPT(3,32))*XMCA*G(1,32) X(10,33)=EXP(SLOPE(3,33)/T(3)+YCEPT(3,33))*XMCA*G(1,33) X(10,34)=EXP(SLOPE(3,34)/T(3)+YCEPT(3,34))*XMCA*G(1,34) X(10,35)=EXP(SLOPE(3,35)/T(3)+YCEPT(3,35))*XMCA*G(1,35) X(10,36)=EXP(SLOPE(3,36)/T(3)+YCEPT(3,36))*XMCA*G(1,36) X(10,37)=EXP(SLOPE(3,37)/T(3)+YCEPT(3,37))*XMCA*G(1,37) X(10,21)=(SLOPE(3,21)*DTI+YCEPT(3,21))*ZRLIQ X(10,22)=(SLOPE(3,22)*DTI+YCEPT(3,22))*BLIQ X(10,38)=(SLOPE(3,38)*DTI+YCEPT(3,38))*FLIQ X(10,39)=(SLOPE(3,39)*DTI+YCEPT(3,39))*TALIQ C WRITE(*,*) ' SLOPE HF ',' YCEPT HF ' C WRITE(*,*) SLOPE(3,38), YCEPT(3,38) C WRITE(*,*) ' ZRLIQ ' ,' BLIQ ',' FLIQ ' ,' TALIQ ' C WRITE(*,*) ZRLIQ,BLIQ,FLIQ,TALIQ C WRITE(*,*) G(1,1), G(1,2) DO 567 I=3,NT X(10,2)=X(10,2)-X(10,I) C WRITE(*,*) I, DTI, X(10,I),G(1,I) 567 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C PLAGIOCLASE LIQUIDUS CALCULATION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 259 DP=500. 260 ALKAL=(2*CLC+SODA+POTA)/SIO2 FRAL=ALUM-SODA-POTA ELT=G(1,6)*G(1,3)**2*G(1,4)**2 AN=EXP(4053./T(4)-2.0551*ALKAL-12.3235* :SIO2-11.9235*FRAL+7.8628)*ELT NA=EXP(8942./T(4)-1.0583*ALKAL+1.1986*SIO2+0.3451* :FRAL-6.3112)*G(1,1)*G(1,4)**3 KSP=.153*G(1,5)*G(1,4)**3 FEPLAG=0.02*G(1,9) MGPLAG=.015*G(1,2) D=.2-NA-AN-KSP-FEPLAG-MGPLAG T(4)=T(4)-DP*D IF(ABS(D).GT..0001) GO TO 260 AN=.2-NA-KSP-FEPLAG-MGPLAG WRITE(*,265)AN*500.,T(4)-273. 265 FORMAT(' PLAGIOCLASE ',F7.2,' ANORTHITE ',F7.1,' TEMP') DO 270 I=14,NT PLTR(I)=EXP(SLOPE(4,I)/T(4)+YCEPT(4,I))*G(1,I) 270 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C SPINEL LIQUIDUS CALCULATION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC TER=2. CNT=0. IF(START.EQ. 0.) TER=10. D=1. 275 SM=(-895./T(5)+.89)*(G(1,2)/(G(1,2)+G(1,9))) SF=1./3.-SM+ST-SN-SMN ST=(EXP(25000./T(5)-18.6)*G(1,7)*G(1,2)**2)/(SM**2) S3=EXP(8800./T(5)-4.95)*G(1,10) IF(T(5).GT.2500.) GO TO 9999 IF(T(5).LT.1448.46797) GO TO 280 SA=EXP(-4000./T(5)+5.2)*ALUM**2 GO TO 285 280 SA=EXP(-59718./T(5)+43.66838)*ALUM**2 285 SC=EXP(11261./T(5)-1.8)*G(1,8) SN=EXP(3970./T(5)-1.37)*G(1,11) SMN=EXP(11521./T(5)-9.17)*G(1,13) CNT=CNT+1. IF(CNT.GT.300.) GO TO 9999 B=D D=2./3.-SA-S3-SC-2*ST IF(D.GT.0.) GO TO 295 IF (B.LT.0.) GO TO 290 TER=TER/2. 290 T(5)=T(5)+TER GO TO 310 295 IF(B.GT.0.) GO TO 300 TER=TER/2. 300 T(5)=T(5)-TER 310 IF(ABS(D).GT..0005) GO TO 275 313 SA=1.-ST-S3-SF-SM-SN-SMN-SC WRITE(*,315) SM,SA,SC,T(5)-273 315 FORMAT(' SPINEL ',' MG ',F6.4,' AL ',F6.4,' CR ',F6.4, :' TEMP ',F6.1) DO 320 I=14,NT SPTR(I)=EXP(SLOPE(5,I)/T(5)+YCEPT(5,I))*G(1,I) 320 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ILMENITE LIQUIDUS CALCULATION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMG=G(1,2) ICR=11.*G(1,8) CNT=0. 325 IFE=EXP(8950./T(6)-3.3)*G(1,7)*G(1,9) D=1.-IMG-IFE-ICR T(6)=T(6)-100*D CNT=CNT+1. IF(CNT.GT.100.) GO TO 9999 IF(ABS(D).GT..001) GO TO 325 IFE=1.-IMG-ICR WRITE(*,330)IFE,T(6)-273. 330 FORMAT(' ILMENITE ',F7.4,' FE ',F6.1,' TEMP ') DO 335 I=14,NT ILMTR(I)=EXP(SLOPE(6,I)/T(6)+YCEPT(6,I))*G(1,I) 335 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C APATITE LIQUIDUS DETERMINATION C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC T(8)=(8400.+((SILICA-.5)*26400.))/ :(DLOG(.385/(1.1*PHOS))+(3.1+(12.4*(SILICA-.5)))) TAP=T(8)-273 PRINT *,' APATITE ',TAP PRINT *,' ' PRINT *,' ' SUMAP=0. DO 9335 I=14,NT APTR(I)=EXP(SLOPE(8,I)/T(8)+YCEPT(8,I))*G(1,I) SUMAP=SUMAP+APTR(I) 9335 CONTINUE TLARG=T(1) DO 336 I=2,8 IF (T(I).LT.TLARG) GO TO 336 TLARG=T(I) 336 CONTINUE IF(T(1).EQ.TLARG) GO TO 520 IF(T(2).EQ.TLARG) GO TO 540 IF(T(3).EQ.TLARG) GO TO 560 IF(T(4).EQ.TLARG) GO TO 580 IF(T(5).EQ.TLARG) GO TO 600 IF(T(6).EQ.TLARG) GO TO 620 IF(T(7).EQ.TLARG) GO TO 632 GO TO 638 340 TTEST=TEMP-3. TDIFF=TEMP-5. DO 345 I=1,NT X(10,I)=0 345 CONTINUE FETOT=G(1,9)+G(1,10) IF(DFO2.LE.-5.)GO TO 347 G(1,9)=FETOT/(2.*EXP(RTIO)+1.) G(1,10)=FETOT-G(1,9) 347 IF (TTEST.GT.T(1)) GO TO 650 GO TO 349 650 IF(F(1).GT.0.) GO TO 640 C C ALGORYTHM FOR THE CALCULATION OF FE+3/+2 RATIO C REGRESSION OF KILINC ET AL.(1983) C 349 Z=0. ZTOT=0. DO 350 I=1,13 G(2,I)=G(1,I) ZTOT=ZTOT+G(1,I) C PRINT *,ELEM(I),G(1,I),G(2,I) 350 CONTINUE G(2,3)=G(2,3)/2. G(2,1)=G(2,1)/2. G(2,5)=G(2,5)/2. DO 355 I=1,13 Z=Z+G(2,I) 355 CONTINUE Z=Z-G(2,7)-G(2,8)-G(2,11)-G(2,12)-G(2,13) DO 360 I=1,13 G(2,I)=G(2,I)/Z 360 CONTINUE IF (DFO2.LE.-5.) GO TO 370 FO2=(-25738./TEMP)+9.+DFO2 C PRINT *,'FO2',FO2,'TEMP',TEMP,'DFO2' RTIO=0.50235*FO2+13184.7/TEMP-4.54-2.15*G(2,4) RTIO=RTIO-8.351*G(2,3) RTIO=RTIO-4.494*(G(2,9)+G(2,10))+0.073*G(2,6)-5.4364*G(2,2) RTIO=RTIO+3.5415*G(2,1)+4.187*G(2,5) GO TO 390 370 RTIO=DLOG((G(1,10)/G(1,9))/2.D0) FO2=RTIO-13184.7/TEMP+4.54+2.15*G(2,4)+8.351*G(2,3) FO2=FO2+4.494*(G(2,9)+G(2,10))-0.073*G(2,6)+5.4364*G(2,2) FO2=FO2-3.5415*G(2,1)-4.187*G(2,5) FO2=FO2/0.50235 390 MIX=MIX+C WRITE(*,*)'FRACTION CRYSTALLIZED',F(9) WRITE(*,*)'FRACTION CRYSTALLIZED SINCE LAST OUTPUT',F(18) IF(G(1,2).LT.0.) GO TO 9999 IF(G(1,3).LT.0.) GO TO 9999 IF(G(1,7).LT.0.) GO TO 9999 IF(G(1,2)/ZTOT.LT.0.01) PER=MIX IF(G(1,3)/ZTOT.LT.0.09) PER=MIX IF(G(1,7)/ZTOT.LT.0.001) PER=MIX IF(G(1,9)/ZTOT.LT.0.01) PER=MIX ZTOT=G(1,2)/(G(1,2)+G(1,9)) IF(ZTOT.LT.0.1) PER=MIX WRITE(*,*)' PER ', PER, ' MIX ', MIX, ' MG# ', ZTOT IF(START.EQ.0.) GO TO 410 IF(MIX .GE. PER) GO TO 410 IF(SOLMAG.LT.0.9) GO TO 400 IF(F(18).GT. 0.02) GO TO 410 400 GO TO 100 410 PCTFR=F(9)*100. TEMP=TEMP-273. FETOT=G(1,9)+G(1,10) IF(DFO2.LE.-5.) GO TO 412 G(1,9)=FETOT/(2.*EXP(RTIO)+1.) G(1,10)=FETOT-G(1,9) 412 Z=0. DO 425 I=1,NT G(1,I)=G(1,I)*WT(I) Z=Z+G(1,I) 425 CONTINUE DO 430 I=1,NT G(1,I)=(G(1,I)*100.)/Z 430 CONTINUE G(1,8)=G(1,8)*10000. G(1,11)=G(1,11)*10000. FO=W(1,2)*100./(W(1,9)+W(1,2)+.000000001) CEN=W(3,2)*100./(W(3,2)+W(3,6)+W(3,9)+.000000001) CWO=W(3,6)*100./(W(3,2)+W(3,6)+W(3,9)+.000000001) CFS=W(3,9)*100./(W(3,2)+W(3,6)+W(3,9)+.000000001) OEN=W(2,2)*100./(W(2,2)+W(2,6)+W(2,9)+.000000001) OWO=W(2,6)*100./(W(2,2)+W(2,6)+W(2,9)+.000000001) OFS=W(2,9)*100./(W(2,2)+W(2,6)+W(2,9)+.000000001) AN=W(4,6)*100./(W(4,6)+W(4,1)+W(4,5)+.000000001) DO 440 I=1,8 Z=0. DO 435 J=1,NT W(I,J)=W(I,J)*WT(J) Z=Z+W(I,J) 435 CONTINUE IF(Z.LT.0.1) GO TO 440 DO 450 J=1,NT W(I,J)=(W(I,J)*100.)/Z 450 CONTINUE 440 CONTINUE DO 455 K=14,NT G(1,K)=G(1,K)*10000. 455 CONTINUE DO 460 I=1,8 W(I,11)=W(I,11)*10000. DO 465 K=14,NT W(I,K)=W(I,K)*10000. 465 CONTINUE IF (I.EQ.5) GO TO 460 W(I,8)=W(I,8)*10000. 460 CONTINUE GO TO 800 470 DO 475 J=1,8 CUM(J)=0. DO 480 I=1,NT W(J,I)=0. 480 CONTINUE 475 CONTINUE IF(XCAO(3).GT.0.15) THEN XCAO(3)=0.02 XMGO(3)=0.46 ENDIF IF(XCAO(1).LT.0.1) THEN XCAO(1)=0.2 XMGO(1)=0.2 ENDIF F(18)=0. IC=IC+1 IF(IC.GE.9000) IC=IC/10 C C GENERATION OF GAUSIAN RANDOM NUMBERS USING ABROMOWITZ AND C STEGUN. GAUSIAN NUMBERS ARE GENERATED FROM UNIFORM RANDOM # C RDM(1)=DND(IC+28) 9059 P=RDM(1) SIGN=1.0 IF (P.GT.0.5) GO TO 9033 GO TO 9044 9033 SIGN=-1.0 P=P-0.5 9044 IF (ABS(P).LT.1.E-15) P=.0001 RP=1./(P**2) TSQ=DLOG(RP) TSQ=SQRT(TSQ) A=2.515517 +0.802853*TSQ+0.010328*TSQ**2 B=1. + 1.432788*TSQ+ 0.189269*TSQ**2 + 0.001308*TSQ**3 Z=TSQ - (A/B) + .00004 RDM(1)=(SIGN*Z)/SQRT(2.0) RDM(1)=RDM(1)/2. RDM(1)=RDM(1)+1.0 IF(RDM(1).LE.0.0) RDM(1)=.001 START=1. GO TO 100 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C OLIVINE TABULATION SECTION C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 520 DO 521 I=1,NT X(10,I)=0. 521 CONTINUE X(10,2)=OM X(10,4)=1./3. X(10,6)=.002 X(10,7)=OT X(10,8)=RC X(10,9)=OF X(10,11)=OLNI X(10,13)=MNOL DO 525 I=14,NT X(10,I)=OLTR(I) 525 CONTINUE TEMP=T(1) F(1)=F(1)+C*(1.-ABS(FRAC(1))) CUM(1)=CUM(1)+C*FRAC(1) DO 530 I=1,NT X(1,I)=F(1)*X(10,I) W(1,I)=X(10,I) G(4,I)=G(4,I)-FRAC(1)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 530 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ORTHOPYROXENE TABULATION SECTION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 540 TEMP=T(2) C PRINT *,'ORTHOPYROXENE' F(2)=F(2)+C*(1.-FRAC(2)) CUM(2)=CUM(2)+C*FRAC(2) DO 550 I=1,NT X(2,I)=F(2)*X(10,I) W(2,I)=X(10,I) G(4,I)=G(4,I)-FRAC(2)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 550 CONTINUE DO 551 I=1,10 C PRINT *,ELEM(I),G(1,I) 551 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CLINOPYROXENE TABULATION SECTION C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 560 TEMP=T(3) F(3)=F(3)+C*(1.-ABS(FRAC(3))) CUM(3)=CUM(3)+C*FRAC(3) DO 570 I=1,NT X(3,I)=F(3)*X(10,I) W(3,I)=X(10,I) G(4,I)=G(4,I)-FRAC(3)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 570 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PLAGIOCLASE TABULATION SECTION C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 580 DO 581 I=1,NT X(10,I)=0. 581 CONTINUE X(10,1)=NA X(10,2)=MGPLAG X(10,4)=.4+KSP+NA X(10,6)=AN X(10,5)=KSP X(10,3)=.4-NA-KSP X(10,9)=FEPLAG DO 585 I=14,NT X(10,I)=PLTR(I) 585 CONTINUE TEMP=T(4) F(4)=F(4)+C*(1.-ABS(FRAC(4))) CUM(4)=CUM(4)+C*FRAC(4) DO 590 I=1,NT X(4,I)=F(4)*X(10,I) W(4,I)=X(10,I) G(4,I)=G(4,I)-FRAC(4)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 590 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 C C SPINEL TABULATION SECTION C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 600 DO 601 I=1,NT X(10,I)=0. 601 CONTINUE X(10,2)=SM X(10,10)=S3 DTI=ST/TITAN X(10,7)=ST X(10,8)=SC X(10,9)=SF X(10,11)=SN X(10,13)=SMN X(10,3)=SA IF(SC.LT.0.00001) GO TO 604 IF(SC*C.GT..021*G(1,8)) C=(.02*G(1,8))/SC 604 DO 605 I=14,NT X(10,I)=SPTR(I) 605 CONTINUE DO 3665 I=25,37 3665 X(10,I)=(YCEPT(5,I)*G(1,I))/0.4 IF (DTI.GT.4.0) GO TO 602 X(10,21)=0.026*DTI*ZRLIQ X(10,22)=0.026*DTI*BLIQ X(10,38)=0.024*DTI*FLIQ X(10,39)=0.021*DTI*TALIQ GO TO 3693 602 X(10,21)=(SLOPE(5,21)*DTI+YCEPT(5,21))*ZRLIQ X(10,22)=(SLOPE(5,22)*DTI+YCEPT(5,22))*BLIQ X(10,38)=(SLOPE(5,38)*DTI+YCEPT(5,38))*FLIQ X(10,39)=(SLOPE(5,39)*DTI+YCEPT(5,39))*TALIQ 3693 C=C/5 TEMP=T(5) F(5)=F(5)+C*(1.-ABS(FRAC(5))) CUM(5)=CUM(5)+C*FRAC(5) DO 610 I=1,NT X(5,I)=F(5)*X(10,I) W(5,I)=X(10,I) G(4,I)=G(4,I)-FRAC(5)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 610 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ILMENITE TABULATION SECTION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 620 DO 621 I=1,NT X(10,I)=0. 621 CONTINUE X(10,2)=IMG X(10,7)=.5 X(10,8)=ICR X(10,9)=IFE DO 625 I=14,NT X(10,I)=ILMTR(I) 625 CONTINUE DO 4605 I=25,37 X(10,I)=(YCEPT(6,I)*G(1,I))/0.4 4605 CONTINUE TEMP=T(6) F(6)=F(6)+C*(1.-ABS(FRAC(6))) CUM(6)=CUM(6)+C*FRAC(6) DO 630 I=1,NT X(6,I)=F(6)*X(10,I) W(6,I)=X(10,I) G(4,I)=G(4,I)-FRAC(6)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 630 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C PIGEONITE TABULATION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 632 TEMP=T(7) F(7)=F(7)+C*(1.-ABS(FRAC(7))) CUM(7)=CUM(7)+C*FRAC(7) DO 637 I=1,NT X(7,I)=F(7)*X(10,I) W(7,I)=X(10,I) G(4,I)=G(4,I)-FRAC(7)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 637 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C APATITE TABULATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 638 TEMP=T(8) DO 9640 I=1,NT X(10,I)=0. 9640 CONTINUE C=C*0.1 DO 9665 I=14,NT X(10,I)=APTR(I) 9665 CONTINUE X(10,6)=5./8.-SUMAP X(10,12)=3./8. F(8)=F(8)+C*(1.-ABS(FRAC(8))) CUM(8)=CUM(8)+C*FRAC(8) DO 9637 I=1,NT X(8,I)=F(8)*X(10,I) W(8,I)=X(10,I) G(4,I)=G(4,I)-FRAC(8)*C*X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 9637 CONTINUE F(18)=F(18)+C F(9)=F(9)+C GO TO 340 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C OLIVINE RESORPTION SECTION C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 640 X(10,2)=OM X(10,1)=0. X(10,3)=.0 X(10,5)=0. X(10,10)=0. X(10,4)=1./3. X(10,6)=.002 X(10,7)=OT X(10,8)=RC X(10,9)=OF X(10,11)=OLNI X(10,13)=MNOL DO 645 I=14,NT X(10,I)=OLTR(I) 645 CONTINUE F(1)=F(1)-C DO 660 I=1,NT X(1,I)=F(1)*X(10,I) W(1,I)=X(10,I) G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) 660 CONTINUE F(18)=F(18)+C F(9)=F(9)-C MIX=MIX-C GO TO 349 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C SECTION TO OUTPUT TO BLF.OUT C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 800 WRITE(3,805) PCTFR,SIZEMAG,TEMP,FO2 805 FORMAT(/,/,' PERCENT CRYSTALLIZED ',F7.1,5X,' SYSTEM SIZE ', :F7.3,/,' TEMPERATURE ',F6.1,' LOG FO2 ',F7.2,/) WRITE(3,810) 810 FORMAT(' OLIV OPX AUG PLAG ', :' SPINEL ILMEN PIG APAT') WRITE(3,815)(F(I),I=1,8) 815 FORMAT(' EQUIL XLS',8F8.4) WRITE(3,820)(CUM(I),I=1,8) 820 FORMAT(' FRAC XLS ',8F8.4) WRITE(3,870) 840 FORMAT(F12.6) 870 FORMAT(/,' MELT COMPOSITION (OXIDE WT %) ',/) WRITE(3,875) 875 FORMAT(' NA MG AL SI K CA TI CR FE', :' FE+3 NI P MN') START=1.0 IC=IC+1 IF(IC.GE.9000) IC=IC/5 WRITE(3,880) (G(1,I),I=1,13) 880 FORMAT(7F6.2,F6.0,2F6.2,F6.0,2F6.2,/) DO 890 I=14,NT WRITE(3,885) ELEM(I),G(1,I) 885 FORMAT(A,F12.3) 890 CONTINUE WRITE(2,*)PCTFR,CHAR(09),TEMP,CHAR(09),(G(1,I),CHAR(09),I=1,NT) 894 FORMAT(F12.5) WRITE(3,895) 895 FORMAT(' EQUILIBRIUM MINERALOGY ',/) 900 IF(W(1,4).LT.0.0001) GO TO 920 WRITE(3,905) FO 905 FORMAT(' OLIVINE ',F5.1,' FO ') WRITE(3,910)(W(1,I),I=1,13) 910 FORMAT(7F6.2,F6.0,2F6.2,F6.0,2F6.2) DO 916 K=14,NT WRITE(3,915) ELEM(K),W(1,K) 915 FORMAT(A,F12.2) 916 CONTINUE 920 IF(W(2,4).LT.0.0001) GO TO 940 WRITE(3,925) OEN,OWO,OFS 925 FORMAT(' OPX ',F5.1,' EN ',F5.1,' WO ',F5.1,' FS ') WRITE(3,930)(W(2,I),I=1,13) 930 FORMAT(7F6.2,F6.0,2F6.2,F6.0,2F6.2) DO 936 K=14,NT WRITE(3,935) ELEM(K),W(2,K) 935 FORMAT(A,F12.2) 936 CONTINUE 940 IF(W(3,4).LT.0.0001) GO TO 960 WRITE(3,945) CEN,CWO,CFS 945 FORMAT(' AUGITE ',F5.1,' EN ',F5.1,' WO ',F5.1,' FS ') WRITE(3,950)(W(3,I),I=1,13) 950 FORMAT(7F6.2,F6.0,2F6.2,F6.0,2F6.2) DO 956 K=14,NT WRITE(3,955) ELEM(K),W(3,K) 955 FORMAT(A,F12.2) 956 CONTINUE 960 IF(W(4,4).LT.0.0001) GO TO 980 WRITE(3,965) AN 965 FORMAT(' PLAGIOCLASE ',F6.1,' AN ') WRITE(3,970)(W(4,I),I=1,13) 970 FORMAT(13F6.2) DO 976 K=14,NT WRITE(3,975) ELEM(K),W(4,K) 975 FORMAT(A,F12.2) 976 CONTINUE 980 IF(W(5,9).LT.0.0001) GO TO 1000 WRITE(3,985) 985 FORMAT(' SPINEL ') WRITE(3,990)(W(5,I),I=1,13) 990 FORMAT(10F6.2,F6.0,2F6.2) DO 996 K=14,NT WRITE(3,995) ELEM(K),W(5,K) 995 FORMAT(A,F12.2) 996 CONTINUE 1000 IF(W(6,7).LT.0.0001) GO TO 1940 WRITE(3,1005) 1005 FORMAT(' ILMENITE ') WRITE(3,1010)(W(6,I),I=1,13) 1010 FORMAT(13F6.2) DO 1016 K=14,NT WRITE(3,1015) ELEM(K),W(6,K) 1015 FORMAT(A,F7.2) 1016 CONTINUE 1940 IF(W(7,4).LT.0.0001) GO TO 2002 WRITE(3,1945) XMGO(3)*200.,XCAO(3)*200.,XFEO(3)*200. 1945 FORMAT(' PIG ',F5.1,' EN ',F5.1,' WO ',F5.1,' FS ') WRITE(3,1950)(W(7,I),I=1,13) 1950 FORMAT(7F6.2,F6.0,2F6.2,F6.0,2F6.2) DO 1956 K=14,NT WRITE(3,1955) ELEM(K),W(7,K) 1955 FORMAT(A,F7.2) 1956 CONTINUE 2002 IF(W(8,12).LT.0.0001) GO TO 470 WRITE(3,2005) 2005 FORMAT(' APATITE ') WRITE(3,2010)(W(8,I),I=1,13) 2010 FORMAT(13F6.2) DO 2016 K=14,NT WRITE(3,2015) ELEM(K),W(8,K) 2015 FORMAT(A,F12.2) 2016 CONTINUE GO TO 470 C C PERIODIC MIXING SECTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 9000 SUM=0. SUM1=0. SUM6=0. IF (IC.GE.9000) IC=IC/100 DO 9002 I=2,5 RDM(I)=DND(IC+I*16) 9003 P=RDM(I) SIGN=1.0 IF(P.GT.0.5) GO TO 9004 GO TO 9006 9004 SIGN=-1.0 P=P-0.5 9006 IF(ABS(P).LT.1.E-15) P=.0001 RP=1/(P**2) TSQ=DLOG(SQRT(RP)) A=2.515517+0.802853*TSQ+0.010328*TSQ**2 B=1.+1.432788*TSQ+.189269*TSQ**2+.001308*TSQ**3 Z=TSQ-(A/B)+0.00004 RDM(I)=(SIGN*Z)/SQRT(2.) RDM(I)=RDM(I)/2. RDM(I)=RDM(I)+1.0 IF(RDM(I).LE.0.0) RDM(I)=0.001 WRITE(*,*) I,RDM(I) 9002 CONTINUE PER=PERI IF(RD(2).EQ.1.0) PER=PERI*RDM(2) IF(RD(3).EQ.1.0) RECHG=RECHGI*RDM(3) IF(RD(4).EQ.1.0) ASSIM=ASSIMI*RDM(4) IF(RD(5).EQ.1.0) ERUPT=ERUPTI*RDM(5) IF(PER .GT. 0.85) PER=PERI DO 9005 I=1,NT SUM=SUM+G(4,I) SUM6=SUM6+G(6,I) 9005 CONTINUE IF (SUM.GE.ERUPT*MIX) GO TO 9199 WRITE(*,*) 'SYSTEM PURGE',' MIX ',MIX,' SUM ',SUM PRINT *,'SYSTEM PURGE',' MIX ',MIX,' SUM ',SUM DO 9109 I=1,NT G(4,I)=0.0 9109 CONTINUE SUM=1. 9199 IC=IC+1 SIZEMAG=0. DO 9010 I=1,NT G(1,I)=G(4,I)-X(1,I)-X(2,I)-X(3,I)-X(4,I)-X(5,I) :-X(6,I)-X(7,I)-X(8,I) SUM1=SUM1+G(1,I) G(6,I)=G(6,I)-SOLMAG*G(6,I)/SUM6 :+SOLMAG*MIX*ASSIM*G(3,I)+SOLMAG*G(1,I) SIZEMAG=SIZEMAG+G(6,I) 9010 CONTINUE IF(SIZEMAG.LT..05) GO TO 9999 C C C COMPILE THE FRACTION XLIZED SINCE THE LAST RECHARGE C PCXLSL=PCXLSL+SOLMAG*MIX IF(UT.EQ.0.) GO TO 9035 IF (PCXLSL .LT. UT) GO TO 9035 DO 9036 I=1,NT G(6,I)=G(6,I)-(ERUPT*PCXLSL*(G(6,I)/SIZEMAG))+PCXLSL*RECHG*G(5,I) SIZEMAG=SIZEMAG-(ERUPT*PCXLSL*(G(6,I)/SIZEMAG)) :+PCXLSL*RECHG*G(5,I) 9036 CONTINUE C C RANDOMIZE PERIODICITY OF RECHARGE C IF(RD(1).EQ.1.0) UT=UTI*RDM(1) WRITE(*,519) UT 519 FORMAT(' NEXT RECHARGE IN ',F8.5,' FRACTION XLIZED ') WRITE(3,9016) 9016 FORMAT(' RECHARGE EPISODE ',/) PCXLSL=0 9035 WRITE(3,9015) SIZEMAG 9015 FORMAT(' MIXING EPISODE SYSTEM SIZE = ',F12.8,/) DO 9025 I=1,NT G(4,I)=G(6,I)/SIZEMAG C WRITE(3,9014) ELEM(I),G(4,I),G(6,I) C 9014 FORMAT(A,2F12.8) 9025 CONTINUE MIX=0. START=0. DO 9020 J=1,8 DO 9019 I=1,NT X(J,I)=0.0 9019 CONTINUE F(J)=0. 9020 CONTINUE GO TO 110 9999 CLOSE(UNIT=3) CLOSE(UNIT=2) CLOSE(UNIT=1) STOP END C C SUBROUTINE PXCALC(INDEXA,XL,UL,DULDT,DULDX,TP,H,DZ,R) C C *********************** C CALCULATION OF COTECTIC COMPOSITION AND TEMPERATURE OF C SOLIDS COEXISTING WITH BASALTIC LIQUIDS GOVERNED BY C NIELSEN'S MODEL. --PM DAVIDSON, NBS, MAY 1986 C ************************ C C PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP) DIMENSION H(100),DZ(100),TP(4) DATA MAXFN/80/ H(29)=1.D0 XSET=H(22) YSET=H(23) ICONT=0 H(9)=XL(1)+XL(2)+XL(3) CALL EQUIL(XSET,YSET,INDEXA,MAXFN,ICONT,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) C PRINT *,'AFTER EQUIL' RETURN 3021 STOP END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE EQUIL(XSET,YSET,INDEXA,MAXFN,ICONT) C CALCULATION OF T AND X5(N) COMPOSITIONAL VARIABLES FOR PHASE EQM C GENERALIZED NEWTON-RAPHSON ITERATION TO FIND X,T C P.M. DAVIDSON, NBS, MARCH 86. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE EQUIL(XSET,YSET,INDEXA,MAXFN,ICONT,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) PARAMETER(NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION UJAC,V,SIMUL,DZ PARAMETER(NMAX=5,MMAX=6) C C VARIABLES--X5(1),X5(2),...X5(NMAX-1) MOL FRACTIONS OF COMPONENTS IN C PHASE 2 C X5(NMAX)=T(K) C DZ(1), DEL-MU CONVERGENCE LIMIT C DZ(2), MATRIX SINGULARITY LIMIT C DZ(3), INCREMENT CONVERGENCE C DZ(4), MAX. FRACTIONAL INCREMENT C NMAX, NUMBER OF COMPONENTS(CONSTRAINTS) C MMAX =NMAX+1 DIMENSION X5(MMAX),UJAC(MMAX,MMAX),V(MMAX) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP) DIMENSION TP(4),H(100),DZ(100) C DATA MAXFN/20/ DATA SCALE/1.D0/ C C INITIALIZATION ERCT=1. DO 1 I=1,MMAX X5(I)=0.D0 V(I)=0.D0 DO 1 J=1,MMAX 1 UJAC(I,J)=0.D0 C X5(1)=XSET X5(2)=YSET C C IF((X5(1)+X5(2)).GE.1.) THEN X5(1)=X5(1)/(X5(1)+X5(2)+.01) X5(2)=X5(2)/(X5(1)+X5(2)+.01) END IF C DETERMINE NDIM, NUMBER OF UNKNOWNS AND MDIM 5 NDIM=3 C 8 MDIM=NDIM+1 X5(NDIM)=H(30) IER=0 DO 100 ITER=1,MAXFN IF(IER.GT.5) GO TO 80 IF(H(30).LT.500.) GO TO 80 C CALCULATE JACOBIAN MATRIX UJAC, AND DELMU(I) CALL DERIV(X5,UJAC,INDEXA,NDIM,MDIM,ICONT,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) C PRINT *, 'AFTER DERIV' STAB=UJAC(2,1)*UJAC(3,2)-UJAC(2,2)*UJAC(2,2) C SCALE TEMPERATURE,DIAGONAL SCITER=1.D0+.2D0/(1.D0*ITER) C DO 9 I=1,NDIM UJAC(I,NDIM)=UJAC(I,NDIM)*SCALE UJAC(I,I)=UJAC(I,I)*SCITER 9 CONTINUE C WRITE(*,*) 'ITER ',ITER C WRITE(*,210) (X5(L),L=1,NDIM) C WRITE(*,210) ((UJAC(L,M),M=1,MDIM),L=1,NDIM) C WRITE(*,*) ' STAB= ',STAB C 210 FORMAT(4E12.5) C CHECK DEL-MU CONVERGENCE ITCON=1 DO 10 I=1,NDIM 10 IF(ABS(UJAC(I,MDIM))-DZ(1).GT.0.) ITCON=0 IF (ITCON) 15,15,300 C C FIND INCREMENT TO X5(I) = V(I) 15 INDIC=0 C PRINT *,'BEFORE SIMUL' DETER=SIMUL(NDIM,UJAC,V,INDIC,MDIM) C PRINT *,'AFTER SIMUL' C CHECK FOR SINGULARITY C WRITE(*,*) ' FROM DETER:UJAC,V,NDIM,MDIM ',NDIM,MDIM C WRITE(*,210) ((UJAC(MM,NN),NN=1,MDIM),MM=1,NDIM) C WRITE(*,210) (V(MM),MM=1,NDIM) C WRITE(*,*) ' DETER: ',DETER,' STAB= ',STAB IF (DETER) 20,80,20 C C CHECK INCREMENT SIZE 20 ITCON=2 IF(STAB.LE.0.) THEN ITCON=0 IER=IER+1 IF(INDEXA.EQ.1) THEN X5(1)=X5(1)-.05D0 GO TO 100 END IF IF(INDEXA.EQ.3) THEN X5(2)=0.75D0*X5(2) END IF END IF C DO 30 I=1,NDIM VI=V(I)/X5(I) IF(ABS(VI).GE.DZ(3)) THEN ITCON=0 IF(ABS(VI).GE.DZ(4)) VI=SIGN(DZ(4),VI) V(I)=VI*X5(I) END IF 30 CONTINUE IF(ITCON) 40,40,300 C C TEST NEW VALUES 40 IER=0 IMAX=NDIM-1 XT=0.D0 DO 50 I=1,IMAX XI=X5(I)+V(I) IF(I.EQ.2) THEN IF(XI.GE..5D0) IER=1 ENDIF C ENSURE X5 < 1 IF(XI*XI.GT.XI) IER=1 IF(I.LT.3.OR.ICONT.EQ.0) THEN XT=XT+XI IF(XT*XT.GT.XT) IER=1 END IF X5(I)=XI 50 CONTINUE C IF (ICONT.EQ.0) THEN X5(NDIM)=X5(NDIM)+V(NDIM)*SCALE IF(X5(NDIM).LE.0.) IER=1 ENDIF C C ADJUST X TO BE BETWEEN 0 AND 1 IF(IER.EQ.1) THEN IF(ERCT.GT.100.) GO TO 80 AERCT=1.D0+5.D-03*ERCT ERCT=ERCT*2. DO 60 I=1,NDIM 60 V(I)=-V(I)/AERCT GO TO 40 END IF C 100 CONTINUE C C NO SOLUTION 80 WRITE(*,*)' NO SOLUTION.' WRITE(*,201) (X5(K),K=1,MMAX) 201 FORMAT(7E12.5) C XCON=REAL(ITCON) 300 CONTINUE H(22)=X5(1) H(23)=X5(2) H(30)=X5(3) RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE DERIV(X5,UJAC,INDEXA,NDIM,MDIM,ICONT,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) C C CALCULATION OF ELEMENTS OF JACOBIAN MATRIX UJAC C UJAC(I,J)=D(DELTA MU(I))/DX(J) C PARAMETER(NMAX=5,MMAX=6,NCOMP=10) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION UJAC DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP) DIMENSION TP(4),H(100),DZ(100) DIMENSION X5(MMAX),UJAC(MMAX,MMAX) COMMON DMU1,DMU2,DMU3,GX,GY,DTDX,DTDY,DGDT,GTT,GTX,GTY,GXX,GXY,GYY C IF(ICONT.EQ.0) THEN H(30)=X5(3) ENDIF H(24)=R*H(30) C CALL ULIQ(INDEXA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) C GO TO(10,20,30), INDEXA C C INDEXA=1: LIQ-CPX. 2: LIQ-OPX. 3:LIQ-PIG C 10 IF(X5(1).EQ.0.D0) THEN C INITIALIZE AUG X5(1)=.54 X5(2)=.45 END IF C 30 IF(X5(1).EQ.0.D0) THEN C INITIALIZE PIG X5(1)=.90 X5(2)=.09 END IF C XA=X5(1) YA=MIN(X5(2),.99999D0-XA) IOPX=0 CALL MARGUL(XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R ) CALL MU0 (XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) U0P1=H(26) U0P2=H(27) U0P4=H(28) DTU0P1=DZ(33) DTU0P2=DZ(34) DTU0P4=DZ(36) H(44)=H(43)+H(30)*DZ(21)+H(29)*DZ(22) CALL SITE(XA,YA,TA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) CALL MU1(XA,YA,TA,XL,UL,DULDT,DULDX,TP,H,DZ,A,B,C,D,GT1,R) CALL PTDER(XA,YA,TA,DGDTK,DGDP,GTTK,GPT,GTKX,GTKY, 1GPY,GPX,IOPX,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) GO TO 50 C 20 IF(X5(1).EQ.0.D0) THEN X5(1)=.95 X5(2)=.04 END IF XA=X5(1) YA=MIN(X5(2),.99999D0-XA) IOPX=1 H(31)=H(39)+H(30)*DZ(20) CALL MU0(XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) U0P1=H(40) U0P2=H(41) U0P4=H(42) DTU0P1=DZ(38) DTU0P2=DZ(39) DTU0P4=DZ(41) CALL SITOPX(XA,YA,TA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) CALL MU1OPX(XA,YA,TA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) CALL PTDER(XA,YA,TA,DGDTK,DGDP,GTTK,GPT,GTKX,GTKY, 1GPY,GPX,IOPX,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) GO TO 50 C C 50 UJAC(1,MDIM)=UL(1)-DMU1 1 +U0P1 UJAC(2,MDIM)=UL(1)-UL(2)-GX 1 +U0P1-U0P2 UJAC(3,MDIM)=2.D0*(UL(4)-UL(2))-GY 1 +2.D0*(U0P4-U0P2) UJAC(2,1)=GXX+DTDX*GTX UJAC(2,2)=GXY+DTDX*GTY UJAC(3,1)=UJAC(2,2) UJAC(3,2)=GYY+DTDY*GTY UJAC(1,1)=(1.D0-XA)*UJAC(2,1)-YA*UJAC(3,1) UJAC(1,2)=(1.D0-XA)*UJAC(2,2)-YA*UJAC(3,2) DU1DTK=DGDTK+(1.D0-XA)*(GTKX+GTTK*DTDX)-YA*(GTKY+GTTK*DTDY) UJAC(1,3)=DU1DTK-DULDT(1)-DTU0P1 UJAC(2,3)=GTKX+DTDX*GTTK - DULDT(1)+DULDT(2) 1 -DTU0P1+DTU0P2 UJAC(3,3)=GTKY+DTDY*GTTK-2.D0*(DULDT(4)-DULDT(2)) 1 -2.D0*(DTU0P4-DTU0P2) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C FUNCTIONS AND SUBROUTINES TO DETERMINE FREE ENERGY AND DERIVATIVES C C AND EQUILIBRIUM DEGREE OF ORDER FOR OPX AND CPX C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE MARGUL(XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) C C C EVALUATES FUNCTIONS OF MARGULES PARAMS C UPDATA MAY 86 TO LOAD ALL PYX SOLUTION PARAMS AT P,T C PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) W13=DZ(43)+H(29)*DZ(44) W31=DZ(45)+H(29)*DZ(46) W23=DZ(47)+H(29)*DZ(48) W32=DZ(49)+H(29)*DZ(50) A=W13-W31+W23-W32 B=W13-W31-W23+W32 C=W13+W31-3.D0*W23+W32 D=W13+W31+5.D0*W23-3.D0*W32 RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE MU1(XA,YA,TA,XL,UL,DULDT,DULDX,TP,H,DZ,A,B,C,D,GT1,R) PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) COMMON DMU1,DMU2,DMU3,GX,GY,DTDX,DTDY,DGDT,GTT,GTX,GTY,GXX,GXY,GYY DM1=XA+YA+TA/2.D0 DM2=2.D0*XA-DM1 FE1=1.D0-DM1 C2=2.D0*YA FE2=1.D0-DM2-C2 BB=H(44)*(XA*(1.D0-XA-YA)+YA*TA/2.D0+TA*TA/4.D0)+H(25)*YA*(XA+YA 1+TA/2.D0) 1-H(45)*(YA*(1.D0-XA-YA)+(1.D0-YA)*TA/2.D0) TS=-H(24)*(DM1*DLOG(DM1)+DM2*DLOG(DM2)+FE1*DLOG(FE1)+ 1 FE2*DLOG(FE2)+C2*DLOG(C2)) H(99)=YA*(C*(XA-TA/2.D0)-D*YA+B*(2.D0*YA*TA-4.D0*XA*YA) 1 +A*(XA*XA+3.D0*YA*YA+TA*TA/4.D0-XA*TA)) 1 +.5D0*YA*(D-C-2.D0*(A-B)) 1 +2.D0*YA*(H(28)-H(27))+XA*(H(26)-H(27))+H(27) 1 +H(46)*DM2*FE2/2.D0 1 +H(47)*DM1*FE1/2.D0 GT1=BB-TS+H(99) DBDX=H(44)*(1.D0-2.D0*XA-YA)+H(25)*YA+H(45)*YA DBDY=H(44)*(-XA+TA/2.D0)+H(25)*(XA+2.D0*YA+TA/2.D0)-H(45)*(1.D0 1-XA-2.D0*YA-TA/2.D0) DTSDX=-H(24)*DLOG(DM1*DM2/(FE1*FE2)) DTSDY=-H(24)*DLOG(C2*C2*DM1/(FE1*FE2*DM2)) DHDX=YA*(C+A*(2.D0*XA-TA)-4.D0*B*YA) 1+H(46)*(1.D0-2.D0*XA+TA)/2.D0 1+H(47)*(.5D0-XA-YA-.5D0*TA) 2+H(26)-H(27) DHDY=C*(XA-TA/2.D0)-2.D0*YA*D+B*(4.D0*YA*TA-8.D0*XA*YA) 1 +A*(XA*XA+9.D0*YA*YA+TA*TA/4.D0-XA*TA) 1 +.5D0*(D-C-2.D0*(A-B)) 1 +2.D0*(H(28)-H(27)) 1 +H(46)*(2.D0*YA-1.D0)/2.D0 1+H(47)*(.5D0-XA-YA-.5D0*TA) GX=DBDX+DHDX-DTSDX GY=DBDY+DHDY-DTSDY DMU1=GT1+(1.D0-XA)*GX-YA*GY DMU2=GT1-XA*GX-YA*GY DMU3=GT1-XA*GX+(1.D0-YA)*GY C C C BXX=-2.D0*H(44) BYY=2.D0*(H(25)+H(45)) BXY=H(25)-H(44)+H(45) TSXX=-H(24)*(1.D0/DM1+1.D0/DM2+1.D0/FE1+1.D0/FE2) TSYY=-H(24)*(4.D0/C2+1.D0/DM1+1.D0/FE1+1.D0/FE2+1.D0/DM2) TSXY=-H(24)*(1.D0/DM1-1.D0/DM2+1.D0/FE1+1.D0/FE2) HXX=2.D0*YA*(A)-H(46)-H(47) HYY=-2.D0*D+18.D0*A*YA+B*(4.D0*TA-8.D0*XA) 1+H(46)-H(47) HXY=C+A*(2.D0*XA-TA)-8.D0*B*YA-H(47) GXX=BXX+HXX-TSXX GXY=BXY+HXY-TSXY GYY=BYY+HYY-TSYY DGDT=YA*(A*(TA-2.D0*XA)+4.D0*B*YA-C) 2+H(46)*(XA-.5D0-TA/2.D0)+H(44)*(YA+TA)+H(25)*YA-H(45)*(1.D0-YA) 1+H(47)*(.5D0-XA-YA-.5D0*TA) 3+H(24)*DLOG(DM1*FE2/(FE1*DM2)) GTT=H(44)-H(46)/2.D0+YA*(A)+H(24)*(1.D0/DM1+1.D0/DM2 1+1.D0/FE1+1.D0/FE2)/2.D0 1-H(47)/2.D0 DGDT=DGDT/2.D0 GTT=GTT/2.D0 GTX=YA*(-A)+H(46)/2.D0+H(24)*(1./DM1-1./FE2+1./FE1 1-1./DM2)/2. 1-H(47)/2.D0 GTY=A*(TA-2.D0*XA)+8.D0*B*YA-C 2+H(44)+H(25)+H(45)+H(24)*(1.D0/DM1-1.D0/FE2+1.D0/FE1+1.D0/DM2) 1-H(47) GTY=GTY/2.D0 DTDX=-GTX/GTT DTDY=-GTY/GTT RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SITE(XA,YA,TA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) DATA TOL1,TOL2/1.D-11,1.D-10/ DATA MAXFN/100/ IF(XA.GE..5) B0=2.*(1.D0-XA-YA) IF(XA.LT..5) B0=2.*(XA-YA) C THIS IS RANDOM DISTN A0=-2.*YA*(1.D0-XA-YA)/(1.D0-YA) IF((XA+YA).LT..5) A0=-2.*(XA+YA) IF((XA+YA).GE..5) A0=-2.*(1.D0-XA-YA) IF(A0.LT.B0) GO TO 2 TA=.5D0*(A0-B0) IER=4 GO TO 905 2 TNEW=A0+0.75*(B0-A0) C IF(A0.LT.TA.AND.TA.LT.B0) TNEW=TA 3 EPS=1.D-10 IER=0 C DO 200 I=1,MAXFN IF(EPS.GT.1.D-12) GO TO 1 IER=1 GO TO 905 1 TA=TNEW DM1=XA+YA+.5D0*TA DM2=2.D0*XA-DM1 FE1=1.D0-DM1 FE2=1.D0-DM2-2.D0*YA F=YA*(H(25)+H(44)+H(45)-C+A*(TA-2.D0*XA)+4.D0*B*YA)-H(45) 1 +TA*(H(44)-.5D0*H(46))+H(46)*(XA-.5D0) 1 +H(24)*DLOG(DM1*FE2/(FE1*DM2)) 1 +H(47)*(.5D0-XA-YA-.5D0*TA) FF=YA*A+H(44)-.5D0*H(46)+.5D0*H(24)*(1.D0/DM1+1.D0/DM2+1.D0/FE1+ 11.D0/FE2) -.5D0*H(47) TINC=F/FF TNEW=TA-TINC C SAFEGUARDS 30 IF(TNEW.GT.A0) GO TO 10 TNEW=A0+EPS EPS=EPS*.5D0 10 IF(TNEW.LT.B0) GO TO 20 TNEW=B0-EPS EPS=EPS*.5D0 GO TO 30 20 IF(DABS(TINC).LT.TOL1) GO TO 900 IF(DABS(F).LT.TOL2) GO TO 900 200 CONTINUE 900 TA=TNEW 905 IF(I.GE.MAXFN) IER=2 IF(IER.EQ.0) GO TO 999 999 RETURN END C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE MU1OPX(XA,YA,TA,XL,UL,DULDT,DULDX, 1TP,H,DZ,A,B,C,D,GT1,R) PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) COMMON DMU1,DMU2,DMU3,GX,GY,DTDX,DTDY,DGDT,GTT,GTX,GTY,GXX,GXY,GYY H(44)=H(31) H(45)=H(32) H(47)=H(33) H(46)=H(34) CO=H(35) DOO=H(36) F01O=H(37) DM1=XA+YA+TA/2.D0 DM2=2.D0*XA-DM1 FE1=1.D0-DM1 C2=2.D0*YA FE2=1.D0-DM2-C2 BB=H(44)*(XA*(1.D0-XA-YA)+YA*TA/2.D0+TA*TA/4.D0)+F01O*YA*(XA+YA 1 +TA/2.D0) 1-H(45)*(YA*(1.D0-XA-YA)+(1.D0-YA)*TA/2.D0) C STD STATES NOT OMITTED: UO10,20,30,40 NO LONGER ZERO]] TS=-H(24)*(DM1*DLOG(DM1)+DM2*DLOG(DM2)+FE1*DLOG(FE1)+ 1 FE2*DLOG(FE2)+C2*DLOG(C2)) H(99)=YA*(CO*(XA-TA/2.D0)-DOO*YA) 1 +.5D0*YA*(DOO-CO) 1 +H(46)*DM2*FE2/2.D0+H(47)*DM1*FE1/2.D0 1 +2.D0*YA*(H(42)-H(41))+XA*(H(40)-H(41))+H(41) GT1=BB-TS+H(99) DBDX=H(44)*(1.D0-2.D0*XA-YA)+F01O*YA+H(45)*YA DBDY=H(44)*(-XA+TA/2.D0)+F01O*(XA+2.D0*YA+TA/2.D0)-H(45)*(1.D0 1-XA-2.D0*YA-TA/2.D0) DTSDX=-H(24)*DLOG(DM1*DM2/(FE1*FE2)) DTSDY=-H(24)*DLOG(C2*C2*DM1/(FE1*FE2*DM2)) DHDX=YA*(CO-H(47))-H(47)*TA 1+(H(47)+H(46))*(1.D0-2.D0*XA+TA)/2.D0 2+H(40)-H(41) DHDY=CO*(XA-TA/2.D0)-2.D0*YA*DOO 1 -H(47)*(XA+TA/2.D0) 1 +.5D0*(DOO-CO) 1 +(H(46)-H(47))*(2.D0*YA-1.D0)/2.D0 1 +2.D0*(H(42)-H(41)) GX=DBDX+DHDX-DTSDX GY=DBDY+DHDY-DTSDY DMU1=GT1+(1.D0-XA)*GX-YA*GY DMU2=GT1-XA*GX-YA*GY DMU3=GT1-XA*GX+(1.D0-YA)*GY C C BXX=-2.D0*H(44) BYY=2.D0*(F01O+H(45)) BXY=F01O-H(44)+H(45) TSXX=-H(24)*(1.D0/DM1+1.D0/DM2+1.D0/FE1+1.D0/FE2) TSYY=-H(24)*(4.D0/C2+1.D0/DM1+1.D0/FE1+1.D0/FE2+1.D0/DM2) TSXY=-H(24)*(1.D0/DM1-1.D0/DM2+1.D0/FE1+1.D0/FE2) HXX=-H(46)-H(47) HYY=-2.D0*DOO+H(46)-H(47) HXY=CO-H(47) GXX=BXX+HXX-TSXX GXY=BXY+HXY-TSXY GYY=BYY+HYY-TSYY DGDT=YA*(-H(47)-CO)+H(47)*(+.5D0-XA-.5D0*TA) 2+H(46)*(XA-.5D0-TA/2.D0)+H(44)*(YA+TA)+F01O*YA-H(45)*(1.D0-YA) 3+H(24)*DLOG(DM1*FE2/(FE1*DM2)) GTT=H(44)-H(46)/2.D0-H(47)*.5D0+H(24)*(1.D0/DM1+1.D0/DM2 1+1.D0/FE1+1.D0/FE2)/2.D0 DGDT=DGDT/2.D0 GTT=GTT/2.D0 GTX=-H(47)/2.D0+H(46)/2.D0+H(24)*(1.D0/DM1-1.D0/FE2+1.D0/FE1 1-1.D0/DM2)/2.D0 GTY=-CO-H(47) 2+H(44)+F01O+H(45)+H(24)*(1.D0/DM1-1.D0/FE2+1.D0/FE1+1.D0/DM2) GTY=GTY/2. DTDX=-GTX/GTT 10 DTDY=-GTY/GTT RETURN END C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SITOPX(XA,YA,TA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) DATA TOL1,TOL2/1.D-11,1.D-10/ DATA MAXFN/100/ H(44)=H(31) H(45)=H(32) H(47)=H(33) H(46)=H(34) CO=H(35) DOO=H(36) F01O=H(37) IF(XA.GE..5) B0=2.*(1.D0-XA-YA) IF(XA.LT..5) B0=2.*(XA-YA) C THIS IS RANDOM DISTN A0=-2.*YA*(1.D0-XA-YA)/(1.D0-YA) IF((XA+YA).LT..5) A0=-2.*(XA+YA) IF((XA+YA).GE..5) A0=-2.*(1.D0-XA-YA) IF(A0.LT.B0) GO TO 2 TA=.5D0*(A0-B0) IER=4 GO TO 905 C INITIALIZE TA,EPS,IER 2 TNEW=A0+.75*(B0-A0) C IF(A0.LT.TA.AND.TA.LT.B0) TNEW=TA 3 EPS=1.D-10 IER=0 C DO 200 I=1,MAXFN IF(EPS.GT.1.D-12) GO TO 1 IER=1 GO TO 905 1 TA=TNEW DM1=XA+YA+.5D0*TA DM2=2.D0*XA-DM1 FE1=1.D0-DM1 FE2=1.D0-DM2-2.D0*YA F=YA*(F01O+H(44)+H(45)-CO-H(47))-H(45) 1 +TA*(H(44)-.5D0*H(46)-.5D0*H(47))+(H(46)-H(47))*(XA-.5D0) 1 +H(24)*DLOG(DM1*FE2/(FE1*DM2)) FF=-.5D0*(H(46)+H(47))+.5D0*H(24)*(1.D0/DM1+1.D0/DM2+1.D0/FE1 1+1.D0/FE2)+H(44) TINC=F/FF TNEW=TA-TINC C SAFEGUARDS 30 IF(TNEW.GT.A0) GO TO 10 TNEW=A0+EPS EPS=EPS*.5D0 10 IF(TNEW.LT.B0) GO TO 20 TNEW=B0-EPS EPS=EPS*.5D0 GO TO 30 20 IF(DABS(TINC).LT.TOL1) GO TO 900 IF(DABS(F).LT.TOL2) GO TO 900 200 CONTINUE 900 TA=TNEW 905 IF(I.GE.MAXFN) IER=2 IF(IER.EQ.0) GO TO 999 999 RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION SIMUL(NDIM,AAA,V,INDIC,MDIM) C TAKEN FROM APPLIED NUM. METHODS, P.290-291 PARAMETER(NCOMP=10,NMAX=5,MMAX=6) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IROW(50),JCOL(50),JORD(50),Y(50),AAA(MMAX,MMAX),V(MMAX) DIMENSION DZ(100) MAX=NDIM INDIC=0 IF(INDIC.GE.0) MAX=NDIM+1 IF(NDIM.LE.50) GO TO 5 WRITE(*,200) SIMUL=0. RETURN C ...BEGIN ELIMINATION PROCEDURE... 5 DETER=1. DO 18 K=1,NDIM KM1=K-1 C SEARCH FOR PIVOT ELEMENT PIVOT=0.0 DO 11 I=1,NDIM DO 11 J=1,NDIM C SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS IF(K.EQ.1) GO TO 9 DO 8 ISCAN=1,KM1 DO 8 JSCAN=1,KM1 IF(I.EQ.IROW(ISCAN)) GO TO 11 IF (J.EQ.JCOL(JSCAN)) GO TO 11 8 CONTINUE 9 IF(DABS(AAA(I,J)).LE.DABS(PIVOT)) GO TO 11 PIVOT=AAA(I,J) IROW(K)=I JCOL(K)=J 11 CONTINUE C INSURE THAT SELECTED PIVOT IS LARGER THAN DZ(2)(OLD EPS2) IF (DABS(PIVOT).GT.DZ(2)) GO TO 13 SIMUL=0. RETURN C UPDATE THE DETERMINANT VALUE 13 IROWK=IROW(K) JCOLK=JCOL(K) DETER=DETER*PIVOT C NORMALIZE PIVOT ROW ELEMENTS DO 14 J=1,MAX 14 AAA(IROWK,J)=AAA(IROWK,J)/PIVOT C CARRY OUT ELIMINATION AND DEVELOP INVERSE AAA(IROWK,JCOLK)=1.D0/PIVOT DO 18 I=1,NDIM AIJCK=AAA(I,JCOLK) IF(I.EQ.IROWK) GO TO 18 AAA(I,JCOLK)=-AIJCK/PIVOT DO 17 J=1,MAX 17 IF(J.NE.JCOLK) AAA(I,J)=AAA(I,J)-AIJCK*AAA(IROWK,J) 18 CONTINUE C ORDER SOLUTION VALUES (IF ANY) AND CREATE JORD ARRAY DO 20 I=1,NDIM IROWI=IROW(I) JCOLI=JCOL(I) JORD(IROWI)=JCOLI 20 IF(INDIC.GE.0) V(JCOLI)=AAA(IROWI,MAX) C ADJUST SIGN OF DETERMINANT INTCH=0 NM1=NDIM-1 C PRINT *,'NM1', NM1,'N',N,'NDIM',NDIM DO 22 I=1,NM1 IP1=I+1 IF(DFO2.EQ.-11.) WRITE(*,*) NM1,NDIM DO 22 J=IP1,NDIM IF(JORD(J).GE.JORD(I)) GO TO 22 JTEMP=JORD(J) JORD(J)=JORD(I) JORD(I)=JTEMP INTCH=INTCH+1 22 CONTINUE C PRINT *,'AFTER 22' IF(INTCH/2*2.NE.INTCH) DETER=-DETER C IF INDIC POSITIVE RETURN WITH RESULTS IF (INDIC.LE.0) GO TO 26 SIMUL=DETER RETURN C IF INDIC IS NEGATIVE OR ZERO, UNSCRAMBLE THE INVERSE 26 DO 28 J=1,NDIM DO 27 I=1,NDIM IROWI=IROW(I) JCOLI=JCOL(I) 27 Y(JCOLI)=AAA(IROWI,J) DO 28 I=1,NDIM 28 AAA(I,J)=Y(I) C PRINT *,'AFTER 28' DO 30 I=1,NDIM DO 29 J=1,NDIM IROWJ=IROW(J) JCOLJ=JCOL(J) 29 Y(IROWJ)=AAA(I,JCOLJ) DO 30 J=1,NDIM 30 AAA(I,J)=Y(J) C PRINT *,'AFTER 30' C RETURN FOR INDIC NEGATIVE OR ZERO SIMUL=DETER RETURN 200 FORMAT(10H0N TOO BIG ) END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE PTDER(XA,YA,TA,DGDTK,DGDP,GTTK,GPT,GTKX,GTKY, 1GPY,GPX,IOPX,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) IF(IOPX.EQ.1) THEN DG0DT=DZ(20) DG0DP=0.D0 DF0DT=DZ(6) DF0DP=DZ(8) U01DT=DZ(23) U02DT=DZ(24) U04DT=DZ(25) U01DP=DZ(26) U02DP=DZ(27) U04DP=DZ(28) DADP1=0.D0 DBDP1=0.D0 DCDP1=0.D0 DDDP1=0.D0 A1=0.D0 B1=0.D0 C1=H(35) D1=H(36) ELSE DG0DT=DZ(21) DG0DP=DZ(22) DF0DT=DZ(5) DF0DP=DZ(7) U01DT=DZ(10) U02DT=DZ(11) U04DT=DZ(12) U01DP=DZ(13) U02DP=DZ(14) U04DP=DZ(15) DADP1=DZ(16) DBDP1=DZ(17) DCDP1=DZ(18) DDDP1=DZ(19) A1=A B1=B C1=C D1=D END IF DM1=XA+YA+TA/2.D0 DM2=2.D0*XA-DM1 FE1=1.D0-DM1 C2=2.D0*YA FE2=1.D0-DM2-C2 G00=XA*(1.D0-XA-YA)+YA*TA/2.D0+TA*TA/4.D0 F00=YA*DM1 S=-R*(DM1*DLOG(DM1)+DM2*DLOG(DM2)+FE1*DLOG(FE1)+ 1 FE2*DLOG(FE2)+C2*DLOG(C2)) DGDTK=DG0DT*G00+DF0DT*F00+2.D0*YA*(U04DT-U02DT) 1 -S+XA*(U01DT-U02DT)+U02DT DGDP=DG0DP*G00+DF0DP*F00+2.D0*YA*(U04DP-U02DP) 1 +XA*(U01DP-U02DP)+U02DP 1 +YA*(DCDP1*(XA-TA/2.D0)-DDDP1*YA+DBDP1*(2.D0*YA*TA-4.D0*XA*YA) 1 +DADP1*(XA*XA+3.D0*YA*YA+TA*TA/4.D0-XA*TA)) 1 +.5D0*YA*(DDDP1-DCDP1-2.D0*(DADP1-DBDP1)) GTTK=.5D0*(DG0DT*(YA+TA)+DF0DT*YA+R*DLOG(DM1*FE2/FE1/DM2)) GPT=.5D0*(DG0DP*(YA+TA)+DF0DP*YA+YA*(DADP1*(TA-2.D0*XA) 1 +4.D0*YA*DBDP1-DCDP1)) GTKX=DG0DT*(1.D0-2.D0*XA-YA)+DF0DT*YA+R*DLOG(DM1*DM2/FE1/FE2) 1 +U01DT-U02DT GTKY=DG0DT*(-XA+TA/2.D0)+DF0DT*(YA+DM1) 1 +R*DLOG(C2*C2*DM1/FE1/FE2/DM2) 1 +2.D0*(U04DT-U02DT) GPY=DG0DP*(-XA+TA/2.D0)+DF0DP*(YA+DM1) 1 +DCDP1*(XA-TA/2.D0)-2.D0*YA*DDDP1+DBDP1*(4.D0*YA*TA-8.D0*XA*YA) 1 +DADP1*(XA*XA+9.D0*YA*YA+TA*TA/4.D0-XA*TA) 1 +.5D0*(DDDP1-DCDP1-2.D0*(DADP1-DBDP1)) 1 +2.D0*(U04DP-U02DP) GPX=DG0DP*(1.D0-2.D0*XA-YA)+YA*(DF0DP+DCDP1+DADP1*(2.D0*XA-TA) 1 -4.D0*YA*DBDP1)+U01DP-U02DP RETURN END C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE MU0(XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) PARAMETER (NCOMP=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION UO(4),GU(4),A1(4),B1(4),C1(4),S(4),DU0DT(4),D2U0DT(4), 1 V(4),DELUO(4) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) DATA TR1,TR2,T9/903.15D0,413.15D0,298.15D0/ C FERROSILITE STD STATES REVISED AUG,85 BY DHL DATA GU(1),GU(2),GU(3),GU(4)/-16393.D0,-4183.D0,-34256.D0, 1 -26400.D0/ DATA S(1),S(2),S(3),S(4)/32.4D0,45.8711D0,34.2D0,40.7D0/ DATA A1(1),A1(2),A1(3),A1(4)/49.1D0,55.6866D0,52.87D0,54.81D0/ DATA B1(1),B1(2),B1(3),B1(4)/9.48D-03,7.26688D-03,7.84D-03, 1 8.17D-03/ DATA C1(1),C1(2),C1(3),C1(4)/12.56D05,15.447D05,15.74D05,15.01D05/ DATA FALPHA,FBETA/39.3D-06,1.01D-06/ C DATA GU(1),GU(2),GU(3),GU(4)/-16393.D0,-4834.D0,-34256.D0, C 1 -26400.D0/ C DATA S(1),S(2),S(3),S(4)/32.4D0,45.2D0,34.2D0,40.7D0/ C DATA A1(1),A1(2),A1(3),A1(4)/49.1D0,52.98D0,52.87D0,54.81D0/ C DATA B1(1),B1(2),B1(3),B1(4)/9.48D-03,10.14D-03,7.84D-03,8.17D-03/ C DATA C1(1),C1(2),C1(3),C1(4)/12.56D05,11.1D05,15.74D05,15.01D05/ C] DATA V(1),V(2),V(3),V(4)/6.2875D0,6.63686D0,6.62D0,6.788D0/ DATA V(1),V(2),V(3),V(4)/6.2552D0,6.60160D0,6.609D0,6.827D0/ C UO(I) IS MU0 FOR CEN,CFS,DI,HD: 1,2,3,4. GU(4)=H(48)/4.184D0 DO 10 I=1,4 T=H(30) P2=H(29) IF(I.EQ.1) THEN T=TR1 P2=1.D0 END IF IF(I.EQ.2) THEN C T=TR2 C AUG 85 FERROSILITE IS OFS (WITH COMPRESSIBILITY+THERMAL EXPANS)! P2=1.D0 END IF UO(I)=GU(I)-S(I)*(T-T9)+A1(I)*(T-T9-T*DLOG(T/T9))+ 1 (C1(I)-B1(I)*T*T9*T9)*(T-T9)*(T-T9)/(2.D0*T*T9*T9) UO(I)=UO(I)*4.184D0+V(I)*(P2-1.D0) 10 CONTINUE DO 15 I=2,4 DU0DT(I)=4.184D0*(-S(I)-A1(I)*DLOG(H(30)/T9)+C1(I)*(H(30)*H(30)- 1T9*T9)/(2.D0*T9*T9*H(30)*H(30))-B1(I)*(H(30)-T9)) 15 D2U0DT(I)=4.184D0*(-A1(I)/H(30)-B1(I)+C1(I)/H(30)/H(30)/H(30)) C REM U0'S IN J]]] D2U0DT(1)=-8.368D0*28.765D0/H(30) C D2U0DT(2)=-8.368D0*21.D0/H(30)-9.D-03 S903=S(1)+.386D0+A1(1)*DLOG(TR1/T9)+B1(1)*(TR1-T9)+C1(1)/2.D0 1 *(1.D0/TR1/TR1-1.D0/T9/T9) UO(1)=UO(1)+8.368D0*28.765D0*(H(30)-TR1-H(30)*DLOG(H(30)/TR1)) 1 +6.2592D0*(H(29)-1.D0)-4.184D0*S903*(H(30)-TR1) 1 +3561.D0-1.91D0*H(30)+.0355D0*H(29) VTOFS=V(2)*(1.D0+FALPHA*(H(30)-T9)) UO(2)=UO(2)+VTOFS*((1.D0+FBETA)*(H(29)-1.D0) 1 -.5D0*FBETA*(H(29)*H(29)-1.D0)) 1 +1843.D0-1.6760D0*H(30)+.03686D0*H(29) DU0DT(1)=8.368D0*(28.765D0*DLOG(TR1/H(30))-S903/2.D0)-1.91D0 C DU0DT(2)=4.184D0*(A1(1)*DLOG(TR2/H(30))-B1(1)*(H(30)-TR2)) C 1 -1.6496D0 H(25)=2.D0*(UO(3)-UO(4))+UO(2)-UO(1) C F0V=2*(6.62-6.788)+6.63686-6.2875=.0134 C] H(25)=H(25)+H(29)*.0134D0 CC DZ(5)=2.D0*(DU0DT(3)-DU0DT(4))+DU0DT(2)-DU0DT(1) C ABOVE DOES NOT USE NEW FS STD STATE!--OCT 85 H(37)=H(25)+33552.D0-17.614D0*H(30)-.01002D0*H(29) C FOLLOWING IS FOR ALTERNATE HD-FS MODEL C] H(37)=H(25)+39648.D0-18.622D0*H(30)-.00831D0*H(29) CC DZ(6)=DZ(5)-17.614D0 DZ(9)=2.D0*(D2U0DT(3)-D2U0DT(4))+D2U0DT(2)-D2U0DT(1) C DATA FROM DI-EN,HD-FS DEL-G ENDMEMBERS C DEL MU0'S FROM OPX-CPX; U0'S NO LONGER SCALED TO U0OPX=0 FOR EN,FS,& HD C MU(OPX)=MU(CPX)-(MU(CPX)-MU(OPX)) C DELUO(1)=-3561.D0+1.91D0*H(30)-.0355D0*H(29) DELUO(2)=-1843.D0+1.676D0*H(30)-.03686D0*H(29) DELUO(4)=-5261.D0-.53D0*H(30)-.0951D0*H(29) DELUO(1)=-DELUO(1) DELUO(2)=-DELUO(2) H(40)=UO(1)-DELUO(1) H(41)=UO(2)-DELUO(2) H(42)=UO(4)-DELUO(4) H(26)=UO(1) H(27)=UO(2) H(28)=UO(4) DZ(33)=DU0DT(1) DZ(34)=DU0DT(2) DZ(35)=DU0DT(3) DZ(36)=DU0DT(4) DZ(38)=DU0DT(1)+1.91D0 DZ(39)=DU0DT(2)+1.676D0 DZ(40)=DU0DT(3)-8.16D0 DZ(41)=DU0DT(4)+0.53D0 DZ(37)=2.D0*(DZ(35)-DZ(36))+DZ(34)-DZ(33) DZ(42)=2.D0*(DZ(40)-DZ(41))+DZ(39)-DZ(38) RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ULIQ(INDEXA,XL,UL,DULDT,DULDX, 1 TP,H,DZ,A,B,C,D,GT1,R) PARAMETER(NCOMP=10) C C CALCULATION OF MU LIQUID = DEL MU (INDEXA PHASE-LIQUID) C + RT (LN A),LIQUID FOR SPECIFIED COMPONENTS C NIELSEN-DUNGAN MODEL (1983) UPDATED AUG 86 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XL(NCOMP),UL(NCOMP),DULDT(NCOMP),DULDX(NCOMP),TP(4) DIMENSION H(100),DZ(100) C XMG=XL(1) XFE=XL(2) XCA=XL(3) XSI=XL(4) DULDX(1)=2.D0*H(24)/XMG DULDX(3)=H(24)/XMG DULDX(4)=-H(24)/XFE DULDX(2)=2.*DULDX(4) C GO TO (10,20) INDEXA C C CPX-LIQ 10 RLN1=R*DLOG(XMG*XMG*XSI*XSI) RLN2=R*DLOG(XFE*XFE*XSI*XSI) RLN3=R*DLOG(XMG*XCA*XSI*XSI) RLN4=R*DLOG(XFE*XCA*XSI*XSI) UL(1)=H(30)*RLN1-H(1)-H(30)*DZ(10)-H(29)*DZ(13)-H(12)/H(30) UL(2)=H(30)*RLN2-H(2)-H(30)*DZ(11)-H(29)*DZ(14)-H(13)/H(30) UL(3)=H(30)*RLN3-H(10)-H(30)*DZ(30)-H(29)*DZ(32)-H(14)/H(30) UL(4)=H(30)*RLN4-H(3)-H(30)*DZ(12)-H(29)*DZ(15)-H(15)/H(30) DULDT(1)=-DZ(10)+RLN1+H(12)/H(30)/H(30) DULDT(2)=-DZ(11)+RLN2+H(13)/H(30)/H(30) DULDT(3)=-DZ(30)+RLN3+H(14)/H(30)/H(30) DULDT(4)=-DZ(12)+RLN4+H(15)/H(30)/H(30) RETURN C C OPX-LIQ 20 RLN1=R*DLOG(XMG*XMG*XSI*XSI) RLN2=R*DLOG(XFE*XFE*XSI*XSI) RLN3=R*DLOG(XMG*XCA*XSI*XSI) RLN4=R*DLOG(XFE*XCA*XSI*XSI) UL(1)=H(30)*RLN1-H(5)-H(30)*DZ(23)-H(29)*DZ(26)-H(16)/H(30) UL(2)=H(30)*RLN2-H(6)-H(30)*DZ(24)-H(29)*DZ(27)-H(17)/H(30) UL(3)=H(30)*RLN3-H(11)-H(30)*DZ(29)-H(29)*DZ(31)-H(18)/H(30) UL(4)=H(30)*RLN4-H(7)-H(30)*DZ(25)-H(29)*DZ(28)-H(19)/H(30) DULDT(1)=-DZ(23)+RLN1+H(16)/H(30)/H(30) DULDT(2)=-DZ(24)+RLN2+H(17)/H(30)/H(30) DULDT(3)=-DZ(29)+RLN3+H(18)/H(30)/H(30) DULDT(4)=-DZ(25)+RLN4+H(19)/H(30)/H(30) RETURN END