anchor.diff:       sil(i)    = 50.  *1.0e-3   !Convert observed S. Ocean ave surf SiO2 (umol/L) to model units (mol/m3) 
        phos(i)   = 1.8  *1.0e-3   !Convert observed S. Ocean ave surf PO4  (umol/L) to model units (mol/m3)
        lat(i)    = -60.           !Latitude used to convert model depth to pressure (Saunders, 1981, JPO)
  -     depth(i) = real(i-1) * 1000. ! depth varies from 0 to 5000 m by 1000 m
  +     depth(i)  = real(i-1) * 1000. ! depth varies from 0 to 5000 m by 1000 m
  +     Patm(i)   = 1.0            !Atmospheric pressure (atm)
        N = i
      END DO
   
  
positive.code:    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
               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)
  
negative.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, Patm, depth, lat
  !  "vars" Input options
     CHARACTER(10) :: optCON, optT, optP, optB, optKf, optK1K2
  
  