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

Contents of /MITgcm/pkg/obcs/orlanski_north.F

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

Revision 1.5 - (show annotations) (download)
Mon Sep 20 23:22:58 2004 UTC (19 years, 9 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.4: +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 C $Header: /u/gcmpack/MITgcm/pkg/obcs/orlanski_north.F,v 1.4 2002/07/11 16:22:30 jmc Exp $
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 OBCs 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 OBs.
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 #ifdef ALLOW_OBCS_NORTH
80
81 C == Local variables ==
82 INTEGER I, K, J_obc
83 _RL CL, ab1, ab2, fracCVEL, f1, f2
84
85 ab1 = 1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
86 ab2 = -0.5 _d 0 - abEps
87 /* CMAX is maximum allowable phase speed-CFL for AB-II */
88 /* cvelTimeScale is averaging period for phase speed in sec. */
89
90 fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
91 f1 = fracCVEL /* dont change this. Set cvelTimeScale */
92 f2 = 1.0-fracCVEL /* dont change this. set cvelTimeScale */
93
94 C Northern OB (Orlanski Radiation Condition)
95 DO K=1,Nr
96 DO I=1-Olx,sNx+Olx
97 J_obc=OB_Jn(I,bi,bj)
98 IF (J_obc.ne.0) THEN
99 C uVel
100 IF ((UN_STORE_2(I,K,bi,bj).eq.0.).and.
101 & (UN_STORE_3(I,K,bi,bj).eq.0.)) THEN
102 CL=0.
103 ELSE
104 CL=-(uVel(I,J_obc-1,K,bi,bj)-UN_STORE_1(I,K,bi,bj))/
105 & (ab1*UN_STORE_2(I,K,bi,bj) + ab2*UN_STORE_3(I,K,bi,bj))
106 ENDIF
107 IF (CL.lt.0.) THEN
108 CL=0.
109 ELSEIF (CL.gt.CMAX) THEN
110 CL=CMAX
111 ENDIF
112 CVEL_UN(I,K,bi,bj) = f1*(CL*dyU(I,J_obc-1,bi,bj)/deltaT)+
113 & f2*CVEL_UN(I,K,bi,bj)
114 C update OBC to next timestep
115 OBNu(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)-
116 & CVEL_UN(I,K,bi,bj)*deltaT*recip_dyU(I,J_obc,bi,bj)*
117 & (ab1*(uVel(I,J_obc,K,bi,bj)-uVel(I,J_obc-1,K,bi,bj)) +
118 & ab2*(UN_STORE_4(I,K,bi,bj)-UN_STORE_1(I,K,bi,bj)))
119 C vVel
120 IF ((VN_STORE_2(I,K,bi,bj).eq.0.).and.
121 & (VN_STORE_3(I,K,bi,bj).eq.0.)) THEN
122 CL=0.
123 ELSE
124 CL=-(vVel(I,J_obc-1,K,bi,bj)-VN_STORE_1(I,K,bi,bj))/
125 & (ab1*VN_STORE_2(I,K,bi,bj) + ab2*VN_STORE_3(I,K,bi,bj))
126 ENDIF
127 IF (CL.lt.0.) THEN
128 CL=0.
129 ELSEIF (CL.gt.CMAX) THEN
130 CL=CMAX
131 ENDIF
132 CVEL_VN(I,K,bi,bj) = f1*(CL*dyF(I,J_obc-2,bi,bj)/deltaT)+
133 & f2*CVEL_VN(I,K,bi,bj)
134 C update OBC to next timestep
135 OBNv(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)-
136 & CVEL_VN(I,K,bi,bj)*deltaT*recip_dyF(I,J_obc-1,bi,bj)*
137 & (ab1*(vVel(I,J_obc,K,bi,bj)-vVel(I,J_obc-1,K,bi,bj)) +
138 & ab2*(VN_STORE_4(I,K,bi,bj)-VN_STORE_1(I,K,bi,bj)))
139 C Temperature
140 IF ((TN_STORE_2(I,K,bi,bj).eq.0.).and.
141 & (TN_STORE_3(I,K,bi,bj).eq.0.)) THEN
142 CL=0.
143 ELSE
144 CL=-(theta(I,J_obc-1,K,bi,bj)-TN_STORE_1(I,K,bi,bj))/
145 & (ab1*TN_STORE_2(I,K,bi,bj) + ab2*TN_STORE_3(I,K,bi,bj))
146 ENDIF
147 IF (CL.lt.0.) THEN
148 CL=0.
149 ELSEIF (CL.gt.CMAX) THEN
150 CL=CMAX
151 ENDIF
152 CVEL_TN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
153 & f2*CVEL_TN(I,K,bi,bj)
154 C update OBC to next timestep
155 OBNt(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)-
156 & CVEL_TN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
157 & (ab1*(theta(I,J_obc,K,bi,bj)-theta(I,J_obc-1,K,bi,bj))+
158 & ab2*(TN_STORE_4(I,K,bi,bj)-TN_STORE_1(I,K,bi,bj)))
159 C Salinity
160 IF ((SN_STORE_2(I,K,bi,bj).eq.0.).and.
161 & (SN_STORE_3(I,K,bi,bj).eq.0.)) THEN
162 CL=0.
163 ELSE
164 CL=-(salt(I,J_obc-1,K,bi,bj)-SN_STORE_1(I,K,bi,bj))/
165 & (ab1*SN_STORE_2(I,K,bi,bj) + ab2*SN_STORE_3(I,K,bi,bj))
166 ENDIF
167 IF (CL.lt.0.) THEN
168 CL=0.
169 ELSEIF (CL.gt.CMAX) THEN
170 CL=CMAX
171 ENDIF
172 CVEL_SN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
173 & f2*CVEL_SN(I,K,bi,bj)
174 C update OBC to next timestep
175 OBNs(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)-
176 & CVEL_SN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
177 & (ab1*(salt(I,J_obc,K,bi,bj)-salt(I,J_obc-1,K,bi,bj)) +
178 & ab2*(SN_STORE_4(I,K,bi,bj)-SN_STORE_1(I,K,bi,bj)))
179 C wVel
180 #ifdef ALLOW_NONHYDROSTATIC
181 IF ((WN_STORE_2(I,K,bi,bj).eq.0.).and.
182 & (WN_STORE_3(I,K,bi,bj).eq.0.)) THEN
183 CL=0.
184 ELSE
185 CL=-(wVel(I,J_obc-1,K,bi,bj)-WN_STORE_1(I,K,bi,bj))/
186 & (ab1*WN_STORE_2(I,K,bi,bj)+ab2*WN_STORE_3(I,K,bi,bj))
187 ENDIF
188 IF (CL.lt.0.) THEN
189 CL=0.
190 ELSEIF (CL.gt.CMAX) THEN
191 CL=CMAX
192 ENDIF
193 CVEL_WN(I,K,bi,bj)=f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)
194 & + f2*CVEL_WN(I,K,bi,bj)
195 C update OBC to next timestep
196 OBNw(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)-
197 & CVEL_WN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
198 & (ab1*(wVel(I,J_obc,K,bi,bj)-wVel(I,J_obc-1,K,bi,bj))+
199 & ab2*(WN_STORE_4(I,K,bi,bj)-WN_STORE_1(I,K,bi,bj)))
200 #endif
201 C update/save storage arrays
202 C uVel
203 C copy t-1 to t-2 array
204 UN_STORE_3(I,K,bi,bj)=UN_STORE_2(I,K,bi,bj)
205 C copy (current time) t to t-1 arrays
206 UN_STORE_2(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj) -
207 & uVel(I,J_obc-2,K,bi,bj)
208 UN_STORE_1(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj)
209 UN_STORE_4(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)
210 C vVel
211 C copy t-1 to t-2 array
212 VN_STORE_3(I,K,bi,bj)=VN_STORE_2(I,K,bi,bj)
213 C copy (current time) t to t-1 arrays
214 VN_STORE_2(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj) -
215 & vVel(I,J_obc-2,K,bi,bj)
216 VN_STORE_1(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj)
217 VN_STORE_4(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)
218 C Temperature
219 C copy t-1 to t-2 array
220 TN_STORE_3(I,K,bi,bj)=TN_STORE_2(I,K,bi,bj)
221 C copy (current time) t to t-1 arrays
222 TN_STORE_2(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj) -
223 & theta(I,J_obc-2,K,bi,bj)
224 TN_STORE_1(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj)
225 TN_STORE_4(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)
226 C Salinity
227 C copy t-1 to t-2 array
228 SN_STORE_3(I,K,bi,bj)=SN_STORE_2(I,K,bi,bj)
229 C copy (current time) t to t-1 arrays
230 SN_STORE_2(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj) -
231 & salt(I,J_obc-2,K,bi,bj)
232 SN_STORE_1(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj)
233 SN_STORE_4(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)
234 C wVel
235 #ifdef ALLOW_NONHYDROSTATIC
236 C copy t-1 to t-2 array
237 WN_STORE_3(I,K,bi,bj)=WN_STORE_2(I,K,bi,bj)
238 C copy (current time) t to t-1 arrays
239 WN_STORE_2(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj) -
240 & wVel(I,J_obc-2,K,bi,bj)
241 WN_STORE_1(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj)
242 WN_STORE_4(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)
243 #endif
244 ENDIF
245 ENDDO
246 ENDDO
247
248 #endif /* ALLOW_OBCS_NORTH */
249 #endif /* ALLOW_ORLANSKI */
250 RETURN
251 END
  ViewVC Help
Powered by ViewVC 1.1.22