anchor.diff:    !       -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007)
     !               **** BUT Valid for  9 < T < 33°C and 10 < S < 40.
     !       -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994)
  +  !     -----------
  +  !     optGAS: choose in situ vs. potential fCO2 and pCO2
  +  !     ---------
  +  !       PRESSURE corrections for K0 and the fugacity coefficient (Cf) 
  +  !       -> 'Pzero'   = 'zero order' fCO2 and pCO2 (typical approach, which is flawed)
  +  !                      considers in situ T & only atm pressure (hydrostatic=0)
  +  !       -> 'Ppot'    = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
  +  !                      considers potential T & only atm pressure (hydrostatic press = 0)
  +  !       -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure)
  +  !                      considers in situ T & total pressure (atm + hydrostatic)
  +  !     ---------
   
     !     OUTPUT variables:
     !     =================
     !     ph   = pH on total scale
     !     pco2 = CO2 partial pressure (uatm)
     !     fco2 = CO2 fugacity (uatm)
  -  !     co2  = aqueous CO2 concentration [mol/m^3]
  -  !     hco3 = bicarbonate (HCO3-) concentration [mol/m^3]
  -  !     co3  = carbonate (CO3--) concentration [mol/m^3]
  +  !     co2  = aqueous CO2 concentration in [mol/kg] or [mol/m^3] depending on optCON
  +  !     hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] depending on optCON
  +  !     co3  = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] depending on optCON
     !     OmegaA = Omega for aragonite, i.e., the aragonite saturation state
     !     OmegaC = Omega for calcite, i.e., the   calcite saturation state
     !     BetaD = Revelle factor   dpCO2/pCO2 / dDIC/DIC
  
positive.code:   USE mp80
    USE mrho
    USE msw_temp
    USE mphsolvers
  
    IMPLICIT NONE
  
  
negative.code:   !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010) 
  !f2py character*3 optional, intent(in) :: optK1K2='l'
    CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2
    !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) 
    !! 'Ppot'    - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
    !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf
    !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean
  !f2py character*7 optional, intent(in) :: optGAS='Pinsitu'
    CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS
  
  ! Output variables:
    !> pH on the <b>total scale</b>
  