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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 7953

Last change on this file since 7953 was 7953, checked in by gm, 7 years ago

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