!...Compute rotational convective term in physical space  (X-mom eq.)
!...For staggered grid
!...u2, du2, and du4 are on W nodes, rest on UVP nodes
!...(Fixed k=Nz plane on 1/24)
!...Corrected near wall on 4/14/96 {w(DZ/2)=0.5*w(DZ)}
! sc: better to have 0.25*w(dz), since really parabolic.
subroutine convec (cx,cy,cz,u1,u2,u3,du1d2,du1d3,du2d1,du2d3,du3d1,du3d2,jt)
use types,only:rprec
use param,only:ld,ld_big,nx,nx2,ny,ny2,nz
use fft
implicit none
integer,intent(in)::jt
real(kind=rprec),dimension(ld,ny,nz)::u1,u2,u3,cx,cy,cz
real(kind=rprec),dimension(ld,ny,nz)::du1d2,du1d3,du2d1,du2d3,du3d1,du3d2
real(kind=rprec),dimension(ld_big,ny2,nz)::cc_big,u1_big,u2_big,u3_big
real(kind=rprec),dimension(ld_big,ny2,nz)::vort1_big,vort2_big,vort3_big
! should these be ``common'ed'' to save RAM?
! look elsewhere for possble ``partners''
! vorticity saves RAM and FFT time!
! note the way the vorticity/uvp/w node stuff works out nicely
integer::jz
real(kind=rprec)::ignore_me,const,factor

!..Recall dudz, and dvdz are on UVP node for k=1 only
!...So du2 does not vary from arg2a to arg2b in 1st plane (k=1)

! sc: it would be nice to NOT have to loop through slices here
! Loop through horizontal slices

const=1._rprec/(nx*ny)
!$omp parallel do default(shared) private(jz)		
do jz=1,nz
! use cx,cy,cz for temp storage here! 
   cx(:,:,jz)=const*u1(:,:,jz)
   cy(:,:,jz)=const*u2(:,:,jz)
   cz(:,:,jz)=const*u3(:,:,jz)
! do forward fft on normal-size arrays
   call rfftwnd_f77_one_real_to_complex(forw,cx(:,:,jz),ignore_me)
   call rfftwnd_f77_one_real_to_complex(forw,cy(:,:,jz),ignore_me)
   call rfftwnd_f77_one_real_to_complex(forw,cz(:,:,jz),ignore_me)
! zero pad: padd takes care of the oddballs
   call padd(u1_big(:,:,jz),cx(:,:,jz))
   call padd(u2_big(:,:,jz),cy(:,:,jz))
   call padd(u3_big(:,:,jz),cz(:,:,jz))
! now do the same, but with the vorticity!
   cx(:,:,jz)=const*(du3d2(:,:,jz)-du2d3(:,:,jz))
   cy(:,:,jz)=const*(du1d3(:,:,jz)-du3d1(:,:,jz))
   cz(:,:,jz)=const*(du2d1(:,:,jz)-du1d2(:,:,jz))

! sc: need to modify the vorticity at the wall to accomodate the 
! correction to the finite differencing error
! want to put du1d3(:,:,2)=du1d3(:,:,2)+u_star/vonK/dz*(1-log(3.))
   call rfftwnd_f77_one_real_to_complex(forw,cx(:,:,jz),ignore_me)
   call rfftwnd_f77_one_real_to_complex(forw,cy(:,:,jz),ignore_me)
   call rfftwnd_f77_one_real_to_complex(forw,cz(:,:,jz),ignore_me)
   call padd(vort1_big(:,:,jz),cx(:,:,jz))
   call padd(vort2_big(:,:,jz),cy(:,:,jz))
   call padd(vort3_big(:,:,jz),cz(:,:,jz))

! Back to physical space
! the normalization should be ok...
   call rfftwnd_f77_one_complex_to_real(back_big,u1_big(:,:,jz),ignore_me)
   call rfftwnd_f77_one_complex_to_real(back_big,u2_big(:,:,jz),ignore_me)
   call rfftwnd_f77_one_complex_to_real(back_big,u3_big(:,:,jz),ignore_me)
   call rfftwnd_f77_one_complex_to_real(back_big,vort1_big(:,:,jz),ignore_me)
   call rfftwnd_f77_one_complex_to_real(back_big,vort2_big(:,:,jz),ignore_me)
   call rfftwnd_f77_one_complex_to_real(back_big,vort3_big(:,:,jz),ignore_me)
end do
!$omp end parallel do

! CX
! redefinition of const
const=1._rprec/(nx2*ny2)

