!$Id: bdy_3d.F90,v 1.10 2006-03-01 15:54:08 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE:  bdy_3d - 3D boundary conditions \label{bdy-3d}
!
! !INTERFACE:
   module bdy_3d
!
! !DESCRIPTION:
!  
! Here, the three-dimensional boundary
! conditions for temperature and salinity are handled.
!
! !USES:
   use halo_zones, only : H_TAG,U_TAG,V_TAG
   use domain, only: imin,jmin,imax,jmax,H,az,au,av
   use domain, only: iimin,jjmin,iimax,jjmax,kmax
   use domain, only: nsbv,NWB,NNB,NEB,NSB,bdy_index
   use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
   use variables_3d
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
   public init_bdy_3d, do_bdy_3d
   REALTYPE, public, allocatable       :: S_bdy(:,:),T_bdy(:,:)
!
! !PRIVATE DATA MEMBERS:
!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_bdy_3d - initialising 3D boundary conditions
! \label{sec-init-bdy-3d}
!
! !INTERFACE:
   subroutine init_bdy_3d()
!
! !DESCRIPTION:
!
! Here, the necessary fields {\tt S\_bdy} and {\tt T\_bdy} for
! salinity and temperature, respectively, are allocated.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
   integer                   :: rc,i,j,k,n
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_bdy_3d() # ',Ncall
#endif

   LEVEL2 'init_bdy_3d()'
   allocate(S_bdy(0:kmax,nsbv),stat=rc)
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (S_bdy)'

   allocate(T_bdy(0:kmax,nsbv),stat=rc)
   if (rc /= 0) stop 'init_init_bdy_3d: Error allocating memory (T_bdy)'

#ifdef DEBUG
   write(debug,*) 'Leaving init_bdy_3d()'
   write(debug,*)
#endif
   return
   end subroutine init_bdy_3d
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE:  do_bdy_3d  - updating 3D boundary conditions
! \label{sec-do-bdy-3d}
!
! !INTERFACE:
   subroutine do_bdy_3d(tag,field)
!
! !DESCRIPTION:
!
! Here, the boundary conditions for salinity and temperature are
! copied to the boundary points and relaxed to the near boundary points
! by means of the flow relaxation scheme by \cite{MARTINSENea87}.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in)                 :: tag
!
! !INPUT/OUTPUT PARAMETERS:
   REALTYPE, intent(inout)             :: field(I3DFIELD)
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
   integer                   :: i,j,k,l,n,ii,jj
   REALTYPE                  :: sp(1:4),rat
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_bdy_3d() # ',Ncall
#endif

#if 0
   select case (tag)
      case (1)
      ! Lateral zero-gradient boundary condition (north & south)
         do k=1,kmax
            do i=iimin,iimax
               if (au(i,jjmin) .eq. 3) field(i,jjmin,k)=field(i,jjmin+1,k)
               if (au(i,jjmax) .eq. 3) field(i,jjmax,k)=field(i,jjmax-1,k)
            end do
         end do
      case (2)
      ! Lateral zero-gradient boundary conditions (west & east)
         do k=1,kmax
            do j=jjmin,jjmax
               if (av(iimin,j) .eq. 3) field(iimin,j,k)=field(iimin+1,j,k)
               if (av(iimax,j) .eq. 3) field(iimax,j,k)=field(iimax-1,j,k)
            end do
         end do
      case default
         FATAL 'Non valid tag'
         stop 'do_bdy_3d'
   end select
#endif

#ifndef NO_BAROCLINIC
! Sponge layer factors according to Martinsen and Engedahl, 1987.
   sp(1)=1.0
   sp(2)=0.5625
   sp(3)=0.25
   sp(4)=0.0625

   l = 0
   k = 0
   do n=1,NWB
      l = l+1
      k = bdy_index(l)
      i = wi(n)
      do j=wfj(n),wlj(n)
         do ii=1,4
            if (az(i-1+ii,j).gt.0) then
               S(i-1+ii,j,:) = sp(ii)*S_bdy(:,k)+(1.-sp(ii))*S(i-1+ii,j,:)
               T(i-1+ii,j,:) = sp(ii)*T_bdy(:,k)+(1.-sp(ii))*T(i-1+ii,j,:)
            end if
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i+1,j,:)
            end if
#endif
         end do
         k = k+1
      end do
   end do

   do n = 1,NNB
      l = l+1
      k = bdy_index(l)
      j = nj(n)
      do i = nfi(n),nli(n)
         do jj=1,4
            if (az(i,j+1-jj).gt.0) then
               S(i,j+1-jj,:) = sp(jj)*S_bdy(:,k)+(1.-sp(jj))*S(i,j+1-jj,:)
               T(i,j+1-jj,:) = sp(jj)*T_bdy(:,k)+(1.-sp(jj))*T(i,j+1-jj,:)
            end if
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i,j-1,:)
            end if
#endif
         end do
         k = k+1
      end do
   end do

   do n=1,NEB
      l = l+1
      k = bdy_index(l)
      i = ei(n)
      do j=efj(n),elj(n)
         do ii=1,4
            if (az(i+1-ii,j).gt.0) then
               S(i+1-ii,j,:) = sp(ii)*S_bdy(:,k)+(1.-sp(ii))*S(i+1-ii,j,:)
               T(i+1-ii,j,:) = sp(ii)*T_bdy(:,k)+(1.-sp(ii))*T(i+1-ii,j,:)
            end if
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i-1,j,:)
            end if
#endif
         end do
         k = k+1
      end do
   end do

   do n = 1,NSB
      l = l+1
      k = bdy_index(l)
      j = sj(n)
      do i = sfi(n),sli(n)
         do jj=1,4
            if (az(i,j-1+jj).gt.0) then
               S(i,j-1+jj,:) = sp(jj)*S_bdy(:,k)+(1.-sp(jj))*S(i,j-1+jj,:)
               T(i,j-1+jj,:) = sp(jj)*T_bdy(:,k)+(1.-sp(jj))*T(i,j-1+jj,:)
            end if
#ifdef GETM_BIO
            if ( allocated(cc3d) ) then
               cc3d(:,i,j,:) = cc3d(:,i,j+1,:)
            end if
#endif
         end do
         k = k+1
      end do
   end do

   call mirror_bdy_3d(T,H_TAG)
   call mirror_bdy_3d(S,H_TAG)

#ifdef GETM_BIO
   if ( allocated(cc3d) ) then
      do n=1,size(cc3d,1)
         call mirror_bdy_3d(cc3d(n,:,:,:),H_TAG)
      end do
   end if
#endif


#endif

#ifdef DEBUG
   write(debug,*) 'leaving do_bdy_3d()'
   write(debug,*)
#endif
   return
   end subroutine do_bdy_3d
!EOC

!-----------------------------------------------------------------------

   end module bdy_3d

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Karsten Bolding and Hans Burchard               !
!-----------------------------------------------------------------------