internal_wave/code/obcs_calc.F

42c525bfb4 Alis*0001 #include "OBCS_OPTIONS.h"
                0002
3a86c9b47d Oliv*0003       SUBROUTINE OBCS_CALC( futureTime, futureIter,
2b169ecc4b Jean*0004      &                      uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0005      &                      myThid )
2b169ecc4b Jean*0006 C     *==========================================================*
                0007 C     | SUBROUTINE OBCS_CALC
                0008 C     | o Calculate future boundary data at open boundaries
                0009 C     |   at time = futureTime
                0010 C     *==========================================================*
42c525bfb4 Alis*0011       IMPLICIT NONE
                0012
                0013 C     === Global variables ===
                0014 #include "SIZE.h"
                0015 #include "EEPARAMS.h"
                0016 #include "PARAMS.h"
b9b591469d Jean*0017 #ifdef ALLOW_EXCH2
                0018 # include "W2_EXCH2_SIZE.h"
                0019 #endif /* ALLOW_EXCH2 */
                0020 #include "SET_GRID.h"
2b169ecc4b Jean*0021 #include "GRID.h"
87cd09b3a4 Jean*0022 #include "OBCS_PARAMS.h"
                0023 #include "OBCS_GRID.h"
                0024 #include "OBCS_FIELDS.h"
2b169ecc4b Jean*0025 #include "EOS.h"
42c525bfb4 Alis*0026
                0027 C     == Routine arguments ==
683a18afb1 Jean*0028       INTEGER futureIter
42c525bfb4 Alis*0029       _RL futureTime
                0030       _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0031       _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0032       _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0033       _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0034       _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0035       INTEGER myThid
                0036
                0037 #ifdef ALLOW_OBCS
                0038
                0039 C     == Local variables ==
3a86c9b47d Oliv*0040       INTEGER bi, bj
42c525bfb4 Alis*0041       INTEGER I, J ,K
                0042       _RL obTimeScale,Uinflow,rampTime2
                0043       _RL vertStructWst(Nr)
                0044       _RL mz,strat,kx
                0045       _RL tmpsum
                0046
2b169ecc4b Jean*0047 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0048
                0049 #ifdef ALLOW_DEBUG
                0050       IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid)
                0051 #endif
                0052
42c525bfb4 Alis*0053 C Vertical mode number
e7832df42a Alis*0054       mz=1.0 _d 0
42c525bfb4 Alis*0055 C Stratification
                0056       strat = 1.0 _d -6 / (gravity*tAlpha)
                0057
                0058 C Create a vertical structure function with zero mean
                0059       tmpsum=0.
                0060       do K=1,Nr
                0061        vertStructWst(K)=cos(mz*PI* (rC(K)/rF(Nr+1)) )
                0062        tmpsum=tmpsum+vertStructWst(K)*drF(K)
                0063       enddo
                0064       tmpsum=tmpsum/rF(Nr+1)
                0065       do K=1,Nr
                0066        vertStructWst(K)=vertStructWst(K)-tmpsum
                0067       enddo
                0068 c
e7832df42a Alis*0069       obTimeScale = 44567.0 _d 0
                0070        kx=mz*2. _d 0*pi/400.0 _d 0
                0071      &  *sqrt((2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale)
42c525bfb4 Alis*0072      & - f0*f0)/(1.0 _d -6
e7832df42a Alis*0073      & - 2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale)))
                0074       Uinflow = 0.024 _d 0
                0075 C *NOTE* I have commented out the ramp function below
                0076 C just to speed things up. You will probably want to use it
                0077 C for smoother looking solutions.
                0078       rampTime2 = 4. _d 0*44567.0 _d 0
42c525bfb4 Alis*0079
3a86c9b47d Oliv*0080       DO bj=myByLo(myThid),myByHi(myThid)
                0081       DO bi=myBxLo(myThid),myBxHi(myThid)
42c525bfb4 Alis*0082
                0083 C     Eastern OB
                0084       IF (useOrlanskiEast) THEN
                0085         CALL ORLANSKI_EAST(
2b169ecc4b Jean*0086      &          bi, bj, futureTime,
                0087      &          uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0088      &          myThid )
                0089       ELSE
                0090         DO K=1,Nr
                0091           DO J=1-Oly,sNy+Oly
                0092             OBEu(J,K,bi,bj)=0.
                0093             OBEv(J,K,bi,bj)=0.
                0094             OBEt(J,K,bi,bj)=tRef(K)
                0095             OBEs(J,K,bi,bj)=sRef(K)
                0096 #ifdef ALLOW_NONHYDROSTATIC
                0097             OBEw(J,K,bi,bj)=0.
                0098 #endif
                0099           ENDDO
                0100         ENDDO
                0101       ENDIF
                0102
                0103 C     Western OB
                0104       IF (useOrlanskiWest) THEN
                0105         CALL ORLANSKI_WEST(
2b169ecc4b Jean*0106      &          bi, bj, futureTime,
                0107      &          uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0108      &          myThid )
                0109       ELSE
                0110         DO K=1,Nr
                0111           DO J=1-Oly,sNy+Oly
