tsr=-(Dt)*xkar/(LOG(zt/zot10)-psit(zt/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 )

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

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

      hf=-cpair*rhoa*usr*tsr         !sensible W/m2

        dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr

        dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr

        dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr

      _RL CD,usr,tsr,qsr