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

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

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

Revision 1.7 - (show annotations) (download)
Thu Dec 15 01:12:58 2005 UTC (18 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint58e_post, checkpoint58u_post, checkpoint58w_post, checkpoint60, checkpoint61, checkpoint62, checkpoint58r_post, checkpoint58n_post, checkpoint58x_post, checkpoint58t_post, checkpoint58h_post, checkpoint58q_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint58, checkpoint58f_post, checkpoint58d_post, checkpoint58c_post, checkpoint58a_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62d, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint58p_post, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y, checkpoint58b_post, checkpoint58m_post
Changes since 1.6: +10 -9 lines 
remove unused variables (reduce number of warnings).
1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_sponge.F,v 1.6 2005/04/27 14:10:06 jmc Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5 #include "AD_CONFIG.h"
6
7 CStartOfInterface
8 SUBROUTINE OBCS_SPONGE_U(
9 I iMin, iMax, jMin, jMax,bi,bj,kLev,
10 I myCurrentTime,myThid)
11 C /==========================================================\
12 C | S/R OBCS_SPONGE_U |
13 C | o Contains problem specific forcing for zonal velocity. |
14 C |==========================================================|
15 C | Adds terms to gU for forcing by external sources |
16 C | e.g. wind stress, bottom friction etc.................. |
17 C \==========================================================/
18 IMPLICIT NONE
19
20 C == Global data ==
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "GRID.h"
25 #include "DYNVARS.h"
26 #include "FFIELDS.h"
27 #ifdef ALLOW_OBCS
28 # include "OBCS.h"
29 # ifdef ALLOW_CAL
30 # include "cal.h"
31 # endif
32 #endif
33
34 C == Routine arguments ==
35 C iMin - Working range of tile for applying forcing.
36 C iMax
37 C jMin
38 C jMax
39 C kLev
40 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
41 _RL myCurrentTime
42 INTEGER myThid
43 CEndOfInterface
44
45 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
46 C == Local variables ==
47 C Loop counters
48 INTEGER I, J, Isl, Jsl
49 _RL urelax, lambda_obcs_u
50 #ifndef ALLOW_CAL
51 INTEGER secondsperday
52 PARAMETER (secondsperday=86400)
53 #endif
54
55
56 IF (useOBCSsponge) THEN
57
58 C Northern Open Boundary
59 #ifdef ALLOW_OBCS_NORTH
60 DO i=iMin,iMax
61 IF ((OB_Jn(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
62 DO jsl= 1,spongeThickness
63 j=OB_Jn(i,bi,bj)-jsl
64
65 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
66 urelax=(
67 & float(spongeThickness-jsl)*OBNu(i,kLev,bi,bj)
68 & + float(jsl)*uVel(i,j,kLev,bi,bj) )
69 & / float(spongeThickness)
70
71 lambda_obcs_u = (
72 & float(spongeThickness-jsl)*Vrelaxobcsbound
73 & + float(jsl)*Vrelaxobcsinner)
74 & / float(spongeThickness)
75
76 IF (lambda_obcs_u.ne.0.) THEN
77 lambda_obcs_u = 1. _d 0 / lambda_obcs_u
78 ELSE
79 lambda_obcs_u = 0. _d 0
80 ENDIF
81
82 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
83 & - _maskW(i,j,kLev,bi,bj) * lambda_obcs_u
84 & * ( uVel(i,j,kLev,bi,bj) - urelax )
85 ENDIF
86
87 ENDDO
88 ENDIF
89 ENDDO
90 #endif
91
92 C Southern Open Boundary
93 #ifdef ALLOW_OBCS_SOUTH
94 DO i=iMin,iMax
95 IF ((OB_Js(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
96 DO jsl= 1,spongeThickness
97 j=OB_Js(i,bi,bj)+jsl
98
99 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
100 urelax=(
101 & float(spongeThickness-jsl)*OBSu(i,kLev,bi,bj)
102 & + float(jsl)*uVel(i,j,kLev,bi,bj) )
103 & / float(spongeThickness)
104
105 lambda_obcs_u = (
106 & float(spongeThickness-jsl)*Vrelaxobcsbound
107 & + float(jsl)*Vrelaxobcsinner)
108 & / float(spongeThickness)
109
110 if (lambda_obcs_u.ne.0.) then
111 lambda_obcs_u = 1. _d 0 / lambda_obcs_u
112 else
113 lambda_obcs_u = 0. _d 0
114 endif
115
116 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
117 & - _maskW(i,j,kLev,bi,bj) * lambda_obcs_u
118 & * ( uVel(i,j,kLev,bi,bj) - urelax )
119 ENDIF
120
121 ENDDO
122 ENDIF
123 ENDDO
124 #endif
125
126 C Eastern Open Boundary
127 #ifdef ALLOW_OBCS_EAST
128 DO j=jMin,jMax
129 IF ((OB_Ie(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
130 DO isl= 1,spongeThickness
131 i=OB_Ie(j,bi,bj)-isl
132
133 IF ((i.ge.imin).and.(i.le.imax)) THEN
134 urelax=(
135 & float(spongeThickness-isl)*OBEu(j,kLev,bi,bj)
136 & + float(isl)*uVel(i,j,kLev,bi,bj) )
137 & / float(spongeThickness)
138
139 lambda_obcs_u = (
140 & float(spongeThickness-isl)*Urelaxobcsbound
141 & + float(isl)*Urelaxobcsinner)
142 & / float(spongeThickness)
143
144 if (lambda_obcs_u.ne.0.) then
145 lambda_obcs_u = 1. _d 0 / lambda_obcs_u
146 else
147 lambda_obcs_u = 0. _d 0
148 endif
149
150 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
151 & - _maskW(i,j,kLev,bi,bj) * lambda_obcs_u
152 & * ( uVel(i,j,kLev,bi,bj) - urelax )
153 ENDIF
154
155 ENDDO
156 ENDIF
157 ENDDO
158 #endif
159
160 C Western Open Boundary
161 #ifdef ALLOW_OBCS_WEST
162 DO j=jMin,jMax
163 IF ((OB_Iw(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
164 DO isl= 1,spongeThickness
165 i=OB_Iw(j,bi,bj)+isl+1
166
167 IF ((i.ge.imin).and.(i.le.imax)) THEN
168 urelax=(
169 & float(spongeThickness-isl)*OBWu(j,kLev,bi,bj)
170 & + float(isl)*uVel(i,j,kLev,bi,bj) )
171 & / float(spongeThickness)
172
173 lambda_obcs_u= (
174 & float(spongeThickness-isl)*Urelaxobcsbound
175 & + float(isl)*Urelaxobcsinner)
176 & / float(spongeThickness)
177
178 if (lambda_obcs_u.ne.0.) then
179 lambda_obcs_u = 1. _d 0 / lambda_obcs_u
180 else
181 lambda_obcs_u = 0. _d 0
182 endif
183
184 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
185 & - _maskW(i,j,kLev,bi,bj) * lambda_obcs_u
186 & * ( uVel(i,j,kLev,bi,bj) - urelax )
187 ENDIF
188
189 ENDDO
190 ENDIF
191 ENDDO
192 #endif
193
194 ENDIF
195
196 #endif /* ALLOW_OBCS & ALLOW_OBCS_SPONGE */
197
198 RETURN
199 END
200
201 CStartOfInterface
202 SUBROUTINE OBCS_SPONGE_V(
203 I iMin, iMax, jMin, jMax,bi,bj,kLev,
204 I myCurrentTime,myThid)
205 C /==========================================================\
206 C | S/R OBCS_SPONGE_V |
207 C | o Contains problem specific forcing for merid velocity. |
208 C |==========================================================|
209 C | Adds terms to gV for forcing by external sources |
210 C | e.g. wind stress, bottom friction etc.................. |
211 C \==========================================================/
212 IMPLICIT NONE
213
214 C == Global data ==
215 #include "SIZE.h"
216 #include "EEPARAMS.h"
217 #include "PARAMS.h"
218 #include "GRID.h"
219 #include "DYNVARS.h"
220 #include "FFIELDS.h"
221 #ifdef ALLOW_OBCS
222 # include "OBCS.h"
223 # ifdef ALLOW_CAL
224 # include "cal.h"
225 # endif
226 #endif
227
228 C == Routine arguments ==
229 C iMin - Working range of tile for applying forcing.
230 C iMax
231 C jMin
232 C jMax
233 C kLev
234 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
235 _RL myCurrentTime
236 INTEGER myThid
237 CEndOfInterface
238
239 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
240 C == Local variables ==
241 C Loop counters
242 INTEGER I, J, Isl, Jsl
243 _RL vrelax,lambda_obcs_v
244 #ifndef ALLOW_CAL
245 INTEGER secondsperday
246 PARAMETER (secondsperday=86400)
247 #endif
248
249
250 IF (useOBCSsponge) THEN
251
252 C Northern Open Boundary
253 #ifdef ALLOW_OBCS_NORTH
254 DO i=iMin,iMax
255 IF ((OB_Jn(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
256 DO jsl= 1,spongeThickness
257 j=OB_Jn(i,bi,bj)-jsl
258
259 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
260 vrelax=(
261 & float(spongeThickness-jsl)*OBNv(i,kLev,bi,bj)
262 & + float(jsl)*vVel(i,j,kLev,bi,bj) )
263 & / float(spongeThickness)
264
265 lambda_obcs_v = (
266 & float(spongeThickness-jsl)*Vrelaxobcsbound
267 & + float(jsl)*Vrelaxobcsinner)
268 & / float(spongeThickness)
269
270 IF (lambda_obcs_v.ne.0.) THEN
271 lambda_obcs_v = 1. _d 0 / lambda_obcs_v
272 ELSE
273 lambda_obcs_v = 0. _d 0
274 ENDIF
275
276 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
277 & - _maskS(i,j,kLev,bi,bj) * lambda_obcs_v
278 & * ( vVel(i,j,kLev,bi,bj) - vrelax )
279 ENDIF
280
281 ENDDO
282 ENDIF
283 ENDDO
284 #endif
285
286 C Southern Open Boundary
287 #ifdef ALLOW_OBCS_SOUTH
288 DO i=iMin,iMax
289 IF ((OB_Js(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
290 DO jsl= 1,spongeThickness
291 j=OB_Js(i,bi,bj)+jsl+1
292
293 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
294 vrelax=(
295 & float(spongeThickness-jsl)*OBSv(i,kLev,bi,bj)
296 & + float(jsl)*vVel(i,j,kLev,bi,bj) )
297 & / float(spongeThickness)
298
299 lambda_obcs_v = (
300 & float(spongeThickness-jsl)*Vrelaxobcsbound
301 & + float(jsl)*Vrelaxobcsinner)
302 & / float(spongeThickness)
303
304 if (lambda_obcs_v.ne.0.) then
305 lambda_obcs_v = 1. _d 0 / lambda_obcs_v
306 else
307 lambda_obcs_v = 0. _d 0
308 endif
309
310 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
311 & - _maskS(i,j,kLev,bi,bj) * lambda_obcs_v
312 & * ( vVel(i,j,kLev,bi,bj) - vrelax )
313 ENDIF
314
315 ENDDO
316 ENDIF
317 ENDDO
318 #endif
319
320 C Eastern Open Boundary
321 #ifdef ALLOW_OBCS_EAST
322 DO j=jMin,jMax
323 IF ((OB_Ie(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
324 DO isl= 1,spongeThickness
325 i=OB_Ie(j,bi,bj)-isl
326
327 IF ((i.ge.imin).and.(i.le.imax)) THEN
328 vrelax=(
329 & float(spongeThickness-isl)*OBEv(j,kLev,bi,bj)
330 & + float(isl-1)*vVel(i,j,kLev,bi,bj) )
331 & / float(spongeThickness)
332
333 lambda_obcs_v = (
334 & float(spongeThickness-isl)*Urelaxobcsbound
335 & + float(isl)*Urelaxobcsinner)
336 & / float(spongeThickness)
337
338 if (lambda_obcs_v.ne.0.) then
339 lambda_obcs_v = 1. _d 0 / lambda_obcs_v
340 else
341 lambda_obcs_v = 0. _d 0
342 endif
343
344 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
345 & - _maskS(i,j,kLev,bi,bj) * lambda_obcs_v
346 & * ( vVel(i,j,kLev,bi,bj) - vrelax )
347 ENDIF
348
349 ENDDO
350 ENDIF
351 ENDDO
352 #endif
353
354 C Western Open Boundary
355 #ifdef ALLOW_OBCS_WEST
356 DO j=jMin,jMax
357 IF ((OB_Iw(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
358 DO isl= 1,spongeThickness
359 i=OB_Iw(j,bi,bj)+isl
360
361 IF ((i.ge.imin).and.(i.le.imax)) THEN
362
363 vrelax=(
364 & float(spongeThickness-isl)*OBWv(j,kLev,bi,bj)
365 & + float(isl)*vVel(i,j,kLev,bi,bj) )
366 & / float(spongeThickness)
367
368 lambda_obcs_v = (
369 & float(spongeThickness-isl)*Urelaxobcsbound
370 & + float(isl)*Urelaxobcsinner)
371 & / float(spongeThickness)
372
373 if (lambda_obcs_v.ne.0.) then
374 lambda_obcs_v = 1. _d 0 / lambda_obcs_v
375 else
376 lambda_obcs_v = 0. _d 0
377 endif
378
379 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
380 & - _maskS(i,j,kLev,bi,bj) * lambda_obcs_v
381 & * ( vVel(i,j,kLev,bi,bj) - vrelax )
382 ENDIF
383
384 ENDDO
385 ENDIF
386 ENDDO
387 #endif
388
389 ENDIF
390
391 #endif /* ALLOW_OBCS & ALLOW_OBCS_SPONGE */
392
393 RETURN
394 END
395
396
397 CStartOfInterface
398 SUBROUTINE OBCS_SPONGE_T(
399 I iMin, iMax, jMin, jMax,bi,bj,kLev,
400 I myCurrentTime,myThid)
401 C /==========================================================\
402 C | S/R OBCS_SPONGE_T |
403 C | o Contains problem specific forcing for zonal velocity. |
404 C |==========================================================|
405 C | Adds terms to gT for forcing by external sources |
406 C | e.g. wind stress, bottom friction etc.................. |
407 C \==========================================================/
408 IMPLICIT NONE
409
410 C == Global data ==
411 #include "SIZE.h"
412 #include "EEPARAMS.h"
413 #include "PARAMS.h"
414 #include "GRID.h"
415 #include "DYNVARS.h"
416 #include "FFIELDS.h"
417 #ifdef ALLOW_OBCS
418 # include "OBCS.h"
419 # ifdef ALLOW_CAL
420 # include "cal.h"
421 # endif
422 #endif
423 #ifdef ALLOW_AUTODIFF_TAMC
424 # include "tamc.h"
425 # include "tamc_keys.h"
426 #endif
427
428 C == Routine arguments ==
429 C iMin - Working range of tile for applying forcing.
430 C iMax
431 C jMin
432 C jMax
433 C kLev
434 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
435 _RL myCurrentTime
436 INTEGER myThid
437 CEndOfInterface
438
439 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
440 C == Local variables ==
441 C Loop counters
442 INTEGER I, J, Isl, Jsl
443 _RL trelax, lambda_obcs_t
444 #ifndef ALLOW_CAL
445 INTEGER secondsperday
446 PARAMETER (secondsperday=86400)
447 #endif
448
449
450 IF (useOBCSsponge) THEN
451
452 #ifdef ALLOW_AUTODIFF_TAMC
453 act1 = bi - myBxLo(myThid)
454 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
455 act2 = bj - myByLo(myThid)
456 max2 = myByHi(myThid) - myByLo(myThid) + 1
457 act3 = myThid - 1
458 max3 = nTx*nTy
459 act4 = ikey_dynamics - 1
460 ikey = (act1 + 1) + act2*max1
461 & + act3*max1*max2
462 & + act4*max1*max2*max3
463 kkey = (ikey-1)*Nr + klev
464 #endif /* ALLOW_AUTODIFF_TAMC */
465
466 C Northern Open Boundary
467 #ifdef ALLOW_OBCS_NORTH
468
469 #ifdef ALLOW_AUTODIFF_TAMC
470 CADJ STORE OBNt(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
471 #endif
472
473 DO i=iMin,iMax
474 IF ((OB_Jn(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
475 DO jsl= 1,spongeThickness
476 j=OB_Jn(i,bi,bj)-jsl
477
478 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
479 IF (OBNt(i,klev,bi,bj).ne. 0.d0) then
480 trelax=(
481 & float(spongeThickness-jsl)*OBNt(i,kLev,bi,bj)
482 & + float(jsl)*theta(i,j,kLev,bi,bj) )
483 & / float(spongeThickness)
484 lambda_obcs_t = (
485 & float(spongeThickness-jsl)*Vrelaxobcsbound
486 & + float(jsl)*Vrelaxobcsinner)
487 & / float(spongeThickness)
488
489 IF (lambda_obcs_t.ne.0.) THEN
490 lambda_obcs_t = 1. _d 0 / lambda_obcs_t
491 ELSE
492 lambda_obcs_t = 0. _d 0
493 ENDIF
494
495 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
496 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_t
497 & * ( theta(i,j,kLev,bi,bj) - trelax )
498 endif
499 ENDIF
500
501 ENDDO
502 ENDIF
503 ENDDO
504 #endif
505
506 C Southern Open Boundary
507 #ifdef ALLOW_OBCS_SOUTH
508
509 #ifdef ALLOW_AUTODIFF_TAMC
510 CADJ STORE OBSt(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
511 #endif
512
513 DO i=iMin,iMax
514 IF ((OB_Js(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
515 DO jsl= 1,spongeThickness
516 j=OB_Js(i,bi,bj)+jsl
517
518 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
519 IF (OBSt(i,klev,bi,bj).ne. 0.d0) then
520 trelax=(
521 & float(spongeThickness-jsl)*OBSt(i,kLev,bi,bj)
522 & + float(jsl)*theta(i,j,kLev,bi,bj) )
523 & / float(spongeThickness)
524
525 lambda_obcs_t = (
526 & float(spongeThickness-jsl)*Vrelaxobcsbound
527 & + float(jsl)*Vrelaxobcsinner)
528 & / float(spongeThickness)
529
530 if (lambda_obcs_t.ne.0.) then
531 lambda_obcs_t = 1. _d 0 / lambda_obcs_t
532 else
533 lambda_obcs_t = 0. _d 0
534 endif
535
536 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
537 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_t
538 & * ( theta(i,j,kLev,bi,bj) - trelax )
539 endif
540 ENDIF
541
542 ENDDO
543 ENDIF
544 ENDDO
545 #endif
546
547 C Eastern Open Boundary
548 #ifdef ALLOW_OBCS_EAST
549
550 #ifdef ALLOW_AUTODIFF_TAMC
551 CADJ STORE OBEt(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
552 #endif
553
554 DO j=jMin,jMax
555 IF ((OB_Ie(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
556 DO isl= 1,spongeThickness
557 i=OB_Ie(j,bi,bj)-isl
558
559 IF ((i.ge.imin).and.(i.le.imax)) THEN
560 IF (OBEt(j,klev,bi,bj).ne. 0.d0) then
561 trelax=(
562 & float(spongeThickness-isl)*OBEt(j,kLev,bi,bj)
563 & + float(isl)*theta(i,j,kLev,bi,bj) )
564 & / float(spongeThickness)
565
566 lambda_obcs_t = (
567 & float(spongeThickness-isl)*Urelaxobcsbound
568 & + float(isl)*Urelaxobcsinner)
569 & / float(spongeThickness)
570
571 if (lambda_obcs_t.ne.0.) then
572 lambda_obcs_t = 1. _d 0 / lambda_obcs_t
573 else
574 lambda_obcs_t = 0. _d 0
575 endif
576
577 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
578 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_t
579 & * ( theta(i,j,kLev,bi,bj) - trelax )
580 endif
581 ENDIF
582
583 ENDDO
584 ENDIF
585 ENDDO
586 #endif
587
588 C Western Open Boundary
589 #ifdef ALLOW_OBCS_WEST
590
591 #ifdef ALLOW_AUTODIFF_TAMC
592 CADJ STORE OBWt(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
593 #endif
594
595 DO j=jMin,jMax
596 IF ((OB_Iw(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
597 DO isl= 1,spongeThickness
598 cgg i=OB_Iw(j,bi,bj)+isl+1
599 cgg Needed to fix the coordinate of the tracer open boundary. This is the classic "cut-and-paste" bug.
600 i=OB_Iw(j,bi,bj)+isl
601
602 IF ((i.ge.imin).and.(i.le.imax)) THEN
603 IF (OBWt(j,klev,bi,bj).ne. 0.d0) then
604 trelax=(
605 & float(spongeThickness-isl)*OBWt(j,kLev,bi,bj)
606 & + float(isl)*theta(i,j,kLev,bi,bj) )
607 & / float(spongeThickness)
608
609 lambda_obcs_t= (
610 & float(spongeThickness-isl)*Urelaxobcsbound
611 & + float(isl)*Urelaxobcsinner)
612 & / float(spongeThickness)
613
614 if (lambda_obcs_t .ne. 0.) then
615 lambda_obcs_t = 1. _d 0 / lambda_obcs_t
616 else
617 lambda_obcs_t = 0. _d 0
618 endif
619
620 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
621 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_t
622 & * ( theta(i,j,kLev,bi,bj) - trelax )
623 endif
624 ENDIF
625
626 ENDDO
627 ENDIF
628 ENDDO
629 #endif
630
631 ENDIF
632
633 #endif /* ALLOW_OBCS & ALLOW_OBCS_SPONGE */
634
635 RETURN
636 END
637
638
639 CStartOfInterface
640 SUBROUTINE OBCS_SPONGE_S(
641 I iMin, iMax, jMin, jMax,bi,bj,kLev,
642 I myCurrentTime,myThid)
643 C /==========================================================\
644 C | S/R OBCS_SPONGE_U |
645 C | o Contains problem specific forcing for zonal velocity. |
646 C |==========================================================|
647 C | Adds terms to gU for forcing by external sources |
648 C | e.g. wind stress, bottom friction etc.................. |
649 C \==========================================================/
650 IMPLICIT NONE
651
652 C == Global data ==
653 #include "SIZE.h"
654 #include "EEPARAMS.h"
655 #include "PARAMS.h"
656 #include "GRID.h"
657 #include "DYNVARS.h"
658 #include "FFIELDS.h"
659 #ifdef ALLOW_OBCS
660 # include "OBCS.h"
661 # ifdef ALLOW_CAL
662 # include "cal.h"
663 # endif
664 #endif
665 #ifdef ALLOW_AUTODIFF_TAMC
666 # include "tamc.h"
667 # include "tamc_keys.h"
668 #endif
669
670 C == Routine arguments ==
671 C iMin - Working range of tile for applying forcing.
672 C iMax
673 C jMin
674 C jMax
675 C kLev
676 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
677 _RL myCurrentTime
678 INTEGER myThid
679 CEndOfInterface
680
681 #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))
682 C == Local variables ==
683 C Loop counters
684 INTEGER I, J, Isl, Jsl
685 _RL srelax, lambda_obcs_s
686 #ifndef ALLOW_CAL
687 INTEGER secondsperday
688 PARAMETER (secondsperday=86400)
689 #endif
690
691
692 IF (useOBCSsponge) THEN
693
694 #ifdef ALLOW_AUTODIFF_TAMC
695 act1 = bi - myBxLo(myThid)
696 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
697 act2 = bj - myByLo(myThid)
698 max2 = myByHi(myThid) - myByLo(myThid) + 1
699 act3 = myThid - 1
700 max3 = nTx*nTy
701 act4 = ikey_dynamics - 1
702 ikey = (act1 + 1) + act2*max1
703 & + act3*max1*max2
704 & + act4*max1*max2*max3
705 kkey = (ikey-1)*Nr + klev
706 #endif /* ALLOW_AUTODIFF_TAMC */
707
708 C Northern Open Boundary
709 #ifdef ALLOW_OBCS_NORTH
710
711 #ifdef ALLOW_AUTODIFF_TAMC
712 CADJ STORE OBNs(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
713 #endif
714
715 DO i=iMin,iMax
716 IF ((OB_Jn(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
717 DO jsl= 1,spongeThickness
718 j=OB_Jn(i,bi,bj)-jsl
719
720 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
721 IF (OBNs(i,klev,bi,bj).ne. 0.d0) then
722 srelax=(
723 & float(spongeThickness-jsl)*OBNs(i,kLev,bi,bj)
724 & + float(jsl)*salt(i,j,kLev,bi,bj) )
725 & / float(spongeThickness)
726
727 lambda_obcs_s = (
728 & float(spongeThickness-jsl)*Vrelaxobcsbound
729 & + float(jsl)*Vrelaxobcsinner)
730 & / float(spongeThickness)
731
732 IF (lambda_obcs_s.ne.0.) THEN
733 lambda_obcs_s = 1. _d 0 / lambda_obcs_s
734 ELSE
735 lambda_obcs_s = 0. _d 0
736 ENDIF
737
738 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
739 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_s
740 & * ( salt(i,j,kLev,bi,bj) - srelax )
741 endif
742 ENDIF
743
744 ENDDO
745 ENDIF
746 ENDDO
747 #endif
748
749 C Southern Open Boundary
750 #ifdef ALLOW_OBCS_SOUTH
751
752 #ifdef ALLOW_AUTODIFF_TAMC
753 CADJ STORE OBSs(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
754 #endif
755
756 DO i=iMin,iMax
757 IF ((OB_Js(i,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
758 DO jsl= 1,spongeThickness
759 j=OB_Js(i,bi,bj)+jsl
760
761 IF ((j.ge.jmin).and.(j.le.jmax)) THEN
762 IF (OBSs(i,klev,bi,bj).ne. 0.d0) then
763 srelax=(
764 & float(spongeThickness-jsl)*OBSs(i,kLev,bi,bj)
765 & + float(jsl)*salt(i,j,kLev,bi,bj) )
766 & / float(spongeThickness)
767
768 lambda_obcs_s = (
769 & float(spongeThickness)*Vrelaxobcsbound
770 & + float(jsl)*Vrelaxobcsinner)
771 & / float(spongeThickness)
772
773 if (lambda_obcs_s.ne.0.) then
774 lambda_obcs_s = 1. _d 0 / lambda_obcs_s
775 else
776 lambda_obcs_s = 0. _d 0
777 endif
778
779 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
780 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_s
781 & * ( salt(i,j,kLev,bi,bj) - srelax )
782 endif
783 ENDIF
784
785 ENDDO
786 ENDIF
787 ENDDO
788 #endif
789
790 C Eastern Open Boundary
791 #ifdef ALLOW_OBCS_EAST
792
793 #ifdef ALLOW_AUTODIFF_TAMC
794 CADJ STORE OBEs(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
795 #endif
796
797 DO j=jMin,jMax
798 IF ((OB_Ie(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
799 DO isl= 1,spongeThickness
800 i=OB_Ie(j,bi,bj)-isl
801
802 IF ((i.ge.imin).and.(i.le.imax)) THEN
803 IF (OBEs(j,klev,bi,bj).ne. 0.d0) then
804 srelax=(
805 & float(spongeThickness-isl)*OBEs(j,kLev,bi,bj)
806 & + float(isl)*salt(i,j,kLev,bi,bj) )
807 & / float(spongeThickness)
808
809 lambda_obcs_s = (
810 & float(spongeThickness-isl)*Urelaxobcsbound
811 & + float(isl)*Urelaxobcsinner)
812 & / float(spongeThickness)
813
814 if (lambda_obcs_s.ne.0.) then
815 lambda_obcs_s = 1. _d 0 / lambda_obcs_s
816 else
817 lambda_obcs_s = 0. _d 0
818 endif
819
820 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
821 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_s
822 & * ( salt(i,j,kLev,bi,bj) - srelax )
823 endif
824 ENDIF
825
826 ENDDO
827 ENDIF
828 ENDDO
829 #endif
830
831 C Western Open Boundary
832 #ifdef ALLOW_OBCS_WEST
833
834 #ifdef ALLOW_AUTODIFF_TAMC
835 CADJ STORE OBWs(:,klev,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
836 #endif
837
838 DO j=jMin,jMax
839 IF ((OB_Iw(j,bi,bj).ne.0).and.(spongeThickness.ne.0)) THEN
840 DO isl= 1,spongeThickness
841
842 cgg i=OB_Iw(j,bi,bj)+isl+1
843 cgg Fix the tracer o.b. coordinate.
844 i=OB_Iw(j,bi,bj)+isl
845
846 IF ((i.ge.imin).and.(i.le.imax)) THEN
847 IF (OBWs(j,klev,bi,bj).ne. 0.d0) then
848 srelax=(
849 & float(spongeThickness-isl)*OBWs(j,kLev,bi,bj)
850 & + float(isl)*salt(i,j,kLev,bi,bj) )
851 & / float(spongeThickness)
852
853 lambda_obcs_s= (
854 & float(spongeThickness-isl)*Urelaxobcsbound
855 & + float(isl)*Urelaxobcsinner)
856 & / float(spongeThickness)
857
858 if (lambda_obcs_s.ne.0.) then
859 lambda_obcs_s = 1. _d 0 / lambda_obcs_s
860 else
861 lambda_obcs_s = 0. _d 0
862 endif
863
864 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
865 & - maskC(i,j,kLev,bi,bj) * lambda_obcs_s
866 & * ( salt(i,j,kLev,bi,bj) - srelax )
867 endif
868 ENDIF
869
870 ENDDO
871 ENDIF
872 ENDDO
873 #endif
874
875 ENDIF
876
877 #endif /* ALLOW_OBCS & ALLOW_OBCS_SPONGE */
878
879 RETURN
880 END
  ViewVC Help
Powered by ViewVC 1.1.22