subroutine ph_eqcharge(npt,k,ssum,tco2,th2s,tboh,tnh4,tpo4,mnii, & feii,ca,so4,no3,initflag,debug,hp,co2,co3,hco3, & hs,h2s,h3po4,h2po4,hpo4,po4,nh3,nh4,boh3,boh4) C DESCRIPTION: C This function calculates the H+ concentration and the concentration C of the dissociated species from the INPUT variables, listed later C on the bases of the charge balance equation. C To calculate the pH from H+ use: ph=-log10(hp/1.0E+03) C C USAGE c INTEGER npt c parameter (npt=10) c double precision k(12) c double precision ssum c double precision co2(npt),co3(npt),hco3(npt),hs(npt),h2s(npt) c double precision boh3(npt),boh4(npt),hp(npt) c double precision tnh4(npt),tpo4(npt),mnii(npt),feii(npt),ca(npt) c double precision so4(npt),no3(npt) c double precision h3po4(npt),h2po4(npt),hpo4(npt),po4(npt) c double precision nh3(npt),nh4(npt) c call ph_eqcharge(npt,k,ssum,tco2,th2s,tboh,tnh4,tpo4,mnii, C & feii,ca,so4,no3,initflag,debug,hp,co2,co3,hco3, C & hs,h2s,h3po4,h2po4,hpo4,po4,nh3,nh4,boh3,boh4) C SEE ph_test.f FOR DETAILS OF USING THIS PROGRAM!!!!!!!!!!!!!!! C C INPUT: C k = equilibrium constants from ph_thermo C ssum = sum concentration of the conservative species in seawater C from ph_swconst C tco2 = Total carbonate concentration [mmol/l] C th2s = Total suffate concentration [mmol/l] C tboh = Total borate concentration [mmol/l] C tnh4 = Total amonium concentration [mmol/l] C tpo4 = Total phospate concentration [mmol/l] C mnii = Mn2+ concentration [mmol/l] C feii = Fe2+ concentration [mmol/l] C ca = concentration of total calcite in seawater C so4 = concentration of total sulfate in seawater C no3 = Nitrate concentration [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 h3po4 = Equilibrium concentration of H3PO4 [mmol/l] C h2po4 = Equilibrium concentration of H2PO4 [mmol/l] C hpo4 = Equilibrium concentration of HPO4 [mmol/l] C nh3 = Equilibrium concentration of NH3 [mmol/l] C nh4 = Equilibrium concentration of NH4 [mmol/l] C po4 = Equilibrium concentration of PO4 [mmol/l] C boh3 = Equilibrium concentration of BOH3 [mmol/l] C boh4 = Equilibrium concentration of BOH4 [mmol/l] C C AUTHOR: C Bernard Boudreau(1), Roger Luff(2) and Matthias Haeckel(2) C (1) Dept. Oceanography C Dalhousie University C Halifax, NS B3H 4J1,Canada C bernie@boudreau.ocean.dal.ca C C (2) 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 C REFERENCES: C Boudreau B. P. (1996) A method-of-lines code for carbon and C nutrient diagenesis in aquatic sediments.Computers & Geosciences C 22(5), 479-496. 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 tco2(npt),th2s(npt),tboh(npt),tnh4(npt) double precision tpo4(npt),mnii(npt),feii(npt),ca(npt) double precision so4(npt),no3(npt) double precision co2(npt),co3(npt),hco3(npt),hs(npt),h2s(npt) double precision boh3(npt),boh4(npt),hp(npt),h3po4(npt) double precision h2po4(npt),hpo4(npt),po4(npt),nh3(npt),nh4(npt) C------------------------------------------------------------------- C INTERNAL VARIABLES C------------------------------------------------------------------- integer jakobi parameter (jakobi=15) integer i, j, kk, info double precision tol INTEGER IPVT(jakobi),zeile,spalte,iter,maxiter double precision MATRIX(jakobi,jakobi), eqwork(jakobi),det(2) double precision f(jakobi) double precision dd(jakobi) double precision atoll,ssum,oh C------------------------------------------------------------------- C SET A FIRST GUESS TO THE FINAL EQUILIBRIUM VALUES BASED ON THE C SUM CONCENTRATIONS C------------------------------------------------------------------- if (initflag.eq.1) then print *, 'init ph_eqcharge' DO I=1,npt co2(i) = tco2(i)*0.1 hco3(i) = tco2(i)*0.8 co3(i) = tco2(i)*0.1 h2s(i) = th2s(i)*0.5 hs(i) = th2s(i)*0.5 h3po4(i)= tpo4(i)*0.01 h2po4(i)= tpo4(i)*0.45 hpo4(i) = tpo4(i)*0.45 po4(i) = tpo4(i)*0.09 nh3(i) = tnh4(i)*0.5 nh4(i) = tnh4(i)*0.5 boh3(i) = tboh(i)*0.5 boh4(i) = tboh(i)*0.5 hp(i) = co2(i)*k(2)/hco3(i) oh = k(1)/hp(i) enddo endif C------------------------------------------------------------------- C SET UP LOCAL VARIABLES C------------------------------------------------------------------- maxiter = 50 tol = 1.0D-12 DO I=1,npt iter = 0 C------------------------------------------------------------------- C ITERATION LOOP STARTS HERE C------------------------------------------------------------------- 2000 continue iter = iter + 1 C------------------------------------------------------------------- C SET UP THE JACOBI C------------------------------------------------------------------- do spalte = 1,jakobi do zeile = 1,jakobi MATRIX(zeile,spalte)=0.0 enddo enddo crl line 1 MATRIX(1,1) = -1.0 MATRIX(1,7) = k(2) crl line 2 MATRIX(2,1) = -1.0 MATRIX(2,7) = -hp(i) MATRIX(2,8) = k(3) MATRIX(2,15) = -1.0 crl line 3 MATRIX(3,1) = -1.0 MATRIX(3,8) = -hp(i) MATRIX(3,15) = -2.0 crl line 4 MATRIX(4,2) = -1.0 MATRIX(4,14) = k(10) crl line 5 MATRIX(5,2) = -1.0 MATRIX(5,14) = -hp(i) MATRIX(5,15) = -1.0 crl line 6 MATRIX(6,3) = -1.0 MATRIX(6,10) = k(5) crl line 7 MATRIX(7,3) = -1.0 MATRIX(7,10) = -hp(i) MATRIX(7,11) = k(6) MATRIX(7,15) = -1.0 crl line 8 MATRIX(8,3) = -1.0 MATRIX(8,11) = -hp(i) MATRIX(8,12) = k(7) MATRIX(8,15) = -2.0 crl line 9 MATRIX(9,3) = -1.0 MATRIX(9,12) = -hp(i) MATRIX(9,15) = -3.0 crl line 10 MATRIX(10,4) = -1.0 MATRIX(10,13) = -hp(i) crl line 11 MATRIX(11,4) = -1.0 MATRIX(11,13) = k(9) MATRIX(11,15) = 1.0 crl line 12 MATRIX(12,5) = -1.0 MATRIX(12,9) = k(4) crl line 13 MATRIX(13,5) = -1.0 MATRIX(13,9) = -hp(i) MATRIX(13,15) = -1.0 crl line 14 MATRIX(14,6) = -oh MATRIX(14,7) = -hco3(i) MATRIX(14,8) = -co3(i) MATRIX(14,9) = -boh4(i) MATRIX(14,10) = -h2po4(i) MATRIX(14,11) = -hpo4(i) MATRIX(14,12) = -po4(i) MATRIX(14,13) = -nh3(i) MATRIX(14,14) = -hs(i) MATRIX(14,15) = 1.0 crl line 15 MATRIX(15,6) = -hp(i) MATRIX(15,15) = -1.0 C------------------------------------------------------------------- C INVERT THE JACOBI USING THE LINPACK 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) = tco2(i) - (co2(i)+co3(i)+hco3(i)) F(2) = th2s(i) - (h2s(i)+hs(i)) F(3) = tpo4(i) - (h3po4(i) + h2po4(i) + hpo4(i) + po4(i)) F(4) = tnh4(i) - (nh3(i) + nh4(i)) F(5) = tboh(i) - (boh3(i) + boh4(i)) F(6) = K(1) - hp(i)*oh F(7) = K(2)*co2(i) - hp(i)*hco3(i) f(8) = K(3)*hco3(i) - hp(i)*co3(i) f(9) = k(4)*boh3(i) - hp(i)*boh4(i) f(10)= k(5)*h3po4(i) - hp(i)*h2po4(i) f(11) = k(6)*h2po4(i) - hp(i)*hpo4(i) f(12) = k(7)*hpo4(i) - hp(i)*po4(i) f(13) = k(9)* nh4(i) - hp(i)*nh3(i) f(14) = k(10)*h2s(i) - hp(i)*hs(i) f(15) = SSUM+2.0*ca(i)+hp(i)+nh4(i)+2.0*(MnII(i)+FeII(i))- & hco3(i)-h2po4(i)-hs(i)-boh4(i)-2.0*so4(i)- & 2.0*co3(i)-2.0*hpo4(i)-3.0*po4(i)-oh-no3(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) h3po4(i)=h3po4(i)-dd(6) h2po4(i)=h2po4(i)-dd(7) hpo4(i)=hpo4(i)-dd(8) po4(i)=po4(i)-dd(9) nh3(i)=nh3(i)-dd(10) nh4(i)=nh4(i)-dd(11) boh3(i)=boh3(i)-dd(12) boh4(i)=boh4(i)-dd(13) hp(i) =hp(i)-dd(14) oh =oh-dd(15) 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 2000 C------------------------------------------------------------------- C IN THE CASE THAT THE FUNCTION DOES NOT CONVERGENTE, END WITH AND C ERROR MESSAGE OR WITH A H+ CONCENTRATION OF NaN C------------------------------------------------------------------- if (iter.eq.maxiter) then print *,'Iter = Maxiter' print *, 'pH_eqcharge: 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).le.0.0) stop 'Error CO2 < 0' if (hco3(i).le.0.0) stop 'Error HCO3 < 0' if (co3(i).le.0.0) stop 'Error co3 < 0' if (h2s(i).le.0.0) stop 'Error H2S < 0' if (hs(i).le.0.0) stop 'Error HS < 0' if (h3po4(i).le.0.0) stop 'Error H3PO4 < 0' if (h2po4(i).le.0.0) stop 'Error H2PO4 < 0' if (hpo4(i).le.0.0) stop 'Error HPO4 < 0' if (po4(i).le.0.0) stop 'Error PO4 < 0' if (nh4(i).le.0.0) stop 'Error NH4 < 0' if (nh3(i).le.0.0) stop 'Error NH3 < 0' if (boh3(i).le.0.0) stop 'Error BOH3 < 0' if (boh4(i).le.0.0) stop 'Error BOH4 < 0' if (hp(i).le.0.0) stop 'Error HP < 0' if (oh.le.0.0) stop 'Error OH < 0' endif enddo initflag = 0 return end