c @(#) nbill.f 1.7 (IDS) 1/21/95

C+
C*
C*   This file contains the following routines in the nodal part 
C*   of the host-node program :
C*   
C*   1) Main program (hnbill)
C*   2) subroutine get_bill_params - obtains the inputs for the nodal part
C*   3) subroutine send_bill_params - returns the interpolated results for
C*                                    either the bilinear or bicubic 
C*                                    interpolation.
C*   4) subroutine print_time       - computes the timings of the nodal part.
C*   5) subroutine nbill            - performs the lasagnia decomposition
C*                                    of the inputs among the various nodes
C-

C+
C  The main program initializes the nodal part of the host-node and then
C  enters a loop waiting for a request from the host part of the program. 
C  
C  Current requests are :-
C
C  1) MPBILL for bilinear interpolation
C  2) MPBIC  for bicubic interpolation (as in bic1nu.f)
C
C  The values are defined in slgcm.mp.  Note slgcm.mpp is no longer used.
C  Once the nodal program receives a request for performing an 
C  interpolation,  Node 0 alone issues a series of "gets" to get the
C  input parameters from the host. Then all nodes execute a call to
C  nbill to perform the interpolation.  Finally node 0 issues two sends
C  to return the two interpolated sets of values.  Finally print_time
C  is called to update timings.  The nodal program terminates when the
C  host program terminates.  The host program does not currently issue
C  MPQUIT to cause the nodal program to quit the loop.  Even if the loop
C  were exited the nodal program does not terminate till the host does.
C
C  Irrespective of the interpolation, the same parameters are inputs and
C  the same parameters are outputs.  [The bicubic interpolation performed
C  in bicpnu.f is still done by the host (serial code).]
C
C
C  Inputs : wina corresponds to ust (yuupt)
C           winb    "        "  vst (yvvpt)
C           rlatitude "      "  spuh (sphi)
C           rlongitude "     "  sluh (slam)
C           faclama   
C           faclamb    
C           facphia
C           facphib 
C           ja, jb      - range of indices
C           interpol    - bilinear (bicubic)
C  
C  Outputs : woua corresponds to c1a (yuupts or yvvpts)
C            woub     "       "  c2a (yvupts or yuvpts)
C  
C  The items in parantheses refers to the bicubic interpolation. For
C  further details about the quantities one is referred to Yong Li's 
C  notes.
C-

      program hnbill

      include "slgcm.slgprm"
c1      include "slgcm.mpp"  - Not used any more !!
      include "slgcm.mp"

      implicit none

        integer host, node, status, count

        dimension woua ( imp3 ,jm   )   !first output field
        dimension woub ( imp3 ,jm   )   !second output field
        real woua, woub
         dimension wina (imp3, jmp3)
         dimension winb (imp3, jmp3)
         dimension rlatitude ( imp3, jm)
         dimension rlongitude ( imp3,jm)
         real wina, winb
         real rlatitude, rlongitude

         integer dk, dj

C1         real   wina ( imp3 ,jmp3 ) !first input field
C1         real   winb ( imp3 ,jmp3 ) !second input field
C1         real   rlatitude  ( imp3 ,jm   ) ! arrays with coordinates on
C1        real   rlongitude  ( imp3 ,jm   ) ! which the subroutine interpolates
         real   faclama
         real   faclamb
         real   facphia
         real   facphib
         integer ja, jb
C+
C            cabtime - accumulated total cpu (busy) time
C            gcbtime - accumulated cpu (busy) time for global communication
C            lcbtime - accumulated cpu (busy) time for local communication
C            llitime -     "        "          "    "  calc. ii(i,j), jj(i,j)
C            brtime  -     "        "          "    "  calc. degree of circ.
C            cstime  -     "        "          "    "  creating southern part
C                                                       of segment
C            cntime  -     "        "          "    "  creating northern part
C                                                       of segment
C            reatime -     "        "          "    "  calc. local indices
C            fbtime  -     "        "          "    "   find_bc_jj 
C                                                      subroutine overhead
C            intime  -     "        "          "    "  cubit subroutine
C                                                      overhead
C            loitime - NOT USED
C-
         real*8 gcbtime, lcbtime, cabtime, llitime, loitime,
     &          brtime, cstime, cntime, reatime, fbtime, intime
         common /ntime/ cabtime, gcbtime, lcbtime, llitime, 
     &          loitime, brtime, cstime, cntime, reatime, fbtime,
     &          intime

         real niter, siter
         common /citer/ niter, siter

         integer tmp

         real*8 mp_busy_time
         integer mp_timer_clear
