subroutine forcing(jt)
use types,only:rprec
use param,only:use_bldg,nx,ny,nz,DT,dns_bc,z_i,L_z,dx,dy,&
     jx_s,x_relax,face_avg,inflow
use sim_param,only:u,v,w,dpdx,dpdy,dpdz,RHSx,RHSy,RHSz,RHSx_f,RHSy_f,RHSz_f
use immersedbc
implicit none
real(kind=rprec)::force,RHS
integer,intent(in)::jt
integer::px,py,lx,ly,lz
integer :: jx,jy,jz,i,isize,jsize
real(kind=rprec),dimension(4)::y,y1,y2,y12
! this sets pressure grad's inside bldg's
if (use_bldg) then
      do i=1,n_bldg
         px=bldg_pts(1,i)
         py=bldg_pts(2,i)
         lx=bldg_pts(3,i)
         ly=bldg_pts(4,i)
         lz=bldg_pts(5,i)
         dpdx(px:px+lx,py:py+ly-1,1:1+lz)=0._rprec
         dpdy(px:px+lx-1,py:py+ly,1:1+lz)=0._rprec
         dpdz(px:px+lx-1,py:py+ly-1,1:1+lz)=0._rprec
      end do
   end if
! Finalize right hand side for explicit time advancement
   RHSx=RHSx-dpdx
   RHSy=RHSy-dpdy
   RHSz=RHSz-dpdz
! Start first step using new C as former C
! should be able to move this out of time loop
   if ((jt == 1) .and. (RHSx_f(1,1,1) == 0._rprec)) then 
      RHSx_f=RHSx
      RHSy_f=RHSy
      RHSz_f=RHSz
   end if

   force=z_i/L_z
! start calculation of body forces (fx,fy,fz)
! 'force' is the mean pressure gradient
   if (use_bldg) then
      do i=1,n_bldg
         px=bldg_pts(1,i)
         py=bldg_pts(2,i)
         lx=bldg_pts(3,i)
         ly=bldg_pts(4,i)
         lz=bldg_pts(5,i)
         do jz=1,1+lz
         do jy=py,py+ly-1
         do jx=px,px+lx
            RHS=1.5_rprec*RHSx(jx,jy,jz)-0.5_rprec*RHSx_f(jx,jy,jz)+force
            fx(jx,jy,jz) = (u_des(jx,jy,jz)-u(jx,jy,jz))/dt - RHS
!            RHS=1.5_rprec*RHSy(jx,jy,jz)-0.5_rprec*RHSy_f(jx,jy,jz)
!            fy(jx,jy,jz) = (v_des(jx,jy,jz)-v(jx,jy,jz))/dt - RHS
!            RHS=1.5_rprec*RHSz(jx,jy,jz)-0.5_rprec*RHSz_f(jx,jy,jz)
!            fz(jx,jy,jz) = (w_des(jx,jy,jz)-w(jx,jy,jz))/dt - RHS
         end do
         end do
         do jy=py,py+ly
         do jx=px,px+lx-1
            RHS=1.5_rprec*RHSy(jx,jy,jz)-0.5_rprec*RHSy_f(jx,jy,jz)
            fy(jx,jy,jz) = (v_des(jx,jy,jz)-v(jx,jy,jz))/dt - RHS
         end do
         end do
         do jy=py,py+ly
         do jx=px,px+lx
            RHS=1.5_rprec*RHSz(jx,jy,jz)-0.5_rprec*RHSz_f(jx,jy,jz)
            fz(jx,jy,jz) = (w_des(jx,jy,jz)-w(jx,jy,jz))/dt - RHS
         end do
         end do
         end do
      end do
! end calculation of forces
   endif
! this is the mean pressure gradient forcing
   do jz=1,nz-1
   do jy=1,ny
   do jx=1,nx
      RHS=1.5_rprec*RHSx(jx,jy,jz)-0.5_rprec*RHSx_f(jx,jy,jz)+force
      u(jx,jy,jz) = (u(jx,jy,jz) + dt*(RHS+fx(jx,jy,jz)))
      RHS=1.5_rprec*RHSy(jx,jy,jz)-0.5_rprec*RHSy_f(jx,jy,jz)
      v(jx,jy,jz) = (v(jx,jy,jz) + dt*(RHS+fy(jx,jy,jz)))
   end do
   end do
   end do
!TS ADD UNIFORM FLOW
   if(inflow)then
   do jy=0,x_relax-1
      jx=jx_s-jy
      u(jx,1:ny,1:nz-1)=face_avg+&
      (u(jx,1:ny,1:nz-1)-face_avg)*real(jy,kind=rprec)/real(x_relax,kind=rprec)
      v(jx,1:ny,1:nz-1)=(v(jx,1:ny,1:nz-1)-0._rprec)*&
           real(jy,kind=rprec)/real(x_relax,kind=rprec)
      w(jx,1:ny,1:nz)=(w(jx,1:ny,1:nz)-0._rprec)*&
           real(jy,kind=rprec)/real(x_relax,kind=rprec)
   end do
   endif

! no-stress top
   u(:,:,nz)=u(:,:,nz-1)
! no-stress top
   v(:,:,nz)=v(:,:,nz-1)
   if (dns_bc) then
      do jz=2,nz-1
      do jy=1,ny
      do jx=1,nx
         RHS=1.5_rprec*RHSz(jx,jy,jz)-0.5_rprec*RHSz_f(jx,jy,jz)
         w(jx,jy,jz) = (w(jx,jy,jz) + dt*(RHS+fz(jx,jy,jz)))
      end do
      end do
      end do
      w(:,:,1)=0._rprec
      w(:,:,nz)=0._rprec
   else
      call stepbl_w(w, RHSz, RHSz_f, fz)              
   end if
 end subroutine forcing
