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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

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