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

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

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


Revision 1.10 - (hide annotations) (download)
Tue Sep 18 20:09:17 2012 UTC (11 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint64, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, HEAD
Changes since 1.9: +3 -3 lines
use new parameter OB_indexNone for null index value (instead of hard-coded 0)

1 jmc 1.10 C $Header: /u/gcmpack/MITgcm/pkg/obcs/orlanski_south.F,v 1.9 2011/05/24 14:31:14 jmc Exp $
2 adcroft 1.3 C $Name: $
3 adcroft 1.2
4     #include "OBCS_OPTIONS.h"
5    
6 jmc 1.7 SUBROUTINE ORLANSKI_SOUTH( bi, bj, futureTime,
7     I uVel, vVel, wVel, theta, salt,
8 adcroft 1.2 I myThid )
9     C /==========================================================\
10     C | SUBROUTINE ORLANSKI_SOUTH |
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 jmc 1.9 #include "OBCS_PARAMS.h"
25     #include "OBCS_GRID.h"
26     #include "OBCS_FIELDS.h"
27 adcroft 1.2 #include "ORLANSKI.h"
28    
29 jmc 1.7 C SPK 6/2/00: Added radiative OBCs for salinity.
30     C SPK 6/6/00: Changed calculation of OB*w. When K=1, the
31 adcroft 1.2 C upstream value is used. For example on the eastern OB:
32     C IF (K.EQ.1) THEN
33     C OBEw(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
34     C ENDIF
35 jmc 1.7 C
36 adcroft 1.2 C SPK 7/7/00: 1) Removed OB*w fix (see above).
37     C 2) Added variable CMAX. Maximum diagnosed phase speed is now
38 jmc 1.7 C clamped to CMAX. For stability of AB-II scheme (CFL) the
39 adcroft 1.2 C (non-dimensional) phase speed must be <0.5
40     C 3) (Sonya Legg) Changed application of uVel and vVel.
41 jmc 1.7 C uVel on the western OB is actually applied at I_obc+1
42 adcroft 1.2 C while vVel on the southern OB is applied at J_obc+1.
43 adcroft 1.3 C 4) (Sonya Legg) Added templates for forced OBs.
44 adcroft 1.2 C
45     C SPK 7/17/00: Non-uniform resolution is now taken into account in diagnosing
46     C phase speeds and time-stepping OB values. CL is still the
47     C non-dimensional phase speed; CVEL is the dimensional phase
48 jmc 1.7 C speed: CVEL = CL*(dx or dy)/dt, where dx and dy is the
49     C appropriate grid spacings. Note that CMAX (with which CL
50 adcroft 1.2 C is compared) remains non-dimensional.
51 jmc 1.7 C
52     C SPK 7/18/00: Added code to allow filtering of phase speed following
53     C Blumberg and Kantha. There is now a separate array
54 adcroft 1.2 C CVEL_**, where **=Variable(U,V,T,S,W)Boundary(E,W,N,S) for
55 jmc 1.7 C the dimensional phase speed. These arrays are initialized to
56 adcroft 1.2 C zero in ini_obcs.F. CVEL_** is filtered according to
57     C CVEL_** = fracCVEL*CVEL(new) + (1-fracCVEL)*CVEL_**(old).
58     C fracCVEL=1.0 turns off filtering.
59     C
60 jmc 1.7 C SPK 7/26/00: Changed code to average phase speed. A new variable
61 adcroft 1.2 C 'cvelTimeScale' was created. This variable must now be
62 jmc 1.7 C specified. Then, fracCVEL=deltaT/cvelTimeScale.
63     C Since the goal is to smooth out the 'singularities' in the
64     C diagnosed phase speed, cvelTimeScale could be picked as the
65     C duration of the singular period in the unfiltered case. Thus,
66     C for a plane wave cvelTimeScale might be the time take for the
67 adcroft 1.2 C wave to travel a distance DX, where DX is the width of the region
68     C near which d(phi)/dx is small.
69    
70     C == Routine arguments ==
71     INTEGER bi, bj
72     _RL futureTime
73     _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
74     _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
75     _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
76     _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
77     _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
78     INTEGER myThid
79    
80     #ifdef ALLOW_ORLANSKI
81 heimbach 1.5 #ifdef ALLOW_OBCS_SOUTH
82 adcroft 1.2
83     C == Local variables ==
84     INTEGER I, K, J_obc
85     _RL CL, ab1, ab2, fracCVEL, f1, f2
86    
87     ab1 = 1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
88     ab2 = -0.5 _d 0 - abEps
89     /* CMAX is maximum allowable phase speed-CFL for AB-II */
90     /* cvelTimeScale is averaging period for phase speed in sec. */
91    
92     fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
93 adcroft 1.3 f1 = fracCVEL /* dont change this. Set cvelTimeScale */
94     f2 = 1.0-fracCVEL /* dont change this. set cvelTimeScale */
95 adcroft 1.2
96     C Southern OB (Orlanski Radiation Condition)
97     DO K=1,Nr
98 jmc 1.10 DO I=1-OLx,sNx+OLx
99 adcroft 1.2 J_obc=OB_Js(I,bi,bj)
100 jmc 1.10 IF ( J_obc.NE.OB_indexNone ) THEN
101 adcroft 1.2 C uVel
102     IF ((US_STORE_2(I,K,bi,bj).eq.0.).and.
103     & (US_STORE_3(I,K,bi,bj).eq.0.)) THEN
104     CL=0.
105     ELSE
106     CL=(uVel(I,J_obc+1,K,bi,bj)-US_STORE_1(I,K,bi,bj))/
107     & (ab1*US_STORE_2(I,K,bi,bj) + ab2*US_STORE_3(I,K,bi,bj))
108     ENDIF
109     IF (CL.lt.0.) THEN
110     CL=0.
111     ELSEIF (CL.gt.CMAX) THEN
112     CL=CMAX
113     ENDIF
114     CVEL_US(I,K,bi,bj) = f1*(CL*dyU(I,J_obc+2,bi,bj)/deltaT)+
115     & f2*CVEL_US(I,K,bi,bj)
116     C update OBC to next timestep
117     OBSu(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)+
118 jmc 1.4 & CVEL_US(I,K,bi,bj)*deltaT*recip_dyU(I,J_obc+1,bi,bj)*
119 adcroft 1.2 & (ab1*(uVel(I,J_obc+1,K,bi,bj)-uVel(I,J_obc,K,bi,bj)) +
120     & ab2*(US_STORE_1(I,K,bi,bj)-US_STORE_4(I,K,bi,bj)))
121     C vVel (to be applied at J_obc+1)
122     IF ((VS_STORE_2(I,K,bi,bj).eq.0.).and.
123     & (VS_STORE_3(I,K,bi,bj).eq.0.)) THEN
124     CL=0.
125     ELSE
126     CL=(vVel(I,J_obc+2,K,bi,bj)-VS_STORE_1(I,K,bi,bj))/
127     & (ab1*VS_STORE_2(I,K,bi,bj) + ab2*VS_STORE_3(I,K,bi,bj))
128     ENDIF
129     IF (CL.lt.0.) THEN
130     CL=0.
131     ELSEIF (CL.gt.CMAX) THEN
132     CL=CMAX
133     ENDIF
134     CVEL_VS(I,K,bi,bj) = f1*(CL*dyF(I,J_obc+2,bi,bj)/deltaT)+
135     & f2*CVEL_VS(I,K,bi,bj)
136     C update OBC to next timestep
137     OBSv(I,K,bi,bj)=vVel(I,J_obc+1,K,bi,bj)+
138 jmc 1.4 & CVEL_VS(I,K,bi,bj)*deltaT*recip_dyF(I,J_obc+1,bi,bj)*
139 adcroft 1.2 & (ab1*(vVel(I,J_obc+2,K,bi,bj)-vVel(I,J_obc+1,K,bi,bj))+
140     & ab2*(VS_STORE_1(I,K,bi,bj)-VS_STORE_4(I,K,bi,bj)))
141     C Temperature
142     IF ((TS_STORE_2(I,K,bi,bj).eq.0.).and.
143     & (TS_STORE_3(I,K,bi,bj).eq.0.)) THEN
144     CL=0.
145     ELSE
146     CL=(theta(I,J_obc+1,K,bi,bj)-TS_STORE_1(I,K,bi,bj))/
147     & (ab1*TS_STORE_2(I,K,bi,bj) + ab2*TS_STORE_3(I,K,bi,bj))
148     ENDIF
149     IF (CL.lt.0.) THEN
150     CL=0.
151     ELSEIF (CL.gt.CMAX) THEN
152     CL=CMAX
153     ENDIF
154     CVEL_TS(I,K,bi,bj) = f1*(CL*dyC(I,J_obc+2,bi,bj)/deltaT)+
155     & f2*CVEL_TS(I,K,bi,bj)
156     C update OBC to next timestep
157     OBSt(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)+
158 jmc 1.4 & CVEL_TS(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc+1,bi,bj)*
159 adcroft 1.2 & (ab1*(theta(I,J_obc+1,K,bi,bj)-theta(I,J_obc,K,bi,bj))+
160     & ab2*(TS_STORE_1(I,K,bi,bj)-TS_STORE_4(I,K,bi,bj)))
161     C Salinity
162     IF ((SS_STORE_2(I,K,bi,bj).eq.0.).and.
163     & (SS_STORE_3(I,K,bi,bj).eq.0.)) THEN
164     CL=0.
165     ELSE
166     CL=(salt(I,J_obc+1,K,bi,bj)-SS_STORE_1(I,K,bi,bj))/
167     & (ab1*SS_STORE_2(I,K,bi,bj) + ab2*SS_STORE_3(I,K,bi,bj))
168     ENDIF
169     IF (CL.lt.0.) THEN
170     CL=0.
171     ELSEIF (CL.gt.CMAX) THEN
172     CL=CMAX
173     ENDIF
174     CVEL_SS(I,K,bi,bj) = f1*(CL*dyC(I,J_obc+2,bi,bj)/deltaT)+
175     & f2*CVEL_SS(I,K,bi,bj)
176     C update OBC to next timestep
177     OBSs(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)+
178 jmc 1.4 & CVEL_SS(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc+1,bi,bj)*
179 adcroft 1.2 & (ab1*(salt(I,J_obc+1,K,bi,bj)-salt(I,J_obc,K,bi,bj)) +
180     & ab2*(SS_STORE_1(I,K,bi,bj)-SS_STORE_4(I,K,bi,bj)))
181 jmc 1.8 #ifdef ALLOW_NONHYDROSTATIC
182     IF ( nonHydrostatic ) THEN
183 adcroft 1.2 C wVel
184 jmc 1.8 IF ((WS_STORE_2(I,K,bi,bj).eq.0.).and.
185     & (WS_STORE_3(I,K,bi,bj).eq.0.)) THEN
186     CL=0.
187     ELSE
188     CL=(wVel(I,J_obc+1,K,bi,bj)-WS_STORE_1(I,K,bi,bj))/
189     & (ab1*WS_STORE_2(I,K,bi,bj)+ab2*WS_STORE_3(I,K,bi,bj))
190     ENDIF
191     IF (CL.lt.0.) THEN
192     CL=0.
193     ELSEIF (CL.gt.CMAX) THEN
194     CL=CMAX
195     ENDIF
196     CVEL_WS(I,K,bi,bj)=f1*(CL*dyC(I,J_obc+2,bi,bj)/deltaT)
197 adcroft 1.2 & + f2*CVEL_WS(I,K,bi,bj)
198 jmc 1.8 C update OBC to next timestep
199     OBSw(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)+
200     & CVEL_WS(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc+1,bi,bj)*
201     & (ab1*(wVel(I,J_obc+1,K,bi,bj)-wVel(I,J_obc,K,bi,bj))+
202     & ab2*(WS_STORE_1(I,K,bi,bj)-WS_STORE_4(I,K,bi,bj)))
203     ENDIF
204     #endif /* ALLOW_NONHYDROSTATIC */
205 adcroft 1.2 C update/save storage arrays
206     C uVel
207     C copy t-1 to t-2 array
208     US_STORE_3(I,K,bi,bj)=US_STORE_2(I,K,bi,bj)
209     C copy (current time) t to t-1 arrays
210     US_STORE_2(I,K,bi,bj)=uVel(I,J_obc+2,K,bi,bj) -
211     & uVel(I,J_obc+1,K,bi,bj)
212     US_STORE_1(I,K,bi,bj)=uVel(I,J_obc+1,K,bi,bj)
213     US_STORE_4(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)
214     C vVel
215     C copy t-1 to t-2 array
216     VS_STORE_3(I,K,bi,bj)=VS_STORE_2(I,K,bi,bj)
217     C copy (current time) t to t-1 arrays
218     VS_STORE_2(I,K,bi,bj)=vVel(I,J_obc+3,K,bi,bj) -
219     & vVel(I,J_obc+2,K,bi,bj)
220     VS_STORE_1(I,K,bi,bj)=vVel(I,J_obc+2,K,bi,bj)
221     VS_STORE_4(I,K,bi,bj)=vVel(I,J_obc+1,K,bi,bj)
222     C Temperature
223     C copy t-1 to t-2 array
224     TS_STORE_3(I,K,bi,bj)=TS_STORE_2(I,K,bi,bj)
225     C copy (current time) t to t-1 arrays
226     TS_STORE_2(I,K,bi,bj)=theta(I,J_obc+2,K,bi,bj) -
227     & theta(I,J_obc+1,K,bi,bj)
228     TS_STORE_1(I,K,bi,bj)=theta(I,J_obc+1,K,bi,bj)
229     TS_STORE_4(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)
230     C Salinity
231     C copy t-1 to t-2 array
232     SS_STORE_3(I,K,bi,bj)=SS_STORE_2(I,K,bi,bj)
233     C copy (current time) t to t-1 arrays
234     SS_STORE_2(I,K,bi,bj)=salt(I,J_obc+2,K,bi,bj) -
235     & salt(I,J_obc+1,K,bi,bj)
236     SS_STORE_1(I,K,bi,bj)=salt(I,J_obc+1,K,bi,bj)
237     SS_STORE_4(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)
238 jmc 1.8 #ifdef ALLOW_NONHYDROSTATIC
239     IF ( nonHydrostatic ) THEN
240 adcroft 1.2 C wVel
241     C copy t-1 to t-2 array
242     WS_STORE_3(I,K,bi,bj)=WS_STORE_2(I,K,bi,bj)
243     C copy (current time) t to t-1 arrays
244     WS_STORE_2(I,K,bi,bj)=wVel(I,J_obc+2,K,bi,bj) -
245     & wVel(I,J_obc+1,K,bi,bj)
246     WS_STORE_1(I,K,bi,bj)=wVel(I,J_obc+1,K,bi,bj)
247     WS_STORE_4(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)
248 jmc 1.8 ENDIF
249     #endif /* ALLOW_NONHYDROSTATIC */
250 adcroft 1.2 ENDIF
251     ENDDO
252     ENDDO
253    
254 jmc 1.6 #endif /* ALLOW_OBCS_SOUTH */
255 adcroft 1.2 #endif /* ALLOW_ORLANSKI */
256     RETURN
257     END

  ViewVC Help
Powered by ViewVC 1.1.22