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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 8226

Last change on this file since 8226 was 8226, checked in by clem, 7 years ago

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