#include "cppdefs.h"
MODULE AC3dbc_mod
#ifdef SOLVE3D
!
!svn $Id: t3dbc_im.F 732 2008-09-07 01:55:51Z jcwarner $
!=======================================================================
! Copyright (c) 2002-2019 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.txt Hernan G. Arango !
! John C. Warner !
! !
! This subroutine sets lateral boundary conditions for the !
! wave action density field. !
! !
!=======================================================================
!
implicit none
PRIVATE
PUBLIC :: AC3dbc_tile
CONTAINS
!
!***********************************************************************
SUBROUTINE AC3dbc (ng, tile, nout)
!***********************************************************************
!
USE mod_param
USE mod_inwave_vars
USE mod_ocean
USE mod_stepping
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, nout
!
! Local variable declarations.
!
# include "tile.h"
!
CALL AC3dbc_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nstp(ng),nout, &
& WAVEP(ng)% AC)
RETURN
END SUBROUTINE AC3dbc
!
!***********************************************************************
SUBROUTINE AC3dbc_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nstp,nout, &
& AC)
!***********************************************************************
USE mod_param
USE mod_inwave_params
USE mod_inwave_bound
USE mod_inwave_vars
USE mod_boundary
USE mod_grid
USE mod_scalars
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: nstp, nout
!
# ifdef ASSUMED_SHAPE
real(r8), intent(inout) :: AC(LBi:,LBj:,:,:)
# else
real(r8), intent(inout) :: AC(LBi:UBi,LBj:UBj,ND,3)
# endif
!
! Local variable declarations.
!
integer :: i, j, k, d_bnd
real(r8), parameter :: eps =1.0E-20_r8
real(r8) :: Ce, Cx, cff, dTde, dTdt, dTdx, tau, ramp
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
# include "set_bounds.h"
# ifdef RAMP_INWAVE
ramp=TANH((tdays(ng)-dstart)/0.14_r8)
# else
ramp=1.0_r8
# endif
!
!-----------------------------------------------------------------------
! Lateral boundary conditions at the western edge.
!-----------------------------------------------------------------------
!
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
!
! Western edge, implicit upstream radiation condition.
!
IF (LBC(iwest,isAC3d,ng)%radiation) THEN
DO k=1,ND
DO j=Jstr,Jend+1
grad(Istr-1,j)=AC(Istr-1,j ,k,nstp)- &
& AC(Istr-1,j-1,k,nstp)
# ifdef MASKING
grad(Istr-1,j)=grad(Istr-1,j)* &
& GRID(ng)%vmask(Istr-1,j)
# endif
grad(Istr ,j)=AC(Istr ,j ,k,nstp)- &
& AC(Istr ,j-1,k,nstp)
# ifdef MASKING
grad(Istr ,j)=grad(Istr ,j)* &
& GRID(ng)%vmask(Istr ,j)
# endif
END DO
DO j=Jstr,Jend
IF (LBC_apply(ng)%west(j)) THEN
dTdt=AC(Istr,j,k,nstp)-AC(Istr ,j,k,nout)
dTdx=AC(Istr,j,k,nout)-AC(Istr+1,j,k,nout)
IF ((dTdt*dTdx).lt.0.0_r8) dTdt=0.0_r8
IF ((dTdt*(grad(Istr,j)+grad(Istr,j+1))).gt.0.0_r8) THEN
dTde=grad(Istr,j )
ELSE
dTde=grad(Istr,j+1)
END IF
cff=MAX(dTdx*dTdx+dTde*dTde,eps)
Cx=dTdt*dTdx
Ce=MIN(cff,MAX(dTdt*dTde,-cff))
AC(Istr-1,j,k,nout)=(cff*AC(Istr-1,j,k,nstp)+ &
& Cx *AC(Istr ,j,k,nout)- &
& MAX(Ce,0.0_r8)* &
& grad(Istr-1,j )- &
& MIN(Ce,0.0_r8)* &
& grad(Istr-1,j+1))/ &
& (cff+Cx)
# ifdef MASKING
AC(Istr-1,j,k,nout)=AC(Istr-1,j,k,nout)* &
& GRID(ng)%rmask(Istr-1,j)
# endif
END IF
END DO
END DO
!
! Western edge, clamped boundary condition.
!
ELSE IF (LBC(iwest,isAC3d,ng)%clamped) THEN
DO k=1,ND
DO d_bnd=1,WAVEB(ng)%ND_bnd
IF(WAVEB(ng)%WD_bnd(d_bnd).EQ.WAVEG(ng)%wd(k))then
DO j=Jstr,Jend
IF (LBC_apply(ng)%west(j)) THEN
AC(Istr-1,j,k,nout)=WAVEB(ng)%AC_west(j,d_bnd)*ramp
# ifdef MASKING
AC(Istr-1,j,k,nout)=AC(Istr-1,j,k,nout)* &
& GRID(ng)%rmask(Istr-1,j)
# endif
END IF
END DO
END IF
END DO
END DO
!
! Western edge, gradient boundary condition.
!
ELSE IF (LBC(iwest,isAC3d,ng)%gradient) THEN
DO k=1,ND
DO j=Jstr,Jend
IF (LBC_apply(ng)%west(j)) THEN
AC(Istr-1,j,k,nout)=AC(Istr,j,k,nout)
# ifdef MASKING
AC(Istr-1,j,k,nout)=AC(Istr-1,j,k,nout)* &
& GRID(ng)%rmask(Istr-1,j)
# endif
END IF
END DO
END DO
!
! Western edge, closed boundary condition.
!
ELSE IF (LBC(iwest,isAC3d,ng)%closed) THEN
DO k=1,ND
DO j=Jstr,Jend
IF (LBC_apply(ng)%west(j)) THEN
AC(Istr-1,j,k,nout)=AC(Istr,j,k,nout)
# ifdef MASKING
AC(Istr-1,j,k,nout)=AC(Istr-1,j,k,nout)* &
& GRID(ng)%rmask(Istr-1,j)
# endif
END IF
END DO
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Lateral boundary conditions at the eastern edge.
!-----------------------------------------------------------------------
!
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
!
! Eastern edge, implicit upstream radiation condition.
!
IF (LBC(ieast,isAC3d,ng)%radiation) THEN
DO k=1,ND
DO j=Jstr,Jend+1
grad(Iend ,j)=AC(Iend ,j ,k,nstp)- &
& AC(Iend ,j-1,k,nstp)
# ifdef MASKING
grad(Iend ,j)=grad(Iend ,j)* &
& GRID(ng)%vmask(Iend ,j)
# endif
grad(Iend+1,j)=AC(Iend+1,j ,k,nstp)- &
& AC(Iend+1,j-1,k,nstp)
# ifdef MASKING
grad(Iend+1,j)=grad(Iend+1,j)* &
& GRID(ng)%vmask(Iend+1,j)
# endif
END DO
DO j=Jstr,Jend
IF (LBC_apply(ng)%east(j)) THEN
dTdt=AC(Iend,j,k,nstp)-AC(Iend ,j,k,nout)
dTdx=AC(Iend,j,k,nout)-AC(Iend-1,j,k,nout)
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=MIN(cff,MAX(dTdt*dTde,-cff))
AC(Iend+1,j,k,nout)=(cff*AC(Iend+1,j,k,nstp)+ &
& Cx *AC(Iend ,j,k,nout)- &
& MAX(Ce,0.0_r8)* &
& grad(Iend+1,j )- &
& MIN(Ce,0.0_r8)* &
& grad(Iend+1,j+1))/ &
& (cff+Cx)
# ifdef MASKING
AC(Iend+1,j,k,nout)=AC(Iend+1,j,k,nout)* &
& GRID(ng)%rmask(Iend+1,j)
# endif
END IF
END DO
END DO
!
! Eastern edge, clamped boundary condition.
!
ELSE IF (LBC(ieast,isAC3d,ng)%clamped) THEN
DO d_bnd=1,WAVEB(ng)%ND_bnd
DO k=1,ND
IF(WAVEB(ng)%WD_bnd(d_bnd).EQ.WAVEG(ng)%wd(k)) THEN
DO j=Jstr,Jend
IF (LBC_apply(ng)%east(j)) THEN
AC(Iend+1,j,k,nout)=WAVEB(ng)%AC_east(j,d_bnd)*ramp
# ifdef MASKING
AC(Iend+1,j,k,nout)=AC(Iend+1,j,k,nout)* &
& GRID(ng)%rmask(Iend+1,j)
# endif
END IF
END DO
END IF
END DO
END DO
!
! Eastern edge, gradient boundary condition.
!
ELSE IF (LBC(ieast,isAC3d,ng)%gradient) THEN
DO k=1,ND
DO j=Jstr,Jend
IF (LBC_apply(ng)%east(j)) THEN
AC(Iend+1,j,k,nout)=AC(Iend,j,k,nout)
# ifdef MASKING
AC(Iend+1,j,k,nout)=AC(Iend+1,j,k,nout)* &
& GRID(ng)%rmask(Iend+1,j)
# endif
END IF
END DO
END DO
!
! Eastern edge, closed boundary condition.
!
ELSE IF (LBC(ieast,isAC3d,ng)%closed) THEN
DO k=1,ND
DO j=Jstr,Jend
IF (LBC_apply(ng)%east(j)) THEN
AC(Iend+1,j,k,nout)=AC(Iend,j,k,nout)
# ifdef MASKING
AC(Iend+1,j,k,nout)=AC(Iend+1,j,k,nout)* &
& GRID(ng)%rmask(Iend+1,j)
# endif
END IF
END DO
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Lateral boundary conditions at the southern edge.
!-----------------------------------------------------------------------
!
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
!
! Southern edge, implicit upstream radiation condition.
!
IF (LBC(isouth,isAC3d,ng)%radiation) THEN
DO k=1,ND
DO i=Istr,Iend+1
grad(i,Jstr )=AC(i ,Jstr ,k,nstp)- &
& AC(i-1,Jstr ,k,nstp)
# ifdef MASKING
grad(i,Jstr )=grad(i,Jstr )* &
& GRID(ng)%umask(i,Jstr )
# endif
grad(i,Jstr-1)=AC(i ,Jstr-1,k,nstp)- &
& AC(i-1,Jstr-1,k,nstp)
# ifdef MASKING
grad(i,Jstr-1)=grad(i,Jstr-1)* &
& GRID(ng)%umask(i,Jstr-1)
# endif
END DO
DO i=Istr,Iend
IF (LBC_apply(ng)%south(i)) THEN
dTdt=AC(i,Jstr,k,nstp)-AC(i,Jstr ,k,nout)
dTde=AC(i,Jstr,k,nout)-AC(i,Jstr+1,k,nout)
IF ((dTdt*dTde).lt.0.0_r8) dTdt=0.0_r8
IF ((dTdt*(grad(i,Jstr)+grad(i+1,Jstr))).gt.0.0_r8) THEN
dTdx=grad(i ,Jstr)
ELSE
dTdx=grad(i+1,Jstr)
END IF
cff=MAX(dTdx*dTdx+dTde*dTde,eps)
Cx=MIN(cff,MAX(dTdt*dTdx,-cff))
Ce=dTdt*dTde
AC(i,Jstr-1,k,nout)=(cff*AC(i,Jstr-1,k,nstp)+ &
& Ce *AC(i,Jstr ,k,nout)- &
& MAX(Cx,0.0_r8)* &
& grad(i ,Jstr-1)- &
& MIN(Cx,0.0_r8)* &
& grad(i+1,Jstr-1))/ &
& (cff+Ce)
# ifdef MASKING
AC(i,Jstr-1,k,nout)=AC(i,Jstr-1,k,nout)* &
& GRID(ng)%rmask(i,Jstr-1)
# endif
END IF
END DO
END DO
!
! Southern edge, clamped boundary condition.
!
ELSE IF (LBC(isouth,isAC3d,ng)%clamped) THEN
DO d_bnd=1,WAVEB(ng)%ND_bnd
DO k=1,ND
IF(WAVEB(ng)%WD_bnd(d_bnd).EQ.WAVEG(ng)%wd(k))then
DO i=Istr,Iend
IF (LBC_apply(ng)%south(i)) THEN
AC(i,Jstr-1,k,nout)=WAVEB(ng)%AC_south(i,d_bnd)*ramp
# ifdef MASKING
AC(i,Jstr-1,k,nout)=AC(i,Jstr-1,k,nout)* &
& GRID(ng)%rmask(i,Jstr-1)
# endif
END IF
END DO
END IF
END DO
END DO
!
! Southern edge, gradient boundary condition.
!
ELSE IF (LBC(isouth,isAC3d,ng)%gradient) THEN
DO k=1,ND
DO i=Istr,Iend
IF (LBC_apply(ng)%south(i)) THEN
AC(i,Jstr-1,k,nout)=AC(i,Jstr,k,nout)
# ifdef MASKING
AC(i,Jstr-1,k,nout)=AC(i,Jstr-1,k,nout)* &
& GRID(ng)%rmask(i,Jstr-1)
# endif
END IF
END DO
END DO
!
! Southern edge, closed boundary condition.
!
ELSE IF (LBC(isouth,isAC3d,ng)%closed) THEN
DO k=1,ND
DO i=Istr,Iend
IF (LBC_apply(ng)%south(i)) THEN
AC(i,Jstr-1,k,nout)=AC(i,Jstr,k,nout)
# ifdef MASKING
AC(i,Jstr-1,k,nout)=AC(i,Jstr-1,k,nout)* &
& GRID(ng)%rmask(i,Jstr-1)
# endif
END IF
END DO
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Lateral boundary conditions at the northern edge.
!-----------------------------------------------------------------------
!
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
!
! Northern edge, implicit upstream radiation condition.
!
IF (LBC(inorth,isAC3d,ng)%radiation) THEN
DO k=1,ND
DO i=Istr,Iend+1
grad(i,Jend )=AC(i ,Jend ,k,nstp)- &
& AC(i-1,Jend ,k,nstp)
# ifdef MASKING
grad(i,Jend )=grad(i,Jend )* &
& GRID(ng)%umask(i,Jend )
# endif
grad(i,Jend+1)=AC(i ,Jend+1,k,nstp)- &
& AC(i-1,Jend+1,k,nstp)
# ifdef MASKING
grad(i,Jend+1)=grad(i,Jend+1)* &
& GRID(ng)%umask(i,Jend+1)
# endif
END DO
DO i=Istr,Iend
IF (LBC_apply(ng)%north(i)) THEN
dTdt=AC(i,Jend,k,nstp)-AC(i,Jend ,k,nout)
dTde=AC(i,Jend,k,nout)-AC(i,Jend-1,k,nout)
IF ((dTdt*dTde).lt.0.0_r8) dTdt=0.0_r8
IF ((dTdt*(grad(i,Jend)+grad(i+1,Jend))).gt.0.0_r8) THEN
dTdx=grad(i ,Jend)
ELSE
dTdx=grad(i+1,Jend)
END IF
cff=MAX(dTdx*dTdx+dTde*dTde,eps)
Cx=MIN(cff,MAX(dTdt*dTdx,-cff))
Ce=dTdt*dTde
AC(i,Jend+1,k,nout)=(cff*AC(i,Jend+1,k,nstp)+ &
& Ce *AC(i,Jend ,k,nout)- &
& MAX(Cx,0.0_r8)* &
& grad(i ,Jend+1)- &
& MIN(Cx,0.0_r8)* &
& grad(i+1,Jend+1))/ &
& (cff+Ce)
# ifdef MASKING
AC(i,Jend+1,k,nout)=AC(i,Jend+1,k,nout)* &
& GRID(ng)%rmask(i,Jend+1)
# endif
END IF
END DO
END DO
!
! Northern edge, clamped boundary condition.
!
ELSE IF (LBC(inorth,isAC3d,ng)%clamped) THEN
DO d_bnd=1,WAVEB(ng)%ND_bnd
DO k=1,ND
IF(WAVEB(ng)%WD_bnd(d_bnd).EQ.WAVEG(ng)%wd(k))then
DO i=Istr,Iend
IF (LBC_apply(ng)%north(i)) THEN
AC(i,Jend+1,k,nout)=WAVEB(ng)%AC_north(i,d_bnd)*ramp
# ifdef MASKING
AC(i,Jend+1,k,nout)=AC(i,Jend+1,k,nout)* &
& GRID(ng)%rmask(i,Jend+1)
# endif
END IF
END DO
END IF
END DO
END DO
!
! Northern edge, gradient boundary condition.
!
ELSE IF (LBC(inorth,isAC3d,ng)%gradient) THEN
DO k=1,ND
DO i=Istr,Iend
IF (LBC_apply(ng)%north(i)) THEN
AC(i,Jend+1,k,nout)=AC(i,Jend,k,nout)
# ifdef MASKING
AC(i,Jend+1,k,nout)=AC(i,Jend+1,k,nout)* &
& GRID(ng)%rmask(i,Jend+1)
# endif
END IF
END DO
END DO
!
! Northern edge, closed boundary condition.
!
ELSE IF (LBC(inorth,isAC3d,ng)%closed) THEN
DO k=1,ND
DO i=Istr,Iend
IF (LBC_apply(ng)%north(i)) THEN
AC(i,Jend+1,k,nout)=AC(i,Jend,k,nout)
# ifdef MASKING
AC(i,Jend+1,k,nout)=AC(i,Jend+1,k,nout)* &
& GRID(ng)%rmask(i,Jend+1)
# endif
END IF
END DO
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Boundary corners.
!-----------------------------------------------------------------------
!
IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
IF (LBC_apply(ng)%south(Istr-1).and. &
& LBC_apply(ng)%west (Jstr-1)) THEN
DO k=1,ND
AC(Istr-1,Jstr-1,k,nout)=0.5_r8* &
& (AC(Istr ,Jstr-1,k,nout)+ &
& AC(Istr-1,Jstr ,k,nout))
END DO
END IF
END IF
IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
IF (LBC_apply(ng)%south(Iend+1).and. &
& LBC_apply(ng)%east (Jstr-1)) THEN
DO k=1,ND
AC(Iend+1,Jstr-1,k,nout)=0.5_r8* &
& (AC(Iend ,Jstr-1,k,nout)+ &
& AC(Iend+1,Jstr ,k,nout))
END DO
END IF
END IF
IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
IF (LBC_apply(ng)%north(Istr-1).and. &
& LBC_apply(ng)%west (Jend+1)) THEN
DO k=1,ND
AC(Istr-1,Jend+1,k,nout)=0.5_r8* &
& (AC(Istr-1,Jend ,k,nout)+ &
& AC(Istr ,Jend+1,k,nout))
END DO
END IF
END IF
IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
IF (LBC_apply(ng)%north(Iend+1).and. &
& LBC_apply(ng)%east (Jend+1)) THEN
DO k=1,ND
AC(Iend+1,Jend+1,k,nout)=0.5_r8* &
& (AC(Iend+1,Jend ,k,nout)+ &
& AC(Iend ,Jend+1,k,nout))
END DO
END IF
END IF
END IF
RETURN
END SUBROUTINE AC3dbc_tile
#endif
END MODULE AC3dbc_mod