c %Z% %M% %I% (IDS) %G%

      program hnbill

      include "slgcm.slgprm"
c1      include "slgcm.mpp"
      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

        integer host_stat
        integer mp_host_addr
        integer mp_self_addr
        integer mp_msg_pending
        
        call mp_enable_host()
        host = mp_host_addr()
        node = mp_self_addr()

        call mp_init()
D        print *, "In nbill 1st time, Node : ", node

        count = 0
        status = host_stat()
        do while (status .ne. MPQUIT)
           count = count + 1
D           print *, "Hello I am bill at count = ", count
           call mp_suspend_server()
           
           if (node .eq. 0) then
              call get_bill_params(wina, faclama, facphia,   
     *                winb, faclamb, facphib, 
     *                rlatitude,rlongitude, ja, jb)

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

      
      
           call mp_sync()

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

C!
C!           Uncomment if nodes are using UNIX system calls
           call mp_suspend_server()

C1           if (node .eq. 0) then
C1              call send_bill_params(woua, woub)
C1           end if

C1           do while (mp_msg_pending(host,0) .eq. 0)
C1              count = count + 1
C1           end do

C1           status = host_stat() 
        end do

      end

      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

C1         host = cmmd_host_node()
         host = mp_host_addr()

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

         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

      subroutine send_bill_params(woua, woub)
        include "slgcm.slgprm"
c1      include "slgcm.mpp"
        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

C1         host = cmmd_host_node()
         host = mp_host_addr()

C1         status = cmmd_send(host, 0, woua, imp3 * jm * 4)
C1         status = cmmd_send(host, 0, woub, imp3 * jm * 4)
         status = mp_send_rblock(host, 0, woua, imp3, jm)
         status = mp_send_rblock(host, 0, woub, imp3, jm)
         
      end
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
      integer dk, dj

      integer nx, nx0, ny, ny0, nynth, nysth, nxest, nxwst
      common /commvars/ nx,nx0,ny,ny0, nynth, nysth, nxest, nxwst

      integer count,
     &        lower,
     &        upper,
     &        i, tmp,
     &        range,
     &        j, k, dny0

      integer mp_rcv_int
      integer mp_send_int
      integer mp_rcv_rblock
      integer mp_send_rblock

      parameter (count = imp3)

       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                      		          ! the fields wina and winb


C     make the IO independent

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


       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"
        lower = 1
        upper = max_col

        do k = 2, ny -1
           lower = upper + 1
           upper = 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)
C1                 print 107, k, lower, wina(i, lower + j-1), wina_d(i,j)    
               enddo
            enddo
D            print *, ny0, " : The first part"

            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, wina_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
C3
C3        Debug
D          print *, ny0, " : The first part"
C3        Debug
C3          
          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)
C1                 print 107, ny-1, lower, 
C1     &            wina(i, lower + j-1), wina_d(i,j)                 
             end do
          end do

c
c     Need only the bounds of the decomposed wina, winb arrays
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, wina_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)


C3
C3        Debug
D          print *, ny0, " : The Second part"
C3        Debug
C3          
            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)
C1                 print 107, wina(i,j), wina_d(i,j)

               enddo
            enddo
C3
C3        Debug
D          print *, ny0, " : The Third part"
C3        Debug
C3          

       else
           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)
       endif

C1       call cmmd_sync_with_nodes()
       call mp_sync()
C3
C3        Debug
D          print *,  ny0, " : The Fourth part"
C3        Debug
C3          

C     Interpolate
      if (interpol .eq. MPBILL) then
        call bill_mp(wina_d, woua_d,faclama, facphia,
     &               winb_d, woub_d,faclamb, facphib, 
     &               rlatitude_d,rlongitude_d, ja, jb,
     &               lower, upper)
      else
        call bic_mp(wina_d, woua_d,faclama, facphia,
     &               winb_d, woub_d,faclamb, facphib, 
     &               rlatitude_d,rlongitude_d, ja, jb,
     &               lower, upper)
      endif
       
       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)
C3
C3        Debug
D          print *,  ny0, " : The fifth part"
C3        Debug
C3          
       else
          lower = 1
	  upper = max_col
	  range = upper - lower + 1
	  
	  do k = 1, ny
	    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
	    
	    lower = upper + 1
	    if (k .eq. ny) then 
	      upper = jm
	    else
              upper = upper + max_col
	    endif
	
	    range = upper - lower + 1
	    
	    if (k .ne. 1) then
	      tmp = mp_rcv_rblock(k -1, k -1, woua_d,
     &                            imp3, max_col_p3)
      	      tmp = mp_rcv_rblock(k -1, k -1, woub_d,
     &                            imp3, max_col_p3)
	    end if
	 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

      

      subroutine bill_mp(wina_d, woua_d,faclama, facphia,
     *                   winb_d, woub_d,faclamb, facphib, 
     *                   rlatitude_d, rlongitude_d, ja, jb,
     *                   lo, up)

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

      implicit none

      integer nx, nx0, ny, ny0, nynth, nysth, nxest, nxwst
      common /commvars/ nx,nx0,ny,ny0, nynth, nysth, nxest, nxwst

      integer lo, up
      integer ja, jb
      real    faclama,
     &         faclamb,
     &         facphia,
     &         facphib

