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_NORTH( bi, bj, futureTime,
                0004      I                      uVel, vVel, wVel, theta, salt,
42c525bfb4 Alis*0005      I                      myThid )
9e0fdb4065 Garr*0006 C     *==========================================================*
42c525bfb4 Alis*0007 C     | SUBROUTINE OBCS_RADIATE                                  |
                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.
                0064 
                0065 C     == Routine arguments ==
                0066       INTEGER bi, bj
                0067       _RL futureTime
                0068       _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0069       _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0070       _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0071       _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0072       _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0073       INTEGER myThid
                0074 
                0075 #ifdef ALLOW_ORLANSKI
96bbd4e2a5 Patr*0076 #ifdef ALLOW_OBCS_NORTH
42c525bfb4 Alis*0077 C     == Local variables ==
                0078       INTEGER  I, K, J_obc
                0079       _RL CL, ab1, ab2, fracCVEL, f1, f2
9e0fdb4065 Garr*0080       _RL denom
                0081 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
42c525bfb4 Alis*0082 
                0083       ab1   =  1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
                0084       ab2   = -0.5 _d 0 - abEps
                0085       /* CMAX is maximum allowable phase speed-CFL for AB-II */
                0086       /* cvelTimeScale is averaging period for phase speed in sec. */
                0087 
                0088       fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
e9b27c9813 Alis*0089       f1 = fracCVEL /* dont change this. Set cvelTimeScale */
                0090       f2 = 1.0-fracCVEL   /* dont change this. set cvelTimeScale */
42c525bfb4 Alis*0091 
                0092 C     Northern OB (Orlanski Radiation Condition)
                0093        DO K=1,Nr
74019f026d Jean*0094          DO I=1-OLx,sNx+OLx
42c525bfb4 Alis*0095             J_obc=OB_Jn(I,bi,bj)
74019f026d Jean*0096             IF ( J_obc.NE.OB_indexNone ) THEN
9e0fdb4065 Garr*0097 C-             uVel
                0098                denom =
42c525bfb4 Alis*0099      &          (ab1*UN_STORE_2(I,K,bi,bj) + ab2*UN_STORE_3(I,K,bi,bj))
9e0fdb4065 Garr*0100                IF ( denom .NE. 0 _d 0 ) THEN
                0101                  CL =-(uVel(I,J_obc-1,K,bi,bj)-UN_STORE_1(I,K,bi,bj))
                0102      &              / denom
                0103                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0104                ELSE
                0105                  CL = 0. _d 0
42c525bfb4 Alis*0106                ENDIF
                0107                CVEL_UN(I,K,bi,bj) = f1*(CL*dyU(I,J_obc-1,bi,bj)/deltaT)+
                0108      &                f2*CVEL_UN(I,K,bi,bj)
                0109 C              update OBC to next timestep
                0110                OBNu(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)-
9f591cc1fa Jean*0111      &           CVEL_UN(I,K,bi,bj)*deltaT*recip_dyU(I,J_obc,bi,bj)*
42c525bfb4 Alis*0112      &           (ab1*(uVel(I,J_obc,K,bi,bj)-uVel(I,J_obc-1,K,bi,bj)) +
                0113      &           ab2*(UN_STORE_4(I,K,bi,bj)-UN_STORE_1(I,K,bi,bj)))
9e0fdb4065 Garr*0114 
                0115 C-             vVel
                0116                denom =
42c525bfb4 Alis*0117      &          (ab1*VN_STORE_2(I,K,bi,bj) + ab2*VN_STORE_3(I,K,bi,bj))
9e0fdb4065 Garr*0118                IF ( denom .NE. 0 _d 0 ) THEN
                0119                  CL =-(vVel(I,J_obc-1,K,bi,bj)-VN_STORE_1(I,K,bi,bj))
                0120      &              / denom
                0121                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0122                ELSE
                0123                  CL = 0. _d 0
42c525bfb4 Alis*0124                ENDIF
                0125                CVEL_VN(I,K,bi,bj) = f1*(CL*dyF(I,J_obc-2,bi,bj)/deltaT)+
                0126      &                f2*CVEL_VN(I,K,bi,bj)
                0127 C              update OBC to next timestep
                0128                OBNv(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)-
9f591cc1fa Jean*0129      &           CVEL_VN(I,K,bi,bj)*deltaT*recip_dyF(I,J_obc-1,bi,bj)*
42c525bfb4 Alis*0130      &           (ab1*(vVel(I,J_obc,K,bi,bj)-vVel(I,J_obc-1,K,bi,bj)) +
a24437ad80 Jean*0131      &           ab2*(VN_STORE_4(I,K,bi,bj)-VN_STORE_1(I,K,bi,bj)))
9e0fdb4065 Garr*0132 
                0133 C-             Temperature
                0134                denom =
42c525bfb4 Alis*0135      &          (ab1*TN_STORE_2(I,K,bi,bj) + ab2*TN_STORE_3(I,K,bi,bj))
9e0fdb4065 Garr*0136                IF ( denom .NE. 0 _d 0 ) THEN
                0137                  CL =-(theta(I,J_obc-1,K,bi,bj)-TN_STORE_1(I,K,bi,bj))
                0138      &              / denom
                0139                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0140                ELSE
                0141                  CL = 0. _d 0
42c525bfb4 Alis*0142                ENDIF
                0143                CVEL_TN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
                0144      &                f2*CVEL_TN(I,K,bi,bj)
                0145 C              update OBC to next timestep
                0146                OBNt(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)-
9f591cc1fa Jean*0147      &          CVEL_TN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
42c525bfb4 Alis*0148      &          (ab1*(theta(I,J_obc,K,bi,bj)-theta(I,J_obc-1,K,bi,bj))+
a24437ad80 Jean*0149      &          ab2*(TN_STORE_4(I,K,bi,bj)-TN_STORE_1(I,K,bi,bj)))
9e0fdb4065 Garr*0150 
                0151 C-             Salinity
                0152                denom =
