Back to home page

MITgcm

 
 

    


File indexing completed on 2024-05-14 05:10:40 UTC

view on githubraw file Latest commit 9e0fdb40 on 2024-05-13 20:05:25 UTC
42c525bfb4 Alis*0001 #include "OBCS_OPTIONS.h"
                0002 
a24437ad80 Jean*0003       SUBROUTINE ORLANSKI_EAST( bi, bj, futureTime,
                0004      I                      uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0005      I                      myThid )
9e0fdb4065 Garr*0006 C     *==========================================================*
42c525bfb4 Alis*0007 C     | SUBROUTINE ORLANSKI_EAST                                 |
                0008 C     | o Calculate future boundary data at open boundaries      |
                0009 C     |   at time = futureTime by applying Orlanski radiation    |
                0010 C     |   conditions.                                            |
9e0fdb4065 Garr*0011 C     *==========================================================*
42c525bfb4 Alis*0012       IMPLICIT NONE
                0013 
                0014 C     === Global variables ===
                0015 #include "SIZE.h"
                0016 #include "EEPARAMS.h"
                0017 #include "PARAMS.h"
                0018 #include "GRID.h"
9b4f2a04e2 Jean*0019 #include "OBCS_PARAMS.h"
                0020 #include "OBCS_GRID.h"
                0021 #include "OBCS_FIELDS.h"
42c525bfb4 Alis*0022 #include "ORLANSKI.h"
                0023 
a24437ad80 Jean*0024 C SPK 6/2/00: Added radiative OBCs for salinity.
                0025 C SPK 6/6/00: Changed calculation of OB*w. When K=1, the
42c525bfb4 Alis*0026 C             upstream value is used. For example on the eastern OB:
                0027 C                IF (K.EQ.1) THEN
                0028 C                   OBEw(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
                0029 C                ENDIF
a24437ad80 Jean*0030 C
42c525bfb4 Alis*0031 C SPK 7/7/00: 1) Removed OB*w fix (see above).
                0032 C             2) Added variable CMAX. Maximum diagnosed phase speed is now
a24437ad80 Jean*0033 C                clamped to CMAX. For stability of AB-II scheme (CFL) the
42c525bfb4 Alis*0034 C                (non-dimensional) phase speed must be <0.5
                0035 C             3) (Sonya Legg) Changed application of uVel and vVel.
a24437ad80 Jean*0036 C                uVel on the western OB is actually applied at I_obc+1
42c525bfb4 Alis*0037 C                while vVel on the southern OB is applied at J_obc+1.
e9b27c9813 Alis*0038 C             4) (Sonya Legg) Added templates for forced OBs.
42c525bfb4 Alis*0039 C
                0040 C SPK 7/17/00: Non-uniform resolution is now taken into account in diagnosing
                0041 C              phase speeds and time-stepping OB values. CL is still the
                0042 C              non-dimensional phase speed; CVEL is the dimensional phase
a24437ad80 Jean*0043 C              speed: CVEL = CL*(dx or dy)/dt, where dx and dy is the
                0044 C              appropriate grid spacings. Note that CMAX (with which CL
42c525bfb4 Alis*0045 C              is compared) remains non-dimensional.
a24437ad80 Jean*0046 C
                0047 C SPK 7/18/00: Added code to allow filtering of phase speed following
                0048 C              Blumberg and Kantha. There is now a separate array
42c525bfb4 Alis*0049 C              CVEL_**, where **=Variable(U,V,T,S,W)Boundary(E,W,N,S) for
a24437ad80 Jean*0050 C              the dimensional phase speed. These arrays are initialized to
42c525bfb4 Alis*0051 C              zero in ini_obcs.F. CVEL_** is filtered according to
                0052 C              CVEL_** = fracCVEL*CVEL(new) + (1-fracCVEL)*CVEL_**(old).
                0053 C              fracCVEL=1.0 turns off filtering.
                0054 C
a24437ad80 Jean*0055 C SPK 7/26/00: Changed code to average phase speed. A new variable
42c525bfb4 Alis*0056 C              'cvelTimeScale' was created. This variable must now be
a24437ad80 Jean*0057 C              specified. Then, fracCVEL=deltaT/cvelTimeScale.
                0058 C              Since the goal is to smooth out the 'singularities' in the
                0059 C              diagnosed phase speed, cvelTimeScale could be picked as the
                0060 C              duration of the singular period in the unfiltered case. Thus,
                0061 C              for a plane wave cvelTimeScale might be the time take for the
42c525bfb4 Alis*0062 C              wave to travel a distance DX, where DX is the width of the region
                0063 C              near which d(phi)/dx is small.
