c %Z% %M% %I% (IDS) %G%
C+
C*
C*   This file contains the following routines in the nodal part 
C*   of the host-node program :
C*   
C*   1) subroutine lobill - Performs a bilinear interpolation 
C*   2) subroutine get_lsd_bl_rng - Finds the degree of circulation
C*   3) subroutine find_bc_jj - not used 
C*                              
C-


C+ 
C  Subroutine - lobill
C  This routine first calls get_lsd_bl_rng to find the degree of
C  circulation in the North-South direction.  It then calls              
C  comm_south and comm_north to build a segment. Further details
C  are given in SCCS 694: Case Studies on Parallelization of NASA's 
C  Weather Code : Results on Interpolation Techniques. It then 
C  generates the set of local indices. Each node then performs
C  a local bilinear inerpolation The other major difference
C   is all the points of interpolation of a component are first
C   computed and stored in an array before the actual interpolation 
C   unlike the sequential
C   code, where the interpolation is performed after computing 
C   a single point. See also bill.f (original code)
C
C
C  Inputs : 
C
C           win_d         - decomposed wind component
C           delx          - relative x dist from origin of grid square.
C           dely          - relative y dist from origin of grid square.
C           ja, jb        - range of latitude indices
C           lower, upper  - partition bounds on each node
C           ii, jj        - longitude, latitude indices
C  
C  Outputs : wou_d        - interpolated results
C  
C     
C-


       subroutine lobill(win_d, delx, dely,  wou_d, ii,
     &                    jj, lower, upper, ja, jb)


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

      implicit none

C
C     Common block which contains information about the logical
C     layout of nodes for the algorithm.
C
      integer nx, nx0, ny, ny0, nynth, nysth, nxest, nxwst
      common /commvars/ nx,nx0,ny,ny0, nynth, nysth, nxest, nxwst

C
C     Decomposed arrays.
C
      real      win_d(imp3,max_col_p3),
     &          delx(imp3,max_col_p3),
     &          dely(imp3,max_col_p3),
     &          wou_d(imp3,max_col_p3)

      integer  ii(imp3,max_col_p3), 
     &          jj(imp3,max_col_p3)

C
C    Partition bounds - lower, upper.  ja, jb are limits on
C    latitude indices.
C
      integer  lower, 
     &          upper,
     &         ja, jb
      integer  jjp1, jjp2, jjp3, jjp4


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


C     Segment Buffer
      real win_n( imp3, jmp3)

C
C     Temporary locations
      integer i, j, k 
C1    ,l, m
      integer max_iter(2), range
      integer lo, up
C1      integer st, fi, rem, lim
      integer dk, dj

C
C     Variables used in generation of local indices.
C
      integer wr_south, wr_north, width
      integer wr_lower, co_lower, wrapd
      integer njj(imp3, max_col_p3)

C****
C***  Debug Variables
C1      integer dlower(6), dupper(6), lindex, comm_s, dwin(4,30) ?
C1      common /dvars/ dlower, dupper, lindex, comm_s, dwin

C     Bilinear interpolation function prototype      
      real bli

         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

         integer tmp

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

D      print *, "Starting"

      call mp_init() 

      call mp_sync()

      range = upper - lower + 1

C
C     Determination of latitude indices for computation
C
      lo = 1
      up = range
      if (ny0 .eq. 0) then
         lo = ja
      else if (ny0 .eq. ny - 1) then
CT         up = jb
         up = jb - lower + 1
      endif

C
C     Initialization of segment - need it only for debugging.
C
      do k = 1, imp3
         do j = 1, jmp3
            win_n(k, j) = -1
         end do
      end do

C1      open (unit = 13, file = 'debug.out', status = 'scratch' )
C1      close(13)
C1      open (unit = 13, file = 'debug.out', status = 'new' )


C
C     Determination of range of circulation.  Please see SCCS 694
C     for detailed explanation.  max_iter(1) - has the maximum
C     inter-nodal distance in northern direction and max_iter(2)
C     has for the southern.
C
      call get_lsd_bl_rng(max_col, lo, up, jj, max_iter,
     *                        ny, ny0, nysth)


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


