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

C+ 
C  Subroutine - bic_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  'locubit' to perform the interpolation. The same is repeated
C   for the next component.  This code is essentially the same
C   as the sequential version, except the wrap subroutine is
C   called in the host program.  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 bic1nu.f (original code)
C
C
C  Inputs : wina_d corresponds to yuupt
C           winb_d    "        "  yvvpt
C           rlatitude_d "      "  sphi
C           rlongitude_d "     "  slam
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 yuupts or yvvpts
C            woub     "       "  yvupts or yuvpts
C  
C     
C-

      subroutine bic_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
CT    Timing Code 
C1      real*8   etime,
C1     &         btime,
C1     &         itime
CT    Timing Code 

      real    facdlma
      real    facdlmb, r

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


CD
 108  format ( I10, F16.8)
 109  format ( I10, I4)

D      print *, "I am in subroutine bic mp"
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


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

C
C       Perform interpolation 
C
D        print *, "Calling Locubit"
        call locubit(wina_d, tem1_d, tem2_d,  woua_d,
     &               ii_d,  jj_d, svdlo, svdup, 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 
    
        do j = lo1, 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

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 

C*** Debug Code      
CD        comm_a = 0
C
C       Perform interpolation 
C
        call mp_sync()
D        print *, "Calling Locubit", ny0
        call locubit(winb_d, tem1_d, tem2_d,  woub_d,
     &               ii_d,  jj_d, svdlo, svdup, ja, jb)
           

 107    format ( I10, F16.8, F16.8)


C       End of routine bic_mp
      return
      end




