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 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

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