C1       real wina_d(imp3, max_col)
C1       real woua_d(imp3, max_col)
C1       real winb_d(imp3, max_col)
C1       real woub_d(imp3, max_col)
C1       real rlatitude_d(imp3, max_col)
C1       real rlongitude_d(imp3, max_col)
       
C1      dimension tem1_d(imp3, max_col)
C1      dimension tem2_d(imp3, max_col)
C1      dimension core1_d(imp3, max_col)
C1      dimension ii_d(imp3, max_col)
C1      dimension jj_d(imp3, max_col)

       real wina_d(imp3, max_col_p3)
       real woua_d(imp3, max_col_p3)
       real winb_d(imp3, max_col_p3)
       real woub_d(imp3, max_col_p3)
       real rlatitude_d(imp3, max_col_p3)
       real rlongitude_d(imp3, max_col_p3)
       
      dimension tem1_d(imp3, max_col_p3)
      dimension tem2_d(imp3, max_col_p3)
      dimension core1_d(imp3, max_col_p3)
      dimension ii_d(imp3, max_col_p3)
      dimension jj_d(imp3, max_col_p3)

      real tem1_d,
     &      tem2_d,
     &      core1_d
      integer ii_d,
     &         jj_d
      integer  lo1, up1, tmp, range
      integer  i, j,  svdlo, svdup
C***
C      Debug Code
      integer dny0
C***  

C***    Timing Code 
      real*8   etime,
     &         btime,
     &         itime
C***    Timing Code 

      real    facdlma
      real    facdlmb, r

C***    Timing Code 
      integer mp_timer_clear
      integer mp_timer_start
      integer mp_timer_stop
      real*8 mp_elapsed_time
      real*8 mp_busy_time
      real*8 mp_idle_time
C***    Timing Code 

 107  format ( I4, F16.8)
 
              range = up -lo + 1
