!
! 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 assign_node_subdomain ()

use mem_grid
use node_mod

implicit none

  integer :: ngr
  integer :: jnode, ncols

  ! Decompose all grids into subdomains

  do ngr = 1,ngrids
    ! Decompose the domain into rectangular sub-domains that have horizontal
    ! areas (ie, number of columns) match as closely as possible.
    CALL par_decomp_domain (nnxp(ngr),nnyp(ngr),nmachs,ixb(1,ngr),ixe(1,ngr) &
                           ,iyb(1,ngr),iye(1,ngr))

    if (nmachs .gt. 1) then
      print*,'!------------------------------------------------'
      print*,'!            Domain decomposition'
      print*,'!   grid# node# x-beg x-end y-beg y-end #cols'
      print*,'!------------------------------------------------'
      do jnode = 1,nmachs
         ncols= (1+ixe(jnode,ngr)-ixb(jnode,ngr))  &
               *(1+iye(jnode,ngr)-iyb(jnode,ngr))
         print '('' ! '',7i6,f12.4)', ngr,jnode,ixb(jnode,ngr) &
               ,ixe(jnode,ngr),iyb(jnode,ngr),iye(jnode,ngr),ncols
      enddo
      print*,'!------------------------------------------------'
      print*,''
    endif

    ! Figure out and record which sub domains are in the corners
    ! of the full domain
    do jnode = 1, nmachs
      if ((ixb(jnode,ngr) .eq.             2) .and. (iyb(jnode,ngr) &
         .eq.             2)) fd_sw_num(ngr) = jnode
      if ((ixb(jnode,ngr) .eq.             2) .and. (iye(jnode,ngr) &
         .eq. (nnyp(ngr)-1))) fd_nw_num(ngr) = jnode
      if ((ixe(jnode,ngr) .eq. (nnxp(ngr)-1)) .and. (iyb(jnode,ngr) &
         .eq.             2)) fd_se_num(ngr) = jnode
      if ((ixe(jnode,ngr) .eq. (nnxp(ngr)-1)) .and. (iye(jnode,ngr) &
         .eq. (nnyp(ngr)-1))) fd_ne_num(ngr) = jnode
    enddo
  enddo

return
END SUBROUTINE assign_node_subdomain

!##############################################################################
Subroutine node_init ()

use mem_grid
use node_mod

implicit none

!DEBUG
!  integer :: ngr, inode, itype, info
!DEBUG

  ! finish off the domain decomposition by filling in boundary information
  CALL par_decomp_bounds ()

  ! figure out which nodes need to talk to each other and record these paths
  CALL par_node_paths (maxgrds,ngrids,nxpmax,nypmax,nnxp,nnyp  &
   ,nxtnest,my_rams_num,maxmach,nmachs,ipm,jpm,ixb,ixe,iyb,iye,ipaths &
   ,iget_paths,ibnd,jbnd,not_a_node)

!DEBUG
!do ngr = 1,ngrids
!   do inode=1,nmachs
!      print'(a,i0,a)', 'NODE', my_rams_num, ':'
!      print'(a,i0,a)', 'NODE', my_rams_num, ':  ngr idn itype info ipaths '
!      do itype=1,7
!        print'(a,i0,a)', 'NODE', my_rams_num, ':'
!         do info=1,5
!            write(6,332) 'NODE', my_rams_num, ': ',ngr,inode,itype,info &
!               ,ipaths(info,itype,ngr,inode)
! 332        format(a,i0,a,5i6)
!         enddo
!      enddo
!   enddo
!enddo
!
!do ngr = 1,ngrids
!   do inode = 1,nmachs
!      print'(a,i0,a)', 'NODE', my_rams_num, ':'
!      print'(a,i0,a)', 'NODE', my_rams_num, ':  ngr isn itype iget_paths '
!      do itype = 1,6
!         write(6,334) 'NODE', my_rams_num, ': ',ngr,inode,itype &
!             ,iget_paths(itype,ngr,inode)
! 334     format(a,i0,a,4i6)
!      enddo
!   enddo
!enddo
!DEBUG


