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_SOUTH( 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_SOUTH                                |
                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_SOUTH
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     Southern OB (Orlanski Radiation Condition)
                0093        DO K=1,Nr
74019f026d Jean*0094          DO I=1-OLx,sNx+OLx
42c525bfb4 Alis*0095             J_obc=OB_Js(I,bi,bj)
74019f026d Jean*0096             IF ( J_obc.NE.OB_indexNone ) THEN
9e0fdb4065 Garr*0097 C-             uVel
                0098                 denom =
42c525bfb4 Alis*0099      &          (ab1*US_STORE_2(I,K,bi,bj) + ab2*US_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)-US_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_US(I,K,bi,bj) = f1*(CL*dyU(I,J_obc+2,bi,bj)/deltaT)+
                0108      &                f2*CVEL_US(I,K,bi,bj)
                0109 C              update OBC to next timestep
                0110                OBSu(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)+
9f591cc1fa Jean*0111      &          CVEL_US(I,K,bi,bj)*deltaT*recip_dyU(I,J_obc+1,bi,bj)*
42c525bfb4 Alis*0112      &          (ab1*(uVel(I,J_obc+1,K,bi,bj)-uVel(I,J_obc,K,bi,bj)) +
                0113      &         ab2*(US_STORE_1(I,K,bi,bj)-US_STORE_4(I,K,bi,bj)))
9e0fdb4065 Garr*0114 
                0115 C-             vVel (to be applied at J_obc+1)
                0116                 denom =
42c525bfb4 Alis*0117      &          (ab1*VS_STORE_2(I,K,bi,bj) + ab2*VS_STORE_3(I,K,bi,bj))
9e0fdb4065 Garr*0118                IF ( denom .NE. 0 _d 0 ) THEN
                0119                  CL = (vVel(I,J_obc+2,K,bi,bj)-VS_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_VS(I,K,bi,bj) = f1*(CL*dyF(I,J_obc+2,bi,bj)/deltaT)+
                0126      &                f2*CVEL_VS(I,K,bi,bj)
                0127 C              update OBC to next timestep
                0128                OBSv(I,K,bi,bj)=vVel(I,J_obc+1,K,bi,bj)+
9f591cc1fa Jean*0129      &          CVEL_VS(I,K,bi,bj)*deltaT*recip_dyF(I,J_obc+1,bi,bj)*
42c525bfb4 Alis*0130      &          (ab1*(vVel(I,J_obc+2,K,bi,bj)-vVel(I,J_obc+1,K,bi,bj))+
                0131      &          ab2*(VS_STORE_1(I,K,bi,bj)-VS_STORE_4(I,K,bi,bj)))
9e0fdb4065 Garr*0132 
                0133 C-             Temperature
                0134                 denom =
42c525bfb4 Alis*0135      &          (ab1*TS_STORE_2(I,K,bi,bj) + ab2*TS_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)-TS_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_TS(I,K,bi,bj) = f1*(CL*dyC(I,J_obc+2,bi,bj)/deltaT)+
                0144      &                f2*CVEL_TS(I,K,bi,bj)
                0145 C              update OBC to next timestep
                0146                OBSt(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)+
9f591cc1fa Jean*0147      &          CVEL_TS(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc+1,bi,bj)*
42c525bfb4 Alis*0148      &          (ab1*(theta(I,J_obc+1,K,bi,bj)-theta(I,J_obc,K,bi,bj))+
                0149      &          ab2*(TS_STORE_1(I,K,bi,bj)-TS_STORE_4(I,K,bi,bj)))
9e0fdb4065 Garr*0150 
                0151 C-             Salinity
                0152                 denom =
42c525bfb4 Alis*0153      &          (ab1*SS_STORE_2(I,K,bi,bj) + ab2*SS_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)-SS_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_SS(I,K,bi,bj) = f1*(CL*dyC(I,J_obc+2,bi,bj)/deltaT)+
                0162      &                f2*CVEL_SS(I,K,bi,bj)
                0163 C              update OBC to next timestep
                0164                OBSs(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)+
9f591cc1fa Jean*0165      &          CVEL_SS(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc+1,bi,bj)*
42c525bfb4 Alis*0166      &          (ab1*(salt(I,J_obc+1,K,bi,bj)-salt(I,J_obc,K,bi,bj)) +
                0167      &         ab2*(SS_STORE_1(I,K,bi,bj)-SS_STORE_4(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*WS_STORE_2(I,K,bi,bj)+ab2*WS_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)-WS_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_WS(I,K,bi,bj)=f1*(CL*dyC(I,J_obc+2,bi,bj)/deltaT)
42c525bfb4 Alis*0182      &                   + f2*CVEL_WS(I,K,bi,bj)
41a255859f Jean*0183 C              update OBC to next timestep
                0184                OBSw(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)+
                0185      &           CVEL_WS(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc+1,bi,bj)*
                0186      &           (ab1*(wVel(I,J_obc+1,K,bi,bj)-wVel(I,J_obc,K,bi,bj))+
                0187      &           ab2*(WS_STORE_1(I,K,bi,bj)-WS_STORE_4(I,K,bi,bj)))
                0188              ENDIF
                0189 #endif /* ALLOW_NONHYDROSTATIC */
9e0fdb4065 Garr*0190 
42c525bfb4 Alis*0191 C              update/save storage arrays
                0192 C              uVel
                0193 C              copy t-1 to t-2 array
                0194                US_STORE_3(I,K,bi,bj)=US_STORE_2(I,K,bi,bj)
                0195 C              copy (current time) t to t-1 arrays
                0196                US_STORE_2(I,K,bi,bj)=uVel(I,J_obc+2,K,bi,bj) -
                0197      &         uVel(I,J_obc+1,K,bi,bj)
                0198                US_STORE_1(I,K,bi,bj)=uVel(I,J_obc+1,K,bi,bj)
                0199                US_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                VS_STORE_3(I,K,bi,bj)=VS_STORE_2(I,K,bi,bj)
                0203 C              copy (current time) t to t-1 arrays
                0204                VS_STORE_2(I,K,bi,bj)=vVel(I,J_obc+3,K,bi,bj) -
                0205      &         vVel(I,J_obc+2,K,bi,bj)
                0206                VS_STORE_1(I,K,bi,bj)=vVel(I,J_obc+2,K,bi,bj)
                0207                VS_STORE_4(I,K,bi,bj)=vVel(I,J_obc+1,K,bi,bj)
                0208 C              Temperature
                0209 C              copy t-1 to t-2 array
                0210                TS_STORE_3(I,K,bi,bj)=TS_STORE_2(I,K,bi,bj)
                0211 C              copy (current time) t to t-1 arrays
                0212                TS_STORE_2(I,K,bi,bj)=theta(I,J_obc+2,K,bi,bj) -
                0213      &         theta(I,J_obc+1,K,bi,bj)
                0214                TS_STORE_1(I,K,bi,bj)=theta(I,J_obc+1,K,bi,bj)
                0215                TS_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                SS_STORE_3(I,K,bi,bj)=SS_STORE_2(I,K,bi,bj)
                0219 C              copy (current time) t to t-1 arrays
                0220                SS_STORE_2(I,K,bi,bj)=salt(I,J_obc+2,K,bi,bj) -
                0221      &         salt(I,J_obc+1,K,bi,bj)
                0222                SS_STORE_1(I,K,bi,bj)=salt(I,J_obc+1,K,bi,bj)
                0223                SS_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                WS_STORE_3(I,K,bi,bj)=WS_STORE_2(I,K,bi,bj)
                0229 C              copy (current time) t to t-1 arrays
                0230                WS_STORE_2(I,K,bi,bj)=wVel(I,J_obc+2,K,bi,bj) -
                0231      &         wVel(I,J_obc+1,K,bi,bj)
                0232                WS_STORE_1(I,K,bi,bj)=wVel(I,J_obc+1,K,bi,bj)
                0233                WS_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 
e9741b789e Jean*0240 #endif /* ALLOW_OBCS_SOUTH */
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!