!
! Copyright (C) 1991-2004  ; All Rights Reserved ; Colorado State University
! Colorado State University Research Foundation ; ATMET, LLC
! 
! This file is free software; you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software 
! Foundation; either version 2 of the License, or (at your option) any later version.
! 
! This software is distributed in the hope that it will be useful, but WITHOUT ANY 
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 
! PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this 
! code; if not, write to the Free Software Foundation, Inc., 
! 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!======================================================================================

Subroutine node_sendnbc (ifm,icm)

use mem_grid
use node_mod
use var_tables
use mem_basic

implicit none

integer :: ifm,icm
integer :: nm,i1,i2,j1,j2,k1,k2,ng,itype,mtp,iptr,nv
real, allocatable, save :: buffnest(:)
integer, save :: ncall=0, membuff,membuff_extra,nvar

itype=5

!______________________
!
!   First, before we send anything, let's post the receives.

do nm=1,nmachs
   irecv_req(nm)=0
   if (iget_paths(itype,ifm,nm).ne.not_a_node) then
      CALL par_get_noblock (node_buffs(nm)%lbc_recv_buff(1)  &
          ,node_buffs(nm)%nrecv ,5000+icm,machnum(nm),irecv_req(nm) )
   endif
enddo

! Send coarse grid points necessary for fine grid boundary interpolation
!   to fine grid nodes. Note that even though coarse grid points are sent,
!   ipaths is referenced by the fine grid, since all nests only have one
!   parent, not vice versa.


! Compute size of buffer needed and allocate if necessary
if(ncall == 0) then

   ncall=1
   membuff_extra=nvar*2+100
   membuff=0      

   do ng=1,ngrids
      do nm=1,nmachs
         if(ipaths(5,itype,ng,nm).ne.not_a_node) then
            i1=ipaths(1,itype,ng,nm)
            i2=ipaths(2,itype,ng,nm)
            j1=ipaths(3,itype,ng,nm)
            j2=ipaths(4,itype,ng,nm)
            k1=1
            k2=nnzp(ng)

            nvar=4 + num_scalar(ng)
            mtp=(i2-i1+1)*(j2-j1+1)*(k2-k1+1)
            membuff=max(membuff,mtp*nvar)
         endif
      enddo
   enddo

   membuff=membuff+membuff_extra
   allocate (buffnest(membuff))
   if(iprntstmt>=1)print*,'sending nesting condition:' &
      ,'Allocate buffer for:',my_rams_num,membuff,nvar
endif


do nm=1,nmachs

   isend_req(nm)=0

   if(ipaths(5,itype,ifm,nm).ne.not_a_node) then

      i1=ipaths(1,itype,ifm,nm)
      i2=ipaths(2,itype,ifm,nm)
      j1=ipaths(3,itype,ifm,nm)
      j2=ipaths(4,itype,ifm,nm)
      k1=1
      k2=nnzp(icm)

      mtp=(i2-i1+1)*(j2-j1+1)*(k2-k1+1)

! Put variables into buffer. All need coarse grid density weighting first.

      iptr=0
      CALL mknest_buff (1,basic_g(icm)%uc(1,1,1),buffnest(1+iptr)  &
          ,basic_g(icm)%dn0(1,1,1),mmzp(icm),mmxp(icm),mmyp(icm)  &
          ,mi0(icm),mj0(icm),i1,i2,j1,j2,k1,k2)
      iptr=iptr+mtp
      CALL mknest_buff (2,basic_g(icm)%vc(1,1,1),buffnest(1+iptr)  &
          ,basic_g(icm)%dn0(1,1,1),mmzp(icm),mmxp(icm),mmyp(icm)  &
          ,mi0(icm),mj0(icm),i1,i2,j1,j2,k1,k2)
      iptr=iptr+mtp
      CALL mknest_buff (3,basic_g(icm)%wc(1,1,1),buffnest(1+iptr)  &
          ,basic_g(icm)%dn0(1,1,1),mmzp(icm),mmxp(icm),mmyp(icm)  &
          ,mi0(icm),mj0(icm),i1,i2,j1,j2,k1,k2)
      iptr=iptr+mtp
      CALL mknest_buff (4,basic_g(icm)%pc(1,1,1),buffnest(1+iptr)  &
          ,basic_g(icm)%dn0(1,1,1),mmzp(icm),mmxp(icm),mmyp(icm)  &
          ,mi0(icm),mj0(icm),i1,i2,j1,j2,k1,k2)
      iptr=iptr+mtp

      do nv=1,num_scalar(ifm)
         CALL mknest_buff (5,scalar_tab(nv,icm)%var_p,buffnest(1+iptr)  &
             ,basic_g(icm)%dn0(1,1,1),mmzp(icm),mmxp(icm),mmyp(icm)  &
             ,mi0(icm),mj0(icm),i1,i2,j1,j2,k1,k2)         
         iptr=iptr+mtp
      enddo


      CALL par_init_put (node_buffs(nm)%lbc_send_buff(1)  &
                       ,node_buffs(nm)%nsend )

      CALL par_put_int (i1,1)
      CALL par_put_int (i2,1)
      CALL par_put_int (j1,1)
      CALL par_put_int (j2,1)
      CALL par_put_int (k1,1)
      CALL par_put_int (k2,1)
      CALL par_put_int (my_rams_num,1)
      CALL par_put_int (nvar,1)
      CALL par_put_int (iptr,1)

      CALL par_put_float (buffnest,iptr)

      CALL par_send_noblock (machnum(ipaths(5,itype,ifm,nm)),5000+icm,isend_req(nm))

   endif

