tsr=-(Dt)*xkar/(LOG(zt/zot10)-psit(zt/L10))

      qsr=-(Dq)*xkar/(LOG(zq/zot10)-psit(zq/L10))

       pst=psit(zt/L)

       tsr=-(Dt)*xkar/(LOG(zt/zot)-psit(zt/L))

       qsr=-(Dq)*xkar/(LOG(zq/zoq)-psit(zq/L))

      q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L))

     &                *xkar/(LOG(zt/zot10)-psit(zt/L10))

      dEvapdTs = rhoa*usr*( xkar/(LOG(zq/zoq)-psit(zq/L)) )

      _RL FUNCTION psit(zL)

       psit = (1.-f)*psik+f*psic

       psit = -( (1.+2.*zL/3.)**1.5

      _RL psiu,psit,zot10,Ct10,CC,Ribu