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

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

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

Revision 1.5 - (hide 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 heimbach 1.5 C $Header: /u/gcmpack/MITgcm/pkg/obcs/orlanski_north.F,v 1.4 2002/07/11 16:22:30 jmc Exp $
2 adcroft 1.3 C $Name: $
3 adcroft 1.2
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 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 heimbach 1.5 #ifdef ALLOW_OBCS_NORTH
80 adcroft 1.2
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 adcroft 1.3 f1 = fracCVEL /* dont change this. Set cvelTimeScale */
92     f2 = 1.0-fracCVEL /* dont change this. set cvelTimeScale */
93 adcroft 1.2
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 jmc 1.4 & CVEL_UN(I,K,bi,bj)*deltaT*recip_dyU(I,J_obc,bi,bj)*
117 adcroft 1.2 & (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 jmc 1.4 & CVEL_VN(I,K,bi,bj)*deltaT*recip_dyF(I,J_obc-1,bi,bj)*
137 adcroft 1.2 & (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 jmc 1.4 & CVEL_TN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
157 adcroft 1.2 & (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 jmc 1.4 & CVEL_SN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
177 adcroft 1.2 & (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 jmc 1.4 & CVEL_WN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
198 adcroft 1.2 & (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 heimbach 1.5 #endif /* ALLOW_OBCS_NORTH */
249 adcroft 1.2 #endif /* ALLOW_ORLANSKI */
250     RETURN
251     END
  ViewVC Help
Powered by ViewVC 1.1.22