f2f2d1c632 Alis*0064 C
                0065 C JBG 4/10/03: Fixed phase speed at western boundary (as suggested by
3daafce20b Jean*0066 C              Dale Durran in his MWR paper). Fixed value (in m/s) is
f2f2d1c632 Alis*0067 C              passed in as variable CFIX in data.obcs.
                0068 C              also now allow choice of Orlanski or fixed wavespeed
                0069 C              (by means of new booleans useFixedCEast and
                0070 C              useFixedCWest) without having to recompile each time
                0071 C
42c525bfb4 Alis*0072 
                0073 C     == Routine arguments ==
                0074       INTEGER bi, bj
                0075       _RL futureTime
                0076       _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0077       _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0078       _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0079       _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0080       _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0081       INTEGER myThid
                0082 
                0083 #ifdef ALLOW_ORLANSKI
96bbd4e2a5 Patr*0084 #ifdef ALLOW_OBCS_EAST
42c525bfb4 Alis*0085 C     == Local variables ==
                0086       INTEGER J, K, I_obc
                0087       _RL CL, ab1, ab2, fracCVEL, f1, f2
9e0fdb4065 Garr*0088       _RL denom
                0089 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
42c525bfb4 Alis*0090 
                0091       ab1   =  1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
                0092       ab2   = -0.5 _d 0 - abEps
                0093       /* CMAX is maximum allowable phase speed-CFL for AB-II */
                0094       /* cvelTimeScale is averaging period for phase speed in sec. */
                0095 
                0096       fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
e9b27c9813 Alis*0097       f1 = fracCVEL /* dont change this. Set cvelTimeScale */
                0098       f2 = 1.0-fracCVEL   /* dont change this. set cvelTimeScale */
42c525bfb4 Alis*0099 
                0100 C     Eastern OB (Orlanski Radiation Condition)
                0101       DO K=1,Nr
74019f026d Jean*0102          DO J=1-OLy,sNy+OLy
42c525bfb4 Alis*0103             I_obc=OB_Ie(J,bi,bj)
74019f026d Jean*0104             IF ( I_obc.NE.OB_indexNone ) THEN
9e0fdb4065 Garr*0105 C-             uVel
                0106                denom =
