New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_opa_sponge.F90 in trunk/NEMO/NST_SRC – NEMO

source: trunk/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 1146

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago

Add svn Id (first try), see ticket #210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 15.1 KB
Line 
1#define SPONGE
2
3Module agrif_opa_sponge
4#if defined key_agrif
5   USE par_oce
6   USE oce
7   USE dom_oce
8   USE in_out_manager
9   USE agrif_oce
10
11   IMPLICIT NONE
12   PRIVATE
13
14   PUBLIC Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptn, interpsn, interpun, interpvn
15
16
17   CONTAINS
18
19   SUBROUTINE Agrif_Sponge_Tra
20      !!---------------------------------------------
21      !!   *** ROUTINE Agrif_Sponge_Tra ***
22      !!---------------------------------------------
23#include "domzgr_substitute.h90"
24
25      INTEGER :: ji,jj,jk
26      INTEGER :: spongearea
27      REAL(wp) :: timecoeff
28      REAL(wp) :: zta, zsa, zabe1, zabe2, zbtr
29      REAL(wp), DIMENSION(jpi,jpj) :: localviscsponge
30      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tbdiff, sbdiff
31      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu ,ztv, zsu ,zsv
32      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztab
33
34#if defined SPONGE
35
36      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
37
38      Agrif_SpecialValue=0.
39      Agrif_UseSpecialValue = .TRUE.
40      ztab = 0.e0
41      CALL Agrif_Bc_Variable(ztab, ta,calledweight=timecoeff,procname=interptn)
42      Agrif_UseSpecialValue = .FALSE.
43
44      tbdiff(:,:,:) = tb(:,:,:) - ztab(:,:,:)
45
46      ztab = 0.e0
47      Agrif_SpecialValue=0.
48      Agrif_UseSpecialValue = .TRUE.
49      CALL Agrif_Bc_Variable(ztab, sa,calledweight=timecoeff,procname=interpsn)
50      Agrif_UseSpecialValue = .FALSE.
51
52      sbdiff(:,:,:) = sb(:,:,:) - ztab(:,:,:)
53
54
55      spongearea = 2 + 2 * Agrif_irhox()
56
57      localviscsponge = 0.
58     
59      IF (.NOT. spongedoneT) THEN
60         spe1ur(:,:) = 0.
61         spe2vr(:,:) = 0.
62
63      IF ((nbondi == -1).OR.(nbondi == 2)) THEN
64         DO ji = 2, spongearea
65            localviscsponge(ji,:) = visc_tra * (spongearea-ji)/real(spongearea-2)
66         ENDDO
67   
68    spe1ur(2:spongearea-1,:)=0.5 * (localviscsponge(2:spongearea-1,:) + localviscsponge(3:spongearea,:)) &
69          * e2u(2:spongearea-1,:) / e1u(2:spongearea-1,:)
70
71         spe2vr(2:spongearea,1:jpjm1) = 0.5 * (localviscsponge(2:spongearea,1:jpjm1) + &
72             localviscsponge(2:spongearea,2:jpj)) &
73           * e1v(2:spongearea,1:jpjm1) / e2v(2:spongearea,1:jpjm1)
74      ENDIF
75
76      IF ((nbondi == 1).OR.(nbondi == 2)) THEN
77         DO ji = nlci-spongearea + 1,nlci-1
78            localviscsponge(ji,:) = visc_tra * (ji - (nlci-spongearea+1))/real(spongearea-2)
79         ENDDO
80   
81    spe1ur(nlci-spongearea + 1:nlci-2,:)=0.5 * (localviscsponge(nlci-spongearea + 1:nlci-2,:) + &
82           localviscsponge(nlci-spongearea + 2:nlci-1,:)) &
83          * e2u(nlci-spongearea + 1:nlci-2,:) / e1u(nlci-spongearea + 1:nlci-2,:)
84
85         spe2vr(nlci-spongearea + 1:nlci-1,1:jpjm1) = 0.5 * (localviscsponge(nlci-spongearea + 1:nlci-1,1:jpjm1) &
86              + localviscsponge(nlci-spongearea + 1:nlci-1,2:jpj)) &
87           * e1v(nlci-spongearea + 1:nlci-1,1:jpjm1) / e2v(nlci-spongearea + 1:nlci-1,1:jpjm1)
88      ENDIF
89
90
91      IF ((nbondj == -1).OR.(nbondj == 2)) THEN
92         DO jj = 2, spongearea
93            localviscsponge(:,jj) = visc_tra * (spongearea-jj)/real(spongearea-2)
94         ENDDO
95   
96    spe1ur(1:jpim1,2:spongearea)=0.5 * (localviscsponge(1:jpim1,2:spongearea) + &
97           localviscsponge(2:jpi,2:spongearea)) &
98          * e2u(1:jpim1,2:spongearea) / e1u(1:jpim1,2:spongearea)
99
100         spe2vr(:,2:spongearea-1) = 0.5 * (localviscsponge(:,2:spongearea-1) + &
101             localviscsponge(:,3:spongearea)) &
102           * e1v(:,2:spongearea-1) / e2v(:,2:spongearea-1)
103      ENDIF
104
105      IF ((nbondj == 1).OR.(nbondj == 2)) THEN
106         DO jj = nlcj-spongearea + 1,nlcj-1
107            localviscsponge(:,jj) = visc_tra * (jj - (nlcj-spongearea+1))/real(spongearea-2)
108         ENDDO
109   
110    spe1ur(1:jpim1,nlcj-spongearea + 1:nlcj-1)=0.5 * (localviscsponge(1:jpim1,nlcj-spongearea + 1:nlcj-1) + &
111            localviscsponge(2:jpi,nlcj-spongearea + 1:nlcj-1)) &
112          * e2u(1:jpim1,nlcj-spongearea + 1:nlcj-1) / e1u(1:jpim1,nlcj-spongearea + 1:nlcj-1)
113
114         spe2vr(:,nlcj-spongearea + 1:nlcj-2) = 0.5 * (localviscsponge(:,nlcj-spongearea + 1:nlcj-2) + &
115            localviscsponge(:,nlcj-spongearea + 2:nlcj-1)) &
116           * e1v(:,nlcj-spongearea + 1:nlcj-2) / e2v(:,nlcj-spongearea + 1:nlcj-2)
117      ENDIF
118     
119         spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:))
120
121         spongedoneT = .TRUE.
122      ENDIF
123
124      DO jk = 1, jpkm1
125         DO jj = 1, jpjm1
126            DO ji = 1, jpim1
127#if defined key_zco
128               zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj)
129               zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj)
130#else
131               zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk)
132               zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk)
133#endif
134               ztu(ji,jj,jk) = zabe1 * ( tbdiff(ji+1,jj  ,jk) - tbdiff(ji,jj,jk) )
135               zsu(ji,jj,jk) = zabe1 * ( sbdiff(ji+1,jj  ,jk) - sbdiff(ji,jj,jk) )
136               ztv(ji,jj,jk) = zabe2 * ( tbdiff(ji  ,jj+1,jk) - tbdiff(ji,jj,jk) )
137               zsv(ji,jj,jk) = zabe2 * ( sbdiff(ji  ,jj+1,jk) - sbdiff(ji,jj,jk) )
138            ENDDO
139         ENDDO
140
141         DO jj = 2,jpjm1
142            DO ji = 2,jpim1
143#if defined key_zco
144               zbtr = spbtr2(ji,jj)
145#else
146               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
147#endif
148               ! horizontal diffusive trends
149               zta = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
150                  &          + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
151               zsa = zbtr * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk)   &
152                  &          + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
153               ! add it to the general tracer trends
154               ta(ji,jj,jk) = (ta(ji,jj,jk) + zta)
155               sa(ji,jj,jk) = (sa(ji,jj,jk) + zsa)
156            END DO
157         END DO
158
159      ENDDO
160
161#endif
162
163   END SUBROUTINE Agrif_Sponge_Tra
164
165   SUBROUTINE Agrif_Sponge_dyn
166      !!---------------------------------------------
167      !!   *** ROUTINE Agrif_Sponge_dyn ***
168      !!---------------------------------------------
169#include "domzgr_substitute.h90"
170
171      INTEGER :: ji,jj,jk
172      INTEGER :: spongearea
173      REAL(wp) :: timecoeff
174      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
175      REAL(wp), DIMENSION(jpi,jpj) :: localviscsponge
176      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztab, ubdiff, vbdiff,rotdiff,hdivdiff
177
178#if defined SPONGE
179
180      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
181
182      Agrif_SpecialValue=0.
183      Agrif_UseSpecialValue = ln_spc_dyn
184      ztab = 0.e0
185      CALL Agrif_Bc_Variable(ztab, ua,calledweight=timecoeff,procname=interpun)
186      Agrif_UseSpecialValue = .FALSE.
187
188      ubdiff(:,:,:) = (ub(:,:,:) - ztab(:,:,:))*umask(:,:,:)
189
190      ztab = 0.e0
191      Agrif_SpecialValue=0.
192      Agrif_UseSpecialValue = ln_spc_dyn
193      CALL Agrif_Bc_Variable(ztab, va,calledweight=timecoeff,procname=interpvn)
194      Agrif_UseSpecialValue = .FALSE.
195
196      vbdiff(:,:,:) = (vb(:,:,:) - ztab(:,:,:))*vmask(:,:,:)
197
198      spongearea = 2 + 2 * Agrif_irhox()
199
200      localviscsponge = 0.
201
202      IF (.NOT. spongedoneU) THEN
203         spe1ur2(:,:) = 0.
204         spe2vr2(:,:) = 0.
205
206      IF ((nbondi == -1).OR.(nbondi == 2)) THEN
207         DO ji = 2, spongearea
208            localviscsponge(ji,:) = visc_dyn * (spongearea-ji)/real(spongearea-2)
209         ENDDO
210   
211    spe1ur2(2:spongearea-1,:)=0.5 * (localviscsponge(2:spongearea-1,:) + localviscsponge(3:spongearea,:))
212
213         spe2vr2(2:spongearea,1:jpjm1) = 0.5 * (localviscsponge(2:spongearea,1:jpjm1) + &
214             localviscsponge(2:spongearea,2:jpj))
215      ENDIF
216
217      IF ((nbondi == 1).OR.(nbondi == 2)) THEN
218         DO ji = nlci-spongearea + 1,nlci-1
219            localviscsponge(ji,:) = visc_dyn * (ji - (nlci-spongearea+1))/real(spongearea-2)
220         ENDDO
221   
222    spe1ur2(nlci-spongearea + 1:nlci-2,:)=0.5 * (localviscsponge(nlci-spongearea + 1:nlci-2,:) + &
223           localviscsponge(nlci-spongearea + 2:nlci-1,:))
224
225         spe2vr2(nlci-spongearea + 1:nlci-1,1:jpjm1) = 0.5 * (localviscsponge(nlci-spongearea + 1:nlci-1,1:jpjm1) &
226              + localviscsponge(nlci-spongearea + 1:nlci-1,2:jpj))
227      ENDIF
228
229
230      IF ((nbondj == -1).OR.(nbondj == 2)) THEN
231         DO jj = 2, spongearea
232            localviscsponge(:,jj) = visc_dyn * (spongearea-jj)/real(spongearea-2)
233         ENDDO
234   
235    spe1ur2(1:jpim1,2:spongearea)=0.5 * (localviscsponge(1:jpim1,2:spongearea) + &
236           localviscsponge(2:jpi,2:spongearea))
237
238         spe2vr2(:,2:spongearea-1) = 0.5 * (localviscsponge(:,2:spongearea-1) + &
239             localviscsponge(:,3:spongearea))
240      ENDIF
241
242      IF ((nbondj == 1).OR.(nbondj == 2)) THEN
243         DO jj = nlcj-spongearea + 1,nlcj-1
244            localviscsponge(:,jj) = visc_dyn * (jj - (nlcj-spongearea+1))/real(spongearea-2)
245         ENDDO
246   
247    spe1ur2(1:jpim1,nlcj-spongearea + 1:nlcj-1)=0.5 * (localviscsponge(1:jpim1,nlcj-spongearea + 1:nlcj-1) + &
248            localviscsponge(2:jpi,nlcj-spongearea + 1:nlcj-1))
249
250         spe2vr2(:,nlcj-spongearea + 1:nlcj-2) = 0.5 * (localviscsponge(:,nlcj-spongearea + 1:nlcj-2) + &
251            localviscsponge(:,nlcj-spongearea + 2:nlcj-1))
252      ENDIF
253
254         spongedoneU = .TRUE.
255   
256    spbtr3(:,:) = 1./( e1f(:,:) * e2f(:,:))
257      ENDIF
258     
259      IF (.NOT. spongedoneT) THEN
260        spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:))     
261      ENDIF
262     
263      DO jk=1,jpkm1
264      ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:)
265      vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:)
266      ENDDO
267     
268      hdivdiff = 0.
269      rotdiff = 0.
270
271      DO jk = 1, jpkm1                                 ! Horizontal slab
272         !                                             ! ===============
273
274         !                                             ! --------
275         ! Horizontal divergence                       !   div
276         !                                             ! --------
277         DO jj = 2, jpjm1
278            DO ji = 2, jpim1   ! vector opt.
279#if defined key_zco
280               zbtr = spbtr2(ji,jj)
281               hdivdiff(ji,jj,jk) = (  e2u(ji,jj) * ubdiff(ji,jj,jk) &
282                  - e2u(ji-1,jj  ) * ubdiff(ji-1,jj  ,jk)      &
283                  &               + e1v(ji,jj) * vbdiff(ji,jj,jk) - &
284                  &              e1v(ji  ,jj-1) * vbdiff(ji  ,jj-1,jk)  ) * zbtr
285#else
286               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
287               hdivdiff(ji,jj,jk) =   &
288                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * & 
289                  ubdiff(ji,jj,jk) - e2u(ji-1,jj  )* &
290                  fse3u(ji-1,jj  ,jk)  * ubdiff(ji-1,jj  ,jk)       &
291                  + e1v(ji,jj)*fse3v(ji,jj,jk) * &
292                  vbdiff(ji,jj,jk) - e1v(ji  ,jj-1)* &
293                  fse3v(ji  ,jj-1,jk)  * vbdiff(ji  ,jj-1,jk)  ) * zbtr
294#endif
295            END DO
296         END DO
297
298         DO jj = 1, jpjm1
299            DO ji = 1, jpim1   ! vector opt.
300#if defined key_zco     
301               zbtr = spbtr3(ji,jj)
302               rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    &
303                  &              - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
304                  &           * fmask(ji,jj,jk) * zbtr
305#else
306               zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk)
307               rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    &
308                  &              - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
309                  &           * fmask(ji,jj,jk) * zbtr
310#endif       
311            END DO
312         END DO
313
314      ENDDO
315
316      !                                                ! ===============
317      DO jk = 1, jpkm1                                 ! Horizontal slab
318         !                                             ! ===============
319         DO jj = 2, jpjm1
320            DO ji = 2, jpim1   ! vector opt.
321#if defined key_zco
322               ! horizontal diffusive trends
323               ze2u = rotdiff (ji,jj,jk)
324               ze1v = hdivdiff(ji,jj,jk)
325               zua = - (                ze2u                  - &
326                  rotdiff (ji,jj-1,jk) ) / e2u(ji,jj)   &
327                  + ( hdivdiff(ji+1,jj,jk) -     &
328                  ze1v                  ) / e1u(ji,jj)
329
330               zva = + (                ze2u                  - &
331                  rotdiff (ji-1,jj,jk) ) / e1v(ji,jj)   &
332                  + ( hdivdiff(ji,jj+1,jk) -       &
333                  ze1v                  ) / e2v(ji,jj)
334#else
335               ze2u = rotdiff (ji,jj,jk)
336               ze1v = hdivdiff(ji,jj,jk)
337               ! horizontal diffusive trends
338               zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
339                  + ( hdivdiff(ji+1,jj,jk) - ze1v      &
340                  ) / e1u(ji,jj)
341
342               zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
343                  + ( hdivdiff(ji,jj+1,jk) - ze1v    &
344                  ) / e2v(ji,jj)
345#endif
346
347               ! add it to the general momentum trends
348               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
349               va(ji,jj,jk) = va(ji,jj,jk) + zva
350            END DO
351         END DO
352         !                                             ! ===============
353      END DO                                           !   End of slab
354      !                                                ! ===============
355
356#endif
357
358   END SUBROUTINE Agrif_Sponge_dyn
359
360   SUBROUTINE interptn(tabres,i1,i2,j1,j2,k1,k2)
361      !!---------------------------------------------
362      !!   *** ROUTINE interptn ***
363      !!---------------------------------------------
364#  include "domzgr_substitute.h90"       
365     
366      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
367      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
368
369      tabres(i1:i2,j1:j2,k1:k2) = tn(i1:i2,j1:j2,k1:k2)
370
371   END SUBROUTINE interptn
372
373   SUBROUTINE interpsn(tabres,i1,i2,j1,j2,k1,k2)
374      !!---------------------------------------------
375      !!   *** ROUTINE interpsn ***
376      !!---------------------------------------------
377#  include "domzgr_substitute.h90"       
378     
379      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
380      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
381
382      tabres(i1:i2,j1:j2,k1:k2) = sn(i1:i2,j1:j2,k1:k2)
383
384   END SUBROUTINE interpsn
385
386   SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
387      !!---------------------------------------------
388      !!   *** ROUTINE interpun ***
389      !!---------------------------------------------
390#  include "domzgr_substitute.h90"       
391     
392      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
393      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
394
395      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
396
397   END SUBROUTINE interpun
398
399   SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
400      !!---------------------------------------------
401      !!   *** ROUTINE interpvn ***
402      !!---------------------------------------------
403#  include "domzgr_substitute.h90"       
404     
405      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
406      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
407
408      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
409
410   END SUBROUTINE interpvn
411
412#else
413CONTAINS
414
415   SUBROUTINE agrif_opa_sponge_empty
416      !!---------------------------------------------
417      !!   *** ROUTINE agrif_OPA_sponge_empty ***
418      !!---------------------------------------------
419      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
420   END SUBROUTINE agrif_opa_sponge_empty
421#endif
422
423END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.