_RL    zu, zt, zq

     &       zu, zt, zq,

     I           deltaTtracer, zu, useFluxLimit,

     I           deltaTtracer, zu, useFluxLimit,

     I           deltaTtracer, zu, useFluxLimit,

      u10 = Du*LOG(10. _d 0/zo)/LOG(zu/zo)

      Cd = (xkar/LOG(zu/zo10))**2

      Ribcu=-zu/(zi*0.004 _d 0*xBeta**3)  ! Saturation or plateau Rib

      Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2)

      L10=zu/zetu                       ! MO length

      usr= Du*xkar/(LOG(zu/zo10)-psiu(zu/L10))

     &    ' L,zu,zL,zt    =', L, zu, zL, zt

     &    ' ln(zu/zo),psu,diff,zL*=', LOG(zu/zo), psu, LOG(zu/zo)-psu,

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

       L=zu/zL

       psu=psiu(zu/L)

       usr=Du*xkar/(LOG(zu/zo)-psiu(zu/L))

              hu=zu

      zu=10. _d 0

      zt=zu

      zq=zu

     &   dsolms,dsolmn,zu,zt,zq,

     &                 )*2. _d 0 *zu

          wvel(i,j,bi,bj) = localW(i,j)/2. _d 0 /zu

     I           deltaTtracer, zu,

     &        +1. _d 0/2. _d 0/zu*fVer(i,j)

      _RL deltaTtracer, zu