subroutine ic()
!...Log profile that is modified to flatten at z=z_i
use types,only:rprec
!only:ld,nx,ny,nz,z_i,wt_s,num_patch,zo1,zo2,dz,vonk,s_flag,u_star,theta_s1
use param
use sim_param,only:u,v,w
use bottombc,only:num_patch,zo1,zo2,theta_s1
implicit none
!      real*8 theta(ld,ny,nz),q(ld,ny,nz)
real(kind=rprec),dimension(nz)::ubar
real(kind=rprec)::rms, noise, arg, arg2
real(kind=rprec)::z,zo_avg,w_star,T_star,q_star,ran3
integer::jx,jy,jz,seed

w_star=(9.81_rprec/theta_s1*wt_s*z_i)**(1._rprec/3._rprec)
!      T_star=wt_s/w_star
!      q_star=T_star
seed=-80

! included by jan
if (num_patch.gt.1) then
   zo_avg=(zo1+zo2)/2._rprec
else
   zo_avg=zo1
end if

print *,'Modified Log Profile for IC'
do jz=1,nz
   z=(jz-.5_rprec)*dz*z_i
! IC in equilibrium with rough surface (rough dominates in effective zo)
   arg2=z/zo_avg
   arg=(1._rprec/vonk)*log(arg2)!-1./(2.*vonk*z_i*z_i)*z*z
   ubar(jz)=arg
! sc: I changed around some parenthesis here
   if ((s_flag).and.(z.gt.(.4_rprec*z_i))) ubar(jz)=ubar(jz-1)
end do
do jz=1,nz
   print *,'k, ubar:',jz,ubar(jz)
end do

rms = 3._rprec
do jz=1,nz
do jy=1,ny
do jx=1,nx
!...Ran3 returns uniform RV between 0 and 1. (sigma_rv=0.289)
!...Taking std dev of vel as 1 at all heights
   z=(real(jz)-.5_rprec)*dz*z_i
   if (z.le.z_i) then
      noise=rms/.289_rprec*(ran3(seed)-.5_rprec)
      u(jx,jy,jz)=noise*(1._rprec-z/z_i)*w_star/u_star+ubar(jz)
      noise=rms/.289_rprec*(ran3(seed)-0.5_rprec)
      v(jx,jy,jz)=noise*(1._rprec-z/z_i)*w_star/u_star !noise
      noise=rms/.289_rprec*(ran3(seed)-.5_rprec)
      w(jx,jy,jz)=noise*(1._rprec-z/z_i)*w_star/u_star
!	      theta(jx,jy,jz)=(theta_s1+10.*noise*(1-z/z_i)*T_star)/T_scale
!	      q(jx,jy,jz)=q_mix+50.*noise*(1-z/z_i)*q_star
   else
      noise=rms/.289_rprec*(ran3(seed)-.5_rprec)
      u(jx,jy,jz)=noise*w_star/u_star*.01_rprec+ubar(jz)
      noise=rms/.289_rprec*(ran3(seed)-0.5_rprec)
      v(jx,jy,jz)=noise*w_star/u_star*.01_rprec
      noise=rms/.289_rprec*(ran3(seed)-0.5_rprec)
      w(jx,jy,jz)=noise*w_star/u_star*.01_rprec
      noise=ran3(seed)-.5_rprec
! is this right for the neutral case? can't we put just theta=0., q=0. ?
!	      theta(jx,jy,jz)=(theta_s1+(z-z_i)*0.003)/T_scale
!	      noise=ran3(seed)-0.5
!	      q(jx,jy,jz)=q_mix+10.*noise*(1-z/z_i)*q_star
   end if
end do
end do
end do

!...BC for W
do jy=1,ny
do jx=1,nx
   w(jx,jy,1)=0._rprec
   w(jx,jy,nz)=0._rprec
end do
end do

!...BC for U, V
do jy=1,ny
do jx=1,nx
   u(jx,jy,nz)=u(jx,jy,nz-1)
   v(jx,jy,nz)=v(jx,jy,nz-1)
end do
end do
end subroutine ic
