C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_calc.F,v 1.29 2010/12/06 14:25:36 jmc Exp $
C $Name: $
#include "OBCS_OPTIONS.h"
SUBROUTINE OBCS_CALC( futureTime, futureIter,
& uVel, vVel, wVel, theta, salt,
& myThid )
C *==========================================================*
C | SUBROUTINE OBCS_CALC |
C | o Calculate future boundary data at open boundaries |
C | at time = futureTime |
C *==========================================================*
C | |
C *==========================================================*
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "OBCS.h"
#ifdef ALLOW_PTRACERS
# include "PTRACERS_SIZE.h"
# include "PTRACERS_PARAMS.h"
# include "PTRACERS_FIELDS.h"
# include "OBCS_PTRACERS.h"
#endif /* ALLOW_PTRACERS */
#ifdef ALLOW_NEST_CHILD
# include "NEST_CHILD.h"
#endif /* ALLOW_NEST_CHILD */
C == Routine arguments ==
INTEGER futureIter
_RL futureTime
_RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
INTEGER myThid
#ifdef ALLOW_OBCS
C == Local variables ==
C bi, bj :: tile indices
C I,J,K :: loop indices
C I_obc, J_obc :: local index of open boundary
C msgBuf :: Informational/error message buffer
INTEGER bi, bj
INTEGER I, J, K, I_obc, J_obc
#ifdef ALLOW_OBCS_BALANCE
_RL Tr_T, Ar_T, Ar
#endif /* ALLOW_OBCS_BALANCE */
#ifdef ALLOW_PTRACERS
CHARACTER*(MAX_LEN_MBUF) msgBuf
INTEGER iTracer
#endif /* ALLOW_PTRACERS */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid)
#endif
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_NEST_CHILD
IF ( useNEST_CHILD ) THEN
IF ( PASSI.LT.2 ) THEN
CALL NEST_CHILD_RECV ( myThid )
ENDIF
ENDIF
#endif /* ALLOW_NEST_CHILD */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_OBCS_EAST
C Eastern OB
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: East',myThid)
#endif
IF (useOrlanskiEast) THEN
#ifdef ALLOW_ORLANSKI
CALL ORLANSKI_EAST(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
#endif
#ifdef ALLOW_NEST_CHILD
ELSEIF ( useNEST_CHILD ) THEN
DO k=1,Nr
DO j=1-Oly,sNy+Oly
I_obc=OB_Ie(j,bi,bj)
IF (I_obc.ne.0) THEN
OBEu(j,k,bi,bj)= U_F1(j,k,2)
OBEv(j,k,bi,bj)= V_F1(j,k,2)
OBEt(j,k,bi,bj)= T_F1(j,k,2)
OBEs(j,k,bi,bj)= S_F1(j,k,2)
#ifdef NONLIN_FRSURF
OBEeta(j,bi,bj)= ETA_F1(j,1,2)
#endif
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_NEST_CHILD */
ELSE
DO K=1,Nr
DO J=1-Oly,sNy+Oly
I_obc=OB_Ie(J,bi,bj)
IF (I_obc.ne.0) THEN
OBEu(J,K,bi,bj)=0.
OBEv(J,K,bi,bj)=0.
OBEt(J,K,bi,bj)=tRef(K)
OBEs(J,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBEw(J,K,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
OBEeta(J,bi,bj)=0.
#endif
ENDIF
ENDDO
ENDDO
ENDIF
#endif /* ALLOW_OBCS_EAST */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_OBCS_WEST
C Western OB
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: West',myThid)
#endif
IF (useOrlanskiWest) THEN
#ifdef ALLOW_ORLANSKI
CALL ORLANSKI_WEST(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
#endif
#ifdef ALLOW_NEST_CHILD
ELSEIF ( useNEST_CHILD ) THEN
DO k=1,Nr
DO j=1-Oly,sNy+Oly
I_obc=OB_Iw(j,bi,bj)
IF (I_obc.ne.0) THEN
OBWu(j,k,bi,bj)= U_F1(j,k,1)
OBWv(j,k,bi,bj)= V_F1(j,k,1)
OBWt(j,k,bi,bj)= T_F1(j,k,1)
OBWs(j,k,bi,bj)= S_F1(j,k,1)
#ifdef NONLIN_FRSURF
OBWeta(j,bi,bj)= ETA_F1(j,1,1)
#endif
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_NEST_CHILD */
ELSE
DO K=1,Nr
DO J=1-Oly,sNy+Oly
I_obc=OB_Iw(J,bi,bj)
IF (I_obc.ne.0) THEN
OBWu(J,K,bi,bj)=0.
OBWv(J,K,bi,bj)=0.
OBWt(J,K,bi,bj)=tRef(K)
OBWs(J,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBWw(J,K,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
OBWeta(J,bi,bj)=0.
#endif
ENDIF
ENDDO
ENDDO
ENDIF
#endif /* ALLOW_OBCS_WEST */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_OBCS_NORTH
C Northern OB
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: North',myThid)
#endif
IF (useOrlanskiNorth) THEN
#ifdef ALLOW_ORLANSKI
CALL ORLANSKI_NORTH(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
#endif
ELSE
DO K=1,Nr
DO I=1-Olx,sNx+Olx
J_obc=OB_Jn(I,bi,bj)
IF (J_obc.ne.0) THEN
OBNv(I,K,bi,bj)=0.
OBNu(I,K,bi,bj)=0.
OBNt(I,K,bi,bj)=tRef(K)
OBNs(I,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBNw(I,K,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
OBNeta(I,bi,bj)=0.
#endif
ENDIF
ENDDO
ENDDO
ENDIF
#endif /* ALLOW_OBCS_NORTH */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_OBCS_SOUTH
C Southern OB
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: South',myThid)
#endif
IF (useOrlanskiSouth) THEN
#ifdef ALLOW_ORLANSKI
CALL ORLANSKI_SOUTH(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
#endif
ELSE
DO K=1,Nr
DO I=1-Olx,sNx+Olx
J_obc=OB_Js(I,bi,bj)
IF (J_obc.ne.0) THEN
OBSu(I,K,bi,bj)=0.
OBSv(I,K,bi,bj)=0.
OBSt(I,K,bi,bj)=tRef(K)
OBSs(I,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBSw(I,K,bi,bj)=0.
#endif
#ifdef NONLIN_FRSURF
OBSeta(I,bi,bj)=0.
#endif
ENDIF
ENDDO
ENDDO
ENDIF
#endif /* ALLOW_OBCS_SOUTH */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_PTRACERS
IF ( usePTRACERS ) THEN
C
C Calculate some default open boundary conditions for passive tracers:
C The default is a homogeneous v.Neumann conditions, that is, the
C tracer gradient across the open boundary is nearly zero;
C only nearly, because the boundary conditions are applied throughout
C the time step during which the interior field does change; therefore
C we have to use the values from the previous time step here. If you
C really want exact v.Neumann conditions, you have to modify
C obcs_apply_ptracer directly.
C
# ifdef ALLOW_OBCS_EAST
C Eastern OB
# ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC: East, pTracers',myThid)
# endif
IF (useOrlanskiEast) THEN
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: useOrlanskiEast Rad OBC with'
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: pTracers not yet implemented'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R OBCS_CALC'
ELSE
DO iTracer=1,PTRACERS_numInUse
DO K=1,Nr
DO J=1-Oly,sNy+Oly
I_obc=OB_Ie(J,bi,bj)
IF (I_obc.ne.0) THEN
OBEptr(J,K,bi,bj,iTracer) =
& pTracer(I_obc-1,J,K,bi,bj,iTracer)
& *_maskW(I_obc,J,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
# endif /* ALLOW_OBCS_EAST */
C ------------------------------------------------------------------------------
# ifdef ALLOW_OBCS_WEST
C Western OB
# ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC: West, pTracers',myThid)
# endif
IF (useOrlanskiWest) THEN
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: useOrlanskiWest Rad OBC with'
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: pTracers not yet implemented'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R OBCS_CALC'
ELSE
DO iTracer=1,PTRACERS_numInUse
DO K=1,Nr
DO J=1-Oly,sNy+Oly
I_obc=OB_Iw(J,bi,bj)
IF (I_obc.ne.0) THEN
OBWptr(J,K,bi,bj,iTracer) =
& pTracer(I_obc+1,J,K,bi,bj,iTracer)
& *_maskW(I_obc+1,J,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
# endif /* ALLOW_OBCS_WEST */
C ------------------------------------------------------------------------------
# ifdef ALLOW_OBCS_NORTH
C Northern OB
# ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC: North, pTracers',myThid)
# endif
IF (useOrlanskiNorth) THEN
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: useOrlanskiNorth Rad OBC with'
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: pTracers not yet implemented'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R OBCS_CALC'
ELSE
DO iTracer=1,PTRACERS_numInUse
DO K=1,Nr
DO I=1-Olx,sNx+Olx
J_obc=OB_Jn(I,bi,bj)
IF (J_obc.ne.0) THEN
OBNptr(I,K,bi,bj,iTracer) =
& pTracer(I,J_obc-1,K,bi,bj,iTracer)
& *_maskS(I,J_obc,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
# endif /* ALLOW_OBCS_NORTH */
C ------------------------------------------------------------------------------
# ifdef ALLOW_OBCS_SOUTH
C Southern OB
# ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC: South, pTracers',myThid)
#endif
IF (useOrlanskiSouth) THEN
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: useOrlanskiSouth Rad OBC with'
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(A)')
& 'OBCS_CALC: ERROR: pTracers not yet implemented'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R OBCS_CALC'
ELSE
DO iTracer=1,PTRACERS_numInUse
DO K=1,Nr
DO I=1-Olx,sNx+Olx
J_obc=OB_Js(I,bi,bj)
IF (J_obc.ne.0) THEN
OBSptr(I,K,bi,bj,iTracer) =
& pTracer(I,J_obc+1,K,bi,bj,iTracer)
& *_maskS(I,J_obc+1,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
# endif /* ALLOW_OBCS_SOUTH */
C end if (usePTracers)
ENDIF
#endif /* ALLOW_PTRACERS */
C-- end bi,bj loops.
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_OBCS_PRESCRIBE
IF (useOBCSprescribe) THEN
C-- Calculate future values on open boundaries
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('OBCS_PRESCRIBE_READ',myThid)
#endif
CALL OBCS_PRESCRIBE_READ( futureTime, futureIter, myThid )
ENDIF
#endif /* ALLOW_OBCS_PRESCRIBE */
C ------------------------------------------------------------------------------
#ifndef ALLOW_AUTODIFF_TAMC
C I do not think that we want this to be differentiated for now
#ifdef ALLOW_OBCS_STEVENS
C The Stevens (1990) boundary conditions come after reading data from
C files, because they use the data to compute a mix of simplified
C Orlanski and prescribed boundary conditions
IF (useStevensNorth.OR.useStevensSouth.OR.
& useStevensEast.OR.useStevensWest) THEN
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('OBCS_CALC_STEVENS',myThid)
#endif
CALL OBCS_CALC_STEVENS( futureTime, futureIter, myThid )
ENDIF
#endif /* ALLOW_OBCS_STEVENS */
#endif /* ndef ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_OBCS_BALANCE
IF ( useOBCSbalance) THEN
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_MSG('useOBCSbalance=.TRUE.',myThid)
#endif
#ifdef ALLOW_OBCS_EAST
Tr_T = 0. _d 0
Ar_T = 0. _d 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO J=1,sNy
I_obc=OB_Ie(J,bi,bj)
IF (I_obc.ne.0) THEN
Ar = drF(k)*hFacC(I_obc,j,k,bi,bj)*dyG(I_obc,j,bi,bj)
Ar_T = Ar_T + Ar
Tr_T = Tr_T + Ar * OBEu(J,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
_GLOBAL_SUM_RL( Ar_T , myThid )
IF ( Ar_T .GT. 0. _d 0 ) THEN
_GLOBAL_SUM_RL( Tr_T , myThid )
Tr_T = (0. _d 0 - Tr_T)/Ar_T
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO J=1-Oly,sNy+Oly
I_obc=OB_Ie(J,bi,bj)
IF (I_obc.ne.0) THEN
OBEu(J,K,bi,bj) = OBEu(J,K,bi,bj) + Tr_T
c OBEv(J,K,bi,bj) = 0.
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
#endif
#ifdef ALLOW_OBCS_WEST
Tr_T = 0. _d 0
Ar_T = 0. _d 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO J=1,sNy
I_obc=OB_Iw(J,bi,bj)
IF (I_obc.ne.0) THEN
Ar = drF(k)*hFacC(I_obc,j,k,bi,bj)*dyG(I_obc,j,bi,bj)
Ar_T = Ar_T + Ar
Tr_T = Tr_T + Ar * OBWu(J,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
_GLOBAL_SUM_RL( Ar_T , myThid )
IF ( Ar_T .GT. 0. _d 0 ) THEN
_GLOBAL_SUM_RL( Tr_T , myThid )
Tr_T = (0. _d 0 - Tr_T)/Ar_T
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO J=1-Oly,sNy+Oly
I_obc=OB_Iw(J,bi,bj)
IF (I_obc.ne.0) THEN
OBWu(J,K,bi,bj) = OBWu(J,K,bi,bj) + Tr_T
c OBWv(J,K,bi,bj) = 0.
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
#endif
#ifdef ALLOW_OBCS_NORTH
Tr_T = 0. _d 0
Ar_T = 0. _d 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO I=1,sNx
J_obc=OB_Jn(I,bi,bj)
IF (J_obc.ne.0) THEN
Ar = drF(k)*hFacC(i,J_obc,k,bi,bj)*dxG(i,J_obc,bi,bj)
Ar_T = Ar_T + Ar
Tr_T = Tr_T + Ar * OBNv(I,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
_GLOBAL_SUM_RL( Ar_T , myThid )
IF ( Ar_T .GT. 0. _d 0 ) THEN
_GLOBAL_SUM_RL( Tr_T , myThid )
Tr_T = (0. _d 0 - Tr_T)/Ar_T
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO I=1-Olx,sNx+Olx
J_obc=OB_Jn(I,bi,bj)
IF (J_obc.ne.0) THEN
c OBNu(I,K,bi,bj) = 0.
OBNv(I,K,bi,bj) = OBNv(I,K,bi,bj) + Tr_T
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
#endif
#ifdef ALLOW_OBCS_SOUTH
Tr_T = 0. _d 0
Ar_T = 0. _d 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO I=1,sNx
J_obc=OB_Js(I,bi,bj)
IF (J_obc.ne.0) THEN
Ar = drF(k)*hFacC(i,J_obc,k,bi,bj)*dxG(i,J_obc,bi,bj)
Ar_T = Ar_T + Ar
Tr_T = Tr_T + Ar * OBSv(I,K,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
_GLOBAL_SUM_RL( Ar_T , myThid )
IF ( Ar_T .GT. 0. _d 0 ) THEN
_GLOBAL_SUM_RL( Tr_T , myThid )
Tr_T = (0. _d 0 - Tr_T)/Ar_T
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO K=1,Nr
DO I=1-Olx,sNx+Olx
J_obc=OB_Js(I,bi,bj)
IF (J_obc.ne.0) THEN
c OBSu(I,K,bi,bj) = 0.
OBSv(I,K,bi,bj) = OBSv(I,K,bi,bj) + Tr_T
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
#endif
ENDIF
#endif /* ALLOW_OBCS_BALANCE */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid)
#endif
#endif /* ALLOW_OBCS */
RETURN
END