anchor.diff:  !!    total alkalinity (ALK), dissolved inorganic carbon (DIC),
   !!    silica and phosphate concentrations (all 1-D arrays)
   SUBROUTINE vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis,  &
  -                temp, sal, alk, dic, sil, phos, depth, lat, N,                          &
  +                temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,                          &
                   optCON, optT, optP, optB, optK1K2, optKf                                  )
   
     !   Purpose:
  
positive.code:         hco3(i) = REAL(cb * drho)
          co3(i)  = REAL(cc * drho)
  
  !       Determine CO2 pressure and fugacity (in microatm)
  !       NOTE: equation below for pCO2 requires CO2* in mol/kg
          dfCO2 = cu   * 1.e6_r8/K0
          B=(-1636.75d0+12.0408d0*tk-0.0327957d0*(tk*tk)+0.0000316528d0* &
               (tk*tk*tk))*1e-6
  
          dpCO2= dfCO2 / EXP(100000.0d0 &
               *(B+ 2.0d0*(57.7d0 - 0.118d0*tk)*1e-6_r8) &
               / (8.314d0*tk) &
               )
  
          fCO2(i) = REAL(dfCO2)
          pCO2(i) = REAL(dpCO2)
  
negative.code:   !> phosphate concentration in <b>[mol/m^3]</b> (when optCON = 'mol/m3') OR in <b>[mol/kg]</b> (when optCON = 'mol/kg')
    REAL(kind=r4), INTENT(in), DIMENSION(N) :: phos
  !f2py optional , depend(sal) :: n=len(sal)
    !> atmospheric pressure <b>[atm]</b>
    REAL(kind=r4), INTENT(in), DIMENSION(N) :: Patm
    !> depth in \b meters (when optP='m') or \b decibars (when optP='db')
    REAL(kind=r4), INTENT(in),    DIMENSION(N) :: depth
    !> latitude <b>[degrees north]</b>
  