MODULE co2calc !----------------------------------------------------------------------------- ! based upon OCMIP2 co2calc ! ! CVS:$Id: co2calc.F90 941 2006-05-12 21:36:48Z klindsay $ ! CVS:$Name$ !----------------------------------------------------------------------------- USE constants USE blocks, ONLY : nx_block, ny_block, block, get_block USE domain, ONLY : blocks_clinic USE domain_size, ONLY : max_blocks_clinic USE kinds_mod USE state_mod, ONLY : ref_pressure USE io_types, ONLY : stdout USE time_management, ONLY : nsteps_run #ifdef CCSMCOUPLED !*** ccsm USE shr_vmath_mod #endif IMPLICIT NONE !----------------------------------------------------------------------------- ! public/private declarations !----------------------------------------------------------------------------- PRIVATE PUBLIC :: co2calc_row, comp_CO3terms, comp_co3_sat_vals !----------------------------------------------------------------------------- ! module parameters !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! The current setting of xacc, a tolerance critera, will result in co2star ! being 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). !----------------------------------------------------------------------------- REAL(KIND=r8), PARAMETER :: xacc = 1e-10_r8 INTEGER(KIND=int_kind), PARAMETER :: max_bracket_grow_it = 3 INTEGER(KIND=int_kind), PARAMETER :: maxit = 100 !----------------------------------------------------------------------------- ! declarations for function coefficients & species concentrations !----------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block,max_blocks_clinic) :: & kw, kb, ks, kf, k1p, k2p, k3p, ksi, & bt, st, ft, dic, ta, pt, sit !***************************************************************************** CONTAINS !***************************************************************************** SUBROUTINE co2calc_row(iblock, j, mask, locmip_k1_k2_bug_fix, t, s, & dic_in, ta_in, pt_in, sit_in, phlo, phhi, ph, xco2_in, atmpres, & co2star, dco2star, pCO2surf, dpco2) !--------------------------------------------------------------------------- ! SUBROUTINE co2calc_row ! ! PURPOSE : Calculate delta co2*, etc. from total alkalinity, total CO2, ! temp, salinity (s), etc. !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: iblock, j LOGICAL(KIND=log_kind), DIMENSION(nx_block), INTENT(IN) :: mask LOGICAL(KIND=log_kind), INTENT(IN) :: locmip_k1_k2_bug_fix REAL(KIND=r8), DIMENSION(nx_block), INTENT(IN) :: & t, & ! temperature (degrees C) s, & ! salinity (PSU) dic_in, & ! total inorganic carbon (nmol/cm^3) ta_in, & ! total alkalinity (neq/cm^3) pt_in, & ! inorganic phosphate (nmol/cm^3) sit_in, & ! inorganic silicate (nmol/cm^3) xco2_in, & ! atmospheric mole fraction CO2 in dry air (ppmv) atmpres ! atmospheric pressure (atmosphere) !--------------------------------------------------------------------------- ! input/output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(INOUT) :: & phlo, & ! lower limit of pH range phhi ! upper limit of pH range !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(OUT) :: & ph, & ! computed ph values, for initial guess on next time step co2star, & ! CO2*water (nmol/cm^3) dco2star, & ! delta CO2 (nmol/cm^3) pco2surf, & ! oceanic pCO2 (ppmv) dpco2 ! Delta pCO2, i.e, pCO2ocn - pCO2atm (ppmv) !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- INTEGER(KIND=int_kind) :: i INTEGER(KIND=int_kind) :: k REAL(KIND=r8) :: & mass_to_vol, & ! (mol/kg) -> (mmol/m^3) vol_to_mass, & ! (mmol/m^3) -> (mol/kg) co2starair, & ! co2star saturation htotal2 REAL(KIND=r8), DIMENSION(nx_block) :: & xco2, & ! atmospheric CO2 (atm) htotal, & ! free concentration of H ion k0,k1,k2, & ! equilibrium constants for CO2 species ff ! fugacity of CO2 !--------------------------------------------------------------------------- ! check for existence of ocean points !--------------------------------------------------------------------------- IF (COUNT(mask) == 0) THEN ph = c0 co2star = c0 dco2star = c0 pCO2surf = c0 dpCO2 = c0 RETURN END IF !--------------------------------------------------------------------------- ! set unit conversion factors !--------------------------------------------------------------------------- mass_to_vol = 1e6_r8 * rho_sw vol_to_mass = c1 / mass_to_vol !--------------------------------------------------------------------------- ! convert tracer units to per mass & xco2 from uatm to atm !--------------------------------------------------------------------------- DO i = 1,nx_block IF (mask(i)) THEN dic(i,iblock) = dic_in(i) * vol_to_mass xco2(i) = xco2_in(i) * 1e-6_r8 END IF ! if mask END DO ! i loop !--------------------------------------------------------------------------- ! compute htotal !--------------------------------------------------------------------------- k = 1 CALL comp_htotal_kvals(iblock, j, k, mask, t, s, dic_in, ta_in, pt_in, sit_in, & phlo, phhi, htotal, k0, k1, k2, ff, & k1_k2_pH_tot=locmip_k1_k2_bug_fix) !--------------------------------------------------------------------------- ! Compute co2starair !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! 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) !--------------------------------------------------------------------------- DO i = 1,nx_block IF (mask(i)) THEN htotal2 = htotal(i) ** 2 co2star(i) = dic(i,iblock) * htotal2 / & (htotal2 + k1(i) * htotal(i) + k1(i) * k2(i)) co2starair = xco2(i) * ff(i) * atmpres(i) dco2star(i) = co2starair - co2star(i) ph(i) = -LOG10(htotal(i)) !--------------------------------------------------------------------- ! Add two output arguments for storing pCO2surf ! Should we be using K0 or ff for the solubility here? !--------------------------------------------------------------------- pCO2surf(i) = co2star(i) / ff(i) dpCO2(i) = pCO2surf(i) - xco2(i) * atmpres(i) !--------------------------------------------------------------------- ! Convert units of output arguments ! Note: pCO2surf and dpCO2 are calculated in atm above. !--------------------------------------------------------------------- co2star(i) = co2star(i) * mass_to_vol dco2star(i) = dco2star(i) * mass_to_vol pCO2surf(i) = pCO2surf(i) * 1e6_r8 dpCO2(i) = dpCO2(i) * 1e6_r8 ELSE ! if mask ph(i) = c0 co2star(i) = c0 dco2star(i) = c0 pCO2surf(i) = c0 dpCO2(i) = c0 END IF ! if mask END DO ! i loop END SUBROUTINE co2calc_row !***************************************************************************** SUBROUTINE comp_CO3terms(iblock, k, mask, t, s, dic_in, ta_in, pt_in, sit_in, & phlo, phhi, ph, H2CO3, HCO3, CO3) !--------------------------------------------------------------------------- ! SUBROUTINE comp_CO3terms ! ! PURPOSE : Calculate H2CO3, HCO3, CO3 from ! total alkalinity, total CO2, temp, salinity (s), etc. !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: iblock INTEGER(KIND=int_kind), INTENT(IN) :: k LOGICAL(KIND=log_kind), DIMENSION(nx_block, ny_block), INTENT(IN) :: mask REAL(KIND=r8), DIMENSION(nx_block,ny_block), INTENT(IN) :: & t, & ! temperature (degrees C) s, & ! salinity (PSU) dic_in, & ! total inorganic carbon (nmol/cm^3) ta_in, & ! total alkalinity (neq/cm^3) pt_in, & ! inorganic phosphate (nmol/cm^3) sit_in ! inorganic silicate (nmol/cm^3) !--------------------------------------------------------------------------- ! input/output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block,ny_block), INTENT(INOUT) :: & phlo, & ! lower limit of pH range phhi ! upper limit of pH range !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block, ny_block), INTENT(OUT) :: & pH, & ! computed ph values, for initial guess on next time step H2CO3, & ! Carbonic Acid Concentration HCO3, & ! Bicarbonate Ion Concentration CO3 ! Carbonate Ion Concentration !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- INTEGER(KIND=int_kind) :: i, j REAL(KIND=r8) :: & mass_to_vol, & ! (mol/kg) -> (mmol/m^3) vol_to_mass, & ! (mmol/m^3) -> (mol/kg) htotal2, denom REAL(KIND=r8), DIMENSION(nx_block) :: & htotal, & ! free concentration of H ion k0,k1,k2, & ! equilibrium constants for CO2 species ff ! fugacity of CO2 !--------------------------------------------------------------------------- ! check for existence of ocean points !--------------------------------------------------------------------------- IF (COUNT(mask) == 0) THEN ph = c0 H2CO3 = c0 HCO3 = c0 CO3 = c0 RETURN END IF !--------------------------------------------------------------------------- ! set unit conversion factors !--------------------------------------------------------------------------- mass_to_vol = 1e6_r8 * rho_sw vol_to_mass = c1 / mass_to_vol DO j = 1,ny_block !------------------------------------------------------------------------ ! check for existence of ocean points on this row !------------------------------------------------------------------------ IF (COUNT(mask(:,j)) == 0) THEN ph(:,j) = c0 H2CO3(:,j) = c0 HCO3(:,j) = c0 CO3(:,j) = c0 CYCLE END IF !------------------------------------------------------------------------ ! convert tracer units to per mass !------------------------------------------------------------------------ DO i = 1,nx_block IF (mask(i,j)) THEN dic(i,iblock) = dic_in(i,j) * vol_to_mass END IF ! if mask END DO ! i loop !------------------------------------------------------------------------ ! compute htotal !------------------------------------------------------------------------ CALL comp_htotal_kvals(iblock, j, k, mask(:,j), t(:,j), s(:,j), dic_in(:,j), & ta_in(:,j), pt_in(:,j), sit_in(:,j), & phlo(:,j), phhi(:,j), htotal, k0, k1, k2, ff) !------------------------------------------------------------------------ ! 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-51) !------------------------------------------------------------------------ DO i = 1,nx_block IF (mask(i,j)) THEN htotal2 = htotal(i) ** 2 denom = c1 / (htotal2 + k1(i) * htotal(i) + k1(i) * k2(i)) H2CO3(i,j) = dic(i,iblock) * htotal2 * denom HCO3(i,j) = dic(i,iblock) * k1(i) * htotal(i) * denom CO3(i,j) = dic(i,iblock) * k1(i) * k2(i) * denom ph(i,j) = -LOG10(htotal(i)) !------------------------------------------------------------------ ! Convert units of output arguments !------------------------------------------------------------------ H2CO3(i,j) = H2CO3(i,j) * mass_to_vol HCO3(i,j) = HCO3(i,j) * mass_to_vol CO3(i,j) = CO3(i,j) * mass_to_vol ELSE ! if mask ph(i,j) = c0 H2CO3(i,j) = c0 HCO3(i,j) = c0 CO3(i,j) = c0 END IF ! if mask END DO ! i loop END DO ! j loop END SUBROUTINE comp_CO3terms !***************************************************************************** SUBROUTINE comp_co3_sat_vals(k, mask, t, s, co3_sat_calc, co3_sat_arag) !--------------------------------------------------------------------------- ! SUBROUTINE comp_co3_sat_vals ! ! PURPOSE : Calculate co3 concentration at calcite and aragonite saturation ! from temp, salinity (s), press !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: k LOGICAL(KIND=log_kind), DIMENSION(nx_block, ny_block), INTENT(IN) :: mask REAL(KIND=r8), DIMENSION(nx_block,ny_block), INTENT(IN) :: & t, & ! temperature (degrees C) s ! salinity (PSU) !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block, ny_block), INTENT(OUT) :: & co3_sat_calc,&! co3 concentration at calcite saturation co3_sat_arag ! co3 concentration at aragonite saturation !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- INTEGER(KIND=int_kind) :: i, j REAL(KIND=r8) :: & mass_to_vol, & ! (mol/kg) -> (mmol/m^3) press_bar ! pressure at level k [bars] REAL(KIND=r8), DIMENSION(nx_block) :: & tk, & ! temperature (K) log10tk,invtk,sqrts,s15,invRtk,arg,& K_calc, & ! thermodynamic constant for calcite K_arag, & ! thermodynamic constant for aragonite deltaV,Kappa, & ! pressure correction terms lnKfac,Kfac, & ! pressure correction terms inv_Ca ! inverse of Calcium concentration (mol/kg) !--------------------------------------------------------------------------- ! check for existence of ocean points !--------------------------------------------------------------------------- IF (COUNT(mask) == 0) THEN co3_sat_calc = c0 co3_sat_arag = c0 RETURN END IF !--------------------------------------------------------------------------- ! set unit conversion factors !--------------------------------------------------------------------------- mass_to_vol = 1e6_r8 * rho_sw press_bar = ref_pressure(k) DO j = 1,ny_block !------------------------------------------------------------------------ ! check for existence of ocean points on this row !------------------------------------------------------------------------ IF (COUNT(mask(:,j)) == 0) THEN co3_sat_calc(:,j) = c0 co3_sat_arag(:,j) = c0 CYCLE END IF tk = T0_Kelvin + t(:,j) #ifdef CCSMCOUPLED CALL shr_vmath_log(tk, log10tk, nx_block) #else log10tk = LOG(tk) #endif log10tk = log10tk/LOG(c10) invtk = c1 / tk invRtk = (c1 / 83.1451_r8) * invtk #ifdef CCSMCOUPLED CALL shr_vmath_sqrt(s(:,j), sqrts, nx_block) #else sqrts = SQRT(s(:,j)) #endif s15 = sqrts * s(:,j) !------------------------------------------------------------------------ ! Compute K_calc, K_arag, apply pressure factor ! Mucci, Amer. J. of Science 283:781-799, 1983 & Millero 1979 !------------------------------------------------------------------------ arg = -171.9065_r8 - 0.077993_r8 * tk + 2839.319_r8 * invtk + 71.595_r8 * log10tk + & (-0.77712_r8 + 0.0028426_r8 * tk + 178.34_r8 * invtk) * sqrts - & 0.07711_r8 * s(:,j) + 0.0041249_r8 * s15 arg = LOG(c10) * arg #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, K_calc, nx_block) #else K_calc = EXP(arg) #endif IF (k > 1) THEN deltaV = -48.76_r8 + 0.5304_r8 * t(:,j) Kappa = (-11.76_r8 + 0.3692_r8 * t(:,j)) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif K_calc = K_calc * Kfac END IF arg = -171.945_r8 - 0.077993_r8 * tk + 2903.293_r8 * invtk + 71.595_r8 * log10tk + & (-0.068393_r8 + 0.0017276_r8 * tk + 88.135_r8 * invtk) * sqrts - & 0.10018_r8 * s(:,j) + 0.0059415_r8 * s15 arg = LOG(c10) * arg #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, K_arag, nx_block) #else K_arag = EXP(arg) #endif IF (k > 1) THEN deltaV = deltaV + 2.8_r8 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif K_arag = K_arag * Kfac END IF WHERE ( mask(:,j) ) !------------------------------------------------------------------ ! Compute CO3 concentration at calcite & aragonite saturation !------------------------------------------------------------------ inv_Ca = (35.0_r8 / 0.01028_r8) / s(:,j) co3_sat_calc(:,j) = K_calc * inv_Ca co3_sat_arag(:,j) = K_arag * inv_Ca !------------------------------------------------------------------ ! Convert units of output arguments !------------------------------------------------------------------ co3_sat_calc(:,j) = co3_sat_calc(:,j) * mass_to_vol co3_sat_arag(:,j) = co3_sat_arag(:,j) * mass_to_vol ELSEWHERE co3_sat_calc(:,j) = c0 co3_sat_arag(:,j) = c0 END WHERE END DO ! j loop END SUBROUTINE comp_co3_sat_vals !***************************************************************************** SUBROUTINE comp_htotal_kvals(iblock, j, k, mask, t, s, dic_in, ta_in, pt_in, sit_in, & phlo, phhi, htotal, k0, k1, k2, ff, k1_k2_pH_tot) !--------------------------------------------------------------------------- ! SUBROUTINE comp_htotal_kvals ! ! PURPOSE : Calculate htotal, k0, k1, k2, ff from total alkalinity, total CO2, ! temp, salinity (s), etc. !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: iblock, j INTEGER(KIND=int_kind), INTENT(IN) :: k LOGICAL(KIND=log_kind), DIMENSION(nx_block), INTENT(IN) :: mask REAL(KIND=r8), DIMENSION(nx_block), INTENT(IN) :: & t, & ! temperature (degrees C) s, & ! salinity (PSU) dic_in, & ! total inorganic carbon (nmol/cm^3) ta_in, & ! total alkalinity (neq/cm^3) pt_in, & ! inorganic phosphate (nmol/cm^3) sit_in ! inorganic silicate (nmol/cm^3) LOGICAL(KIND=log_kind), INTENT(IN), OPTIONAL :: k1_k2_pH_tot !--------------------------------------------------------------------------- ! input/output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(INOUT) :: & phlo, & ! lower limit of pH range phhi ! upper limit of pH range !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(OUT) :: & htotal, & ! free concentration of H ion k0,k1,k2, & ! equilibrium constants for CO2 species ff ! fugacity of CO2 !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- INTEGER(KIND=int_kind) :: i REAL(KIND=r8) :: & mass_to_vol, & ! (mol/kg) -> (mmol/m^3) vol_to_mass ! (mmol/m^3) -> (mol/kg) REAL(KIND=r8), DIMENSION(nx_block) :: & x1, x2 ! bounds on htotal for solver LOGICAL(KIND=log_kind) :: k1_k2_pH_tot_loc !--------------------------------------------------------------------------- ! check for existence of ocean points !--------------------------------------------------------------------------- IF (COUNT(mask) == 0) THEN htotal = c0 RETURN END IF !--------------------------------------------------------------------------- ! set unit conversion factors !--------------------------------------------------------------------------- mass_to_vol = 1e6_r8 * rho_sw vol_to_mass = c1 / mass_to_vol !--------------------------------------------------------------------------- ! which pH scale is to be used for k1 & k2 constants !--------------------------------------------------------------------------- IF (PRESENT(k1_k2_pH_tot)) THEN k1_k2_pH_tot_loc = k1_k2_pH_tot ELSE k1_k2_pH_tot_loc = .true. END IF !--------------------------------------------------------------------------- ! convert tracer units to per mass !--------------------------------------------------------------------------- DO i = 1,nx_block IF (mask(i)) THEN dic(i,iblock) = dic_in(i) * vol_to_mass ta(i,iblock) = ta_in(i) * vol_to_mass pt(i,iblock) = max(pt_in(i),c0) * vol_to_mass sit(i,iblock) = max(sit_in(i),c0) * vol_to_mass x1(i) = c10 ** (-phhi(i)) x2(i) = c10 ** (-phlo(i)) END IF ! if mask END DO ! i loop CALL comp_co3_coeffs(iblock, k, mask, t, s, k0, k1, k2, ff, k1_k2_pH_tot_loc) !--------------------------------------------------------------------------- ! 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. !--------------------------------------------------------------------------- CALL drtsafe_row(iblock, j, k, mask, k1, k2, x1, x2, xacc, htotal) END SUBROUTINE comp_htotal_kvals !***************************************************************************** SUBROUTINE comp_co3_coeffs(iblock, k, mask, t, s, k0, k1, k2, ff, k1_k2_pH_tot) !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: iblock INTEGER(KIND=int_kind), INTENT(IN) :: k LOGICAL(KIND=log_kind), DIMENSION(nx_block), INTENT(IN) :: mask REAL(KIND=r8), DIMENSION(nx_block), INTENT(IN) :: & t, & ! temperature (degrees C) s ! salinity (PSU) LOGICAL(KIND=log_kind), INTENT(IN) :: k1_k2_pH_tot !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(OUT) :: & k0,k1,k2, & ! equilibrium constants for CO2 species ff ! fugacity of CO2 !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- INTEGER(KIND=int_kind) :: i REAL(KIND=r8) :: & press_bar ! pressure at level k [bars] REAL(KIND=r8), DIMENSION(nx_block) :: & tk, & ! temperature (K) is, & ! ionic strength scl, & ! chlorinity tk100, tk1002, invtk, dlogtk, is2, sqrtis, & s2, sqrts, s15, invRtk, arg, & deltaV,Kappa,lnKfac,Kfac, & ! pressure correction terms log_1_m_1p005em3_s, & log_1_p_tot_sulfate_div_ks !--------------------------------------------------------------------------- press_bar = ref_pressure(k) !--------------------------------------------------------------------------- ! 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 !--------------------------------------------------------------------------- tk = T0_Kelvin + t tk100 = tk * 1e-2_r8 tk1002 = tk100 * tk100 invtk = c1 / tk #ifdef CCSMCOUPLED CALL shr_vmath_log(tk, dlogtk, nx_block) #else dlogtk = LOG(tk) #endif invRtk = (c1 / 83.1451_r8) * invtk is = 19.924_r8 * s / (c1000 - 1.005_r8 * s) is2 = is * is #ifdef CCSMCOUPLED CALL shr_vmath_sqrt(is, sqrtis, nx_block) CALL shr_vmath_sqrt(s, sqrts, nx_block) #else sqrtis = SQRT(is) sqrts = SQRT(s) #endif s2 = s * s scl = s / 1.80655_r8 arg = c1 - 0.001005_r8 * s #ifdef CCSMCOUPLED CALL shr_vmath_log(arg, log_1_m_1p005em3_s, nx_block) #else log_1_m_1p005em3_s = LOG(arg) #endif !--------------------------------------------------------------------------- ! f = k0(1-pH2O)*correction term for non-ideality ! Weiss & Price (1980, Mar. Chem., 8, 347-359; ! Eq 13 with table 6 values) !--------------------------------------------------------------------------- arg = -162.8301_r8 + 218.2968_r8 / tk100 + & 90.9241_r8 * (dlogtk + LOG(1e-2_r8)) - 1.47696_r8 * tk1002 + & s * (.025695_r8 - .025225_r8 * tk100 + 0.0049867_r8 * tk1002) #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, ff, nx_block) #else ff = EXP(arg) #endif !--------------------------------------------------------------------------- ! K0 from Weiss 1974 !--------------------------------------------------------------------------- arg = 93.4517_r8 / tk100 - 60.2409_r8 + 23.3585_r8 * (dlogtk + LOG(1e-2_r8)) + & s * (.023517_r8 - 0.023656_r8 * tk100 + 0.0047036_r8 * tk1002) #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, k0, nx_block) #else k0 = EXP(arg) #endif !--------------------------------------------------------------------------- ! k1 = [H][HCO3]/[H2CO3] ! k2 = [H][CO3]/[HCO3] ! if k1_k2_pH_tot == .true., then use ! Lueker, Dickson, Keeling (2000) using Mehrbach et al. data on total scale ! otherwise, use ! Millero p.664 (1995) using Mehrbach et al. data on seawater scale ! this is only present to be consistent w/ OCMIP2 code ! it should not be used for new runs ! the only reason to use it is to be compatible with prior ! long spun up runs that had used it ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS !--------------------------------------------------------------------------- IF (k1_k2_pH_tot) THEN ! total pH scale arg = 3633.86_r8 * invtk - 61.2172_r8 + & 9.67770_r8 * dlogtk - 0.011555_r8 * s + & 0.0001152_r8 * s2 ELSE ! seawater pH scale, see comment above arg = 3670.7_r8 * invtk - 62.008_r8 + & 9.7944_r8 * dlogtk - 0.0118_r8 * s + & 0.000116_r8 * s2 END IF arg = -LOG(c10) * arg #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, k1, nx_block) #else k1 = EXP(arg) #endif IF (k > 1) THEN deltaV = -25.5_r8 + 0.1271_r8 * t Kappa = (-3.08_r8 + 0.0877_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif k1 = k1 * Kfac END IF IF (k1_k2_pH_tot) THEN ! total pH scale arg = 471.78_r8 * invtk + 25.9290_r8 - & 3.16967_r8 * dlogtk - 0.01781_r8 * s + 0.0001122_r8 * s2 ELSE ! seawater pH scale, see comment above arg = 1394.7_r8 * invtk + 4.777_r8 - & 0.0184_r8 * s + 0.000118_r8 * s2 END IF arg = -LOG(c10) * arg #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, k2, nx_block) #else k2 = EXP(arg) #endif IF (k > 1) THEN deltaV = -15.82_r8 - 0.0219_r8 * t Kappa = (1.13_r8 - 0.1475_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif k2 = k2 * Kfac END IF !--------------------------------------------------------------------------- ! kb = [H][BO2]/[HBO2] ! Millero p.669 (1995) using data from Dickson (1990) ! CO2SYS states that this in on total pH scale ! pressure correction from Millero 1979, p. 1657 ! omitting salinity contribution !--------------------------------------------------------------------------- arg = (-8966.90_r8 - 2890.53_r8 * sqrts - & 77.942_r8 * s + 1.728_r8 * s * sqrts - & 0.0996_r8 * s2) * invtk + & (148.0248_r8 + 137.1942_r8 * sqrts + 1.62142_r8 * s) + & (-24.4344_r8 - 25.085_r8 * sqrts - 0.2474_r8 * s) * dlogtk + & 0.053105_r8 * sqrts * tk #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, kb(:,iblock), nx_block) #else kb(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -29.48_r8 + (0.1622_r8 - 0.002608_r8 * t) * t Kappa = -2.84_r8 * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif kb(:,iblock) = kb(:,iblock) * Kfac END IF !--------------------------------------------------------------------------- ! k1p = [H][H2PO4]/[H3PO4] ! DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS !--------------------------------------------------------------------------- arg = -4576.752_r8 * invtk + 115.525_r8 - & 18.453_r8 * dlogtk + & (-106.736_r8 * invtk + 0.69171_r8) * sqrts + & (-0.65643_r8 * invtk - 0.01844_r8) * s #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, k1p(:,iblock), nx_block) #else k1p(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -14.51_r8 + (0.1211_r8 - 0.000321_r8 * t) * t Kappa = (-2.67_r8 + 0.0427_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif k1p(:,iblock) = k1p(:,iblock) * Kfac END IF !--------------------------------------------------------------------------- ! k2p = [H][HPO4]/[H2PO4] ! DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)) ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS !--------------------------------------------------------------------------- arg = -8814.715_r8 * invtk + 172.0883_r8 - & 27.927_r8 * dlogtk + & (-160.340_r8 * invtk + 1.3566_r8) * sqrts + & (0.37335_r8 * invtk - 0.05778_r8) * s #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, k2p(:,iblock), nx_block) #else k2p(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -23.12_r8 + (0.1758_r8 - 0.002647_r8 * t) * t Kappa = (-5.15_r8 + 0.09_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif k2p(:,iblock) = k2p(:,iblock) * Kfac END IF !--------------------------------------------------------------------------- ! k3p = [H][PO4]/[HPO4] ! DOE(1994) eq 7.2.26 with footnote using data from Millero (1974) ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS !--------------------------------------------------------------------------- arg = -3070.75_r8 * invtk - 18.141_r8 + & (17.27039_r8 * invtk + 2.81197_r8) * sqrts + & (-44.99486_r8 * invtk - 0.09984_r8) * s #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, k3p(:,iblock), nx_block) #else k3p(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -26.57_r8 + (0.202_r8 - 0.003042_r8 * t) * t Kappa = (-4.08_r8 + 0.0714_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif k3p(:,iblock) = k3p(:,iblock) * Kfac END IF !--------------------------------------------------------------------------- ! ksi = [H][SiO(OH)3]/[Si(OH)4] ! Millero p.671 (1995) using data from Yao and Millero (1995) ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS ! apply boric acid values !--------------------------------------------------------------------------- arg = -8904.2_r8 * invtk + 117.385_r8 - & 19.334_r8 * dlogtk + & (-458.79_r8 * invtk + 3.5913_r8) * sqrtis + & (188.74_r8 * invtk - 1.5998_r8) * is + & (-12.1652_r8 * invtk + 0.07871_r8) * is2 + & log_1_m_1p005em3_s #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, ksi(:,iblock), nx_block) #else ksi(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -29.48_r8 + (0.1622_r8 - 0.002608_r8 * t) * t Kappa = -2.84_r8 * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif ksi(:,iblock) = ksi(:,iblock) * Kfac END IF !--------------------------------------------------------------------------- ! kw = [H][OH] ! Millero p.670 (1995) using composite data ! following DOE Handbook, 0.015 substracted from constant to ! approximately convert from SWS pH scale to total pH scale ! pressure correction from Millero 1983 ! note that deltaV coeffs in Millero 1995 are those actually ! freshwater deltaV coeffs from Millero 1983 !--------------------------------------------------------------------------- arg = -13847.26_r8 * invtk + 148.9652_r8 - 23.6521_r8 * dlogtk + & (118.67_r8 * invtk - 5.977_r8 + 1.0495_r8 * dlogtk) * sqrts - & 0.01615_r8 * s #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, kw(:,iblock), nx_block) #else kw = EXP(arg) #endif IF (k > 1) THEN deltaV = -20.02_r8 + (0.1119_r8 - 0.001409_r8 * t) * t Kappa = (-5.13_r8 + 0.0794_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif kw(:,iblock) = kw(:,iblock) * Kfac END IF !--------------------------------------------------------------------------- ! ks = [H][SO4]/[HSO4], free pH scale ! Dickson (1990, J. chem. Thermodynamics 22, 113) ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS !--------------------------------------------------------------------------- arg = -4276.1_r8 * invtk + 141.328_r8 - 23.093_r8 * dlogtk + & (-13856.0_r8 * invtk + 324.57_r8 - 47.986_r8 * dlogtk) * sqrtis + & (35474.0_r8 * invtk - 771.54_r8 + 114.723_r8 * dlogtk) * is - & 2698.0_r8 * invtk * is * sqrtis + & 1776.0_r8 * invtk * is2 + & log_1_m_1p005em3_s #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, ks(:,iblock), nx_block) #else ks(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -18.03_r8 + (0.0466_r8 + 0.000316_r8 * t) * t Kappa = (-4.53_r8 + 0.09_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif ks(:,iblock) = ks(:,iblock) * Kfac END IF !--------------------------------------------------------------------- ! kf = [H][F]/[HF] ! Dickson and Riley (1979) -- change pH scale to total ! pressure correction from Millero 1995, p. 675 ! w/ typo corrections from CO2SYS !--------------------------------------------------------------------- arg = c1 + (0.1400_r8 / 96.062_r8) * (scl) / ks(:,iblock) #ifdef CCSMCOUPLED CALL shr_vmath_log(arg, log_1_p_tot_sulfate_div_ks, nx_block) #else log_1_p_tot_sulfate_div_ks = LOG(arg) #endif arg = 1590.2_r8 * invtk - 12.641_r8 + 1.525_r8 * sqrtis + & log_1_m_1p005em3_s + log_1_p_tot_sulfate_div_ks #ifdef CCSMCOUPLED CALL shr_vmath_exp(arg, kf(:,iblock), nx_block) #else kf(:,iblock) = EXP(arg) #endif IF (k > 1) THEN deltaV = -9.78_r8 - (0.009_r8 + 0.000942_r8 * t) * t Kappa = (-3.91_r8 + 0.054_r8 * t) * p001 lnKfac = (-deltaV + p5 * Kappa * press_bar) * press_bar * invRtk #ifdef CCSMCOUPLED CALL shr_vmath_exp(lnKfac, Kfac, nx_block) #else Kfac = EXP(lnKfac) #endif kf(:,iblock) = kf(:,iblock) * Kfac END IF !--------------------------------------------------------------------- ! Calculate concentrations for borate, sulfate, and fluoride ! bt : Uppstrom (1974) ! st : Morris & Riley (1966) ! ft : Riley (1965) !--------------------------------------------------------------------- bt(:,iblock) = 0.000232_r8 / 10.811_r8 * scl st(:,iblock) = 0.14_r8 / 96.062_r8 * scl ft(:,iblock) = 0.000067_r8 / 18.9984_r8 * scl ! DO i = 1,nx_block ! IF (mask(i)) THEN ! END IF ! if mask ! END DO ! i loop END SUBROUTINE comp_co3_coeffs !***************************************************************************** SUBROUTINE talk_row(iblock, mask, k1, k2, x, fn, df) !--------------------------------------------------------------------------- ! This routine computes 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. !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: iblock LOGICAL(KIND=log_kind), DIMENSION(nx_block), INTENT(IN) :: mask REAL(KIND=r8), DIMENSION(nx_block), INTENT(IN) :: k1, k2, x !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(OUT) :: fn, df !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- INTEGER(KIND=int_kind) :: i REAL(KIND=r8) :: & x1, x1_r, x2, x2_r, x3, k12, k12p, k123p, & a, a_r, a2_r, da, b, b_r, b2_r, db, c, c_r, & kb_p_x1_r, ksi_p_x1_r, c1_p_c_ks_x1_r_r, c1_p_kf_x1_r_r !--------------------------------------------------------------------------- DO i = 1,nx_block IF (mask(i)) THEN x1 = x(i) x1_r = c1 / x1 x2 = x1 * x1 x2_r = x1_r * x1_r x3 = x2 * x1 k12 = k1(i) * k2(i) k12p = k1p(i,iblock) * k2p(i,iblock) k123p = k12p * k3p(i,iblock) a = x3 + k1p(i,iblock) * x2 + k12p * x1 + k123p a_r = c1 / a a2_r = a_r * a_r da = c3 * x2 + c2 * k1p(i,iblock) * x1 + k12p b = x2 + k1(i) * x1 + k12 b_r = c1 / b b2_r = b_r * b_r db = c2 * x1 + k1(i) c = c1 + st(i,iblock) / ks(i,iblock) c_r = c1 / c kb_p_x1_r = c1 / (kb(i,iblock) + x1) ksi_p_x1_r = c1 / (ksi(i,iblock) + x1) c1_p_c_ks_x1_r_r = c1 / (c1 + c * ks(i,iblock) * x1_r) c1_p_kf_x1_r_r = c1 / (c1 + kf(i,iblock) * x1_r) !--------------------------------------------------------------------- ! fn = hco3+co3+borate+oh+hpo4+2*po4+silicate-hfree-hso4-hf-h3po4-ta !--------------------------------------------------------------------- fn(i) = k1(i) * dic(i,iblock) * x1 * b_r & + c2 * dic(i,iblock) * k12 * b_r & + bt(i,iblock) * kb(i,iblock) * kb_p_x1_r & + kw(i,iblock) * x1_r & + pt(i,iblock) * k12p * x1 * a_r & + c2 * pt(i,iblock) * k123p * a_r & + sit(i,iblock) * ksi(i,iblock) * ksi_p_x1_r & - x1 * c_r & - st(i,iblock) * c1_p_c_ks_x1_r_r & - ft(i,iblock) * c1_p_kf_x1_r_r & - pt(i,iblock) * x3 * a_r & - ta(i,iblock) !--------------------------------------------------------------------- ! df = d(fn)/dx !--------------------------------------------------------------------- df(i) = k1(i) * dic(i,iblock) * (b - x1 * db) * b2_r & - c2 * dic(i,iblock) * k12 * db * b2_r & - bt(i,iblock) * kb(i,iblock) * kb_p_x1_r * kb_p_x1_r & - kw(i,iblock) * x2_r & + (pt(i,iblock) * k12p * (a - x1 * da)) * a2_r & - c2 * pt(i,iblock) * k123p * da * a2_r & - sit(i,iblock) * ksi(i,iblock) * ksi_p_x1_r * ksi_p_x1_r & - c1 * c_r & - st(i,iblock) * c1_p_c_ks_x1_r_r * c1_p_c_ks_x1_r_r * (c * ks(i,iblock) * x2_r) & - ft(i,iblock) * c1_p_kf_x1_r_r * c1_p_kf_x1_r_r * kf(i,iblock) * x2_r & - pt(i,iblock) * x2 * (c3 * a - x1 * da) * a2_r END IF ! if mask END DO ! i loop END SUBROUTINE talk_row !***************************************************************************** SUBROUTINE drtsafe_row(iblock, j, k, mask_in, k1, k2, x1, x2, xacc, soln) !--------------------------------------------------------------------------- ! Vectorized version of drtsafe, which was a modified version of ! Numerical Recipes algorithm. ! Keith Lindsay, Oct 1999 ! ! Algorithm comment : ! Iteration from Newtons method is used unless it leaves ! bracketing interval or the dx is > 0.5 the previous dx. ! In that case, bisection method is used. !--------------------------------------------------------------------------- #ifdef CCSMCOUPLED USE shr_sys_mod, ONLY : shr_sys_abort #endif !--------------------------------------------------------------------------- ! input arguments !--------------------------------------------------------------------------- INTEGER(KIND=int_kind), INTENT(IN) :: iblock, j, k LOGICAL(KIND=log_kind), DIMENSION(nx_block), INTENT(IN) :: mask_in REAL(KIND=r8), DIMENSION(nx_block), INTENT(IN) :: k1, k2 REAL(KIND=r8), INTENT(IN) :: xacc !--------------------------------------------------------------------------- ! input/output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(INOUT) :: x1, x2 !--------------------------------------------------------------------------- ! output arguments !--------------------------------------------------------------------------- REAL(KIND=r8), DIMENSION(nx_block), INTENT(OUT) :: soln !--------------------------------------------------------------------------- ! local variable declarations !--------------------------------------------------------------------------- LOGICAL(KIND=log_kind) :: leave_bracket, dx_decrease LOGICAL(KIND=log_kind), DIMENSION(nx_block) :: mask INTEGER(KIND=int_kind) :: i, it REAL(KIND=r8) :: temp REAL(KIND=r8), DIMENSION(nx_block) :: xlo, xhi, flo, fhi, f, df, dxold, dx TYPE (block) :: this_block !--------------------------------------------------------------------------- ! bracket root at each location and set up first iteration !--------------------------------------------------------------------------- mask = mask_in it = 0 DO CALL talk_row(iblock, mask, k1, k2, x1, flo, df) CALL talk_row(iblock, mask, k1, k2, x2, fhi, df) WHERE ( mask ) mask = (flo > c0 .AND. fhi > c0) .OR. & (flo < c0 .AND. fhi < c0) END WHERE IF (.NOT. ANY(mask)) EXIT it = it + 1 DO i = 1,nx_block IF (mask(i)) THEN this_block = get_block(blocks_clinic(iblock), iblock) WRITE(stdout,*) '(co2calc.F90:drtsafe_row) ', & 'i_glob = ', this_block%i_glob(i), & ', j_glob = ', this_block%j_glob(j), ', k = ', k, & ', nsteps_run = ', nsteps_run, ', it = ', it WRITE(stdout,*) '(co2calc.F90:drtsafe_row) ', & ' x1,f = ', x1(i), flo(i) WRITE(stdout,*) '(co2calc.F90:drtsafe_row) ', & ' x2,f = ', x2(i), fhi(i) END IF END DO IF (it > max_bracket_grow_it) THEN CALL shr_sys_abort('bounding bracket for pH solution not found') END IF WHERE ( mask ) dx = sqrt(x2 / x1) x2 = x2 * dx x1 = x1 / dx END WHERE END DO mask = mask_in DO i = 1,nx_block IF (mask(i)) THEN IF (flo(i) .LT. c0) THEN xlo(i) = x1(i) xhi(i) = x2(i) ELSE xlo(i) = x2(i) xhi(i) = x1(i) temp = flo(i) flo(i) = fhi(i) fhi(i) = temp END IF soln(i) = p5 * (xlo(i) + xhi(i)) dxold(i) = ABS(xlo(i) - xhi(i)) dx(i) = dxold(i) END IF END DO CALL talk_row(iblock, mask, k1, k2, soln, f, df) !--------------------------------------------------------------------------- ! perform iterations, zeroing mask when a location has converged !--------------------------------------------------------------------------- DO it = 1,maxit DO i = 1,nx_block IF (mask(i)) THEN leave_bracket = ((soln(i) - xhi(i)) * df(i) - f(i)) * & ((soln(i) - xlo(i)) * df(i) - f(i)) .GE. 0 dx_decrease = ABS(c2 * f(i)) .LE. ABS(dxold(i) * df(i)) IF (leave_bracket .OR. .NOT. dx_decrease) THEN dxold(i) = dx(i) dx(i) = p5 * (xhi(i) - xlo(i)) soln(i) = xlo(i) + dx(i) IF (xlo(i) .EQ. soln(i)) mask(i) = .FALSE. ELSE dxold(i) = dx(i) dx(i) = -f(i) / df(i) temp = soln(i) soln(i) = soln(i) + dx(i) IF (temp .EQ. soln(i)) mask(i) = .FALSE. END IF IF (ABS(dx(i)) .LT. xacc) mask(i) = .FALSE. END IF END DO IF (.NOT. ANY(mask)) RETURN CALL talk_row(iblock, mask, k1, k2, soln, f, df) DO i = 1,nx_block IF (mask(i)) THEN IF (f(i) .LT. c0) THEN xlo(i) = soln(i) flo(i) = f(i) ELSE xhi(i) = soln(i) fhi(i) = f(i) END IF END IF END DO END DO ! iteration loop #ifdef CCSMCOUPLED CALL shr_sys_abort('lack of convergence in drtsafe_row') #endif END SUBROUTINE drtsafe_row !***************************************************************************** END MODULE co2calc