! the cc's contain the normalization factor for the upcoming fft's
cc_big(:,:,1)=const*(u2_big(:,:,1)*(-vort3_big(:,:,1))&
     +0.5_rprec*u3_big(:,:,2)*(vort2_big(:,:,2)))
! according to the comment at the beginning, dudz, dvdz are on the UVP nodes
! for jz=1.  If that's so, then why do we use the 0.5?
! whats going on at du1d3_big(:,:,1), etc.?
! I don't think these quantites are 0! Or, if don't want them, then why 0.5?

!$omp parallel do default(shared) private(jz)
do jz=2,nz-1
   cc_big(:,:,jz)=const*(u2_big(:,:,jz)*(-vort3_big(:,:,jz))&
        +0.5_rprec*(u3_big(:,:,jz+1)*(vort2_big(:,:,jz+1))&
        +u3_big(:,:,jz)*(vort2_big(:,:,jz))))
end do
!$omp end parallel do

cc_big(:,:,nz)=const*(u2_big(:,:,nz)*(-vort3_big(:,:,nz))&
     +u3_big(:,:,nz)*(vort2_big(:,:,nz)))
! doesn't this sort of assume a no stress BC on top?
! so what happens with the sponge thing?

! Loop through horizontal slices
!$omp parallel do default(shared) private(jz)	
do jz=1,nz
   call rfftwnd_f77_one_real_to_complex(forw_big,cc_big(:,:,jz),ignore_me)
! un-zero pad
! note: cc_big is going into cx!!!!
   call unpadd(cx(:,:,jz),cc_big(:,:,jz))
! Back to physical space
   call rfftwnd_f77_one_complex_to_real(back,cx(:,:,jz),ignore_me)

! sc: why aren't we using this?
! call filt_da(cx)
end do
!$omp end parallel do

! CY
! const should be 1./(nx2*ny2) here
! the cc's contain the normalization factor for the upcoming fft's
cc_big(:,:,1)=const*(u1_big(:,:,1)*(vort3_big(:,:,1))&
     +0.5_rprec*u3_big(:,:,2)*(-vort1_big(:,:,2)))
! again, what's the deal with the 0.5 and not incl (:,:,1)?
 
!$omp parallel do default(shared) private(jz)
do jz=2,nz-1
   cc_big(:,:,jz)=const*(u1_big(:,:,jz)*(vort3_big(:,:,jz))&
        +0.5_rprec*(u3_big(:,:,jz+1)*(-vort1_big(:,:,jz+1))&
        +u3_big(:,:,jz)*(-vort1_big(:,:,jz))))
end do
!$omp end parallel do
cc_big(:,:,nz)=const*(u1_big(:,:,nz)*(vort3_big(:,:,nz))&
     +u3_big(:,:,nz)*(-vort1_big(:,:,nz)))

!$omp parallel do default(shared) private(jz)		
do jz=1,nz
   call rfftwnd_f77_one_real_to_complex(forw_big,cc_big(:,:,jz),ignore_me)
 ! un-zero pad
! note: cc_big is going into cy!!!!
   call unpadd(cy(:,:,jz),cc_big(:,:,jz))

! Back to physical space
   call rfftwnd_f77_one_complex_to_real(back,cy(:,:,jz),ignore_me)

!! sc: why aren't we using this?
! call filt_da(cy)
end do
!$omp end parallel do

! CZ
! There is no convective acceleration of w at wall or at top.
cc_big(:,:,1)=0.
!$omp parallel do default(shared) private(jz)
do jz=2,nz-1
   cc_big(:,:,jz)=const*0.5_rprec*(&
        (u1_big(:,:,jz)+u1_big(:,:,jz-1))*(-vort2_big(:,:,jz))&
        +(u2_big(:,:,jz)+u2_big(:,:,jz-1))*(vort1_big(:,:,jz))&
         )
!	end if
end do
!$omp end parallel do
cc_big(:,:,nz)=0._rprec ! according to JDA paper p.242

! Loop through horizontal slices
!$omp parallel do default(shared) private(jz)		
do jz=1,nz
   call rfftwnd_f77_one_real_to_complex(forw_big,cc_big(:,:,jz),ignore_me)

! un-zero pad
! note: cc_big is going into cz!!!!
   call unpadd(cz(:,:,jz),cc_big(:,:,jz))

! Back to physical space
   call rfftwnd_f77_one_complex_to_real(back,cz(:,:,jz),ignore_me)

! sc: why aren't we using this?
! call filt_da(cz)
end do
!$omp end parallel do
end subroutine convec
