_RL  dayfrac, yday, delta

        delta = (0.006918 - (0.399912*cos(yday))      !cosine zenith angle

          sun1 = -sin(delta)/cos(delta) * sin(lat)/cos(lat)

          cosz = ( sin(delta)*sin(lat)+              !average zenith angle

     &            (cos(delta)*cos(lat)*sin(dayhrs)/dayhrs) )

      _KPP_RL delta

            delta = (hbl(i) + zgrid(ki)) / (zgrid(ki) - zgrid(ki+1))

               dstar         = (1.- delta)**2 * dkm1(i,md)

     &                       + delta**2 * dkmp5

               blmc(i,ki,md) = (1.- delta)*diffus(i,ki,md)

     &                       + delta*dstar