CT
C
C       function prototypes
C
        integer host_stat
        integer mp_host_addr
        integer mp_self_addr
        integer mp_msg_pending
C 
C       Initialize the host-node mode        
C
        call mp_enable_host()
        host = mp_host_addr()
        node = mp_self_addr()

        call mp_init()
D        print *, "In nbill 1st time, Node : ", node
CT
CT      Initialize timing variables
CT 
        cabtime = 0.0
        gcbtime = 0.0
        lcbtime = 0.0
        llitime = 0.0 
        loitime = 0.0
        brtime  = 0.0
        cstime  = 0.0
        cntime  = 0.0 
        reatime = 0.0 
        fbtime  = 0.0
        intime  = 0.0
CT

        count = 0
C
C       Receive request from host and enter loop.
C
        status = host_stat()
        do while (status .ne. MPQUIT)
C
C          Keep a count of the number of interpolations
C
           count = count + 1
D           print *, "Hello I am bill at count = ", count
C
C         Request host to stop serving as a server for serial i/o
C         Communication with host can only take place after this 
C         is done.
C
           call mp_suspend_server()
           
C
C          Get the input parameters
C
           if (node .eq. 0) then
              call get_bill_params(wina, faclama, facphia,   
     *                winb, faclamb, facphib, 
     *                rlatitude,rlongitude, ja, jb)

C           Debug code
C 
C1 107          format ( F16.8)
C1 
C1              print *, '+WINA'
C1              do dj = 1, jmp3
C1                 do dk = 1, imp3
C1                    print 107, wina(dk,dj)
C1                 end do
C1              end do
C1              print *, '+WINA OVER'
           endif
D           print *, "In nbill - 2nd time, Node : ", node

      
C
C          Perform the lasagnia decomposition and interpolate
C      
           call mp_sync()

D           print *, "I am in nodal main with ", status
              call nbill(wina, woua,faclama, facphia,
     *                winb, woub,faclamb, facphib, 
     *                rlatitude,rlongitude, ja, jb, status)

C
C          Update timings for the interpolation
C              
           if (node .eq. 5) then
             call print_time(count)
           endif
           call mp_suspend_server()

C
C         Return the interpolated results to the host.
C
           if (node .eq. 0) then
              call send_bill_params(woua, woub)
           end if

           status = host_stat() 
        end do

      end


C+
C  Subroutine get_bill_params
C  It issues a "get" or "receive" for every input from the host.
C  But conventionally these are the output parameters of this 
C  subroutine.
C
C     Input Parameters : None
C  
C     Output Parameters :
C           wina corresponds to ust (yuupt)
C           winb    "        "  vst (yvvpt)
C           rlatitude "      "  spuh (sphi)
C           rlongitude "     "  sluh (slam)
C           faclama   
C           faclamb    
C           facphia
C           facphib 
C           ja, jb      - range of indices
C           interpol    - bilinear (bicubic)
C  
C-

      subroutine get_bill_params(wina, faclama, facphia,   
     *                winb, faclamb, facphib, 
     *                rlatitude,rlongitude, ja, jb)

        include "slgcm.slgprm"
