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

source: branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 8762

Last change on this file since 8762 was 8762, checked in by jchanut, 6 years ago

AGRIF + vvl: Final changes, update SETTE tests (these are ok except SAS) - #1965
  • Property svn:keywords set to Id
File size: 16.5 KB
Line 
1#define SPONGE && define SPONGE_TOP
2
3MODULE agrif_opa_sponge
4   !!======================================================================
5   !!                ***  MODULE agrif_opa_update  ***
6   !! AGRIF :   
7   !!======================================================================
8   !! History : 
9   !!----------------------------------------------------------------------
10#if defined key_agrif
11   USE par_oce
12   USE oce
13   USE dom_oce
14   USE in_out_manager
15   USE agrif_oce
16   USE wrk_nemo 
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
23   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
24
25   !!----------------------------------------------------------------------
26   !! NEMO/NST 3.7 , NEMO Consortium (2015)
27   !! $Id$
28   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
29   !!----------------------------------------------------------------------
30CONTAINS
31
32   SUBROUTINE Agrif_Sponge_Tra
33      !!---------------------------------------------
34      !!   *** ROUTINE Agrif_Sponge_Tra ***
35      !!---------------------------------------------
36      REAL(wp) :: timecoeff
37      !!---------------------------------------------
38      !
39#if defined SPONGE
40!!      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
41!! Assume persistence:
42      timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
43
44      CALL Agrif_Sponge
45      Agrif_SpecialValue=0.
46      Agrif_UseSpecialValue = .TRUE.
47      tabspongedone_tsn = .FALSE.
48
49      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
50
51      Agrif_UseSpecialValue = .FALSE.
52#endif
53      !
54   END SUBROUTINE Agrif_Sponge_Tra
55
56
57   SUBROUTINE Agrif_Sponge_dyn
58      !!---------------------------------------------
59      !!   *** ROUTINE Agrif_Sponge_dyn ***
60      !!---------------------------------------------
61      REAL(wp) :: timecoeff
62      !!---------------------------------------------
63
64#if defined SPONGE
65!!      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
66!! Assume persistence:
67      timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
68
69      Agrif_SpecialValue=0.
70      Agrif_UseSpecialValue = ln_spc_dyn
71
72      tabspongedone_u = .FALSE.
73      tabspongedone_v = .FALSE.         
74      CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
75
76      tabspongedone_u = .FALSE.
77      tabspongedone_v = .FALSE.
78      CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
79
80      Agrif_UseSpecialValue = .FALSE.
81#endif
82      !
83   END SUBROUTINE Agrif_Sponge_dyn
84
85
86   SUBROUTINE Agrif_Sponge
87      !!---------------------------------------------
88      !!   *** ROUTINE  Agrif_Sponge ***
89      !!---------------------------------------------
90      INTEGER  :: ji,jj,jk
91      INTEGER  :: ispongearea, ilci, ilcj
92      LOGICAL  :: ll_spdone
93      REAL(wp) :: z1spongearea, zramp
94      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
95
96#if defined SPONGE || defined SPONGE_TOP
97      ll_spdone=.TRUE.
98      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
99         ! Define ramp from boundaries towards domain interior
100         ! at T-points
101         ! Store it in ztabramp
102         ll_spdone=.FALSE.
103
104         CALL wrk_alloc( jpi, jpj, ztabramp )
105
106         ispongearea  = 2 + nn_sponge_len * Agrif_irhox()
107         ilci = nlci - ispongearea
108         ilcj = nlcj - ispongearea 
109         z1spongearea = 1._wp / REAL( ispongearea - 2 )
110
111         ztabramp(:,:) = 0._wp
112
113         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
114            DO jj = 1, jpj
115               IF ( umask(2,jj,1) == 1._wp ) THEN
116                 DO ji = 2, ispongearea                 
117                    ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
118                 END DO
119               ENDIF
120            ENDDO
121         ENDIF
122
123         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
124            DO jj = 1, jpj
125               IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
126                  DO ji = ilci+1,nlci-1
127                     zramp = (ji - (ilci+1) ) * z1spongearea
128                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
129                  ENDDO
130               ENDIF
131            ENDDO
132         ENDIF
133
134         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
135            DO ji = 1, jpi
136               IF ( vmask(ji,2,1) == 1._wp ) THEN
137                  DO jj = 2, ispongearea
138                     zramp = ( ispongearea-jj ) * z1spongearea
139                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
140                  END DO
141               ENDIF
142            ENDDO
143         ENDIF
144
145         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
146            DO ji = 1, jpi
147               IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
148                  DO jj = ilcj+1,nlcj-1
149                     zramp = (jj - (ilcj+1) ) * z1spongearea
150                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
151                  END DO
152               ENDIF
153            ENDDO
154         ENDIF
155
156      ENDIF
157
158      ! Tracers
159      IF( .NOT. spongedoneT ) THEN
160         fsaht_spu(:,:) = 0._wp
161         fsaht_spv(:,:) = 0._wp
162         DO jj = 2, jpjm1
163            DO ji = 2, jpim1   ! vector opt.
164               fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj  ))
165               fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji  ,jj+1))
166            END DO
167         END DO
168
169         CALL lbc_lnk( fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
170         CALL lbc_lnk( fsaht_spv, 'V', 1. )
171         spongedoneT = .TRUE.
172      ENDIF
173
174      ! Dynamics
175      IF( .NOT. spongedoneU ) THEN
176         fsahm_spt(:,:) = 0._wp
177         fsahm_spf(:,:) = 0._wp
178         DO jj = 2, jpjm1
179            DO ji = 2, jpim1   ! vector opt.
180               fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
181               fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) &
182                                                     &  +ztabramp(ji,jj) + ztabramp(ji+1,jj  ) )
183            END DO
184         END DO
185
186         CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
187         CALL lbc_lnk( fsahm_spf, 'F', 1. )
188         spongedoneU = .TRUE.
189      ENDIF
190      !
191      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
192      !
193#endif
194      !
195   END SUBROUTINE Agrif_Sponge
196
197
198   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
199      !!---------------------------------------------
200      !!   *** ROUTINE interptsn_sponge ***
201      !!---------------------------------------------
202      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
203      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
204      LOGICAL, INTENT(in) :: before
205      !
206      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
207      INTEGER  ::   iku, ikv
208      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
209      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
210      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
211      !
212      IF( before ) THEN
213         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsb(i1:i2,j1:j2,k1:k2,n1:n2)
214      ELSE   
215         !
216         tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)   
217         DO jn = 1, jpts           
218            DO jk = 1, jpkm1
219               DO jj = j1,j2-1
220                  DO ji = i1,i2-1
221                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
222                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
223                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
224                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
225                  END DO
226               END DO
227               !
228               IF( ln_zps ) THEN      ! set gradient at partial step level
229                  DO jj = j1,j2-1
230                     DO ji = i1,i2-1
231                        ! last level
232                        iku = mbku(ji,jj)
233                        ikv = mbkv(ji,jj)
234                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
235                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
236                     END DO
237                  END DO
238               ENDIF
239            END DO
240            !
241            DO jk = 1, jpkm1
242               DO jj = j1+1,j2-1
243                  DO ji = i1+1,i2-1
244                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
245                        zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
246                        ! horizontal diffusive trends
247                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  )
248                        ! add it to the general tracer trends
249                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
250                     ENDIF
251                  END DO
252               END DO
253            END DO
254            !
255         END DO
256         !
257         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
258         !
259      ENDIF
260      !
261   END SUBROUTINE interptsn_sponge
262
263
264   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
265      !!---------------------------------------------
266      !!   *** ROUTINE interpun_sponge ***
267      !!---------------------------------------------   
268      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
269      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
270      LOGICAL, INTENT(in) :: before
271
272      INTEGER :: ji,jj,jk
273
274      ! sponge parameters
275      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
276      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
277      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
278      INTEGER :: jmax
279      !!---------------------------------------------   
280      !
281      IF( before ) THEN
282         tabres = un(i1:i2,j1:j2,:)
283      ELSE
284         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
285         !
286         DO jk = 1, jpkm1                                 ! Horizontal slab
287            !                                             ! ===============
288
289            !                                             ! --------
290            ! Horizontal divergence                       !   div
291            !                                             ! --------
292            DO jj = j1,j2
293               DO ji = i1+1,i2   ! vector opt.
294                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
295                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
296                                     &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
297               END DO
298            END DO
299
300            DO jj = j1,j2-1
301               DO ji = i1,i2   ! vector opt.
302                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
303                  rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
304                                       +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) & 
305                                    & ) * fmask(ji,jj,jk) * zbtr 
306               END DO
307            END DO
308         END DO
309         !
310         DO jj = j1+1, j2-1
311            DO ji = i1+1, i2-1   ! vector opt.
312
313               IF (.NOT. tabspongedone_u(ji,jj)) THEN
314                  DO jk = 1, jpkm1                                 ! Horizontal slab
315                     ze2u = rotdiff (ji,jj,jk)
316                     ze1v = hdivdiff(ji,jj,jk)
317                     ! horizontal diffusive trends
318                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
319                           + ( hdivdiff(ji+1,jj,jk) - ze1v  ) / e1u(ji,jj)
320
321                     ! add it to the general momentum trends
322                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
323
324                  END DO
325               ENDIF
326
327            END DO
328         END DO
329
330         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
331
332         jmax = j2-1
333         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
334
335         DO jj = j1+1, jmax
336            DO ji = i1+1, i2   ! vector opt.
337
338               IF (.NOT. tabspongedone_v(ji,jj)) THEN
339                  DO jk = 1, jpkm1                                 ! Horizontal slab
340                     ze2u = rotdiff (ji,jj,jk)
341                     ze1v = hdivdiff(ji,jj,jk)
342
343                     ! horizontal diffusive trends
344                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
345                           + ( hdivdiff(ji,jj+1,jk) - ze1v  ) / e2v(ji,jj)
346
347                     ! add it to the general momentum trends
348                     va(ji,jj,jk) = va(ji,jj,jk) + zva
349                  END DO
350               ENDIF
351               !
352            END DO
353         END DO
354         !
355         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
356         !
357      ENDIF
358      !
359   END SUBROUTINE interpun_sponge
360
361
362   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
363      !!---------------------------------------------
364      !!   *** ROUTINE interpvn_sponge ***
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      LOGICAL, INTENT(in) :: before
369      INTEGER, INTENT(in) :: nb , ndir
370      !
371      INTEGER  ::   ji, jj, jk
372      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
373      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
374      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
375      INTEGER :: imax
376      !!---------------------------------------------
377
378      IF( before ) THEN
379         tabres = vn(i1:i2,j1:j2,:)
380      ELSE
381         !
382         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
383         !
384         DO jk = 1, jpkm1                                 ! Horizontal slab
385            !                                             ! ===============
386
387            !                                             ! --------
388            ! Horizontal divergence                       !   div
389            !                                             ! --------
390            DO jj = j1+1,j2
391               DO ji = i1,i2   ! vector opt.
392                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
393                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
394                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
395               END DO
396            END DO
397            DO jj = j1,j2
398               DO ji = i1,i2-1   ! vector opt.
399                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
400                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
401                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
402               END DO
403            END DO
404         END DO
405
406         !                                                ! ===============
407         !                                               
408
409         imax = i2-1
410         IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-3)
411
412         DO jj = j1+1, j2
413            DO ji = i1+1, imax   ! vector opt.
414               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
415                  DO jk = 1, jpkm1
416                     ua(ji,jj,jk) = ua(ji,jj,jk)                                                               &
417                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  &
418                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
419                  END DO
420               ENDIF
421            END DO
422         END DO
423         !
424         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
425         !
426         DO jj = j1+1, j2-1
427            DO ji = i1+1, i2-1   ! vector opt.
428               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
429                  DO jk = 1, jpkm1
430                     va(ji,jj,jk) = va(ji,jj,jk)                                                                  &
431                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
432                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
433                  END DO
434               ENDIF
435            END DO
436         END DO
437         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
438      ENDIF
439      !
440   END SUBROUTINE interpvn_sponge
441
442#else
443CONTAINS
444   SUBROUTINE agrif_opa_sponge_empty
445      !!---------------------------------------------
446      !!   *** ROUTINE agrif_OPA_sponge_empty ***
447      !!---------------------------------------------
448      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
449   END SUBROUTINE agrif_opa_sponge_empty
450#endif
451
452   !!======================================================================
453END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.