Skip to content

p4zche

Overview

The p4zche module implements seawater chemistry calculations for the PISCES (Pelagic Interaction Scheme for Carbon and Ecosystem Studies) biogeochemical model, following OCMIP (Ocean Carbon-Cycle Model Intercomparison Project) protocols. This module computes chemical equilibrium constants and solubilities essential for simulating ocean carbon cycle dynamics.

Module Information

  • Primary Standards: MOCSY (Marine Carbonate System) standards (since v3.6)

Dependencies

  • oce_trc: Shared ocean-tracer variables
  • trc: Passive tracer common variables
  • sms_pisces: PISCES source-sink variables
  • lib_mpp: MPP (Massively Parallel Processing) library
  • eosbn2: Equation of state (specifically neos parameter)

Public Interface

Public Subroutines

  • p4z_che: Main chemistry computation routine
  • p4z_che_alloc: Memory allocation
  • ahini_for_at: Initial pH guess for alkalinity solver
  • solve_at_general: Universal pH solver

Public Variables

Chemical Equilibrium Arrays

3D arrays (horizontal + vertical)

  • sio3eq: Silicon chemistry equilibrium
  • fekeq: Iron chemistry equilibrium
  • chemc: CO2 solubility coefficients (3 components)
  • chemo2: O2 solubility
  • akfe2ox: Fe(II) oxidation rate
  • salinprac: Practical salinity
  • tempis: In situ temperature

Carbonate System Constants

  • akb3: Boric acid dissociation constant
  • akw3: Water dissociation constant
  • akf3: Hydrofluoric acid dissociation constant
  • aks3: Bisulfate dissociation constant
  • ak1p3, ak2p3, ak3p3: Phosphoric acid dissociation constants
  • aksi3: Silicic acid dissociation constant

Concentration Variables

  • borat: Total borate concentration
  • fluorid: Total fluoride concentration
  • sulfat: Total sulfate concentration
  • fesol: Iron solubility (5 species)

Physical Constants

Atmospheric and Gas Constants

  • atcox = 0.20946: Atmospheric oxygen fraction (atm)
  • o2atm = 1/(1000 × 0.20946): Oxygen conversion factor
  • rgas = 83.14472: Universal gas constant (cm³·bar·K⁻¹·mol⁻¹)
  • oxyco = 1/22.4144: Ideal gas to moles conversion

Pressure Correction Coefficients

The module includes extensive pressure correction coefficients (devk series) for various equilibrium constants following Millero (1995):

  • devk10 through devk110: First-order temperature coefficients
  • devk20 through devk210: Second-order temperature coefficients
  • devk30 through devk310: Third-order temperature coefficients
  • devk40 through devk410: Linear pressure terms
  • devk50 through devk510: Quadratic pressure terms

Main Subroutines

p4z_che

Main routine computing seawater chemistry for all grid points.

Inputs: Kbb, Kmm - Time level indices

Process:

  1. Salinity Conversion: Converts absolute salinity to practical salinity when TEOS-10 equation of state is used: 

    salinprac = sal × 35.0 / 35.16504  (if neos = -1)

  2. In Situ Temperature Calculation: Approximates in situ temperature from potential temperature with pressure correction (accuracy <0.04°C): 

    tempis = θ - P × (a1 - a2 × P)

  3. Surface Layer CO2 Solubility: Computes using Weiss (1974) coefficients in mol/(L·µatm)

  4. Deep Ocean O2 Solubility: Uses Gordon & Garcia scaled temperature approach

  5. Deep Ocean Chemical Constants:

  6. Pressure calculation using Saunders (1980) formulation

  7. Dissociation constants for:

    • Carbonic acid (K1, K2) - Mehrbach (1973) refit by Millero (1995)
    • Boric acid (Kb)
    • Water (Kw) - Millero (1995)
    • Phosphoric acid (K1p, K2p, K3p)
    • Silicic acid (Ksi)
    • Sulfuric acid (Ks) - Dickson (1990)
    • Hydrofluoric acid (Kf) - Dickson & Riley (1979)

  8. pH Scale Conversions: Converts between Free, Seawater (SWS), and Total pH scales

  9. Pressure Corrections: Applies depth-dependent corrections to all constants

  10. Calcite Solubility: Computes Ksp using Mucci (1983) formulation

  11. Iron Chemistry:

    • Fe(II) oxidation kinetics
    • Solubility of 5 iron species using Liu & Millero (1999)

  12. Silicon Saturation: Computes Si(OH)4 equilibrium concentration

