program ph_test C------------------------------------------------------------------- C DESCRIPTION: C Test program for the ph-library C - Calculate the density of sewater C - Calculate the total borate concentration from seawater C - Calculate the thermodynamic constants C - Calculate the species distribtions with ph_eqalk C - Calculate the mole fractions C - Calculate the species distribtions with ph_eqcharge C - Print the results to screen C C AUTHOR: C Roger Luff 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. implicit none INTEGER npt parameter (npt=10) INTEGER initflag,debug C------------------------------------------------------------------- c used by all solvers C------------------------------------------------------------------- double precision k(12) double precision tbohw,so4w,ssum,caw double precision RSW,S,T,p double precision co2(npt),co3(npt),hco3(npt),hs(npt),h2s(npt) double precision boh3(npt),boh4(npt),hp(npt) C------------------------------------------------------------------- c used by alkalinity conservation and charge conservation solver C------------------------------------------------------------------- double precision alk(npt),tco2(npt),th2s(npt),tboh(npt) C------------------------------------------------------------------- c used by charge conservation solver C------------------------------------------------------------------- double precision tnh4(npt),tpo4(npt),mnii(npt),feii(npt),ca(npt) double precision so4(npt),no3(npt) double precision h3po4(npt),h2po4(npt),hpo4(npt),po4(npt) double precision nh3(npt),nh4(npt) C------------------------------------------------------------------- c used by advancement solver C------------------------------------------------------------------- double precision co2i(npt),co3i(npt),hco3i(npt),hsi(npt) double precision h2si(npt),boh3i(npt),boh4i(npt) C------------------------------------------------------------------- c mole fractions with respect to the sum concentrations C------------------------------------------------------------------- double precision fco2(npt),fco3(npt),fhco3(npt) double precision fh2s(npt),fhs(npt),fboh3(npt),fboh4(npt) double precision fahco3(npt),faco3(npt),faboh4(npt),fahs(npt) C------------------------------------------------------------------- c local variables C------------------------------------------------------------------- integer i double precision ph C------------------------------------------------------------------- c Parameters of the environment C------------------------------------------------------------------- s = 35.0 t = 4.0 p = 250.0 C------------------------------------------------------------------- C Calculate the density of sewater C------------------------------------------------------------------- call ph_density(RSW,S,T,P) C------------------------------------------------------------------- C Calculate the total borate concentration from seawater C------------------------------------------------------------------- call ph_swconst(s,rsw,tbohw,so4w,ssum,caw) C------------------------------------------------------------------- C Calculate the thermodynamic constants C------------------------------------------------------------------- call ph_tdconst(k,t,s,p) C------------------------------------------------------------------- C Initialize parameters for ph_eqalk C------------------------------------------------------------------- print *, 'The concentration of TH2S increases with depth!' do i=1,npt tco2(i) = 2.2 th2s(i) = float(i)/10.0 alk(i) = 2.5 tboh(i) = tbohw enddo initflag = 1 debug = 0 C------------------------------------------------------------------- C Calculate the species distribtions with ph_eqalk C------------------------------------------------------------------- call ph_eqalk(npt,k,alk,tco2,th2s,tboh,initflag,debug,hp, & co2,co3,hco3,hs,h2s,boh3,boh4) C------------------------------------------------------------------- C Calculate the mole fractions C------------------------------------------------------------------- call ph_eqfac (npt,1,k,hp,alk,hco3,hs,boh4,fco2,fco3,fhco3, & fh2s,fhs,fboh3,fboh4,fahco3,faco3,faboh4,fahs) C------------------------------------------------------------------- C Print the results to screen C------------------------------------------------------------------- print *, 'Results from ph_eqalk:' do i=1,npt ph = -LOG10(HP(i)/1.0D+03) print *, 'pH in layer', i, ' = ',ph,' co2fac = ',fco2(i) enddo C------------------------------------------------------------------- C Initialize parameters for ph_eqcharge C------------------------------------------------------------------- do i=1,npt tco2(i) = 2.2 th2s(i) = float(i)/10.0 tboh(i) = tbohw so4(i) = so4w ca(i) = caw tnh4(i) = 0.0 tpo4(i) = 0.0 no3(i) = 0.0 mnii(i) = 0.0 feii(i) = 0.0 enddo initflag = 1 debug = 0 C------------------------------------------------------------------- C Calculate the species distribtions with ph_eqcharge C------------------------------------------------------------------- call 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------------------------------------------------------------------- C Calculate the mole fractions C------------------------------------------------------------------- call ph_eqfac (npt,3,k,hp,alk,hco3,hs,boh4,fco2,fco3,fhco3, & fh2s,fhs,fboh3,fboh4,fahco3,faco3,faboh4,fahs) C------------------------------------------------------------------- C Print the results to screen C------------------------------------------------------------------- print *, 'Results from ph_eqcharge:' do i=1,npt ph = -LOG10(HP(i)/1.0D+03) print *, 'pH in layer', i, ' = ',ph,' co2fac = ',fco2(i) enddo C------------------------------------------------------------------- C Initialize parameters for ph_eqadvan C------------------------------------------------------------------- do i=1,npt co2i(i) = tco2(i)*0.1+float(i)/10.0 hco3i(i) = tco2(i)*0.8+float(i)/10.0 co3i(i) = tco2(i)*0.1+float(i)/10.0 h2si(i) = th2s(i)*0.5 hsi(i) = th2s(i)*0.5 boh3i(i) = boh3(i) boh4i(i) = boh4(i) enddo initflag = 1 debug = 0 C------------------------------------------------------------------- C Calculate the species distribtions with ph_eqadvan C------------------------------------------------------------------- call ph_eqadvan(npt,k,co2i,co3i,hco3i,hsi,h2si,boh3i, & boh4i,initflag,debug,hp,co2,co3,hco3,hs,h2s,boh3,boh4) C------------------------------------------------------------------- C Calculate the mole fractions C------------------------------------------------------------------- call ph_eqfac (npt,2,k,hp,alk,hco3,hs,boh4,fco2,fco3,fhco3, & fh2s,fhs,fboh3,fboh4,fahco3,faco3,faboh4,fahs) C------------------------------------------------------------------- C Print the results to screen C------------------------------------------------------------------- print *, 'Results from ph_eqadvan:' do i=1,npt ph = -LOG10(HP(i)/1.0D+03) print *, 'pH in layer', i, ' = ',ph,' co2fac = ',fco2(i) enddo print *, 'Note, resultes above represent different scenarios!' end