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

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

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

Revision 1.8 - (hide annotations) (download)
Mon Sep 20 23:22:58 2004 UTC (19 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint57m_post, checkpoint57g_pre, checkpoint57s_post, checkpoint58b_post, checkpoint57b_post, checkpoint57g_post, checkpoint56b_post, checkpoint57y_post, checkpoint57r_post, checkpoint57d_post, checkpoint57i_post, checkpoint59, checkpoint58, checkpoint57, checkpoint56, checkpoint58f_post, checkpoint57n_post, checkpoint58d_post, checkpoint58a_post, checkpoint57z_post, checkpoint58y_post, checkpoint58t_post, checkpoint55i_post, checkpoint58m_post, checkpoint57l_post, checkpoint57t_post, checkpoint55c_post, checkpoint57v_post, checkpoint57f_post, checkpoint60, checkpoint61, checkpoint57a_post, checkpoint57h_pre, checkpoint58w_post, checkpoint57h_post, checkpoint57y_pre, checkpoint55g_post, checkpoint58o_post, checkpoint57c_post, checkpoint58p_post, checkpoint58q_post, checkpoint55d_post, checkpoint58e_post, checkpoint55d_pre, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, checkpoint55h_post, checkpoint58n_post, checkpoint57e_post, checkpoint55b_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint55f_post, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, eckpoint57e_pre, checkpoint58k_post, checkpoint58v_post, checkpoint56a_post, checkpoint58l_post, checkpoint57h_done, checkpoint57j_post, checkpoint57f_pre, checkpoint61f, checkpoint58g_post, checkpoint58x_post, checkpoint61n, checkpoint59j, checkpoint58h_post, checkpoint56c_post, checkpoint58j_post, checkpoint57a_pre, checkpoint55a_post, checkpoint57o_post, checkpoint61q, checkpoint57k_post, checkpoint57w_post, checkpoint61e, checkpoint58i_post, checkpoint57x_post, checkpoint58c_post, checkpoint58u_post, checkpoint58s_post, checkpoint55e_post, checkpoint61g, checkpoint61d, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p
Changes since 1.7: +3 -1 lines 
o merged code to
  * prescribe/read time-dependent open boundaries
    (works in conjunction with exf, cal)
  * sponge layer code for open boundaries
  * each boundary N/S/E/W now has its own CPP option
    (healthy for the adjoint)
1 heimbach 1.8 C $Header: /u/gcmpack/MITgcm/pkg/obcs/orlanski_west.F,v 1.7 2004/07/08 17:03:29 adcroft Exp $
2 adcroft 1.6 C $Name: $
3 adcroft 1.5 cc
4 adcroft 1.2
5     #include "OBCS_OPTIONS.h"
6    
7     SUBROUTINE ORLANSKI_WEST( bi, bj, futureTime,
8     I uVel, vVel, wVel, theta, salt,
9     I myThid )
10     C /==========================================================\
11     C | SUBROUTINE ORLANSKI_WEST |
12     C | o Calculate future boundary data at open boundaries |
13     C | at time = futureTime by applying Orlanski radiation |
14     C | conditions. |
15     C |==========================================================|
16     C | |
17     C \==========================================================/
18     IMPLICIT NONE
19    
20     C === Global variables ===
21     #include "SIZE.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "GRID.h"
25     #include "OBCS.h"
26     #include "ORLANSKI.h"
27    
28 adcroft 1.3 C SPK 6/2/00: Added radiative OBCs for salinity.
29 adcroft 1.2 C SPK 6/6/00: Changed calculation of OB*w. When K=1, the
30     C upstream value is used. For example on the eastern OB:
31     C IF (K.EQ.1) THEN
32     C OBEw(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
33     C ENDIF
34     C
35     C SPK 7/7/00: 1) Removed OB*w fix (see above).
36     C 2) Added variable CMAX. Maximum diagnosed phase speed is now
37     C clamped to CMAX. For stability of AB-II scheme (CFL) the
38     C (non-dimensional) phase speed must be <0.5
39     C 3) (Sonya Legg) Changed application of uVel and vVel.
40     C uVel on the western OB is actually applied at I_obc+1
41     C while vVel on the southern OB is applied at J_obc+1.
42 adcroft 1.3 C 4) (Sonya Legg) Added templates for forced OBs.
43 adcroft 1.2 C
44     C SPK 7/17/00: Non-uniform resolution is now taken into account in diagnosing
45     C phase speeds and time-stepping OB values. CL is still the
46     C non-dimensional phase speed; CVEL is the dimensional phase
47     C speed: CVEL = CL*(dx or dy)/dt, where dx and dy is the
48     C appropriate grid spacings. Note that CMAX (with which CL
49     C is compared) remains non-dimensional.
50     C
51     C SPK 7/18/00: Added code to allow filtering of phase speed following
52     C Blumberg and Kantha. There is now a separate array
53     C CVEL_**, where **=Variable(U,V,T,S,W)Boundary(E,W,N,S) for
54     C the dimensional phase speed. These arrays are initialized to
55     C zero in ini_obcs.F. CVEL_** is filtered according to
56     C CVEL_** = fracCVEL*CVEL(new) + (1-fracCVEL)*CVEL_**(old).
57     C fracCVEL=1.0 turns off filtering.
58     C
59     C SPK 7/26/00: Changed code to average phase speed. A new variable
60     C 'cvelTimeScale' was created. This variable must now be
61     C specified. Then, fracCVEL=deltaT/cvelTimeScale.
62     C Since the goal is to smooth out the 'singularities' in the
63     C diagnosed phase speed, cvelTimeScale could be picked as the
64     C duration of the singular period in the unfiltered case. Thus,
65     C for a plane wave cvelTimeScale might be the time take for the
66     C wave to travel a distance DX, where DX is the width of the region
67     C near which d(phi)/dx is small.
68 adcroft 1.5 C
69     C JBG 3/24/03: Fixed phase speed at western boundary (as suggested by
70     C Dale Durran's MWR paper). Fixed value (in m/s) is
71     C passed in as variable CFIX in data.obcs.
72     C
73     C JBG 4/10/03: allow choice of Orlanski or fixed wavespeed (by means of new
74     C booleans useFixedCEast and useFixedCWest) without
75     C having to recompile each time
76     c SAL 1/7/03: Fixed bug: implementation for salinity was incomplete
77     C
78 adcroft 1.2
79     C == Routine arguments ==
80     INTEGER bi, bj
81     _RL futureTime
82     _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
83     _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
84     _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
85     _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
86     _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
87     INTEGER myThid
88    
89     #ifdef ALLOW_ORLANSKI
90 heimbach 1.8 #ifdef ALLOW_OBCS_WEST
91 adcroft 1.2
92     C == Local variables ==
93     INTEGER J, K, I_obc
94     _RL CL, ab1, ab2, fracCVEL, f1, f2
95    
96     ab1 = 1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
97     ab2 = -0.5 _d 0 - abEps
98     /* CMAX is maximum allowable phase speed-CFL for AB-II */
99     /* cvelTimeScale is averaging period for phase speed in sec. */
100    
101     fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
102 adcroft 1.3 f1 = fracCVEL /* dont change this. Set cvelTimeScale */
103     f2 = 1.0-fracCVEL /* dont change this. set cvelTimeScale */
104 adcroft 1.2
105     C Western OB (Orlanski Radiation Condition)
106     DO K=1,Nr
107     DO J=1-Oly,sNy+Oly
108     I_obc=OB_Iw(J,bi,bj)
109     IF (I_obc.ne.0) THEN
110     C uVel (to be applied at I_obc+1)
111     IF ((UW_STORE_2(J,K,bi,bj).eq.0.).and.
112     & (UW_STORE_3(J,K,bi,bj).eq.0.)) THEN
113     CL=0.
114     ELSE
115     CL=(uVel(I_obc+2,J,K,bi,bj)-UW_STORE_1(J,K,bi,bj))/
116     & (ab1*UW_STORE_2(J,K,bi,bj) + ab2*UW_STORE_3(J,K,bi,bj))
117     ENDIF
118     IF (CL.lt.0.) THEN
119     CL=0.
120     ELSEIF (CL.gt.CMAX) THEN
121     CL=CMAX
122     ENDIF
123 adcroft 1.5 IF (useFixedCWest) THEN
124     C Fixed phase speed (ignoring all of that painstakingly
125     C saved data...)
126     CVEL_UW(J,K,bi,bj) = CFIX
127     ELSE
128 adcroft 1.6 CVEL_UW(J,K,bi,bj) = f1*(CL*dxF(I_obc+2,J,bi,bj)/deltaT
129     & )+f2*CVEL_UW(J,K,bi,bj)
130 adcroft 1.5 ENDIF
131 adcroft 1.2 C update OBC to next timestep
132     OBWu(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)+
133 adcroft 1.7 & CVEL_UW(J,K,bi,bj)*(deltaT*recip_dxF(I_obc+1,J,bi,bj))*
134 adcroft 1.2 & (ab1*(uVel(I_obc+2,J,K,bi,bj)-uVel(I_obc+1,J,K,bi,bj))+
135     & ab2*(UW_STORE_1(J,K,bi,bj)-UW_STORE_4(J,K,bi,bj)))
136     C vVel
137     IF ((VW_STORE_2(J,K,bi,bj).eq.0.).and.
138     & (VW_STORE_3(J,K,bi,bj).eq.0.)) THEN
139     CL=0.
140     ELSE
141     CL=(vVel(I_obc+1,J,K,bi,bj)-VW_STORE_1(J,K,bi,bj))/
142     & (ab1*VW_STORE_2(J,K,bi,bj) + ab2*VW_STORE_3(J,K,bi,bj))
143     ENDIF
144     IF (CL.lt.0.) THEN
145     CL=0.
146     ELSEIF (CL.gt.CMAX) THEN
147     CL=CMAX
148     ENDIF
149 adcroft 1.5 IF (useFixedCWest) THEN
150     C Fixed phase speed (ignoring all of that painstakingly
151     C saved data...)
152     CVEL_VW(J,K,bi,bj) = CFIX
153     ELSE
154 adcroft 1.6 CVEL_VW(J,K,bi,bj) = f1*(CL*dxV(I_obc+2,J,bi,bj)/deltaT
155     & )+f2*CVEL_VW(J,K,bi,bj)
156 adcroft 1.5 ENDIF
157 adcroft 1.2 C update OBC to next timestep
158     OBWv(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)+
159 adcroft 1.7 & CVEL_VW(J,K,bi,bj)*(deltaT*recip_dxV(I_obc+1,J,bi,bj))*
160 adcroft 1.2 & (ab1*(vVel(I_obc+1,J,K,bi,bj)-vVel(I_obc,J,K,bi,bj))+
161     & ab2*(VW_STORE_1(J,K,bi,bj)-VW_STORE_4(J,K,bi,bj)))
162     C Temperature
163     IF ((TW_STORE_2(J,K,bi,bj).eq.0.).and.
164     & (TW_STORE_3(J,K,bi,bj).eq.0.)) THEN
165     CL=0.
166     ELSE
167     CL=(theta(I_obc+1,J,K,bi,bj)-TW_STORE_1(J,K,bi,bj))/
168     & (ab1*TW_STORE_2(J,K,bi,bj) + ab2*TW_STORE_3(J,K,bi,bj))
169     ENDIF
170     IF (CL.lt.0.) THEN
171     CL=0.
172     ELSEIF (CL.gt.CMAX) THEN
173     CL=CMAX
174     ENDIF
175 adcroft 1.5 IF (useFixedCWest) THEN
176     C Fixed phase speed (ignoring all of that painstakingly
177     C saved data...)
178     CVEL_TW(J,K,bi,bj) = CFIX
179     ELSE
180 adcroft 1.6 CVEL_TW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT
181     & )+f2*CVEL_TW(J,K,bi,bj)
182 adcroft 1.5 ENDIF
183 adcroft 1.2 C update OBC to next timestep
184     OBWt(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)+
185 adcroft 1.7 & CVEL_TW(J,K,bi,bj)*(deltaT*recip_dxC(I_obc+1,J,bi,bj))*
186 adcroft 1.2 & (ab1*(theta(I_obc+1,J,K,bi,bj)-theta(I_obc,J,K,bi,bj))+
187     & ab2*(TW_STORE_1(J,K,bi,bj)-TW_STORE_4(J,K,bi,bj)))
188     C Salinity
189     IF ((SW_STORE_2(J,K,bi,bj).eq.0.).and.
190     & (SW_STORE_3(J,K,bi,bj).eq.0.)) THEN
191     CL=0.
192     ELSE
193     CL=(salt(I_obc+1,J,K,bi,bj)-SW_STORE_1(J,K,bi,bj))/
194     & (ab1*SW_STORE_2(J,K,bi,bj) + ab2*SW_STORE_3(J,K,bi,bj))
195     ENDIF
196     IF (CL.lt.0.) THEN
197     CL=0.
198     ELSEIF (CL.gt.CMAX) THEN
199     CL=CMAX
200     ENDIF
201 adcroft 1.5 IF (useFixedCWest) THEN
202     C Fixed phase speed (ignoring all of that painstakingly
203     C saved data...)
204     CVEL_SW(J,K,bi,bj) = CFIX
205     ELSE
206 adcroft 1.6 CVEL_SW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT
207     & )+f2*CVEL_SW(J,K,bi,bj)
208 adcroft 1.5 ENDIF
209 adcroft 1.2 C update OBC to next timestep
210     OBWs(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)+
211 adcroft 1.7 & CVEL_SW(J,K,bi,bj)*(deltaT*recip_dxC(I_obc+1,J,bi,bj))*
212 adcroft 1.2 & (ab1*(salt(I_obc+1,J,K,bi,bj)-salt(I_obc,J,K,bi,bj))+
213     & ab2*(SW_STORE_1(J,K,bi,bj)-SW_STORE_4(J,K,bi,bj)))
214     C wVel
215     #ifdef ALLOW_NONHYDROSTATIC
216     IF ((WW_STORE_2(J,K,bi,bj).eq.0.).and.
217     & (WW_STORE_3(J,K,bi,bj).eq.0.)) THEN
218     CL=0.
219     ELSE
220     CL=(wVel(I_obc+1,J,K,bi,bj)-WW_STORE_1(J,K,bi,bj))/
221     & (ab1*WW_STORE_2(J,K,bi,bj)+ab2*WW_STORE_3(J,K,bi,bj))
222     ENDIF
223     IF (CL.lt.0.) THEN
224     CL=0.
225     ELSEIF (CL.gt.CMAX) THEN
226     CL=CMAX
227     ENDIF
228 adcroft 1.5 IF (useFixedCWest) THEN
229     C Fixed phase speed (ignoring all of that painstakingly
230     C saved data...)
231     CVEL_WW(J,K,bi,bj) = CFIX
232     ELSE
233     CVEL_WW(J,K,bi,bj)=f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)
234 adcroft 1.2 & + f2*CVEL_WW(J,K,bi,bj)
235 adcroft 1.5 ENDIF
236 adcroft 1.2 C update OBC to next timestep
237     OBWw(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)+
238 adcroft 1.7 & CVEL_WW(J,K,bi,bj)*(deltaT*recip_dxC(I_obc+1,J,bi,bj))*
239 adcroft 1.2 & (ab1*(wVel(I_obc+1,J,K,bi,bj)-wVel(I_obc,J,K,bi,bj))+
240     & ab2*(WW_STORE_1(J,K,bi,bj)-WW_STORE_4(J,K,bi,bj)))
241     #endif
242     C update/save storage arrays
243     C uVel
244     C copy t-1 to t-2 array
245     UW_STORE_3(J,K,bi,bj)=UW_STORE_2(J,K,bi,bj)
246     C copy (current time) t to t-1 arrays
247     UW_STORE_2(J,K,bi,bj)=uVel(I_obc+3,J,K,bi,bj) -
248     & uVel(I_obc+2,J,K,bi,bj)
249     UW_STORE_1(J,K,bi,bj)=uVel(I_obc+2,J,K,bi,bj)
250     UW_STORE_4(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)
251     C vVel
252     C copy t-1 to t-2 array
253     VW_STORE_3(J,K,bi,bj)=VW_STORE_2(J,K,bi,bj)
254     C copy (current time) t to t-1 arrays
255     VW_STORE_2(J,K,bi,bj)=vVel(I_obc+2,J,K,bi,bj) -
256     & vVel(I_obc+1,J,K,bi,bj)
257     VW_STORE_1(J,K,bi,bj)=vVel(I_obc+1,J,K,bi,bj)
258     VW_STORE_4(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)
259     C Temperature
260     C copy t-1 to t-2 array
261     TW_STORE_3(J,K,bi,bj)=TW_STORE_2(J,K,bi,bj)
262     C copy (current time) t to t-1 arrays
263     TW_STORE_2(J,K,bi,bj)=theta(I_obc+2,J,K,bi,bj) -
264     & theta(I_obc+1,J,K,bi,bj)
265     TW_STORE_1(J,K,bi,bj)=theta(I_obc+1,J,K,bi,bj)
266     TW_STORE_4(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)
267 adcroft 1.5 c Salinity
268     C copy t-1 to t-2 array
269     SW_STORE_3(J,K,bi,bj)=SW_STORE_2(J,K,bi,bj)
270     C copy (current time) t to t-1 arrays
271     SW_STORE_2(J,K,bi,bj)=salt(I_obc+2,J,K,bi,bj) -
272     & salt(I_obc+1,J,K,bi,bj)
273     SW_STORE_1(J,K,bi,bj)=salt(I_obc+1,J,K,bi,bj)
274     SW_STORE_4(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)
275 adcroft 1.2 C wVel
276     #ifdef ALLOW_NONHYDROSTATIC
277     C copy t-1 to t-2 array
278     WW_STORE_3(J,K,bi,bj)=WW_STORE_2(J,K,bi,bj)
279     C copy (current time) t to t-1 arrays
280     WW_STORE_2(J,K,bi,bj)=wVel(I_obc+2,J,K,bi,bj) -
281     & wVel(I_obc+1,J,K,bi,bj)
282     WW_STORE_1(J,K,bi,bj)=wVel(I_obc+1,J,K,bi,bj)
283     WW_STORE_4(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)
284     #endif
285     ENDIF
286     ENDDO
287     ENDDO
288    
289 heimbach 1.8 #endif
290 adcroft 1.2 #endif /* ALLOW_ORLANSKI */
291     RETURN
292     END
  ViewVC Help
Powered by ViewVC 1.1.22