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