      MODULE t3dbc_mod
!
!svn $Id: t3dbc_im.F 294 2009-01-09 21:37:26Z arango $
!=======================================================================
!  Copyright (c) 2002-2009 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                           Hernan G. Arango   !
!========================================== Alexander F. Shchepetkin ===
!                                                                      !
!  This subroutine sets lateral boundary conditions for the ITRC-th    !
!  tracer field.                                                       !
!                                                                      !
!=======================================================================
!
      implicit none
      PRIVATE
      PUBLIC  :: t3dbc_tile
      CONTAINS
!
!***********************************************************************
      SUBROUTINE t3dbc (ng, tile, nout, itrc)
!***********************************************************************
!
      USE mod_param
      USE mod_ocean
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, nout, itrc
!
!  Local variable declarations.
!
      integer :: IminS, ImaxS, JminS, JmaxS
      integer :: LBi, UBi, LBj, UBj, LBij, UBij
!
!  Set horizontal starting and ending indices for automatic private storage
!  arrays.
!
      IminS=BOUNDS(ng)%Istr(tile)-3
      ImaxS=BOUNDS(ng)%Iend(tile)+3
      JminS=BOUNDS(ng)%Jstr(tile)-3
      JmaxS=BOUNDS(ng)%Jend(tile)+3
!
!  Determine array lower and upper bounds in the I- and J-directions.
!
      LBi=BOUNDS(ng)%LBi(tile)
      UBi=BOUNDS(ng)%UBi(tile)
      LBj=BOUNDS(ng)%LBj(tile)
      UBj=BOUNDS(ng)%UBj(tile)
!
!  Set array lower and upper bounds for MIN(I,J)- and MAX(I,J)-directions.
!
      LBij=BOUNDS(ng)%LBij
      UBij=BOUNDS(ng)%UBij
!
      CALL t3dbc_tile (ng, tile, itrc,                                  &
     &                 LBi, UBi, LBj, UBj, N(ng), NT(ng),               &
     &                 IminS, ImaxS, JminS, JmaxS,                      &
     &                 nstp(ng), nout,                                  &
     &                 OCEAN(ng)% t)
      RETURN
      END SUBROUTINE t3dbc
!
!***********************************************************************
      SUBROUTINE t3dbc_tile (ng, tile, itrc,                            &
     &                       LBi, UBi, LBj, UBj, UBk, UBt,              &
     &                       IminS, ImaxS, JminS, JmaxS,                &
     &                       nstp, nout,                                &
     &                       t)
!***********************************************************************
!
      USE mod_param
      USE mod_boundary
      USE mod_grid
      USE mod_scalars
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, itrc
      integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: nstp, nout
!
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
!
!  Local variable declarations.
!
      integer :: i, j, k
      real(r8), parameter :: eps =1.0E-20_r8
      real(r8) :: Ce, Cx, cff, dTde, dTdt, dTdx, tau
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
!
!-----------------------------------------------------------------------
!  Set lower and upper tile bounds and staggered variables bounds for
!  this horizontal domain partition.  Notice that if tile=-1, it will
!  set the values for the global grid.
!-----------------------------------------------------------------------
!
      integer :: Istr, IstrR, IstrT, IstrU, Iend, IendR, IendT
      integer :: Jstr, JstrR, JstrT, JstrV, Jend, JendR, JendT
!
      Istr =BOUNDS(ng)%Istr (tile)
      IstrR=BOUNDS(ng)%IstrR(tile)
      IstrT=BOUNDS(ng)%IstrT(tile)
      IstrU=BOUNDS(ng)%IstrU(tile)
      Iend =BOUNDS(ng)%Iend (tile)
      IendR=BOUNDS(ng)%IendR(tile)
      IendT=BOUNDS(ng)%IendT(tile)
      Jstr =BOUNDS(ng)%Jstr (tile)
      JstrR=BOUNDS(ng)%JstrR(tile)
      JstrT=BOUNDS(ng)%JstrT(tile)
      JstrV=BOUNDS(ng)%JstrV(tile)
      Jend =BOUNDS(ng)%Jend (tile)
      JendR=BOUNDS(ng)%JendR(tile)
      JendT=BOUNDS(ng)%JendT(tile)
!
!-----------------------------------------------------------------------
!  Lateral boundary conditions at the western edge.
!-----------------------------------------------------------------------
!
      IF (Istr.eq.1) THEN
!
!  Western edge, closed boundary condition.
!
        DO k=1,N(ng)
          DO j=Jstr,Jend
            t(Istr-1,j,k,nout,itrc)=t(Istr,j,k,nout,itrc)
            t(Istr-1,j,k,nout,itrc)=t(Istr-1,j,k,nout,itrc)*            &
     &                              GRID(ng)%rmask(Istr-1,j)
          END DO
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Lateral boundary conditions at the eastern edge.
!-----------------------------------------------------------------------
!
      IF (Iend.eq.Lm(ng)) THEN
