!-----Version: 11.02.2015 ! *********** ! * THERIAQ * ! *********** ! ! Program written by Christian de Capitani ! at the Department of Geological Sciences, ! University of British Columbia, Vancouver, B.C. Canada ! (May 1984 - Sept 1987) ! ! revision: April 1987 ! minor changes: December 1987 ! revision: July 1993 ! addition of aqueous phases: July 1995 ! ! for details of algorithm see: ! de Capitani C. and Brown T.H. : The computation of chemical ! equilibrium in complex systems containing non-ideal solutions. ! Geochim. Cosmochim. Acta 51(1987):2639-2652 ! ! Any suggestions, complaints or comments are greatly appreciated ! by the author and should be sent to: ! Christian de Capitani ! Mineralogisch-Petrographisches Institut ! Universitaet Basel ! Bernoullistrasse 30 ! CH-4056 BASEL ! ! IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !-----END OF COMMON VARIABLES INTEGER*4 MAXBL PARAMETER (MAXBL=10) !----- CHARACTER*32 CH16 REAL*8 X16,BULBL(MAXBL,COMAX),TCBL(MAXBL),PBL(MAXBL), & FLUBL(MAXBL,COMAX),SOLBL(MAXBL,COMAX),KGWBL(MAXBL) INTEGER*4 IS,NBEO,FLOWMOD,NBL,IBL CHARACTER*25 USEBL(MAXBL) !-----END OF VARIABLES FOR MEGAMAT REAL*8 XKGH2O,F1,F2 INTEGER*4 AQFAIL,NAQLOOP,VERYSPEZ !----- LOGICAL*4 MORE INTEGER*4 I001,I002,I,II,I1,COMAY REAL*8 FF,SOLBUL(COMAX),LABUL(COMAX),INFLOW(COMAX) CHARACTER*500 CH001,CH002,SYREC,CHIN(3),ZEITSTRING CHARACTER*16 TEXT(COMAX) CHARACTER*8 CH8 LOGICAL*4 ISASOLID(COMAX) !----- REAL*8 TCBK,PBK,FLUBK(COMAX),SOLBK(COMAX),KGWBK,PHBK,PEBK CHARACTER*25 USEBK !---- END OF VARIABLES FOR AQ-EQUIL. !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- INTEGER*4 NRLOOP,DHPAR REAL*8 NRFIDI,NRMUZ,NRRED COMMON /NERAIN/ NRLOOP,DHPAR COMMON /NERARE/ NRFIDI,NRMUZ,NRRED !---- INTEGER*4 ierr,j progname='THERIAQ' vers='11.02.2015' task='"Computation of aqueous equilibriua"' call initialize('$THERIAQ-FILES',ierr) if(ierr.ne.0) STOP !***** DO 400,I=1,3 400 CHIN(I)=' ' !***** !------------------ ! Open UNIT=log !------------------ j=log line=filename(j)(1:fnl(j))//ext(j) path=wpath akzess=' ' state=' ' call openfile(j,ierr) if(ierr.ne.0) STOP !----- READ (UNIT=log,FMT='(I6,2X,E12.5,2X,E12.5,2X,E12.5,I3)',END=411) & NRLOOP,NRFIDI,NRMUZ,NRRED,DHPAR DO 410,I=1,3 410 READ (UNIT=log,FMT='(A500)',END=411) CHIN(I) 411 CONTINUE WRITE (6,100) 100 FORMAT (/ & ' -------------------'/ & ' database definition'/ & ' -------------------') CALL LABLA(CHIN(1),I002) IF (I002.EQ.0) I002=1 CH002=' Enter [ "?" | CR | "files" | database filename ] <'// & CHIN(1)(1:I002)//'>?' !----- 412 CONTINUE CALL PUST (6,CH002) READ (5,FMT='(A500)') CH001 IF (CH001.EQ.'?') THEN CALL helpme('$THK-START') GOTO 412 END IF IF (VERGL(CH001,'files')) THEN CALL listfiles GOTO 412 END IF IF (CH001.EQ.' ') THEN CH001=CHIN(1) I001=I002 ELSE CHIN(1)=CH001 END IF CH002=CH001 CALL TAXI(CH002,DBNAME) CALL LABLA(DBNAME,I001) !------------------ ! Open UNIT=out !------------------ j=out line=filename(j)(1:fnl(j))//ext(j) path=wpath akzess=' ' state=' ' call openfile(j,ierr) if(ierr.ne.0) STOP !------------------ ! open UNIT=dbs !------------------ j=dbs line=DBNAME path=wpath akzess=' ' state='old' call openfile(j,ierr) if(ierr.ne.0) STOP !------------------ ! open UNIT=dat !------------------ j=dat line=filename(j)(1:fnl(j))//ext(j) path=wpath akzess=' ' state='old' call openfile(j,ierr) if(ierr.ne.0) STOP !----- WRITE (scr,101) DBNAME(1:I001) WRITE (out,101) DBNAME(1:I001) 101 FORMAT (/' database for this run: ',A/) !***** COMAY=COMAX CALL PROREAD(SYREC) Call LABLA(filename(dat),I1) WRITE (scr,102) filename(dat)(1:I1) WRITE (out,102) filename(dat)(1:I1) 102 FORMAT (/,' Input from file',1x,a) write(scr,104) ('-',I=1,I1+16) write(out,104) ('-',I=1,I1+16) 104 format(1x,130a1,:) write(scr,106) TC,P write(out,106) TC,P 106 format(' T =',F8.2,' C P =',F9.2,' Bar') CALL PUST (scr,' '//SYREC) CALL PUST (out,' '//SYREC) !----- !-----READ PRTCOD, FORMUL AND USE FROM SYREC !-----SET UP FIRST NUN COLUMNS OF MATRIX CALL GELI(SYREC,FF) PRTCOD=IDINT(FF) IF (PRTCOD.EQ.0) TEST=DABS(TEST) DO 650,I=1,11 650 PRTLOG(I)=.FALSE. IF (PRTCOD.LE.-2) PRTLOG(1)=.TRUE. IF (PRTCOD.EQ.-1) THEN DO 652,I=2,5 652 PRTLOG(I)=.TRUE. END IF IF (PRTCOD.EQ.0) THEN DO 654,I=5,6 654 PRTLOG(I)=.TRUE. END IF IF (PRTCOD.GE.1) THEN DO 656,I=2,8 656 PRTLOG(I)=.TRUE. END IF ! IF (PRTLOG(1)) STOP ! CLOSE(UNIT=dat) IF (PRTLOG(5)) THEN WRITE (UNIT=scr,FMT=140) TEST,LO1MAX,EQUALX,DXMIN WRITE (UNIT=out,FMT=140) TEST,LO1MAX,EQUALX,DXMIN 140 FORMAT (/,' TEST =',1PE11.4,8X,'LO1MAX =',I4,13X,'EQUALX =', & 1PE11.4,6X,'DELXMIN =',1PE11.4) WRITE (UNIT=scr,FMT=141) DXSCAN,DXSTAR,STPSTA,STPMAX,GCMAX WRITE (UNIT=out,FMT=141) DXSCAN,DXSTAR,STPSTA,STPMAX,GCMAX 141 FORMAT (' DELXSCAN =',1PE11.4,4X,'DELXSTAR =',1PE11.4,4X, & 'STEPSTAR =',I4,11X,'STEPMAX =',I4,12X,'GCMAX =',I5) END IF !----- CALL TAXI(SYREC,FORMUL) CALL TAXI(SYREC,USE) !+++++ CALL CHEMIE(COMAY,NC,OXYDE,OXANZ,FORMUL,CHEM) CALL LABLA(USE,LUSE) !----------------------------------------------------------------- ! NVAR=1 ! VARIS(1)='loop' !================================================================= !------------------ ! Open UNIT=ibk !------------------ j=ibk line=filename(j)(1:fnl(j))//ext(j) path=wpath akzess=' ' state=' ' call openfile(j,ierr) if(ierr.ne.0) STOP DO 810,I=1,NC 810 INFLOW(I)=0.0D0 10 READ (UNIT=ibk,FMT='(A500)',END=888) CH001 CALL TAXI(CH001,CH8) CALL GELI(CH001,FF) I001=0 DO 800,I=1,NC IF (CH8.EQ.OXYDE(I)) THEN I001=I GOTO 11 END IF 800 CONTINUE 11 IF (I001.EQ.0) THEN WRITE (UNIT=scr,FMT='(''Unknown element:'',A8)') CH8 STOP END IF INFLOW(I001)=FF GOTO 10 888 CONTINUE ! WRITE (UNIT=scr,FMT=160) (OXYDE(I),INFLOW(I),I=1,NC) ! 160 FORMAT (A8,1PE12.5) !================================================================= !================================================================= !---- input in file "blocks" !---- FLOWMOD: 0: one flow direction, [1: flow in both directions (not implemented)] !---- for each block: TCBL: Temperature, PBL: Pressure, KGWBL: Kg of water !---- bulk: BULBL(i), i=1,NC, USEBL: use-variable NBL=0 OPEN (UNIT=40,FILE='blocks',STATUS='UNKNOWN') READ (UNIT=40,FMT='(A500)',END=77) CH001 !C CALL PUST (6,CH001) CALL GELI(CH001,FF) FLOWMOD=IDINT(FF) 7 READ (UNIT=40,FMT='(A500)',END=77) CH001 !C CALL PUST (6,CH001) CALL GELI(CH001,F1) CALL GELI(CH001,F2) IF (F1.EQ.0.0D0.OR.F2.EQ.0.0D0) GOTO 77 NBL=NBL+1 TCBL(NBL)=F1 PBL(NBL)=F2 CALL GELI(CH001,FF) KGWBL(NBL)=FF READ (UNIT=40,FMT='(A500)',END=77) CH001 !C CALL PUST (6,CH001) CALL GELI(CH001,F1) !C WRITE (*,*) F1 CALL TAXI(CH001,FORMUL) !C CALL PUST (6,FORMUL) CALL TAXI(CH001,CH) !C WRITE (*,*) CH USEBL(NBL)=CH CALL CHEMIE(COMAY,NC,OXYDE,OXANZ,FORMUL,CHEM) DO 700,I=1,NC BULBL(NBL,I)=CHEM(I) !C WRITE (*,*) OXYDE(I),CHEM(I) 700 CONTINUE GOTO 7 77 CONTINUE CLOSE (UNIT=40) DO 702,II=1,NBL WRITE (scr,1000) II WRITE (out,1000) II 1000 FORMAT (//' block:',I3/' ---------') WRITE (scr,1004) TCBL(II),PBL(II),USEBL(II),KGWBL(II) WRITE (out,1004) TCBL(II),PBL(II),USEBL(II),KGWBL(II) 1004 FORMAT (' T =',F8.2,' C',' P =',F9.2' Bar', & ' use =',A8,' KgW =',1PE15.8) DO I=1,NC IF (BULBL(II,I).NE.0.0D0) THEN WRITE (scr,1006) OXYDE(I),BULBL(II,I) WRITE (out,1006) OXYDE(I),BULBL(II,I) 1006 FORMAT (1X,A10,F11.5) END IF END DO 702 CONTINUE !================================================================= !---- input: NAQLOOP: number of main loops, XKGH2O: Kg water flowing per loop VERYSPEZ=0 WRITE (UNIT=6,FMT='('' '')') CALL LABLA(CHIN(2),I002) IF (I002.EQ.0) I002=1 CH002='number of loops ? ('//CHIN(2)(1:I002)//'):' CALL PUST (6,CH002) READ (5,FMT='(A500)') CH001 IF (CH001.EQ.' ') THEN CH001=CHIN(2) I001=I002 ELSE CHIN(2)=CH001 END IF CALL GELI(CH001,FF) NAQLOOP=IDINT(FF) IF (NAQLOOP.EQ.0) VERYSPEZ=1 CALL LABLA(CHIN(3),I002) IF (I002.EQ.0) I002=1 CH002='kg of water per loop ? ('//CHIN(3)(1:I002)//'):' CALL PUST (6,CH002) READ (5,FMT='(A500)') CH001 IF (CH001.EQ.' ') THEN CH001=CHIN(3) I001=I002 ELSE CHIN(3)=CH001 END IF CALL GELI(CH001,XKGH2O) !-----store terminal input CLOSE (UNIT=log) ! OPEN (UNIT=log,FILE='theriaq.last',STATUS='OLD') !------------------ ! Open UNIT=log !------------------ j=log line=filename(j)(1:fnl(j))//ext(j) path=wpath akzess=' ' state='OLD' call openfile(j,ierr) if(ierr.ne.0) STOP WRITE (UNIT=log,FMT='(I6,2X,E12.5,2X,E12.5,2X,E12.5,I3)') & NRLOOP,NRFIDI,NRMUZ,NRRED,DHPAR DO 145,I=1,3 145 CALL PUST(log,CHIN(I)) CLOSE (UNIT=log) !================================================================= !CCC PRTLOG(9)=.TRUE. PRTLOG(9)=.FALSE. PRTLOG(2)=.TRUE. NVARTBL=0 NROWTBL=0 !--- IF (VERYSPEZ.EQ.1) THEN CALL SPEZSUB(AQFAIL) GOTO 711 END IF !---- SOLBL(ibl,): Bulk of solids in blodk !---- FLUBL(ibl,): Bulk in fluid in blodk !---- initial: all FLUBL(ibl,i)=0.0D0 DO 707,IBL=1,NBL DO 707,I=1,NC SOLBL(IBL,I)=BULBL(IBL,I) 707 FLUBL(IBL,I)=0.0D0 !----------------------------------------------------------------- !--- start main aqueous loop !----------------------------------------------------------------- DO 710,NBEO=1,NAQLOOP !----------------------------------------------------------------- !--- start block loop !----------------------------------------------------------------- DO 720,IBL=1,NBL WRITE (UNIT=6,FMT='(/'' aqueous loop:'',I4)') NBEO WRITE (UNIT=out,FMT='(/'' aqueous loop:'',I4)') NBEO WRITE (UNIT=6,FMT='('' block: '',I4)') IBL WRITE (UNIT=out,FMT='('' block: '',I4)') IBL !C DO 722,I=1,NC !C 722 FLUBL(IBL,I)=0.0D0 !----- NROWTBL=NROWTBL+1 IF (NROWTBL.GT.MAXTBL) THEN WRITE (UNIT=6,FMT='(''too many rows: calculation stopped'')') WRITE (UNIT=out,FMT='(''too many rows: calculation stopped'')') STOP END IF DO 600,I=1,NVARTBL 600 OUTTBL(NROWTBL,I)=0.0D0 !----- !----- CH16='loop' X16=DBLE(NBEO) CALL SETTBL(CH16,X16) CH16='block' X16=DBLE(IBL) CALL SETTBL(CH16,X16) !----- !============================================================= !---- for one block: SOLBK: bulk of solids, FLUBK: bulk in fluid !---- TCBK: temperatuer, PBK: pressure, KGWBK: kg water !---- PEBK: pE, PHBK: pH, USEBK: USE DO 723,I=1,NC SOLBK(I)=SOLBL(IBL,I) 723 FLUBK(I)=FLUBL(IBL,I) TCBK=TCBL(IBL) PBK=PBL(IBL) KGWBK=KGWBL(IBL) PEBK=0.0D0 PHBK=0.0D0 USEBK=USEBL(IBL) CALL SEAQUIL(TCBK,PBK,FLUBK,SOLBK,USEBK,KGWBK,PHBK,PEBK) DO 725,I=1,NC SOLBL(IBL,I)=SOLBK(I) FLUBL(IBL,I)=FLUBK(I) 725 BULBL(IBL,I)=SOLBK(I)+FLUBK(I) !C STOP !C GOTO 710 !============================================================= !----- 720 CONTINUE !----------------------------------------------------------------- !--- end block loop !----------------------------------------------------------------- !----- mix water !----- DO 750,IBL=NBL,1,-1 WRITE (6,2010) IBL WRITE (out,2010) IBL 2010 FORMAT (/' fluid composition (FLUBL, BULBL) block:',I5) DO 752,I=1,NC IF (FLUBL(IBL,I).NE.0.0D0.OR.BULBL(IBL,I).NE.0.0D0) THEN WRITE (6,2012) OXYDE(I),FLUBL(IBL,I),BULBL(IBL,I) WRITE (out,2012) OXYDE(I),FLUBL(IBL,I),BULBL(IBL,I) 2012 FORMAT (1X,A8,2(2X,1PE15.8)) END IF 752 CONTINUE !== IF (IBL.GT.1) THEN DO 754,I=1,NC IF (OXYDE(I).NE.'O '.AND.OXYDE(I).NE.'H '.AND.OXYDE(I).NE.'E ') & THEN FF=FLUBL(IBL,I)/KGWBL(IBL) FLUBL(IBL,I)=FLUBL(IBL,I)-XKGH2O*FF IF (FLUBL(IBL,I).LT.0.0D0) FLUBL(IBL,I)=0.0D0 FF=FLUBL(IBL-1,I)/KGWBL(IBL-1) FLUBL(IBL,I)=FLUBL(IBL,I)+XKGH2O*FF !CCC IF (BULBL(IBL,I).LT.0.0D0) BULBL(IBL,I)=0.0D0 END IF 754 CONTINUE END IF !== IF (IBL.EQ.1) THEN DO 755,I=1,NC IF (OXYDE(I).NE.'O '.AND.OXYDE(I).NE.'H '.AND.OXYDE(I).NE.'E ') & THEN FF=FLUBL(IBL,I)/KGWBL(IBL) FLUBL(IBL,I)=FLUBL(IBL,I)-XKGH2O*FF IF (FLUBL(IBL,I).LT.0.0D0) FLUBL(IBL,I)=0.0D0 END IF 755 CONTINUE END IF !== WRITE (6,2014) WRITE (out,2014) 2014 FORMAT (/' new FLUBL') DO 756,I=1,NC IF (FLUBL(IBL,I).NE.0.0D0) THEN WRITE (6,2016) OXYDE(I),FLUBL(IBL,I) WRITE (out,2016) OXYDE(I),FLUBL(IBL,I) 2016 FORMAT (1X,A8,(2X,1PE15.8)) END IF 756 CONTINUE 750 CONTINUE !----- 710 CONTINUE !----------------------------------------------------------------- !--- end main aqueous loop !----------------------------------------------------------------- !---- the last call to DBREAD is apparently not necessary (??) 711 CALL PRTTBL !------------------ ! Open UNIT=lbk !------------------ j=lbk line=filename(j)(1:fnl(j))//ext(j) path=wpath akzess=' ' state=' ' call openfile(j,ierr) if(ierr.ne.0) STOP WRITE (UNIT=lbk,FMT='(''*** LAST BULK COMPOSITION ***'')') DO 770,I=1,NUN WRITE (UNIT=lbk,FMT=1114) CHNAME(I),LABUL(I) 1114 FORMAT (1X,A8,1X,1PE12.5) 770 CONTINUE WRITE (UNIT=lbk,FMT='(''*** AQUEOUS COMPOSITION ***'')') DO 772,I=1,NUN WRITE (UNIT=lbk,FMT=1115) OXYDE(CHMCOD(I)),QSUM(I) 1115 FORMAT (2X,A8,1PE12.5) 772 CONTINUE CLOSE(UNIT=lbk) !----------------------------------------------------------------- !----------------------------------------------------------------- !-----AT THIS POINT THE PROGRAM MAY BE CHANGED TO PERFORM REPEATED !-----CALCULATIONS WITH VARYING T,P OR BULK COMPOSITION !----------------------------------------------------------------- !-----CALCULATE THE EQUILIBRIUM ASSEMBLAGE: !***** ! CALL CALSTR ! CALL PRININ (if you want anything printed before calculation) ! CALL THERIA !***** !----------------------------------------------------------------- !-----CHANGE T OR P !***** ! TC=..... (Temperature in deg. C) ! P=..... (Pressure in Bars) ! CALL NURVONPT !***** !----------------------------------------------------------------- !-----CHANGE THE BULK COMPOSITION: !-----(THIS IS SLIGHTLY MORE COMPLICATED. THE FOLLOWING IS AN EXAMPLE.) !***** ! READ (UNIT=5,FMT='(A170)') FORMUL ! CALL CHEMIE(COMAY,NC,OXYDE,OXANZ,FORMUL,CHE) ! MORE=.FALSE. ! DO 601,I=1,NC ! IF (CHE(I).EQ.0.0D0.NEQV.CHEM(I).EQ.0.0D0) MORE=.TRUE. ! CHEM(I)=CHE(I) ! 601 CONTINUE ! IF (MORE) THEN ! CALL DBREAD ! CALL NURVONPT ! ELSE ! DO 602,I=1,NUN ! 602 BULK(I)=CHE(CHMCOD(I)) ! END IF !***** !----------------------------------------------------------------- !-----CONTROL THE AMOUNT OF OUTPUT PRODUCED: !-----IF PRTLOG(n) IS SET .FALSE. THEN THE CORRESPONDING OUTPUT !-----IS OMITTED. !***** !*****BEFORE THE CALCULATION: !***** ! PRTLOG(1)=.TRUE. (stop after reading database) ! PRTLOG(2)=.TRUE. (Print bulk composition) ! PRTLOG(3)=.TRUE. (Print list of considered phases) ! PRTLOG(4)=.TRUE. (Print a summary of the solution models) ! PRTLOG(5)=.TRUE. (Print the parameters) !***** !*****AFTER THE CALCULATION: !***** ! PRTLOG(6)=.TRUE.OR.PRTLOG(7)=.TRUE.OR.PRTLOG(8)=.TRUE. ! (Print stable assemblage) ! PRTLOG(6)=.TRUE. (Print volumes and densities) ! PRTLOG(7)=.TRUE. (Print compositions of all stable phases) ! PRTLOG(8)=.TRUE. (Print activities of all considered phases) ! PRTLOG(9)=.TRUE. (used to print table in loop) ! PRTLOG(10)=.TRUE. (used to print image for thermap) ! PRTLOG(11)=.TRUE. (used to print short table (e.g. theriaq)) !***** !----------------------------------------------------------------- ! CALL NURVONPT ! CALL CALSTR ! IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ ! CALL THERIA CALL CPUTIME(ZEITSTRING) CALL LABLA(ZEITSTRING,I001) WRITE (scr,150) ZEITSTRING(1:I001) WRITE (out,150) ZEITSTRING(1:I001) 150 FORMAT (/,'exit THERIAK',/,A) END !----- !******************************** SUBROUTINE SEAQUIL3(TCBK,PBK,FLUBK,SOLBK,USEBK,KGWBK,PHBK,PEBK) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !-----END OF COMMON VARIABLES !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- !---- REAL*8 TCBK,PBK,FLUBK(COMAX),SOLBK(COMAX),KGWBK,PHBK,PEBK, & BMDIS(COMAX),AFUER(PHMAX,PHMAX),CFUER(PHMAX),DFUER(PHMAX), & NEWNN(COMAX),BMDISNEW(COMAX) CHARACTER*16 TEXT CHARACTER*25 USEBK CHARACTER*32 CH16 CHARACTER*125 AS0,AS1 !---- END OF VARIABLES FOR AQ-EQUIL. REAL*8 FF,TOTBK(COMAX),AQPOT0(COMAX),POTDIFF,FBULK(COMAX), & BULK0(COMAX),NNPOT(COMAX),MFBULK(COMAX),PROJG,X16,MOLPKG, & XFLU(COMAX),MUEFLU(COMAX),GFLU,SUM,AWAT,IONSUM,AWAT0, & MUWAT,MUWAT0,MUKOR INTEGER*4 I,II,AQFAIL,AQPHAS(COMAX),NDLOOP,NDLMAX,IDPOT(COMAX), & NFUER,ERFUER,K,ISNUL,IH,IO,IH2O,INH,INO,INH2O,IXH2O LOGICAL*4 FIXPOT(COMAX),FINIS !----- MOLPKG=55.508435062D0 PRTEST=.FALSE. FINIS=.FALSE. NDLMAX=10 AS0=' ' AS1=' ' IF (.NOT.PRTEST) THEN DO 200,I=1,11 200 PRTLOG(I)=.FALSE. END IF !============================================================= !---- first: TOTBK: total bulk of solids plus fluids !---- AQPHAS: for each component a phase with G=0 (same as "ELEMENTS") !---- AQPHAS(i)=NUMMER of phase for component i !---- for ndloop=1: aqueous phases: NULL=.true. (not used) !---- SOLBK(i) solid bulk of block (1 to NC) !---- FLUBK(i) fluid bulk of block (1 to NC) !---- TOTBK(i) total bulk of Block (solid+fluid) (1 to NC) !---- !---- BULK(i) bulk used for THERIAK, is initially set in DBREAD from CHEM (1 to NUN) !---- BULK0(i) total bulk of of block (1 to NUN) !---- FBULK(i) fluid bulk of block (1 to NUN) !---- MFBULK(i) maximum bulk for fluid (usually = BULK0) AWAT0=1.0D0 MUWAT0=0.0D0 DO 300,I=1,NC TOTBK(I)=SOLBK(I)+FLUBK(I) FLUBK(I)=0.0D0 300 CHEM(I)=TOTBK(I) USE=USEBK CALL LABLA(USE,LUSE) CALL DBREAD CALL POTPHAS(AQPHAS) DO 302,I=1,NPHA IF (PHASID(I).EQ.'AQU'.OR.PHASID(I).EQ.'AQP') THEN NULL(I)=.TRUE. END IF 302 CONTINUE DO 304,I=1,NUN BMDIS(I)=0.0D0 IDPOT(I)=0 IF (CHNAME(I).EQ.'H ') BULK(I)=2.0D0*MOLPKG BULK0(I)=BULK(I) MFBULK(I)=BULK(I) FBULK(I)=0.0D0 NNPOT(I)=0.0D0 FIXPOT(I)=.FALSE. 304 AQPOT0(I)=0.0D0 INH=0 INO=0 INH2O=0 DO 310,I=1,NPHA IF (ABK(I).EQ.'H2O') INH2O=I IF (NAME(I).EQ.'"H"') INH=I IF (NAME(I).EQ.'"O"') INO=I 310 CONTINUE NDLOOP=0 !==================================================================== !---- ndloop=1: call THERIAK with total bulk (solid + fluid) !==================================================================== WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (equil.)'')') NDLOOP TC=TCBK P=PBK CALL NURVONPT CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ !================ CALL THERIA !================ CALL ANASS(AS0) WRITE (6,1000) WRITE (out,1000) 1000 FORMAT (/' initial assemblage (total bulk):') CALL PUST(6,' '//AS0) CALL PUST(out,' '//AS0) NPOT=0 PROJG=0.0D0 POTDIFF=0.0D0 !==================================================================== 1 CONTINUE NDLOOP=NDLOOP+1 WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (aq. 1)'')') NDLOOP !====================== CALL AQEQ(AQFAIL,IONSUM) IF (AQFAIL.EQ.1) RETURN !====================== !==================================================================== !---- set NULL of water to .true. ?? !---- define fluid solution in calc structure !==================================================================== IH=0 IO=0 IH2O=0 DO 350,I=1,NMAX IF (NUMMER(I).GT.0) THEN IF (NUMMER(I).EQ.INH2O) IH2O=I IF (NUMMER(I).EQ.INH) IH=I IF (NUMMER(I).EQ.INO) IO=I END IF 350 CONTINUE MUKOR=G(IH2O)-2.0D0*G(IH)-G(IO) WRITE (UNIT=6,FMT='(/,'' g of H2O'',2X,1PE12.5)') G(IH2O) WRITE (UNIT=6,FMT='('' g of H '',2X,1PE12.5)') G(IH) WRITE (UNIT=6,FMT='('' g of O '',2X,1PE12.5)') G(IO) !C NULL(NUMMER(IH2O))=.TRUE. !----- change GGK to calculated potential WRITE (UNIT=6,FMT='(/,'' calculated m'')') SUM=0.0D0 DO 400,I=1,NUN XFLU(I)=QSUM(I) IF (CHNAME(I).EQ.'H ') XFLU(I)=XFLU(I)+2.0D0*MOLPKG IF (CHNAME(I).EQ.'O ') XFLU(I)=XFLU(I)+MOLPKG DO 401,K=1,NMAX IF (NUMMER(K).EQ.I) MUEFLU(I)=G(K) 401 CONTINUE WRITE (6,2000) CHNAME(I),XFLU(I),MUEFLU(I) WRITE (out,2000) CHNAME(I),XFLU(I),MUEFLU(I) 2000 FORMAT (1X,A8,4(2X,1PE12.5)) IF (XFLU(I).GT.0.0D0) SUM=SUM+XFLU(I) 400 CONTINUE !--- WRITE (UNIT=6,FMT='(/,'' calculated x'')') GFLU=0.0D0 DO 405,I=1,NUN IF (XFLU(I).GT.0.0D0) THEN !c XFLU(I)=XFLU(I)/SUM ELSE XFLU(I)=0.0D0 END IF GFLU=GFLU+XFLU(I)*MUEFLU(I) WRITE (6,2002) CHNAME(I),XFLU(I),MUEFLU(I) WRITE (out,2002) CHNAME(I),XFLU(I),MUEFLU(I) 2002 FORMAT (1X,A8,4(2X,1PE12.5)) 405 CONTINUE WRITE (6,2004) GFLU WRITE (out,2004) GFLU 2004 FORMAT (/,' g of fluid:',4(2X,1PE12.5)) CALL ADDAQ(XFLU,NDLOOP) AWAT=1.0D0-0.017D0*IONSUM/55.508435062D0 AWAT=1.0D0-IONSUM/(IONSUM+55.508435062D0) MUWAT=RT*DLOG(AWAT) WRITE (UNIT=6,FMT='(/,'' a of water'',2X,1PE12.5)') AWAT GGK(NPHA)=-GFLU !C GGK(NPHA)=-GFLU-2000.0d0 !C GGK(NPHA)=-GFLU-MUKOR !C GGK(NPHA)=GGK(NPHA)+MUKOR*MOLPKG FF=MUWAT WRITE (UNIT=6,FMT='(/,'' mue of water'',2X,1PE12.5)') FF FF=MUWAT*MOLPKG WRITE (UNIT=6,FMT='('' g correction'',2X,1PE12.5)') FF CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ !================ CALL THERIA !================ WRITE (UNIT=6,FMT='(/,'' 1=continue, 0=end'')') READ (*,*) I IF (I.EQ.1) GOTO 1 !----------------------------------------------------------------- !***** END !----- !******************************** SUBROUTINE SEAQUIL(TCBK,PBK,FLUBK,SOLBK,USEBK,KGWBK,PHBK,PEBK) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !-----END OF COMMON VARIABLES !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- REAL*8 TCBK,PBK,FLUBK(COMAX),SOLBK(COMAX),KGWBK,PHBK,PEBK, & BMDIS(COMAX),AFUER(PHMAX,PHMAX),CFUER(PHMAX),DFUER(PHMAX), & NEWNN(COMAX),BMDISNEW(COMAX) CHARACTER*16 TEXT CHARACTER*25 USEBK CHARACTER*32 CH16 CHARACTER*125 AS0,AS1 !---- END OF VARIABLES FOR AQ-EQUIL. REAL*8 FF,TOTBK(COMAX),AQPOT0(COMAX),POTDIFF,FBULK(COMAX), & BULK0(COMAX),NNPOT(COMAX),MFBULK(COMAX),PROJG,X16,IONSUM INTEGER*4 I,II,AQFAIL,AQPHAS(COMAX),NDLOOP,NDLMAX,IDPOT(COMAX), & NFUER,ERFUER,K,ISNUL LOGICAL*4 FIXPOT(COMAX),FINIS !----- PRTEST=.FALSE. FINIS=.FALSE. NDLMAX=10 AS0=' ' AS1=' ' IF (.NOT.PRTEST) THEN DO 200,I=1,11 200 PRTLOG(I)=.FALSE. PRTLOG(2)=.TRUE. END IF !============================================================= !---- first: TOTBK: total bulk of solids plus fluids !---- AQPHAS: for each component a phase with G=0 (same as "ELEMENTS") !---- AQPHAS(i)=NUMMER of phase for component i !---- for ndloop=1: aqueous phases: NULL=.true. (not used) !---- SOLBK(i) solid bulk of block (1 to NC) !---- FLUBK(i) fluid bulk of block (1 to NC) !---- TOTBK(i) total bulk of Block (solid+fluid) (1 to NC) !---- !---- BULK(i) bulk used for THERIAK, is initially set in DBREAD from CHEM (1 to NUN) !---- BULK0(i) total bulk of of block (1 to NUN) !---- FBULK(i) fluid bulk of block (1 to NUN) !---- MFBULK(i) maximum bulk for fluid (usually = BULK0) DO 300,I=1,NC TOTBK(I)=SOLBK(I)+FLUBK(I) FLUBK(I)=0.0D0 300 CHEM(I)=TOTBK(I) USE=USEBK CALL LABLA(USE,LUSE) CALL DBREAD CALL POTPHAS(AQPHAS) DO 302,I=1,NPHA IF (PHASID(I).EQ.'AQU'.OR.PHASID(I).EQ.'AQP') THEN NULL(I)=.TRUE. END IF 302 CONTINUE DO 304,I=1,NUN BMDIS(I)=0.0D0 IDPOT(I)=0 BULK0(I)=BULK(I) MFBULK(I)=BULK(I) FBULK(I)=0.0D0 NNPOT(I)=0.0D0 FIXPOT(I)=.FALSE. 304 AQPOT0(I)=0.0D0 NDLOOP=0 !==================================================================== !---- ndloop=1: call THERIAK with total bulk (solid + fluid) !==================================================================== WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (equil.)'')') NDLOOP TC=TCBK P=PBK WRITE (UNIT=6,FMT='(/,'' TC= '',F8.2,'' P ='',F8.2)') TC,P CALL NURVONPT CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ !================ CALL THERIA !================ CALL ANASS(AS0) WRITE (6,1000) WRITE (out,1000) 1000 FORMAT (/' initial assemblage (total bulk):') CALL PUST(6,' '//AS0) CALL PUST(out,' '//AS0) IF (.NOT.PRTEST) CALL PRTSHORT NPOT=0 PROJG=0.0D0 POTDIFF=0.0D0 !-- DO I=1,NMAX DO II=1,NUN IF (NUMMER(I).EQ.AQPHAS(II)) THEN AQPOT(II)=G(I) END IF END DO END DO !-- 1 CONTINUE IF (NDLOOP.GT.NDLMAX) GOTO 999 !==================================================================== DO II=1,NUN AQPOT0(II)=AQPOT(II) END DO NDLOOP=NDLOOP+1 WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (aq. 1)'')') NDLOOP !====================== CALL AQEQ(AQFAIL,IONSUM) !====================== !==================================================================== !----- check if dissolved > bulk !==================================================================== IF (PRTEST) WRITE (6,2000) KGWBK IF (PRTEST) WRITE (out,2000) KGWBK 2000 FORMAT (/' bulk and aqueous compositions (moles)', & 27X,'mass of water: ',F12.6,' Kg', & /15X,'bulk',10X,'fluid',9X,'max.fl',4X,'(fl-max)/max') NPOT=0 DO 400,I=1,NUN FF=QSUM(I)*KGWBK IF (CHNAME(I).NE.'O'.AND.CHNAME(I).NE.'H'.AND.CHNAME(I).NE.'E') & THEN IF (PRTEST) WRITE (6,2002) CHNAME(I),BULK(I),FF,MFBULK(I), & (FF-MFBULK(I))/MFBULK(I) IF (PRTEST) WRITE (out,2002) CHNAME(I),BULK(I),FF,MFBULK(I), & (FF-MFBULK(I))/MFBULK(I) 2002 FORMAT (1X,A8,4(2X,1PE12.5)) END IF IF ((FF-MFBULK(I))/MFBULK(I).GT.1D-5.AND.CHNAME(I).NE.'E') THEN NPOT=NPOT+1 POT(NPOT)=I CPOT(NPOT)=MFBULK(I) END IF 400 CONTINUE !==================================================================== !----- if dissolved > bulk fix compositions call AQEQ,THERIA !==================================================================== IF (NPOT.GT.0) THEN WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (aq. 2)'')') NDLOOP !====================== CALL AQEQ(AQFAIL,IONSUM) !====================== IF (AQFAIL.EQ.1) RETURN !----- change GGK to calculated potential DO 500,I=1,NMAX DO 502,II=1,NPOT IF (NUMMER(I).EQ.AQPHAS(POT(II))) THEN NULL(AQPHAS(POT(II)))=.FALSE. GGK(AQPHAS(POT(II)))=GGK(AQPHAS(POT(II)))-G(I)-MUEG(POT(II)) AQPOT0(POT(II))=0.0D0 FIXPOT(POT(II))=.TRUE. WRITE (6,2010) NAME(AQPHAS(POT(II))),GGK(AQPHAS(POT(II))) WRITE (out,2010) NAME(AQPHAS(POT(II))),GGK(AQPHAS(POT(II))) 2010 FORMAT (/,1X,A16,' GGK is set to:',2X,1PE15.8) END IF 502 CONTINUE 500 CONTINUE !----- calculate fluid fbulk DO 510,I=1,NUN FF=QSUM(I)*KGWBK IF (CHNAME(I).NE.'E') THEN FBULK(I)=FF END IF 510 CONTINUE !----- define new solid bulk (with AQPOT) II=0 WRITE (UNIT=6,FMT='(/''BULK,FBULK,NNPOT'')') DO 512,I=1,NUN IF (FIXPOT(I)) THEN BULK(I)=BULK0(I) ELSE BULK(I)=BULK0(I)-FBULK(I) IF (BULK(I).LE.0.0D0) BULK(I)=0.0D0 END IF ! ! IF (PRTEST) WRITE (6,2015) CHNAME(I),BULK(I),FBULK(I),NNPOT(I) WRITE (6,2015) CHNAME(I),BULK(I),FBULK(I),NNPOT(I) WRITE (out,2015) CHNAME(I),BULK(I),FBULK(I),NNPOT(I) 2015 FORMAT (9X,A8,3(2X,1PE15.8)) 512 CONTINUE CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ !================ CALL THERIA !================ CALL ANASS(AS0) WRITE (6,2020) WRITE (out,2020) 2020 FORMAT (/'new initial assemblage (with AQPOT):') CALL PUST(6,AS0) CALL PUST(out,AS0) IF (.NOT.PRTEST) CALL PRTSHORT NPOT=0 PROJG=0.0D0 POTDIFF=0.0D0 !-- DO I=1,NMAX DO II=1,NUN IF (NUMMER(I).EQ.AQPHAS(II)) THEN AQPOT0(II)=G(I) END IF END DO END DO !-- WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (aq. 3)'')') NDLOOP !====================== CALL AQEQ(AQFAIL,IONSUM) !====================== !----- END IF !==================================================================== !----- calculate fbulk and substract from bulk0 !==================================================================== DO 610,I=1,NUN FF=QSUM(I)*KGWBK IF (CHNAME(I).NE.'E') THEN FBULK(I)=FF END IF 610 CONTINUE !----- define new solid bulk WRITE (UNIT=6,FMT='(/'' BULK,FBULK'')') DO 612,I=1,NUN IF (FIXPOT(I)) THEN BULK(I)=BULK0(I) ELSE BULK(I)=BULK0(I)-FBULK(I) IF (BULK(I).LE.0.0D0) BULK(I)=0.0D0 END IF ! IF (PRTEST) WRITE (6,3000) CHNAME(I),BULK(I),FBULK(I) WRITE (6,3000) CHNAME(I),BULK(I),FBULK(I) WRITE (out,3000) CHNAME(I),BULK(I),FBULK(I) 3000 FORMAT (9X,A8,3(2X,1PE15.8)) 612 CONTINUE CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ !================ CALL THERIA !================ CALL ANASS(AS1) WRITE (6,3010) WRITE (out,3010) 3010 FORMAT (/' assemblages:') CALL PUST(6,' '//AS0) CALL PUST(out,' '//AS0) CALL PUST(6,' '//AS1) CALL PUST(out,' '//AS1) IF (.NOT.PRTEST) CALL PRTSHORT !-- DO I=1,NMAX DO II=1,NUN IF (NUMMER(I).EQ.AQPHAS(II)) THEN AQPOT(II)=G(I) END IF END DO END DO !-- !==================================================================== !----- check if ass0 = ass1 !==================================================================== IF (AS0.NE.AS1) THEN WRITE (6,3012) WRITE (out,3012) 3012 FORMAT (' assemblage changes after dissolution') GOTO 999 END IF !==================================================================== !----- check chemical potentials !==================================================================== PROJG=0.0D0 DO I=1,NUN PROJG=PROJG+AQPOT(I)*BULK0(I) END DO WRITE (6,3020) PROJG WRITE (out,3020) PROJG 3020 FORMAT (/,' projected G = ',1PE15.8) WRITE (6,3022) WRITE (out,3022) 3022 FORMAT (/' chemical potentials',/ & 23X,'previous',12X,'now',9X,'difference') POTDIFF=0.0D0 DO I=1,NUN FF=DABS(AQPOT(I)-AQPOT0(I)) IF (FF.GT.POTDIFF) POTDIFF=FF WRITE (6,3024) CHNAME(I),AQPOT0(I),AQPOT(I), & AQPOT(I)-AQPOT0(I) WRITE (out,3024) CHNAME(I),AQPOT0(I),AQPOT(I), & AQPOT(I)-AQPOT0(I) 3024 FORMAT (9X,A8,3(2X,1PE15.8)) END DO WRITE (6,3026) POTDIFF WRITE (out,3026) POTDIFF 3026 FORMAT (/,' potdiff = ',1PE15.8) !================ ! WRITE (UNIT=6,FMT='(/,'' 1=continue, 0=end'')') ! READ (*,*) I ! IF (I.EQ.1) GOTO 1 IF (POTDIFF.GT.1D-5) GOTO 1 !----- print last equilibrium 999 CONTINUE PRTLOG(6)=.TRUE. PRTLOG(8)=.TRUE. PRTLOG(11)=.TRUE. CALL PRTCAL PRTLOG(6)=.FALSE. PRTLOG(8)=.FALSE. PRTLOG(11)=.FALSE. PRTEST=.TRUE. CALL AQEQ(AQFAIL,IONSUM) PRTEST=.FALSE. WRITE (6,1008) WRITE (out,1008) 1008 FORMAT (/' equilibrium compositions',/ & 24X,'solids',11X,'fluid',12X,'total') DO I=1,NUN WRITE (6,1010) CHNAME(I),BULK(I),FBULK(I),BULK0(I) WRITE (out,1010) CHNAME(I),BULK(I),FBULK(I),BULK0(I) 1010 FORMAT (9X,A8,3(2X,1PE15.8)) CH16='(blk)'//CHNAME(I)(1:8) X16=BULK(I) CALL SETTBL(CH16,X16) END DO !----- before returning: set FLUBK and SOLBK DO 330,I=1,NC SOLBK(I)=0.0D0 330 FLUBK(I)=0.0D0 DO 332,I=1,NUN SOLBK(CHMCOD(I))=BULK(I) 332 FLUBK(CHMCOD(I))=FBULK(I) !***** END !----- !******************************** SUBROUTINE SEAQUIL2(TCBK,PBK,FLUBK,SOLBK,USEBK,KGWBK,PHBK,PEBK) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !-----END OF COMMON VARIABLES !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- REAL*8 TCBK,PBK,FLUBK(COMAX),SOLBK(COMAX),KGWBK,PHBK,PEBK CHARACTER*25 USEBK CHARACTER*32 CH16 CHARACTER*125 AS0,AS1 !---- END OF VARIABLES FOR AQ-EQUIL. REAL*8 FF,TOTBK(COMAX),AQPOT0(COMAX),POTDIFF,FBULK(COMAX), & BULK0(COMAX),NNPOT(COMAX),MFBULK(COMAX),PROJG,X16,IONSUM INTEGER*4 I,II,AQFAIL,AQPHAS(COMAX),NDLOOP,NDLMAX,IDPOT(COMAX) LOGICAL*4 FIXPOT(COMAX),FINIS !----- PRTEST=.TRUE. FINIS=.FALSE. NDLMAX=10 AS0=' ' AS1=' ' IF (.NOT.PRTEST) THEN DO 200,I=1,11 200 PRTLOG(I)=.FALSE. END IF !============================================================= !---- first: TOTBK: total bulk of solids plus fluids !---- AQPHAS: for each component a phase with G=0 (same as "ELEMENTS") !---- AQPHAS(i)=NUMMER of phase for component i !---- for ndloop=1: aqueous phases: NULL=.true. (not used) !---- SOLBK(i) solid bulk of block (1 to NC) !---- FLUBK(i) fluid bulk of block (1 to NC) !---- TOTBK(i) total bulk of Block (solid+fluid) (1 to NC) !---- !---- BULK(i) bulk used for THERIAK, is initially set in DBREAD from CHEM (1 to NUN) !---- BULK0(i) total bulk of of block (1 to NUN) !---- FBULK(i) fluid bulk of block (1 to NUN) !---- MFBULK(i) maximum bulk for fluid (usually = BULK0) DO 300,I=1,NC TOTBK(I)=SOLBK(I)+FLUBK(I) FLUBK(I)=0.0D0 300 CHEM(I)=TOTBK(I) USE=USEBK CALL LABLA(USE,LUSE) CALL DBREAD CALL POTPHAS(AQPHAS) DO 302,I=1,NPHA IF (PHASID(I).EQ.'AQU'.OR.PHASID(I).EQ.'AQP') THEN NULL(I)=.TRUE. END IF 302 CONTINUE DO 304,I=1,NUN IDPOT(I)=0 BULK0(I)=BULK(I) MFBULK(I)=BULK(I) FBULK(I)=0.0D0 NNPOT(I)=0.0D0 FIXPOT(I)=.FALSE. 304 AQPOT0(I)=0.0D0 NDLOOP=0 TC=TCBK P=PBK CALL NURVONPT !============================================================= !---- ndloop=1: call THERIAK with total bulk (solid + fluid) !==== hier zurueck 1 CONTINUE NDLOOP=NDLOOP+1 IF (NDLOOP.GT.NDLMAX) RETURN WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (equil.)'')') NDLOOP WRITE (UNIT=out,FMT='(/,'' loop = '',I3,'' (equil.)'')') NDLOOP IF (PRTEST) WRITE (6,3004) IF (PRTEST) WRITE (out,3004) 3004 FORMAT (/' pre-theriak compositions',/ & 24X,'solids',11X,'fluid',12X,'nnpot') DO 305,I=1,NUN IF (FIXPOT(I)) THEN BULK(I)=BULK0(I) ELSE BULK(I)=BULK0(I)-FBULK(I) END IF IF (PRTEST) WRITE (6,3006) CHNAME(I),BULK(I),FBULK(I),NNPOT(I) IF (PRTEST) WRITE (out,3006) CHNAME(I),BULK(I),FBULK(I),NNPOT(I) 3006 FORMAT (9X,A8,3(2X,1PE15.8)) 305 CONTINUE CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ !================ CALL THERIA !================ AS0=AS1 CALL ANASS(AS1) WRITE (6,3010) WRITE (out,3010) 3010 FORMAT (/'assemblages:') CALL PUST(6,' '//AS0) CALL PUST(out,' '//AS0) CALL PUST(6,' '//AS1) CALL PUST(out,' '//AS1) NPOT=0 DO 310,I=1,NMAX DO 310,II=1,NUN IF (NUMMER(I).EQ.AQPHAS(II)) THEN AQPOT(II)=G(I) NNPOT(II)=NN(I) IDPOT(II)=I IF (NNPOT(II).NE.0.0D0) MFBULK(II)=NN(I) END IF 310 CONTINUE PROJG=0.0D0 DO 312,I=1,NUN 312 PROJG=PROJG+AQPOT(I)*BULK0(I) WRITE (6,3008) PROJG WRITE (out,3008) PROJG 3008 FORMAT (/,' projected G = ',1PE15.8) IF (PRTEST) WRITE (6,1000) IF (PRTEST) WRITE (out,1000) 1000 FORMAT (/' chemical potentials',/ & 23X,'previous',12X,'now',9X,'difference') POTDIFF=0.0D0 DO 314,I=1,NUN FF=DABS(AQPOT(I)-AQPOT0(I)) IF (FF.GT.POTDIFF) POTDIFF=FF IF (PRTEST) WRITE (6,1002) CHNAME(I),AQPOT0(I),AQPOT(I), & AQPOT(I)-AQPOT0(I) IF (PRTEST) WRITE (out,1002) CHNAME(I),AQPOT0(I),AQPOT(I), & AQPOT(I)-AQPOT0(I) 1002 FORMAT (9X,A8,3(2X,1PE15.8)) !++ ! IF (NDLOOP.GE.2.AND.DABS(AQPOT(I)-AQPOT0(I)).GT.1D-3) THEN ! GGK(AQPHAS(I))=GGK(AQPHAS(I))-(AQPOT(I)-AQPOT0(I))/2.0D0 ! NULL(AQPHAS(I))=.FALSE. !CC AQPOT0(I)=0.0D0 ! FIXPOT(I)=.TRUE. ! END IF !++ AQPOT0(I)=AQPOT(I) IF ((FBULK(I)-MFBULK(I))/MFBULK(I).GT.1D-4) THEN NPOT=NPOT+1 POT(NPOT)=I CPOT(NPOT)=MFBULK(I) END IF 314 CONTINUE IF (PRTEST) THEN WRITE (6,1004) WRITE (out,1004) 1004 FORMAT (/' temporary compositions',/ & 24X,'solids',11X,'fluid',12X,'max.fl') DO 316,I=1,NUN WRITE (6,1006) CHNAME(I),BULK(I),FBULK(I),MFBULK(I) WRITE (out,1006) CHNAME(I),BULK(I),FBULK(I),MFBULK(I) 1006 FORMAT (9X,A8,3(2X,1PE15.8)) 316 CONTINUE END IF !------------- calculation will stop after next aq.eq. IF (POTDIFF.LT.1D-4.AND.NPOT.EQ.0 & .OR.NDLOOP.GE.6.AND.NPOT.EQ.0) THEN IF (POTDIFF.GE.1D-4) THEN WRITE (6,1009) WRITE (out,1009) 1009 FORMAT (/' change in assemblage') END IF FINIS=.TRUE. PRTLOG(6)=.TRUE. PRTLOG(11)=.TRUE. CALL PRTCAL PRTLOG(6)=.FALSE. PRTLOG(11)=.FALSE. PRTEST=.TRUE. END IF !----------------------------------------------------------------- WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (aq. 1)'')') NDLOOP !====================== CALL AQEQ(AQFAIL,IONSUM) !====================== IF (AQFAIL.EQ.1) RETURN !----- change GGK to calculated potential DO 500,I=1,NMAX DO 502,II=1,NPOT IF (NUMMER(I).EQ.AQPHAS(POT(II))) THEN NULL(AQPHAS(POT(II)))=.FALSE. GGK(AQPHAS(POT(II)))=GGK(AQPHAS(POT(II)))-G(I)-MUEG(POT(II)) AQPOT0(POT(II))=0.0D0 FIXPOT(POT(II))=.TRUE. WRITE (6,2015) NAME(AQPHAS(POT(II))),GGK(AQPHAS(POT(II))) WRITE (out,2015) NAME(AQPHAS(POT(II))),GGK(AQPHAS(POT(II))) 2015 FORMAT (/,1X,A16,' GGK is set to:',2X,1PE15.8) END IF 502 CONTINUE 500 CONTINUE !----------------------------------------------------------------- !----- if dissolved > bulk fix compositions and call AQEQ IF (PRTEST) WRITE (6,2012) KGWBK IF (PRTEST) WRITE (out,2012) KGWBK 2012 FORMAT (/' bulk and aqueous compositions (moles)', & 27X,'mass of water: ',F12.6,' Kg', & /15X,'bulk',10X,'fluid',9X,'max.fl',4X,'(fl-max)/max') NPOT=0 DO 766,I=1,NUN FF=QSUM(I)*KGWBK IF (CHNAME(I).NE.'O'.AND.CHNAME(I).NE.'H'.AND.CHNAME(I).NE.'E') & THEN IF (PRTEST) WRITE (6,2013) CHNAME(I),BULK(I),FF,MFBULK(I), & (FF-MFBULK(I))/MFBULK(I) IF (PRTEST) WRITE (out,2013) CHNAME(I),BULK(I),FF,MFBULK(I), & (FF-MFBULK(I))/MFBULK(I) 2013 FORMAT (1X,A8,4(2X,1PE12.5)) END IF IF ((FF-MFBULK(I))/MFBULK(I).GT.1D-5.AND.CHNAME(I).NE.'E') THEN NPOT=NPOT+1 POT(NPOT)=I CPOT(NPOT)=MFBULK(I) END IF 766 CONTINUE IF (NPOT.GT.0) THEN WRITE (UNIT=6,FMT='(/,'' loop = '',I3,'' (aq. 2)'')') NDLOOP !====================== CALL AQEQ(AQFAIL,IONSUM) !====================== IF (AQFAIL.EQ.1) RETURN END IF !----------------------------------------------------------------- !----- calculate effective fbulk DO 780,I=1,NUN FF=QSUM(I)*KGWBK IF (CHNAME(I).NE.'E') THEN FBULK(I)=FF END IF 780 CONTINUE !----- change GGK to calculated potential DO 510,I=1,NMAX DO 512,II=1,NPOT IF (NUMMER(I).EQ.AQPHAS(POT(II))) THEN NULL(AQPHAS(POT(II)))=.FALSE. GGK(AQPHAS(POT(II)))=GGK(AQPHAS(POT(II)))-G(I)-MUEG(POT(II)) AQPOT0(POT(II))=0.0D0 FIXPOT(POT(II))=.TRUE. WRITE (6,2025) NAME(AQPHAS(POT(II))),GGK(AQPHAS(POT(II))) WRITE (out,2025) NAME(AQPHAS(POT(II))),GGK(AQPHAS(POT(II))) 2025 FORMAT (/,1X,A16,' GGK is set to:',2X,1PE15.8) END IF 512 CONTINUE 510 CONTINUE !===== IF (FINIS) THEN DO 318,I=1,NUN BULK(I)=BULK0(I)-FBULK(I) IF (BULK(I).LT.0.0D0) BULK(I)=0.0D0 318 CONTINUE WRITE (6,1008) WRITE (out,1008) 1008 FORMAT (/' equilibrium compositions',/ & 24X,'solids',11X,'fluid',12X,'total') DO 320,I=1,NUN WRITE (6,1010) CHNAME(I),BULK(I),FBULK(I),BULK0(I) WRITE (out,1010) CHNAME(I),BULK(I),FBULK(I),BULK0(I) 1010 FORMAT (9X,A8,3(2X,1PE15.8)) CH16='(blk)'//CHNAME(I)(1:8) X16=BULK(I) CALL SETTBL(CH16,X16) 320 CONTINUE DO 330,I=1,NC SOLBK(I)=0.0D0 330 FLUBK(I)=0.0D0 DO 332,I=1,NUN SOLBK(CHMCOD(I))=BULK(I) 332 FLUBK(CHMCOD(I))=FBULK(I) RETURN END IF !===== GOTO 1 !----- RETURN END !----- !******************************** SUBROUTINE AQEQ(AQFAIL,IONSUM) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !-----END OF COMMON VARIABLES !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- REAL*8 LOGI(PHMAX),MOLW(PHMAX),DUM3(PHMAX),XM(PHMAX), & IONST,IONSUM,DUM1,DUM2,XRES(PHMAX),CC0(PHMAX), & JACOB(PHMAX,PHMAX),DX0(PHMAX),DXG(PHMAX), & LNMOL(PHMAX),LNXM(PHMAX),FF,AH2O,FINDIFF,AVONP,BVONP,TFURP,F1 INTEGER*4 I,IR,IC,N0,IH2O,ERRC,AQFAIL,IMEGA REAL*8 RH2O,AL,BE,DALDT,DBEDT,DBEDP,EE,DEDP,DEDT,DETT CHARACTER*32 CH16 CHARACTER*132 CH132 !----- INTEGER*4 NRLOOP,DHPAR REAL*8 NRFIDI,NRMUZ,NRRED,ZWERT COMMON /NERAIN/ NRLOOP,DHPAR COMMON /NERARE/ NRFIDI,NRMUZ,NRRED !---- DO 510,I=1,NPOT WRITE (6,1001) CHNAME(POT(I)),CPOT(I) WRITE (out,1001) CHNAME(POT(I)),CPOT(I) 1001 FORMAT (/,1X,A8,'set to ',1PE15.8) 510 CONTINUE WRITE (UNIT=CH132,FMT='(1X,131A1)') ('-',I=1,131) !---- AQFAIL=0 ZWERT=0.0D0 ! FINDIFF=1.0D-8 FINDIFF=NRFIDI IZ=0 ICH=0 NROW=NUN NCOL=0 DO 500,IR=1,NROW RNAME(IR)=CHNAME(IR) IF (RNAME(IR)(1:1).EQ.' ') THEN CH16=RNAME(IR)(2:) RNAME(IR)=CH16 END IF IF (RNAME(IR).EQ.'E') IZ=IR 500 CONTINUE IF (IZ.EQ.0) THEN WRITE (UNIT=6,FMT='(/''no component E'')') WRITE (UNIT=out,FMT='(/''no component E'')') RETURN END IF !+++++ !+++++Copy all information into AA etc. !+++++ DO 550,IC=1,NUN IF (SUGG(IC).LE.NUN) THEN WRITE (6,1000) WRITE (out,1000) 1000 FORMAT (/' Assemblage not buffered, no aqueous equilibrium ', & 'is calculated') RETURN !CC STOP END IF 550 CONTINUE !----- IF (PRTEST) WRITE (6,1026) CH132 IF (PRTEST) WRITE (out,1026) CH132 1026 FORMAT (/,A,/ & ' aqueous equilibria (only species with G<1E5)'/ & ' ------------------') IH2O=0 IEP=0 DO 600,IC=1,NMAX IF (NUMMER(IC).NE.0) THEN IF (ABK(NUMMER(IC)).EQ.'H2O') IH2O=IC IF (ABK(NUMMER(IC)).EQ.'E') IEP=IC IF ((PHASID(NUMMER(IC)).EQ.'AQU' & .OR.PHASID(NUMMER(IC)).EQ.'AQP')) THEN IF (G(IC).LT.1D5.OR.NAME(NUMMER(IC)).EQ.'H+') THEN !CC >.AND.G(IC).LT.1D5) THEN !dC >.AND.G(IC).LT.1D20) THEN !----- NCOL=NCOL+1 CNAME(NCOL)=NAME(NUMMER(IC)) IF (CNAME(NCOL).EQ.'H+') ICH=NCOL !DC GW(NCOL)=G(IC) GV(NCOL)=G(IC) ZZ(NCOL)=X(IC,IZ) PHRA(NCOL)=INDIVI(NUMMER(IC),1) PHRB(NCOL)=INDIVI(NUMMER(IC),2) DO 610,IR=1,NROW AA(IR,NCOL)=X(IC,IR) 610 CONTINUE !----- ELSE IF (PRTEST) WRITE (6,1028) NAME(NUMMER(IC)),G(IC) IF (PRTEST) WRITE (out,1028) NAME(NUMMER(IC)),G(IC) 1028 FORMAT (' not used: ',A16,' G = ',1PE15.8) END IF END IF END IF 600 CONTINUE IF (ICH.EQ.0) THEN WRITE (6,1030) WRITE (out,1030) 1030 FORMAT (/'no H+, no calculations') RETURN END IF IF (IH2O.EQ.0) THEN WRITE (6,1032) WRITE (out,1032) 1032 FORMAT (/'no H2O, no calculations') RETURN ELSE AH2O=DEXP(-G(IH2O)/RT) END IF !+++++ AVONP=0.0D0 BVONP=0.0D0 TFURP=T CALL ABPHREEQ(TFURP,AVONP,BVONP) !+++++ L10=DLOG(10.0D0) CALL RHOETC(P,T,RH2O,AL,BE,DALDT,DBEDT,DBEDP) CALL JN91(T,RH2O,BE,AL,DALDT,EE,DEDP,DEDT,DETT) AGA=1.824829238D6*DSQRT(RH2O)/DSQRT((EE*T)**3) BGA=50.29158649D8*DSQRT(RH2O)/DSQRT(EE*T)*1D-8 IF (PRTEST) WRITE (6,1050) IF (PRTEST) WRITE (out,1050) 1050 FORMAT (/ & ' -----------------'/ & ' properties of H2O'/ & ' -----------------') IF (PRTEST) WRITE (UNIT=6,FMT=1051) P,PGAS,TC,T IF (PRTEST) WRITE (UNIT=out,FMT=1051) P,PGAS,TC,T 1051 FORMAT (' P =',F9.2,' bar',7X,'P(Gas) =',F9.2,' bar T =', & F8.2,' C =',F8.2,' K') IF (PRTEST) WRITE (6,1052) AH2O,RH2O,AL,BE,EE IF (PRTEST) WRITE (out,1052) AH2O,RH2O,AL,BE,EE 1052 FORMAT (/ & ' activity: ',7X,1PE12.5,7X,'density:',1PE12.5,' g/ml'/ & ' alpha: ',7X,1PE12.5,7X,'beta: ',1PE12.5/ & ' diel.const:',7X,0PF12.6) IF (PRTEST) WRITE (6,1053) AGA,BGA IF (PRTEST) WRITE (out,1053) AGA,BGA 1053 FORMAT (/ & ' Debye-Hueckel constants (SUPCRT): A =',1PE12.5,5X, & 'B =',1PE12.5) IF (PRTEST) WRITE (6,1054) AVONP,BVONP IF (PRTEST) WRITE (out,1054) AVONP,BVONP 1054 FORMAT ( & ' Debye-Hueckel constants (PHREEQ): A =',1PE12.5,5X, & 'B =',1PE12.5) IF (DHPAR.EQ.1) THEN AGA=AVONP BGA=BVONP END IF IF (PRTEST) WRITE (6,1055) AGA,BGA IF (PRTEST) WRITE (out,1055) AGA,BGA 1055 FORMAT (/ & ' Debye-Hueckel constants (USED) : A =',1PE12.5,5X, & 'B =',1PE12.5) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! nrow: number of components (=NUN) ! ncol: number of species (only phasid=AQU, and G<1D5) ! AA(nrow,ncol)=X(phase,comp) !!!!transpose!!! ! GV(i): chemical potential of species after THERIAK ! GW(i): chemical potential of species after AQEQ ! ZZ(i): charge of species ! IZ: index of component "E" ! IEP: index of phase "E" (in CALSTR) ! IH2O: index of phase "H2O" (in CALSTR) ! ICH: index of species "H+" ! variables: Ln(mi) [i=1,ncol] ! MUZ (correction of chem.pot. of charge) ! MUEG(i) (correction of chem.pot. of components with fixed concentrations) ! GW(i)=GV(i)+summe(AA(IR,IC)*MUEG(IR)) ! ! function F: i=1,ncol: ! F(i)=(GW(i)+ZZ(i)*MUZ)/RT+Ln(mi)+Ln(gammai) ! F(i)=(GW(I)+ZZ(I)*MUZ)/RT+DLOG(XM(I))+LOGI(I)*L10 ! i=ncol+1: ! total charge - value: ZSUM - ZWERT ! summe(mj*zj) (j=1,ncol) ! i>ncol+1: ! conc. of component - value: QSUM(j) - CPOT(j) ! summe(mj*AA(IR,IC)) (j=1,ncol) ! Newton-Raphson: ! F(xi)=0 ! F(xi + dx) = F(xi) + J*dx Jij=dFi/dxj ! solve: J*dx = -F ! xnew = xold + dx ! F: CC0(n0) ! J: JACOB(n0,n0) ! dx: DX0(n0) ! ! dFi/dxj: ! i,j=1 to ncol: finite difference of CC0 (subroutine RESID) ! ! i=ncol+1, j=1 to ncol: d(ZSUM-ZWERT)/dLn(xj) = ZZ(j)*MOLW(j) ! i=ncol+1, j=ncol+1: d(ZSUM-ZWERT)/dMUZ = 0 ! i=ncol+1, j=ncol+2 to n0: d(ZSUM-ZWERT)/dMUEGj = 0 ! ! i=ncol+1+ii, j=1 to ncol: d(QSUM-CPOT)ii/dLn(xj) = AA(POT(ii),j)*MOLW(j) ! i=ncol+1+ii, j=ncol+1: d(QSUM-CPOT)ii/dMUZ = 0 ! i=ncol+1+ii, j=ncol+2 to n0: d(QSUM-CPOT)ii/dMUEGj = 0 ! ! i=1 to ncol, j=ncol+1: dFi/dMUZ = ZZ(i)/RT ! i=1 to ncol, j=ncol+1+ii: dFi/dMUEGii = AA(POT(ii),i)/RT ! !---- The initial guess ! MUZ=0.0D0 MUZ=NRMUZ DO 410,IR=1,NUN 410 MUEG(IR)=0.0D0 DO 411,IC=1,NCOL 411 GW(IC)=GV(IC) DO 700,I=1,NCOL LNMOL(I)=(-GW(I)-ZZ(I)*MUZ)/RT MOLW(I)=DEXP(LNMOL(I)) 700 CONTINUE !===== start the megaloop DO 900,IMEGA=1,NRLOOP MEGA=IMEGA !----- the constant !----- CC0(1 to ncol)=(GW(I)+ZZ(I)*MUZ)/RT+DLOG(XM(I))+LOGI(I)*L10 DO 705,I=1,NCOL LNXM(I)=LNMOL(I) 705 XM(I)=DEXP(LNMOL(I)) CALL RESID(XM,IONST,IONSUM,LOGI,XRES) ZSUM=0.0D0 DO 710,I=1,NCOL CC0(I)=XRES(I) ZSUM=ZSUM+ZZ(I)*MOLW(I) 710 CONTINUE N0=NCOL+1 !----- CC0(ncol+1)=ZSUM-ZWERT CC0(N0)=ZSUM-ZWERT DO 715,IR=1,NROW QSUM(IR)=0.0D0 DO 720,IC=1,NCOL 720 QSUM(IR)=QSUM(IR)+AA(IR,IC)*MOLW(IC) 715 CONTINUE !----- CC0(ncol+1+i)=QSUM(POT(I))-CPOT(I) DO 721,I=1,NPOT N0=N0+1 CC0(N0)=QSUM(POT(I))-CPOT(I) 721 CONTINUE !DC ! WRITE (6,1042) MEGA,(CC0(I),I=1,N0) ! 1042 FORMAT ('CC0',I4,100(2X,1PE12.5)) ! WRITE (6,1043) MEGA,(GW(I),I=1,NCOL) ! WRITE (out,1043) MEGA,(GW(I),I=1,NCOL) ! 1043 FORMAT ('GW ',I4,100(2X,1PE12.5)) !---- test convergence FERTIG=0 DO 725,I=1,NCOL FF=DABS(CC0(I)*RT) IF (FF.GT.1D-3.AND.MOLW(I).GT.1D-50) THEN FERTIG=FERTIG+1 UNF(FERTIG)=I END IF !... DUM3(I)=FF 725 CONTINUE IF (DABS(ZSUM).GT.1.0D-8) THEN FERTIG=FERTIG+1 UNF(FERTIG)=NCOL+1 END IF DO 726,I=1,NPOT FF=DABS(CC0(NCOL+I+1)/CPOT(I)) IF (FF.GT.1D-5) THEN FERTIG=FERTIG+1 UNF(FERTIG)=NCOL+I+1 END IF 726 CONTINUE !... DUM3(NCOL+1)=ZSUM !... WRITE (UNIT=6,FMT='(''tst '',20E12.4)') (DUM3(I),I=1,NCOL) !CC CALL AQPRT(MOLW,LOGI,IONST,IONSUM) IF (FERTIG.EQ.0) GOTO 901 !----- !----- the derivatives !----- columns 1 to ncol !----- d(XRES)/d(Ln(m) DO 730,IC=1,NCOL LNXM(IC)=LNXM(IC)+FINDIFF XM(IC)=DEXP(LNXM(IC)) CALL RESID(XM,DUM1,DUM2,DUM3,XRES) DO 735,IR=1,NCOL 735 JACOB(IR,IC)=(XRES(IR)-CC0(IR))/(FINDIFF) !----- d(total charge)/d(Ln(mi) JACOB(NCOL+1,IC)=ZZ(IC)*MOLW(IC) !----- d(component)/d(Ln(mi) DO 736,I=1,NPOT JACOB(NCOL+I+1,IC)=AA(POT(I),IC)*MOLW(IC) 736 CONTINUE LNXM(IC)=LNMOL(IC) XM(IC)=MOLW(IC) 730 CONTINUE !----- column ncol+1 !----- d(XRES)/d(MUZ) DO 740,IR=1,NCOL 740 JACOB(IR,NCOL+1)=ZZ(IR)/RT !----- columns ncol+1 to n0 !----- d(XRES)/d(MUEGi) DO 750,IR=1,NCOL DO 750,I=1,NPOT JACOB(IR,NCOL+I+1)=AA(POT(I),IR)/RT 750 CONTINUE !----- columns ncol+1 to n0, and rows ncol+1 to n0 !----- d(total charge)/d(MUZ) and (components)/d(MUEGi) DO 753,IR=NCOL+1,N0 DO 753,IC=NCOL+1,N0 753 JACOB(IR,IC)=0.0D0 !---- print jacob ! DO 752,IR=1,N0 ! WRITE (UNIT=6,FMT='(''jac '',20E12.4)') (JACOB(IR,I),I=1,N0) ! 752 CONTINUE !---- !C CALL LINSOFT(N0,JACOB,CC0,DX0,ERRC) CALL LINSOFT2(N0,JACOB,CC0,DX0,ERRC) !----- IF (ERRC.EQ.1) THEN AQFAIL=1 DO 743,IR=1,N0 WRITE (UNIT=6,FMT='(''jac '',20E12.4)') (JACOB(IR,I),I=1,N0) WRITE (UNIT=out,FMT='(''jac '',20E12.4)') (JACOB(IR,I),I=1,N0) 743 CONTINUE END IF !---- !----check solution DO 741,IR=1,N0 DUM3(IR)=0.0D0 DO 741 IC=1,N0 DUM3(IR)=DUM3(IR)+DX0(IC)*JACOB(IR,IC) 741 CONTINUE DO 742,I=1,N0 FF=DABS(DUM3(I)-CC0(I)) IF (FF.GT.1D-8) THEN WRITE (6,2000) I,FF WRITE (out,2000) I,FF 2000 FORMAT (/' LINSOFT test>1D-8: index =',I4,' value =',1PE15.4) END IF 742 CONTINUE !----- IF (ERRC.EQ.1) THEN CALL AQPRT(MOLW,LOGI,IONST,IONSUM) RETURN END IF !----- DO 745,I=1,NCOL LNMOL(I)=LNMOL(I)-DX0(I)/NRRED IF (LNMOL(I).GT.10.0D0) LNMOL(I)=10.0D0 IF (LNMOL(I).LT.-200.0D0) LNMOL(I)=-200.0D0 MOLW(I)=DEXP(LNMOL(I)) 745 CONTINUE MUZ=MUZ-DX0(NCOL+1) DO 746,I=1,NPOT MUEG(POT(I))=MUEG(POT(I))-DX0(NCOL+I+1)/NRRED 746 CONTINUE !DC ! WRITE (6,1056) MEGA,(MUEG(I),I=1,NCOL) ! WRITE (out,1056) MEGA,(MUEG(I),I=1,NCOL) ! 1056 FORMAT ('MUG',I4,100(2X,1PE12.5)) DO 747,IC=1,NCOL GW(IC)=GV(IC) DO 747,IR=1,NROW 747 GW(IC)=GW(IC)+AA(IR,IC)*MUEG(IR) ! CALL AQPRT(MOLW,LOGI,IONST,IONSUM) 900 CONTINUE !===== end of megaloop 901 CONTINUE !C IF (PRTEST) CALL AQPRT(MOLW,LOGI,IONST,IONSUM) !C IF (PRTEST) WRITE (UNIT=6,FMT='(A132)') CH132 CALL AQPRT(MOLW,LOGI,IONST,IONSUM) WRITE (UNIT=6,FMT='(A132)') CH132 !----- RETURN END !----- !******************************** SUBROUTINE RESID(XM,IONST,IONSUM,LOGI,XRES) IMPLICIT NONE INCLUDE 'theriak.cmn' !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- REAL*8 XM(PHMAX),LOGI(PHMAX),IONST,IONSUM,XRES(PHMAX) INTEGER*4 I !---- IONST=0.0D0 IONSUM=0.0D0 DO 500,I=1,NCOL IONSUM=IONSUM+XM(I) IONST=IONST+XM(I)*ZZ(I)**2 500 CONTINUE IONST=IONST/2.0D0 DO 600,I=1,NCOL IF (ZZ(I).EQ.0.0D0) THEN ! LOGI(I)=0.0D0 IF (PHRB(I).EQ.0.0D0) THEN LOGI(I)=0.1D0*IONST ELSE LOGI(I)=PHRB(I)*IONST END IF ELSE !==== !---- Davies equation ? ! LOGI(I)=-AGA*ZZ(I)**2*(DSQRT(IONST) & ! /(1.0D0+DSQRT(IONST))+0.2D0*IONST) !---- Davies equation as in PHREEQ IF (PHRB(I).EQ.0.0D0.AND.PHRA(I).EQ.0.0D0) THEN LOGI(I)=-AGA*ZZ(I)**2*(DSQRT(IONST) & /(1.0D0+DSQRT(IONST))+0.3D0*IONST) ELSE LOGI(I)=-AGA*ZZ(I)**2*DSQRT(IONST) & /(1.0D0+BGA*PHRA(I)*DSQRT(IONST))+PHRB(I)*IONST END IF !---- equation 2 ! LOGI(I)=-AGA*ZZ(I)**2*DSQRT(IONST) ! >/(1.0D0+DSQRT(IONST)) !---- equation 3 !1 LOGI(I)=LOGI(I) !1 >-DLOG(1.0D0+0.0180153D0*IONSUM)/L10 !==== END IF XRES(I)=(GW(I)+ZZ(I)*MUZ)/RT+DLOG(XM(I))+LOGI(I)*L10 600 CONTINUE RETURN END !----- !******************************** SUBROUTINE AQPRT(MOLW,LOGI,IONST,IONSUM) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !-----Common block CHARACTER*32 CH16 REAL*8 X16,FF !-----END OF VARIABLES FOR MEGAMAT !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- REAL*8 LOGI(PHMAX),MOLW(PHMAX),GCORR(CALMAX), & IONST,IONSUM,ACT,LACT,EXPO,FF1,FF2, & LOGMI,GWCALC,PHW,ACC,LOGAI,XACC,PEW,AWAT INTEGER*4 I,II,I1,I2,ITOP(PHMAX) !---- WRITE (6,1050) WRITE (out,1050) 1050 FORMAT (/ & ' ----------------------------'/ & ' Composition of aqueous phase'/ & ' ----------------------------') WRITE (6,1003) MEGA,FERTIG WRITE (out,1003) MEGA,FERTIG 1003 FORMAT (/' number of iterations:',I4,' obvious errors =',I4) DO 300,I=1,FERTIG IF (UNF(I).GT.NCOL+1) THEN WRITE (UNIT=6,FMT='('' not converging: pot'')') WRITE (UNIT=out,FMT='('' not converging: pot'')') ELSE WRITE (UNIT=6,FMT='('' not converging: '',A16)') CNAME(UNF(I)) WRITE (UNIT=out,FMT='('' not converging: '',A16)') CNAME(UNF(I)) END IF 300 CONTINUE IF (PRTEST) WRITE (6,1005) ZSUM,MUZ,IONST,IONSUM IF (PRTEST) WRITE (out,1005) ZSUM,MUZ,IONST,IONSUM 1005 FORMAT (/ & ' total charge =',1PE12.5,10X, & 'chemical potential of charge =',0PF15.4/ & ' ionic strength =',1PE12.5,10X, & 'sum of molalities = ',1PE12.5) PHW=-DLOG(MOLW(ICH))/L10-LOGI(ICH) PEW=(-G(IEP)-MUZ)/RT/L10 IF (PRTEST) WRITE (6,1006) PHW,PEW IF (PRTEST) WRITE (out,1006) PHW,PEW 1006 FORMAT ( & ' pH =',F12.5,/ & ' pe =',F12.5) !----- CH16='pH' X16=PHW CALL SETTBL(CH16,X16) CH16='pe' X16=PEW CALL SETTBL(CH16,X16) !----- ITOP(1)=1 DO 400,I=2,NCOL DO 410,I1=1,I-1 IF (MOLW(ITOP(I1)).LT.MOLW(I)) THEN DO 420,I2=I,I1+1,-1 420 ITOP(I2)=ITOP(I2-1) ITOP(I1)=I GOTO 411 END IF 410 CONTINUE ITOP(I)=I 411 CONTINUE 400 CONTINUE !----- ! WRITE (6,1010) ! WRITE (out,1010) !1010 FORMAT (/ ! >' species',10X,'z(i)',10X,'G',14X,'G(calc)', ! >7X,'log g',9X,'log m',10X,'m(i)',10X,'a(i)'/ ! >' -------',10X,'----',10X,'-',14X,'-------', ! >7X,'-----',9X,'-----',10X,'----',10X,'----') ! DO 500,I=1,NCOL ! LOGMI=DLOG(MOLW(ITOP(I)))/L10 ! ACC=MOLW(ITOP(I))*DEXP(LOGI(ITOP(I))*L10) ! GWCALC=-RT*(DLOG(MOLW(ITOP(I)))+LOGI(ITOP(I))*L10)-ZZ(ITOP(I))*MUZ ! WRITE (6,1000) CNAME(ITOP(I)),ZZ(ITOP(I)),GW(ITOP(I)),GWCALC, ! >LOGI(ITOP(I)),LOGMI,MOLW(ITOP(I)),ACC ! WRITE (out,1000) CNAME(ITOP(I)),ZZ(ITOP(I)),GW(ITOP(I)),GWCALC, ! >LOGI(ITOP(I)),LOGMI,MOLW(ITOP(I)),ACC !1000 FORMAT (2X,A16,F4.1,2X,F15.4,2X,F15.4, ! >2X,1PE12.5,2X,1PE12.5,2X,1PE12.5,2X,1PE12.5) ! 500 CONTINUE !----- WRITE (6,1011) WRITE (out,1011) 1011 FORMAT (/ & ' species',10X,'z(i)',6X,'m(i)',9X,'log m', & 9X,'log g',10X,'a(i)',9X,'log a',8X,'x(i)', & 13X,'G',10X,'G(calc)'/ & ' -------',10X,'----',6X,'----',9X,'-----', & 9X,'-----',10X,'----',9X,'-----',8X,'----', & 13X,'-',10X,'-------') DO 501,I=1,NCOL LOGMI=DLOG(MOLW(ITOP(I)))/L10 ACC=MOLW(ITOP(I))*DEXP(LOGI(ITOP(I))*L10) GWCALC=-RT*(DLOG(MOLW(ITOP(I)))+LOGI(ITOP(I))*L10)-ZZ(ITOP(I)) & *MUZ LOGAI=DLOG(ACC)/L10 XACC=MOLW(ITOP(I))/(IONSUM+55.508435062D0) WRITE (6,1001) CNAME(ITOP(I)),ZZ(ITOP(I)),MOLW(ITOP(I)),LOGMI, & LOGI(ITOP(I)),ACC,LOGAI,XACC,GW(ITOP(I)),GWCALC WRITE (out,1001) CNAME(ITOP(I)),ZZ(ITOP(I)),MOLW(ITOP(I)),LOGMI, & LOGI(ITOP(I)),ACC,LOGAI,XACC,GW(ITOP(I)),GWCALC 1001 FORMAT (2X,A16,F4.1,2X,1PE12.5,2X,1PE12.5, & 2X,1PE12.5,2X,1PE12.5,2X,1PE12.5,2X,1PE10.3, & 1X,0PF13.2,1X,0PF13.2) !----- CH16='log_m.'//CNAME(ITOP(I))(1:10) X16=LOGMI CALL SETTBL(CH16,X16) !----- 501 CONTINUE !----- AWAT=1.0D0-0.017D0*IONSUM/55.508435062D0 IF (PRTEST) WRITE (6,1015) AWAT IF (PRTEST) WRITE (out,1015) AWAT 1015 FORMAT (/,'activity of water (Garrels & Christ, 1965): ',F10.6) AWAT=1.0D0-IONSUM/(IONSUM+55.508435062D0) IF (PRTEST) WRITE (6,1016) AWAT IF (PRTEST) WRITE (out,1016) AWAT 1016 FORMAT ('concentration of water ( 1 - sum(x) ) : ',F10.6) !----- !----- ends after 700 IF (PRTEST) THEN FF1=55.508435062D0 FF2=2*FF1 WRITE (6,1021) FF1,FF2 WRITE (out,1021) FF1,FF2 1021 FORMAT (/,' m(O in H2O)',5X,1PE12.5,/ & ' m(H in H2O)',5X,1PE12.5) WRITE (6,1020) WRITE (out,1020) 1020 FORMAT (/ & ' element',14X,'m(i)',10X,'g/Kg',24X,'G(orig)',7X,'dG(corr)'/ & ' -------',14X,'----',10X,'----',24X,'-------',7X,'--------') FF1=0.0D0 FF2=0.0D0 DO 600,II=1,NROW ! FF=0.0D0 ! DO 605,I=1,NCOL ! 605 FF=FF+AA(II,I)*MOLW(I) FF=0.0D0 DO 607,I=1,NMAX IF (NUMMER(I).EQ.II) FF=G(I) 607 CONTINUE FF1=FF1+QSUM(II) FF2=FF2+QSUM(II)*MOLWT(II) WRITE (6,1025) RNAME(II),QSUM(II),QSUM(II)*MOLWT(II), & CHNAME(II),FF,MUEG(II) WRITE (out,1025) RNAME(II),QSUM(II),QSUM(II)*MOLWT(II), & CHNAME(II),FF,MUEG(II) 1025 FORMAT (2X,A16,2(1PE12.5,2X),A16,1PE12.5,2X,1PE15.8) CH16='(aq)'//RNAME(II) X16=QSUM(II) CALL SETTBL(CH16,X16) 600 CONTINUE WRITE (6,1026) FF1,FF2 WRITE (out,1026) FF1,FF2 1026 FORMAT (18X,'------------',2X,'------------',/ & 18X,2(1PE12.5,2X)) ! !----- DO 750,I=1,NMAX GCORR(I)=G(I) DO 751,II=1,NUN 751 GCORR(I)=GCORR(I)+X(I,II)*MUEG(II) 750 CONTINUE WRITE (6,1030) WRITE (out,1030) 1030 FORMAT (/ & ' --------------------'/ & ' activities of phases'/ & ' --------------------'/ & /2X,'phase',19X,'N',13X,'G',9X,'Activity',7X,'Log(Act)', & /2X,'-----',19X,'-',13X,'-',9X,'--------',7X,'--------') DO 700,I=1,NMAX IF (NUMMER(I).GT.NUN) THEN IF ((PHASID(NUMMER(I)).NE.'AQU' & .AND.PHASID(NUMMER(I)).NE.'AQP')) THEN !DC EXPO=-G(I)/RT EXPO=-GCORR(I)/RT LACT=EXPO/2.302585093D0 IF (EXPO.LT.-150.0D0) THEN ACT=0.0D0 ELSE ACT=DEXP(DMIN1(150.0D0,EXPO)) END IF !DC WRITE (scr,1035) NAME(NUMMER(I)),NN(I),G(I),ACT,LACT !DC WRITE (out,1035) NAME(NUMMER(I)),NN(I),G(I),ACT,LACT WRITE (scr,1035) NAME(NUMMER(I)),NN(I),GCORR(I),ACT,LACT WRITE (out,1035) NAME(NUMMER(I)),NN(I),GCORR(I),ACT,LACT 1035 FORMAT (2X,A16,2X,1PE12.5,2X,1PE12.5,2X,1PE12.5,2X,1PE12.5) ! !----- CH16='log_a.'//NAME(NUMMER(I))(1:10) X16=LACT CALL SETTBL(CH16,X16) !----- END IF END IF 700 CONTINUE END IF RETURN END !----- !******************************** SUBROUTINE LINSOFT2(N0,ARR0,CRR0,DX0,ERRC) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !.....solves a linear system of n equations with n unknowns !..... !..... x(1)*A(1,1) + x(2)*A(1,2) + ... = C(1) !..... x(1)*A(2,1) + x(2)*A(2,2) + ... = C(2) !..... ... !..... x(1)*A(n,1) + x(2)*A(n,2) + ... = C(n) !..... !.....N0,N: number of equations !.....ARR0,ARR: coefficients !.....CRR0,CRR: constants (right hand side) !.....DX0: solution !.....ERRC: 0:OK, 1:rank not N0 REAL*8 ARR0(PHMAX,PHMAX),CRR0(PHMAX),DX0(PHMAX),COLSUM(PHMAX), & FF(COMAX),WERT,F1,FC,RTEST INTEGER*4 N0,IC,IR,ERRC,IIX,IIY,PIV,II,I !-----common blocks REAL*8 ARR(PHMAX,PHMAX),CRR(PHMAX) INTEGER*4 N,COLNR(PHMAX),RANG !-----common blocks !-----end of common RTEST=0.0D0 ERRC=0 N=N0 DO 300,IR=1,N COLNR(IR)=IR CRR(IR)=CRR0(IR) DO 300,IC=1,N ARR(IR,IC)=ARR0(IR,IC) 300 CONTINUE !----- !===== !===== reduc RANG=0 DO 600,PIV=1,N !----- mach colsum DO 400,IC=PIV,N COLSUM(IC)=0.0D0 DO 400,IR=1,N 400 COLSUM(IC)=COLSUM(IC)+DABS(ARR(IR,IC)) !----- end mach colsum IIY=0 IIX=0 WERT=0.0D0 !--- suche max DO 610,IR=PIV,N DO 620,IC=PIV,N IF (COLSUM(IC).LE.0.0D0) GOTO 609 F1=DABS(ARR(IR,IC)/COLSUM(IC)) IF (F1.GT.WERT) THEN WERT=F1 IIX=IR IIY=IC END IF 609 CONTINUE 620 CONTINUE 610 CONTINUE !--- end suche max IF (IIX.EQ.0.OR.IIY.EQ.0) GOTO 601 IF (IIX.NE.PIV) THEN !---- switch rows FC=CRR(IIX) DO 310,I=1,N 310 FF(I)=ARR(IIX,I) CRR(IIX)=CRR(PIV) DO 312,I=1,N 312 ARR(IIX,I)=ARR(PIV,I) CRR(PIV)=FC DO 314,I=1,N 314 ARR(PIV,I)=FF(I) !---- end switch rows END IF IF (IIY.NE.PIV) THEN !---- switch columns II=COLNR(IIY) DO 320,I=1,N 320 FF(I)=ARR(I,IIY) COLNR(IIY)=COLNR(PIV) DO 322,I=1,N 322 ARR(I,IIY)=ARR(I,PIV) COLNR(PIV)=II DO 324,I=1,N 324 ARR(I,PIV)=FF(I) !---- end switch columns END IF F1=ARR(PIV,PIV) RANG=PIV !---- lindiv CRR(PIV)=CRR(PIV)/F1 DO 330,I=1,N 330 ARR(PIV,I)=ARR(PIV,I)/F1 !---- end lindiv DO 500,II=1,N IF (II.NE.PIV) THEN F1=ARR(II,PIV) !---- linsub CRR(II)=CRR(II)-CRR(PIV)*F1 DO 340,I=1,N ARR(II,I)=ARR(II,I)-ARR(PIV,I)*F1 IF (RTEST.NE.0.0D0) THEN IF (DABS(ARR(II,I)).LE.RTEST) ARR(II,I)=0.0D0 END IF 340 CONTINUE !---- end linsub END IF 500 CONTINUE 600 CONTINUE !===== end reduc 601 CONTINUE !===== IF (RANG.NE.N) THEN WRITE (6,1000) RANG,N WRITE (out,1000) RANG,N 1000 FORMAT ('LINSOFT failed: rank =',I4,' N =',I4) ERRC=1 RETURN END IF !---- write result into DX0 DO 700,IC=1,N 700 DX0(COLNR(IC))=CRR(IC) !===== ! DO 800,IR=1,N ! WRITE (6,2000) IR,DX0(IR) ! WRITE (out,2000) IR,DX0(IR) !2000 FORMAT ('linsoft: X =',I5,F20.8) ! 800 CONTINUE !===== RETURN END !----- !******************************** SUBROUTINE LINSOFT(N0,AA0,CC0,XX,ERRC) include 'files.cmn' INTEGER*4 PHMAX PARAMETER (PHMAX=500) !.....solves a linear system of n equations with n unknowns !..... !..... x(1)*A(1,1) + x(2)*A(1,2) + ... = C(1) !..... x(1)*A(2,1) + x(2)*A(2,2) + ... = C(2) !..... ... !..... x(1)*A(n,1) + x(2)*A(n,2) + ... = C(n) !..... !.....N0,N: number of equations !.....AA0,AA: coefficients !.....CC0,CC: constants (right hand side) !.....XX: solution REAL*8 AA0(PHMAX,PHMAX),CC0(*),XX(*) INTEGER*4 N0,IC,IR,ERRC !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common ERRC=0 N=N0 DO 500,IR=1,N ROWNR(IR)=IR COLNR(IR)=IR CC(IR)=CC0(IR) DO 500,IC=1,N AA(IR,IC)=AA0(IR,IC) 500 CONTINUE !----- ! CALL LINPRT CALL REDUC IF (RANG.NE.N) THEN WRITE (6,1000) RANG,N WRITE (out,1000) RANG,N 1000 FORMAT ('LINSOFT failed: rank =',I4,' N =',I4) CALL LINPRT ERRC=1 RETURN END IF ! CALL LINPRT DO 550,IC=1,N 550 XX(COLNR(IC))=CC(IC) !===== ! DO 800,IR=1,N ! WRITE (6,2000) IR,XX(IR) ! WRITE (out,2000) IR,XX(IR) !2000 FORMAT ('linsoft: X =',I5,F20.8) ! 800 CONTINUE !===== RETURN END !----- !******************************** SUBROUTINE REDUC IMPLICIT NONE INTEGER*4 PHMAX PARAMETER (PHMAX=500) !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common INTEGER*4 PIV REAL*8 COLSUM(PHMAX),WERT !---- INTEGER*4 IR,IC,IX,IY,II,I2 REAL*8 FF,F1 RANG=0 DO 600,PIV=1,N DO 400,IC=PIV,N COLSUM(IC)=0.0D0 DO 400,IR=1,N 400 COLSUM(IC)=COLSUM(IC)+DABS(AA(IR,IC)) IY=0 IX=0 WERT=0.0D0 DO 610,IR=PIV,N DO 620,IC=PIV,N IF (COLSUM(IC).LE.0.0D0) GOTO 609 FF=DABS(AA(IR,IC)/COLSUM(IC)) IF (FF.GT.WERT) THEN WERT=FF IX=IR IY=IC END IF 609 CONTINUE 620 CONTINUE ! IF (WERT.NE.0.0D0) GOTO 611 610 CONTINUE 611 CONTINUE IF (IX.EQ.0.OR.IY.EQ.0) RETURN IF (IX.NE.PIV) CALL LINROW(IX,PIV) IF (IY.NE.PIV) CALL LINCOL(IY,PIV) FF=AA(PIV,PIV) RANG=PIV CALL LINDIV(PIV,FF) DO 500,II=1,N IF (II.NE.PIV) THEN F1=AA(II,PIV) I2=II CALL LINSUB(PIV,F1,I2) END IF 500 CONTINUE 600 CONTINUE RETURN END !----- !******************************** SUBROUTINE LINPRT IMPLICIT NONE include 'files.cmn' INTEGER*4 PHMAX PARAMETER (PHMAX=500) !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common INTEGER*4 I,II !---- WRITE (UNIT=6,FMT='(10(20X,10I10/))') (COLNR(I),I=1,N) WRITE (UNIT=out,FMT='(10(20X,10I10/))') (COLNR(I),I=1,N) DO 500,II=1,N WRITE (UNIT=6,FMT='(/I10,11E10.3,10(/10X,10E10.3))') & ROWNR(II),CC(II),(AA(II,I),I=1,N) WRITE (UNIT=out,FMT='(/I10,11E10.3,10(/10X,10E10.3))') & ROWNR(II),CC(II),(AA(II,I),I=1,N) 500 CONTINUE RETURN END !----- !******************************** SUBROUTINE LINDIV(I1,FF) IMPLICIT NONE INTEGER*4 PHMAX PARAMETER (PHMAX=500) !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common INTEGER*4 I1,I REAL*8 FF CC(I1)=CC(I1)/FF DO 500,I=1,N 500 AA(I1,I)=AA(I1,I)/FF RETURN END !---- !----- !******************************** SUBROUTINE LINCOL(I1,I2) IMPLICIT NONE INTEGER*4 PHMAX PARAMETER (PHMAX=500) !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common INTEGER*4 I,II,I1,I2 REAL*8 FF(PHMAX) !---- II=COLNR(I1) DO 500,I=1,N 500 FF(I)=AA(I,I1) COLNR(I1)=COLNR(I2) DO 510,I=1,N 510 AA(I,I1)=AA(I,I2) COLNR(I2)=II DO 520,I=1,N 520 AA(I,I2)=FF(I) RETURN END !----- !******************************** SUBROUTINE LINROW(I1,I2) IMPLICIT NONE INTEGER*4 PHMAX PARAMETER (PHMAX=500) !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common INTEGER*4 I,II,I1,I2 REAL*8 FF(PHMAX),FC !---- II=ROWNR(I1) FC=CC(I1) DO 500,I=1,N 500 FF(I)=AA(I1,I) ROWNR(I1)=ROWNR(I2) CC(I1)=CC(I2) DO 510,I=1,N 510 AA(I1,I)=AA(I2,I) ROWNR(I2)=II CC(I2)=FC DO 520,I=1,N 520 AA(I2,I)=FF(I) RETURN END !----- !******************************** SUBROUTINE LINSUB(I1,F1,I2) IMPLICIT NONE INTEGER*4 PHMAX PARAMETER (PHMAX=500) !-----common blocks REAL*8 AA(PHMAX,PHMAX),CC(PHMAX) INTEGER*4 N,COLNR(PHMAX),ROWNR(PHMAX),RANG COMMON /LINRE/ AA,CC COMMON /LININ/ N,COLNR,ROWNR,RANG !-----end of common INTEGER*4 I,I1,I2 REAL*8 F1 !---- CC(I2)=CC(I2)-CC(I1)*F1 DO 500,I=1,N AA(I2,I)=AA(I2,I)-AA(I1,I)*F1 ! IF (DABS(AA(I2,I)).LE.1.0D-10) AA(I2,I)=0.0D0 500 CONTINUE RETURN END !----- !******************************** SUBROUTINE ABPHREEQ(T,A,B) IMPLICIT NONE REAL*8 T,A,B,tc,tk,s1,s2,s3,c1 ! ! Lines stolen from PHREEQ file: model.c ! ! compute temperature dependence of a and b for debye-huckel ! tc=T-273.15D0 tk=T s1=374.11-tc !C WRITE (UNIT=6,FMT='('' s1 ='',1PE15.8)') s1 ! s2=pow(s1,1.0/3.0); ! s2=s1**(1.0/3.0) !DC IF (s1.LT.0.0D0) THEN s2=-(-s1)**(1.0/3.0) ELSE s2=s1**(1.0/3.0) END IF !C WRITE (UNIT=6,FMT='('' s2 ='',1PE15.8)') s2 s3=1.0+0.1342489*s2-3.946263e-03*s1 !C WRITE (UNIT=6,FMT='('' s3 ='',1PE15.8)') s3 s3=s3/(3.1975-0.3151548*s2-1.203374e-03*s1+ & 7.48908e-13*(s1*s1*s1*s1)) !CCC WRITE (UNIT=6,FMT='('' s3a='',1PE15.8)') s3 s3=DSQRT(s3) !CCC WRITE (UNIT=6,FMT='('' s3b='',1PE15.8)') s3 if (tk.GT.373.15) THEN c1=5321.0/tk+233.76-tk*(tk*(8.292e-07*tk-1.417e-03)+0.9297) !C WRITE (UNIT=6,FMT='('' c1a='',1PE15.8)') c1 ELSE ! /* replaced by wateq4f expression ! c1=87.74-tc_x*(tc_x*(1.41e-06*tc_x-9.398e-04)+0.4008); ! */ c1 = 2727.586+0.6224107*tk-466.9151*DLOG(tk)-52000.87/tk !C WRITE (UNIT=6,FMT='('' c1b='',1PE15.8)') c1 END IF c1=sqrt(c1*tk) !CCC WRITE (UNIT=6,FMT='('' c1c='',1PE15.8)') c1 ! /* replaced by wateq4f expressions ! a=1824600.0*s3/(c1 * c1 * c1); ! b=50.29*s3/c1; ! */ A = 1824827.7 * s3 / (c1 * c1 * c1) !C WRITE (UNIT=6,FMT='('' A ='',1PE15.8)') A B = 50.2905 * s3 / c1 !C WRITE (UNIT=6,FMT='('' B ='',1PE15.8)') B RETURN END !----- !******************************** !--- SUBROUTINE POTPHAS(AQPHAS) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !----- INTEGER*4 I,IN,AQPHAS(COMAX) !----- DO 500,IN=1,NUN !===== NPHA=NPHA+1 AQPHAS(IN)=NPHA SUGNR=SUGNR+1 NAME(NPHA)='aq-pot.'//CHNAME(IN) ABK(NPHA)='a-p.'//CHNAME(IN) PHASID(NPHA)='POT' DO 600,I=1,NUN 600 XX(NPHA,I)=0.0D0 XX(NPHA,IN)=1.0D0 EMSOL(NPHA)=0 EMNR(NPHA)=0 NULL(NPHA)=.TRUE. G0R=0.0D0 H0R=0.0D0 S0R=0.0D0 V0R=0.0D0 AA0=0.0D0 AAT=0.0D0 BB0=0.0D0 BBT=0.0D0 K1=0.0D0 K2=0.0D0 K3=0.0D0 K4=0.0D0 K5=0.0D0 K6=0.0D0 K7=0.0D0 K8=0.0D0 K9=0.0D0 NLANDA=0 VOLVO=0 TRTYP=0 NLANDA=0 NCOM=0 DO 523,I=1,10 ICOM(I)=0 523 FFCOM(I)=0.0D0 DO 517,I=1,4 ASPK(I)=0.0D0 BSPK(I)=0.0D0 TQ1B(I)=0.0D0 TEQ(I)=0.0D0 DVDT(I)=0.0D0 DVDP(I)=0.0D0 TRE(I)=0.0D0 DHTR(I)=0.0D0 517 DVTR(I)=0.0D0 TD0=0.0D0 TDMAX=0.0D0 VADJ=0.0D0 D1=0.0D0 D2=0.0D0 D3=0.0D0 D4=0.0D0 D5=0.0D0 D6=0.0D0 D7=0.0D0 D8=0.0D0 D9=0.0D0 VTA=0.0D0 VTB=0.0D0 VPA=0.0D0 VPB=0.0D0 VAA=0.0D0 VAB=0.0D0 VB=0.0D0 VL0=0.0D0 VLA=0.0D0 VLN=0.0D0 VL2=0.0D0 TKRI=0.0D0 SMA=0.0D0 NAT=0.0D0 FALL=' ' RDK=.FALSE. VDW=.FALSE. SPC=.FALSE. COM=.TRUE. DIS=.FALSE. VO1=.TRUE. VO2=.FALSE. VO3=.FALSE. AQU=.FALSE. FIX=.FALSE. TL1=.FALSE. CALL DASAVE(NPHA) !===== 500 CONTINUE !----- RETURN END !----- !******************************** !--- SUBROUTINE ADDAQ(XES,NDLOOP) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !----- INTEGER*4 I,IN,AQPHAS(COMAX),NDLOOP REAL*8 XES(COMAX),THEG CHARACTER*8 TEXT !----- !===== WRITE (UNIT=TEXT,FMT='(''fluid'',I3.3)') NDLOOP IF (NAME(NPHA)(1:5).NE.'fluid') NPHA=NPHA+1 SUGNR=SUGNR+1 ! NAME(NPHA)='fluid' NAME(NPHA)=TEXT ABK(NPHA)='fl' PHASID(NPHA)='POT' DO 600,I=1,NUN 600 XX(NPHA,I)=XES(I) EMSOL(NPHA)=0 EMNR(NPHA)=0 NULL(NPHA)=.FALSE. G0R=0.0D0 H0R=0.0D0 S0R=0.0D0 V0R=0.0D0 AA0=0.0D0 AAT=0.0D0 BB0=0.0D0 BBT=0.0D0 K1=0.0D0 K2=0.0D0 K3=0.0D0 K4=0.0D0 K5=0.0D0 K6=0.0D0 K7=0.0D0 K8=0.0D0 K9=0.0D0 NLANDA=0 VOLVO=0 TRTYP=0 NLANDA=0 NCOM=0 DO 523,I=1,10 ICOM(I)=0 523 FFCOM(I)=0.0D0 DO 517,I=1,4 ASPK(I)=0.0D0 BSPK(I)=0.0D0 TQ1B(I)=0.0D0 TEQ(I)=0.0D0 DVDT(I)=0.0D0 DVDP(I)=0.0D0 TRE(I)=0.0D0 DHTR(I)=0.0D0 517 DVTR(I)=0.0D0 TD0=0.0D0 TDMAX=0.0D0 VADJ=0.0D0 D1=0.0D0 D2=0.0D0 D3=0.0D0 D4=0.0D0 D5=0.0D0 D6=0.0D0 D7=0.0D0 D8=0.0D0 D9=0.0D0 VTA=0.0D0 VTB=0.0D0 VPA=0.0D0 VPB=0.0D0 VAA=0.0D0 VAB=0.0D0 VB=0.0D0 VL0=0.0D0 VLA=0.0D0 VLN=0.0D0 VL2=0.0D0 TKRI=0.0D0 SMA=0.0D0 FALL=' ' RDK=.FALSE. VDW=.FALSE. SPC=.FALSE. COM=.TRUE. DIS=.FALSE. VO1=.TRUE. VO2=.FALSE. VO3=.FALSE. AQU=.FALSE. FIX=.FALSE. TL1=.FALSE. CALL DASAVE(NPHA) !===== !----- RETURN END !******************************** !--- SUBROUTINE ADDWATER(INH2O,IXH2O) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !----- INTEGER*4 I,IN,AQPHAS(COMAX),NDLOOP,INH2O,IXH2O REAL*8 XES(COMAX),THEG,MUWAT CHARACTER*8 TEXT !----- !===== IXH2O=0 DO 500,I=1,NPHA IF (NAME(I).EQ.'H2Opot') IXH2O=I 500 CONTINUE IF (IXH2O.NE.0) RETURN NPHA=NPHA+1 SUGNR=SUGNR+1 NAME(NPHA)='H2Opot' NAME(NPHA)=TEXT ABK(NPHA)='H2Op' PHASID(NPHA)='POT' DO 600,I=1,NUN 600 XX(NPHA,I)=XX(INH2O,I) EMSOL(NPHA)=0 EMNR(NPHA)=0 NULL(NPHA)=.FALSE. G0R=0.0D0 H0R=0.0D0 S0R=0.0D0 V0R=0.0D0 AA0=0.0D0 AAT=0.0D0 BB0=0.0D0 BBT=0.0D0 K1=0.0D0 K2=0.0D0 K3=0.0D0 K4=0.0D0 K5=0.0D0 K6=0.0D0 K7=0.0D0 K8=0.0D0 K9=0.0D0 NLANDA=0 VOLVO=0 TRTYP=0 NLANDA=0 NCOM=0 DO 523,I=1,10 ICOM(I)=0 523 FFCOM(I)=0.0D0 DO 517,I=1,4 ASPK(I)=0.0D0 BSPK(I)=0.0D0 TQ1B(I)=0.0D0 TEQ(I)=0.0D0 DVDT(I)=0.0D0 DVDP(I)=0.0D0 TRE(I)=0.0D0 DHTR(I)=0.0D0 517 DVTR(I)=0.0D0 TD0=0.0D0 TDMAX=0.0D0 VADJ=0.0D0 D1=0.0D0 D2=0.0D0 D3=0.0D0 D4=0.0D0 D5=0.0D0 D6=0.0D0 D7=0.0D0 D8=0.0D0 D9=0.0D0 VTA=0.0D0 VTB=0.0D0 VPA=0.0D0 VPB=0.0D0 VAA=0.0D0 VAB=0.0D0 VB=0.0D0 VL0=0.0D0 VLA=0.0D0 VLN=0.0D0 VL2=0.0D0 TKRI=0.0D0 SMA=0.0D0 FALL=' ' RDK=.FALSE. VDW=.FALSE. SPC=.FALSE. COM=.TRUE. DIS=.FALSE. VO1=.TRUE. VO2=.FALSE. VO3=.FALSE. AQU=.FALSE. FIX=.FALSE. TL1=.FALSE. CALL DASAVE(NPHA) !===== !----- RETURN END !----- !******************************** SUBROUTINE PRTSHORT IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !----- WRITE (UNIT=scr,FMT=1002) WRITE (UNIT=out,FMT=1002) 1002 FORMAT (/ & ' ----------------'/ & ' stable assemblge'/ & ' ----------------'/ & 9X,'phase',19X,'N',13X,'G',9X,'activity') CALL PRTSTR(1,NUN) !----- RETURN END !----- !******************************** SUBROUTINE ANASS(CURAS) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !----- CHARACTER*125 CURAS INTEGER*4 I,I1,IL CHARACTER*16 TEXT !----- !+++++ CURAS=' ' I1=1 DO 500,I=1,NUN2 IF (NUMMER(I).LE.0) THEN TEXT=SOLNAM(EMCODE(I)) ELSE TEXT=NAME(NUMMER(I)) END IF CALL LABLA(TEXT,IL) CURAS(I1:)=TEXT I1=I1+IL+1 500 CONTINUE !+++++ !----- RETURN END !----- !******************************** !--- SUBROUTINE SPEZSUB(AQFAIL) IMPLICIT NONE INCLUDE 'theriak.cmn' include 'files.cmn' !---- !-----Common block REAL*8 AA(COMAX,PHMAX),GW(PHMAX),ZZ(PHMAX),AGA,BGA,L10, & PHRA(PHMAX),PHRB(PHMAX),GV(PHMAX),CPOT(COMAX),AQPOT(COMAX), & QSUM(COMAX),MUEG(COMAX),MUZ,ZSUM INTEGER*4 NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF(0:500),IEP, & NPOT,POT(COMAX) CHARACTER*16 CNAME(PHMAX),RNAME(COMAX) LOGICAL*4 PRTEST COMMON /MARE/ AA,GW,ZZ,AGA,BGA,L10,PHRA,PHRB,GV,CPOT,AQPOT,QSUM, & MUEG,MUZ,ZSUM COMMON /MAIN/ NCOL,NROW,IZ,ICH,MEGA,FERTIG,UNF,IEP,NPOT,POT COMMON /MACH/ CNAME,RNAME COMMON /MALO/ PRTEST !---- REAL*8 XKGH2O,XXMIN,XXMAX,XXVAL,XXDELTA,IONSUM INTEGER*4 I,AQFAIL,NBEO,NAQLOOP CHARACTER*16 XXVAR !---- INTEGER*4 NRLOOP,DHPAR REAL*8 NRFIDI,NRMUZ,NRRED COMMON /NERAIN/ NRLOOP,DHPAR COMMON /NERARE/ NRFIDI,NRMUZ,NRRED !----- XXVAR='Ox-Buffer' XXMIN=-35.0D0 XXMAX=-45.0D0 NAQLOOP=50+1 XXDELTA=(XXMAX-XXMIN)/DBLE(NAQLOOP) XXVAL=XXMIN-XXDELTA !----- DO 4500,NBEO=1,NAQLOOP WRITE (UNIT=6,FMT='(/''spezial loop:'',I4)') NBEO WRITE (UNIT=out,FMT='(/''spezial loop:'',I4)') NBEO CALL DBREAD !----------------------------------------------------------------- XXVAL=XXVAL+XXDELTA DO 400,I=1,NPHA IF (NAME(I).EQ.XXVAR) THEN REDATA(3,I)=-R*XXVAL*2.302585093D0 END IF 400 CONTINUE !----------------------------------------------------------------- DO 672,I=1,NPHA IF (PHASID(I).EQ.'AQU'.OR.PHASID(I).EQ.'AQP') THEN NULL(I)=.TRUE. END IF 672 CONTINUE CALL NURVONPT CALL CALSTR IF (PRTLOG(2).OR.PRTLOG(3).OR.PRTLOG(4)) CALL PRININ CALL THERIA !----------------------------------------------------------------- !----------------------------------------------------------------- !---- CALL AQEQ(AQFAIL,IONSUM) ! DO 750,I=1,NC ! SOLBUL(I)=0.0D0 ! LABUL(I)=BULK(I) ! 750 CHE(I)=0.0D0 ! ! ! DO 754,I=1,NUN2 ! ISASOLID(I)=.TRUE. ! IF (NUMMER(I).LE.0) THEN ! TEXT(I)=SOLNAM(EMCODE(I)) ! IS=EMCODE(I) ! DO 755,II=1,NEND(IS) ! IF (PHASID(EM(IS,II)).NE.'MIN') ISASOLID(I)=.FALSE. ! 755 CONTINUE ! ELSE ! TEXT(I)=NAME(NUMMER(I)) ! IF (PHASID(NUMMER(I)).NE.'MIN') ISASOLID(I)=.FALSE. ! END IF ! 754 CONTINUE !C----- ! DO 757,I=1,NUN ! IF (OXYDE(CHMCOD(I)).NE.'O '.AND.OXYDE(CHMCOD(I)).NE.'H '. ! >AND.OXYDE(CHMCOD(I)).NE.'E ') THEN ! CHE(CHMCOD(I))=BULK(I)-XKGH2O*55.0D+1*QSUM(I) ! ELSE ! CHE(CHMCOD(I))=BULK(I) ! END IF ! IF (CHE(CHMCOD(I)).LE.0.0D0) CHE(CHMCOD(I))=0.0D0 ! 757 CONTINUE ! MORE=.FALSE. ! DO 760,I=1,NC ! CHE(I)=CHE(I)+XKGH2O*55.0D+1*INFLOW(I) ! IF (CHE(I).EQ.0.0D0.NEQV.CHEM(I).EQ.0.0D0) MORE=.TRUE. ! CHEM(I)=CHE(I) ! 760 CONTINUE ! IF (MORE) THEN ! CALL DBREAD ! ELSE ! DO 765,I=1,NUN ! 765 BULK(I)=CHE(CHMCOD(I)) ! END IF 4500 CONTINUE !---- the last call to DBREAD is apparently not necessary (??) !----------------------------------------------------------------- !----------------------------------------------------------------- RETURN END