c1      include "slgcm.mpp"
        include "slgcm.mp"

         implicit none

         real   wina ( imp3 ,jmp3 ) !first input field
         real   winb ( imp3 ,jmp3 ) !second input field
         real   rlatitude  ( imp3 ,jm   ) ! arrays with coordinates on
         real   rlongitude  ( imp3 ,jm   ) ! which the subroutine interpolates
         real   faclama
         real   faclamb
         real   facphia
         real   facphib
         integer ja, jb

         integer host, status

        integer mp_host_addr
        integer mp_rcv_rblock
        integer mp_rcv_real
        integer mp_rcv_int

         host = mp_host_addr()

         status = mp_rcv_rblock(host, 0, wina, imp3, jmp3)
         status = mp_rcv_rblock(host, 0, winb, imp3, jmp3)
         status = mp_rcv_rblock(host, 0, rlatitude, imp3, jm)
         status = mp_rcv_rblock(host, 0, rlongitude, imp3, jm)
         status = mp_rcv_real(host, 0, faclama)
         status = mp_rcv_real(host, 0, faclamb)
         status = mp_rcv_real(host, 0, facphia)
         status = mp_rcv_real(host, 0, facphib)
         status = mp_rcv_int(host, 0, ja)
         status = mp_rcv_int(host, 0, jb)

      end


C+
C  Subroutine send_bill_params
C  It issues a "send" for the two interpolated results from the 
C  nodes. By convention these are the input parameters for this 
C  subroutine.
C
C     Input Parameters : 
C            woua corresponds to c1a (yuupts or yvvpts)
C            woub     "       "  c2a (yvupts or yuvpts)
C  
C     Output Parameters : None
C  
C-
      subroutine send_bill_params(woua, woub)
        include "slgcm.slgprm"
        include "slgcm.mp"

         implicit none
         integer host, status

         real   woua ( imp3 ,jm ) !first input field
         real   woub ( imp3 ,jm ) !second input field

        integer mp_host_addr
        integer mp_send_rblock

         host = mp_host_addr()

         status = mp_send_rblock(host, 0, woua, imp3, jm)
         status = mp_send_rblock(host, 0, woub, imp3, jm)
         
      end
c
c


C+
C  Subroutine print_time
C  This routine prints the updated timings for each request for 
C  bilinear (bicubic) interpolation.  The variable count
C  keeps track of the number of times the interpolation is
C  performed.  The latest timings are written to the file
C  'timings.dat'. The older timings are wiped out in the process.
C
C     Input Parameters : 
C            count   - number of times the interpolations was 
C                      performed.
C  
C     Output Parameters : None
C  
C-
      subroutine print_time(count)

      integer count

         real*8 gcbtime, lcbtime, cabtime, llitime, loitime,
     &          brtime, cstime, cntime, reatime, fbtime, intime
         common /ntime/ cabtime, gcbtime, lcbtime, llitime, 
     &          loitime, brtime, cstime, cntime, reatime, fbtime,
     &          intime

         real niter, siter
         common /citer/ niter, siter

      include "slgcm.slgprm"
      include "slgcm.mp"

        real*8 cobtime, ttime
        real*8 scabtime, sgcbtime, slcbtime
        real titer

D        print *, "Total time : ", cabtime
C
C      Find the average calculation and communication 
C      times.
C
        scabtime  = cabtime - lcbtime - gcbtime
        scabtime = scabtime /  count
        sgcbtime = gcbtime /  count
        slcbtime = lcbtime /  count
        
C
C      Find averaged total commn. time and
C      avgd. total commn. + calcn. time
C
        cobtime = slcbtime + sgcbtime
        ttime = cobtime + scabtime

 200      format (A, E22.15.3)
 201      format (A, F12.7)
D           print *, "+Writing Timings"

C
C          Overwrite contents of timings.dat with new timings.
C
           open (unit = 20, file = 'timings.dat', status = 'scratch' )
           close(20)
           open (unit = 20, file = 'timings.dat', status = 'new' )

           write (20, *), " IM : ", im, " JM : ", jm, " Count ", 
     &           count
           write (20, 200), " Global Comm. Time : ", sgcbtime
           write (20, 200), " Other Comm. Time : ",  slcbtime
           write (20, 201), " Total Comm. Time : ",  cobtime
           write (20, 201), " Calculation Time : ",  scabtime
           write (20, 201), " Total Time : ", ttime
           ttime = llitime / count
           write (20, 201), "Lat Long Calcns : ", ttime
           ttime = brtime / count
           write (20,201), "Degree of Circ. : ", ttime
           ttime = cstime / count
           write (20,201), "Comm. South : ", ttime
           ttime = cntime / count
           write (20,201), "Comm. North : ", ttime
           ttime = reatime / count
           write (20,201), "Local Indices : ", ttime
           ttime = fbtime / count
           write (20,201), "Find BC JJ overhead : ", ttime
           ttime = intime / count
           write (20,201), "Interpolation overhead : ", ttime
           ttime = llitime + brtime + cstime + cntime + reatime
     &             + fbtime + intime
           ttime = ttime / count
           write (20, 201), " Total Time (no fn. overhead) : ", ttime

