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 bill_mp - Calculates the indices at which interpolation
C*                          is desired.
C-

C+ 
C  Subroutine - bill_mp 
C  This routine computes for each of the two components, the points   
C  where interpolation is desired. After computing the points of         
C  interpolation for the first component, the subroutine calls 
C  'lobill' to perform the interpolation. The same is repeated
C   for the next component.  The 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 : wina_d corresponds to ust
C           winb_d    "        "  vst
C           rlatitude_d "      "  spuh
C           rlongitude_d "     "  sluh
C           faclama   
C           faclamb    
C           facphia
C           facphib 
C           ja, jb      - range of latitude indices
C           lo, up,     - partition bounds on each node
C  
C  Outputs : woua corresponds to c1a
C            woub     "       "  c2a
C  
C     
C-

      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"
      include "slgcm.mp"

      implicit none

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

C***  Debug Variables
C1          integer comm_a
C1          common /switch/ comm_a

C
C     Timing Code
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

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


       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)          ! temporary buffers
      dimension tem2_d(imp3, max_col_p3)          !      "
      dimension core1_d(imp3, max_col_p3)         !      "
      dimension ii_d(imp3, max_col_p3)            ! longitude indices
      dimension jj_d(imp3, max_col_p3)            ! latitude indices

      real tem1_d,
     &      tem2_d,
     &      core1_d
      integer ii_d,
     &         jj_d
      integer  lo1, up1, tmp, range               ! computational range
      integer  i, j,  svdlo, svdup                ! temporary storage

CD
      integer dny0, drange
CD
C***    Timing Code 
C1      real*8   etime,
C1     &         btime,
C1     &         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 


CD
 108  format ( I10, F16.8)
 109  format ( I10, I4)
      
D      print *, "In bill mp"
C1              drange = up -lo + 1
C1              if (ny0 .eq. 0) then
C1                 print *, '+WINA'
C1              endif
C1
C1              do j = 1, drange
C1                 do i = 1, imp3
C1                    dny0 = 100000 + ny0 * 1000 + j * 10 + i
C1                    print 108, 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                 drange = jm - lo + 1
C1              endif
C1              do j = 1, drange
C1                 do i = 1, imp3
C1                    dny0 = 100000 + ny0 * 1000 + j * 10 + i
C1                    print 108, dny0, rlatitude_d(i,j)
C1                 end do
C1              end do
C1              call mp_sync()

CT    Timing Code for measuring ii(i,j), jj(i,j) calculation 
	tmp = mp_timer_clear(3)
        if (tmp .eq. -1) then
           print *, 'Failure in clearing timer'
           stop
        end if
        tmp = mp_timer_start(3)
        if (tmp .eq. -1) then
           print *, 'Failure in starting timer'
           stop
        end if
CT    Timing Code 

C
C     Find the local indices which form the bound for
C     iterations on each node
C

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

      range = up1 - lo + 1
      

C1      call mp_init()

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

CD
CD        print '(I4, "Lo", I4, " Up", I4, " Lo1", I4, " Up1", 
CD     &          I4, " Range"
CD     &          I4, " Ja", I4, " Jb", I4)', ny0, lo, up, lo1, up1,
CD     &          range, ja, jb
CD

C        do j = 1, range 
        do j = lo1, 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 
CD
CD          if (ny0 .eq. 1) then
CD             print '(" ITER : J ", I4, " I ", I4)', j, i
CD             print '(" rlon", F16.8, " facdlma", F4.2,
CD     &       " r", F16.8)', rlongitude_d(i,j), facdlma, r
CD          endif
CD
          tem1_d(i,j) = (rlongitude_d(i,j)-facdlma)*r   
          ii_d(i,j) =  int(tem1_d(i,j))
