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

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

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


Revision 1.3.4.1 - (hide annotations) (download)
Thu Jul 11 16:38:37 2002 UTC (21 years, 10 months ago) by jmc
Branch: release1
CVS Tags: release1_p13_pre, release1_p13, release1_p8, release1_p9, release1_p5, release1_p6, release1_p7, release1_p12, release1_p10, release1_p11, release1_p16, release1_p17, release1_p14, release1_p15, release1_p12_pre
Branch point for: release1_50yr
Changes since 1.3: +6 -6 lines
use recip_dx*,recip_dy* instead of /dx*,/dy*
 (Pb reported by Arne, Thanks !)

1 jmc 1.3.4.1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/orlanski_east.F,v 1.3 2001/09/27 18:13:13 adcroft Exp $
2 adcroft 1.3 C $Name: $
3 adcroft 1.2
4     #include "OBCS_OPTIONS.h"
5    
6     SUBROUTINE ORLANSKI_EAST( bi, bj, futureTime,
7     I uVel, vVel, wVel, theta, salt,
8     I myThid )
9     C /==========================================================\
10     C | SUBROUTINE ORLANSKI_EAST |
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 adcroft 1.3 C SPK 6/2/00: Added radiative OBCs for salinity.
28 adcroft 1.2 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 adcroft 1.3 C 4) (Sonya Legg) Added templates for forced OBs.
42 adcroft 1.2 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 adcroft 1.3 f1 = fracCVEL /* dont change this. Set cvelTimeScale */
91     f2 = 1.0-fracCVEL /* dont change this. set cvelTimeScale */
92 adcroft 1.2
93     C Eastern OB (Orlanski Radiation Condition)
94     DO K=1,Nr
95     DO J=1-Oly,sNy+Oly
96     I_obc=OB_Ie(J,bi,bj)
97     IF (I_obc.ne.0) THEN
98     C uVel
99     IF ((UE_STORE_2(J,K,bi,bj).eq.0.).and.
100     & (UE_STORE_3(J,K,bi,bj).eq.0.)) THEN
101     CL=0.
102     ELSE
103     CL=-(uVel(I_obc-1,J,K,bi,bj)-UE_STORE_1(J,K,bi,bj))/
104     & (ab1*UE_STORE_2(J,K,bi,bj) + ab2*UE_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_UE(J,K,bi,bj) = f1*(CL*dxF(I_obc-2,J,bi,bj)/deltaT)+
112     & f2*CVEL_UE(J,K,bi,bj)
113     C update OBC to next timestep
114     OBEu(J,K,bi,bj)=uVel(I_obc,J,K,bi,bj)-
115 jmc 1.3.4.1 & CVEL_UE(J,K,bi,bj)*deltaT*recip_dxF(I_obc-1,J,bi,bj)*
116 adcroft 1.2 & (ab1*(uVel(I_obc,J,K,bi,bj)-uVel(I_obc-1,J,K,bi,bj)) +
117     & ab2*(UE_STORE_4(J,K,bi,bj)-UE_STORE_1(J,K,bi,bj)))
118     C vVel
119     IF ((VE_STORE_2(J,K,bi,bj).eq.0.).and.
120     & (VE_STORE_3(J,K,bi,bj).eq.0.)) THEN
121     CL=0.
122     ELSE
123     CL=-(vVel(I_obc-1,J,K,bi,bj)-VE_STORE_1(J,K,bi,bj))/
124     & (ab1*VE_STORE_2(J,K,bi,bj) + ab2*VE_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_VE(J,K,bi,bj) = f1*(CL*dxV(I_obc-1,J,bi,bj)/deltaT)+
132     & f2*CVEL_VE(J,K,bi,bj)
133     C update OBC to next timestep
134     OBEv(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)-
135 jmc 1.3.4.1 & CVEL_VE(J,K,bi,bj)*deltaT*recip_dxV(I_obc,J,bi,bj)*
136 adcroft 1.2 & (ab1*(vVel(I_obc,J,K,bi,bj)-vVel(I_obc-1,J,K,bi,bj)) +
137     & ab2*(VE_STORE_4(J,K,bi,bj)-VE_STORE_1(J,K,bi,bj)))
138     C Temperature
139     IF ((TE_STORE_2(J,K,bi,bj).eq.0.).and.
140     & (TE_STORE_3(J,K,bi,bj).eq.0.)) THEN
141     CL=0.
142     ELSE
143     CL=-(theta(I_obc-1,J,K,bi,bj)-TE_STORE_1(J,K,bi,bj))/
144     & (ab1*TE_STORE_2(J,K,bi,bj) + ab2*TE_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_TE(J,K,bi,bj) = f1*(CL*dxC(I_obc-1,J,bi,bj)/deltaT)+
152     & f2*CVEL_TE(J,K,bi,bj)
153     C update OBC to next timestep
154     OBEt(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)-
155 jmc 1.3.4.1 & CVEL_TE(J,K,bi,bj)*deltaT*recip_dxC(I_obc,J,bi,bj)*
156 adcroft 1.2 & (ab1*(theta(I_obc,J,K,bi,bj)-theta(I_obc-1,J,K,bi,bj))+
157     & ab2*(TE_STORE_4(J,K,bi,bj)-TE_STORE_1(J,K,bi,bj)))
158     C Salinity
159     IF ((SE_STORE_2(J,K,bi,bj).eq.0.).and.
160     & (SE_STORE_3(J,K,bi,bj).eq.0.)) THEN
161     CL=0.
162     ELSE
163     CL=-(salt(I_obc-1,J,K,bi,bj)-SE_STORE_1(J,K,bi,bj))/
164     & (ab1*SE_STORE_2(J,K,bi,bj) + ab2*SE_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_SE(J,K,bi,bj) = f1*(CL*dxC(I_obc-1,J,bi,bj)/deltaT)+
172     & f2*CVEL_SE(J,K,bi,bj)
173     C update OBC to next timestep
174     OBEs(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)-
175 jmc 1.3.4.1 & CVEL_SE(J,K,bi,bj)*deltaT*recip_dxC(I_obc,J,bi,bj)*
176 adcroft 1.2 & (ab1*(salt(I_obc,J,K,bi,bj)-salt(I_obc-1,J,K,bi,bj))+
177     & ab2*(SE_STORE_4(J,K,bi,bj)-SE_STORE_1(J,K,bi,bj)))
178     C wVel
179     #ifdef ALLOW_NONHYDROSTATIC
180     IF ((WE_STORE_2(J,K,bi,bj).eq.0.).and.
181     & (WE_STORE_3(J,K,bi,bj).eq.0.)) THEN
182     CL=0.
183     ELSE
184     CL=-(wVel(I_obc-1,J,K,bi,bj)-WE_STORE_1(J,K,bi,bj))/
185     & (ab1*WE_STORE_2(J,K,bi,bj)+ab2*WE_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_WE(J,K,bi,bj)=f1*(CL*dxC(I_obc-1,J,bi,bj)/deltaT)
193     & + f2*CVEL_WE(J,K,bi,bj)
194     C update OBC to next timestep
195     OBEw(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)-
196 jmc 1.3.4.1 & CVEL_WE(J,K,bi,bj)*deltaT*recip_dxC(I_obc,J,bi,bj)*
197 adcroft 1.2 & (ab1*(wVel(I_obc,J,K,bi,bj)-wVel(I_obc-1,J,K,bi,bj))+
198     & ab2*(WE_STORE_4(J,K,bi,bj)-WE_STORE_1(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     UE_STORE_3(J,K,bi,bj)=UE_STORE_2(J,K,bi,bj)
204     C copy (current time) t to t-1 arrays
205     UE_STORE_2(J,K,bi,bj)=uVel(I_obc-1,J,K,bi,bj) -
206     & uVel(I_obc-2,J,K,bi,bj)
207     UE_STORE_1(J,K,bi,bj)=uVel(I_obc-1,J,K,bi,bj)
208     UE_STORE_4(J,K,bi,bj)=uVel(I_obc,J,K,bi,bj)
209     C vVel
210     C copy t-1 to t-2 array
211     VE_STORE_3(J,K,bi,bj)=VE_STORE_2(J,K,bi,bj)
212     C copy (current time) t to t-1 arrays
213     VE_STORE_2(J,K,bi,bj)=vVel(I_obc-1,J,K,bi,bj) -
214     & vVel(I_obc-2,J,K,bi,bj)
215     VE_STORE_1(J,K,bi,bj)=vVel(I_obc-1,J,K,bi,bj)
216     VE_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     TE_STORE_3(J,K,bi,bj)=TE_STORE_2(J,K,bi,bj)
220     C copy (current time) t to t-1 arrays
221     TE_STORE_2(J,K,bi,bj)=theta(I_obc-1,J,K,bi,bj) -
222     & theta(I_obc-2,J,K,bi,bj)
223     TE_STORE_1(J,K,bi,bj)=theta(I_obc-1,J,K,bi,bj)
224     TE_STORE_4(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)
225     C Salinity
226     C copy t-1 to t-2 array
227     SE_STORE_3(J,K,bi,bj)=SE_STORE_2(J,K,bi,bj)
228     C copy (current time) t to t-1 arrays
229     SE_STORE_2(J,K,bi,bj)=salt(I_obc-1,J,K,bi,bj) -
230     & salt(I_obc-2,J,K,bi,bj)
231     SE_STORE_1(J,K,bi,bj)=salt(I_obc-1,J,K,bi,bj)
232     SE_STORE_4(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)
233     C wVel
234     #ifdef ALLOW_NONHYDROSTATIC
235     C copy t-1 to t-2 array
236     WE_STORE_3(J,K,bi,bj)=WE_STORE_2(J,K,bi,bj)
237     C copy (current time) t to t-1 arrays
238     WE_STORE_2(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj) -
239     & wVel(I_obc-2,J,K,bi,bj)
240     WE_STORE_1(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
241     WE_STORE_4(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)
242     #endif
243     ENDIF
244     ENDDO
245     ENDDO
246    
247     #endif /* ALLOW_ORLANSKI */
248     RETURN
249     END

  ViewVC Help
Powered by ViewVC 1.1.22