C
C          Find the ratio of calculation to communication and 
C          write it to the file
C
           ttime = ttime / cobtime
           write (20, 201), " Adjusted Ratio : ", ttime
           ttime = scabtime / cobtime
           write (20, 201), " Ratio with overhead : ", ttime

           titer = niter / count
           write (20, 201), "North Iter : ", titer
           titer = siter / count
           write (20, 201), "South Iter : ", titer
           titer = (siter + niter) / count
           write (20, 201), "Total Iter : ", titer
           close(20)

      end


C+
C  Subroutine nbill
C  It performs a lasagnia decomposition on the input arrays      
C  wina, winb, woua, woub, rlongitude, rlatitude on Node 0.  
C  The width of lasagnia strip is precomputed and entered
C  into the header file slgcm.mp.  This is because it is
C  not possible to define the size of an array at run time
C  in f77.  The width depends on the resolution of the 
C  input arrays wina (or winb) and the 
C  number of nodes being used on a particular machine.
C  It then ships the partitions to the other nodes.  It 
C  calls the subroutine bill_mp or bic_mp depending on 
C  whether it is a bilinear or bicubic interpolation.
C  The other nodes then return the interpolated results
C  to node 0 after performing the required interpolation.
C
C
C  Inputs : wina corresponds to ust (yuupt)
C           winb    "        "  vst (yvvpt)
C           rlatitude "      "  spuh (sphi)
C           rlongitude "     "  sluh (slam)
C           faclama   
C           faclamb    
C           facphia
C           facphib 
C           ja, jb      - range of indices
C           interpol    - bilinear (bicubic)
C  
C  Outputs : woua corresponds to c1a (yuupts or yvvpts)
C            woub     "       "  c2a (yvupts or yuvpts)
C  
C-

      subroutine nbill(wina, woua,faclama, facphia,
     *                winb, woub,faclamb, facphib, 
     *                rlatitude,rlongitude, ja, jb, interpol)

      include "slgcm.slgprm"
c1      include "slgcm.mpp"
      include "slgcm.mp"

      implicit none

      real   wina ( imp3 ,jmp3 )   !first input field
      real   woua ( imp3 ,jm   )   !first output field
      real   winb ( imp3 ,jmp3 )   !second input field
      real   woub ( imp3 ,jm   )   !second output field
      real   rlatitude  ( imp3 ,jm   ) ! arrays with coordinates on
      real   rlongitude  ( imp3 ,jm   ) ! which the subroutine interpolates
      real   faclama
      real   faclamb
      real   facphia
      real   facphib
      integer ja, jb, interpol

C
C     variables used for subscripting arrays in debug mode
C
      integer dk, dj

C
C     Common block 
C     nx - number of nodes in the x (ii) axis - 0 for lasagna
C     nx0 - node address in x dirn (not used)
C     ny  - number of nodes in the y (jj) axis - 32 for our CM5
C     ny0 - node address in y dirn 
C     nynth - logical node address of northern neighbor
C     nysth - logical node address of southern neighbor
C     nxest - logical node address of eastern neighbor (not used)
C     nxwst - logical node address of western neighbor (not used)
C
      integer nx, nx0, ny, ny0, nynth, nysth, nxest, nxwst
      common /commvars/ nx,nx0,ny,ny0, nynth, nysth, nxest, nxwst
