C $Header: /u/gcmpack/MITgcm/verification/internal_wave/code/obcs_calc.F,v 1.9 2011/12/12 19:04:25 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 *==========================================================*
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#ifdef ALLOW_EXCH2
# include "W2_EXCH2_SIZE.h"
#endif /* ALLOW_EXCH2 */
#include "SET_GRID.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "EOS.h"
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 ==
INTEGER bi, bj
INTEGER I, J ,K
_RL obTimeScale,Uinflow,rampTime2
_RL vertStructWst(Nr)
_RL mz,strat,kx
_RL tmpsum
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid)
#endif
C Vertical mode number
mz=1.0 _d 0
C Stratification
strat = 1.0 _d -6 / (gravity*tAlpha)
C Create a vertical structure function with zero mean
tmpsum=0.
do K=1,Nr
vertStructWst(K)=cos(mz*PI* (rC(K)/rF(Nr+1)) )
tmpsum=tmpsum+vertStructWst(K)*drF(K)
enddo
tmpsum=tmpsum/rF(Nr+1)
do K=1,Nr
vertStructWst(K)=vertStructWst(K)-tmpsum
enddo
c
obTimeScale = 44567.0 _d 0
kx=mz*2. _d 0*pi/400.0 _d 0
& *sqrt((2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale)
& - f0*f0)/(1.0 _d -6
& - 2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale)))
Uinflow = 0.024 _d 0
C *NOTE* I have commented out the ramp function below
C just to speed things up. You will probably want to use it
C for smoother looking solutions.
rampTime2 = 4. _d 0*44567.0 _d 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
C Eastern OB
IF (useOrlanskiEast) THEN
CALL ORLANSKI_EAST(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=1,Nr
DO J=1-Oly,sNy+Oly
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
ENDDO
ENDDO
ENDIF
C Western OB
IF (useOrlanskiWest) THEN
CALL ORLANSKI_WEST(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=1,Nr
DO J=1-Oly,sNy+Oly
OBWu(J,K,bi,bj)=0. _d 0
& +Uinflow
& *vertStructWst(K)
& *sin(2. _d 0*PI*futureTime/obTimeScale)
c & *(exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
& *cos(kx*(3. _d 0-2. _d 0-0.5 _d 0)*delX(1))
OBWv(J,K,bi,bj)=0. _d 0
& +Uinflow
& *f0/(2.0 _d 0*PI/obTimeScale)
& *vertStructWst(K)
& *cos(2. _d 0*PI*futureTime/obTimeScale )
& * (exp(futureTime/rampTime2)
& - exp(-futureTime/rampTime2))
& /(exp(futureTime/rampTime2)
& + exp(-futureTime/rampTime2))
OBWt(J,K,bi,bj)=tRef(K)
& + Uinflow*sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
& * sin(2.0 _d 0*PI*futureTime/obTimeScale)
& *sqrt(strat/(tAlpha*gravity))
& *sqrt(2.0 _d 0*PI/obTimeScale*2.0*PI/obTimeScale - f0*f0)
& /(2.0 _d 0*PI/obTimeScale)
c & * (exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
#ifdef ALLOW_NONHYDROSTATIC
OBWw(J,K,bi,bj)=-Uinflow
& *sqrt(2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale - f0*f0)
& /sqrt(strat*strat -
& 2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale)
& *sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
& *cos(2. _d 0*PI*futureTime/obTimeScale)
c & *(exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
#endif
ENDDO
ENDDO
ENDIF
C Northern OB, template for forcing
IF (useOrlanskiNorth) THEN
CALL ORLANSKI_NORTH(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=1,Nr
DO I=1-Olx,sNx+Olx
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
ENDDO
ENDDO
ENDIF
C Southern OB, template for forcing
IF (useOrlanskiSouth) THEN
CALL ORLANSKI_SOUTH(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=1,Nr
DO I=1-Olx,sNx+Olx
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
ENDDO
ENDDO
ENDIF
C-- end bi,bj loops.
ENDDO
ENDDO
#ifdef ALLOW_OBCS_BALANCE
IF ( useOBCSbalance ) THEN
CALL OBCS_BALANCE_FLOW( futureTime, futureIter, myThid )
ENDIF
#endif /* ALLOW_OBCS_BALANCE */
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid)
#endif
#endif /* ALLOW_OBCS */
RETURN
END