C1             comm_s = 1           ! Not used anymore
C1             lindex = ny0 + 1     ! Not used anymore
C
C       Build the segment in the southern direction.       
C
             call comm_south(win_d, upper, lower, max_iter, 
     &                       win_n, wr_south, wr_lower, co_lower, 
     &                       wrapd, width)  

C1             comm_s = 0           ! Not used anymore
C1             lindex = ny0 + 1     ! Not used anymore

C
C       Build the segment in the northern direction.       
C
             call comm_north(win_d, upper, lower, max_iter, 
     &                       win_n, wr_north, width)  

C
C       Generate local indices
C
             call readjust_jj(jj, wr_south, wr_north, width,
     &            wr_lower, co_lower, wrapd, njj)  

CD             if (ny0 .eq. 0) then
CD                print *, "Going to interpolate "
CD             endif
109      format ( 4F8.4 )
C1             print *, "Node", ny0
C1             print 109, ((win_n(k,j), k = 1, imp3), j = 1, jmp3)

C1          write (13, 103) ny0, wr_south, wr_north
C1          write (13, 104) width, wr_lower, co_lower, wrapd 
C1          write (13, 105) ((jj(k,j), k = 1, imp3), j = 1, range)
C1          write (13, 106) ((njj(k,j), k = 1, imp3), j = 1, range)
C1          close(13)

C
C     Determine index ranges for interpolation.
C
      lo = 1

      if (ny0 .eq. 0) then
         range = upper - lower + 1
         lo = ja
      else if (ny0 .eq. ny - 1) then
         range = jb - lower + 1
      endif
         

CT    Timing Code for measuring bli overhead
            tmp = mp_timer_clear(5)
            if (tmp .eq. -1) then
               print *, 'Failure in clearing timer'
               stop
            end if
         
C
C     Interpolate
C
      do j = lo, range
         do i = 1, imp3

CC            call find_bl_jj(njj(i,j), jjp1, jjp2, jjp3, 
CC     *                            jjp4, width)

            wou_d(i,j) =  bli ( delx(i,j) , dely(i,j),
     &                win_n(ii(i,j)  ,njj(i,j)+2), 
     &                win_n(ii(i,j) , njj(i,j)+3),
     &                win_n(ii(i,j)+1,njj(i,j)+2), 
     &                win_n(ii(i,j)+1,njj(i,j)+3) )

         end do
      end do

