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

     &          ( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )

     &    ' tsr,tta,q,qsr  =', tsr, tta, q, qsr

       zL=xkar*gravity*zu*( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )

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

       Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr)

      ef=-lath*rhoa*usr*qsr          !latent W/m2

      evap=-rhoa*usr*qsr

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

      _RL CD,usr,tsr,qsr