anchor.diff:     END DO
   
      call vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis,         &  ! OUTPUT
  -             temp, sal, alk, dic, sil, phos, depth, lat, N,                                   &  ! INPUT
  +             temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,                             &  ! INPUT
                optCON='mol/kg', optT='Tinsitu', optP='db', optB='l10', optK1K2='m10', optKf='dg')  ! INPUT OPTIONS
   !  Print out results (typical for data: concentration units differ)
      write(*,50)
  
positive.code:    REAL(kind=r4), DIMENSION(100) :: ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis
  !  "vars" Input variables
     INTEGER :: N
     REAL(kind=r4), DIMENSION(100) :: temp, sal, alk, dic, sil, phos, depth, lat
  !  "vars" Input options
     CHARACTER(10) :: optCON, optT, optP, optB, optKf, optK1K2
  
  
negative.code: !> Call mocsy's main subroutine to compute carbonate system variables: pH, pCO2, fCO2, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R
  !> FROM temperature, salinity, total alkalinity, dissolved inorganic carbon, silica, phosphate, depth (or pressure) (1-D arrays)
     call vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis,  &  ! OUTPUT
               temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,                      &  ! INPUT
               optCON, optT, optP)
  !            optCON, optT, optP, optB, optK1K2, optKf)                                    ! INPUT OPTIONS
  !  Print out results
     write(*,50)
     write(*,100)
  