enddo

return
END SUBROUTINE node_sendnbc

!##############################################################################
Subroutine mknest_buff (ivarn,ac,acs,den,m1,m2,m3,i0,j0  &
        ,i1,i2,j1,j2,k1,k2)

implicit none

integer :: ivarn,m1,m2,m3,i0,j0,i1,i2,j1,j2,k1,k2
real :: ac(m1,m2,m3),acs(0:k2-k1,0:i2-i1,0:j2-j1),den(m1,m2,m3)
integer :: i,j,k

!     ivarn = variable types 1- u
!                            2- v
!                            3- w
!                            4- p
!                            5- scalar

if(ivarn.eq.5) then
   do j=j1,j2
      do i=i1,i2
         do k=k1,k2
            acs(k-k1,i-i1,j-j1)=ac(k,i-i0,j-j0)*den(k,i-i0,j-j0)
         enddo
      enddo
   enddo
elseif(ivarn.eq.1) then
   do j=j1,j2
      do i=i1,i2
         do k=k1,k2
            acs(k-k1,i-i1,j-j1)=ac(k,i-i0,j-j0)  &
                 *((den(k,i-i0,j-j0)+den(k,i+1-i0,j-j0))*.5)
         enddo
      enddo
   enddo
elseif(ivarn.eq.2) then
   do j=j1,j2
      do i=i1,i2
         do k=k1,k2
            acs(k-k1,i-i1,j-j1)=ac(k,i-i0,j-j0)  &
                 *((den(k,i-i0,j-j0)+den(k,i-i0,j+1-j0))*.5)
         enddo
      enddo
   enddo
elseif(ivarn.eq.3) then
   do j=j1,j2
      do i=i1,i2
         do k=k1,k2-1
            acs(k-k1,i-i1,j-j1)=ac(k,i-i0,j-j0)  &
                 *((den(k,i-i0,j-j0)+den(k+1,i-i0,j-j0))*.5)
         enddo
      enddo
   enddo
elseif(ivarn.eq.4) then
   do j=j1,j2
      do i=i1,i2
         do k=k1,k2
            acs(k-k1,i-i1,j-j1)=ac(k,i-i0,j-j0)
         enddo
      enddo
   enddo
endif

return
END SUBROUTINE mknest_buff

!##############################################################################
Subroutine node_getnbc (ifm)

use mem_grid
use node_mod
use var_tables
use mem_scratch
use mem_basic
use mem_nestb

implicit none

integer :: ifm
integer, dimension(maxmach) :: i1c,i2c,j1c,j2c,k1c,k2c,iptv,iptc
integer :: itype,nm,ibytes,msgid,ihostnum,iptr,machf,nvar,nwords,nv  &
          ,nxc,nyc,nzc,mtp
real, allocatable, save :: buffnest(:)
integer, save :: nbuff_save=0

itype=5

!_____________________________________________________________________
!
!  First, let's make sure our sends are all finished and de-allocated

do nm=1,nmachs
   if(ipaths(5,itype,ifm,nm).ne.not_a_node)then
      CALL par_wait (isend_req(nm),ibytes,msgid,ihostnum)
   endif
enddo
!_____________________________________________________________________
!
!  Now, let's wait on our receives

do nm=1,nmachs
   if(iget_paths(itype,ifm,nm).ne.not_a_node)then
      CALL par_wait (irecv_req(nm),ibytes,msgid,ihostnum)
   endif
enddo
!_____________________________________________________________________
!


! Compute size of buffer needed and allocate if necessary

if(nbuff_nest > nbuff_save) then
   nbuff_save = nbuff_nest
   allocate (buffnest(nbuff_nest))
endif

!     From the fine grid nodes, get the coarse grid buffers,
!      interpolate the boundaries, and put them in the "b" array.

iptr=0
do nm=1,nmachs

   if(iget_paths(itype,ifm,nm).ne.not_a_node) then

      CALL par_assoc_buff (node_buffs(nm)%lbc_recv_buff(1)  &
                         ,node_buffs(nm)%nrecv)

      CALL par_get_int (i1c(nm),1)
      CALL par_get_int (i2c(nm),1)
      CALL par_get_int (j1c(nm),1)
      CALL par_get_int (j2c(nm),1)
      CALL par_get_int (k1c(nm),1)
      CALL par_get_int (k2c(nm),1)
      CALL par_get_int (machf,1)
      CALL par_get_int (nvar,1)
      CALL par_get_int (nwords,1)

      CALL par_get_float (buffnest(1+iptr),nwords)
      iptc(nm)=1+iptr
      iptv(nm)=0

      iptr=iptr+nwords
   endif