42c525bfb4 Alis*0107      &          (ab1*UE_STORE_2(J,K,bi,bj) + ab2*UE_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0108                IF ( denom .NE. 0 _d 0 ) THEN
                0109                  CL =-(uVel(I_obc-1,J,K,bi,bj)-UE_STORE_1(J,K,bi,bj))
                0110      &              / denom
                0111                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0112                ELSE
                0113                  CL = 0. _d 0
42c525bfb4 Alis*0114                ENDIF
f2f2d1c632 Alis*0115                IF (useFixedCEast) THEN
41a255859f Jean*0116 C                Fixed phase speed (ignoring all of that painstakingly
                0117 C                saved data...)
f2f2d1c632 Alis*0118                  CVEL_UE(J,K,bi,bj) = CFIX
                0119                ELSE
725711d378 Alis*0120                  CVEL_UE(J,K,bi,bj) = f1*(CL*dxF(I_obc-2,J,bi,bj)/deltaT
                0121      &              )+f2*CVEL_UE(J,K,bi,bj)
f2f2d1c632 Alis*0122                ENDIF
42c525bfb4 Alis*0123 C              update OBC to next timestep
                0124                OBEu(J,K,bi,bj)=uVel(I_obc,J,K,bi,bj)-
077d6102ae Alis*0125      &           CVEL_UE(J,K,bi,bj)*(deltaT*recip_dxF(I_obc-1,J,bi,bj))*
42c525bfb4 Alis*0126      &           (ab1*(uVel(I_obc,J,K,bi,bj)-uVel(I_obc-1,J,K,bi,bj)) +
                0127      &           ab2*(UE_STORE_4(J,K,bi,bj)-UE_STORE_1(J,K,bi,bj)))
9e0fdb4065 Garr*0128 
                0129 C-             vVel
                0130                denom =
42c525bfb4 Alis*0131      &          (ab1*VE_STORE_2(J,K,bi,bj) + ab2*VE_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0132                IF ( denom .NE. 0 _d 0 ) THEN
                0133                  CL =-(vVel(I_obc-1,J,K,bi,bj)-VE_STORE_1(J,K,bi,bj))
                0134      &              / denom
                0135                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0136                ELSE
                0137                  CL = 0. _d 0
42c525bfb4 Alis*0138                ENDIF
f2f2d1c632 Alis*0139                IF (useFixedCEast) THEN
41a255859f Jean*0140 C                Fixed phase speed (ignoring all of that painstakingly
                0141 C                saved data...)
f2f2d1c632 Alis*0142                  CVEL_VE(J,K,bi,bj) = CFIX
                0143                ELSE
                0144                  CVEL_VE(J,K,bi,bj) = f1*(CL*dxV(I_obc-1,J,bi,bj)
                0145      $                 /deltaT)+f2*CVEL_VE(J,K,bi,bj)
                0146                ENDIF
42c525bfb4 Alis*0147 C              update OBC to next timestep
                0148                OBEv(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)-
077d6102ae Alis*0149      &           CVEL_VE(J,K,bi,bj)*(deltaT*recip_dxV(I_obc,J,bi,bj))*
a24437ad80 Jean*0150      &           (ab1*(vVel(I_obc,J,K,bi,bj)-vVel(I_obc-1,J,K,bi,bj)) +
077d6102ae Alis*0151      &           ab2*(VE_STORE_4(J,K,bi,bj)-VE_STORE_1(J,K,bi,bj)))
9e0fdb4065 Garr*0152 
                0153 C-             Temperature
                0154                denom =
42c525bfb4 Alis*0155      &          (ab1*TE_STORE_2(J,K,bi,bj) + ab2*TE_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0156                IF ( denom .NE. 0 _d 0 ) THEN
                0157                  CL =-(theta(I_obc-1,J,K,bi,bj)-TE_STORE_1(J,K,bi,bj))
                0158      &              / denom
                0159                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0160                ELSE
                0161                  CL = 0. _d 0
42c525bfb4 Alis*0162                ENDIF
f2f2d1c632 Alis*0163                IF (useFixedCEast) THEN
41a255859f Jean*0164 C                Fixed phase speed (ignoring all of that painstakingly
                0165 C                saved data...)
f2f2d1c632 Alis*0166                  CVEL_TE(J,K,bi,bj) = CFIX
                0167                ELSE
                0168                  CVEL_TE(J,K,bi,bj) = f1*(CL*dxC(I_obc-1,J,bi,bj)
                0169      $                 /deltaT)+f2*CVEL_TE(J,K,bi,bj)
                0170                ENDIF
42c525bfb4 Alis*0171 C              update OBC to next timestep
                0172                OBEt(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)-
077d6102ae Alis*0173      &           CVEL_TE(J,K,bi,bj)*(deltaT*recip_dxC(I_obc,J,bi,bj))*
                0174      &           (ab1*(theta(I_obc,J,K,bi,bj)-theta(I_obc-1,J,K,bi,bj))+
42c525bfb4 Alis*0175      &           ab2*(TE_STORE_4(J,K,bi,bj)-TE_STORE_1(J,K,bi,bj)))
9e0fdb4065 Garr*0176 
                0177 C-             Salinity
                0178                denom =
42c525bfb4 Alis*0179      &          (ab1*SE_STORE_2(J,K,bi,bj) + ab2*SE_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0180                IF ( denom .NE. 0 _d 0 ) THEN
                0181                  CL =-(salt(I_obc-1,J,K,bi,bj)-SE_STORE_1(J,K,bi,bj))
                0182      &              / denom
                0183                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0184                ELSE
                0185                  CL = 0. _d 0
42c525bfb4 Alis*0186                ENDIF
f2f2d1c632 Alis*0187                IF (useFixedCEast) THEN
41a255859f Jean*0188 C                Fixed phase speed (ignoring all of that painstakingly
                0189 C                saved data...)
f2f2d1c632 Alis*0190                  CVEL_SE(J,K,bi,bj) = CFIX
                0191                ELSE
                0192                  CVEL_SE(J,K,bi,bj) = f1*(CL*dxC(I_obc-1,J,bi,bj)
                0193      $                 /deltaT)+f2*CVEL_SE(J,K,bi,bj)
                0194                ENDIF
42c525bfb4 Alis*0195 C              update OBC to next timestep
                0196                OBEs(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)-
077d6102ae Alis*0197      &           CVEL_SE(J,K,bi,bj)*(deltaT*recip_dxC(I_obc,J,bi,bj))*
42c525bfb4 Alis*0198      &           (ab1*(salt(I_obc,J,K,bi,bj)-salt(I_obc-1,J,K,bi,bj))+
a24437ad80 Jean*0199      &           ab2*(SE_STORE_4(J,K,bi,bj)-SE_STORE_1(J,K,bi,bj)))
9e0fdb4065 Garr*0200 
42c525bfb4 Alis*0201 #ifdef ALLOW_NONHYDROSTATIC
41a255859f Jean*0202              IF ( nonHydrostatic ) THEN
9e0fdb4065 Garr*0203 C-             wVel
                0204                denom =
41a255859f Jean*0205      &          (ab1*WE_STORE_2(J,K,bi,bj)+ab2*WE_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0206                IF ( denom .NE. 0 _d 0 ) THEN
                0207                  CL =-(wVel(I_obc-1,J,K,bi,bj)-WE_STORE_1(J,K,bi,bj))
                0208      &              / denom
                0209                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0210                ELSE
                0211                  CL = 0. _d 0
41a255859f Jean*0212                ENDIF
f2f2d1c632 Alis*0213                IF (useFixedCEast) THEN
41a255859f Jean*0214 C                Fixed phase speed (ignoring all of that painstakingly
                0215 C                saved data...)
f2f2d1c632 Alis*0216                  CVEL_WE(J,K,bi,bj) = CFIX
                0217                ELSE
                0218                  CVEL_WE(J,K,bi,bj)=f1*(CL*dxC(I_obc-1,J,bi,bj)/deltaT)
42c525bfb4 Alis*0219      &                   + f2*CVEL_WE(J,K,bi,bj)
f2f2d1c632 Alis*0220                ENDIF
41a255859f Jean*0221 C              update OBC to next timestep
                0222                OBEw(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)-
                0223      &           CVEL_WE(J,K,bi,bj)*(deltaT*recip_dxC(I_obc,J,bi,bj))*
                0224      &           (ab1*(wVel(I_obc,J,K,bi,bj)-wVel(I_obc-1,J,K,bi,bj))+
                0225      &           ab2*(WE_STORE_4(J,K,bi,bj)-WE_STORE_1(J,K,bi,bj)))
                0226              ENDIF
                0227 #endif /* ALLOW_NONHYDROSTATIC */
9e0fdb4065 Garr*0228 
                0229 C-             update/save storage arrays
42c525bfb4 Alis*0230 C              uVel
                0231 C              copy t-1 to t-2 array
                0232                UE_STORE_3(J,K,bi,bj)=UE_STORE_2(J,K,bi,bj)
                0233 C              copy (current time) t to t-1 arrays
                0234                UE_STORE_2(J,K,bi,bj)=uVel(I_obc-1,J,K,bi,bj) -
                0235      &         uVel(I_obc-2,J,K,bi,bj)
                0236                UE_STORE_1(J,K,bi,bj)=uVel(I_obc-1,J,K,bi,bj)
                0237                UE_STORE_4(J,K,bi,bj)=uVel(I_obc,J,K,bi,bj)
                0238 C              vVel
                0239 C              copy t-1 to t-2 array
                0240                VE_STORE_3(J,K,bi,bj)=VE_STORE_2(J,K,bi,bj)
                0241 C              copy (current time) t to t-1 arrays
                0242                VE_STORE_2(J,K,bi,bj)=vVel(I_obc-1,J,K,bi,bj) -
                0243      &         vVel(I_obc-2,J,K,bi,bj)
                0244                VE_STORE_1(J,K,bi,bj)=vVel(I_obc-1,J,K,bi,bj)
                0245                VE_STORE_4(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)
                0246 C              Temperature
                0247 C              copy t-1 to t-2 array
                0248                TE_STORE_3(J,K,bi,bj)=TE_STORE_2(J,K,bi,bj)
                0249 C              copy (current time) t to t-1 arrays
                0250                TE_STORE_2(J,K,bi,bj)=theta(I_obc-1,J,K,bi,bj) -
                0251      &         theta(I_obc-2,J,K,bi,bj)
                0252                TE_STORE_1(J,K,bi,bj)=theta(I_obc-1,J,K,bi,bj)
                0253                TE_STORE_4(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)
                0254 C              Salinity
                0255 C              copy t-1 to t-2 array
                0256                SE_STORE_3(J,K,bi,bj)=SE_STORE_2(J,K,bi,bj)
                0257 C              copy (current time) t to t-1 arrays
                0258                SE_STORE_2(J,K,bi,bj)=salt(I_obc-1,J,K,bi,bj) -
                0259      &         salt(I_obc-2,J,K,bi,bj)
                0260                SE_STORE_1(J,K,bi,bj)=salt(I_obc-1,J,K,bi,bj)
                0261                SE_STORE_4(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)
                0262 #ifdef ALLOW_NONHYDROSTATIC
41a255859f Jean*0263              IF ( nonHydrostatic ) THEN
                0264 C              wVel
42c525bfb4 Alis*0265 C              copy t-1 to t-2 array
                0266                WE_STORE_3(J,K,bi,bj)=WE_STORE_2(J,K,bi,bj)
                0267 C              copy (current time) t to t-1 arrays
                0268                WE_STORE_2(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj) -
                0269      &         wVel(I_obc-2,J,K,bi,bj)
                0270                WE_STORE_1(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
                0271                WE_STORE_4(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)
41a255859f Jean*0272              ENDIF
                0273 #endif /* ALLOW_NONHYDROSTATIC */
42c525bfb4 Alis*0274             ENDIF
                0275          ENDDO
                0276       ENDDO
                0277 
9e0fdb4065 Garr*0278 #endif /* ALLOW_OBCS_EAST */
42c525bfb4 Alis*0279 #endif /* ALLOW_ORLANSKI */
                0280       RETURN
                0281       END