subroutine derivwall()
! calculates deriv at wall as suggested in Moeng '84'
use types,only:rprec
use param,only:vonk,dz,ld,nx,ny,nz
use sim_param,only:dudz,dvdz,u,v
use bottombc,only:zo,phi_m,psi_m
implicit none
integer::jx,jy, method
real(kind=rprec),dimension(nx,ny)::ustar,u_avg

!TS Remove the following line when obukhov.f is used
psi_m=0._rprec
phi_m=1._rprec

method = 1
if (method == 1) then
   u_avg(1:nx,1:ny)=sqrt(u(1:nx,1:ny,1)**2+v(1:nx,1:ny,1)**2)
!TS ADD for non-neutral case
   ustar=u_avg*vonK/(log((0.5_rprec*dz)/zo)-psi_m)
   do jy=1,ny
   do jx=1,nx
      !this is as in Moeng 84
      dudz(jx,jy,1)=ustar(jx,jy)/(0.5_rprec*dz*vonK)*u(jx,jy,1)/u_avg(jx,jy)&
!TS ADD for non-neutral case
           *phi_m(jx,jy)
! sc: i dont know if this is correct yet
      dvdz(jx,jy,1)=ustar(jx,jy)/(0.5_rprec*dz*vonK)*v(jx,jy,1)/u_avg(jx,jy)&
!TS ADD for non-neutral case
           *phi_m(jx,jy)
      dudz(jx,jy,1)=merge(0._rprec,dudz(jx,jy,1),u(jx,jy,1).eq.0._rprec)
      dvdz(jx,jy,1)=merge(0._rprec,dvdz(jx,jy,1),v(jx,jy,1).eq.0._rprec)
   end do
   end do
elseif (method == 2) then
   dudz(1:nx,:,1)=dudz(1:nx,:,2)*2._rprec
   dvdz(1:nx,:,1)=dvdz(1:nx,:,2)/2._rprec
end if
end subroutine derivwall
