module bottombc
use types,only:rprec
use param,only:nx,ny,ld
implicit none
! num_patch= numbr of patches, zo?=surface roughness for the patches types
! ptypes=number of pacthes types to be used, usually we use 2
integer,parameter::num_patch=1,ptypes=2
real(kind=rprec),parameter::zo1=.025_rprec,zo2=.25_rprec
integer,parameter::square=0
real(kind=rprec),parameter::zo_out=.025_rprec, zo_square=.25_rprec
integer,parameter::Nx_sq=5,Ny_sq=3
real(kind=rprec),parameter::t_scale=300._rprec
real(kind=rprec),parameter::theta_s1=300._rprec,theta_s2=300._rprec
real(kind=rprec),parameter::q_s1=11._rprec,q_s2=13._rprec
real(kind=rprec),dimension(nx,ny)::zo,T_s,q_s
!TS add for non-neutral case
real(kind=rprec),dimension(nx,ny)::phi_m,psi_m
integer,dimension(nx,ny)::patch
integer,dimension(num_patch)::patchnum

contains
subroutine patches()
! looks like this sets the roughness for the different patches
! and fills the lookup tables patch and patchnum
!use types,only:rprec
use param,only:z_i
implicit none

integer::i, j, patchtype, begini, endi, type
!integer,dimension(ptypes),intent(inout)::patchnum
!integer,dimension(nx,ny),intent(inout)::patch

! sc: 
! without this, some compilers may give total junk
! this was cause the u_avg to be WAY off, and so txz was WAY off, etc.
patchnum=0
!
if (square==1) then!for square patch
   zo(:,:)=zo_out/z_i
   patch(:,:)= 1
   zo(Nx/2-(Nx_sq-1)/2:Nx/2+(Nx_sq-1)/2,&
        Ny/2-(Ny_sq-1)/2:Ny/2+(Ny_sq-1)/2) = zo_square/z_i
   patch(Nx/2-(Nx_sq-1)/2:Nx/2+(Nx_sq-1)/2,&
        Ny/2-(Ny_sq-1)/2:Ny/2+(Ny_sq-1)/2) = 2
!TS   patchnum(2)=Nx_sq*Ny_sq
!TS   patchnum(1)=Nx*Ny-patchnum(2)
   patchnum=Nx*Ny
else!for regular patches
   do j=1,ny
      endi = 0
      type=1
      do patchtype=1,num_patch
         begini=endi+1
         endi=(nx*patchtype)/num_patch
         do i=begini,endi
            if (type.eq.1) then
               zo(i,j)=zo1/z_i
               T_s(i,j)=theta_s1/t_scale
               q_s(i,j)=q_s1
               patch(i,j)=type
            end if
            if(type.eq.2) then
               zo(i,j)=zo2/z_i
               T_s(i,j)=theta_s2/t_scale
               q_s(i,j)=q_s2
               patch(i,j)=type
            end if
            patchnum(type)=patchnum(type)+1
          end do
          if (type.eq.1) then
             type=2
          elseif (type.eq.2) then
             type=1
          end if
       end do       
      end do
end if
end subroutine patches

subroutine avgpatch(u_avg)
! computes the averaged value of a variable (at the wall) over a patch
! and assigns it to an nx X ny array

! sc: 
! note: this is inefficient: should calculate directly to u_avg
! --no maybe its good
!use types,only:rprec
use sim_param,only:u
!use bottombc,only:ptypes
implicit none
real(kind=rprec),dimension(nx,ny),intent(inout)::u_avg
!integer,dimension(ptypes),intent(in)::patchnum
!integer,dimension(nx,ny),intent(in)::patch
integer::i,j,k
real(kind=rprec),dimension(ptypes)::temp
!real(kind=rprec),dimension(ld,ny,1),intent(in)::u

temp=0._rprec
do j=1,ny
do i=1,nx
do k=1,ptypes
   if (patch(i,j).eq.k) then
      temp(patch(i,j))=temp(patch(i,j))+u(i,j,1)/real(patchnum(patch(i,j)))
   end if
end do
end do
end do

do j=1,ny
do i=1,nx
do k=1,ptypes
   if (patch(i,j).eq.k) then
      u_avg(i,j)=real(temp(k))
   end if
end do
end do
end do
end subroutine avgpatch

end module bottombc
