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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 3750

Last change on this file since 3750 was 3698, checked in by rblod, 11 years ago

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