CT
CT       Timing code
CT 
         real*8 gcbtime, lcbtime, cabtime, llitime, loitime,
     &          brtime, cstime, cntime, reatime, fbtime, intime
         common /ntime/ cabtime, gcbtime, lcbtime, llitime, 
     &          loitime, brtime, cstime, cntime, reatime, fbtime,
     &          intime

CT
      integer count,               ! not used anymore
     &        lower,               ! lower bound of partition
     &        upper,               ! upper bound of partition
     &        i, tmp,              ! temporary variables
     &        range,               ! width of partition
     &        j, k,                ! counters
     &        dny0                 ! not used anymore

C
C     Function prototypes
C
      integer mp_rcv_int           
      integer mp_send_int
      integer mp_rcv_rblock
      integer mp_send_rblock
      integer mp_rcv_real
      integer mp_send_real

      parameter (count = imp3)    ! not used anymore

C
C      buffers to hold partitions. Array sections are not
C      allowed in f77.  
C
       dimension wina_d(imp3, max_col_p3)        
       dimension woua_d(imp3, max_col_p3)
       dimension winb_d(imp3, max_col_p3)
       dimension woub_d(imp3, max_col_p3)
       dimension rlatitude_d(imp3, max_col_p3)
       dimension rlongitude_d(imp3, max_col_p3)
       real wina_d, woua_d, winb_d, woub_d,
     &       rlongitude_d, rlatitude_d
C                      		         
CT
CT    function prototypes for the timing code.
CT  
      integer mp_timer_clear
      integer mp_timer_start
      integer mp_timer_stop
      real*8  mp_busy_time
CT

D      print *, "Starting"
C1 107   format ( I4, I4, 2(F16.8))

C     make the IO independent
       call mp_init()
D      print *, "In subrt. nbill"


C
C      Set the common block variables such as the address of
C      northern, southern neighbors etc.
C
       call set_lsd_commons(nx,nx0,ny,ny0, nynth, 
     &                      nysth, nxest, nxwst)

D          print *, "After lsd commons : ny0", ny0
       if (ny0 .eq. 0) then

D        print *, ny0, " : The first part"
C
C       Initialize with the bounds of the first partition (on node 0)
C
        lower = 1
        upper = max_col

C
C     Copy the partition into temporary buffers.  
C     This is done as array sections cannot be
C     passed in f77 as an argument.  The last 
C     node is dealt with separately, as the
C     input arrays cannot be split into 32
C     equal partitions.  
C
        do k = 2, ny -1                 ! for node 1 to 30
           lower = upper + 1            ! increment bounds for next partition
           upper = upper + max_col
           range = upper - lower + 1    

C
C          Copy into temporary buffers
C
           do j = 1, range              
              do i = 1, imp3
                 wina_d(i,j) = wina(i, lower +j-1)
                 winb_d(i,j) = winb(i, lower +j-1)
                 rlatitude_d(i,j) = rlatitude(i, lower +j-1)
                 rlongitude_d(i,j) = rlongitude(i, lower +j-1)
               enddo
            enddo
D            print *, ny0, " : The first part"

C
C   Ship the bounds of the partition followed by the
C   buffers themselves.
C
            tmp = mp_send_int(k - 1, 0, lower)
            tmp = mp_send_int(k - 1, 0, upper)
            tmp = mp_send_rblock(k - 1, 0, wina_d, imp3, max_col_p3)
            tmp = mp_send_rblock(k - 1, 0, winb_d, imp3, max_col_p3)
            tmp = mp_send_rblock(k - 1, 0, rlatitude_d, imp3, 
     &                            max_col_p3)
            tmp = mp_send_rblock(k - 1, 0, rlongitude_d, imp3, 
     &                            max_col_p3)

          enddo

D          print *, ny0, " : The first part"


C
C        For the last node the partitions of rlatitude
C        and rlongitude are of the same size, but differ
C        from that of the wina, winb partitions.
C
          lower = upper + 1
          upper = jm
          range = upper - lower + 1
            
          do j = 1, range
             do i = 1, imp3
                rlatitude_d(i,j) = rlatitude(i, lower +j-1)
                rlongitude_d(i,j) = rlongitude(i, lower +j-1)
             end do
          end do

          upper = jmp3
          range = upper - lower + 1

          do j = 1, range
             do i = 1, imp3
                wina_d(i,j) = wina(i, lower +j-1)
                winb_d(i,j) = winb(i, lower +j-1)
             end do
          end do

