MODULE asmpco2bal !Contains subroutine to convert a global pco2 increment to independent DIC and alkalinity increments USE par_kind, ONLY: wp USE par_oce, ONLY: jpi, jpj, jpk USE dom_oce, ONLY: tmask #if defined key_hadocc USE had_bgc_const, ONLY: RHO_WATER_SI #else USE phycst, ONLY: RHO_WATER_SI => rau0 #endif IMPLICIT NONE PRIVATE PUBLIC asm_pco2_bal CONTAINS SUBROUTINE asm_pco2_bal(pco2inc, tco2_tracer_bkg, alk_tracer_bkg, & & tem_tracer_bkg, sal_tracer_bkg, & & tco2_tracer_inc, alk_tracer_inc) !Converts a global pco2 increment to independent DIC and alkalinity increments !Calculates the solution to the equation [DDIC, DAlk]=(DpCO2 / |Grad (pCO2)|)* Grad(pCO2); where Grad(pco2)=[dpCO2/dDIC, dpCO2/dalk] !Gradient code pco2_grad is a modification of assimcd1.mf77 by Ian Totterdell. !James While 06/09 IMPLICIT NONE REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pco2inc !pCO2 increment REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg, & !In: tCO2 (DIC) background alk_tracer_bkg !In alkalinity background REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg, & !In: temperature background sal_tracer_bkg !In salinity background REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: tco2_tracer_inc, & !Out: tCO2 (DIC) increment alk_tracer_inc !Out alkalinity !local variables REAL(wp), DIMENSION(jpi,jpj) :: Dpco2_tco2, Dpco2_alk REAL(wp) :: wgt_fct INTEGER :: ji,jj !calculate the gradients with respect to DIC and Alkalinity CALL pco2_grad(tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, Dpco2_tco2, Dpco2_alk) !loop over the ocean grid calculating the increments to tCO2 (DIC) and alkalinity DO jj=1, jpj DO ji=1, jpi wgt_fct=pco2inc(ji,jj)/( (Dpco2_tco2(ji,jj)*Dpco2_tco2(ji,jj)) + (Dpco2_alk(ji,jj)*Dpco2_alk(ji,jj)) ) tco2_tracer_inc(ji,jj)=wgt_fct * Dpco2_tco2(ji,jj) * tmask(ji,jj,1) alk_tracer_inc(ji,jj)=wgt_fct * Dpco2_alk(ji,jj) * tmask(ji,jj,1) ENDDO ENDDO END SUBROUTINE asm_pco2_bal SUBROUTINE pco2_grad( tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, DPDC, DPDA ) ! !-------M---------M---------M---------M---------M---------M---------M-- ! ! Ian J Totterdell ! 2007 August 08 ! ! Code for calculating derivatives of pCO2 wrt DIC and Alk ! ! J While ! June 2009 ! Modiefied for use in NEMO !-------M---------M---------M---------M---------M---------M---------M-- IMPLICIT NONE ! ! REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg, & !In: tCO2 (DIC) background alk_tracer_bkg !In alkalinity background REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg, & !In: temperature background sal_tracer_bkg !In salinity background ! REAL(wp), INTENT(OUT), DIMENSION(jpi,jpj) :: & DPDC, & DPDA ! ! ! Local variables. ! INTEGER :: & mxiter, & ji,jj, & iter ! REAL(wp) :: & a1_co2, a2_co2, a3_co2, & b1_co2, b2_co2, b3_co2, & a1_k1, a2_k1, a3_k1, & b1_k1, b2_k1, b3_k1, b4_k1, b5_k1, & a1_k2, a2_k2, a3_k2, & b1_k2, b2_k2, b3_k2, b4_k2, b5_k2, & a1_kb, a2_kb, a3_kb, a4_kb, a5_kb, & b1_kb, b2_kb, b3_kb, & c1_kb, c2_kb, c3_kb, d_kb, & a1_kw, a2_kw, a3_kw, & b1_kw, b2_kw, b3_kw, c_kw ! REAL(wp) :: & tk, logtk, ttk, & sal, sqsal, & kone, ktwo, kboric, kwater, & alkprime, qa,qb,qc, & fxa, xtrma, xtrmb, & dxda, dxdc, & dgda, dgdc ! REAL(wp), DIMENSION(jpi,jpj) :: & co2_sol, & rek, & rg1, & rg2, & rg3, & boron, & tco2, & alk_t, & xs, & x1, & x2, & x2p, & x3, & xx, & xf ! ! Set parameter and constant values. ! PARAMETER ( mxiter = 10 ) ! ! CO2 solubility, constants for the calculation of. ! From Weiss (1974) ! DATA & a1_co2 / -60.2409 /, a2_co2 / 93.4517 /, a3_co2 / 23.3585 /, & b1_co2 / 0.023517 /, b2_co2 / -0.023656 /, b3_co2 / 0.0047036 / ! ! kone, constants for the calculation of. ! From Roy et al (1993) ! DATA a1_k1 /-2307.1266/, a2_k1 / 2.83655 /, a3_k1 /-1.5529413/ & , b1_k1 / -4.0484 /, b2_k1 /-0.20760841/, b3_k1 /0.08468345/ & , b4_k1 /-0.00654208/, b5_k1 / -0.001005 / ! ! ktwo, constants for the calculation of. ! From Roy et al (1993) ! DATA a1_k2 /-3351.6106/, a2_k2 / -9.226508 /, a3_k2 /-0.2005743/ & , b1_k2 / -23.9722 /, b2_k2 /-0.106901773/,b3_k2 /0.1130822/ & , b4_k2 /-0.00846934/, b5_k2 / -0.001005 / ! ! kboric, constants for the calculation of. ! Using data from Dickson (1990); ! DATA a1_kb / -8966.90 /, a2_kb / -2890.53 /, a3_kb / -77.942 / & , a4_kb / 1.728 /, a5_kb /-0.0996/ & , b1_kb / 148.0248 /, b2_kb / 137.1942 /, b3_kb / 1.62142 / & , c1_kb / -24.4344 /, c2_kb / -25.085 /, c3_kb / -0.2474 / & , d_kb / 0.053105 / ! ! kwater, constants for the calculation of. ! From Millero (1995, p670), using composite data. ! DATA a1_kw / -13847.26 /, a2_kw / 148.9652 /, a3_kw / -23.6521 / & , b1_kw / 118.67 /, b2_kw / -5.977 /, b3_kw / 1.0495 / & , c_kw / -0.01615 / ! ! Executable code begins here. ! ! Calculate the carbon system constants. ! DO ji=1,jpi DO jj=1,jpj ! tk = 273.15 + MAX( -2.0, MIN( 40.0, tem_tracer_bkg(ji,jj) ) ) logtk = LOG(tk) ttk = MAX( 2.71, 0.01 * ( 273.15 + tem_tracer_bkg(ji,jj) ) ) sal = MAX( 0.0, MIN( 45.0, sal_tracer_bkg(ji,jj) ) ) sqsal = SQRT(sal) ! ! Calculate the CO2 solubility from Weiss (1974) ! Units are mol/(kg.soln)/atm . ! Check value (Zeebe and Wolf-G; tc=25,sal=35): ln(co2_sol) = -3.5617 ! co2_sol(ji,jj) = EXP( & a1_co2 + a2_co2 / ttk + a3_co2 * LOG(ttk) + & sal * ( b1_co2 + ttk * ( b2_co2 + b3_co2 * ttk ) ) ) ! ! Calculate the 1st constant of carbonic acid from Roy et al. (1993) ! Units are mol/(kg.soln) on the total hydrogen ion scale. ! K1 = (HCO3-).(H+)/(CO2(g)) ! Check value (Dickson and Goyet; tc=25,sal=35): ln(K1) = -13.4847 ! kone & = 1.0e-7 *( 1.0 + b5_k1*sal ) * EXP( a1_k1/tk + a2_k1 + a3_k1*logtk & + ( b1_k1/tk + b2_k1 )*sqsal + b3_k1*sal + b4_k1*sal*sqsal ) ! ! Calculate the 2nd constant of carbonic acid from Roy et al. (1993) ! Units are mol/(kg.soln) on the total hydrogen ion scale. ! K2 = (CO3--).(H+)/(HCO3-) ! Check value (Dickson and Goyet; tc=25,sal=35): ln(K2) = -20.5504 ! ktwo & = 1.0e-7 *( 1.0 + b5_k2*sal ) * EXP( a1_k2/tk + a2_k2 + a3_k2*logtk & + ( b1_k2/tk + b2_k2 )*sqsal + b3_k2*sal + b4_k2*sal*sqsal ) ! ! Calculate the equilibrium constant of boric acid from Dickson (1990). ! Units are mol/(kg.soln) on the total hydrogen ion scale. ! Kb = (B(OH)4-).(H+)/(B(OH)3) (with (Bt) = (B(OH)3) + (B(OH)4-) ). ! Check value (Dickson and Goyet; tc=25,sal=35): ln(Kb) = -19.7964 ! kboric & = 1.0e-7 * EXP( ( a1_kb + a2_kb*sqsal + a3_kb*sal + a4_kb*sal*sqsal & + a5_kb*sal*sal )/tk + ( b1_kb + b2_kb*sqsal + b3_kb*sal ) & + ( c1_kb + c2_kb*sqsal + c3_kb*sal )*logtk + d_kb*sqsal*tk ) ! ! Calculate the equilibrium constant of water from Millero (1994). ! Units are ( mol/(kg.soln) )^2 on the total hydrogen ion scale. ! Kw = (OH-).(H+) ! Check value (Dickson and Goyet; tc=25,sal=35): ln(Kw) = -30.434 ! kwater & = 1.0e-14 *EXP( ( a1_kw/tk + a2_kw + a3_kw*logtk ) & + ( b1_kw/tk + b2_kw + b3_kw*logtk )*sqsal + c_kw*sal ) ! ! Calculate the combinations (as in PPCO2). ! rek(ji,jj) = SQRT( kone / ktwo ) rg3(ji,jj) = SQRT( kone * ktwo ) rg1(ji,jj) = rg3(ji,jj) / kboric rg2(ji,jj) = kwater / rg3(ji,jj) ! ! Set Boron concentration according to the surface salinity ! (Boron accounts for about 10% of the total Alkalinity). ! Use relationship used by Peng et al, Tellus 39b 1987, ! originally from Culkin 1965: ! ! boron(ji,jj) = 4.16e-4 * sal / 35.0 ! Correct value but not in UM boron(ji,jj) = 4.106e-4 * sal / 35.0 ! ! tco2 and alk_t are got from the modelled tracers at the surface ! (converted to units of mol/kg). ! tco2(ji,jj) = MAX( 0.0, tco2_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI ) alk_t(ji,jj) = MAX( 0.0, alk_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI ) ! ! Set initial values for iterative method. ! xs(ji,jj) = 1.0 x2p(ji,jj) = 0.0 x3(ji,jj) = 0.0 xx(ji,jj) = 0.0 xf(ji,jj) = xs(ji,jj) ! ENDDO ENDDO ! ! Use the secant method of similar triangles to calculate the H+ ion ! concentration, [H] (though strictly xf refers to rg3/[H] ). ! To aid vectorisation this has been written to iterate a set ! number of times (mxiter), large enough to bring every example to ! within the 1e-5 fractional change criterion. ! ! The method follows that of Bacastow (1990) in separating the ! alkalinity into two parts: the major part is the carbonate alk ! which uses the [H] that we are trying to find, while the minor ! part (xw) is sum of the borate and water dissociation terms which ! use the previous approximation to [H]. ! This leaves an easy-to-solve quadratic. ! DO iter=1,mxiter ! DO ji=1,jpi DO jj=1,jpj ! IF((abs(xf(ji,jj)-xx(ji,jj))/xf(ji,jj))<1.e-5) CYCLE !prevent iteration collapsing to 0.0/0.0 x1(ji,jj) = x2p(ji,jj) x2(ji,jj) = x3(ji,jj) x2p(ji,jj) = xx(ji,jj) x3(ji,jj) = xf(ji,jj) IF (iter<=2) THEN xx(ji,jj) = xf(ji,jj) ELSE xx(ji,jj) = x2p(ji,jj) + ( x3(ji,jj) - x2p(ji,jj) )*( x1(ji,jj) - x2p(ji,jj) ) & /( x1(ji,jj) - x2(ji,jj) - x2p(ji,jj) + x3(ji,jj) ) ENDIF ! ! All the above has been setting up the similar triangles for the ! secant method and using it to find an improved interpolated ! estimate; it has been completely independent of the equations ! we are trying to solve. ! ! The next bit, however, IS dependent on the equations. ! alkprime is the carbonate alkalinity, which is calculated ! from the total alkalinity by subtracting the other, smaller, terms, ! which use the interpolated previous approximation to xf (xx). ! alkprime = alk_t(ji,jj) - rg2(ji,jj)*xx(ji,jj) + rg3(ji,jj)/xx(ji,jj) & - boron(ji,jj)/( 1.0 + rg1(ji,jj)/xx(ji,jj) ) ! ! Set up the quadratic terms. ! qa = 2.0*tco2(ji,jj) - alkprime qb = -rek(ji,jj)*( alkprime - tco2(ji,jj) ) qc = -alkprime xf(ji,jj) = ( -qb + SQRT( qb*qb - 4.0*qa*qc ) )*0.5/qa ENDDO ENDDO ! ENDDO ! DO ji=1,jpi DO jj=1,jpj ! xtrma = 1.0/( 1.0 + rek(ji,jj)*xf(ji,jj) + xf(ji,jj)*xf(ji,jj) ) xtrmb = rg2(ji,jj) + rg3(ji,jj)/(xf(ji,jj)*xf(ji,jj)) & + rg1(ji,jj)*boron(ji,jj)/(( rg1(ji,jj) + xf(ji,jj) )*( rg1(ji,jj) + xf(ji,jj) )) ! dxda = 1.0/( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) & + rek(ji,jj)*xf(ji,jj)*xf(ji,jj) )*xtrma*xtrma + xtrmb ) dxdc = -dxda*xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma ! dgdc = xtrma - dxdc*tco2(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma*xtrma dgda = ( 1.0 - dxda*( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) )*xtrma & + xtrmb ) )/(xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )) ! fxa = 1.0e-3*RHO_WATER_SI*co2_sol(ji,jj) DPDC(ji,jj) = dgdc/fxa DPDA(ji,jj) = dgda/fxa ! ENDDO ENDDO RETURN END SUBROUTINE pco2_grad END MODULE asmpco2bal