anchor.diff:    !> 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>
  
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:         tempis(i) = 0.99975*tempis68 + 0.0002
  !       Note: parts (a) and (c) above are tiny corrections;
  !             part  (b) is a big correction for deep waters (but zero at surface)
       ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN
  !       When optT = 'Tinsitu', tempis is input & output (no tempot needed)
          tempis(i) = temp(i)
          tempis68 = (temp(i) - 0.0002) / 0.99975
  !       tempis(i) = 0.99975*tempis68 + 0.0002
       ELSE
          PRINT *,"optT must be either 'Tpot' or 'Tinsitu'"
          PRINT *,"you specified optT =", trim(optT) 
          STOP
       ENDIF
  
  