C     
C     Accumulate cubit overhead.
C     
            intime = intime + mp_busy_time(5) 

 101  format( 1X, '+', 4F16.8 //, '+' 3F16.8)
C1
C1       print *, '+WOU'
C1       print 101, ((wou_d(dk, dj), dj = 1, max_col_p3), dk = 1, imp3)
C1       print *, '+WOU OVER'

C***    Timing Code 
C        btime = mp_busy_time(1)
C        etime = mp_elapsed_time(1)

C       print 200, ">Node : ", ny0, " , Comm. Time:",
C     &             etime
C       print 200, ">Node : ", ny0, " , Comm. Time:",
C     &             btime
C 200   format (A, I4, A, E22.15.3)
C***    Timing Code 


C1 102      format( 4I4)
C1 103      format ( "Ny0", I4, "South Wrap", I4, "North Wrap", I4)
C1 104      format ( "Wdth", I4, "Wr Lr", I4, "Co Lr", I4,  "Wrd", I4 )     
C1 105      format ( "JJ" //, 4(I4))
C1 106      format ( "NJJ" //, 4(I4))
C1 108      format ( "I", I4, " J", I4, " II", I4, " JJ", I4 //, 
C1     &             "DELX", F16.8, " DELY", F16.8 //,
C1     &             "W1", F16.8, "W2", F16.8, "W3", F16.8, 
C1     &             "W4", F16.8 //,
C1     &             "W5", F16.8, "W6", F16.8, "W7", F16.8, 
C1     &             "W8", F16.8 //,
C1     &             "W9", F16.8, "W10", F16.8, "W11", F16.8, 
C1     &             "W12", F16.8 //,
C1     &             "W13", F16.8, "W14", F16.8, "W15", F16.8, 
C1     &             "W16", F16.8)

      end

C**  end of subroutine lobill



C+ 
C  Subroutine - get_lsd_bl_rng
C  This routine first determines the owners of the two extreme  
C  points for a given value of jj(i,j).  One of the owners is            
C  assumed to have jj(i,j) and the other jj(i,j)+3.  Then it   
C  establishes whether the owner can be reached in a shorter hop    
C  in the northern direction (or southern).  Then the current  
C  northern (or southern) maximum is updated if the current hop
C  is greater than the current maximum.                     
C
C  Inputs : 
C
C           width         - partition width
C           lo1           - lower bound for computation (not partition)
C           up1           - upper  "     "     "          "    "
C           jj            - computed latitude indices
C           ny            - Number of nodes used for the N-S dimension
C           ny0           - address of node
C           nysth         - address of southern neighbor
C  
C  Outputs : max_iter(2)   - Northern & Southern Maxima  
C                            Note : Index value = 1 is used for Northern
C                                   Index value = 2 "   "    "  Southern
C     
C-
      subroutine get_lsd_bl_rng(width, lo1, up1, jj, max_iter,
     *                            ny, ny0, nysth)

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

      integer lo1,
     &        up1,
     &        width,
     &        max_iter(2),
c
c    !!! Change Dimensions of jj to jj(imp3, max_col_p3)
c
C     &        jj(4,4),         ! used in sequential simulation
     &         jj(imp3, max_col_p3),
     &        ny, ny0, nysth


CT
CT      Function Prototypes for Timing Code
        integer mp_timer_clear
        integer mp_timer_start
        integer mp_timer_stop
C        real*8 mp_elapsed_time
C        real*8   etime
        real*8 mp_busy_time
CT
D         real*8 gcbtime, lcbtime, cabtime
D         common /ntime/ cabtime, gcbtime, lcbtime
         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
CT
CT       For recording avg number of iterations.
         real niter, siter
         common /citer/ niter, siter

      implicit none

        integer lmax_iter(2),               ! local maxima 1 - North, 2 - South
     &           dist,                      ! distance
     &           ownp0, ownp3,              ! owner nodes
     &           i, j, tmp                  ! counter, status
  
CT    Timing Code for measuring calc. for degree of circulation.
	tmp = mp_timer_clear(6)
        if (tmp .eq. -1) then
           print *, 'Failure in clearing timer'
           stop
        end if
        tmp = mp_timer_start(6)
        if (tmp .eq. -1) then
           print *, 'Failure in starting timer'
           stop
        end if
CT    Timing Code 

C
C       Initialize Variables
C
	lmax_iter(1) = 0
        lmax_iter(2) = 0
        dist      = 0
  
          do j = lo1, up1
           do i = 1, imp3
C
C             Find owner nodes 
C
              ownp0 = int((jj(i,j)) / width)
CC              ownp0 = int((jj(i,j)+2) / width)
              ownp3 = int((jj(i,j)+3) / width)

C
C             If jj(i,j) is exactly divisible by partition width,
C             or owner address greater than the maximum valid address
C             adjust owner address.
C
              tmp = mod(jj(i,j), width)
              if (tmp .eq. 0) then
                 ownp0 = ownp0 -1
              else if (ownp0 .ge. ny) then
                 ownp0 = ny -1
              endif
                 
              if (mod(jj(i,j)+3, width) .eq. 0) then
                 ownp3 = ownp3 -1
              else if (ownp3 .ge. ny) then
                 ownp3 = ny -1
              endif

C
C	      Shift address of current node ny0 to 0.  
C             Thus ownp0, ownp3 will be shifted accordingly.
C             Then try to find the shortest path to the further
C             of the two in either the southerly or northerly 
C             direction.
C
		ownp3 = ownp3 - ny0
		ownp0 = ownp0 - ny0

		if (ownp3 .lt. 0) ownp3 = ny + ownp3
		if (ownp0 .lt. 0) ownp0 = ny + ownp0

                if (ownp0 .gt. ny / 2) then
C
C     The shortest distance is by the southern direction, as the 
C     inter-nodal distance is greater than half the number of 
C     nodes in the northern direction.
C
    		  dist = ny - ownp0
C
C     Update the local southern maximum.
C
		  if (lmax_iter(2) .lt. dist) lmax_iter(2) = dist  
                else 
		  dist = ownp0
		  if (lmax_iter(1) .lt. dist) lmax_iter(1) = dist
		endif

                if (ownp3 .gt. ny / 2) then
    		  dist = ny - ownp3
		  if (lmax_iter(2) .lt. dist) lmax_iter(2) = dist  
                else 
		  dist = ownp3
		  if (lmax_iter(1) .lt. dist) lmax_iter(1) = dist
		endif

C
C     Debug code.  The number of hops should not exceed 5 for 
C     this data-set in either direction.
C
D                if (lmax_iter(1) .gt. 5 ) then
D                   print *, ">Nor NY0 : ", ny0, " JJ ", jj(i,j), " J ",
D     &                  j,  " width : ", width, " iter ", lmax_iter(1)
D                end if
D                if (lmax_iter(2) .gt. 5 ) then
D                   print *, ">Sou NY0 : ", ny0, " JJ ", jj(i,j), " J ",
D     &                  j,  " width : ", width, " iter ", lmax_iter(2)
D                end if
           enddo
        enddo

CT    Timing Code 
        tmp = mp_timer_stop(6)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
           stop
        end if
C
C       Accumulate times for measuring calc. for degree of circulation.
C
        brtime = brtime + mp_busy_time(6) 
CT    Timing Code 

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

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

        max_iter(1) = cmmd_reduce_int(lmax_iter(1), CMMD_combiner_max) 
        max_iter(2) = cmmd_reduce_int(lmax_iter(2), CMMD_combiner_max) 

C
C     Record maxima
C
        niter = niter + max_iter(1)
        siter = siter + max_iter(2)
        
C
C     Debug code
C
C1        max_iter(1) = lmax_iter(1) 
C1        max_iter(2) = lmax_iter(2) 

C***    Timing Code 
        tmp = mp_timer_stop(0)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
           stop
        end if
        gcbtime = gcbtime + mp_busy_time(0)
C
C         print 200, ">Node : ", ny0, " , Global Comm. Time: ",
C     &             etime
 200    format (A, I4, A, E22.15.3)

C***    Timing Code 
D        print *, "L Max Iter is : ", lmax_iter(1), " (N), ", 
D     &            lmax_iter(2), " (S)", " Node : ", ny0
 
D         if (ny0 .eq. 0) then
D            print *, "> North Iter", max_iter(1), "  South Iter", 
D     &      max_iter(2)
D         endif


      end

C** end of subroutine get_lsd_bl_rng


C+ 
C  Subroutine - find_bl_jj
C  This routine is not used
C-
      subroutine find_bl_jj(jj, jjp1, jjp2, jjp3, 
     *                      jjp4, width)

        integer jj, jjp1, jjp2
        integer jjp3, jjp4, width

        jjp1 = jj + 1
        if (jjp1 .gt. width) then
           jjp1 = jjp1 - width
        endif

        jjp2 = jj + 2
        if (jjp2 .gt. width) then
           jjp2 = jjp2 - width
        endif

        jjp3 = jj + 3
        if (jjp3 .gt. width) then
           jjp3 = jjp3 - width
        endif

        jjp4 = jj + 4
        if (jjp4 .gt. width) then
           jjp4 = jjp4 - width
        endif

      end

C** end of subroutine find_bl_jj

C+ 
C  Function - bli
C  This function performs the bilinear interpolation and returns
C  the value.  Please refer to SCCS 694 for further details.
C
C  Inputs : 
C
C           delx, dely    - relative distance from lower right hand corner
C                           of the grid-square in x, y axes
C           w1, w2, w3, w4 - values at the four corners of the grid-square
C  
C  Outputs : 
C           None 
C
C  Return Value : 
C           bli - interpolated value
C-
      function bli(delx,dely,w1,w2,w3,w4) 

      implicit none
      
      real bli           
      
      real delx, dely
      real w1, w2, w3, w4

C***    Timing Code 
      integer mp_timer_start
      integer mp_timer_stop
      integer tmp
C***    Timing Code 

CT
            tmp = mp_timer_start(5)
            if (tmp .eq. -1) then
               print *, 'Failure in starting timer'
               stop
            end if
CT
        bli = (1.-delx)*((1.-dely)*w1+dely*w2) +
     >        delx *((1.-dely)*w3+dely*w4)

CT    Timing Code 
            tmp = mp_timer_stop(5)
            if (tmp .eq. -1) then
               print *, 'Failure in stopping timer'
               stop
            end if
CT

      end
