k12 = k1local*k2local

            b = x2 + k1local*x + k12

            db = 2.0*x + k1local

            fn = k1local*x*diclocal/b +

            df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -

     I                       k1local,k2local,

              k12 = k1local*k2local

              b = x2 + k1local*x + k12

              db = 2.0*x + k1local

              fn = k1local*x*diclocal/b +

        co2star=diclocal*htotal2/(htotal2 + k1local*htotal

     &            + k1local*k2local)

     I                       k1local,k2local,

        _RL  k1local, k2local

        stuff = (1.0-gamm)*(1.0-gamm)*k1local*k1local

     &          - 4.0*k1local*k2local*(1.0-2.0*gamm)

        hnew  = 0.5*( (gamm-1.0)*k1local + sqrt(stuff) )

     &   (1.0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))

        co3local = k1local*k2local*diclocal /

     &         (hnew*hnew + k1local*hnew + k1local*k2local)

     I                       k1local,k2local,

        _RL  k1local, k2local

        stuff = (1.0 _d 0-gamm)*(1.0 _d 0-gamm)*k1local*k1local

     &          - 4.0 _d 0*k1local*k2local*(1.0 _d 0-2.0 _d 0*gamm)

        hnew  = 0.5 _d 0*( (gamm-1.0 _d 0)*k1local + sqrt(stuff) )

     &   (1.0 _d 0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))

        _RL  k1local, k2local