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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

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