subroutine pH_eqalk(npt,k,alk,tco2,th2s,tboh,initflag,debug,hp, & co2,co3,hco3,hs,h2s,boh3,boh4) C------------------------------------------------------------------- C PH_EQALK $Revision: 1.1 $ $Date: 2001/01/07 $ C Copyright (C) GEOMAR, Roger Luff 1998. C C USAGE: C parameter (npt=100) C INTEGER npt,initflag,debug C double precision k(12) C double precision alk(npt),tco2(npt),th2s(npt),tboh(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_eqalk(npt,k,alk,tco2,th2s,tboh,initflag,debug,hp, C & co2,co3,hco3,hs,h2s,boh3,boh4) C SEE ph_test.f FOR DETAILS OF USING THIS PROGRAM!!!!!!!!!!!!!!! C C DESCRIPTION: C This subroutine calculates the H+ concentration and the concentration C of the dissociated species co2,co3,hco3,hs,h2s,boh3,boh4 from C alkalinity, total carbonate, total sulfide and total borate C concentration and the equilibrium constants. To calculate the pH from H+ C use: ph(i)=-log10(hp(i)/1.0E+03) C C INPUT: C npt = number of layers in the model C k = equilibrium constants from ph_thermo C alk = Alkalinity [meq/l] C tco2 = Total carbonate concentration [mmol/l] C th2s = Total suffate concentration [mmol/l] C tboh = Total borate 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 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@geomar.de) C GEOMAR Research Centre for Marine Geosciences C Wischhofstr. 1-3, D-24148 Kiel, Germany C rluff@geomar.de, mhaeckel@geomar.de C C DISCLAIMER: C This software is provided "as it is" without warranty of any kind. C C REFERENCES: C Luff R., Haeckel M. and Wallmann K. (2001) "Robust and fast FORTRAN C and MATLAB libraries to calculate pH distributions in marine systems", C Computers & Geosciences(27), 157-169. C------------------------------------------------------------------- implicit none integer jakobi parameter (jakobi=2) C------------------------------------------------------------------- C External variables C------------------------------------------------------------------- INTEGER npt,initflag,debug double precision k(12) double precision alk(npt),tco2(npt),th2s(npt),tboh(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------------------------------------------------------------------- double precision k3k2,helptb,helpts,r,ran1 double precision co3i,co2i integer i, j, info,randini double precision tol,atoll, relaxfac,alktco2,th2stco2,rfacc INTEGER IPVT(jakobi),iter,maxiter double precision MATRIX(jakobi,jakobi) double precision DETT,MTMP double precision f(jakobi),ddco2,ddco3,ddhco3 C------------------------------------------------------------------- C INITIALIZE THE CARBONAT SPECIES C------------------------------------------------------------------- if (initflag.eq.1) then print *, 'init ph_eqalk' DO I=1,npt co3(i) = 0.09*tco2(i) co2(i) = 0.11*tco2(i) hco3(i) = 0.8*tco2(i) enddo randini = -1 r = ran1(randini) endif c relaxfac = Relaxationsfactor: improvement of the linear c convergence of the N-R algorithm, analogously to the method c of successive overrelaxation (SOR). C------------------------------------------------------------------- C SET UP LOCAL VARIABLES C------------------------------------------------------------------- maxiter = 1000 tol = 1.0D-12 k3k2 = k(3)/k(2) do i=1,npt iter = 0 rfacc = 0.0 alktco2=alk(i)/tco2(i) th2stco2=th2s(i)/tco2(i) C------------------------------------------------------------------- C alktco2 > 1 : CALCULATE EQUILIBRIUM USING CO3 AND HCO3 AS UNKNOWN C alktco2 < 1 : CALCULATE EQUILIBRIUM USING CO2 AND HCO3 AS UNKNOWN C------------------------------------------------------------------- if (alktco2 .gt.1.0) then C------------------------------------------------------------------- C CHECK IF FUNCTION CAN CALCULATE A SOLUTION, SEE REFERENCE C------------------------------------------------------------------- if (alktco2 .lt. 2.0+th2stco2) then C------------------------------------------------------------------- C ITERATION LOOP STARTS HERE FOR CO3 AND HCO3 AS DEPENDENT VARIABLES C------------------------------------------------------------------- 1000 continue rfacc = rfacc + 0.1 relaxfac = min(1.0,0.3+rfacc) iter = iter + 1 co2i = k3k2*hco3(i)**2/co3(i) helptb = k(3)/k(4)*hco3(i)/co3(i) helpts = k(3)/k(10)*hco3(i)/co3(i) C------------------------------------------------------------------- C SET UP THE JACOBI-MATRIX C------------------------------------------------------------------- crl Line 1 MATRIX(1,1) = 1.0+2.0*k3k2*hco3(i)/co3(i) MATRIX(1,2) = 1.0-tboh(i)/(1.0+helptb)**2*helptb/HCO3(i) & -th2s(i)/(1.0+helpts)**2*helpts/HCO3(i) crl Line 2 MATRIX(2,1) = 1.0-k3k2*hco3(i)**2/co3(i)**2 MATRIX(2,2) = 2.0+tboh(i)/(1.0+helptb)**2*helptb/CO3(i) & +th2s(i)/(1.0+helpts)**2*helpts/CO3(i) C------------------------------------------------------------------- C INVERT THE JACOBI C First calculate the Determinant C------------------------------------------------------------------- Dett = MATRIX(1,1)*MATRIX(2,2)-(MATRIX(1,2)*MATRIX(2,1)) C------------------------------------------------------------------- crl secondly we need 1/DET C------------------------------------------------------------------- Dett = 1.0/Dett C------------------------------------------------------------------- crl third calculate the inverse now C------------------------------------------------------------------- MTMP = Dett*MATRIX(2,2) MATRIX(1,2) = -Dett*MATRIX(1,2) MATRIX(2,1) = -Dett*MATRIX(2,1) MATRIX(2,2) = Dett*MATRIX(1,1) MATRIX(1,1) = MTMP C------------------------------------------------------------------- C SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED C------------------------------------------------------------------- f(1)=(co2i+hco3(i)+co3(i))-tco2(i) f(2)=(hco3(i)+2.0*co3(i)+tboh(i)/(1.0+helptb) & +th2s(i)/(1.0+helpts))-alk(i) C------------------------------------------------------------------- C CALCULATE THE DISTANCE TO THE NEXT INTERATION STEP C------------------------------------------------------------------- ddhco3 = MATRIX(1,1)*f(1) + MATRIX(2,1)*f(2) ddco3 = MATRIX(1,2)*f(1) + MATRIX(2,2)*f(2) C------------------------------------------------------------------- C CALCULATE THE NEW DISTRIBUTION C------------------------------------------------------------------- hco3(i)=hco3(i)-ddhco3*relaxfac co3(i)=co3(i)-ddco3*relaxfac C------------------------------------------------------------------- C if Newton-Raphson is going in the wrong direction try it C again with a new choice of first guess concentration. C------------------------------------------------------------------- if (hco3(i) .le. 0.0 .or. co3(i) .le. 0.0) then r = ran1(randini) co3(i)=tco2(i)*r hco3(i)=tco2(i)*(1.0-r) ddhco3=2.0 rfacc = 0.0 endif C------------------------------------------------------------------- C CHECK THE CONVERGENCE C------------------------------------------------------------------- atoll = abs(ddhco3)+abs(ddco3) if ((atoll.gt.tol).and.(iter.lt.maxiter)) goto 1000 C------------------------------------------------------------------- C CALCULATE THE DIRECT DEPENDING VARIABLES C------------------------------------------------------------------- co2(i) = k3k2*hco3(i)**2/co3(i) helptb = k(3)/k(4)*hco3(i)/co3(i) helpts = k(3)/k(10)*hco3(i)/co3(i) boh4(i) = tboh(i)/(1.0+helptb) hs(i) = th2s(i)/(1.0+helpts) hp(i) = k(3)*hco3(i)/co3(i) else C------------------------------------------------------------------- C PRINT A WARNING WHEN NO SOLUTION CAN BE CALCULATED C------------------------------------------------------------------- print *, 'No calculation because of bad input values' endif else C------------------------------------------------------------------- C ITERATION LOOP STARTS HERE FOR CO2 AND HCO3 AS DEPENDENT VARIABLES C------------------------------------------------------------------- 2000 continue rfacc = rfacc + 0.1 relaxfac = min(1.0,0.3+rfacc) iter = iter + 1 co3i = k3k2*hco3(i)**2/co2(i) helptb = k(2)/k(4)*co2(i)/hco3(i) helpts = k(2)/k(10)*co2(i)/hco3(i) C------------------------------------------------------------------- C SET UP THE JACOBI C------------------------------------------------------------------- crl Line 1 MATRIX(1,1) = 1.0+2.0*k3k2*hco3(i)/co2(i) MATRIX(1,2) = 1.0+4.0*k3k2*hco3(i)/co2(i) & +tboh(i)/(1.0+helptb)**2*helptb/hco3(i) & +th2s(i)/(1.0+helpts)**2*helpts/hco3(i) crl Line 2 MATRIX(2,1) = 1.0-k3k2*hco3(i)**2/co2(i)**2 MATRIX(2,2) = -2.0*k3k2*hco3(i)**2/co2(i)**2 & -tboh(i)/(1.0*helptb)**2*helptb/co2(i) & -th2s(i)/(1.0*helpts)**2*helpts/co2(i) C------------------------------------------------------------------- C INVERT THE JACOBI C First calculate the Determinant C------------------------------------------------------------------- Dett = MATRIX(1,1)*MATRIX(2,2)-(MATRIX(1,2)*MATRIX(2,1)) C------------------------------------------------------------------- crl secondly we need 1/DET C------------------------------------------------------------------- Dett = 1.0/Dett C------------------------------------------------------------------- crl third calculate the inverse now C------------------------------------------------------------------- MTMP = Dett*MATRIX(2,2) MATRIX(1,2) = -Dett*MATRIX(1,2) MATRIX(2,1) = -Dett*MATRIX(2,1) MATRIX(2,2) = Dett*MATRIX(1,1) MATRIX(1,1) = MTMP C------------------------------------------------------------------- C SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED C------------------------------------------------------------------- f(1)=(co2(i)+hco3(i)+co3i)-tco2(i) f(2)=(hco3(i)+2.0*co3i+tboh(i)/(1.0+helptb) & +th2s(i)/(1.0+helpts))-alk(i) C------------------------------------------------------------------- C CALCULATE THE DISTANCE TO THE NEXT INTERATION STEP C------------------------------------------------------------------- ddhco3 = MATRIX(1,1)*f(1) + MATRIX(2,1)*f(2) ddco2 = MATRIX(1,2)*f(1) + MATRIX(2,2)*f(2) C------------------------------------------------------------------- C CALCULATE THE NEW DISTRIBUTION C------------------------------------------------------------------- hco3(i)=hco3(i)-ddhco3*relaxfac co2(i)=co2(i)-ddco2*relaxfac C------------------------------------------------------------------- C CHECK THE CONVERGENCE C------------------------------------------------------------------- atoll = abs(ddhco3)+abs(ddco2) if ((atoll.gt.tol).and.(iter.lt.maxiter)) goto 2000 C------------------------------------------------------------------- C CALCULATE THE DIRECT DEPENDING VARIABLES C------------------------------------------------------------------- co3(i) = k3k2*hco3(i)**2/co2(i) helptb = k(2)/k(4)*co2(i)/hco3(i) helpts = k(2)/k(10)*co2(i)/hco3(i) boh4(i) = tboh(i)/(1.0+helptb) hs(i) = th2s(i)/(1.0+helpts) hp(i) = k(2)/hco3(i)*co2(i) endif boh3(i) = tboh(i)-boh4(i) h2s(i) = th2s(i)-hs(i) crl ph(i)=-LOG10(HP(i)/1.0D+03) 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 (hp(i).lt.0.0) then print *, 'Error C<0 at HP, layer', hp(i), i stop 'In Alkalinity' endif if (co2(i).lt.0.0) then print *, 'Error C<0 at co2, layer', co2(i), i stop 'In Alkalinity' endif if (co3(i).lt.0.0) then print *, 'Error C<0 at co3, layer', co3(i), i stop 'In Alkalinity' endif if (hco3(i).lt.0.0) then print *, 'Error C<0 at hco3, layer', hco3(i), i stop 'In Alkalinity' endif if (h2s(i).lt.0.0) then print *, 'Error C<0 at H2S, layer', h2s(i), i print *, co2(i), hco3(i), co3(i), th2s(i) stop 'In Alkalinity' endif if (hs(i).lt.0.0) then print *, 'Error C<0 at HS, layer', hs(i), i stop 'In Alkalinity' endif if (boh4(i).lt.0.0) then print *, 'Error C<0 at HS, layer', hs(i), i stop 'In Alkalinity' endif if (boh3(i).lt.0.0) then print *, 'Error C<0 at HS, layer', hs(i), i stop 'In Alkalinity' endif endif enddo crl call eqfac initflag = 0 return end