Back to home page

MITgcm

 
 

    


File indexing completed on 2024-05-14 05:10:41 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_WEST( 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_WEST                                 |
                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 3/24/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
                0069 C JBG 4/10/03: allow choice of Orlanski or fixed wavespeed (by means of new
                0070 C              booleans useFixedCEast and useFixedCWest) without
                0071 C              having to recompile each time
                0072 c SAL 1/7/03:  Fixed bug: implementation for salinity was incomplete
                0073 C
42c525bfb4 Alis*0074 
                0075 C     == Routine arguments ==
                0076       INTEGER bi, bj
                0077       _RL futureTime
                0078       _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0079       _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0080       _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0081       _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0082       _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0083       INTEGER myThid
                0084 
                0085 #ifdef ALLOW_ORLANSKI
96bbd4e2a5 Patr*0086 #ifdef ALLOW_OBCS_WEST
42c525bfb4 Alis*0087 C     == Local variables ==
                0088       INTEGER J, K, I_obc
                0089       _RL CL, ab1, ab2, fracCVEL, f1, f2
9e0fdb4065 Garr*0090       _RL denom
                0091 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
42c525bfb4 Alis*0092 
                0093       ab1   =  1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
                0094       ab2   = -0.5 _d 0 - abEps
                0095       /* CMAX is maximum allowable phase speed-CFL for AB-II */
                0096       /* cvelTimeScale is averaging period for phase speed in sec. */
                0097 
                0098       fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
e9b27c9813 Alis*0099       f1 = fracCVEL /* dont change this. Set cvelTimeScale */
                0100       f2 = 1.0-fracCVEL   /* dont change this. set cvelTimeScale */
42c525bfb4 Alis*0101 
                0102 C      Western OB (Orlanski Radiation Condition)
                0103        DO K=1,Nr
74019f026d Jean*0104          DO J=1-OLy,sNy+OLy
42c525bfb4 Alis*0105             I_obc=OB_Iw(J,bi,bj)
74019f026d Jean*0106             IF ( I_obc.NE.OB_indexNone ) THEN
9e0fdb4065 Garr*0107 C-             uVel (to be applied at I_obc+1)
                0108                 denom =
42c525bfb4 Alis*0109      &          (ab1*UW_STORE_2(J,K,bi,bj) + ab2*UW_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0110                IF ( denom .NE. 0 _d 0 ) THEN
                0111                  CL = (uVel(I_obc+2,J,K,bi,bj)-UW_STORE_1(J,K,bi,bj))
                0112      &              / denom
                0113                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0114                ELSE
                0115                  CL = 0. _d 0
42c525bfb4 Alis*0116                ENDIF
f2f2d1c632 Alis*0117                IF (useFixedCWest) THEN
41a255859f Jean*0118 C                Fixed phase speed (ignoring all of that painstakingly
                0119 C                saved data...)
f2f2d1c632 Alis*0120                  CVEL_UW(J,K,bi,bj) = CFIX
                0121                ELSE
725711d378 Alis*0122                  CVEL_UW(J,K,bi,bj) = f1*(CL*dxF(I_obc+2,J,bi,bj)/deltaT
                0123      &              )+f2*CVEL_UW(J,K,bi,bj)
f2f2d1c632 Alis*0124                ENDIF
42c525bfb4 Alis*0125 C              update OBC to next timestep
                0126                OBWu(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)+
077d6102ae Alis*0127      &          CVEL_UW(J,K,bi,bj)*(deltaT*recip_dxF(I_obc+1,J,bi,bj))*
42c525bfb4 Alis*0128      &          (ab1*(uVel(I_obc+2,J,K,bi,bj)-uVel(I_obc+1,J,K,bi,bj))+
                0129      &          ab2*(UW_STORE_1(J,K,bi,bj)-UW_STORE_4(J,K,bi,bj)))
9e0fdb4065 Garr*0130 
                0131 C-             vVel
                0132                 denom =
42c525bfb4 Alis*0133      &          (ab1*VW_STORE_2(J,K,bi,bj) + ab2*VW_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0134                IF ( denom .NE. 0 _d 0 ) THEN
                0135                  CL = (vVel(I_obc+1,J,K,bi,bj)-VW_STORE_1(J,K,bi,bj))
                0136      &              / denom
                0137                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0138                ELSE
                0139                  CL = 0. _d 0
42c525bfb4 Alis*0140                ENDIF
f2f2d1c632 Alis*0141                IF (useFixedCWest) THEN
41a255859f Jean*0142 C                Fixed phase speed (ignoring all of that painstakingly
                0143 C                saved data...)
f2f2d1c632 Alis*0144                  CVEL_VW(J,K,bi,bj) = CFIX
                0145                ELSE
725711d378 Alis*0146                  CVEL_VW(J,K,bi,bj) = f1*(CL*dxV(I_obc+2,J,bi,bj)/deltaT
                0147      &              )+f2*CVEL_VW(J,K,bi,bj)
f2f2d1c632 Alis*0148                ENDIF
42c525bfb4 Alis*0149 C              update OBC to next timestep
                0150                OBWv(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)+
077d6102ae Alis*0151      &           CVEL_VW(J,K,bi,bj)*(deltaT*recip_dxV(I_obc+1,J,bi,bj))*
42c525bfb4 Alis*0152      &           (ab1*(vVel(I_obc+1,J,K,bi,bj)-vVel(I_obc,J,K,bi,bj))+
                0153      &           ab2*(VW_STORE_1(J,K,bi,bj)-VW_STORE_4(J,K,bi,bj)))
9e0fdb4065 Garr*0154 
                0155 C-             Temperature
                0156                 denom =
42c525bfb4 Alis*0157      &          (ab1*TW_STORE_2(J,K,bi,bj) + ab2*TW_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0158                IF ( denom .NE. 0 _d 0 ) THEN
                0159                  CL = (theta(I_obc+1,J,K,bi,bj)-TW_STORE_1(J,K,bi,bj))
                0160      &              / denom
                0161                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0162                ELSE
                0163                  CL = 0. _d 0
42c525bfb4 Alis*0164                ENDIF
f2f2d1c632 Alis*0165                IF (useFixedCWest) THEN
41a255859f Jean*0166 C                Fixed phase speed (ignoring all of that painstakingly
                0167 C                saved data...)
f2f2d1c632 Alis*0168                  CVEL_TW(J,K,bi,bj) = CFIX
                0169                ELSE
725711d378 Alis*0170                  CVEL_TW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT
                0171      &              )+f2*CVEL_TW(J,K,bi,bj)
f2f2d1c632 Alis*0172                ENDIF
42c525bfb4 Alis*0173 C              update OBC to next timestep
                0174                OBWt(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)+
077d6102ae Alis*0175      &          CVEL_TW(J,K,bi,bj)*(deltaT*recip_dxC(I_obc+1,J,bi,bj))*
42c525bfb4 Alis*0176      &          (ab1*(theta(I_obc+1,J,K,bi,bj)-theta(I_obc,J,K,bi,bj))+
                0177      &          ab2*(TW_STORE_1(J,K,bi,bj)-TW_STORE_4(J,K,bi,bj)))
9e0fdb4065 Garr*0178 
                0179 C-             Salinity
                0180                 denom =
42c525bfb4 Alis*0181      &          (ab1*SW_STORE_2(J,K,bi,bj) + ab2*SW_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0182                IF ( denom .NE. 0 _d 0 ) THEN
                0183                  CL = (salt(I_obc+1,J,K,bi,bj)-SW_STORE_1(J,K,bi,bj))
                0184      &              / denom
                0185                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0186                ELSE
                0187                  CL = 0. _d 0
42c525bfb4 Alis*0188                ENDIF
f2f2d1c632 Alis*0189                IF (useFixedCWest) THEN
41a255859f Jean*0190 C                Fixed phase speed (ignoring all of that painstakingly
                0191 C                saved data...)
f2f2d1c632 Alis*0192                  CVEL_SW(J,K,bi,bj) = CFIX
                0193                ELSE
725711d378 Alis*0194                  CVEL_SW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT
                0195      &              )+f2*CVEL_SW(J,K,bi,bj)
f2f2d1c632 Alis*0196                ENDIF
42c525bfb4 Alis*0197 C              update OBC to next timestep
                0198                OBWs(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)+
077d6102ae Alis*0199      &           CVEL_SW(J,K,bi,bj)*(deltaT*recip_dxC(I_obc+1,J,bi,bj))*
42c525bfb4 Alis*0200      &           (ab1*(salt(I_obc+1,J,K,bi,bj)-salt(I_obc,J,K,bi,bj))+
                0201      &           ab2*(SW_STORE_1(J,K,bi,bj)-SW_STORE_4(J,K,bi,bj)))
9e0fdb4065 Garr*0202 
42c525bfb4 Alis*0203 #ifdef ALLOW_NONHYDROSTATIC
41a255859f Jean*0204              IF ( nonHydrostatic ) THEN
9e0fdb4065 Garr*0205 C-             wVel
                0206                 denom =
41a255859f Jean*0207      &          (ab1*WW_STORE_2(J,K,bi,bj)+ab2*WW_STORE_3(J,K,bi,bj))
9e0fdb4065 Garr*0208                IF ( denom .NE. 0 _d 0 ) THEN
                0209                  CL = (wVel(I_obc+1,J,K,bi,bj)-WW_STORE_1(J,K,bi,bj))
                0210      &              / denom
                0211                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0212                ELSE
                0213                  CL = 0. _d 0
41a255859f Jean*0214                ENDIF
f2f2d1c632 Alis*0215                IF (useFixedCWest) THEN
41a255859f Jean*0216 C                Fixed phase speed (ignoring all of that painstakingly
                0217 C                saved data...)
f2f2d1c632 Alis*0218                  CVEL_WW(J,K,bi,bj) = CFIX
                0219                ELSE
                0220                  CVEL_WW(J,K,bi,bj)=f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)
42c525bfb4 Alis*0221      &                   + f2*CVEL_WW(J,K,bi,bj)
f2f2d1c632 Alis*0222                ENDIF
41a255859f Jean*0223 C              update OBC to next timestep
                0224                OBWw(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)+
077d6102ae Alis*0225      &           CVEL_WW(J,K,bi,bj)*(deltaT*recip_dxC(I_obc+1,J,bi,bj))*
42c525bfb4 Alis*0226      &            (ab1*(wVel(I_obc+1,J,K,bi,bj)-wVel(I_obc,J,K,bi,bj))+
                0227      &            ab2*(WW_STORE_1(J,K,bi,bj)-WW_STORE_4(J,K,bi,bj)))
41a255859f Jean*0228              ENDIF
                0229 #endif /* ALLOW_NONHYDROSTATIC */
9e0fdb4065 Garr*0230 
                0231 C-             update/save storage arrays
42c525bfb4 Alis*0232 C              uVel
                0233 C              copy t-1 to t-2 array
                0234                UW_STORE_3(J,K,bi,bj)=UW_STORE_2(J,K,bi,bj)
                0235 C              copy (current time) t to t-1 arrays
                0236                UW_STORE_2(J,K,bi,bj)=uVel(I_obc+3,J,K,bi,bj) -
                0237      &         uVel(I_obc+2,J,K,bi,bj)
                0238                UW_STORE_1(J,K,bi,bj)=uVel(I_obc+2,J,K,bi,bj)
                0239                UW_STORE_4(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)
                0240 C              vVel
                0241 C              copy t-1 to t-2 array
                0242                VW_STORE_3(J,K,bi,bj)=VW_STORE_2(J,K,bi,bj)
                0243 C              copy (current time) t to t-1 arrays
                0244                VW_STORE_2(J,K,bi,bj)=vVel(I_obc+2,J,K,bi,bj) -
                0245      &         vVel(I_obc+1,J,K,bi,bj)
                0246                VW_STORE_1(J,K,bi,bj)=vVel(I_obc+1,J,K,bi,bj)
                0247                VW_STORE_4(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)
                0248 C              Temperature
                0249 C              copy t-1 to t-2 array
                0250                TW_STORE_3(J,K,bi,bj)=TW_STORE_2(J,K,bi,bj)
                0251 C              copy (current time) t to t-1 arrays
                0252                TW_STORE_2(J,K,bi,bj)=theta(I_obc+2,J,K,bi,bj) -
                0253      &         theta(I_obc+1,J,K,bi,bj)
                0254                TW_STORE_1(J,K,bi,bj)=theta(I_obc+1,J,K,bi,bj)
                0255                TW_STORE_4(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)
f2f2d1c632 Alis*0256 c              Salinity
                0257 C              copy t-1 to t-2 array
41a255859f Jean*0258                SW_STORE_3(J,K,bi,bj)=SW_STORE_2(J,K,bi,bj)
f2f2d1c632 Alis*0259 C              copy (current time) t to t-1 arrays
                0260                SW_STORE_2(J,K,bi,bj)=salt(I_obc+2,J,K,bi,bj) -
                0261      &         salt(I_obc+1,J,K,bi,bj)
                0262                SW_STORE_1(J,K,bi,bj)=salt(I_obc+1,J,K,bi,bj)
41a255859f Jean*0263                SW_STORE_4(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)
42c525bfb4 Alis*0264 #ifdef ALLOW_NONHYDROSTATIC
41a255859f Jean*0265              IF ( nonHydrostatic ) THEN
                0266 C              wVel
42c525bfb4 Alis*0267 C              copy t-1 to t-2 array
                0268                WW_STORE_3(J,K,bi,bj)=WW_STORE_2(J,K,bi,bj)
                0269 C              copy (current time) t to t-1 arrays
                0270                WW_STORE_2(J,K,bi,bj)=wVel(I_obc+2,J,K,bi,bj) -
                0271      &         wVel(I_obc+1,J,K,bi,bj)
                0272                WW_STORE_1(J,K,bi,bj)=wVel(I_obc+1,J,K,bi,bj)
                0273                WW_STORE_4(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)
41a255859f Jean*0274              ENDIF
                0275 #endif /* ALLOW_NONHYDROSTATIC */
42c525bfb4 Alis*0276             ENDIF
                0277          ENDDO
                0278       ENDDO
                0279 
9e0fdb4065 Garr*0280 #endif /* ALLOW_OBCS_WEST */
42c525bfb4 Alis*0281 #endif /* ALLOW_ORLANSKI */
                0282       RETURN
                0283       END