enddo

!   We have all the coarse grid info. Start looping through each variable.

do nv=1,nvar

!            First, construct coarse grid variable in scr1.

   CALL azero (ubound(scratch%scr1,1),scratch%scr1(1))
   CALL azero (ubound(scratch%scr2,1),scratch%scr2(1))
   do nm=1,nmachs
      if(iget_paths(itype,ifm,nm).ne.not_a_node) then
         nzc=k2c(nm)-k1c(nm)+1
         nxc=i2c(nm)-i1c(nm)+1
         nyc=j2c(nm)-j1c(nm)+1
         mtp=nzc*nxc*nyc

         CALL unmkbuff (scratch%scr1(1),buffnest(iptc(nm)+iptv(nm))  &
              ,maxnzp,maxnxp,maxnyp,nzc,nxc,nyc  &
              ,i1c(nm),i2c(nm),j1c(nm),j2c(nm),k1c(nm),k2c(nm))

         iptv(nm)=iptv(nm)+mtp
      endif
   enddo

!            Do the actual interpolation and put stuff into the "b" array

   if(nv.eq.1) then
      CALL par_bintp (scratch%scr1(1),scratch%scr2(1)  &
           ,basic_g(ifm)%dn0u(1,1,1)  &
           ,maxnzp,maxnxp,maxnyp,nnzp(ifm)  &
           ,mmzp(ifm),mmxp(ifm),mmyp(ifm)  &
           ,ifm,1,mi0(ifm),mj0(ifm),mibcon(ifm)  &
           ,nbounds(ifm)%bux(1,1,1),nbounds(ifm)%buy(1,1,1)  &
           ,nbounds(ifm)%buz(1,1,1))
   elseif(nv.eq.2) then
      CALL par_bintp (scratch%scr1(1),scratch%scr2(1)  &
           ,basic_g(ifm)%dn0v(1,1,1)  &
           ,maxnzp,maxnxp,maxnyp,nnzp(ifm)  &
           ,mmzp(ifm),mmxp(ifm),mmyp(ifm)  &
           ,ifm,2,mi0(ifm),mj0(ifm),mibcon(ifm)  &
           ,nbounds(ifm)%bvx(1,1,1),nbounds(ifm)%bvy(1,1,1)  &
           ,nbounds(ifm)%bvz(1,1,1))
   elseif(nv.eq.3) then
      CALL par_bintp (scratch%scr1(1),scratch%scr2(1)  &
           ,basic_g(ifm)%dn0(1,1,1)  &
           ,maxnzp,maxnxp,maxnyp,nnzp(ifm)  &
           ,mmzp(ifm),mmxp(ifm),mmyp(ifm)  &
           ,ifm,3,mi0(ifm),mj0(ifm),mibcon(ifm)  &
           ,nbounds(ifm)%bwx(1,1,1),nbounds(ifm)%bwy(1,1,1)  &
           ,nbounds(ifm)%bwz(1,1,1))
   elseif(nv.eq.4) then
      CALL par_bintp (scratch%scr1(1),scratch%scr2(1)  &
           ,basic_g(ifm)%dn0(1,1,1)  &
           ,maxnzp,maxnxp,maxnyp,nnzp(ifm)  &
           ,mmzp(ifm),mmxp(ifm),mmyp(ifm)  &
           ,ifm,4,mi0(ifm),mj0(ifm),mibcon(ifm)  &
           ,nbounds(ifm)%bpx(1,1,1),nbounds(ifm)%bpy(1,1,1)  &
           ,nbounds(ifm)%bpz(1,1,1))
   else
      CALL par_bintp (scratch%scr1(1),scratch%scr2(1)  &
           ,basic_g(ifm)%dn0(1,1,1)  &
           ,maxnzp,maxnxp,maxnyp,nnzp(ifm)  &
           ,mmzp(ifm),mmxp(ifm),mmyp(ifm)  &
           ,ifm,5,mi0(ifm),mj0(ifm),mibcon(ifm)  &
           ,nbounds(ifm)%bsx(1,1,1,nv-4),nbounds(ifm)%bsy(1,1,1,nv-4)  &
           ,nbounds(ifm)%bsz(1,1,1,nv-4))
   endif


enddo

return
END SUBROUTINE node_getnbc

!##############################################################################
Subroutine unmkbuff (ac,buff,max1,max2,max3,m1,m2,m3,i1,i2,j1,j2,k1,k2)

implicit none

integer :: max1,max2,max3,m1,m2,m3,i1,i2,j1,j2,k1,k2
real :: ac(max1,max2,max3),buff(0:m1-1,0:m2-1,0:m3-1)
integer :: i,j,k

do j=j1,j2
   do i=i1,i2
      do k=k1,k2
         ac(k,i,j)=buff(k-k1,i-i1,j-j1)
      enddo
   enddo
enddo

return
END SUBROUTINE unmkbuff

