anchor.diff:       ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN
   !       When optT = 'Tinsitu', tempis is input & output (no tempot needed)
           tempis = temp(i)
  +        tempis68 = (temp(i) - 0.0002) / 0.99975
  +        dtempot = sw_ptmp(DBLE(sal(i)), DBLE(tempis68), DBLE(p), 0.0d0)
        ELSE
           PRINT *,"optT must be either 'Tpot' or 'Tinsitu'"
           PRINT *,"you specified optT =", trim(optT) 
  
positive.code:         Kb(i)  = Kb(i)*total2SWS_0p
  
  !       Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw
  !       K1p(i) = K1p(i)*total2SWS_0p 
  !       K2p(i) = K2p(i)*total2SWS_0p
  !       K3p(i) = K3p(i)*total2SWS_0p
  !       Ksi(i) = Ksi(i)*total2SWS_0p
  !       Kw(i)  = Kw(i)*total2SWS_0p
  
  !       Keep other contants on same scale:
  !          - Ks (already on Free scale)
  !          - Kspc, Kspa (no scale)
  
  !       Perform actual pressure correction (on seawater scale)
          K1(i)   = K1(i)*EXP(lnkpok0(1))
  
negative.code:                      St, Ft, Bt,                              &
                       temp, sal, Patm,                         &
                       depth, lat, N,                           &
                       optT, optP, optB, optK1K2, optKf, optGAS)
  
    !   Purpose:
    !     Compute thermodynamic constants
  