/[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.4 - (hide annotations) (download)
Thu Jul 11 16:22:30 2002 UTC (21 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint51k_post, checkpoint47e_post, checkpoint52l_pre, hrcube4, hrcube5, checkpoint46l_post, checkpoint46g_pre, checkpoint47c_post, checkpoint50c_post, checkpoint46f_post, checkpoint52d_pre, checkpoint48e_post, checkpoint50c_pre, checkpoint46b_post, checkpoint52j_pre, checkpoint51o_pre, checkpoint51l_post, checkpoint48i_post, checkpoint46l_pre, checkpoint52l_post, checkpoint52k_post, checkpoint54, checkpoint51, checkpoint50, checkpoint53, checkpoint52, checkpoint50d_post, checkpoint52f_post, checkpoint50b_pre, checkpoint51f_post, checkpoint48b_post, checkpoint51d_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint51t_post, checkpoint51n_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint47a_post, checkpoint48d_pre, checkpoint51j_post, checkpoint47i_post, checkpoint52e_pre, checkpoint52e_post, checkpoint51n_pre, checkpoint47d_post, checkpoint53d_post, checkpoint46d_pre, checkpoint48d_post, checkpoint48f_post, checkpoint45d_post, checkpoint52b_pre, checkpoint46j_pre, checkpoint51l_pre, checkpoint52m_post, checkpoint48h_post, checkpoint51q_post, checkpoint51b_pre, checkpoint46a_post, checkpoint47g_post, checkpoint52b_post, checkpoint52c_post, checkpoint46j_post, checkpoint51h_pre, checkpoint46k_post, checkpoint48a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint52f_pre, checkpoint47j_post, checkpoint54a_pre, checkpoint53c_post, branch-exfmods-tag, branchpoint-genmake2, checkpoint54a_post, checkpoint46e_pre, checkpoint51r_post, checkpoint48c_post, checkpoint46b_pre, checkpoint51i_post, checkpoint51b_post, checkpoint51c_post, checkpoint46c_pre, checkpoint53a_post, checkpoint46, checkpoint47b_post, checkpoint46h_pre, checkpoint52d_post, checkpoint53g_post, checkpoint46m_post, checkpoint46a_pre, checkpoint50g_post, checkpoint46g_post, checkpoint52a_pre, checkpoint50h_post, checkpoint52i_post, checkpoint50e_pre, checkpoint50i_post, checkpoint51i_pre, checkpoint52h_pre, checkpoint53f_post, checkpoint52j_post, checkpoint47f_post, checkpoint50e_post, checkpoint46i_post, checkpoint46c_post, branch-netcdf, checkpoint50d_pre, checkpoint52n_post, checkpoint53b_pre, checkpoint46e_post, checkpoint51e_post, checkpoint47, checkpoint48, checkpoint49, checkpoint46h_post, checkpoint51o_post, checkpoint51f_pre, checkpoint48g_post, checkpoint53b_post, checkpoint47h_post, checkpoint52a_post, checkpoint51g_post, ecco_c52_e35, checkpoint46d_post, checkpoint50b_post, checkpoint51m_post, checkpoint53d_pre, checkpoint51a_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-exfmods-curt, branch-genmake2, branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.3: +6 -6 lines 
use recip_dx*,recip_dy* instead of /dx*,/dy*
  (Pb reported by Arne, Thanks !)
1 jmc 1.4 C $Header: /u/gcmpack/MITgcm/pkg/obcs/orlanski_west.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_WEST( bi, bj, futureTime,
7     I uVel, vVel, wVel, theta, salt,
8     I myThid )
9     C /==========================================================\
10     C | SUBROUTINE ORLANSKI_WEST |
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 Western OB (Orlanski Radiation Condition)
94     DO K=1,Nr
95     DO J=1-Oly,sNy+Oly
96     I_obc=OB_Iw(J,bi,bj)
97     IF (I_obc.ne.0) THEN
98     C uVel (to be applied at I_obc+1)
99     IF ((UW_STORE_2(J,K,bi,bj).eq.0.).and.
100     & (UW_STORE_3(J,K,bi,bj).eq.0.)) THEN
101     CL=0.
102     ELSE
103     CL=(uVel(I_obc+2,J,K,bi,bj)-UW_STORE_1(J,K,bi,bj))/
104     & (ab1*UW_STORE_2(J,K,bi,bj) + ab2*UW_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_UW(J,K,bi,bj) = f1*(CL*dxF(I_obc+2,J,bi,bj)/deltaT)+
112     & f2*CVEL_UW(J,K,bi,bj)
113     C update OBC to next timestep
114     OBWu(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)+
115 jmc 1.4 & CVEL_UW(J,K,bi,bj)*deltaT*recip_dxF(I_obc+1,J,bi,bj)*
116 adcroft 1.2 & (ab1*(uVel(I_obc+2,J,K,bi,bj)-uVel(I_obc+1,J,K,bi,bj))+
117     & ab2*(UW_STORE_1(J,K,bi,bj)-UW_STORE_4(J,K,bi,bj)))
118     C vVel
119     IF ((VW_STORE_2(J,K,bi,bj).eq.0.).and.
120     & (VW_STORE_3(J,K,bi,bj).eq.0.)) THEN
121     CL=0.
122     ELSE
123     CL=(vVel(I_obc+1,J,K,bi,bj)-VW_STORE_1(J,K,bi,bj))/
124     & (ab1*VW_STORE_2(J,K,bi,bj) + ab2*VW_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_VW(J,K,bi,bj) = f1*(CL*dxV(I_obc+2,J,bi,bj)/deltaT)+
132     & f2*CVEL_VW(J,K,bi,bj)
133     C update OBC to next timestep
134     OBWv(J,K,bi,bj)=vVel(I_obc,J,K,bi,bj)+
135 jmc 1.4 & CVEL_VW(J,K,bi,bj)*deltaT*recip_dxV(I_obc+1,J,bi,bj)*
136 adcroft 1.2 & (ab1*(vVel(I_obc+1,J,K,bi,bj)-vVel(I_obc,J,K,bi,bj))+
137     & ab2*(VW_STORE_1(J,K,bi,bj)-VW_STORE_4(J,K,bi,bj)))
138     C Temperature
139     IF ((TW_STORE_2(J,K,bi,bj).eq.0.).and.
140     & (TW_STORE_3(J,K,bi,bj).eq.0.)) THEN
141     CL=0.
142     ELSE
143     CL=(theta(I_obc+1,J,K,bi,bj)-TW_STORE_1(J,K,bi,bj))/
144     & (ab1*TW_STORE_2(J,K,bi,bj) + ab2*TW_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_TW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)+
152     & f2*CVEL_TW(J,K,bi,bj)
153     C update OBC to next timestep
154     OBWt(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)+
155 jmc 1.4 & CVEL_TW(J,K,bi,bj)*deltaT*recip_dxC(I_obc+1,J,bi,bj)*
156 adcroft 1.2 & (ab1*(theta(I_obc+1,J,K,bi,bj)-theta(I_obc,J,K,bi,bj))+
157     & ab2*(TW_STORE_1(J,K,bi,bj)-TW_STORE_4(J,K,bi,bj)))
158     C Salinity
159     IF ((SW_STORE_2(J,K,bi,bj).eq.0.).and.
160     & (SW_STORE_3(J,K,bi,bj).eq.0.)) THEN
161     CL=0.
162     ELSE
163     CL=(salt(I_obc+1,J,K,bi,bj)-SW_STORE_1(J,K,bi,bj))/
164     & (ab1*SW_STORE_2(J,K,bi,bj) + ab2*SW_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_SW(J,K,bi,bj) = f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)+
172     & f2*CVEL_SW(J,K,bi,bj)
173     C update OBC to next timestep
174     OBWs(J,K,bi,bj)=salt(I_obc,J,K,bi,bj)+
175 jmc 1.4 & CVEL_SW(J,K,bi,bj)*deltaT*recip_dxC(I_obc+1,J,bi,bj)*
176 adcroft 1.2 & (ab1*(salt(I_obc+1,J,K,bi,bj)-salt(I_obc,J,K,bi,bj))+
177     & ab2*(SW_STORE_1(J,K,bi,bj)-SW_STORE_4(J,K,bi,bj)))
178     C wVel
179     #ifdef ALLOW_NONHYDROSTATIC
180     IF ((WW_STORE_2(J,K,bi,bj).eq.0.).and.
181     & (WW_STORE_3(J,K,bi,bj).eq.0.)) THEN
182     CL=0.
183     ELSE
184     CL=(wVel(I_obc+1,J,K,bi,bj)-WW_STORE_1(J,K,bi,bj))/
185     & (ab1*WW_STORE_2(J,K,bi,bj)+ab2*WW_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_WW(J,K,bi,bj)=f1*(CL*dxC(I_obc+2,J,bi,bj)/deltaT)
193     & + f2*CVEL_WW(J,K,bi,bj)
194     C update OBC to next timestep
195     OBWw(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)+
196 jmc 1.4 & CVEL_WW(J,K,bi,bj)*deltaT*recip_dxC(I_obc+1,J,bi,bj)*
197 adcroft 1.2 & (ab1*(wVel(I_obc+1,J,K,bi,bj)-wVel(I_obc,J,K,bi,bj))+
198     & ab2*(WW_STORE_1(J,K,bi,bj)-WW_STORE_4(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     UW_STORE_3(J,K,bi,bj)=UW_STORE_2(J,K,bi,bj)
204     C copy (current time) t to t-1 arrays
205     UW_STORE_2(J,K,bi,bj)=uVel(I_obc+3,J,K,bi,bj) -
206     & uVel(I_obc+2,J,K,bi,bj)
207     UW_STORE_1(J,K,bi,bj)=uVel(I_obc+2,J,K,bi,bj)
208     UW_STORE_4(J,K,bi,bj)=uVel(I_obc+1,J,K,bi,bj)
209     C vVel
210     C copy t-1 to t-2 array
211     VW_STORE_3(J,K,bi,bj)=VW_STORE_2(J,K,bi,bj)
212     C copy (current time) t to t-1 arrays
213     VW_STORE_2(J,K,bi,bj)=vVel(I_obc+2,J,K,bi,bj) -
214     & vVel(I_obc+1,J,K,bi,bj)
215     VW_STORE_1(J,K,bi,bj)=vVel(I_obc+1,J,K,bi,bj)
216     VW_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     TW_STORE_3(J,K,bi,bj)=TW_STORE_2(J,K,bi,bj)
220     C copy (current time) t to t-1 arrays
221     TW_STORE_2(J,K,bi,bj)=theta(I_obc+2,J,K,bi,bj) -
222     & theta(I_obc+1,J,K,bi,bj)
223     TW_STORE_1(J,K,bi,bj)=theta(I_obc+1,J,K,bi,bj)
224     TW_STORE_4(J,K,bi,bj)=theta(I_obc,J,K,bi,bj)
225     C wVel
226     #ifdef ALLOW_NONHYDROSTATIC
227     C copy t-1 to t-2 array
228     WW_STORE_3(J,K,bi,bj)=WW_STORE_2(J,K,bi,bj)
229     C copy (current time) t to t-1 arrays
230     WW_STORE_2(J,K,bi,bj)=wVel(I_obc+2,J,K,bi,bj) -
231     & wVel(I_obc+1,J,K,bi,bj)
232     WW_STORE_1(J,K,bi,bj)=wVel(I_obc+1,J,K,bi,bj)
233     WW_STORE_4(J,K,bi,bj)=wVel(I_obc,J,K,bi,bj)
234     #endif
235     ENDIF
236     ENDDO
237     ENDDO
238    
239     #endif /* ALLOW_ORLANSKI */
240     RETURN
241     END
  ViewVC Help
Powered by ViewVC 1.1.22