42c525bfb4 Alis*0153      &          (ab1*SN_STORE_2(I,K,bi,bj) + ab2*SN_STORE_3(I,K,bi,bj))
9e0fdb4065 Garr*0154                IF ( denom .NE. 0 _d 0 ) THEN
                0155                  CL =-(salt(I,J_obc-1,K,bi,bj)-SN_STORE_1(I,K,bi,bj))
                0156      &              / denom
                0157                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0158                ELSE
                0159                  CL = 0. _d 0
42c525bfb4 Alis*0160                ENDIF
                0161                CVEL_SN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
                0162      &                f2*CVEL_SN(I,K,bi,bj)
                0163 C              update OBC to next timestep
                0164                OBNs(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)-
9f591cc1fa Jean*0165      &          CVEL_SN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
a24437ad80 Jean*0166      &          (ab1*(salt(I,J_obc,K,bi,bj)-salt(I,J_obc-1,K,bi,bj)) +
42c525bfb4 Alis*0167      &          ab2*(SN_STORE_4(I,K,bi,bj)-SN_STORE_1(I,K,bi,bj)))
9e0fdb4065 Garr*0168 
42c525bfb4 Alis*0169 #ifdef ALLOW_NONHYDROSTATIC
41a255859f Jean*0170              IF ( nonHydrostatic ) THEN
9e0fdb4065 Garr*0171 C-             wVel
                0172                denom =
41a255859f Jean*0173      &          (ab1*WN_STORE_2(I,K,bi,bj)+ab2*WN_STORE_3(I,K,bi,bj))
9e0fdb4065 Garr*0174                IF ( denom .NE. 0 _d 0 ) THEN
                0175                  CL =-(wVel(I,J_obc-1,K,bi,bj)-WN_STORE_1(I,K,bi,bj))
                0176      &              / denom
                0177                  CL = MIN( MAX( CL, zeroRL ), CMAX )
                0178                ELSE
                0179                  CL = 0. _d 0
41a255859f Jean*0180                ENDIF
                0181                CVEL_WN(I,K,bi,bj)=f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)
                0182      &                + f2*CVEL_WN(I,K,bi,bj)
                0183 C              update OBC to next timestep
                0184                OBNw(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)-
                0185      &           CVEL_WN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
                0186      &           (ab1*(wVel(I,J_obc,K,bi,bj)-wVel(I,J_obc-1,K,bi,bj))+
                0187      &           ab2*(WN_STORE_4(I,K,bi,bj)-WN_STORE_1(I,K,bi,bj)))
                0188              ENDIF
                0189 #endif /* ALLOW_NONHYDROSTATIC */
9e0fdb4065 Garr*0190 
                0191 C-             update/save storage arrays
42c525bfb4 Alis*0192 C              uVel
                0193 C              copy t-1 to t-2 array
                0194                UN_STORE_3(I,K,bi,bj)=UN_STORE_2(I,K,bi,bj)
                0195 C              copy (current time) t to t-1 arrays
                0196                UN_STORE_2(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj) -
                0197      &         uVel(I,J_obc-2,K,bi,bj)
                0198                UN_STORE_1(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj)
                0199                UN_STORE_4(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)
                0200 C              vVel
                0201 C              copy t-1 to t-2 array
                0202                VN_STORE_3(I,K,bi,bj)=VN_STORE_2(I,K,bi,bj)
                0203 C              copy (current time) t to t-1 arrays
                0204                VN_STORE_2(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj) -
                0205      &         vVel(I,J_obc-2,K,bi,bj)
                0206                VN_STORE_1(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj)
                0207                VN_STORE_4(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)
                0208 C              Temperature
                0209 C              copy t-1 to t-2 array
                0210                TN_STORE_3(I,K,bi,bj)=TN_STORE_2(I,K,bi,bj)
                0211 C              copy (current time) t to t-1 arrays
                0212                TN_STORE_2(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj) -
                0213      &         theta(I,J_obc-2,K,bi,bj)
                0214                TN_STORE_1(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj)
                0215                TN_STORE_4(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)
                0216 C              Salinity
                0217 C              copy t-1 to t-2 array
                0218                SN_STORE_3(I,K,bi,bj)=SN_STORE_2(I,K,bi,bj)
                0219 C              copy (current time) t to t-1 arrays
                0220                SN_STORE_2(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj) -
                0221      &         salt(I,J_obc-2,K,bi,bj)
                0222                SN_STORE_1(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj)
                0223                SN_STORE_4(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)
                0224 #ifdef ALLOW_NONHYDROSTATIC
41a255859f Jean*0225              IF ( nonHydrostatic ) THEN
                0226 C              wVel
42c525bfb4 Alis*0227 C              copy t-1 to t-2 array
                0228                WN_STORE_3(I,K,bi,bj)=WN_STORE_2(I,K,bi,bj)
                0229 C              copy (current time) t to t-1 arrays
                0230                WN_STORE_2(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj) -
                0231      &         wVel(I,J_obc-2,K,bi,bj)
                0232                WN_STORE_1(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj)
                0233                WN_STORE_4(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)
41a255859f Jean*0234              ENDIF
                0235 #endif /* ALLOW_NONHYDROSTATIC */
42c525bfb4 Alis*0236             ENDIF
                0237          ENDDO
                0238        ENDDO
                0239 
96bbd4e2a5 Patr*0240 #endif /* ALLOW_OBCS_NORTH */
42c525bfb4 Alis*0241 #endif /* ALLOW_ORLANSKI */
                0242       RETURN
                0243       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!