subroutine ph_eqadvan(npt,k,co2i,co3i,hco3i,hsi,h2si,boh3i, & boh4i,initflag,debug,hp,co2,co3,hco3,hs,h2s,boh3,boh4) C------------------------------------------------------------------- C PH_EQADVAN $Revision: 1.0 $ $Date: 1998/12/07 $ C Copyright (C) GEOMAR, Roger Luff, Matthias Haeckel 1998. C C USAGE: C INTEGER npt,initflag,debug C double precision k(12) C double precision co2i(npt),co3i(npt),hco3i(npt),hsi(npt) C double precision h2si(npt),boh3i(npt),boh4i(npt) C double precision co2(npt),co3(npt),hco3(npt),hs(npt),h2s(npt) C double precision boh3(npt),boh4(npt),hp(npt) C call ph_eqadvan(npt,k,co2i,co3i,hco3i,hsi,h2si,boh3i, C & boh4i,initflag,debug,hp,co2,co3,hco3,hs,h2s,boh3,boh4) C SEE ph_test.f FOR DETAILS OF USING THIS PROGRAM!!!!!!!!!!!!!!! C C DESCRIPTION: C This function calculates the H+ concentration and the concentration C of the dissociated species to the closest equilibrium of co2,co3,hco3,hs, C h2s,boh3,boh4 from the dissociated species not in equilibrium co2,co3, C hco3,hs,h2s,boh3,boh4, and the equilibrium constants. To calculate C the pH from H+ use: ph=-log10(hp/1.0E+03) C C INPUT: C k = equilibrium constants from ph_thermo C hpi = Initial concentration H+ [mmol/l] C co2i = Initial concentration of CO2 [mmol/l] C co3i = Initial concentration of CO3[mmol/l] C hco3i = Initial concentration of HCO3 [mmol/l] C hsi = Initial concentration of HS [mmol/l] C h2si = Initial concentration of H2S [mmol/l] C boh3i = Initial concentration of BOH3 [mmol/l] C boh4i = Initial concentration of BOH4 [mmol/l] C initflag = Initialization flag 1=initialize carbonate species C debug = debug flag 1=debug on (will slow solver) C C OUTPUT: C hp = Equilibrium concentration H+ [mmol/l] C co2 = Equilibrium concentration of CO2 [mmol/l] C co3 = Equilibrium concentration of CO3[mmol/l] C hco3 = Equilibrium concentration of HCO3 [mmol/l] C hs = Equilibrium concentration of HS [mmol/l] C h2s = Equilibrium concentration of H2S [mmol/l] C boh3 = Equilibrium concentration of BOH3 [mmol/l] C boh4 = Equilibrium concentration of BOH4 [mmol/l] C C AUTHOR: C Roger Luff and Matthias Haeckel 98-12-07 (rluff@gmx.de) C GEOMAR Research Centre for Marine Geosciences C Wischhofstr. 1-3, D-24148 Kiel, Germany C visumod@technologist.com C C DISCLAIMER: C This software is provided "as it is" without warranty of any kind, C see ph_info for details. C C REFERENCES: C Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN C and MATLAB libraries to calculate pH distributions in C marine systems", Computers & Geosciences. C------------------------------------------------------------------- implicit none C------------------------------------------------------------------- C External variables C------------------------------------------------------------------- INTEGER npt,initflag,debug double precision k(12) double precision co2i(npt),co3i(npt),hco3i(npt),hsi(npt) double precision h2si(npt),boh3i(npt),boh4i(npt) double precision co2(npt),co3(npt),hco3(npt),hs(npt),h2s(npt) double precision boh3(npt),boh4(npt),hp(npt) C------------------------------------------------------------------- C Internal variables C----------------------------------------------------------------- integer jakobi parameter (jakobi=12) integer i, j, kk, info double precision tol INTEGER IPVT(jakobi),iter,maxiter INTEGER zeile,spalte double precision MATRIX(jakobi,jakobi),eqwork(jakobi),det(2) double precision f(jakobi),hpi,delta1,delta2,delta3,delta4 double precision dd(jakobi) double precision atoll C------------------------------------------------------------------- C INITIALIZE THE EQUILIBRIUM VARIABLES WITH THE INITIAL VALUES C------------------------------------------------------------------- if (initflag.eq.1) then print *, 'init ph_eqadvan' DO I=1,npt co2(i) = co2i(i) hco3(i) = hco3i(i) co3(i) = co3i(i) if (co2i(i).le.0.1) then hp(i) = hco3i(i)*K(3)/co3i(i) else hp(i) = co2i(i)*K(2)/hco3i(i) endif H2S(i) = h2si(i) HS(i) = hsi(i) BOH3(i) = boh3i(i) BOH4(i) = boh4i(i) enddo endif C------------------------------------------------------------------- C SET UP LOCAL VARIABLES C------------------------------------------------------------------- maxiter = 150 tol = 1.0D-11 DO I=1,npt delta1 = 0.0 delta2 = 0.0 delta3 = 0.0 delta4 = 0.0 iter = 0 hpi = hp(i) C------------------------------------------------------------------- C ITERATION LOOP STARTS HERE FOR CO3 AND HCO3 AS DEPENDENT VARIABLES C------------------------------------------------------------------- 1000 continue iter = iter + 1 do spalte = 1,jakobi do zeile = 1,jakobi MATRIX(zeile,spalte)=0.0 enddo enddo C------------------------------------------------------------------- C SET UP THE JACOBI-MATRIX C------------------------------------------------------------------- crl line 1 MATRIX(1,1) = -1.0 MATRIX(1,7) = -k(2) crl line 2 MATRIX(2,2) = -1.0 MATRIX(2,7) = hp(i) MATRIX(2,8) = -k(3) crl line 3 MATRIX(3,3) = -1.0 MATRIX(3,8) = hp(i) crl line 4 MATRIX(4,4) = -1.0 MATRIX(4,9) = -k(10) crl line 5 MATRIX(5,5) = -1.0 MATRIX(5,9) = hp(i) crl line 6 MATRIX(6,6) = -1.0 MATRIX(6,7) = hco3(i) MATRIX(6,8) = co3(i) MATRIX(6,9) = hs(i) MATRIX(6,12) = boh4(i) crl line 7 MATRIX(7,1) = -1.0 MATRIX(7,2) = 1.0 MATRIX(7,6) = 1.0 crl line 8 MATRIX(8,2) = -1.0 MATRIX(8,3) = 1.0 MATRIX(8,6) = 1.0 crl line 9 MATRIX(9,4) = -1.0 MATRIX(9,5) = 1.0 MATRIX(9,6) = 1.0 crl line 10 MATRIX(10,6) = 1.0 MATRIX(10,10) = -1.0 MATRIX(10,11) = 1.0 crl line 11 MATRIX(11,10) = -1.0 MATRIX(11,12) = -k(4) crl line 12 MATRIX(12,11) = -1.0 MATRIX(12,12) = hp(i) C------------------------------------------------------------------- C INVERT THE JACOBI USING LINPAC LIBRARY C------------------------------------------------------------------- call DGEFA(MATRIX,jakobi,jakobi,ipvt,info) if (info.eq.0) then call dgedi(MATRIX,jakobi,jakobi,ipvt,det,eqwork,01) else print *,'ERROR: Jakobian contains a 0 on the diagonal' print *,'pH-eqalk solver stopped calculation' stop endif C------------------------------------------------------------------- C SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED C------------------------------------------------------------------- f(1) = co2i(i)-delta1-co2(i) f(2) = hco3i(i)+delta1-delta2-hco3(i) f(3) = co3i(i)+delta2-co3(i) f(4) = h2si(i)-delta3-h2s(i) f(5) = hsi(i)+delta3-hs(i) f(6) = hpi+delta1+delta2+delta3+delta4-hp(i) f(7) = hco3(i)*hp(i)-k(2)*co2(i) f(8) = co3(i)*hp(i)-k(3)*hco3(i) f(9) = hs(i)*hp(i)-k(10)*h2s(i) f(10)= boh3i(i)-delta4-boh3(i) f(11)= boh4i(i)+delta4-boh4(i) f(12)= boh4(i)*hp(i)-k(4)*boh3(i) do j=1,jakobi dd(j) = 0.0 enddo C------------------------------------------------------------------- C CALCULATE THE DISTANCE TO THE NEXT INTERATION STEP C------------------------------------------------------------------- do j=1,jakobi do kk=1,jakobi dd(kk) =dd(kk) +MATRIX(j,kk)*f(j) enddo enddo C------------------------------------------------------------------- C CALCULATE THE NEW DISTRIBUTION C------------------------------------------------------------------- co2(i)=co2(i)-dd(1) hco3(i)=hco3(i)-dd(2) co3(i)=co3(i)-dd(3) h2s(i)=h2s(i)-dd(4) hs(i)=hs(i)-dd(5) hp(i) =hp(i)-dd(6) delta1=delta1-dd(7) delta2=delta2-dd(8) delta3=delta3-dd(9) delta4=delta4-dd(10) boh3(i)=boh3(i)-dd(11) boh4(i)=boh4(i)-dd(12) atoll = 0.0 C------------------------------------------------------------------- C CHECK THE CONVERGENCE C------------------------------------------------------------------- do j=1,jakobi atoll = atoll + abs(dd(j)) enddo if((atoll.gt.tol).and.(iter.lt.maxiter)) goto 1000 C------------------------------------------------------------------- C IN THE CASE THAT THE FUNCTION DOES NOT CONVERGENTE, END WITH AND C ERROR MESSAGE C------------------------------------------------------------------- if (iter.ge.maxiter) then print *,'Iter = Maxiter' print *, 'pH_eqalk: too much iterations' stop endif C------------------------------------------------------------------- C DEBUG PART, PROGRAM WILL STOP IN CASE OF ERROR. THIS CHECK C SLOWS THE SOLVER SIGNIFICANTLY C------------------------------------------------------------------- if (debug.eq.1) then if (co2(i).lt.0.0) then print *, 'Error C<0 at co2, layer',co2(i),i stop 'In EQUILIB1' endif if (co3(i).lt.0.0) then print *, 'Error C<0 at co3, layer',co3(i),i stop 'In EQUILIB1' endif if (hco3(i).lt.0.0) then print *, 'Error C<0 at hco3, layer',hco3(i),i stop 'In EQUILIB1' endif if (boh3(i).lt.0.0) then print *, 'Error C<0 at BOH3, layer',boh3(i),i stop 'In EQUILIB1' endif if (boh4(i).lt.0.0) then print *, 'Error C<0 at BOH4, layer',boh4(i),i stop 'In EQUILIB1' endif if (h2s(i).lt.0.0) then print *, 'Error C<0 at h2s, layer',h2s(i),i stop 'In EQUILIB1' endif if (hs(i).lt.0.0) then print *, 'Error C<0 at hs, layer',hs(i),i stop 'In EQUILIB1' endif if (hp(i).lt.0.0) then print *, 'Error C<0 at hp, layer',hp(i),i stop 'In EQUILIB1' endif endif enddo initflag = 0 return end