CD
CD          if (ny0 .eq. 1) print '(" tem1", F16.8, " II", I4)',
CD     &      tem1_d(i,j), ii_d(i,j)
CD
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
CD
CD         if (ny0 .eq. 1) PRINT *, 'Second Time ii ', ii_d(i,j)
CD
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))
CD
CD          if (ny0 .eq. 1) print '(" tem1", F16.8, " II ", I4)',
CD     &      tem1_d(i,j), ii_d(i,j)
CD
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)
CD
CD          if (ny0 .eq. 1) print '(" tem1", F16.8, " II", I4)',
CD     &      tem1_d(i,j), ii_d(i,j)
CD
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
       enddo

        svdlo = lo
        svdup = up

CT    Timing Code 
        tmp = mp_timer_stop(3)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
           stop
        end if
C
C       Accumulate ii(i,j), jj(i,j) calculation times for bicubic.
C
        llitime = llitime + mp_busy_time(3) 
CT    Timing Code 

C1        call circulate(wina_d, tem1_d,  tem2_d, 
C1     &                 woua_d, ii_d,    jj_d,
C1     &                 svdlo,     svdup,     
C1     &                 ja,     jb)
  
C*** Debug Code      
CD        comm_a  = 1

CD              drange = up -lo + 1
CD              if (ny0 .eq. 0) then
CD                 print *, '+II'
CD              endif
CD              call mp_sync()
CD
CD              do j = lo1, range
CD                 do i = 1, imp3
CD                    dny0 = 100000 + ny0 * 1000 + j * 10 + i
CD                    print 109, dny0, ii_d(i,j)
CD                 end do
CD              end do
CD              call mp_sync()
CD
CD              if (ny0 .eq. 0) then
CD                 print *, '+JJ'
CD              endif
CD              call mp_sync()
CD
CD              do j = lo1, range
CD                 do i = 1, imp3
CD                    dny0 = 200000 + ny0 * 1000 + j * 10 + i
CD                    print 109, dny0, jj_d(i,j)
CD                 end do
CD              end do
CD              call mp_sync()
CD
CD              if (ny0 .eq. 0) then
CD                 print *, '+TEM1D'
CD              endif
CD
CD              do j = lo1, range
CD                 do i = 1, imp3
CD                    dny0 = 300000 + ny0 * 1000 + j * 10 + i
CD                    print 108, dny0, tem1_d(i,j)
CD                 end do
CD              end do
CD
CD              call mp_sync()
CD              if (ny0 .eq. 0) then
CD                 print *, '+TEM2D'
CD              endif
CD
CD              do j = lo1, range
CD                 do i = 1, imp3
CD                    dny0 = 400000 + ny0 * 1000 + j * 10 + i
CD                    print 108, dny0, tem2_d(i,j)
CD                 end do
CD              end do
CD
CD              call mp_sync()


C
C       Perform interpolation 
C
        call lobill(wina_d, tem1_d, tem2_d,  woua_d,
     &               ii_d,  jj_d, svdlo, svdup, ja, jb)



CD        print '(I4, "Lo", I4, " Up", I4, " Lo1", I4, " Up1", 
CD     &          I4, " Range",
CD     &          I4, " Ja", I4, " Jb", I4)', ny0, lo, up, lo1, up1,
CD     &          range, ja, jb


CT    Timing Code for measuring ii(i,j), jj(i,j) calculations
	tmp = mp_timer_clear(3)
        if (tmp .eq. -1) then
           print *, 'Failure in clearing timer'
           stop
        end if
        tmp = mp_timer_start(3)
        if (tmp .eq. -1) then
           print *, 'Failure in starting timer'
           stop
        end if
CT    Timing Code 

           
C1        do j = 1, range 
        do j = lo1, range 
	 do i = 1,imp3
c
c
c the comments are the same as above