e7832df42a Alis*0112           OBWu(J,K,bi,bj)=0. _d 0
42c525bfb4 Alis*0113      &       +Uinflow
                0114      &       *vertStructWst(K)
e7832df42a Alis*0115      &       *sin(2. _d 0*PI*futureTime/obTimeScale)
                0116 c    &       *(exp(futureTime/rampTime2)
                0117 c    &   - exp(-futureTime/rampTime2))
                0118 c    &   /(exp(futureTime/rampTime2)
                0119 c    &  + exp(-futureTime/rampTime2))
                0120      &   *cos(kx*(3. _d 0-2. _d 0-0.5 _d 0)*delX(1))
                0121           OBWv(J,K,bi,bj)=0. _d 0
42c525bfb4 Alis*0122      &       +Uinflow
e7832df42a Alis*0123      &       *f0/(2.0 _d 0*PI/obTimeScale)
42c525bfb4 Alis*0124      &       *vertStructWst(K)
e7832df42a Alis*0125      &       *cos(2. _d 0*PI*futureTime/obTimeScale )
42c525bfb4 Alis*0126      & * (exp(futureTime/rampTime2)
                0127      &   - exp(-futureTime/rampTime2))
                0128      &   /(exp(futureTime/rampTime2)
                0129      &  + exp(-futureTime/rampTime2))
                0130           OBWt(J,K,bi,bj)=tRef(K)
e7832df42a Alis*0131      & + Uinflow*sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
                0132      & * sin(2.0 _d 0*PI*futureTime/obTimeScale)
42c525bfb4 Alis*0133      & *sqrt(strat/(tAlpha*gravity))
e7832df42a Alis*0134      & *sqrt(2.0 _d 0*PI/obTimeScale*2.0*PI/obTimeScale - f0*f0)
                0135      & /(2.0 _d 0*PI/obTimeScale)
                0136 c    & * (exp(futureTime/rampTime2)
                0137 c    &   - exp(-futureTime/rampTime2))
                0138 c    &   /(exp(futureTime/rampTime2)
                0139 c    &  + exp(-futureTime/rampTime2))
42c525bfb4 Alis*0140 #ifdef ALLOW_NONHYDROSTATIC
                0141           OBWw(J,K,bi,bj)=-Uinflow
e7832df42a Alis*0142      & *sqrt(2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale - f0*f0)
                0143      & /sqrt(strat*strat -
                0144      &          2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale)
                0145      & *sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
                0146      &       *cos(2. _d 0*PI*futureTime/obTimeScale)
                0147 c    &       *(exp(futureTime/rampTime2)
                0148 c    &   - exp(-futureTime/rampTime2))
                0149 c    &   /(exp(futureTime/rampTime2)
                0150 c    &  + exp(-futureTime/rampTime2))
42c525bfb4 Alis*0151 #endif
                0152           ENDDO
                0153         ENDDO
                0154       ENDIF
                0155
                0156 C         Northern OB, template for forcing
                0157       IF (useOrlanskiNorth) THEN
                0158         CALL ORLANSKI_NORTH(
2b169ecc4b Jean*0159      &          bi, bj, futureTime,
                0160      &          uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0161      &          myThid )
                0162       ELSE
                0163         DO K=1,Nr
                0164           DO I=1-Olx,sNx+Olx
                0165             OBNv(I,K,bi,bj)=0.
                0166             OBNu(I,K,bi,bj)=0.
                0167             OBNt(I,K,bi,bj)=tRef(K)
                0168             OBNs(I,K,bi,bj)=sRef(K)
                0169 #ifdef ALLOW_NONHYDROSTATIC
                0170             OBNw(I,K,bi,bj)=0.
                0171 #endif
                0172           ENDDO
                0173         ENDDO
                0174       ENDIF
                0175
                0176 C         Southern OB, template for forcing
2b169ecc4b Jean*0177       IF (useOrlanskiSouth) THEN
42c525bfb4 Alis*0178         CALL ORLANSKI_SOUTH(
2b169ecc4b Jean*0179      &          bi, bj, futureTime,
                0180      &          uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0181      &          myThid )
                0182       ELSE
                0183         DO K=1,Nr
                0184           DO I=1-Olx,sNx+Olx
                0185             OBSu(I,K,bi,bj)=0.
                0186             OBSv(I,K,bi,bj)=0.
                0187             OBSt(I,K,bi,bj)=tRef(K)
                0188             OBSs(I,K,bi,bj)=sRef(K)
                0189 #ifdef ALLOW_NONHYDROSTATIC
                0190             OBSw(I,K,bi,bj)=0.
                0191 #endif
                0192           ENDDO
                0193         ENDDO
                0194       ENDIF
                0195
3a86c9b47d Oliv*0196 C--   end bi,bj loops.
                0197       ENDDO
2b169ecc4b Jean*0198       ENDDO
                0199
                0200 #ifdef ALLOW_DEBUG
                0201       IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid)
                0202 #endif
42c525bfb4 Alis*0203 #endif /* ALLOW_OBCS */
2b169ecc4b Jean*0204
42c525bfb4 Alis*0205       RETURN
                0206       END