!
!  Eastern edge, implicit upstream radiation condition.
!
        DO k=1,N(ng)
          DO j=Jstr,Jend+1
           grad(Iend  ,j)=t(Iend  ,j  ,k,nstp,itrc)-                    &
     &                    t(Iend  ,j-1,k,nstp,itrc)
           grad(Iend  ,j)=grad(Iend  ,j)*                               &
     &                    GRID(ng)%vmask(Iend  ,j)
           grad(Iend+1,j)=t(Iend+1,j  ,k,nstp,itrc)-                    &
     &                    t(Iend+1,j-1,k,nstp,itrc)
           grad(Iend+1,j)=grad(Iend+1,j)*                               &
     &                    GRID(ng)%vmask(Iend+1,j)
          END DO
          DO j=Jstr,Jend
            dTdt=t(Iend,j,k,nstp,itrc)-t(Iend  ,j,k,nout,itrc)
            dTdx=t(Iend,j,k,nout,itrc)-t(Iend-1,j,k,nout,itrc)
            tau=Tobc_out(itrc,ng,ieast)
            IF ((dTdt*dTdx).lt.0.0_r8) tau=Tobc_in(itrc,ng,ieast)
            tau=tau*dt(ng)
            IF ((dTdt*dTdx).lt.0.0_r8) dTdt=0.0_r8
            IF ((dTdt*(grad(Iend,j)+grad(Iend,j+1))).gt.0.0_r8) THEN
              dTde=grad(Iend,j  )
            ELSE
              dTde=grad(Iend,j+1)
            END IF
            cff=MAX(dTdx*dTdx+dTde*dTde,eps)
            Cx=dTdt*dTdx
            Ce=0.0_r8
            t(Iend+1,j,k,nout,itrc)=(cff*t(Iend+1,j,k,nstp,itrc)+       &
     &                               Cx *t(Iend  ,j,k,nout,itrc)-       &
     &                               MAX(Ce,0.0_r8)*                    &
     &                                  grad(Iend+1,j  )-               &
     &                               MIN(Ce,0.0_r8)*                    &
     &                                  grad(Iend+1,j+1))/              &
     &                              (cff+Cx)
            t(Iend+1,j,k,nout,itrc)=t(Iend+1,j,k,nout,itrc)+            &
     &                              tau*(BOUNDARY(ng)%t_east(j,k,itrc)- &
     &                                   t(Iend+1,j,k,nstp,itrc))
            t(Iend+1,j,k,nout,itrc)=t(Iend+1,j,k,nout,itrc)*            &
     &                              GRID(ng)%rmask(Iend+1,j)
          END DO
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Lateral boundary conditions at the southern edge.
!-----------------------------------------------------------------------
!
      IF (Jstr.eq.1) THEN
!
!  Southern edge, closed boundary condition.
!
        DO k=1,N(ng)
          DO i=Istr,Iend
            t(i,Jstr-1,k,nout,itrc)=t(i,Jstr,k,nout,itrc)
            t(i,Jstr-1,k,nout,itrc)=t(i,Jstr-1,k,nout,itrc)*            &
     &                              GRID(ng)%rmask(i,Jstr-1)
          END DO
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Lateral boundary conditions at the northern edge.
!-----------------------------------------------------------------------
!
      IF (Jend.eq.Mm(ng)) THEN
!
!  Northern edge, closed boundary condition.
!
        DO k=1,N(ng)
          DO i=Istr,Iend
            t(i,Jend+1,k,nout,itrc)=t(i,Jend,k,nout,itrc)
            t(i,Jend+1,k,nout,itrc)=t(i,Jend+1,k,nout,itrc)*            &
     &                              GRID(ng)%rmask(i,Jend+1)
          END DO
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Boundary corners.
!-----------------------------------------------------------------------
!
      IF ((Jstr.eq.1).and.(Istr.eq.1)) THEN
        DO k=1,N(ng)
          t(Istr-1,Jstr-1,k,nout,itrc)=0.5_r8*                          &
     &                                 (t(Istr  ,Jstr-1,k,nout,itrc)+   &
     &                                  t(Istr-1,Jstr  ,k,nout,itrc))
        END DO
      END IF
      IF ((Jstr.eq.1).and.(Iend.eq.Lm(ng))) THEN
        DO k=1,N(ng)
          t(Iend+1,Jstr-1,k,nout,itrc)=0.5_r8*                          &
     &                                 (t(Iend  ,Jstr-1,k,nout,itrc)+   &
     &                                  t(Iend+1,Jstr  ,k,nout,itrc))
        END DO
      END IF
      IF ((Jend.eq.Mm(ng)).and.(Istr.eq.1)) THEN
        DO k=1,N(ng)
          t(Istr-1,Jend+1,k,nout,itrc)=0.5_r8*                          &
     &                                 (t(Istr-1,Jend  ,k,nout,itrc)+   &
     &                                  t(Istr  ,Jend+1,k,nout,itrc))
        END DO
      END IF
      IF ((Jend.eq.Mm(ng)).and.(Iend.eq.Lm(ng))) THEN
        DO k=1,N(ng)
          t(Iend+1,Jend+1,k,nout,itrc)=0.5_r8*                          &
     &                                 (t(Iend+1,Jend  ,k,nout,itrc)+   &
     &                                  t(Iend  ,Jend+1,k,nout,itrc))
        END DO
      END IF
      RETURN
      END SUBROUTINE t3dbc_tile
      END MODULE t3dbc_mod
