anchor.diff:  
     !     INPUT variables:
     !     ================
  +  !     Patm    = atmospheric pressure [atm]
     !     depth   = depth [m]     (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure)
     !             = pressure [db] (with optP='db')
     !     lat     = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m')
  
positive.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)
    !> 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>
  
negative.code:         invtk=1.0d0/tk
          dlogtk=LOG(tk)
  
  !       Atmospheric pressure
          Patmd = DBLE(Patm(i))
  
  !       Pressure effect (prb is in bars)
          prb = DBLE(p(i)) / 10.0d0
  
  