return
END SUBROUTINE node_init

!##############################################################################
Subroutine par_set_comm_buff_sizes ()

use mem_grid
use node_mod
use var_tables
use micro_prm, only: nkr

implicit none

  integer :: inode,ng,mp_nzp,icm,ifm,nestvar,nf,nc,nv
  integer :: isn,idn,itype
  integer :: i1,i2,j1,j2
  integer :: ixy,ixyz,numbuff,num_lbc_buff,num_nest_buff,num_feed_buff,memf
  integer :: npvar2,npvar3,npvar7
  integer :: nivar2,nivar3,nivars4,nivarw4,nivars3,nivarb4
  integer :: num_init_buff

  integer, dimension(nmachs) :: recv_buff_sizes

  ! Compute  send and receive buffer sizes. These will be maximum of
  ! long timestep, initialization, turbulence, nest boundaries, and nest feedback.
  !
  ! Small timestep will use same buffers as they are always smaller.
  !
  ! node_buffs(*)%nsend - sizes when I am the source node
  ! node_buffs(*)%nrecv - sizes when I am the destination node

  do inode=1,nmachs
     node_buffs(inode)%nsend=0
     node_buffs(inode)%nrecv=0
  enddo

  mp_nzp=0
  do ng=1,ngrids
     mp_nzp=max(mp_nzp,mmzp(ng))
  enddo

  
  do ng=1,ngrids
  
  !  Find number of nested variables to be communicated.
     icm=ng
     ifm=ng
     if(ng /= 1) icm=nxtnest(ifm)
     nestvar=4
     do nf=1,num_scalar(ifm)
        do nc=1,num_scalar(icm)
           if(scalar_tab(nf,ifm)%name==scalar_tab(nc,icm)%name)  &
                    nestvar=nestvar+1
        enddo
     enddo
  
  !  Find number of lbc variables to be communicated.
     npvar3=0 ; npvar2=0 ; npvar7=0
     do nv = 1,num_var(ng)
        if(vtab_r(nv,ng)%impt1 == 1 ) then
           if (vtab_r(nv,ng)%idim_type==2) npvar2=npvar2+1
           if (vtab_r(nv,ng)%idim_type==3) npvar3=npvar3+1
           if (vtab_r(nv,ng)%idim_type==7) npvar7=npvar7+1
        endif
     enddo
  
  !  Find number of init variables to be communicated.
     nivar2 = 0
     nivar3 = 0
     nivars4 = 0
     nivarw4 = 0
     nivars3 = 0
     nivarb4 = 0
     do nv = 1,num_var(ng)
        if(vtab_r(nv,ng)%impti == 1 ) then
           if (vtab_r(nv,ng)%idim_type==2) nivar2=nivar2+1
           if (vtab_r(nv,ng)%idim_type==3) nivar3=nivar3+1
           if (vtab_r(nv,ng)%idim_type==4) nivars4=nivars4+1
           if (vtab_r(nv,ng)%idim_type==5) nivarw4=nivarw4+1
           if (vtab_r(nv,ng)%idim_type==6) nivars3=nivars3+1
           if (vtab_r(nv,ng)%idim_type==7) nivarb4=nivarb4+1
        endif
     enddo

     ! Figure out buffer sizes when I am the sending node 
     do idn=1,nmachs
        num_lbc_buff=0
        num_init_buff=0
        num_nest_buff=0
        num_feed_buff=0
  
        itype=1
        i1=ipaths(1,itype,ng,idn)
        i2=ipaths(2,itype,ng,idn)
        j1=ipaths(3,itype,ng,idn)
        j2=ipaths(4,itype,ng,idn)
        if(i1.ne.0) then
           ixy=(i2-i1+1)*(j2-j1+1)
           ixyz=(i2-i1+1)*(j2-j1+1)*(mp_nzp)
           num_lbc_buff=ixyz*npvar3+ixy*npvar2+ixyz*nkr*npvar7  &
                +2*(npvar3+npvar2+npvar7+100)
           num_init_buff = (ixy*nivar2) + (ixyz*nivar3) + (nivars4*ixy*nzg*npatch) + &
                           (nivarw4*ixy*nzs*npatch) + (nivars3*ixy*npatch) + &
                           (nivarb4*ixyz*nkr) + &
                           2*(nivar2+nivar3+nivars4+nivarw4+nivars3+nivarb4+100)
        endif
  
        itype=5
        i1=ipaths(1,itype,ng,idn)
        i2=ipaths(2,itype,ng,idn)
        j1=ipaths(3,itype,ng,idn)
        j2=ipaths(4,itype,ng,idn)
        if(i1.ne.0) then
           ixyz=(i2-i1+1)*(j2-j1+1)*(mp_nzp)
           num_nest_buff=ixyz*nestvar+2*(nestvar+100)
        endif
  
        itype=6
        i1=ipaths(1,itype,ng,idn)
        i2=ipaths(2,itype,ng,idn)
        j1=ipaths(3,itype,ng,idn)
        j2=ipaths(4,itype,ng,idn)
        if(i1.ne.0) then
           ixyz=(i2-i1+1)*(j2-j1+1)*(mp_nzp)
           num_feed_buff=ixyz*nestvar+2*(nestvar+100)
        endif

        node_buffs(idn)%nsend = max(node_buffs(idn)%nsend  &
             ,num_lbc_buff,num_init_buff,num_nest_buff,num_feed_buff)
     enddo
  enddo

  ! At this point we have placed the max buffer sizes for sending data to 
  ! other nodes considering all grids. Ie, one send size will account for 
  ! what is needed for all grids.

  ! Now figure out buffer sizes when this node is the destination node.
  !
  ! We need information from the other nodes in order to accomplish
  ! this since our copy of ipaths only accounts for when we are the
  ! sending node (and when another node is sending to us, the buffer
  ! size could be different). Can use an MPI gather to collect the
  ! send buffer sizes from the other nodes and load them (in correct
  ! order) into our receive buffer size vars.
  !
  ! The basic idea is that for node 1, the receive buffer sizes have been 
  ! computed in the first entry of all the other nodes' send buffer size 
  ! arrays (since they are sending to node 1). For node 2, we want to 
  ! collect the second entry of all the other nodes' send buffer size 
  ! arrays. For node 3, collect the third entries, and so on.
  !
  ! The algorithm then becomes:
  !   For each node number
  !     Issue a gather sending your individual send buffer size and 
  !     collecting into a buffer large enough to hold all the numbers 
  !     from the other nodes
  !   Once the gather loop is finished, dole out the numbers in your 
  !   buffer into your receive buffer size entries.
  !
  ! Every node must issue a gather to all nodes including itself. This 
  ! makes the gathers match up across the collection of processes and 
  ! each node receives the entries it needs.
  !
  ! Note that this algorithm effectively transposes the send matrix 
  ! (formed by the send arrays from all the nodes) into the receive matrix.
    
  do isn = 1, nmachs
    CALL par_gather_ints (node_buffs(isn)%nsend, recv_buff_sizes, 1, machnum(isn))
  enddo

  do isn = 1, nmachs
    node_buffs(isn)%nrecv = recv_buff_sizes(isn)
  enddo