CD          if (ny0 .eq. 1) then
CD             print '(" ITER : J ", I4, " I ", I4)', j, i
CD             print '(" rlon", F16.8, " facdlma", F4.2,
CD     &       " r", F16.8)', rlongitude_d(i,j), facdlma, r
CD          endif

          tem1_d(i,j) = (rlongitude_d(i,j)-facdlmb)*r  
          ii_d(i,j) =  int(tem1_d(i,j))

CD
CD          if (ny0 .eq. 1) print '(" tem1", F16.8, " II", I4)',
CD     &      tem1_d(i,j), ii_d(i,j)
CD
	  if(tem1_d(i,j) .lt. 0.) ii_d(i,j) = ii_d(i,j) -1

CD          if (ny0 .eq. 1) PRINT *, 'Second Time ii ', ii_d(i,j)

	  tem1_d(i,j) = (tem1_d(i,j) - ii_d(i,j))*im+two
          ii_d(i,j) = int(tem1_d(i,j))

CD          if (ny0 .eq. 1) print '(" tem1", F16.8, " II ", I4)',
CD     &      tem1_d(i,j), 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


CD          if (ny0 .eq. 1) print '(" tem1", F16.8, " II", I4)',
CD     &      tem1_d(i,j), ii_d(i,j)

c same thing for the latitudes 
CD
CD          if (ny0 .eq. 1) then
CD              print '(" rlat", F16.8, " facphib", F16.8)',
CD     &        rlatitude_d(i,j), facphib
CD              print '(" rdph", F16.8, " pio2", F16.8)', rdph, pio2
CD            core1 = (rlatitude(i,j) + pio2) * rdph
CD            print '(" core1", F16.8)', core1_d(i,j)
CD         endif
CD
          tem2_d(i,j) = core1_d(i,j) + facphib
          jj_d(i,j) = int(tem2_d(i,j))
CD
CD          if (ny0 .eq. 1) print '(" tem2", F16.8, 
CD     &     " JJ ", I4)', tem2_d(i,j), jj_d(i,j)
CD
          tem2_d(i,j) = tem2_d(i,j) - jj_d(i,j)

CD
CD          if (ny0 .eq. 1) print '(" tem2", F16.8)',
CD     &     tem2_d(i,j)
CD
         enddo
        enddo

        svdlo = lo
        svdup = up

CT    Timing Code 
        tmp = mp_timer_stop(3)
        if (tmp .eq. -1) then
           print *, 'Failure in stopping timer'
           stop
        end if
C
C       Accumulate ii(i,j), jj(i,j) calculation times for bicubic.
C
        llitime = llitime + mp_busy_time(3) 
CT    Timing Code 

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

C*** Debug Code      
CD        comm_a = 0

C
C       Perform interpolation 
C
        call mp_sync()
        call lobill(winb_d, tem1_d, tem2_d,  woub_d,
     &               ii_d,  jj_d, svdlo, svdup, ja, jb)
           

 107    format ( I10, F16.8, F16.8)
C1
CD             range = up - lo + 1
CD              if (ny0 .eq. 0) then
CD                 print *, '+WOU'
CD              endif
CD              if (ny0 .eq. (ny-1)) then
CD                 range = jm - lo + 1
CD              endif
CD              do j = 1, range
CD                 do i = 1, imp3
CD                    dny0 = 100000 + ny0 * 1000 + j * 10 + i
CD                    print 107, dny0, woua_d(i,j), woub_d(i,j)
CD                 end do
CD              end do
CD              call mp_sync()

C***    Timing Code - Part 2
        
C1        tmp = mp_timer_stop(0)
C1        if (tmp .eq. -1) then
C1           print *, 'Failure in stopping timer'
C1           stop
C1        end if
C1
C1        etime = mp_elapsed_time(0)
C1        btime = mp_busy_time(0)
C1        itime = mp_idle_time(0)
C1
C1      if (ny0 .eq. 15) then
C1        print 200, ">Node : ",ny0, " , Circul elapse time: ", etime
C1      endif
C1 200    format (A, I3, A, E22.15.3)
C1
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

