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 branches/dev_001_SBC/NEMO/NST_SRC – NEMO

source: branches/dev_001_SBC/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 811

Last change on this file since 811 was 811, checked in by ctlod, 16 years ago

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