-
Notifications
You must be signed in to change notification settings - Fork 237
/
orlanski_north.F
243 lines (229 loc) · 10.5 KB
/
orlanski_north.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
#include "OBCS_OPTIONS.h"
SUBROUTINE ORLANSKI_NORTH( bi, bj, futureTime,
I uVel, vVel, wVel, theta, salt,
I myThid )
C *==========================================================*
C | SUBROUTINE OBCS_RADIATE |
C | o Calculate future boundary data at open boundaries |
C | at time = futureTime by applying Orlanski radiation |
C | conditions. |
C *==========================================================*
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "ORLANSKI.h"
C SPK 6/2/00: Added radiative OBCs for salinity.
C SPK 6/6/00: Changed calculation of OB*w. When K=1, the
C upstream value is used. For example on the eastern OB:
C IF (K.EQ.1) THEN
C OBEw(J,K,bi,bj)=wVel(I_obc-1,J,K,bi,bj)
C ENDIF
C
C SPK 7/7/00: 1) Removed OB*w fix (see above).
C 2) Added variable CMAX. Maximum diagnosed phase speed is now
C clamped to CMAX. For stability of AB-II scheme (CFL) the
C (non-dimensional) phase speed must be <0.5
C 3) (Sonya Legg) Changed application of uVel and vVel.
C uVel on the western OB is actually applied at I_obc+1
C while vVel on the southern OB is applied at J_obc+1.
C 4) (Sonya Legg) Added templates for forced OBs.
C
C SPK 7/17/00: Non-uniform resolution is now taken into account in diagnosing
C phase speeds and time-stepping OB values. CL is still the
C non-dimensional phase speed; CVEL is the dimensional phase
C speed: CVEL = CL*(dx or dy)/dt, where dx and dy is the
C appropriate grid spacings. Note that CMAX (with which CL
C is compared) remains non-dimensional.
C
C SPK 7/18/00: Added code to allow filtering of phase speed following
C Blumberg and Kantha. There is now a separate array
C CVEL_**, where **=Variable(U,V,T,S,W)Boundary(E,W,N,S) for
C the dimensional phase speed. These arrays are initialized to
C zero in ini_obcs.F. CVEL_** is filtered according to
C CVEL_** = fracCVEL*CVEL(new) + (1-fracCVEL)*CVEL_**(old).
C fracCVEL=1.0 turns off filtering.
C
C SPK 7/26/00: Changed code to average phase speed. A new variable
C 'cvelTimeScale' was created. This variable must now be
C specified. Then, fracCVEL=deltaT/cvelTimeScale.
C Since the goal is to smooth out the 'singularities' in the
C diagnosed phase speed, cvelTimeScale could be picked as the
C duration of the singular period in the unfiltered case. Thus,
C for a plane wave cvelTimeScale might be the time take for the
C wave to travel a distance DX, where DX is the width of the region
C near which d(phi)/dx is small.
C == Routine arguments ==
INTEGER bi, bj
_RL futureTime
_RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
INTEGER myThid
#ifdef ALLOW_ORLANSKI
#ifdef ALLOW_OBCS_NORTH
C == Local variables ==
INTEGER I, K, J_obc
_RL CL, ab1, ab2, fracCVEL, f1, f2
_RL denom
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ab1 = 1.5 _d 0 + abEps /* Adams-Bashforth coefficients */
ab2 = -0.5 _d 0 - abEps
/* CMAX is maximum allowable phase speed-CFL for AB-II */
/* cvelTimeScale is averaging period for phase speed in sec. */
fracCVEL = deltaT/cvelTimeScale /* fraction of new phase speed used*/
f1 = fracCVEL /* dont change this. Set cvelTimeScale */
f2 = 1.0-fracCVEL /* dont change this. set cvelTimeScale */
C Northern OB (Orlanski Radiation Condition)
DO K=1,Nr
DO I=1-OLx,sNx+OLx
J_obc=OB_Jn(I,bi,bj)
IF ( J_obc.NE.OB_indexNone ) THEN
C- uVel
denom =
& (ab1*UN_STORE_2(I,K,bi,bj) + ab2*UN_STORE_3(I,K,bi,bj))
IF ( denom .NE. 0 _d 0 ) THEN
CL =-(uVel(I,J_obc-1,K,bi,bj)-UN_STORE_1(I,K,bi,bj))
& / denom
CL = MIN( MAX( CL, zeroRL ), CMAX )
ELSE
CL = 0. _d 0
ENDIF
CVEL_UN(I,K,bi,bj) = f1*(CL*dyU(I,J_obc-1,bi,bj)/deltaT)+
& f2*CVEL_UN(I,K,bi,bj)
C update OBC to next timestep
OBNu(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)-
& CVEL_UN(I,K,bi,bj)*deltaT*recip_dyU(I,J_obc,bi,bj)*
& (ab1*(uVel(I,J_obc,K,bi,bj)-uVel(I,J_obc-1,K,bi,bj)) +
& ab2*(UN_STORE_4(I,K,bi,bj)-UN_STORE_1(I,K,bi,bj)))
C- vVel
denom =
& (ab1*VN_STORE_2(I,K,bi,bj) + ab2*VN_STORE_3(I,K,bi,bj))
IF ( denom .NE. 0 _d 0 ) THEN
CL =-(vVel(I,J_obc-1,K,bi,bj)-VN_STORE_1(I,K,bi,bj))
& / denom
CL = MIN( MAX( CL, zeroRL ), CMAX )
ELSE
CL = 0. _d 0
ENDIF
CVEL_VN(I,K,bi,bj) = f1*(CL*dyF(I,J_obc-2,bi,bj)/deltaT)+
& f2*CVEL_VN(I,K,bi,bj)
C update OBC to next timestep
OBNv(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)-
& CVEL_VN(I,K,bi,bj)*deltaT*recip_dyF(I,J_obc-1,bi,bj)*
& (ab1*(vVel(I,J_obc,K,bi,bj)-vVel(I,J_obc-1,K,bi,bj)) +
& ab2*(VN_STORE_4(I,K,bi,bj)-VN_STORE_1(I,K,bi,bj)))
C- Temperature
denom =
& (ab1*TN_STORE_2(I,K,bi,bj) + ab2*TN_STORE_3(I,K,bi,bj))
IF ( denom .NE. 0 _d 0 ) THEN
CL =-(theta(I,J_obc-1,K,bi,bj)-TN_STORE_1(I,K,bi,bj))
& / denom
CL = MIN( MAX( CL, zeroRL ), CMAX )
ELSE
CL = 0. _d 0
ENDIF
CVEL_TN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
& f2*CVEL_TN(I,K,bi,bj)
C update OBC to next timestep
OBNt(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)-
& CVEL_TN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
& (ab1*(theta(I,J_obc,K,bi,bj)-theta(I,J_obc-1,K,bi,bj))+
& ab2*(TN_STORE_4(I,K,bi,bj)-TN_STORE_1(I,K,bi,bj)))
C- Salinity
denom =
& (ab1*SN_STORE_2(I,K,bi,bj) + ab2*SN_STORE_3(I,K,bi,bj))
IF ( denom .NE. 0 _d 0 ) THEN
CL =-(salt(I,J_obc-1,K,bi,bj)-SN_STORE_1(I,K,bi,bj))
& / denom
CL = MIN( MAX( CL, zeroRL ), CMAX )
ELSE
CL = 0. _d 0
ENDIF
CVEL_SN(I,K,bi,bj) = f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)+
& f2*CVEL_SN(I,K,bi,bj)
C update OBC to next timestep
OBNs(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)-
& CVEL_SN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
& (ab1*(salt(I,J_obc,K,bi,bj)-salt(I,J_obc-1,K,bi,bj)) +
& ab2*(SN_STORE_4(I,K,bi,bj)-SN_STORE_1(I,K,bi,bj)))
#ifdef ALLOW_NONHYDROSTATIC
IF ( nonHydrostatic ) THEN
C- wVel
denom =
& (ab1*WN_STORE_2(I,K,bi,bj)+ab2*WN_STORE_3(I,K,bi,bj))
IF ( denom .NE. 0 _d 0 ) THEN
CL =-(wVel(I,J_obc-1,K,bi,bj)-WN_STORE_1(I,K,bi,bj))
& / denom
CL = MIN( MAX( CL, zeroRL ), CMAX )
ELSE
CL = 0. _d 0
ENDIF
CVEL_WN(I,K,bi,bj)=f1*(CL*dyC(I,J_obc-1,bi,bj)/deltaT)
& + f2*CVEL_WN(I,K,bi,bj)
C update OBC to next timestep
OBNw(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)-
& CVEL_WN(I,K,bi,bj)*deltaT*recip_dyC(I,J_obc,bi,bj)*
& (ab1*(wVel(I,J_obc,K,bi,bj)-wVel(I,J_obc-1,K,bi,bj))+
& ab2*(WN_STORE_4(I,K,bi,bj)-WN_STORE_1(I,K,bi,bj)))
ENDIF
#endif /* ALLOW_NONHYDROSTATIC */
C- update/save storage arrays
C uVel
C copy t-1 to t-2 array
UN_STORE_3(I,K,bi,bj)=UN_STORE_2(I,K,bi,bj)
C copy (current time) t to t-1 arrays
UN_STORE_2(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj) -
& uVel(I,J_obc-2,K,bi,bj)
UN_STORE_1(I,K,bi,bj)=uVel(I,J_obc-1,K,bi,bj)
UN_STORE_4(I,K,bi,bj)=uVel(I,J_obc,K,bi,bj)
C vVel
C copy t-1 to t-2 array
VN_STORE_3(I,K,bi,bj)=VN_STORE_2(I,K,bi,bj)
C copy (current time) t to t-1 arrays
VN_STORE_2(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj) -
& vVel(I,J_obc-2,K,bi,bj)
VN_STORE_1(I,K,bi,bj)=vVel(I,J_obc-1,K,bi,bj)
VN_STORE_4(I,K,bi,bj)=vVel(I,J_obc,K,bi,bj)
C Temperature
C copy t-1 to t-2 array
TN_STORE_3(I,K,bi,bj)=TN_STORE_2(I,K,bi,bj)
C copy (current time) t to t-1 arrays
TN_STORE_2(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj) -
& theta(I,J_obc-2,K,bi,bj)
TN_STORE_1(I,K,bi,bj)=theta(I,J_obc-1,K,bi,bj)
TN_STORE_4(I,K,bi,bj)=theta(I,J_obc,K,bi,bj)
C Salinity
C copy t-1 to t-2 array
SN_STORE_3(I,K,bi,bj)=SN_STORE_2(I,K,bi,bj)
C copy (current time) t to t-1 arrays
SN_STORE_2(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj) -
& salt(I,J_obc-2,K,bi,bj)
SN_STORE_1(I,K,bi,bj)=salt(I,J_obc-1,K,bi,bj)
SN_STORE_4(I,K,bi,bj)=salt(I,J_obc,K,bi,bj)
#ifdef ALLOW_NONHYDROSTATIC
IF ( nonHydrostatic ) THEN
C wVel
C copy t-1 to t-2 array
WN_STORE_3(I,K,bi,bj)=WN_STORE_2(I,K,bi,bj)
C copy (current time) t to t-1 arrays
WN_STORE_2(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj) -
& wVel(I,J_obc-2,K,bi,bj)
WN_STORE_1(I,K,bi,bj)=wVel(I,J_obc-1,K,bi,bj)
WN_STORE_4(I,K,bi,bj)=wVel(I,J_obc,K,bi,bj)
ENDIF
#endif /* ALLOW_NONHYDROSTATIC */
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_OBCS_NORTH */
#endif /* ALLOW_ORLANSKI */
RETURN
END