C1              dny0 = ny0 * 100 + 1
C1              print *, '+WINA'
C1              if (ny0 .eq. 0) then
C1                 print *, '+WINA'
C1              endif
C1
C1              do j = 1, range
C1                 do i = 1, imp3
C1                    dny0 = 1000 + ny0 * 10 + i
C1                    print 107, dny0, wina_d(i,j)
C1                 end do
C1              end do
C1              call mp_sync()
C1
C1              if (ny0 .eq. 0) then
C1                 print *, '+RLAT'
C1              endif
C1              if (ny0 .eq. (ny-1)) then
C1                 range = jm - lo + 1
C1              endif
C1              do j = 1, range
C1                 do i = 1, imp3
C1                    dny0 = 2000 + ny0 * 10 + i
C1                    print 107, dny0, rlatitude_d(i,j)
C1                 end do
C1              end do
C1              call mp_sync()
C1              print *, '+WINA OVER'
C1              print *, '+WINB'
C1              print 101, ((winb(dk, dj), dj = 1, jmp3), dk = 1, imp3)
C1              print *, '+WINB OVER'
C1              print *, '+RLAT'
C1              print 102, ((rlatitude(dk, dj), dj = 1, jm), dk = 1, imp3)
C1              print *, '+RLAT OVER'
C1              print *, '+RLON'
C1              print 102, ((rlongitude(dk, dj), dj = 1, jm), 
C1     &             dk = 1, imp3)
C1              print *, '+RLON OVER'
C1              print 103, faclama, faclamb, facphia, facphib

 101  format( 1X, 8(4F16.8 //), 3F16.8 //)
 102  format( 1X, 8(4F16.8 //))
 103  format( 1X, '+', 'Faclama : ', F16.8, 'Faclamb : ', 
     &         F16.8, 'Facphia : ', F16.8, 'Facphib : ', F16.8 //)


      lo1 = lo
      up1 = up

      if (ny0 .eq. 0) then
         lo1 = ja
      else if (ny0 .eq. ny - 1) then
         up1 = jb
      endif

      range = up1 -lo1 + 1
C Debug         
D      print 201, ny0, lo1, up1
 201  format(3(I4))
C Debug
        call mp_init()


C1      allocate(tem1_d(1:imp3, 1:up-lo + 1))
C1      allocate (tem2_d(1:imp3, 1:up-lo + 1))
C1      allocate (ii_d(1:imp3, 1:up-lo + 1))
C1      allocate (jj_d(1:imp3, 1:up-lo + 1))
C1      allocate (core1_d(1:imp3, 1:up-lo + 1))


C***    Timing Code - Part 1
	tmp = mp_timer_clear(0)
        if (tmp .eq. -1) then
           print *, 'Failure in clearing timer'
C1           stop
        end if
        tmp = mp_timer_start(0)
        if (tmp .eq. -1) then
           print *, 'Failure in starting timer'
C1           stop
        end if
C***    Timing Code - Part 1


c r transforms radians into "turns" 
	r=one/twopi
c dlm is 2*pi/im
c u and v quantities are offset differently from the equator.
        facdlma = (2.-faclama)*dlm 
        facdlmb = (2.-faclamb)*dlm

C1	do j = lo1,up1    
        do j = 1, range 
	 do i = 1,imp3

c add enough turns as to have a positive lon
c and compute the index number of the grid point to the
c east of the lon 
          tem1_d(i,j) = (rlongitude_d(i,j)-facdlma)*r   
          ii_d(i,j) =  int(tem1_d(i,j))
C1          PRINT *, 'First time ii', ii_d(i,j)
c if started with a negative give a turn more to get into positive 
c angles
	 if(tem1_d(i,j) .lt. 0.) ii_d(i,j) = ii_d(i,j)-1
C1          PRINT *, 'Second Time ii', ii_d(i,j)
c get the actual gridpoint index in the lon direction
	   tem1_d(i,j) = (tem1_d(i,j) - ii_d(i,j))*im+two
           ii_d(i,j) = int(tem1_d(i,j))
C1          PRINT *, 'Third Time ii', ii_d(i,j)
           tem1_d(i,j) = tem1_d(i,j) - ii_d(i,j)
           if(ii_d(i,j).eq.1)ii_d(i,j) = imp1
C1          PRINT *, 'Third Time ii', ii_d(i,j)
c same thing for the latitudes 
c the computation for core1 can be used for the second field
         core1_d(i,j) = (  rlatitude_d(i,j) + pio2)*rdph 
         tem2_d(i,j) = core1_d(i,j) + facphia
         jj_d(i,j) = int(tem2_d(i,j))
         tem2_d(i,j) = tem2_d(i,j) - jj_d(i,j)
        enddo
C Debug
D        print 202, "For wina : ", ny0, j
 202    format(A, 2(I4))
C Debug
       enddo

        svdlo = lo
        svdup = up
        call circulate(wina_d, tem1_d,  tem2_d, 
     &                 woua_d, ii_d,    jj_d,
     &                 svdlo,     svdup,     
     &                 ja,     jb)
        
           
C1	do j = lo1,up1    
        do j = 1, range 
	 do i = 1,imp3
c
c
c the comments are the same as above
          tem1_d(i,j) = (rlongitude_d(i,j)-facdlmb)*r  
          ii_d(i,j) =  int(tem1_d(i,j))
	  if(tem1_d(i,j) .lt. 0.) ii_d(i,j) = ii_d(i,j) -1
	  tem1_d(i,j) = (tem1_d(i,j) - ii_d(i,j))*im+two
          ii_d(i,j) = int(tem1_d(i,j))
          tem1_d(i,j) = tem1_d(i,j) - ii_d(i,j)
          if(ii_d(i,j) .eq. 1) ii_d(i,j) = imp1
c same thing for the latitudes 
          tem2_d(i,j) = core1_d(i,j) + facphib
          jj_d(i,j) = int(tem2_d(i,j))
          tem2_d(i,j) = tem2_d(i,j) - jj_d(i,j)
         enddo
        enddo

        svdlo = lo
        svdup = up
        call circulate(winb_d, tem1_d,  tem2_d, 
     &                 woub_d, ii_d,    jj_d,
     &                 svdlo,     svdup,      
     &                 ja,     jb)


C***    Timing Code - Part 2
        
        tmp = mp_timer_stop(0)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
C1           stop
        end if

        etime = mp_elapsed_time(0)
        btime = mp_busy_time(0)
        itime = mp_idle_time(0)

C1      if (ny0 .eq. 15) then
C1        print 200, ">Node : ",ny0, " , Circul elapse time: ", etime
C1      endif
 200    format (A, I3, A, E22.15.3)

C1        print 200, ">Node : ",ny0, " , Circul elapse time: ", etime
C1        print 200, ">Node : ",ny0, " , Circul calcn. time: ", btime
C***    Timing Code - Part 2

C       End of routine bill_mp
      return
      end




