c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE pathwbc(ip,jp,kp,ij,kco0,kco) 

c     ******************************************************************
c
c     This routine determines the zonal and meridional faces through 
c     which the western boundary current flows 
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Argument variables

      integer	ip,jp,kp,ij,kco0,kco


c---- Local variables

      integer	n

      character*12 stringx




c     -----------------------------------------
c     Determine WBC zonal and meridional faces
c
c     Note
c     ----
c     1) The algorithm is based on a Lagrangian 
c        approach whereby the path of a fluid
c        particle along the western boundary is 
c        determined
c     2) The particle always tries to `stick' 
c        on the western boundary by first
c        looking at the possibility to flow 
c        right of its trajectory 
c     ------------------------------------------      

c---- Initialization

      do n=1,n_o
         i_obs(n) = 0
      enddo


c---- First face (prescribed in arguments) 

c     Contraint #

      kco = kco0 


c     WBC increases west of 40W in NH

      if (lat(jp)    .gt.0. and.
     >    lon(ip)*deg.gt.-40.) then
          b_wbc0 = b_wbc0_up
      else
          b_wbc0 = b_wbc0_do
      endif


c     Corresponding element in b of Ax=b

      if (ij.lt.0) b_wbc(kco) = -b_wbc0
      if (ij.gt.0) b_wbc(kco) = +b_wbc0


c     i_obs  () facilitates subsequent definition of A of Ax=b
c     str_wbc() facilitates subsequent definition of x (cf. x_prior.f)

      if (abs(ij).eq.1) then
        do n=1,n_o
           if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco
        enddo
        str_wbc(kco) = stringx('UU',ip,jp,kp)
      else
        do n=1,n_o
           if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco
        enddo
        str_wbc(kco) = stringx('VV',ip,jp,kp)
      endif 


c---- Subsequent faces: Start of loop

100   continue 


c---- Flag flag_vor(,,)

      if (iabs(ij).eq.1) then
        flag_vor(ip-1,jp,kp) = .false.
        flag_vor(ip  ,jp,kp) = .false.
      endif

      if (iabs(ij).eq.2) then
        flag_vor(ip,jp-1,kp) = .false.
        flag_vor(ip,jp  ,kp) = .false.
      endif


c---- WBC increases west of 40W in NH

      if (lat(jp)    .gt.0. and.
     >    lon(ip)*deg.gt.-40.) then
          b_wbc0 = b_wbc0_up
      else
          b_wbc0 = b_wbc0_do
      endif


c---- Particle path 

c--
      if (ij.eq.+1) then	! Particle is on zonal face and moves east
c--

      if (flag_V(ip,jp,kp)) then 

         ip         = ip
         jp         = jp
         ij         = -2
         kco        = kco + 1
         b_wbc(kco) = -b_wbc0 
         do n=1,n_o
            if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco 
         enddo 
         str_wbc(kco) = stringx('VV',ip,jp,kp)

      else if (flag_U(ip+1,jp,kp)) then

         ip         = ip + 1
         jp         = jp
         ij         = +1 
         kco        = kco + 1
         b_wbc(kco) = b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('UU',ip,jp,kp)

      else if (flag_V(ip,jp+1,kp)) then

         ip         = ip 
         jp         = jp + 1
         ij         = +2
         kco        = kco + 1
         b_wbc(kco) = b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('VV',ip,jp,kp)

      else

         b_wbc(kco) = 0.d0 
         ij         = -1

      endif

c--
      else if (ij.eq.-1) then	! Particle is on zonal face and moves west
c--

      if (flag_V(ip-1,jp+1,kp)) then

         ip         = ip - 1
         jp         = jp + 1
         ij         = +2
         kco        = kco + 1
         b_wbc(kco) = b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('VV',ip,jp,kp)

      else if (flag_U(ip-1,jp,kp)) then

         ip         = ip - 1
         jp         = jp
         ij         = -1
         kco        = kco + 1
         b_wbc(kco) = -b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('UU',ip,jp,kp)

      else if (flag_V(ip-1,jp,kp)) then

         ip         = ip - 1 
         jp         = jp 
         ij         = -2
         kco        = kco + 1
         b_wbc(kco) = -b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco)= stringx('VV',ip,jp,kp)

      else

         b_wbc(kco) = 0.d0
         ij         = +1

      endif

c--
      else if (ij.eq.+2) then	! Particle is on meridional face and moves north
c--

      if (flag_U(ip+1,jp,kp)) then

         ip         = ip + 1
         jp         = jp 
         ij         = +1
         kco        = kco + 1
         b_wbc(kco) = b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('UU',ip,jp,kp)

      else if (flag_V(ip,jp+1,kp)) then

         ip         = ip 
         jp         = jp + 1
         ij         = +2 
         kco        = kco + 1
         b_wbc(kco) = b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('VV',ip,jp,kp)

      else if (flag_U(ip,jp,kp)) then

         ip         = ip 
         jp         = jp 
         ij         = -1
         kco        = kco + 1
         b_wbc(kco) = -b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco  
         enddo
         str_wbc(kco) = stringx('UU',ip,jp,kp)

      else

         b_wbc(kco) = 0.d0
         ij         = -2 

      endif

c--
      else if (ij.eq.-2) then   ! Particle is on meridional face and moves south 
c--

      if (flag_U(ip,jp-1,kp)) then

         ip         = ip 
         jp         = jp - 1 
         ij         = -1
         kco        = kco + 1
         b_wbc(kco) = -b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('UU',ip,jp,kp)

      else if (flag_V(ip,jp-1,kp)) then

         ip         = ip 
         jp         = jp - 1
         ij         = -2
         kco        = kco + 1
         b_wbc(kco) = -b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('VV',ip,jp,kp)) i_obs(n) = kco 
         enddo
         str_wbc(kco) = stringx('VV',ip,jp,kp)

      else if (flag_U(ip+1,jp-1,kp)) then

         ip         = ip + 1 
         jp         = jp - 1 
         ij         = +1
         kco        = kco + 1
         b_wbc(kco) = b_wbc0
         do n=1,n_o
            if (string(n).eq.stringx('UU',ip,jp,kp)) i_obs(n) = kco 
         enddo
        str_wbc(kco) = stringx('UU',ip,jp,kp)

      else

         b_wbc(kco) = 0.d0
         ij         = +2

      endif

c--
      endif
c--


c---- Test

c     WBC only in NH
c
c      if (lat(jp).ge.0.d0) goto 100


c     WBC in NH and SH

      if (abs(ij).eq.2 .and. jp.eq.1) goto 101
      goto 100
101   continue


c---- Total number of faces through which WBC flows =
c     Total number of linear observational constraints
c
c      n_wbc = kco




c     ------------------
c     End of the routine
c     ------------------

      return

      end

c     ******************************************************************

