C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_calc_stevens.F,v 1.14 2015/01/12 16:06:38 mlosch Exp $
C $Name: $
#include "OBCS_OPTIONS.h"
#undef CHECK_BALANCE
C-- File obcs_calc_stevens.F:
C-- Contents
C-- o OBCS_CALC_STEVENS
C-- o OBCS_STEVENS_CALC_TRACER_EAST
C-- o OBCS_STEVENS_CALC_TRACER_WEST
C-- o OBCS_STEVENS_CALC_TRACER_NORTH
C-- o OBCS_STEVENS_CALC_TRACER_SOUTH
C-- o OBCS_STEVENS_SAVE_TRACER
CBOP
C !ROUTINE: OBCS_CALC_STEVENS
C !INTERFACE:
SUBROUTINE OBCS_CALC_STEVENS(
I futureTime, futureIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE OBCS_CALC_STEVENS
C | o Calculate future boundary data at open boundaries
C | at time = futureTime
C | from input data following Stevens(1990), and some
C | MOM3 legacy code
C |
C | o the code works like this
C | - the "barotropic" (= vertically averaged) velocity
C | normal to the boundary is assumed to be in
C | OBE/W/N/Su/v (normal) when this routine is entered
C | - the vertically averaged velocity is corrected
C | by the "baroclinic" (= deviation from vertically
C | averaged velocity) velocity to give a new OB?u/v;
C | the "barolinic" velocity is estimated from the previous
C | time step which makes this boundary condition depend on
C | a restart file. If OBCS_STEVENS_USE_INTERIOR_VELOCITY
C | is defined the velocity is simply copied from the model
C | interior to the boundary, thereby avoiding a restart
C | file or complicated reconstruction, but this solution
C | can give unexpected results.
C | (Note: in this context the terms barotropic and baroclinic
C | are MOM jargon and --- to my mind ---- should not be used)
C | - a wave phase speed is estimated from temporal and
C | horizontal variations of the tracer fields for each
C | tracer individually, this similar to Orlanski BCs,
C | but for simplicity the fields of the previous time step
C | are used, and the time derivative is estimated
C | independently of the time stepping procedure by simple
C | differencing
C | - velocity tangential to the boundary is always zero
C | (although this could be changed)
C | - a new tracer is computed from a local advection equation
C | with an upwind scheme: tracer from the interior is
C | advected out of the domain, and tracer from the boundary
C | is "advected" into the domain by a restoring mechanism
C | - for the advection equation only values from the
C | the current (not the updated) time level are used
C |
C *==========================================================*
C | Feb, 2009: started by Martin Losch (Martin.Losch@awi.de)
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "DYNVARS.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_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
_RL futureTime
INTEGER futureIter
INTEGER myThid
#ifdef ALLOW_OBCS_STEVENS
C !LOCAL VARIABLES:
C == Local variables ==
C I,J,K :: loop indices
C msgBuf :: Informational/error message buffer
C uMer/vZonBar :: vertically averaged velocity at open boundary
C drFBar :: local depth for vertical average
C uMer/vZonPri :: velocity anomalies applied to the open boundaries
C gammat/s :: restoring parameters (1./(T/SrelaxStevens - time scale))
C auxillary variables
C cflMer/Zon :: ratio of grid spacing and time step
C aFac :: switch (0 or 1) that turns on advective contribution
C gFacM/Z :: switch (0 or 1) that turns on restoring boundary condition
C pFac :: switch that turns on/off phase velocity contribution
INTEGER bi, bj
INTEGER I, J, K
c CHARACTER*(MAX_LEN_MBUF) msgBuf
_RL cflMer (1-OLy:sNy+OLy,Nr)
_RL gFacM (1-OLy:sNy+OLy,Nr)
_RL uMerPri(Nr)
_RL uMerBar
_RL cflZon (1-OLx:sNx+OLx,Nr)
_RL gFacZ (1-OLx:sNx+OLx,Nr)
_RL vZonPri(Nr)
_RL vZonBar
_RL drFBar
_RL gammat, gammas, pFac, aFac
#ifdef ALLOW_PTRACERS
c INTEGER iTracer
#endif /* ALLOW_PTRACERS */
#ifdef CHECK_BALANCE
_RL uVelLoc, vVelLoc
_RL vPhase
#endif
CEOP
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC_STEVENS',myThid)
#endif
aFac = 1. _d 0
IF (.NOT. useStevensAdvection ) aFac = 0. _d 0
pFac = 1. _d 0
IF (.NOT. useStevensPhaseVel ) pFac = 0. _d 0
gammat = 0. _d 0
IF (TrelaxStevens .GT. 0. _d 0) gammat = 1./TrelaxStevens
gammas = 0. _d 0
IF (SrelaxStevens .GT. 0. _d 0) gammas = 1./SrelaxStevens
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ikey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_OBCS_EAST
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBEt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBEs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
# endif
IF ( useStevensEast ) THEN
C Eastern OB
#ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC_STEVENS: East',myThid)
#endif
C compute vertical average and deviation from vertical
C average for I_obe
DO J=1-OLy,sNy+OLy
I = OB_Ie(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
C first initialize some fields
drFbar = 0. _d 0
uMerBar = 0. _d 0
DO K=1,Nr
uMerPri(K) = 0. _d 0
ENDDO
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
uMerBar = uMerBar + uVel(I-1,J,K,bi,bj)
#else
uMerBar = uMerBar + OBEuStevens(J,K,bi,bj)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& *drF(K)* _hFacW(I,J,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacW(I,J,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
uMerPri(K) = (uVel(I-1,J,K,bi,bj)-uMerBar)
#else
uMerPri(K) = (OBEuStevens(J,K,bi,bj)-uMerBar)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& * _maskW(I,J,K,bi,bj)
ENDDO
C vertical average of input field
drFbar = 0. _d 0
uMerBar = 0. _d 0
DO K=1,Nr
uMerBar = uMerBar + OBEu(J,K,bi,bj)
& *drF(K)* _hFacW(I,J,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacW(I,J,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
C Now the absolute velocity normal to the boundary is
C uMerBar + uMerPri(K).
DO K=1,Nr
OBEu(J,K,bi,bj) = (uMerBar + uMerPri(K))
& * _maskW(I,J,K,bi,bj)
CML OBEv(J,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
OBEw(J,K,bi,bj)=0.
#endif
ENDDO
ENDIF
ENDDO
#ifdef NONLIN_FRSURF
C this is a bit of a hack
IF ( nonlinFreeSurf.GT.0 ) THEN
DO J=1-OLy,sNy+OLy
I = OB_Ie(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
OBEeta(J,bi,bj) = etaN(I-1,J,bi,bj)
ENDIF
ENDDO
ENDIF
#endif /* NONLIN_FRSURF */
C Next, we compute the phase speed correction, which depends on the
C tracer!
DO K=1,Nr
DO J=1-OLy,sNy+OLy
I = OB_Ie(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
cflMer(J,K) = 0.5 _d 0 * _dxC(I-1,J,bi,bj)/dTtracerLev(K)
CML gFacM(J,K) = 0. _d 0
CML IF ( uVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacM(J,K) = 1. _d 0
gFacM(J,K) = ABS(MIN(SIGN(1.D0,uVel(I,J,K,bi,bj)),0.D0))
ELSE
cflMer(J,K) = 0. _d 0
gFacM (J,K) = 0. _d 0
ENDIF
ENDDO
ENDDO
C theta
CALL OBCS_STEVENS_CALC_TRACER_EAST(
U OBEt,
I OBEtStevens, theta, gammat,
I uVel, cflMer, gFacM, pFac, aFac,
I OB_Ie, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C salinity
CALL OBCS_STEVENS_CALC_TRACER_EAST(
U OBEs,
I OBEsStevens, salt, gammas,
I uVel, cflMer, gFacM, pFac, aFac,
I OB_Ie, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C Template for passive tracers, requires work
CML#ifdef ALLOW_PTRACERS
CMLC ptracers
CML IF ( usePtracers ) THEN
CML DO itracer = 1, PTRACERnumInUse
CML CALL OBCS_STEVENS_CALC_TRACER_EAST(
CML O OBEptr (1-OLy,1,1,1,iTracer),
CML I OBEpStevens (1-OLy,1,1,1,iTracer),
CML I pTracer(1-OLx,1-OLy,1,1,1,iTracer), gammas,
CML I uVel, cflMer, gFacM, pFac, aFac,
CML I OB_Ie, OB_indexNone, bi, bj,
CML I futureTime, futureIter,
CML I myThid )
CML ENDDO
CML ENDIF
CML#endif /* ALLOW_PTRACERS */
C IF ( useStevensEast ) THEN
ENDIF
#endif /* ALLOW_OBCS_EAST */
C ------------------------------------------------------------------------------
#ifdef ALLOW_OBCS_WEST
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBWt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBWs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
# endif
IF ( useStevensWest ) THEN
C Western OB
#ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC_STEVENS: West',myThid)
#endif
C compute vertical average and deviation from vertical
C average for I_obw+1
DO J=1-OLy,sNy+OLy
I = OB_Iw(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
C first initialize some fields
drFBar = 0. _d 0
uMerBar = 0. _d 0
DO K=1,Nr
uMerPri(K) = 0. _d 0
ENDDO
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
uMerBar = uMerBar + uVel(I+2,J,K,bi,bj)
#else
uMerBar = uMerBar + OBWuStevens(J,K,bi,bj)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& *drF(K)* _hFacW(I+1,J,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacW(I+1,J,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
uMerPri(K) = (uVel(I+2,J,K,bi,bj)-uMerBar)
#else
uMerPri(K) = (OBWuStevens(J,K,bi,bj)-uMerBar)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& * _maskW(I+1,J,K,bi,bj)
ENDDO
C vertical average of input field
drFBar = 0. _d 0
uMerBar = 0. _d 0
DO K=1,Nr
uMerBar = uMerBar + OBWu(J,K,bi,bj)
& *drF(K)* _hFacW(I+1,J,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacW(I+1,J,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
C Now the absolute velocity normal to the boundary is
C uMerBar + uMerPri(K).
DO K=1,Nr
OBWu(J,K,bi,bj) = (uMerBar + uMerPri(K))
& * _maskW(I+1,J,K,bi,bj)
CML OBWv(J,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
OBWw(J,K,bi,bj)=0.
#endif
ENDDO
ENDIF
ENDDO
#ifdef NONLIN_FRSURF
C this is a bit of a hack
IF ( nonlinFreeSurf.GT.0 ) THEN
DO J=1-OLy,sNy+OLy
I = OB_Iw(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
OBWeta(J,bi,bj) = etaN(I+1,J,bi,bj)
ENDIF
ENDDO
ENDIF
#endif /* NONLIN_FRSURF */
C Next, we compute the phase speed correction, which depends on the
C tracer!
DO K=1,Nr
DO J=1-OLy,sNy+OLy
I = OB_Iw(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
cflMer(J,K) = 0.5 _d 0 * _dxC(I+2,J,bi,bj)/dTtracerLev(K)
CML gFacM = 0. _d 0
CML IF ( uVel(I+1,J,K,bi,bj) .GT. 0. _d 0 ) gFacM = 1. _d 0
gFacM(J,K) = ABS(MAX(SIGN(1.D0,uVel(I+1,J,K,bi,bj)),0.D0))
ELSE
cflMer(J,K) = 0. _d 0
gFacM (J,K) = 0. _d 0
ENDIF
ENDDO
ENDDO
C theta
CALL OBCS_STEVENS_CALC_TRACER_WEST(
U OBWt,
I OBWtStevens, theta, gammat,
I uVel, cflMer, gFacM, pFac, aFac,
I OB_Iw, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C salinity
CALL OBCS_STEVENS_CALC_TRACER_WEST(
U OBWs,
I OBWsStevens, salt, gammas,
I uVel, cflMer, gFacM, pFac, aFac,
I OB_Iw, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C ptracers
C IF ( useStevensWest ) THEN
ENDIF
#endif /* ALLOW_OBCS_WEST */
C ------------------------------------------------------------------------------
#ifdef ALLOW_OBCS_NORTH
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBNt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBNs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
# endif
IF ( useStevensNorth ) THEN
C Northern OB
#ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC_STEVENS: North',myThid)
#endif
C compute vertical average and deviation from vertical
C average for J_obn
DO I=1-OLx,sNx+OLx
J = OB_Jn(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
C first initialize some fields
drFBar = 0. _d 0
vZonBar = 0. _d 0
DO K=1,Nr
vZonPri(K) = 0. _d 0
ENDDO
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
vZonBar = vZonBar + vVel(I,J-1,K,bi,bj)
#else
vZonBar = vZonBar + OBNvStevens(I,K,bi,bj)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& *drF(K)* _hFacS(I,J,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacS(I,J,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
vZonPri(K) = (vVel(I,J-1,K,bi,bj)-vZonBar)
#else
vZonPri(K) = (OBNvStevens(I,K,bi,bj)-vZonBar)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& * _maskS(I,J,K,bi,bj)
ENDDO
C vertical average of input field
drFBar = 0. _d 0
vZonBar = 0. _d 0
DO K=1,Nr
vZonBar = vZonBar + OBNv(I,K,bi,bj)
& *drF(K)* _hFacS(I,J,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacS(I,J,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
C Now the absolute velocity normal to the boundary is
C vZonBar + vZonPri(K).
DO K=1,Nr
OBNv(I,K,bi,bj) = (vZonBar + vZonPri(K))
& * _maskS(I,J,K,bi,bj)
CML OBNu(I,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
OBNw(I,K,bi,bj)=0.
#endif
ENDDO
ENDIF
ENDDO
#ifdef NONLIN_FRSURF
C this is a bit of a hack
IF ( nonlinFreeSurf.GT.0 ) THEN
DO I=1-OLx,sNx+OLx
J = OB_Jn(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
OBNeta(I,bi,bj) = etaN(I,J-1,bi,bj)
ENDIF
ENDDO
ENDIF
#endif /* NONLIN_FRSURF */
C Next, we compute the phase speed correction, which depends on the
C tracer!
DO K=1,Nr
DO I=1-OLx,sNx+OLx
J = OB_Jn(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
cflZon(I,K) = 0.5 _d 0 * _dyC(I,J-1,bi,bj)/dTtracerLev(K)
CML gFacZ(I,K) = 0. _d 0
CML IF ( vVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacZ(I,K) = 1. _d 0
gFacZ(I,K) = ABS(MIN(SIGN(1.D0,vVel(I,J,K,bi,bj)),0.D0))
ELSE
cflZon(I,K) = 0. _d 0
gFacZ (I,K) = 0. _d 0
ENDIF
ENDDO
ENDDO
C theta
CALL OBCS_STEVENS_CALC_TRACER_NORTH(
U OBNt,
I OBNtStevens, theta, gammat,
I vVel, cflZon, gFacZ, pFac, aFac,
I OB_Jn, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C salinity
CALL OBCS_STEVENS_CALC_TRACER_NORTH(
U OBNs,
I OBNsStevens, salt, gammas,
I vVel, cflZon, gFacZ, pFac, aFac,
I OB_Jn, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C ptracers
C IF ( useStevensNorth ) THEN
ENDIF
#endif /* ALLOW_OBCS_NORTH */
C ------------------------------------------------------------------------------
#ifdef ALLOW_OBCS_SOUTH
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBSt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBSs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
# endif
IF ( useStevensSouth ) THEN
C Southern OB
#ifdef ALLOW_DEBUG
IF (debugMode)
& CALL DEBUG_MSG('OBCS_CALC_STEVENS: South',myThid)
#endif
C compute vertical average and deviation from vertical
C average for J_obs+1
DO I=1-OLx,sNx+OLx
J = OB_Js(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
C first initialize some fields
drFBar = 0. _d 0
vZonBar = 0. _d 0
DO K=1,Nr
vZonPri(K) = 0. _d 0
ENDDO
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
vZonBar = vZonBar + vVel(I,J+2,K,bi,bj)
#else
vZonBar = vZonBar + OBSvStevens(I,K,bi,bj)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& *drF(K)* _hFacS(I,J+1,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacS(I,J+1,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
vZonPri(K) = (vVel(I,J+2,K,bi,bj)-vZonBar)
#else
vZonPri(K) = (OBSvStevens(I,K,bi,bj)-vZonBar)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
& * _maskS(I,J+1,K,bi,bj)
ENDDO
C vertical average of input field
drFBar = 0. _d 0
vZonBar = 0. _d 0
DO K=1,Nr
vZonBar = vZonBar + OBSv(I,K,bi,bj)
& *drF(K)* _hFacS(I,J+1,K,bi,bj)
drFBar = drFBar + drF(K)* _hFacS(I,J+1,K,bi,bj)
ENDDO
IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
C Now the absolute velocity normal to the boundary is
C vZonBar + vZonPri(K).
DO K=1,Nr
OBSv(I,K,bi,bj) = (vZonBar + vZonPri(K))
& * _maskS(I,J+1,K,bi,bj)
CML OBSu(I,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
OBSw(I,K,bi,bj)=0.
#endif
ENDDO
ENDIF
ENDDO
#ifdef NONLIN_FRSURF
C this is a bit of a hack
IF ( nonlinFreeSurf.GT.0 ) THEN
DO I=1-OLx,sNx+OLx
J = OB_Js(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
OBSeta(I,bi,bj) = etaN(I,J+1,bi,bj)
ENDIF
ENDDO
ENDIF
#endif /* NONLIN_FRSURF */
C Next, we compute the phase speed correction, which depends on the
C tracer!
DO K=1,Nr
DO I=1-OLx,sNx+OLx
J = OB_Js(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
cflZon(I,K) = 0.5 _d 0 * _dyC(I,J+2,bi,bj)/dTtracerLev(K)
CML gFacZ = 0. _d 0
CML IF ( vVel(I,J+1,K,bi,bj) .GT. 0. _d 0 ) gFacZ = 1. _d 0
gFacZ(I,K) = ABS(MAX(SIGN(1.D0,vVel(I,J+1,K,bi,bj)),0.D0))
ELSE
cflZon(I,K) = 0. _d 0
gFacZ (I,K) = 0. _d 0
ENDIF
ENDDO
ENDDO
C theta
CALL OBCS_STEVENS_CALC_TRACER_SOUTH(
U OBSt,
I OBStStevens, theta, gammat,
I vVel, cflZon, gFacZ, pFac, aFac,
I OB_Js, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C salinity
CALL OBCS_STEVENS_CALC_TRACER_SOUTH(
U OBSs,
I OBSsStevens, salt, gammas,
I vVel, cflZon, gFacZ, pFac, aFac,
I OB_Js, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C ptracers
C IF ( useStevensSouth ) THEN
ENDIF
#endif /* ALLOW_OBCS_SOUTH */
C end bi/bj-loops
ENDDO
ENDDO
C save the tracer fields of the previous time step for the next time step
CALL OBCS_STEVENS_SAVE_TRACERS(
I futureTime, futureIter,
I myThid )
C ------------------------------------------------------------------------------
#ifdef CHECK_BALANCE
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
uPhase=0.
vPhase=0.
uVelLoc = 0.
DO J=1-OLy,sNy+OLy
uMerBar=0. _d 0
DO K=1,Nr
I = OB_Ie(J,bi,bj)
IF ( I.EQ.OB_indexNone ) I = 1
uPhase = uPhase + OBEu(J,K,bi,bj)
& *drF(k)* _hFacW(I,J,K,bi,bj)*dyG(I,J,bi,bj)
I = OB_Iw(J,bi,bj)
IF ( I.EQ.OB_indexNone ) I = 1
vPhase = vPhase + OBWu(J,K,bi,bj)
& *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
CML uVelLoc = uVelLoc + uMerPri(J,K)
CML & *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
CML uMerBar(J)=uMerBar(J) + uMerPri(J,K)
CML & *drF(k)* _hFacW(I+1,J,K,bi,bj)
ENDDO
CML print *, 'ml-obcs: uBar = ', j,uMerBar(J)
ENDDO
C end bi/bj-loops
ENDDO
ENDDO
_GLOBAL_SUM_RL( uPhase, myThid )
_GLOBAL_SUM_RL( vPhase, myThid )
CML _GLOBAL_SUM_RL( uVelLoc, myThid )
print *, 'ml-obcs: OBE = ', uPhase*1 _d -6, ' Sv'
print *, 'ml-obcs: OBW = ', vPhase*1 _d -6, ' Sv'
CML print *, 'ml-obcs: OBWp = ', uVelLoc*1 _d -6, ' Sv'
#endif /* CHECK_BALANCE */
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC_STEVENS',myThid)
#endif
#endif /* ALLOW_OBCS_STEVENS */
RETURN
END
CBOP
C !ROUTINE: OBCS_STEVENS_CALC_TRACER_EAST
C !INTERFACE:
SUBROUTINE OBCS_STEVENS_CALC_TRACER_EAST(
U OBEf,
I OBE_Stevens, tracer, gammaf,
I uVel, cflMer, gFacM, pFac, aFac,
I OB_I, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE OBCS_STEVENS_CALC_TRACER_EAST
C | Calculate tracer value at the eastern OB location
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id number
C bi, bj :: indices of current tile
_RL futureTime
INTEGER futureIter
INTEGER myThid
INTEGER bi, bj
INTEGER OB_indexNone
INTEGER OB_I (1-OLy:sNy+OLy,nSx,nSy)
_RL cflMer (1-OLy:sNy+OLy,Nr)
_RL gFacM (1-OLy:sNy+OLy,Nr)
_RL gammaf, pFac, aFac
_RL OBEf (1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL OBE_Stevens (1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
#ifdef ALLOW_OBCS_STEVENS
C !LOCAL VARIABLES:
C i,j,k :: loop indices
C uPhase :: estimate of phase velocity for radiation condition
C dtracSpace :: horizontal difference of tracer
C dtracTime :: temporal difference of tracer
INTEGER i,j,k
_RL uPhase
_RL dtracSpace
_RL dTracTime
CEOP
DO K=1,Nr
DO J=1-OLy,sNy+OLy
I = OB_I(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
dTracSpace = (tracer(I-1,J,K,bi,bj)-tracer(I-2,J,K,bi,bj))
& * _maskW(I-1,J,K,bi,bj)
dTracTime = (tracer(I-1,J,K,bi,bj)-OBE_Stevens(J,K,bi,bj))
uPhase = cflMer(J,K) * pFac
IF ( dTracSpace .NE. 0. _d 0 ) THEN
uPhase = MIN( cflMer(J,K),
& MAX( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
& ) * pFac
ENDIF
C Compute the tracer tendency here, the tracer will be updated
C with a simple Euler forward step in S/R obcs_apply_ts
OBEf(J,K,bi,bj) = _maskW(I,J,K,bi,bj) * (
& - ( aFac*MAX(0.D0,uVel(I,J,K,bi,bj)) + uPhase )
& *(tracer(I,J,K,bi,bj)-tracer(I-1,J,K,bi,bj))
& * _recip_dxC(I,J,bi,bj)
& - gFacM(J,K) * gammaf
& * (tracer(I,J,K,bi,bj)-OBEf(J,K,bi,bj)) )
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_OBCS_STEVENS */
RETURN
END
CBOP
C !ROUTINE: OBCS_STEVENS_CALC_TRACER_WEST
C !INTERFACE:
SUBROUTINE OBCS_STEVENS_CALC_TRACER_WEST(
U OBWf,
I OBW_Stevens, tracer, gammaf,
I uVel, cflMer, gFacM, pFac, aFac,
I OB_I, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE OBCS_STEVENS_CALC_TRACER_WEST
C | Calculate tracer value at the western OB location
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id number
C bi, bj :: indices of current tile
_RL futureTime
INTEGER futureIter
INTEGER myThid
INTEGER bi, bj
INTEGER OB_indexNone
INTEGER OB_I (1-OLy:sNy+OLy,nSx,nSy)
_RL cflMer (1-OLy:sNy+OLy,Nr)
_RL gFacM (1-OLy:sNy+OLy,Nr)
_RL gammaf, pFac, aFac
_RL OBWf (1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL OBW_Stevens (1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
#ifdef ALLOW_OBCS_STEVENS
C !LOCAL VARIABLES:
C i,j,k :: loop indices
C uPhase :: estimate of phase velocity for radiation condition
C dtracSpace :: horizontal difference of tracer
C dtracTime :: temporal difference of tracer
INTEGER i,j,k
_RL uPhase
_RL dtracSpace
_RL dTracTime
CEOP
DO K=1,Nr
DO J=1-OLy,sNy+OLy
I = OB_I(J,bi,bj)
IF ( I.NE.OB_indexNone ) THEN
dTracSpace = (tracer(I+2,J,K,bi,bj)-tracer(I+1,J,K,bi,bj))
& * _maskW(I+2,J,K,bi,bj)
dTracTime = (tracer(I+1,J,K,bi,bj)-OBW_Stevens(J,K,bi,bj))
uPhase = -cflMer(J,K) * pFac
IF ( dTracSpace .NE. 0. _d 0 ) THEN
uPhase = MAX( -cflMer(J,K),
& MIN( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
& ) * pFac
ENDIF
C Compute the tracer tendency here, the tracer will be updated
C with a simple Euler forward step in S/R obcs_apply_ts
OBWf(J,K,bi,bj) = _maskW(I+1,J,K,bi,bj) * (
& - ( aFac*MIN(0.D0,uVel(I+1,J,K,bi,bj)) + uPhase )
& *(tracer(I+1,J,K,bi,bj)-tracer(I,J,K,bi,bj))
& * _recip_dxC(I+1,J,bi,bj)
& - gFacM(J,K) * gammaf
& * (tracer(I,J,K,bi,bj)-OBWf(J,K,bi,bj)) )
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_OBCS_STEVENS */
RETURN
END
CBOP
C !ROUTINE: OBCS_STEVENS_CALC_TRACER_NORTH
C !INTERFACE:
SUBROUTINE OBCS_STEVENS_CALC_TRACER_NORTH(
U OBNf,
I OBN_Stevens, tracer, gammaf,
I vVel, cflZon, gFacZ, pFac, aFac,
I OB_J, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE OBCS_STEVENS_CALC_TRACER_NORTH
C | Calculate tracer value at the northern OB location
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id number
C bi, bj :: indices of current tile
_RL futureTime
INTEGER futureIter
INTEGER myThid
INTEGER bi, bj
INTEGER OB_indexNone
INTEGER OB_J (1-OLx:sNx+OLx,nSx,nSy)
_RL cflZon (1-OLx:sNx+OLx,Nr)
_RL gFacZ (1-OLx:sNx+OLx,Nr)
_RL gammaf, pFac, aFac
_RL OBNf (1-Olx:sNx+Olx,Nr,nSx,nSy)
_RL OBN_Stevens (1-Olx:sNx+Olx,Nr,nSx,nSy)
_RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
#ifdef ALLOW_OBCS_STEVENS
C !LOCAL VARIABLES:
C i,j,k :: loop indices
C vPhase :: estimate of phase velocity for radiation condition
C dtracSpace :: horizontal difference of tracer
C dtracTime :: temporal difference of tracer
INTEGER i,j,k
_RL vPhase
_RL dtracSpace
_RL dTracTime
CEOP
DO K=1,Nr
DO I=1-OLx,sNx+OLx
J = OB_J(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
C Theta first:
dTracSpace = (tracer(I,J-1,K,bi,bj)-tracer(I,J-2,K,bi,bj))
& * _maskS(I,J-1,K,bi,bj)
dTracTime = (tracer(I,J-1,K,bi,bj)-OBN_Stevens(I,K,bi,bj))
vPhase = cflZon(I,K) * pFac
IF ( dTracSpace .NE. 0. _d 0 ) THEN
vPhase = MIN( cflZon(I,K),
& MAX( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
& ) * pFac
ENDIF
C Compute the tracer tendency here, the tracer will be updated
C with a simple Euler forward step in S/R obcs_apply_ts
OBNf(I,K,bi,bj) = _maskS(I,J,K,bi,bj) * (
& - ( aFac*MAX(0.D0,vVel(I,J,K,bi,bj)) + vPhase )
& *(tracer(I,J,K,bi,bj)-tracer(I,J-1,K,bi,bj))
& * _recip_dyC(I,J,bi,bj)
& - gFacZ(I,K) * gammaf
& * (tracer(I,J,K,bi,bj)-OBNf(I,K,bi,bj)) )
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_OBCS_STEVENS */
RETURN
END
CBOP
C !ROUTINE: OBCS_STEVENS_CALC_TRACER_SOUTH
C !INTERFACE:
SUBROUTINE OBCS_STEVENS_CALC_TRACER_SOUTH(
U OBSf,
I OBS_Stevens, tracer, gammaf,
I vVel, cflZon, gFacZ, pFac, aFac,
I OB_J, OB_indexNone, bi, bj,
I futureTime, futureIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE OBCS_STEVENS_CALC_TRACER_SOUTH
C | Calculate tracer value at the southern OB location
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id number
C bi, bj :: indices of current tile
_RL futureTime
INTEGER futureIter
INTEGER myThid
INTEGER bi, bj
INTEGER OB_indexNone
INTEGER OB_J (1-OLx:sNx+OLx,nSx,nSy)
_RL cflZon (1-OLx:sNx+OLx,Nr)
_RL gFacZ (1-OLx:sNx+OLx,Nr)
_RL gammaf, pFac, aFac
_RL OBSf (1-Olx:sNx+Olx,Nr,nSx,nSy)
_RL OBS_Stevens (1-Olx:sNx+Olx,Nr,nSx,nSy)
_RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
#ifdef ALLOW_OBCS_STEVENS
C !LOCAL VARIABLES:
C i,j,k :: loop indices
C vPhase :: estimate of phase velocity for radiation condition
C dtracSpace :: horizontal difference of tracer
C dtracTime :: temporal difference of tracer
INTEGER i,j,k
_RL vPhase
_RL dtracSpace
_RL dTracTime
CEOP
DO K=1,Nr
DO I=1-OLx,sNx+OLx
J = OB_J(I,bi,bj)
IF ( J.NE.OB_indexNone ) THEN
dTracSpace = (tracer(I,J+2,K,bi,bj)-tracer(I,J+1,K,bi,bj))
& * _maskS(I,J+2,K,bi,bj)
dTracTime = (tracer(I,J+1,K,bi,bj)-OBS_Stevens(I,K,bi,bj))
vPhase = -cflZon(I,K) * pFac
IF ( dTracSpace .NE. 0. _d 0 ) THEN
vPhase = MAX( -cflZon(I,K),
& MIN( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
& ) * pFac
ENDIF
C Compute the tracer tendency here, the tracer will be updated
C with a simple Euler forward step in S/R obcs_apply_ts
OBSf(I,K,bi,bj) = _maskS(I,J+1,K,bi,bj) * (
& - ( aFac*MIN(0.D0,vVel(I,J+1,K,bi,bj)) + vPhase )
& *(tracer(I,J+1,K,bi,bj)-tracer(I,J,K,bi,bj))
& * _recip_dyC(I,J+1,bi,bj)
& - gFacZ(I,K) * gammaf
& * (tracer(I,J,K,bi,bj)-OBSf(I,K,bi,bj)) )
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_OBCS_STEVENS */
RETURN
END
CBOP
C !ROUTINE: OBCS_STEVENS_SAVE_TRACERS
C !INTERFACE:
SUBROUTINE OBCS_STEVENS_SAVE_TRACERS(
I futureTime, futureIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE OBCS_STEVENS_SAVE_TRACERS
C | Save tracers (of previous time step) at the OB location
C | to be used in the next time step for Stevens boundary
C | conditions
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "DYNVARS.h"
CML#ifdef ALLOW_PTRACERS
CML#include "PTRACERS_SIZE.h"
CML#include "PTRACERS_PARAMS.h"
CML#include "PTRACERS_FIELDS.h"
CML#include "OBCS_PTRACERS.h"
CML#endif /* ALLOW_PTRACERS */
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id number
_RL futureTime
INTEGER futureIter
INTEGER myThid
#ifdef ALLOW_OBCS_STEVENS
C !LOCAL VARIABLES:
C bi, bj :: indices of current tile
C i,j,k :: loop indices
C Iobc,Jobc :: position-index of open boundary
INTEGER bi,bj,i,j,k,Iobc,Jobc
CEOP
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_OBCS_NORTH
IF ( tileHasOBN(bi,bj) .AND. useStevensNorth ) THEN
C Northern boundary
DO i=1-OLx,sNx+OLx
Jobc = OB_Jn(i,bi,bj)
IF ( Jobc.NE.OB_indexNone ) THEN
DO k = 1,Nr
OBNtStevens(i,k,bi,bj) = theta(i,Jobc-1,k,bi,bj)
& *maskC(i,Jobc+1,k,bi,bj)
OBNsStevens(i,k,bi,bj) = salt(i,Jobc-1,k,bi,bj)
& *maskC(i,Jobc+1,k,bi,bj)
ENDDO
ENDIF
ENDDO
ENDIF
#endif /* ALLOW_OBCS_NORTH */
#ifdef ALLOW_OBCS_SOUTH
IF ( tileHasOBS(bi,bj) .AND. useStevensSouth ) THEN
C Southern boundary
DO i=1-OLx,sNx+OLx
Jobc = OB_Js(i,bi,bj)
IF ( Jobc.NE.OB_indexNone ) THEN
DO k = 1,Nr
OBStStevens(i,k,bi,bj) = theta(i,Jobc+1,k,bi,bj)
& *maskC(i,Jobc+1,k,bi,bj)
OBSsStevens(i,k,bi,bj) = salt(i,Jobc+1,k,bi,bj)
& *maskC(i,Jobc+1,k,bi,bj)
ENDDO
ENDIF
ENDDO
ENDIF
#endif /* ALLOW_OBCS_SOUTH */
#ifdef ALLOW_OBCS_EAST
IF ( tileHasOBE(bi,bj) .AND. useStevensEast ) THEN
C Eastern boundary
DO j=1-OLy,sNy+OLy
Iobc = OB_Ie(j,bi,bj)
IF ( Iobc.NE.OB_indexNone ) THEN
DO k = 1,Nr
OBEtStevens(j,k,bi,bj) = theta(Iobc-1,j,k,bi,bj)
& *maskC(Iobc-1,j,k,bi,bj)
OBEsStevens(j,k,bi,bj) = salt(Iobc-1,j,k,bi,bj)
& *maskC(Iobc-1,j,k,bi,bj)
ENDDO
ENDIF
ENDDO
ENDIF
#endif /* ALLOW_OBCS_EAST */
#ifdef ALLOW_OBCS_WEST
IF ( tileHasOBW(bi,bj) .AND. useStevensWest ) THEN
C Western boundary
DO j=1-OLy,sNy+OLy
Iobc = OB_Iw(j,bi,bj)
IF ( Iobc.NE.OB_indexNone ) THEN
DO k = 1,Nr
OBWtStevens(j,k,bi,bj) = theta(Iobc+1,j,k,bi,bj)
& *maskC(Iobc+1,j,k,bi,bj)
OBWsStevens(j,k,bi,bj) = salt(Iobc+1,j,k,bi,bj)
& *maskC(Iobc+1,j,k,bi,bj)
ENDDO
ENDIF
ENDDO
ENDIF
#endif /* ALLOW_OBCS_WEST */
ENDDO
ENDDO
#endif /* ALLOW_OBCS_STEVENS */
RETURN
END