!DEBUG
!print'(a,i0,a,100i12)','DEBUG: NODE', my_rams_num,': send buffers:' &
!  ,(node_buffs(isn)%nsend,isn=1,nmachs)
!print'(a,i0,a,100i12)','DEBUG: NODE', my_rams_num,': recv buffers:' &
!  ,(node_buffs(isn)%nrecv,isn=1,nmachs)
!DEBUG
  
  !       Check nest boundary receive buffer size
  itype=5
  do idn=1,nmachs
     newbuff_nest1(idn)=1
     nbuff_nest1(idn)=0
  
     do ng=1,ngrids
        numbuff=0
        icm=nxtnest(ng)

        i1=ipaths(1,itype,ng,idn)
        i2=ipaths(2,itype,ng,idn)
        j1=ipaths(3,itype,ng,idn)
        j2=ipaths(4,itype,ng,idn)
        memf=(i2-i1+1)*(j2-j1+1)*(mp_nzp)*nestvar
        numbuff=numbuff+memf

        nbuff_nest1(idn)=max(nbuff_nest1(idn),numbuff)
     enddo
  
  enddo


return
END SUBROUTINE par_set_comm_buff_sizes

!##############################################################################
Subroutine par_decomp_bounds ()

use mem_grid
use node_mod
use micro_prm, only: nkr