C
C     Need only the bounds of the decomposed wina, winb arrays
C     Ship the bounds and partitions to node 31 (last node)
C
            tmp = mp_send_int(ny - 1, 0, lower)
            tmp = mp_send_int(ny - 1, 0, upper)  
            tmp = mp_send_rblock(ny - 1, 0, wina_d, imp3, max_col_p3)
            tmp = mp_send_rblock(ny - 1, 0, winb_d, imp3, max_col_p3)
            tmp = mp_send_rblock(ny - 1, 0, rlatitude_d, imp3, 
     &           max_col_p3)
            tmp = mp_send_rblock(ny - 1, 0, rlongitude_d, imp3, 
     &           max_col_p3)


D          print *, ny0, " : The Second part"

C
C           Copy the first partition into a buffer 
C           (for node 0) before calling the 
C           interpolation routine bill_mp (bic_mp) on node 0
C
            lower = 1
            upper = max_col
            range = upper - lower + 1


            do j = 1, range
              do i = 1, imp3
                 wina_d(i,j) = wina(i, lower +j-1)
                 winb_d(i,j) = winb(i, lower +j-1)
                 rlatitude_d(i,j) = rlatitude(i, lower +j-1)
                 rlongitude_d(i,j) = rlongitude(i, lower +j-1)

               enddo
            enddo


D          print *, ny0, " : The Third part"

C
C        Send the other inputs (scalar quantities) to
C        all the remianing nodes.
C
          do k = 1, ny -1
            tmp = mp_send_int(k, 0, ja)
            tmp = mp_send_int(k, 0, jb)
            tmp = mp_send_real(k, 0, faclama)
            tmp = mp_send_real(k, 0, facphia)
            tmp = mp_send_real(k, 0, faclamb)
            tmp = mp_send_real(k, 0, facphib)
         enddo

       else
C
C     If not node 0 receive the various inputs.
C

           tmp = mp_rcv_int(0, 0, lower)
           tmp = mp_rcv_int(0, 0, upper)

           tmp = mp_rcv_rblock(0, 0, wina_d, imp3, max_col_p3)
           tmp = mp_rcv_rblock(0, 0, winb_d, imp3, max_col_p3)
           tmp = mp_rcv_rblock(0, 0, rlatitude_d, imp3, max_col_p3)
           tmp = mp_rcv_rblock(0, 0, rlongitude_d, imp3, max_col_p3)

           tmp = mp_rcv_int(k, 0, ja)
           tmp = mp_rcv_int(k, 0, jb)
           tmp = mp_rcv_real(k, 0, faclama)
           tmp = mp_rcv_real(k, 0, facphia)
           tmp = mp_rcv_real(k, 0, faclamb)
           tmp = mp_rcv_real(k, 0, facphib)

       endif

       call mp_sync()

       
D          print *,  ny0, " : The Fourth part"

C
C     Perform a bilinear (or bicubic) interpolation depending
C     on the request.
C
      if (interpol .eq. MPBILL) then
D         if (ny0 .eq. 0) then
D            print *, "I am going to call bill_mp"
D         endif 
CT
CT    Timing Code for measuring total time for interpolation.
	tmp = mp_timer_clear(2)
        if (tmp .eq. -1) then
           print *, 'Failure in clearing timer'
           stop
        end if
        tmp = mp_timer_start(2)
        if (tmp .eq. -1) then
           print *, 'Failure in starting timer'
           stop
        end if
CT    Timing Code 
CT

        call bill_mp(wina_d, woua_d,faclama, facphia,
     &               winb_d, woub_d,faclamb, facphib, 
     &               rlatitude_d,rlongitude_d, ja, jb,
     &               lower, upper)
CT    Timing Code 
        tmp = mp_timer_stop(2)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
           stop
        end if
