c @(#) bicubic.f 1.5 (IDS) 1/19/95

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


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

      implicit none

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

      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)
      integer  lower, 
     &          upper,
     &         ja, jb
      integer  jjp1, jjp2, jjp3, jjp4


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

      real win_n( imp3, jmp3)

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

      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***    Timing Code 
C      integer mp_timer_clear
C      integer mp_timer_start
C      integer mp_timer_stop
C      real*8 mp_elapsed_time
C      real*8 mp_busy_time
C      real*8 mp_idle_time
C***    Timing Code 

D      print *, "Starting"

      call mp_init() 

      call mp_sync()

      range = upper - lower + 1

      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

      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' )

      call get_lsd_bc_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 ?
C1             lindex = ny0 + 1 ?
             call comm_south(win_d, upper, lower, max_iter, 
     &                       win_n, wr_south, wr_lower, co_lower, 
     &                       wrapd, width)  

C1             comm_s = 0 ?
C1             lindex = ny0 + 1 ?
             call comm_north(win_d, upper, lower, max_iter, 
     &                       win_n, wr_north, width)  

             call readjust_jj(jj, wr_south, wr_north, width,
     &            wr_lower, co_lower, wrapd, njj)  

             if (ny0 .eq. 0) then
                print *, "Going to interpolate "
             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)

      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
         
      do j = lo, range
         do i = 1, imp3

            call find_bc_jj(njj(i,j), jjp1, jjp2, jjp3, 
     *                            jjp4, width)
            call  cubit( delx(i,j),dely(i,j), 
     &           win_n(ii(i,j)  ,jjp2),win_n(ii(i,j),jjp3  ),
     &           win_n(ii(i,j)  ,jjp1),win_n(ii(i,j),jjp4  ),
     &           win_n(ii(i,j)+1,jjp2),win_n(ii(i,j)+1,jjp3),
     &           win_n(ii(i,j)+1,jjp1),win_n(ii(i,j)+1,jjp4),
     &           win_n(ii(i,j)-1,jjp2),win_n(ii(i,j)-1,jjp3),
     &           win_n(ii(i,j)-1,jjp1),win_n(ii(i,j)-1,jjp4),
     &           win_n(ii(i,j)+2,jjp2),win_n(ii(i,j)+2,jjp3),
     &           win_n(ii(i,j)+2,jjp1),win_n(ii(i,j)+2,jjp4), 
     &           wou_d(i,j))
C1            print 108, i, j, ii(i,j), njj(i,j), 
C1     &                  delx(i,j), dely(i,j),
C1     &           win_n(ii(i,j)  ,jjp2),win_n(ii(i,j),jjp3  ),
C1     &           win_n(ii(i,j)  ,jjp1),win_n(ii(i,j),jjp4  ),
C1     &           win_n(ii(i,j)+1,jjp2),win_n(ii(i,j)+1,jjp3),
C1     &           win_n(ii(i,j)+1,jjp1),win_n(ii(i,j)+1,jjp4),
C1     &           win_n(ii(i,j)-1,jjp2),win_n(ii(i,j)-1,jjp3),
C1     &           win_n(ii(i,j)-1,jjp1),win_n(ii(i,j)-1,jjp4),
C1     &           win_n(ii(i,j)+2,jjp2),win_n(ii(i,j)+2,jjp3),
C1     &           win_n(ii(i,j)+2,jjp1),win_n(ii(i,j)+2,jjp4) 
C1                       
         end do
      end do

 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 locubit



      subroutine get_lsd_bc_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),
     &         jj(imp3, max_col_p3),
     &        ny, ny0, nysth

        integer mp_timer_clear
        integer mp_timer_start
        integer mp_timer_stop
C        real*8 mp_elapsed_time
C        real*8   etime
CT
         real*8 gcbtime, lcbtime, cabtime
         common /ntime/ cabtime, gcbtime, lcbtime

         real*8 mp_busy_time
CT

      implicit none

        integer lmax_iter(2),
     &           dist,
     &           ownp0, ownp4, i, j, tmp
  
	lmax_iter(1) = 0
        lmax_iter(2) = 0
        dist      = 0
  
          do j = lo1, up1
           do i = 1, imp3
              ownp0 = int((jj(i,j)) / width)
              ownp4 = int((jj(i,j)+4) / width)

              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)+4, width) .eq. 0) then
                 ownp4 = ownp4 -1
              else if (ownp4 .ge. ny) then
                 ownp4 = ny -1
              endif

C
C	      Shift address of current node ny0 to 0.  
C             Thus ownp0, ownp4 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
		ownp4 = ownp4 - ny0
		ownp0 = ownp0 - ny0

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

                if (ownp0 .gt. ny / 2) then
    		  dist = ny - ownp0
		  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 (ownp4 .gt. ny / 2) then
    		  dist = ny - ownp4
		  if (lmax_iter(2) .lt. dist) lmax_iter(2) = dist  
                else 
		  dist = ownp4
		  if (lmax_iter(1) .lt. dist) lmax_iter(1) = dist
		endif

                if (lmax_iter(1) .gt. 5 ) then
                   print *, ">Nor NY0 : ", ny0, " JJ ", jj(i,j), " J ",
     &                  j,  " width : ", width, " iter ", lmax_iter(1)
                end if
                if (lmax_iter(2) .gt. 5 ) then
                   print *, ">Sou NY0 : ", ny0, " JJ ", jj(i,j), " J ",
     &                  j,  " width : ", width, " iter ", lmax_iter(2)
                end if
           enddo
        enddo

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) 
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
 
         if (ny0 .eq. 0) then
            print *, "> North Iter", max_iter(1), "  South Iter", 
     &      max_iter(2)
         endif


      end

C** end of subroutine get_lsd_bc_rng



      subroutine find_bc_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_bc_jj
