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

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

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

Revision 1.13 - (show annotations) (download)
Tue Sep 18 20:09:17 2012 UTC (11 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: 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, checkpoint64, checkpoint65, checkpoint65h, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e
Changes since 1.12: +79 -77 lines 
use new parameter OB_indexNone for null index value (instead of hard-coded 0)
1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_calc_stevens.F,v 1.12 2012/01/03 11:06:58 mlosch Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5 #undef CHECK_BALANCE
6
7 CBOP
8 C !ROUTINE: OBCS_CALC_STEVENS
9 C !INTERFACE:
10 SUBROUTINE OBCS_CALC_STEVENS(
11 I futureTime, futureIter,
12 I myThid )
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | SUBROUTINE OBCS_CALC_STEVENS
16 C | o Calculate future boundary data at open boundaries
17 C | at time = futureTime
18 C | from input data following Stevens(1990), and some
19 C | MOM3 legacy code
20 C |
21 C | o the code works like this
22 C | - the "barotropic" (= vertically averaged) velocity
23 C | normal to the boundary is assumed to be in
24 C | OBE/W/N/Su/v (normal) when this routine is entered
25 C | - the vertically averaged velocity is corrected
26 C | by the "baroclinic" (= deviation from vertically
27 C | averaged velocity) velocity to give a new OB?u/v;
28 C | the "barolinic" velocity is estimated from the previous
29 C | time step which makes this boundary condition depend on
30 C | a restart file. If OBCS_STEVENS_USE_INTERIOR_VELOCITY
31 C | is defined the velocity is simply copied from the model
32 C | interior to the boundary, thereby avoiding a restart
33 C | file or complicated reconstruction, but this solution
34 C | can give unexpected results.
35 C | (Note: in this context the terms barotropic and baroclinic
36 C | are MOM jargon and --- to my mind ---- should not be used)
37 C | - a wave phase speed is estimated from temporal and
38 C | horizontal variations of the tracer fields for each
39 C | tracer individually, this similar to Orlanski BCs,
40 C | but for simplicity the fields of the previous time step
41 C | are used, and the time derivative is estimated
42 C | independently of the time stepping procedure by simple
43 C | differencing
44 C | - velocity tangential to the boundary is always zero
45 C | (although this could be changed)
46 C | - a new tracer is computed from local advection equation
47 C | with an upwind scheme: tracer from the interior is
48 C | advected out of the domain, and tracer from the boundary
49 C | is "advected" into the domain by a restoring mechanism
50 C | - for the advection equation only values from the
51 C | the current (not the updated) time level are used
52 C |
53 C *==========================================================*
54 C | Feb, 2009: started by Martin Losch (Martin.Losch@awi.de)
55 C *==========================================================*
56 C \ev
57 C !USES:
58 IMPLICIT NONE
59 C === Global variables ===
60 #include "SIZE.h"
61 #include "EEPARAMS.h"
62 #include "PARAMS.h"
63 #include "GRID.h"
64 #include "OBCS_PARAMS.h"
65 #include "OBCS_GRID.h"
66 #include "OBCS_FIELDS.h"
67 #include "DYNVARS.h"
68 #ifdef ALLOW_PTRACERS
69 #include "PTRACERS_SIZE.h"
70 #include "PTRACERS_PARAMS.h"
71 #include "PTRACERS_FIELDS.h"
72 #include "OBCS_PTRACERS.h"
73 #endif /* ALLOW_PTRACERS */
74 #ifdef ALLOW_AUTODIFF_TAMC
75 #include "tamc.h"
76 #include "tamc_keys.h"
77 #endif /* ALLOW_AUTODIFF_TAMC */
78
79 C !INPUT/OUTPUT PARAMETERS:
80 C == Routine arguments ==
81 _RL futureTime
82 INTEGER futureIter
83 INTEGER myThid
84
85 #ifdef ALLOW_OBCS_STEVENS
86
87 C !LOCAL VARIABLES:
88 C == Local variables ==
89 C I,J,K :: loop indices
90 C msgBuf :: Informational/error message buffer
91 C uMer/vZonBar :: vertically averaged velocity at open boundary
92 C drFBar :: local depth for vertical average
93 C uMer/vZonPri :: velocity anomalies applied to the open boundaries
94 C gammat/s :: restoring parameters (1./(T/SrelaxStevens - time scale))
95 C u/vPhase :: estimate of phase velocity for radiation condition
96 C auxillary variables
97 C cflMer/Zon :: ratio of grid spacing and time step
98 C dtracSpace :: horizontal difference of tracer
99 C dtracTime :: temporal difference of tracer
100 C aFac :: switch (0 or 1) that turns on advective contribution
101 C gFacM/Z :: switch (0 or 1) that turns on restoring boundary condition
102 C pFac :: switch that turns on/off phase velocity contribution
103 INTEGER bi, bj
104 INTEGER I, J, K
105 c CHARACTER*(MAX_LEN_MBUF) msgBuf
106 _RL uPhase
107 _RL vPhase
108 _RL dtracSpace
109 _RL dTracTime
110 _RL cflMer (1-OLy:sNy+OLy,Nr)
111 _RL gFacM (1-OLy:sNy+OLy,Nr)
112 _RL uMerPri(Nr)
113 _RL uMerBar
114 _RL cflZon (1-OLx:sNx+OLx,Nr)
115 _RL gFacZ (1-OLx:sNx+OLx,Nr)
116 _RL vZonPri(Nr)
117 _RL vZonBar
118 _RL drFBar
119 _RL gammat, gammas, pFac, aFac
120 #ifdef ALLOW_PTRACERS
121 c INTEGER iTracer
122 #endif /* ALLOW_PTRACERS */
123 #ifdef CHECK_BALANCE
124 _RL uVelLoc, vVelLoc
125 _RL vPhase
126 #endif
127 CEOP
128
129 #ifdef ALLOW_DEBUG
130 IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC_STEVENS',myThid)
131 #endif
132
133 aFac = 1. _d 0
134 IF (.NOT. useStevensAdvection ) aFac = 0. _d 0
135 pFac = 1. _d 0
136 IF (.NOT. useStevensPhaseVel ) pFac = 0. _d 0
137 gammat = 0. _d 0
138 IF (TrelaxStevens .GT. 0. _d 0) gammat = 1./TrelaxStevens
139 gammas = 0. _d 0
140 IF (SrelaxStevens .GT. 0. _d 0) gammas = 1./SrelaxStevens
141
142 DO bj=myByLo(myThid),myByHi(myThid)
143 DO bi=myBxLo(myThid),myBxHi(myThid)
144
145 #ifdef ALLOW_AUTODIFF_TAMC
146 act1 = bi - myBxLo(myThid)
147 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
148 act2 = bj - myByLo(myThid)
149 max2 = myByHi(myThid) - myByLo(myThid) + 1
150 act3 = myThid - 1
151 max3 = nTx*nTy
152 act4 = ikey_dynamics - 1
153 ikey = (act1 + 1) + act2*max1
154 & + act3*max1*max2
155 & + act4*max1*max2*max3
156 #endif /* ALLOW_AUTODIFF_TAMC */
157
158 #ifdef ALLOW_OBCS_EAST
159
160 # ifdef ALLOW_AUTODIFF_TAMC
161 CADJ STORE OBEt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
162 CADJ STORE OBEs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
163 # endif
164 IF ( useStevensEast ) THEN
165 C Eastern OB
166 #ifdef ALLOW_DEBUG
167 IF (debugMode)
168 & CALL DEBUG_MSG('OBCS_CALC_STEVENS: East',myThid)
169 #endif
170 C compute vertical average and deviation from vertical
171 C average for I_obe
172 DO J=1-OLy,sNy+OLy
173 I = OB_Ie(J,bi,bj)
174 IF ( I.NE.OB_indexNone ) THEN
175 C first initialize some fields
176 drFbar = 0. _d 0
177 uMerBar = 0. _d 0
178 DO K=1,Nr
179 uMerPri(K) = 0. _d 0
180 ENDDO
181 DO K=1,Nr
182 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
183 uMerBar = uMerBar + uVel(I-1,J,K,bi,bj)
184 #else
185 uMerBar = uMerBar + OBEuStevens(J,K,bi,bj)
186 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
187 & *drF(K)* _hFacW(I,J,K,bi,bj)
188 drFBar = drFBar + drF(K)* _hFacW(I,J,K,bi,bj)
189 ENDDO
190 IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
191 DO K=1,Nr
192 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
193 uMerPri(K) = (uVel(I-1,J,K,bi,bj)-uMerBar)
194 #else
195 uMerPri(K) = (OBEuStevens(J,K,bi,bj)-uMerBar)
196 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
197 & * _maskW(I,J,K,bi,bj)
198 ENDDO
199 C vertical average of input field
200 drFbar = 0. _d 0
201 uMerBar = 0. _d 0
202 DO K=1,Nr
203 uMerBar = uMerBar + OBEu(J,K,bi,bj)
204 & *drF(K)* _hFacW(I,J,K,bi,bj)
205 drFBar = drFBar + drF(K)* _hFacW(I,J,K,bi,bj)
206 ENDDO
207 IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
208 C Now the absolute velocity normal to the boundary is
209 C uMerBar + uMerPri(K).
210 DO K=1,Nr
211 OBEu(J,K,bi,bj) = (uMerBar + uMerPri(K))
212 & * _maskW(I,J,K,bi,bj)
213 CML OBEv(J,K,bi,bj) = 0. _d 0
214 #ifdef ALLOW_NONHYDROSTATIC
215 OBEw(J,K,bi,bj)=0.
216 #endif
217 ENDDO
218 ENDIF
219 ENDDO
220 #ifdef NONLIN_FRSURF
221 C this is a bit of a hack
222 IF ( nonlinFreeSurf.GT.0 ) THEN
223 DO J=1-OLy,sNy+OLy
224 I = OB_Ie(J,bi,bj)
225 IF ( I.NE.OB_indexNone ) THEN
226 OBEeta(J,bi,bj) = etaN(I-1,J,bi,bj)
227 ENDIF
228 ENDDO
229 ENDIF
230 #endif /* NONLIN_FRSURF */
231 C Next, we compute the phase speed correction, which depends on the
232 C tracer!
233 DO K=1,Nr
234 DO J=1-OLy,sNy+OLy
235 I = OB_Ie(J,bi,bj)
236 IF ( I.NE.OB_indexNone ) THEN
237 cflMer(J,K) = 0.5 _d 0 * _dxC(I-1,J,bi,bj)/dTtracerLev(K)
238 CML gFacM(J,K) = 0. _d 0
239 CML IF ( uVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacM(J,K) = 1. _d 0
240 gFacM(J,K) = ABS(MIN(SIGN(1.D0,uVel(I,J,K,bi,bj)),0.D0))
241 ELSE
242 cflMer(J,K) = 0. _d 0
243 gFacM (J,K) = 0. _d 0
244 ENDIF
245 ENDDO
246 ENDDO
247 DO K=1,Nr
248 DO J=1-OLy,sNy+OLy
249 I = OB_Ie(J,bi,bj)
250 IF ( I.NE.OB_indexNone ) THEN
251 C Theta first:
252 dTracSpace = (theta(I-1,J,K,bi,bj)-theta(I-2,J,K,bi,bj))
253 & * _maskW(I-1,J,K,bi,bj)
254 dTracTime = (theta(I-1,J,K,bi,bj)-OBEtStevens(J,K,bi,bj))
255 uPhase = cflMer(J,K) * pFac
256 IF ( dTracSpace .NE. 0. _d 0 ) THEN
257 uPhase = MIN( cflMer(J,K),
258 & MAX( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
259 & ) * pFac
260 ENDIF
261 C Compute the tracer tendency here, the tracer will be updated
262 C with a simple Euler forward step in S/R obcs_apply_ts
263 OBEt(J,K,bi,bj) = _maskW(I,J,K,bi,bj) * (
264 & - ( aFac*MAX(0.D0,uVel(I,J,K,bi,bj)) + uPhase )
265 & *(theta(I,J,K,bi,bj)-theta(I-1,J,K,bi,bj))
266 & * _recip_dxC(I,J,bi,bj)
267 & - gFacM(J,K) * gammat
268 & * (theta(I,J,K,bi,bj)-OBEt(J,K,bi,bj)) )
269 C Now salinity:
270 dTracSpace = (salt(I-1,J,K,bi,bj)-salt(I-2,J,K,bi,bj))
271 & * _maskW(I-1,J,K,bi,bj)
272 dTracTime = (salt(I-1,J,K,bi,bj)-OBEsStevens(J,K,bi,bj))
273 uPhase = cflMer(J,K) * pFac
274 IF ( dTracSpace .NE. 0. _d 0 ) THEN
275 uPhase = MIN( cflMer(J,K),
276 & MAX( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
277 & ) * pFac
278 ENDIF
279 C Compute the tracer tendency here, the tracer will be updated
280 C with a simple Euler forward step in S/R obcs_apply_ts
281 OBEs(J,K,bi,bj) = _maskW(I,J,K,bi,bj) * (
282 & - ( aFac*MAX(0.D0,uVel(I,J,K,bi,bj)) + uPhase )
283 & *(salt(I,J,K,bi,bj)-salt(I-1,J,K,bi,bj))
284 & * _recip_dxC(I,J,bi,bj)
285 & - gFacM(J,K) * gammas
286 & * (salt(I,J,K,bi,bj)-OBEs(J,K,bi,bj)) )
287 ENDIF
288 ENDDO
289 ENDDO
290 C IF ( useStevensEast ) THEN
291 ENDIF
292 #endif /* ALLOW_OBCS_EAST */
293
294 C ------------------------------------------------------------------------------
295
296 #ifdef ALLOW_OBCS_WEST
297
298 # ifdef ALLOW_AUTODIFF_TAMC
299 CADJ STORE OBWt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
300 CADJ STORE OBWs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
301 # endif
302 IF ( useStevensWest ) THEN
303 C Western OB
304 #ifdef ALLOW_DEBUG
305 IF (debugMode)
306 & CALL DEBUG_MSG('OBCS_CALC_STEVENS: West',myThid)
307 #endif
308 C compute vertical average and deviation from vertical
309 C average for I_obw+1
310 DO J=1-OLy,sNy+OLy
311 I = OB_Iw(J,bi,bj)
312 IF ( I.NE.OB_indexNone ) THEN
313 C first initialize some fields
314 drFBar = 0. _d 0
315 uMerBar = 0. _d 0
316 DO K=1,Nr
317 uMerPri(K) = 0. _d 0
318 ENDDO
319 DO K=1,Nr
320 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
321 uMerBar = uMerBar + uVel(I+2,J,K,bi,bj)
322 #else
323 uMerBar = uMerBar + OBWuStevens(J,K,bi,bj)
324 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
325 & *drF(K)* _hFacW(I+1,J,K,bi,bj)
326 drFBar = drFBar + drF(K)* _hFacW(I+1,J,K,bi,bj)
327 ENDDO
328 IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
329 DO K=1,Nr
330 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
331 uMerPri(K) = (uVel(I+2,J,K,bi,bj)-uMerBar)
332 #else
333 uMerPri(K) = (OBWuStevens(J,K,bi,bj)-uMerBar)
334 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
335 & * _maskW(I+1,J,K,bi,bj)
336 ENDDO
337 C vertical average of input field
338 drFBar = 0. _d 0
339 uMerBar = 0. _d 0
340 DO K=1,Nr
341 uMerBar = uMerBar + OBWu(J,K,bi,bj)
342 & *drF(K)* _hFacW(I+1,J,K,bi,bj)
343 drFBar = drFBar + drF(K)* _hFacW(I+1,J,K,bi,bj)
344 ENDDO
345 IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
346 C Now the absolute velocity normal to the boundary is
347 C uMerBar + uMerPri(K).
348 DO K=1,Nr
349 OBWu(J,K,bi,bj) = (uMerBar + uMerPri(K))
350 & * _maskW(I+1,J,K,bi,bj)
351 CML OBWv(J,K,bi,bj) = 0. _d 0
352 #ifdef ALLOW_NONHYDROSTATIC
353 OBWw(J,K,bi,bj)=0.
354 #endif
355 ENDDO
356 ENDIF
357 ENDDO
358 #ifdef NONLIN_FRSURF
359 C this is a bit of a hack
360 IF ( nonlinFreeSurf.GT.0 ) THEN
361 DO J=1-OLy,sNy+OLy
362 I = OB_Iw(J,bi,bj)
363 IF ( I.NE.OB_indexNone ) THEN
364 OBWeta(J,bi,bj) = etaN(I+1,J,bi,bj)
365 ENDIF
366 ENDDO
367 ENDIF
368 #endif /* NONLIN_FRSURF */
369 C Next, we compute the phase speed correction, which depends on the
370 C tracer!
371 DO K=1,Nr
372 DO J=1-OLy,sNy+OLy
373 I = OB_Iw(J,bi,bj)
374 IF ( I.NE.OB_indexNone ) THEN
375 cflMer(J,K) = 0.5 _d 0 * _dxC(I+2,J,bi,bj)/dTtracerLev(K)
376 CML gFacM = 0. _d 0
377 CML IF ( uVel(I+1,J,K,bi,bj) .GT. 0. _d 0 ) gFacM = 1. _d 0
378 gFacM(J,K) = ABS(MAX(SIGN(1.D0,uVel(I+1,J,K,bi,bj)),0.D0))
379 ELSE
380 cflMer(J,K) = 0. _d 0
381 gFacM (J,K) = 0. _d 0
382 ENDIF
383 ENDDO
384 ENDDO
385 DO K=1,Nr
386 DO J=1-OLy,sNy+OLy
387 I = OB_Iw(J,bi,bj)
388 IF ( I.NE.OB_indexNone ) THEN
389 C Theta first:
390 dTracSpace = (theta(I+2,J,K,bi,bj)-theta(I+1,J,K,bi,bj))
391 & * _maskW(I+2,J,K,bi,bj)
392 dTracTime = (theta(I+1,J,K,bi,bj)-OBWtStevens(J,K,bi,bj))
393 uPhase = -cflMer(J,K) * pFac
394 IF ( dTracSpace .NE. 0. _d 0 ) THEN
395 uPhase = MAX( -cflMer(J,K),
396 & MIN( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
397 & ) * pFac
398 ENDIF
399 C Compute the tracer tendency here, the tracer will be updated
400 C with a simple Euler forward step in S/R obcs_apply_ts
401 OBWt(J,K,bi,bj) = _maskW(I+1,J,K,bi,bj) * (
402 & - ( aFac*MIN(0.D0,uVel(I+1,J,K,bi,bj)) + uPhase )
403 & *(theta(I+1,J,K,bi,bj)-theta(I,J,K,bi,bj))
404 & * _recip_dxC(I+1,J,bi,bj)
405 & - gFacM(J,K) * gammat
406 & * (theta(I,J,K,bi,bj)-OBWt(J,K,bi,bj)) )
407 C Now salinity:
408 dTracSpace = (salt(I+2,J,K,bi,bj)-salt(I+1,J,K,bi,bj))
409 & * _maskW(I+2,J,K,bi,bj)
410 dTracTime = (salt(I+1,J,K,bi,bj)-OBWsStevens(J,K,bi,bj))
411 uPhase = -cflMer(J,K) * pFac
412 IF ( dTracSpace .NE. 0. _d 0 ) THEN
413 uPhase = MAX( -cflMer(J,K),
414 & MIN( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
415 & ) * pFac
416 ENDIF
417 C Compute the tracer tendency here, the tracer will be updated
418 C with a simple Euler forward step in S/R obcs_apply_ts
419 OBWs(J,K,bi,bj) = _maskW(I+1,J,K,bi,bj) * (
420 & - ( aFac*MIN(0.D0,uVel(I+1,J,K,bi,bj)) + uPhase )
421 & *(salt(I+1,J,K,bi,bj)-salt(I,J,K,bi,bj))
422 & * _recip_dxC(I+1,J,bi,bj)
423 & - gFacM(J,K) * gammas
424 & * (salt(I,J,K,bi,bj)-OBWs(J,K,bi,bj)) )
425 ENDIF
426 ENDDO
427 ENDDO
428 C IF ( useStevensWest ) THEN
429 ENDIF
430 #endif /* ALLOW_OBCS_WEST */
431
432 C ------------------------------------------------------------------------------
433
434 #ifdef ALLOW_OBCS_NORTH
435 # ifdef ALLOW_AUTODIFF_TAMC
436 CADJ STORE OBNt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
437 CADJ STORE OBNs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
438 # endif
439 IF ( useStevensNorth ) THEN
440 C Northern OB
441 #ifdef ALLOW_DEBUG
442 IF (debugMode)
443 & CALL DEBUG_MSG('OBCS_CALC_STEVENS: North',myThid)
444 #endif
445 C compute vertical average and deviation from vertical
446 C average for J_obn
447 DO I=1-OLx,sNx+OLx
448 J = OB_Jn(I,bi,bj)
449 IF ( J.NE.OB_indexNone ) THEN
450 C first initialize some fields
451 drFBar = 0. _d 0
452 vZonBar = 0. _d 0
453 DO K=1,Nr
454 vZonPri(K) = 0. _d 0
455 ENDDO
456 DO K=1,Nr
457 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
458 vZonBar = vZonBar + vVel(I,J-1,K,bi,bj)
459 #else
460 vZonBar = vZonBar + OBNvStevens(I,K,bi,bj)
461 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
462 & *drF(K)* _hFacS(I,J,K,bi,bj)
463 drFBar = drFBar + drF(K)* _hFacS(I,J,K,bi,bj)
464 ENDDO
465 IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
466 DO K=1,Nr
467 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
468 vZonPri(K) = (vVel(I,J-1,K,bi,bj)-vZonBar)
469 #else
470 vZonPri(K) = (OBNvStevens(I,K,bi,bj)-vZonBar)
471 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
472 & * _maskS(I,J,K,bi,bj)
473 ENDDO
474 C vertical average of input field
475 drFBar = 0. _d 0
476 vZonBar = 0. _d 0
477 DO K=1,Nr
478 vZonBar = vZonBar + OBNv(I,K,bi,bj)
479 & *drF(K)* _hFacS(I,J,K,bi,bj)
480 drFBar = drFBar + drF(K)* _hFacS(I,J,K,bi,bj)
481 ENDDO
482 IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
483 C Now the absolute velocity normal to the boundary is
484 C vZonBar + vZonPri(K).
485 DO K=1,Nr
486 OBNv(I,K,bi,bj) = (vZonBar + vZonPri(K))
487 & * _maskS(I,J,K,bi,bj)
488 CML OBNu(I,K,bi,bj) = 0. _d 0
489 #ifdef ALLOW_NONHYDROSTATIC
490 OBNw(I,K,bi,bj)=0.
491 #endif
492 ENDDO
493 ENDIF
494 ENDDO
495 #ifdef NONLIN_FRSURF
496 C this is a bit of a hack
497 IF ( nonlinFreeSurf.GT.0 ) THEN
498 DO I=1-OLx,sNx+OLx
499 J = OB_Jn(I,bi,bj)
500 IF ( J.NE.OB_indexNone ) THEN
501 OBNeta(I,bi,bj) = etaN(I,J-1,bi,bj)
502 ENDIF
503 ENDDO
504 ENDIF
505 #endif /* NONLIN_FRSURF */
506 C Next, we compute the phase speed correction, which depends on the
507 C tracer!
508 DO K=1,Nr
509 DO I=1-OLx,sNx+OLx
510 J = OB_Jn(I,bi,bj)
511 IF ( J.NE.OB_indexNone ) THEN
512 cflZon(I,K) = 0.5 _d 0 * _dyC(I,J-1,bi,bj)/dTtracerLev(K)
513 CML gFacZ(I,K) = 0. _d 0
514 CML IF ( vVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacZ(I,K) = 1. _d 0
515 gFacZ(I,K) = ABS(MIN(SIGN(1.D0,vVel(I,J,K,bi,bj)),0.D0))
516 ELSE
517 cflZon(I,K) = 0. _d 0
518 gFacZ (I,K) = 0. _d 0
519 ENDIF
520 ENDDO
521 ENDDO
522 DO K=1,Nr
523 DO I=1-OLx,sNx+OLx
524 J = OB_Jn(I,bi,bj)
525 IF ( J.NE.OB_indexNone ) THEN
526 C Theta first:
527 dTracSpace = (theta(I,J-1,K,bi,bj)-theta(I,J-2,K,bi,bj))
528 & * _maskS(I,J-1,K,bi,bj)
529 dTracTime = (theta(I,J-1,K,bi,bj)-OBNtStevens(I,K,bi,bj))
530 vPhase = cflZon(I,K) * pFac
531 IF ( dTracSpace .NE. 0. _d 0 ) THEN
532 vPhase = MIN( cflZon(I,K),
533 & MAX( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
534 & ) * pFac
535 ENDIF
536 C Compute the tracer tendency here, the tracer will be updated
537 C with a simple Euler forward step in S/R obcs_apply_ts
538 OBNt(I,K,bi,bj) = _maskS(I,J,K,bi,bj) * (
539 & - ( aFac*MAX(0.D0,vVel(I,J,K,bi,bj)) + vPhase )
540 & *(theta(I,J,K,bi,bj)-theta(I,J-1,K,bi,bj))
541 & * _recip_dyC(I,J,bi,bj)
542 & - gFacZ(I,K) * gammat
543 & * (theta(I,J,K,bi,bj)-OBNt(I,K,bi,bj)) )
544 C Now salinity:
545 dTracSpace = (salt(I,J-1,K,bi,bj)-salt(I,J-2,K,bi,bj))
546 & * _maskS(I,J-1,K,bi,bj)
547 dTracTime = (salt(I,J-1,K,bi,bj)-OBNsStevens(I,K,bi,bj))
548 vPhase = cflZon(I,K) * pFac
549 IF ( dTracSpace .NE. 0. _d 0 ) THEN
550 vPhase = MIN( cflZon(I,K),
551 & MAX( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
552 & ) * pFac
553 ENDIF
554 C Compute the tracer tendency here, the tracer will be updated
555 C with a simple Euler forward step in S/R obcs_apply_ts
556 OBNs(I,K,bi,bj) = _maskS(I,J,K,bi,bj) * (
557 & - ( aFac*MAX(0.D0,vVel(I,J,K,bi,bj)) + vPhase )
558 & *(salt(I,J,K,bi,bj)-salt(I,J-1,K,bi,bj))
559 & * _recip_dyC(I,J,bi,bj)
560 & - gFacZ(I,K) * gammas
561 & * (salt(I,J,K,bi,bj)-OBNs(I,K,bi,bj)) )
562 ENDIF
563 ENDDO
564 ENDDO
565 C IF ( useStevensNorth ) THEN
566 ENDIF
567 #endif /* ALLOW_OBCS_NORTH */
568
569 C ------------------------------------------------------------------------------
570
571 #ifdef ALLOW_OBCS_SOUTH
572
573 # ifdef ALLOW_AUTODIFF_TAMC
574 CADJ STORE OBSt(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
575 CADJ STORE OBSs(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
576 # endif
577 IF ( useStevensSouth ) THEN
578 C Southern OB
579 #ifdef ALLOW_DEBUG
580 IF (debugMode)
581 & CALL DEBUG_MSG('OBCS_CALC_STEVENS: South',myThid)
582 #endif
583 C compute vertical average and deviation from vertical
584 C average for J_obs+1
585 DO I=1-OLx,sNx+OLx
586 J = OB_Js(I,bi,bj)
587 IF ( J.NE.OB_indexNone ) THEN
588 C first initialize some fields
589 drFBar = 0. _d 0
590 vZonBar = 0. _d 0
591 DO K=1,Nr
592 vZonPri(K) = 0. _d 0
593 ENDDO
594 DO K=1,Nr
595 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
596 vZonBar = vZonBar + vVel(I,J+2,K,bi,bj)
597 #else
598 vZonBar = vZonBar + OBSvStevens(I,K,bi,bj)
599 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
600 & *drF(K)* _hFacS(I,J+1,K,bi,bj)
601 drFBar = drFBar + drF(K)* _hFacS(I,J+1,K,bi,bj)
602 ENDDO
603 IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
604 DO K=1,Nr
605 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
606 vZonPri(K) = (vVel(I,J+2,K,bi,bj)-vZonBar)
607 #else
608 vZonPri(K) = (OBSvStevens(I,K,bi,bj)-vZonBar)
609 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
610 & * _maskS(I,J+1,K,bi,bj)
611 ENDDO
612 C vertical average of input field
613 drFBar = 0. _d 0
614 vZonBar = 0. _d 0
615 DO K=1,Nr
616 vZonBar = vZonBar + OBSv(I,K,bi,bj)
617 & *drF(K)* _hFacS(I,J+1,K,bi,bj)
618 drFBar = drFBar + drF(K)* _hFacS(I,J+1,K,bi,bj)
619 ENDDO
620 IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
621 C Now the absolute velocity normal to the boundary is
622 C vZonBar + vZonPri(K).
623 DO K=1,Nr
624 OBSv(I,K,bi,bj) = (vZonBar + vZonPri(K))
625 & * _maskS(I,J+1,K,bi,bj)
626 CML OBSu(I,K,bi,bj) = 0. _d 0
627 #ifdef ALLOW_NONHYDROSTATIC
628 OBSw(I,K,bi,bj)=0.
629 #endif
630 ENDDO
631 ENDIF
632 ENDDO
633 #ifdef NONLIN_FRSURF
634 C this is a bit of a hack
635 IF ( nonlinFreeSurf.GT.0 ) THEN
636 DO I=1-OLx,sNx+OLx
637 J = OB_Js(I,bi,bj)
638 IF ( J.NE.OB_indexNone ) THEN
639 OBSeta(I,bi,bj) = etaN(I,J+1,bi,bj)
640 ENDIF
641 ENDDO
642 ENDIF
643 #endif /* NONLIN_FRSURF */
644 C Next, we compute the phase speed correction, which depends on the
645 C tracer!
646 DO K=1,Nr
647 DO I=1-OLx,sNx+OLx
648 J = OB_Js(I,bi,bj)
649 IF ( J.NE.OB_indexNone ) THEN
650 cflZon(I,K) = 0.5 _d 0 * _dyC(I,J+2,bi,bj)/dTtracerLev(K)
651 CML gFacZ = 0. _d 0
652 CML IF ( vVel(I,J+1,K,bi,bj) .GT. 0. _d 0 ) gFacZ = 1. _d 0
653 gFacZ(I,K) = ABS(MAX(SIGN(1.D0,vVel(I,J+1,K,bi,bj)),0.D0))
654 ELSE
655 cflZon(I,K) = 0. _d 0
656 gFacZ (I,K) = 0. _d 0
657 ENDIF
658 ENDDO
659 ENDDO
660 DO K=1,Nr
661 DO I=1-OLx,sNx+OLx
662 J = OB_Js(I,bi,bj)
663 IF ( J.NE.OB_indexNone ) THEN
664 C Theta first:
665 dTracSpace = (theta(I,J+2,K,bi,bj)-theta(I,J+1,K,bi,bj))
666 & * _maskS(I,J+2,K,bi,bj)
667 dTracTime = (theta(I,J+1,K,bi,bj)-OBStStevens(I,K,bi,bj))
668 vPhase = -cflZon(I,K) * pFac
669 IF ( dTracSpace .NE. 0. _d 0 ) THEN
670 vPhase = MAX( -cflZon(I,K),
671 & MIN( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
672 & ) * pFac
673 ENDIF
674 C Compute the tracer tendency here, the tracer will be updated
675 C with a simple Euler forward step in S/R obcs_apply_ts
676 OBSt(I,K,bi,bj) = _maskS(I,J+1,K,bi,bj) * (
677 & - ( aFac*MIN(0.D0,vVel(I,J+1,K,bi,bj)) + vPhase )
678 & *(theta(I,J+1,K,bi,bj)-theta(I,J,K,bi,bj))
679 & * _recip_dyC(I,J+1,bi,bj)
680 & - gFacZ(I,K) * gammat
681 & * (theta(I,J,K,bi,bj)-OBSt(I,K,bi,bj)) )
682 C Now salinity:
683 dTracSpace = (salt(I,J+2,K,bi,bj)-salt(I,J+1,K,bi,bj))
684 & * _maskS(I,J+2,K,bi,bj)
685 dTracTime = (salt(I,J+1,K,bi,bj)-OBSsStevens(I,K,bi,bj))
686 vPhase = -cflZon(I,K) * pFac
687 IF ( dTracSpace .NE. 0. _d 0 ) THEN
688 vPhase = MAX( -cflZon(I,K),
689 & MIN( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
690 & ) * pFac
691 ENDIF
692 C Compute the tracer tendency here, the tracer will be updated
693 C with a simple Euler forward step in S/R obcs_apply_ts
694 OBSs(I,K,bi,bj) = _maskS(I,J+1,K,bi,bj) * (
695 & - ( aFac*MIN(0.D0,vVel(I,J+1,K,bi,bj)) + vPhase )
696 & *(salt(I,J+1,K,bi,bj)-salt(I,J,K,bi,bj))
697 & * _recip_dyC(I,J+1,bi,bj)
698 & - gFacZ(I,K) * gammas
699 & * (salt(I,J,K,bi,bj)-OBSs(I,K,bi,bj)) )
700 ENDIF
701 ENDDO
702 ENDDO
703 C IF ( useStevensSouth ) THEN
704 ENDIF
705 #endif /* ALLOW_OBCS_SOUTH */
706
707 C end bi/bj-loops
708 ENDDO
709 ENDDO
710
711 C save the tracer fields of the previous time step for the next time step
712 CALL OBCS_STEVENS_SAVE_TRACERS(
713 I futureTime, futureIter,
714 I myThid )
715 C ------------------------------------------------------------------------------
716
717 #ifdef CHECK_BALANCE
718 DO bj=myByLo(myThid),myByHi(myThid)
719 DO bi=myBxLo(myThid),myBxHi(myThid)
720 uPhase=0.
721 vPhase=0.
722 uVelLoc = 0.
723 DO J=1-OLy,sNy+OLy
724 uMerBar=0. _d 0
725 DO K=1,Nr
726 I = OB_Ie(J,bi,bj)
727 IF ( I.EQ.OB_indexNone ) I = 1
728 uPhase = uPhase + OBEu(J,K,bi,bj)
729 & *drF(k)* _hFacW(I,J,K,bi,bj)*dyG(I,J,bi,bj)
730 I = OB_Iw(J,bi,bj)
731 IF ( I.EQ.OB_indexNone ) I = 1
732 vPhase = vPhase + OBWu(J,K,bi,bj)
733 & *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
734 CML uVelLoc = uVelLoc + uMerPri(J,K)
735 CML & *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
736 CML uMerBar(J)=uMerBar(J) + uMerPri(J,K)
737 CML & *drF(k)* _hFacW(I+1,J,K,bi,bj)
738 ENDDO
739 CML print *, 'ml-obcs: uBar = ', j,uMerBar(J)
740 ENDDO
741 C end bi/bj-loops
742 ENDDO
743 ENDDO
744 _GLOBAL_SUM_RL( uPhase, myThid )
745 _GLOBAL_SUM_RL( vPhase, myThid )
746 CML _GLOBAL_SUM_RL( uVelLoc, myThid )
747 print *, 'ml-obcs: OBE = ', uPhase*1 _d -6, ' Sv'
748 print *, 'ml-obcs: OBW = ', vPhase*1 _d -6, ' Sv'
749 CML print *, 'ml-obcs: OBWp = ', uVelLoc*1 _d -6, ' Sv'
750 #endif /* CHECK_BALANCE */
751
752 #ifdef ALLOW_DEBUG
753 IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC_STEVENS',myThid)
754 #endif
755
756 #endif /* ALLOW_OBCS_STEVENS */
757 RETURN
758 END
759
760 CBOP
761 C !ROUTINE: OBCS_STEVENS_SAVE_TRACERS
762 C !INTERFACE:
763 SUBROUTINE OBCS_STEVENS_SAVE_TRACERS(
764 I futureTime, futureIter,
765 I myThid )
766 C !DESCRIPTION: \bv
767 C *==========================================================*
768 C | SUBROUTINE OBCS_STEVENS_SAVE_TRACERS
769 C | Save tracers (of previous time step) at the OB location
770 C | to be used in the next time step for Stevens boundary
771 C | conditions
772 C *==========================================================*
773 C \ev
774 C !USES:
775 IMPLICIT NONE
776 C == Global variables ==
777 #include "SIZE.h"
778 #include "EEPARAMS.h"
779 #include "GRID.h"
780 #include "OBCS_PARAMS.h"
781 #include "OBCS_GRID.h"
782 #include "OBCS_FIELDS.h"
783 #include "DYNVARS.h"
784 CML#ifdef ALLOW_PTRACERS
785 CML#include "PTRACERS_SIZE.h"
786 CML#include "PTRACERS_PARAMS.h"
787 CML#include "PTRACERS_FIELDS.h"
788 CML#include "OBCS_PTRACERS.h"
789 CML#endif /* ALLOW_PTRACERS */
790
791 C !INPUT/OUTPUT PARAMETERS:
792 C == Routine arguments ==
793 C myThid :: my Thread Id number
794 _RL futureTime
795 INTEGER futureIter
796 INTEGER myThid
797
798 #ifdef ALLOW_OBCS_STEVENS
799 C !LOCAL VARIABLES:
800 C bi, bj :: indices of current tile
801 C i,j,k :: loop indices
802 C Iobc,Jobc :: position-index of open boundary
803 INTEGER bi,bj,i,j,k,Iobc,Jobc
804 CEOP
805
806 DO bj=myByLo(myThid),myByHi(myThid)
807 DO bi=myBxLo(myThid),myBxHi(myThid)
808 #ifdef ALLOW_OBCS_NORTH
809 IF ( tileHasOBN(bi,bj) .AND. useStevensNorth ) THEN
810 C Northern boundary
811 DO i=1-OLx,sNx+OLx
812 Jobc = OB_Jn(i,bi,bj)
813 IF ( Jobc.NE.OB_indexNone ) THEN
814 DO k = 1,Nr
815 OBNtStevens(i,k,bi,bj) = theta(i,Jobc-1,k,bi,bj)
816 & *maskC(i,Jobc+1,k,bi,bj)
817 OBNsStevens(i,k,bi,bj) = salt(i,Jobc-1,k,bi,bj)
818 & *maskC(i,Jobc+1,k,bi,bj)
819 ENDDO
820 ENDIF
821 ENDDO
822 ENDIF
823 #endif /* ALLOW_OBCS_NORTH */
824 #ifdef ALLOW_OBCS_SOUTH
825 IF ( tileHasOBS(bi,bj) .AND. useStevensSouth ) THEN
826 C Southern boundary
827 DO i=1-OLx,sNx+OLx
828 Jobc = OB_Js(i,bi,bj)
829 IF ( Jobc.NE.OB_indexNone ) THEN
830 DO k = 1,Nr
831 OBStStevens(i,k,bi,bj) = theta(i,Jobc+1,k,bi,bj)
832 & *maskC(i,Jobc+1,k,bi,bj)
833 OBSsStevens(i,k,bi,bj) = salt(i,Jobc+1,k,bi,bj)
834 & *maskC(i,Jobc+1,k,bi,bj)
835 ENDDO
836 ENDIF
837 ENDDO
838 ENDIF
839 #endif /* ALLOW_OBCS_SOUTH */
840 #ifdef ALLOW_OBCS_EAST
841 IF ( tileHasOBE(bi,bj) .AND. useStevensEast ) THEN
842 C Eastern boundary
843 DO j=1-OLy,sNy+OLy
844 Iobc = OB_Ie(j,bi,bj)
845 IF ( Iobc.NE.OB_indexNone ) THEN
846 DO k = 1,Nr
847 OBEtStevens(j,k,bi,bj) = theta(Iobc-1,j,k,bi,bj)
848 & *maskC(Iobc-1,j,k,bi,bj)
849 OBEsStevens(j,k,bi,bj) = salt(Iobc-1,j,k,bi,bj)
850 & *maskC(Iobc-1,j,k,bi,bj)
851 ENDDO
852 ENDIF
853 ENDDO
854 ENDIF
855 #endif /* ALLOW_OBCS_EAST */
856 #ifdef ALLOW_OBCS_WEST
857 IF ( tileHasOBW(bi,bj) .AND. useStevensWest ) THEN
858 C Western boundary
859 DO j=1-OLy,sNy+OLy
860 Iobc = OB_Iw(j,bi,bj)
861 IF ( Iobc.NE.OB_indexNone ) THEN
862 DO k = 1,Nr
863 OBWtStevens(j,k,bi,bj) = theta(Iobc+1,j,k,bi,bj)
864 & *maskC(Iobc+1,j,k,bi,bj)
865 OBWsStevens(j,k,bi,bj) = salt(Iobc+1,j,k,bi,bj)
866 & *maskC(Iobc+1,j,k,bi,bj)
867 ENDDO
868 ENDIF
869 ENDDO
870 ENDIF
871 #endif /* ALLOW_OBCS_WEST */
872 ENDDO
873 ENDDO
874 #endif /* ALLOW_OBCS_STEVENS */
875 RETURN
876 END
  ViewVC Help
Powered by ViewVC 1.1.22