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

Last change on this file since 706 was 699, checked in by smasson, 17 years ago

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