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

source: trunk/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 3319

Last change on this file since 3319 was 3294, checked in by rblod, 12 years ago

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