1 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   2 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   3 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   4 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   5 | c
   6 | c.......................................................................
   7 | c
   8 | c  k-cumulative for surface-surface exchanges
   9 | c
  10 | c  in :  * $\phi$                ---> phcdss
  11 | c        * $\overline k$         ---> cbcdss
  12 | c        * $l$                   ---> alcdss
  13 | c        * $k$                   ---> ccdss
  14 | c
  15 | c  out : * $g^{ss} (k;l)$        ---> cdcdss
  16 | c
  17 | c.......................................................................
  18 |       subroutine cdss(phcdss,cbcdss,alcdss,ccdss,cdcdss)
  19 | c.......................................................................
  20 |       implicit double precision (a-h,o-z)
  21 | c.......................................................................
  22 |       cdss1=dsqrt(1.d0+2.d0/phcdss*cbcdss*alcdss)
  23 |       cdss2=phcdss*cdss1
  24 |       cdss3=cbcdss/cdss1
  25 |       call gicd(cdss2,cdss3,ccdss,cdcdss)
  26 | c.......................................................................
  27 |       return
  28 |       end