ahini_for_at

Provides initial H⁺ concentration guess for alkalinity solver using second-order polynomial approximation.

Inputs: Kbb - Time level index

Outputs: p_hini - Initial H⁺ concentration array

Returns:

  • 1×10⁻³: if alkalinity ≤ 0
  • 1×10⁻¹⁰: if alkalinity ≥ 2×DIC + borate
  • 1×10⁻⁷: if discriminant ≤ 0
  • Calculated root: otherwise

Algorithm: Solves cubic polynomial approximation around local minimum using discriminant test.

solve_at_general

Universal pH solver implementing Newton-Raphson method with adaptive bracketing.

Inputs:

  • p_hini: Initial H⁺ concentration
  • Kbb: Time level index

Outputs: zhi- Converged H⁺ concentration (-1 if failed)

Features:

  • Converges from any initial value
  • Automatic bracket determination using anw_infsup
  • Hybrid Newton-bisection method
  • Convergence criterion: |ΔH⁺/H⁺| < 10⁻⁴
  • Maximum iterations: 20

Chemical Species Handled:

  • Carbonate system (H\(_2\)CO\(_3\), HCO\(_3^{-}\), CO\(_3^{2-}\))
  • Borate (B(OH)\(_3\), B(OH)\(_4^-\))
  • Phosphate (H\(_3\)PO\(_4\), H\(_2\)PO\(^-_4\), HPO\(_4^{2-}\), PO\(_4^{3-}\))
  • Silicate (H\(_4\)SiO\(_4\), H\(_3\)SiO\(_4^-\))
  • Sulfate (HSO\(_4^-\), SO\(_4^{2-}\))
  • Fluoride (HF, F\(^-\))
  • Water (H\(_2\)O, OH\(^-\))

anw_infsup

Computes bounds for non-water-selfionization alkalinity contributions.

Outputs:

  • p_alknw_inf: Lower bound (infimum)
  • p_alknw_sup: Upper bound (supremum)

Numerical Parameters

  • pp_rdel_ah_target = 10⁻⁴: Target relative convergence tolerance
  • pp_ln10 = 2.302585...: Natural logarithm of 10
  • jp_maxniter_atgen = 20: Maximum Newton iterations

Memory Allocation

The p4z_che_alloc() function allocates all module arrays with proper error handling, returning the maximum error code from all allocation attempts.

Performance Considerations

  • Uses DO_3D macros for efficient domain decomposition
  • Implements timing instrumentation (ln_timing)
  • Optimized for vectorization with minimal conditional branching
  • Supports RK3 time-stepping with special handling for online coupling

pH Scale Conventions

The module uses Total pH scale as standard with conversion factors:

  • total2free = 1/(1 + SO4/Ks)
  • free2SWS = 1 + SO4/Ks + F/Kf
  • total2SWS = total2free × free2SWS

Limitations and Assumptions

  • In situ temperature approximation accurate to 0.04°C
  • Iron solubility valid only 5-50°C (Liu & Millero 1999)
  • Assumes constant stoichiometric ratios for nutrients
  • Pressure effects limited to Millero (1995) parameterizations

References

Key publications implemented:

  • Weiss (1974): Gas solubilities
  • Dickson (1990, with Riley 1979): Acid dissociation
  • Millero (1995): Comprehensive equilibria
  • Mucci (1983): Calcite solubility
  • Liu & Millero (1999): Iron chemistry

Notes for Users

Typical Usage Pattern

! 1. Allocate memory
ierr = p4z_che_alloc()

! 2. Compute chemistry at each time step
CALL p4z_che(Kbb, Kmm)

! 3. For pH solving (if needed)
CALL ahini_for_at(hini, Kbb)
CALL solve_at_general(hini, hi, Kbb)

Common Issues

  1. Convergence failures: pH solver returns -1 when convergence not achieved in 20 iterations
  2. Temperature range: Iron chemistry limited to 5-50°C range
  3. Salinity conversion: Automatic conversion between absolute and practical salinity based on EOS

Output Variables

All computed chemical constants and concentrations are stored in module-level arrays accessible throughout the PISCES model for biogeochemical calculations.