/[MITgcm]/MITgcm/pkg/obcs/orlanski_west.F
ViewVC logotype

Diff of /MITgcm/pkg/obcs/orlanski_west.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by adcroft, Tue Jan 30 21:03:01 2001 UTC revision 1.2 by adcroft, Fri Feb 2 21:36:30 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    C $Name$
3    
4    #include "OBCS_OPTIONS.h"
5    
6          SUBROUTINE ORLANSKI_WEST( bi, bj, futureTime,
7         I                      uVel, vVel, wVel, theta, salt,
8         I                      myThid )
9    C     /==========================================================\
10    C     | SUBROUTINE ORLANSKI_WEST                                 |
11    C     | o Calculate future boundary data at open boundaries      |
12    C     |   at time = futureTime by applying Orlanski radiation    |
13    C     |   conditions.                                            |
14    C     |==========================================================|
15    C     |                                                          |
16    C     \==========================================================/
17          IMPLICIT NONE
18    
19    C     === Global variables ===
20    #include "SIZE.h"
21    #include "EEPARAMS.h"
22    #include "PARAMS.h"
23    #include "GRID.h"
24    #include "OBCS.h"
25    #include "ORLANSKI.h"
26    
27    C SPK 6/2/00: Added radiative OBC's for salinity.
28    C SPK 6/6/00: Changed calculation of OB*w. When K=1, the
29    C             upstream value is used. For example on the eastern OB:
30    C                IF (K.EQ.1) THEN
31    C                   OBEw(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
32    C                ENDIF
33    C                        
34    C SPK 7/7/00: 1) Removed OB*w fix (see above).
35    C             2) Added variable CMAX. Maximum diagnosed phase speed is now
36    C                clamped to CMAX. For stability of AB-II scheme (CFL) the
37    C                (non-dimensional) phase speed must be <0.5
38    C             3) (Sonya Legg) Changed application of uVel and vVel.
39    C                uVel on the western OB is actually applied at I_obc+1
40    C                while vVel on the southern OB is applied at J_obc+1.
41    C             4) (Sonya Legg) Added templates for forced OB's.
42    C
43    C SPK 7/17/00: Non-uniform resolution is now taken into account in diagnosing
44    C              phase speeds and time-stepping OB values. CL is still the
45    C              non-dimensional phase speed; CVEL is the dimensional phase
46    C              speed: CVEL = CL*(dx or dy)/dt, where dx and dy is the
47    C              appropriate grid spacings. Note that CMAX (with which CL
48    C              is compared) remains non-dimensional.
49    C
50    C SPK 7/18/00: Added code to allow filtering of phase speed following
51    C              Blumberg and Kantha. There is now a separate array
52    C              CVEL_**, where **=Variable(U,V,T,S,W)Boundary(E,W,N,S) for
53    C              the dimensional phase speed. These arrays are initialized to
54    C              zero in ini_obcs.F. CVEL_** is filtered according to
55    C              CVEL_** = fracCVEL*CVEL(new) + (1-fracCVEL)*CVEL_**(old).
56    C              fracCVEL=1.0 turns off filtering.
57    C
58    C SPK 7/26/00: Changed code to average phase speed. A new variable
59    C              'cvelTimeScale' was created. This variable must now be
60    C              specified. Then, fracCVEL=deltaT/cvelTimeScale.
61    C              Since the goal is to smooth out the 'singularities' in the
62    C              diagnosed phase speed, cvelTimeScale could be picked as the
63    C              duration of the singular period in the unfiltered case. Thus,
64    C              for a plane wave cvelTimeScale might be the time take for the
65    C              wave to travel a distance DX, where DX is the width of the region
66    C              near which d(phi)/dx is small.
67    
68    C     == Routine arguments ==
69          INTEGER bi, bj
70          _RL futureTime
71          _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
72          _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
73          _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
74          _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
75          _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
76          INTEGER myThid
77    
78    #ifdef ALLOW_ORLANSKI
79    
80    C     == Local variables ==
81          INTEGER J, K, I_obc
82          _RL CL, ab1, ab2, fracCVEL, f1, f2
83    
84          ab1   =  1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
85          ab2   = -0.5 _d 0 - abEps
86          /* CMAX is maximum allowable phase speed-CFL for AB-II */
87          /* cvelTimeScale is averaging period for phase speed in sec. */
88    
89          fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
90          f1 = fracCVEL /* don't change this. Set cvelTimeScale */
91          f2 = 1.0-fracCVEL   /* don't change this. set cvelTimeScale */
92    
93    C      Western OB (Orlanski Radiation Condition)
94           DO K=1,Nr
95             DO J=1-Oly,sNy+Oly
96                I_obc=OB_Iw(J,bi,bj)
97                IF (I_obc.ne.0) THEN
98    C              uVel (to be applied at I_obc+1)
99                   IF ((UW_STORE_2(J,K,bi,bj).eq.0.).and.
100         &            (UW_STORE_3(J,K,bi,bj).eq.0.)) THEN
101                      CL=0.
102                   ELSE
103                      CL=(uVel(I_obc+2,J,K,bi,bj)-UW_STORE_1(J,K,bi,bj))/
104         &          (ab1*UW_STORE_2(J,K,bi,bj) + ab2*UW_STORE_3(J,K,bi,bj))
105                   ENDIF
106                   IF (CL.lt.0.) THEN
107                      CL=0.
108                   ELSEIF (CL.gt.CMAX) THEN
109                      CL=CMAX
110                   ENDIF
111                   CVEL_UW(J,K,bi,bj) = f1*(CL*dxF(I_obc+2,J,bi,bj)/deltaT)+
112         &                f2*CVEL_UW(J,K,bi,bj)
113    C              update OBC to next timestep
114                   OBWu(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)+
115         &          CVEL_UW(J,K,bi,bj)*(deltaT/dxF(I_obc+1,J,bi,bj))*
116         &          (ab1*(uVel(I_obc+2,J,K,bi,bj)-uVel(I_obc+1,J,K,bi,bj))+
117         &          ab2*(UW_STORE_1(J,K,bi,bj)-UW_STORE_4(J,K,bi,bj)))
118    C              vVel
119                   IF ((VW_STORE_2(J,K,bi,bj).eq.0.).and.
120         &            (VW_STORE_3(J,K,bi,bj).eq.0.)) THEN
121                      CL=0.
122                   ELSE
123                      CL=(vVel(I_obc+1,J,K,bi,bj)-VW_STORE_1(J,K,bi,bj))/
124         &          (ab1*VW_STORE_2(J,K,bi,bj) + ab2*VW_STORE_3(J,K,bi,bj))
125                   ENDIF
126                   IF (CL.lt.0.) THEN
127                      CL=0.
128                   ELSEIF (CL.gt.CMAX) THEN
129                      CL=CMAX
130                   ENDIF
131                   CVEL_VW(J,K,bi,bj) = f1*(CL*dxV(I_obc+2,J,bi,bj)/deltaT)+
132         &                f2*CVEL_VW(J,K,bi,bj)
133    C              update OBC to next timestep
134                   OBWv(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)+
135         &           CVEL_VW(J,K,bi,bj)*(deltaT/dxV(I_obc+1,J,bi,bj))*
136         &           (ab1*(vVel(I_obc+1,J,K,bi,bj)-vVel(I_obc,J,K,bi,bj))+
137         &           ab2*(VW_STORE_1(J,K,bi,bj)-VW_STORE_4(J,K,bi,bj)))
138    C              Temperature
139                   IF ((TW_STORE_2(J,K,bi,bj).eq.0.).and.
140         &            (TW_STORE_3(J,K,bi,bj).eq.0.)) THEN
141                      CL=0.
142                   ELSE
143                      CL=(theta(I_obc+1,J,K,bi,bj)-TW_STORE_1(J,K,bi,bj))/
144         &          (ab1*TW_STORE_2(J,K,bi,bj) + ab2*TW_STORE_3(J,K,bi,bj))
145                   ENDIF
146                   IF (CL.lt.0.) THEN
147                      CL=0.
148                   ELSEIF (CL.gt.CMAX) THEN
149                      CL=CMAX
150                   ENDIF
151                   CVEL_TW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)+
152         &                f2*CVEL_TW(J,K,bi,bj)
153    C              update OBC to next timestep
154                   OBWt(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)+
155         &          CVEL_TW(J,K,bi,bj)*(deltaT/dxC(I_obc+1,J,bi,bj))*
156         &          (ab1*(theta(I_obc+1,J,K,bi,bj)-theta(I_obc,J,K,bi,bj))+
157         &          ab2*(TW_STORE_1(J,K,bi,bj)-TW_STORE_4(J,K,bi,bj)))
158    C              Salinity
159                   IF ((SW_STORE_2(J,K,bi,bj).eq.0.).and.
160         &            (SW_STORE_3(J,K,bi,bj).eq.0.)) THEN
161                      CL=0.
162                   ELSE
163                      CL=(salt(I_obc+1,J,K,bi,bj)-SW_STORE_1(J,K,bi,bj))/
164         &          (ab1*SW_STORE_2(J,K,bi,bj) + ab2*SW_STORE_3(J,K,bi,bj))
165                   ENDIF
166                   IF (CL.lt.0.) THEN
167                      CL=0.
168                   ELSEIF (CL.gt.CMAX) THEN
169                      CL=CMAX
170                   ENDIF
171                   CVEL_SW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)+
172         &                f2*CVEL_SW(J,K,bi,bj)
173    C              update OBC to next timestep
174                   OBWs(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)+
175         &           CVEL_SW(J,K,bi,bj)*(deltaT/dxC(I_obc+1,J,bi,bj))*
176         &           (ab1*(salt(I_obc+1,J,K,bi,bj)-salt(I_obc,J,K,bi,bj))+
177         &           ab2*(SW_STORE_1(J,K,bi,bj)-SW_STORE_4(J,K,bi,bj)))
178    C              wVel
179    #ifdef ALLOW_NONHYDROSTATIC
180                      IF ((WW_STORE_2(J,K,bi,bj).eq.0.).and.
181         &               (WW_STORE_3(J,K,bi,bj).eq.0.)) THEN
182                         CL=0.
183                      ELSE
184                        CL=(wVel(I_obc+1,J,K,bi,bj)-WW_STORE_1(J,K,bi,bj))/
185         &            (ab1*WW_STORE_2(J,K,bi,bj)+ab2*WW_STORE_3(J,K,bi,bj))
186                      ENDIF
187                      IF (CL.lt.0.) THEN
188                         CL=0.
189                      ELSEIF (CL.gt.CMAX) THEN
190                         CL=CMAX
191                      ENDIF
192                      CVEL_WW(J,K,bi,bj)=f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)
193         &                   + f2*CVEL_WW(J,K,bi,bj)
194    C                 update OBC to next timestep
195                      OBWw(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)+
196         &             CVEL_WW(J,K,bi,bj)*(deltaT/dxC(I_obc+1,J,bi,bj))*
197         &            (ab1*(wVel(I_obc+1,J,K,bi,bj)-wVel(I_obc,J,K,bi,bj))+
198         &            ab2*(WW_STORE_1(J,K,bi,bj)-WW_STORE_4(J,K,bi,bj)))
199    #endif
200    C              update/save storage arrays
201    C              uVel
202    C              copy t-1 to t-2 array
203                   UW_STORE_3(J,K,bi,bj)=UW_STORE_2(J,K,bi,bj)
204    C              copy (current time) t to t-1 arrays
205                   UW_STORE_2(J,K,bi,bj)=uVel(I_obc+3,J,K,bi,bj) -
206         &         uVel(I_obc+2,J,K,bi,bj)
207                   UW_STORE_1(J,K,bi,bj)=uVel(I_obc+2,J,K,bi,bj)
208                   UW_STORE_4(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)
209    C              vVel
210    C              copy t-1 to t-2 array
211                   VW_STORE_3(J,K,bi,bj)=VW_STORE_2(J,K,bi,bj)
212    C              copy (current time) t to t-1 arrays
213                   VW_STORE_2(J,K,bi,bj)=vVel(I_obc+2,J,K,bi,bj) -
214         &         vVel(I_obc+1,J,K,bi,bj)
215                   VW_STORE_1(J,K,bi,bj)=vVel(I_obc+1,J,K,bi,bj)
216                   VW_STORE_4(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)
217    C              Temperature
218    C              copy t-1 to t-2 array
219                   TW_STORE_3(J,K,bi,bj)=TW_STORE_2(J,K,bi,bj)
220    C              copy (current time) t to t-1 arrays
221                   TW_STORE_2(J,K,bi,bj)=theta(I_obc+2,J,K,bi,bj) -
222         &         theta(I_obc+1,J,K,bi,bj)
223                   TW_STORE_1(J,K,bi,bj)=theta(I_obc+1,J,K,bi,bj)
224                   TW_STORE_4(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)
225    C              wVel
226    #ifdef ALLOW_NONHYDROSTATIC
227    C              copy t-1 to t-2 array
228                   WW_STORE_3(J,K,bi,bj)=WW_STORE_2(J,K,bi,bj)
229    C              copy (current time) t to t-1 arrays
230                   WW_STORE_2(J,K,bi,bj)=wVel(I_obc+2,J,K,bi,bj) -
231         &         wVel(I_obc+1,J,K,bi,bj)
232                   WW_STORE_1(J,K,bi,bj)=wVel(I_obc+1,J,K,bi,bj)
233                   WW_STORE_4(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)
234    #endif
235                ENDIF
236             ENDDO
237          ENDDO
238    
239    #endif /* ALLOW_ORLANSKI */
240          RETURN
241          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22