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

Diff of /MITgcm/pkg/obcs/orlanski_north.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:00 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_NORTH( bi, bj, futureTime,
7         I                      uVel, vVel, wVel, theta, salt,
8         I                      myThid )
9    C     /==========================================================\
10    C     | SUBROUTINE OBCS_RADIATE                                  |
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  I, K, J_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     Northern OB (Orlanski Radiation Condition)
94           DO K=1,Nr
95             DO I=1-Olx,sNx+Olx
96                J_obc=OB_Jn(I,bi,bj)
97                IF (J_obc.ne.0) THEN
98    C              uVel
99                   IF ((UN_STORE_2(I,K,bi,bj).eq.0.).and.
100         &            (UN_STORE_3(I,K,bi,bj).eq.0.)) THEN
101                      CL=0.
102                   ELSE
103                      CL=-(uVel(I,J_obc-1,K,bi,bj)-UN_STORE_1(I,K,bi,bj))/
104         &          (ab1*UN_STORE_2(I,K,bi,bj) + ab2*UN_STORE_3(I,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_UN(I,K,bi,bj) = f1*(CL*dyU(I,J_obc-1,bi,bj)/deltaT)+
112         &                f2*CVEL_UN(I,K,bi,bj)
113    C              update OBC to next timestep
114                   OBNu(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)-
115         &           CVEL_UN(I,K,bi,bj)*(deltaT/dyU(I,J_obc,bi,bj))*
116         &           (ab1*(uVel(I,J_obc,K,bi,bj)-uVel(I,J_obc-1,K,bi,bj)) +
117         &           ab2*(UN_STORE_4(I,K,bi,bj)-UN_STORE_1(I,K,bi,bj)))
118    C              vVel
119                   IF ((VN_STORE_2(I,K,bi,bj).eq.0.).and.
120         &            (VN_STORE_3(I,K,bi,bj).eq.0.)) THEN
121                      CL=0.
122                   ELSE
123                      CL=-(vVel(I,J_obc-1,K,bi,bj)-VN_STORE_1(I,K,bi,bj))/
124         &          (ab1*VN_STORE_2(I,K,bi,bj) + ab2*VN_STORE_3(I,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_VN(I,K,bi,bj) = f1*(CL*dyF(I,J_obc-2,bi,bj)/deltaT)+
132         &                f2*CVEL_VN(I,K,bi,bj)
133    C              update OBC to next timestep
134                   OBNv(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)-
135         &           CVEL_VN(I,K,bi,bj)*(deltaT/dyF(I,J_obc-1,bi,bj))*
136         &           (ab1*(vVel(I,J_obc,K,bi,bj)-vVel(I,J_obc-1,K,bi,bj)) +
137         &           ab2*(VN_STORE_4(I,K,bi,bj)-VN_STORE_1(I,K,bi,bj)))      
138    C              Temperature
139                   IF ((TN_STORE_2(I,K,bi,bj).eq.0.).and.
140         &            (TN_STORE_3(I,K,bi,bj).eq.0.)) THEN
141                      CL=0.
142                   ELSE
143                      CL=-(theta(I,J_obc-1,K,bi,bj)-TN_STORE_1(I,K,bi,bj))/
144         &          (ab1*TN_STORE_2(I,K,bi,bj) + ab2*TN_STORE_3(I,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_TN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
152         &                f2*CVEL_TN(I,K,bi,bj)
153    C              update OBC to next timestep
154                   OBNt(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)-
155         &          CVEL_TN(I,K,bi,bj)*(deltaT/dyC(I,J_obc,bi,bj))*
156         &          (ab1*(theta(I,J_obc,K,bi,bj)-theta(I,J_obc-1,K,bi,bj))+
157         &          ab2*(TN_STORE_4(I,K,bi,bj)-TN_STORE_1(I,K,bi,bj)))
158    C              Salinity
159                   IF ((SN_STORE_2(I,K,bi,bj).eq.0.).and.
160         &            (SN_STORE_3(I,K,bi,bj).eq.0.)) THEN
161                      CL=0.
162                   ELSE
163                      CL=-(salt(I,J_obc-1,K,bi,bj)-SN_STORE_1(I,K,bi,bj))/
164         &          (ab1*SN_STORE_2(I,K,bi,bj) + ab2*SN_STORE_3(I,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_SN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
172         &                f2*CVEL_SN(I,K,bi,bj)
173    C              update OBC to next timestep
174                   OBNs(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)-
175         &          CVEL_SN(I,K,bi,bj)*(deltaT/dyC(I,J_obc,bi,bj))*
176         &          (ab1*(salt(I,J_obc,K,bi,bj)-salt(I,J_obc-1,K,bi,bj)) +
177         &          ab2*(SN_STORE_4(I,K,bi,bj)-SN_STORE_1(I,K,bi,bj)))
178    C              wVel
179    #ifdef ALLOW_NONHYDROSTATIC
180                      IF ((WN_STORE_2(I,K,bi,bj).eq.0.).and.
181         &               (WN_STORE_3(I,K,bi,bj).eq.0.)) THEN
182                         CL=0.
183                      ELSE
184                        CL=-(wVel(I,J_obc-1,K,bi,bj)-WN_STORE_1(I,K,bi,bj))/
185         &             (ab1*WN_STORE_2(I,K,bi,bj)+ab2*WN_STORE_3(I,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_WN(I,K,bi,bj)=f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)
193         &                   + f2*CVEL_WN(I,K,bi,bj)
194    C                 update OBC to next timestep
195                      OBNw(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)-
196         &             CVEL_WN(I,K,bi,bj)*(deltaT/dyC(I,J_obc,bi,bj))*
197         &             (ab1*(wVel(I,J_obc,K,bi,bj)-wVel(I,J_obc-1,K,bi,bj))+
198         &             ab2*(WN_STORE_4(I,K,bi,bj)-WN_STORE_1(I,K,bi,bj)))
199    #endif
200    C              update/save storage arrays
201    C              uVel
202    C              copy t-1 to t-2 array
203                   UN_STORE_3(I,K,bi,bj)=UN_STORE_2(I,K,bi,bj)
204    C              copy (current time) t to t-1 arrays
205                   UN_STORE_2(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj) -
206         &         uVel(I,J_obc-2,K,bi,bj)
207                   UN_STORE_1(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj)
208                   UN_STORE_4(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)
209    C              vVel
210    C              copy t-1 to t-2 array
211                   VN_STORE_3(I,K,bi,bj)=VN_STORE_2(I,K,bi,bj)
212    C              copy (current time) t to t-1 arrays
213                   VN_STORE_2(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj) -
214         &         vVel(I,J_obc-2,K,bi,bj)
215                   VN_STORE_1(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj)
216                   VN_STORE_4(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)
217    C              Temperature
218    C              copy t-1 to t-2 array
219                   TN_STORE_3(I,K,bi,bj)=TN_STORE_2(I,K,bi,bj)
220    C              copy (current time) t to t-1 arrays
221                   TN_STORE_2(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj) -
222         &         theta(I,J_obc-2,K,bi,bj)
223                   TN_STORE_1(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj)
224                   TN_STORE_4(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)
225    C              Salinity
226    C              copy t-1 to t-2 array
227                   SN_STORE_3(I,K,bi,bj)=SN_STORE_2(I,K,bi,bj)
228    C              copy (current time) t to t-1 arrays
229                   SN_STORE_2(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj) -
230         &         salt(I,J_obc-2,K,bi,bj)
231                   SN_STORE_1(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj)
232                   SN_STORE_4(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)
233    C              wVel
234    #ifdef ALLOW_NONHYDROSTATIC
235    C              copy t-1 to t-2 array
236                   WN_STORE_3(I,K,bi,bj)=WN_STORE_2(I,K,bi,bj)
237    C              copy (current time) t to t-1 arrays
238                   WN_STORE_2(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj) -
239         &         wVel(I,J_obc-2,K,bi,bj)
240                   WN_STORE_1(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj)
241                   WN_STORE_4(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)
242    #endif
243                ENDIF
244             ENDDO
245           ENDDO
246    
247    #endif /* ALLOW_ORLANSKI */
248          RETURN
249          END

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

  ViewVC Help
Powered by ViewVC 1.1.22