C
C       Accumulate total times for bilinear interpolations.
C
        cabtime = cabtime + mp_busy_time(2) 
CT    Timing Code 

      else
D         if (ny0 .eq. 0) then
D            print *, "I am going to call bic_mp"
D         endif 
CT
CT    Timing Code for measuring total times in bicubic interpolation.
	tmp = mp_timer_clear(2)
        if (tmp .eq. -1) then
           print *, 'Failure in clearing timer'
           stop
        end if
        tmp = mp_timer_start(2)
        if (tmp .eq. -1) then
           print *, 'Failure in starting timer'
           stop
        end if
CT    Timing Code 
CT
        call bic_mp(wina_d, woua_d,faclama, facphia,
     &               winb_d, woub_d,faclamb, facphib, 
     &               rlatitude_d,rlongitude_d, ja, jb,
     &               lower, upper)

CT    Timing Code 
        tmp = mp_timer_stop(2)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
           stop
        end if
CT   Timing Code 
C
C     Accumulate total times for all bicubic interpolations
C
        cabtime = cabtime + mp_busy_time(2) 
      endif
  
C
C     If not node 0, return interpolated results.
C     
       if (ny0 .ne. 0) then
          tmp = mp_send_rblock(0, ny0, woua_d,    
     &                             imp3, max_col_p3)
          tmp = mp_send_rblock(0, ny0, woub_d,    
     &                             imp3, max_col_p3)

D          print *,  ny0, " : The fifth part"

       else
C
C         If node 0, copy the results into buffer woua, woub of 
C         size (imp3.jm). Receive the results from other nodes 
C         and copy them into the buffer at the respective places.
C
          lower = 1
	  upper = max_col
	  range = upper - lower + 1
	  
          do j = 1, range
             do i = 1, imp3
                woua(i, lower + j -1) = woua_d(i,j)
                woub(i, lower + j -1) = woub_d(i,j)
             end do
          end do

	  do k = 1, (ny - 1)             ! for nodes 1 to 31
	    
	    lower = upper + 1

            if (k .eq. (ny-1)) then
	      upper = jm
	    else
              upper = upper + max_col
	    endif
	
	    range = upper - lower + 1
	    
            
            tmp = mp_rcv_rblock(k, k, woua_d,
     &                          imp3, max_col_p3)
            tmp = mp_rcv_rblock(k, k, woub_d,
     &                          imp3, max_col_p3)

C
C           Copy all interpolated results into the original output
C           buffer, before returning to host program.
C
             do j = 1, range
                do i = 1, imp3
                   woua(i, lower + j -1) = woua_d(i,j)
                   woub(i, lower + j -1) = woub_d(i,j)
                end do
             end do
	 end do
	    
D         print *,  ny0, " : Finished the rest"
C1 101  format( 1X, 8('+', 4F16.8 //))

C1       print *, 'Count = ', count
C1       print *, '+WOUA'
C1       print 101, ((woua(dk, dj), dj = 1, jm), dk = 1, imp3)
C1       print *, '+WOUA OVER'
C1       print *, '+WOUB'
C1       print 101, ((woub(dk, dj), dj = 1, jm), dk = 1, imp3)
C1       print *, '+WOUB OVER'

       endif

C       End of routine bill
       return
       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Purpose:
c  Uses a bilinear scheme to interpolate the fields wina and winb to the point
c  (rlatitude,rlongitude) and puts the interpolated value in woua woub.
c  on a c-grid, faclam and facphi must be set as follows depending on
c  whether the points being used are the u v or fi points.
c  u points: faclam = 0.5; facphi = 1.0
c  v points: faclam = 1.0; facphi = 0.5
c  h points: faclam = 1.0; facphi = 1.0
c
c  ja and jb are the beggining and ending latitude indexes
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      




CC      subroutine bill_mp(wina_d, woua_d,faclama, facphia,
CC     &               winb_d, woub_d,faclamb, facphib, 
CC     &               rlatitude_d,rlongitude_d, ja, jb,
CC     &               lower, upper)
CC
CC
CC      end




