!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
! Richard D. Slater
!
!
! John P. Dunne
!
!
!
! Surface fCO2 calculation
!
!
!
! Calculate the fugacity of CO2 at the surface in thermodynamic
! equilibrium with the current alkalinity (Alk) and total dissolved
! inorganic carbon (DIC) at a particular temperature and salinity
! using an initial guess for the total hydrogen
! ion concentration (htotal)
!
!
module ocmip2_co2calc_mod !{
!
!------------------------------------------------------------------
!
! modules
!
!------------------------------------------------------------------
!
!
!------------------------------------------------------------------
!
! Global definitions
!
!------------------------------------------------------------------
!
implicit none
private
public :: ocmip2_co2calc
public :: ocmip2_co2_alpha
character(len=128) :: version = '$Id: ocmip2_co2calc.F90,v 16.0.70.1 2009/10/14 14:37:19 smg Exp $'
character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'
!
!-----------------------------------------------------------------------
!
! Subroutine and function definitions
!
!-----------------------------------------------------------------------
!
contains
!#######################################################################
!
!
!
! Calculate CO2 solubility, alpha, from
! temperature (t) and salinity (s).
!
! INPUT
!
! isd = first i-limit of the arrays with halo
! ied = last i-limit of the arrays with halo
! jsd = first j-limit of the arrays with halo
! jed = last j-limit of the arrays with halo
! isc = first i-limit of the arrays for computation
! iec = last i-limit of the arrays for computation
! jsc = first j-limit of the arrays for computation
! jec = last j-limit of the arrays for computation
!
! t = temperature (degrees C)
!
! s = salinity (PSU)
!
! mask = land mask array (0.0 = land)
!
! OUTPUT
! alpha = Solubility of CO2 for air? (mol/kg/atm unless scaled)
!
! IMPORTANT: Some words about units - (JCO, 4/4/1999)
!
! - Models may carry tracers in mol/m^3 (on a per volume basis)
!
! - Conversely, this routine, which was written by observationalists
! (C. Sabine and R. Key), passes input arguments in umol/kg
! (i.e., on a per mass basis)
!
! - Thus, if the input or output units need to be changed from/to mol/m^3
! then set scale to an appropriate value. For example, if the model
! uses mol/m^3, then scale should be set to something like 1.0/1024.5
! to convert from mol/m^3 to mol/kg.
!
!
subroutine ocmip2_co2_alpha(isd, ied, jsd, jed, isc, iec, jsc, jec, t, s, mask, alpha, scale) !{
implicit none
!
! arguments
!
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: isc
integer, intent(in) :: iec
integer, intent(in) :: jsc
integer, intent(in) :: jec
real, dimension(isd:,jsd:), intent(in) :: t
real, dimension(isd:,jsd:), intent(in) :: s
real, dimension(isd:,jsd:), intent(in) :: mask
real, dimension(isc:,jsc:), intent(out) :: alpha
real, intent(in), optional :: scale
!
! local variables
!
integer :: i
integer :: j
integer :: k
real :: log100
real :: tk
real :: tk100
real :: tk1002
real :: logtk
real :: ff
real :: scale_factor
!
! set the scale factor for unit conversion
!
if (present(scale)) then !{
scale_factor = scale
else !}{
scale_factor = 1.0
endif !}
!
!---------------------------------------------------------------------
!
!***********************************************************************
! Calculate all constants needed to convert between various measured
! carbon species. References for each equation are noted in the code.
! Once calculated, the constants are stored and passed in the common
! block "const". The original version of this code was based on
! the code by Dickson in Version 2 of "Handbook of Methods for the
! Analysis of the Various Parameters of the Carbon Dioxide System
! in Seawater", DOE, 1994 (SOP No. 3, p25-26).
!
! Derive simple terms used more than once
!
log100 = log(100.0)
do j = jsc, jec !{
do i = isc, iec !{
if (mask(i,j) .ne. 0.0) then !{
tk = 273.15 + t(i,j)
tk100 = tk / 100.0
tk1002 = tk100 * tk100
logtk = log(tk)
!
! f = k0(1-pH2O)*correction term for non-ideality
!
! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6
! values)
!
ff = exp(-162.8301 + 218.2968 / tk100 + &
90.9241 * (logtk - log100) - 1.47696 * tk1002 + &
s(i,j) * (0.025695 - 0.025225 * tk100 + 0.0049867 * tk1002))
!
!---------------------------------------------------------------
!c Should we be using K0 or ff for the solubility here?
! convert to output units
!---------------------------------------------------------------
!
alpha(i,j) = ff / scale_factor
else !}{
alpha(i,j) = 0.0
endif !}
enddo !} i
enddo !} j
return
end subroutine ocmip2_co2_alpha !}
! NAME="ocmip2_co2_alpha"
!#######################################################################
!
!
!
! Calculate co2* from total alkalinity and total CO2 at
! temperature (t) and salinity (s).
! It is assumed that init_ocmip2_co2calc has already been called with
! the T and S to calculate the various coefficients.
!
! INPUT
!
! isd = first i-limit of the arrays with halo
! ied = last i-limit of the arrays with halo
! jsd = first j-limit of the arrays with halo
! jed = last j-limit of the arrays with halo
! isc = first i-limit of the arrays for computation
! iec = last i-limit of the arrays for computation
! jsc = first j-limit of the arrays for computation
! jec = last j-limit of the arrays for computation
!
! mask = land mask array (0.0 = land)
!
! t = temperature (degrees C)
!
! s = salinity (PSU)
!
! dic_in = total inorganic carbon (mol/kg unless scaled)
! where 1 T = 1 metric ton = 1000 kg
!
! ta_in = total alkalinity (eq/kg unless scaled)
!
! pt_in = inorganic phosphate (mol/kg unless scaled)
!
! sit_in = inorganic silicate (mol/kg unless scaled)
!
! htotallo = factor to set lower limit of htotal range
!
! htotalhi = factor to set upper limit of htotal range
!
! htotal = H+ concentraion
!
! OUTPUT
! co2star = CO2*water (kg/kg unless scaled)
! alpha = Solubility of CO2 for air? (kg/kg/atm unless scaled)
! pco2surf = oceanic pCO2 (ppmv)
!
! k1 = activity factors for carbonate species
!
! k2 (see below)
!
! invtk = 1/(t+273.15)
!
! IMPORTANT: Some words about units - (JCO, 4/4/1999)
!
! - Models may carry tracers in mol/m^3 (on a per volume basis)
!
! - Conversely, this routine, which was written by observationalists
! (C. Sabine and R. Key), passes input arguments in umol/kg
! (i.e., on a per mass basis)
!
! - Thus, if the input or output units need to be changed from/to mol/m^3
! then set scale to an appropriate value. For example, if the model
! uses mol/m^3, then scale should be set to something like 1.0/1024.5
! to convert from mol/m^3 to mol/kg.
!
!
subroutine ocmip2_co2calc(isd, ied, jsd, jed, isc, iec, jsc, jec, mask, t, s, &
dic_in, ta_in, pt_in, sit_in, htotallo, htotalhi, htotal, &
co2star, co3_ion, alpha, pCO2surf, k1_out, k2_out, invtk_out, scale) !{
implicit none
!
! local parameters
!
real, parameter :: permeg = 1.e-6
real, parameter :: xacc = 1.0e-10
!
! arguments
!
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: isc
integer, intent(in) :: iec
integer, intent(in) :: jsc
integer, intent(in) :: jec
real, dimension(isd:,jsd:), intent(in) :: mask
real, dimension(isd:,jsd:), intent(in) :: t
real, dimension(isd:,jsd:), intent(in) :: s
real, dimension(isd:,jsd:), intent(in) :: dic_in
real, dimension(isd:,jsd:), intent(in) :: ta_in
real, dimension(isd:,jsd:), intent(in) :: pt_in
real, dimension(isd:,jsd:), intent(in) :: sit_in
real, dimension(isc:,jsc:), intent(in) :: htotallo
real, dimension(isc:,jsc:), intent(in) :: htotalhi
real, dimension(isc:,jsc:), intent(inout) :: htotal
real, dimension(isc:,jsc:), intent(out), optional :: co2star
real, dimension(isc:,jsc:), intent(out), optional :: co3_ion
real, dimension(isc:,jsc:), intent(out), optional :: alpha
real, dimension(isc:,jsc:), intent(out), optional :: pCO2surf
real, dimension(isc:,jsc:), intent(out), optional :: invtk_out
real, dimension(isc:,jsc:), intent(out), optional :: k1_out
real, dimension(isc:,jsc:), intent(out), optional :: k2_out
real, intent(in), optional :: scale
!
! local variables
!
integer :: i
integer :: j
real :: log100
real :: pt
real :: sit
real :: ta
real :: dic
real :: tk
real :: invtk
real :: tk100
real :: tk1002
real :: logtk
real :: is
real :: is2
real :: sqrtis
real :: s2
real :: sqrts
real :: s15
real :: scl
real :: logf_of_s
real :: ff
!real :: k0
real :: k1
real :: k2
real :: kb
real :: k1p
real :: k2p
real :: k3p
real :: ksi
real :: kw
real :: ks
real :: kf
real :: bt
real :: st
real :: ft
real :: htotal2
real :: co2star_internal
real :: scale_factor
!
! set the scale factor for unit conversion
!
if (present(scale)) then !{
scale_factor = scale
else !}{
scale_factor = 1.0
endif !}
!
!---------------------------------------------------------------------
!
!***********************************************************************
! Calculate all constants needed to convert between various measured
! carbon species. References for each equation are noted in the code.
! Once calculated, the constants are stored and passed in the common
! block "const". The original version of this code was based on
! the code by Dickson in Version 2 of "Handbook of Methods for the
! Analysis of the Various Parameters of the Carbon Dioxide System
! in Seawater", DOE, 1994 (SOP No. 3, p25-26).
!
! Derive simple terms used more than once
!
log100 = log(100.0)
do j = jsc, jec !{
do i = isc, iec !{
if (mask(i,j) .ne. 0.0) then !{
tk = 273.15 + t(i,j)
tk100 = tk / 100.0
tk1002 = tk100 * tk100
invtk = 1.0 / tk
logtk = log(tk)
is = 19.924 * s(i,j) / (1000.0 - 1.005 * s(i,j))
is2 = is * is
sqrtis = sqrt(is)
s2 = s(i,j) * s(i,j)
sqrts = sqrt(s(i,j))
s15 = sqrts ** 3
scl = s(i,j) / 1.80655
logf_of_s = log(1.0 - 0.001005 * s(i,j))
!
! k0 from Weiss 1974
!
!k0 = exp(93.4517/tk100 - 60.2409 + 23.3585 * log(tk100) + &
!s(i,j) * (0.023517 - 0.023656 * tk100 + 0.0047036 * tk1002))
!
! k1 = [H][HCO3]/[H2CO3]
! k2 = [H][CO3]/[HCO3]
!
! Millero p.664 (1995) using Mehrbach et al. data on seawater scale
!
k1 = 10.0 ** (-(3670.7 * invtk - 62.008 + 9.7944 * logtk - &
0.0118 * s(i,j) + 0.000116 * s2))
k2 = 10.0 ** (-(1394.7 * invtk + 4.777 - 0.0184 * s(i,j) + 0.000118 * s2))
!
! kb = [H][BO2]/[HBO2]
!
! Millero p.669 (1995) using data from Dickson (1990)
!
kb = exp((-8966.90 - 2890.53 * sqrts - 77.942 * s(i,j) + &
1.728 * s15 - 0.0996 * s2) * invtk + &
(148.0248 + 137.1942 * sqrts + 1.62142 * s(i,j)) + &
(-24.4344 - 25.085 * sqrts - 0.2474 * s(i,j)) * logtk + &
0.053105 * sqrts * tk)
!
! k1p = [H][H2PO4]/[H3PO4]
!
! DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
!
k1p = exp(-4576.752 * invtk + 115.525 - 18.453 * logtk + &
(-106.736 * invtk + 0.69171) * sqrts + &
(-0.65643 * invtk - 0.01844) * s(i,j))
!
! k2p = [H][HPO4]/[H2PO4]
!
! DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
!
k2p = exp(-8814.715 * invtk + 172.0883 - 27.927 * logtk + &
(-160.340 * invtk + 1.3566) * sqrts + &
(0.37335 * invtk - 0.05778) * s(i,j))
!
!-----------------------------------------------------------------------
! k3p = [H][PO4]/[HPO4]
!
! DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
!
k3p = exp(-3070.75 * invtk - 18.141 + &
(17.27039 * invtk + 2.81197) * sqrts + &
(-44.99486 * invtk - 0.09984) * s(i,j))
!
!-----------------------------------------------------------------------
! ksi = [H][SiO(OH)3]/[Si(OH)4]
!
! Millero p.671 (1995) using data from Yao and Millero (1995)
!
ksi = exp(-8904.2 * invtk + 117.385 - 19.334 * logtk + &
(-458.79 * invtk + 3.5913) * sqrtis + &
(188.74 * invtk - 1.5998) * is + &
(-12.1652 * invtk + 0.07871) * is2 + logf_of_s)
!
!-----------------------------------------------------------------------
! kw = [H][OH]
!
! Millero p.670 (1995) using composite data
!
kw = exp(-13847.26 * invtk + 148.9652 - 23.6521 * logtk + &
(118.67 * invtk - 5.977 + 1.0495 * logtk) * sqrts - 0.01615 * s(i,j))
!
!-----------------------------------------------------------------------
! ks = [H][SO4]/[HSO4]
!
! Dickson (1990, J. chem. Thermodynamics 22, 113)
!
ks = exp(-4276.1 * invtk + 141.328 - 23.093 * logtk + &
(-13856.0 * invtk + 324.57 - 47.986 * logtk) * sqrtis + &
(35474.0 * invtk - 771.54 + 114.723 * logtk) * is - &
2698.0 * invtk * sqrtis ** 3 + 1776.0 * invtk * is2 + logf_of_s)
!
!-----------------------------------------------------------------------
! kf = [H][F]/[HF]
!
! Dickson and Riley (1979) -- change pH scale to total
!
kf = exp(1590.2 * invtk - 12.641 + 1.525 * sqrtis + logf_of_s + &
log(1.0 + (0.1400 / 96.062) * scl / ks))
!
!-----------------------------------------------------------------------
! Calculate concentrations for borate, sulfate, and fluoride
!
! Uppstrom (1974)
!
bt = 0.000232 / 10.811 * scl
!
! Morris & Riley (1966)
!
st = 0.14 / 96.062 * scl
!
! Riley (1965)
!
ft = 0.000067 / 18.9984 * scl
!
! set some stuff to pass back, if requested
!
if (present(k1_out)) then !{
k1_out(i,j) = k1
endif !}
if (present(k2_out)) then !{
k2_out(i,j) = k2
endif !}
if (present(invtk_out)) then !{
invtk_out(i,j) = invtk
endif !}
!
!-----------------------------------------------------------------------
! Possibly convert input in mol/m^3 -> mol/kg
!---------------------------------------------------------------------
!
sit = sit_in(i,j) * scale_factor
ta = ta_in(i,j) * scale_factor
dic = dic_in(i,j) * scale_factor
pt = pt_in(i,j) * scale_factor
!
!***********************************************************************
!
! Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
! The solution converges to err of xacc. The solution must be within
! the range x1 to x2.
!
! If DIC and TA are known then either a root finding or iterative method
! must be used to calculate htotal. In this case we use the
! Newton-Raphson "safe" method taken from "Numerical Recipes"
! (function "rtsafe.f" with error trapping removed).
!
! As currently set, this procedure iterates about 12 times. The x1
! and x2 values set below will accomodate ANY oceanographic values.
! If an initial guess of the pH is known, then the number of
! iterations can be reduced to about 5 by narrowing the gap between
! x1 and x2. It is recommended that the first few time steps be run
! with x1 and x2 set as below. After that, set x1 and x2 to the
! previous value of the pH +/- ~0.5. The current setting of xacc will
! result in co2star accurate to 3 significant figures (xx.y). Making
! xacc bigger will result in faster convergence also, but this is not
! recommended (xacc of 10**-9 drops precision to 2 significant
! figures).
!
htotal(i,j) = drtsafe(htotalhi(i,j)*htotal(i,j), &
htotallo(i,j)*htotal(i,j), &
dic, ta, pt, sit, k1, k2, k1p, k2p, k3p, &
bt, ft, st, kb, kw, kf, ks, ksi, xacc)
!
! Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2,
! ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49)
! Convert units of output arguments
! Note: co2star is calculated in
! mol/kg within this routine
!
htotal2 = htotal(i,j) * htotal(i,j)
co2star_internal = dic * htotal2 / (htotal2 + k1 * (htotal(i,j) + k2)) / scale_factor
if (present(co2star)) then
co2star(i,j) = co2star_internal
endif
if (present(co3_ion)) then
co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2
endif
!ph = -log10(htotal(i,j))
!
! f = k0(1-pH2O)*correction term for non-ideality
!
! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6
! values)
!
if (present(alpha) .or. present(pco2surf)) then !{
ff = exp(-162.8301 + 218.2968 / tk100 + &
90.9241 * (logtk - log100) - 1.47696 * tk1002 + &
s(i,j) * (0.025695 - 0.025225 * tk100 + 0.0049867 * tk1002))
!
!---------------------------------------------------------------
! Add two output arguments for storing pCO2surf
! Should we be using K0 or ff for the solubility here?
! Convert pCO2surf from atm to uatm
! possibly convert output units
!---------------------------------------------------------------
!
if (present(alpha)) then !{
alpha(i,j) = ff / scale_factor
endif !}
if (present(pco2surf)) then !{
pCO2surf(i,j) = co2star_internal / (ff / scale_factor) / permeg
endif !}
endif !}
else !}{
if (present(co2star)) then
co2star(i,j) = 0.0
endif
if (present(k1_out)) then !{
k1_out(i,j) = 0.0
endif !}
if (present(k2_out)) then !{
k2_out(i,j) =0.0
endif !}
if (present(invtk_out)) then !{
invtk_out(i,j) = 0.0
endif !}
if (present(alpha)) then !{
alpha(i,j) = 0.0
endif !}
if (present(pco2surf)) then !{
pCO2surf(i,j) = 0.0
endif !}
endif !}
enddo !} i
enddo !} j
return
end subroutine ocmip2_co2calc !}
! NAME="ocmip2_co2calc"
!#######################################################################
!
!
!
! File taken from Numerical Recipes. Modified R. M. Key 4/94
!
function drtsafe(x1, x2, dic, ta, pt, sit, k1, k2, k1p, k2p, k3p, &
bt, ft, st, kb, kw, kf, ks, ksi, xacc) !{
implicit none
!
! arguments
!
real :: drtsafe
real, intent(in) :: x1
real, intent(in) :: x2
real, intent(in) :: dic
real, intent(in) :: ta
real, intent(in) :: pt
real, intent(in) :: sit
real, intent(in) :: k1
real, intent(in) :: k2
real, intent(in) :: k1p
real, intent(in) :: k2p
real, intent(in) :: k3p
real, intent(in) :: bt
real, intent(in) :: ft
real, intent(in) :: st
real, intent(in) :: kb
real, intent(in) :: kw
real, intent(in) :: kf
real, intent(in) :: ks
real, intent(in) :: ksi
real, intent(in) :: xacc
!
! local parameters
!
integer, parameter :: maxit = 100
!
! local variables
!
integer :: j
real :: fl
real :: df
real :: fh
real :: swap
real :: xl
real :: xh
real :: dxold
real :: dx
real :: f
real :: temp
call ta_iter_1(x1, fl, df, dic, ta, pt, sit, k1, k2, &
k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi)
call ta_iter_1(x2, fh, df, dic, ta, pt, sit, k1, k2, &
k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi)
if(fl .lt. 0.0) then
xl=x1
xh=x2
else
xh=x1
xl=x2
swap=fl
fl=fh
fh=swap
end if
drtsafe=0.5*(x1+x2)
dxold=abs(x2-x1)
dx=dxold
call ta_iter_1(drtsafe, f, df, dic, ta, pt, sit, k1, k2, &
k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi)
do j=1,maxit !{
if (((drtsafe-xh)*df-f)*((drtsafe-xl)*df-f) .ge. 0.0 .or. &
abs(2.0*f) .gt. abs(dxold*df)) then
dxold=dx
dx=0.5*(xh-xl)
drtsafe=xl+dx
if (xl .eq. drtsafe) then
! write (6,*) 'Exiting drtsafe at A on iteration ', j, ', ph = ', -log10(drtsafe)
return
endif
else
dxold=dx
dx=f/df
temp=drtsafe
drtsafe=drtsafe-dx
if (temp .eq. drtsafe) then
! write (6,*) 'Exiting drtsafe at B on iteration ', j, ', ph = ', -log10(drtsafe)
return
endif
end if
if (abs(dx) .lt. xacc) then
! write (6,*) 'Exiting drtsafe at C on iteration ', j, ', ph = ', -log10(drtsafe)
return
endif
call ta_iter_1(drtsafe, f, df, dic, ta, pt, sit, k1, k2, &
k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi)
if(f .lt. 0.0) then
xl=drtsafe
fl=f
else
xh=drtsafe
fh=f
end if
enddo !} j
return
end function drtsafe !}
! NAME="drtsafe"
!#######################################################################
!
!
!
! This routine expresses TA as a function of DIC, htotal and constants.
! It also calculates the derivative of this function with respect to
! htotal. It is used in the iterative solution for htotal. In the call
! "x" is the input value for htotal, "fn" is the calculated value for TA
! and "df" is the value for dTA/dhtotal
!
subroutine ta_iter_1(x, fn, df, dic, ta, pt, sit, k1, k2, &
k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi) !{
implicit none
!
! arguments
!
real, intent(in) :: x
real, intent(out) :: fn
real, intent(out) :: df
real, intent(in) :: dic
real, intent(in) :: ta
real, intent(in) :: pt
real, intent(in) :: sit
real, intent(in) :: k1
real, intent(in) :: k2
real, intent(in) :: k1p
real, intent(in) :: k2p
real, intent(in) :: k3p
real, intent(in) :: bt
real, intent(in) :: ft
real, intent(in) :: st
real, intent(in) :: kb
real, intent(in) :: kw
real, intent(in) :: kf
real, intent(in) :: ks
real, intent(in) :: ksi
!
! local variables
!
real :: x2
real :: x3
real :: k12
real :: k12p
real :: k123p
real :: c
real :: a
real :: a2
real :: da
real :: b
real :: b2
real :: db
x2 = x*x
x3 = x2*x
k12 = k1*k2
k12p = k1p*k2p
k123p = k12p*k3p
c = 1.0 + st/ks
a = x3 + k1p*x2 + k12p*x + k123p
a2 = a*a
da = 3.0*x2 + 2.0*k1p*x + k12p
b = x2 + k1*x + k12
b2 = b*b
db = 2.0*x + k1
!
! fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree+hso4+hf+h3po4-ta
!
fn = k1*x*dic/b + 2.0*dic*k12/b + bt/(1.0 + x/kb) + kw/x + &
pt*k12p*x/a + 2.0*pt*k123p/a + sit/(1.0 + x/ksi) - &
x/c - st/ (1.0 + ks/x/c) - ft/(1.0 + kf/x) - &
pt*x3/a - ta
!
! df = dfn/dx
!
df = ((k1*dic*b) - k1*x*dic*db)/b2 - 2.0*dic*k12*db/b2 - &
bt/kb/(1.0+x/kb)**2 - kw/x2 + (pt*k12p*(a - x*da))/a2 - &
2.0*pt*k123p*da/a2 - sit/ksi/(1.0+x/ksi)**2 - 1.0/c + &
st*(1.0 + ks/x/c)**(-2)*(ks/c/x2) + ft*(1.0 + kf/x)**(-2)*kf/x2 - &
pt*x2*(3.0*a-x*da)/a2
return
end subroutine ta_iter_1 !}
! NAME="ta_iter_1"
end module ocmip2_co2calc_mod !}