implicit none

  integer :: igrid, inode
  integer :: nxpts, nypts

  !        Compute various subdomain boundary numbers for this node
  !             mxbeg,mybeg,mxend,myend - portions of full domain that node will have
  !                                     - includes overlap region
  !             mi0,mj0  - subdomain offsets relative to full domain
  !             mia,miz,mja,mjz - subdomain "compute" points,
  !                         or normal thermodynamic tendency points (2-nx, 2-ny for
  !                          non-parallel run
  !             mibcon - flag denoting if real boundary is on subdomain
  !                       bit 1=west, bit 2=east, bit 3=south, bit 4=north
  !
  !        Subdomain boundary numbers for parallel HDF5 file IO
  !             mem_read(igrid)%xblock, mem_read(igrid)%yblock     - memory block size for read operation
  !             mem_write(igrid)%xblock, mem_write(igrid)%yblock   - memory block size for write operation
  !             mem_read(igrid)%xoff, mem_read(igrid)%yoff         - memory subdomain offset relative to subdomain for read operation
  !             mem_write(igrid)%xoff, mem_write(igrid)%yoff       - memory subdomain offset relative to subdomain for write operation
  !             file_read(igrid)%xblock, file_read(igrid)%yblock   - file block size for read operation
  !             file_write(igrid)%xblock, file_write(igrid)%yblock - file block size for write operation
  !             file_read(igrid)%xoff, file_read(igrid)%yoff       - file subdomain offset relative to full domain for read operation
  !             file_write(igrid)%xoff, file_write(igrid)%yoff     - file subdomain offset relative to full domain for write operation
  !             file_xchunk, file_ychunk             - chunk size write operation

  ! This routine assumes that boundaries are all 1 point wide

  ! need chunk size to be max(all mmxp) by max(all mmyp)
  do igrid = 1, ngrids
    mibcon(igrid) = 0

    file_xchunk(igrid) = 0
    file_ychunk(igrid) = 0

    do inode = 1, nmachs
      ! Find the size of inode's sub-domain
      nxpts = (ixe(inode,igrid) - ixb(inode,igrid)) + 3
      nypts = (iye(inode,igrid) - iyb(inode,igrid)) + 3

      ! Save these if they are the max so far
      if (nxpts > file_xchunk(igrid)) then
        file_xchunk(igrid) = nxpts
      endif
      if (nypts > file_ychunk(igrid)) then
        file_ychunk(igrid) = nypts
      endif

      ! If we are on the inode matched with this process, then record
      ! and calculate all of the other boundary information
      !
      ! Calculate the numbers as if the sub-domain doesn't share a full
      ! domain boundary. Then check for the full domain boundary and adjust
      ! numbers accordingly.
      !
      ! ixb,ixe,iyb,iye denote the compute region of the subdomain
      if (inode .eq. my_rams_num) then
        mmxp(igrid) = nxpts
        mmyp(igrid) = nypts
        mmzp(igrid) = nnzp(igrid)

        mmxyzp(igrid) = mmxp(igrid) * mmyp(igrid) * mmzp(igrid)
        mmxysp(igrid) = mmxp(igrid) * mmyp(igrid) * (nzg+nzs+3) * npatch
        mmxyp(igrid)  = mmxp(igrid) * mmyp(igrid)
        mmxyzbp(igrid) = mmxyzp(igrid) * nkr

        ! mxbeg, mxend, mybeg, myend extend 1 beyond ixb,ixe,iyb,iye
        mxbeg(igrid) = ixb(inode,igrid) - 1
        mxend(igrid) = ixe(inode,igrid) + 1
        mybeg(igrid) = iyb(inode,igrid) - 1
        myend(igrid) = iye(inode,igrid) + 1

        ! offsets are 2 less than ixb,iyb
        mi0(igrid) = ixb(inode,igrid) - 2
        mj0(igrid) = iyb(inode,igrid) - 2

        ! mia,mja are 2
        ! miz,mjz are 2 + (i[xy]e-i[xy]b)
        mia(igrid) = 2
        miz(igrid) = 2 + (ixe(inode,igrid) - ixb(inode,igrid))
        mja(igrid) = 2
        mjz(igrid) = 2 + (iye(inode,igrid) - iyb(inode,igrid))

        ! generate io sizes and offsets
        ! read operation uses entire sub-domain (compute area plus boundaries)
        mem_read(igrid)%xblock = mmxp(igrid)
        mem_read(igrid)%yblock = mmyp(igrid)
        mem_read(igrid)%xoff   = 0
        mem_read(igrid)%yoff   = 0

        file_read(igrid)%xblock = mem_read(igrid)%xblock
        file_read(igrid)%yblock = mem_read(igrid)%yblock
        file_read(igrid)%xoff   = mi0(igrid)
        file_read(igrid)%yoff   = mj0(igrid)
 
        ! write operations uses sub-domain area minus interior boundaries
        ! set values here as if sub-domain is completely surrounded by interior boundaries
        mem_write(igrid)%xblock = (miz(igrid) - mia(igrid)) + 1
        mem_write(igrid)%yblock = (mjz(igrid) - mja(igrid)) + 1
        mem_write(igrid)%xoff   = 1
        mem_write(igrid)%yoff   = 1

        file_write(igrid)%xblock = mem_write(igrid)%xblock
        file_write(igrid)%yblock = mem_write(igrid)%yblock
        file_write(igrid)%xoff   = mi0(igrid) + 1
        file_write(igrid)%yoff   = mj0(igrid) + 1

        ! check west boundary
        if (ixb(inode,igrid) .eq. 2) then
          mibcon(igrid) = mibcon(igrid) + 1
          ! add west boundary back into the io write region
          mem_write(igrid)%xblock = mem_write(igrid)%xblock + 1
          mem_write(igrid)%xoff = mem_write(igrid)%xoff - 1

          file_write(igrid)%xblock = file_write(igrid)%xblock + 1
          file_write(igrid)%xoff = file_write(igrid)%xoff - 1
        endif
        
        ! check east boundary
        if (ixe(inode,igrid) .eq. (nnxp(igrid)-1)) then
          mibcon(igrid) = mibcon(igrid) + 2
          ! add east boundary back into the io write region
          mem_write(igrid)%xblock = mem_write(igrid)%xblock + 1

          file_write(igrid)%xblock = file_write(igrid)%xblock + 1
        endif
        
        ! check south boundary
        if (iyb(inode,igrid) .eq. 2) then
          mibcon(igrid) = mibcon(igrid) + 4
          ! add south boundary back into the io write region
          mem_write(igrid)%yblock = mem_write(igrid)%yblock + 1
          mem_write(igrid)%yoff = mem_write(igrid)%yoff - 1

          file_write(igrid)%yblock = file_write(igrid)%yblock + 1
          file_write(igrid)%yoff = file_write(igrid)%yoff - 1
        endif
        
        ! check north boundary
        if (iye(inode,igrid) .eq. (nnyp(igrid)-1)) then
          mibcon(igrid) = mibcon(igrid) + 8
          ! add north boundary back into the io write region
          mem_write(igrid)%yblock = mem_write(igrid)%yblock + 1

          file_write(igrid)%yblock = file_write(igrid)%yblock + 1
        endif
        
      endif
    enddo
  enddo

return
END SUBROUTINE par_decomp_bounds
