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

Last change on this file since 4528 was 4153, checked in by cetlod, 7 years ago

dev_LOCEAN_2013: merge in trunk changes between r3940 and r4028, see ticket #1169

  • Property svn:keywords set to Id
File size: 16.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      LOGICAL  :: ll_spdone
188      REAL(wp) :: z1spongearea, zramp
189      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
190
191#if defined SPONGE || defined SPONGE_TOP
192      ll_spdone=.TRUE.
193      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
194         ! Define ramp from boundaries towards domain interior
195         ! at T-points
196         ! Store it in ztabramp
197         ll_spdone=.FALSE.
198
199         CALL wrk_alloc( jpi, jpj, ztabramp )
200
201         ispongearea  = 2 + 2 * Agrif_irhox()
202         ilci = nlci - ispongearea
203         ilcj = nlcj - ispongearea 
204         z1spongearea = 1._wp / REAL( ispongearea - 2 )
205         spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
206
207         ztabramp(:,:) = 0.
208
209         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
210            DO jj = 1, jpj
211               IF ( umask(2,jj,1) == 1._wp ) THEN
212                 DO ji = 2, ispongearea                 
213                    ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
214                 END DO
215               ENDIF
216            ENDDO
217         ENDIF
218
219         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
220            DO jj = 1, jpj
221               IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
222                  DO ji = ilci+1,nlci-1
223                     zramp = (ji - (ilci+1) ) * z1spongearea
224                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
225                  ENDDO
226               ENDIF
227            ENDDO
228         ENDIF
229
230         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
231            DO ji = 1, jpi
232               IF ( vmask(ji,2,1) == 1._wp ) THEN
233                  DO jj = 2, ispongearea
234                     zramp = ( ispongearea-jj ) * z1spongearea
235                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
236                  END DO
237               ENDIF
238            ENDDO
239         ENDIF
240
241         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
242            DO ji = 1, jpi
243               IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
244                  DO jj = ilcj+1,nlcj-1
245                     zramp = (jj - (ilcj+1) ) * z1spongearea
246                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
247                  END DO
248               ENDIF
249            ENDDO
250         ENDIF
251
252      ENDIF
253
254      ! Tracers
255      IF( .NOT. spongedoneT ) THEN
256         spe1ur(:,:) = 0.
257         spe2vr(:,:) = 0.
258
259         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
260            spe1ur(2:ispongearea-1,:       ) = visc_tra                                        &
261               &                             *    0.5 * (  ztabramp(2:ispongearea-1,:      )   &
262               &                                         + ztabramp(3:ispongearea  ,:      ) ) &
263               &                             * e2u(2:ispongearea-1,:) / e1u(2:ispongearea-1,:)
264
265            spe2vr(2:ispongearea  ,1:jpjm1 ) = visc_tra                                        &
266               &                             *    0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1)   &
267               &                                         + ztabramp(2:ispongearea,2  :jpj  ) ) &
268               &                             * e1v(2:ispongearea,1:jpjm1) / e2v(2:ispongearea,1:jpjm1)
269         ENDIF
270
271         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
272            spe1ur(ilci+1:nlci-2,:        ) = visc_tra                                   &
273               &                            * 0.5 * (  ztabramp(ilci+1:nlci-2,:      )   & 
274               &                                     + ztabramp(ilci+2:nlci-1,:      ) ) &
275               &                            * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:)
276
277            spe2vr(ilci+1:nlci-1,1:jpjm1  )  = visc_tra                                  &
278               &                            * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1)   & 
279               &                                     + ztabramp(ilci+1:nlci-1,2:jpj  ) ) & 
280               &                            * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1)
281         ENDIF
282
283         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
284            spe1ur(1:jpim1,2:ispongearea  ) = visc_tra                                     &
285               &                            * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  )   & 
286               &                                     + ztabramp(2:jpi  ,2:ispongearea  ) ) &
287               &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea)
288   
289            spe2vr(:      ,2:ispongearea-1) = visc_tra                                     &
290               &                            * 0.5 * (  ztabramp(:      ,2:ispongearea-1)   &
291               &                                     + ztabramp(:      ,3:ispongearea  ) ) &
292               &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1)
293         ENDIF
294
295         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
296            spe1ur(1:jpim1,ilcj+1:nlcj-1) = visc_tra                                   &
297               &                          * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1)   &
298               &                                   + ztabramp(2:jpi  ,ilcj+1:nlcj-1) ) &
299               &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1)
300
301            spe2vr(:      ,ilcj+1:nlcj-2) = visc_tra                                   &
302               &                          * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2)   &
303               &                                   + ztabramp(:      ,ilcj+2:nlcj-1) ) &
304               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2)
305         ENDIF
306         spongedoneT = .TRUE.
307      ENDIF
308
309      ! Dynamics
310      IF( .NOT. spongedoneU ) THEN
311         spe1ur2(:,:) = 0.
312         spe2vr2(:,:) = 0.
313
314         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
315            spe1ur2(2:ispongearea-1,:      ) = visc_dyn                                   &
316               &                             * 0.5 * (  ztabramp(2:ispongearea-1,:      ) &
317               &                                      + ztabramp(3:ispongearea  ,:      ) )
318            spe2vr2(2:ispongearea  ,1:jpjm1) = visc_dyn                                   &
319               &                             * 0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1) &
320               &                                      + ztabramp(2:ispongearea  ,2:jpj  ) ) 
321         ENDIF
322
323         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
324            spe1ur2(ilci+1:nlci-2  ,:      ) = visc_dyn                                   &
325               &                             * 0.5 * (  ztabramp(ilci+1:nlci-2, :       ) &
326               &                                      + ztabramp(ilci+2:nlci-1, :       ) )                     
327            spe2vr2(ilci+1:nlci-1  ,1:jpjm1) = visc_dyn                                   &
328               &                             * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1  ) &
329               &                                      + ztabramp(ilci+1:nlci-1,2:jpj    ) ) 
330         ENDIF
331
332         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
333            spe1ur2(1:jpim1,2:ispongearea  ) = visc_dyn                                   & 
334               &                             * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  ) &
335               &                                      + ztabramp(2:jpi  ,2:ispongearea  ) ) 
336            spe2vr2(:      ,2:ispongearea-1) = visc_dyn                                   &
337               &                             * 0.5 * (  ztabramp(:      ,2:ispongearea-1) &
338               &                                      + ztabramp(:      ,3:ispongearea  ) )
339         ENDIF
340
341         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
342            spe1ur2(1:jpim1,ilcj+1:nlcj-1  ) = visc_dyn                                   &
343               &                             * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1  ) &
344               &                                      + ztabramp(2:jpi  ,ilcj+1:nlcj-1  ) ) 
345            spe2vr2(:      ,ilcj+1:nlcj-2  ) = visc_dyn                                   &
346               &                             * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2  ) &
347               &                                      + ztabramp(:      ,ilcj+2:nlcj-1  ) )
348         ENDIF
349         spongedoneU = .TRUE.
350         spbtr3(:,:) = 1. / ( e1f(:,:) * e2f(:,:) )
351      ENDIF
352      !
353      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
354      !
355#endif
356
357   END SUBROUTINE Agrif_Sponge
358
359   SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2)
360      !!---------------------------------------------
361      !!   *** ROUTINE interptsn ***
362      !!---------------------------------------------
363      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
364      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
365
366      tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
367
368   END SUBROUTINE interptsn
369
370   SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
371      !!---------------------------------------------
372      !!   *** ROUTINE interpun ***
373      !!---------------------------------------------
374      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
375      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
376
377      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
378
379   END SUBROUTINE interpun
380
381   SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
382      !!---------------------------------------------
383      !!   *** ROUTINE interpvn ***
384      !!---------------------------------------------
385      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
386      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
387
388      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
389
390   END SUBROUTINE interpvn
391
392#else
393CONTAINS
394
395   SUBROUTINE agrif_opa_sponge_empty
396      !!---------------------------------------------
397      !!   *** ROUTINE agrif_OPA_sponge_empty ***
398      !!---------------------------------------------
399      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
400   END SUBROUTINE agrif_opa_sponge_empty
401#endif
402
403END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.