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 @ 8214

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

#1883 (HPC-09) - correction of minor issues

  • 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) ) * fmask(ji,jj,jk) * zbtr 
308               END DO
309            END DO
310         END DO
311         !
312         DO jj = j1+1, j2-1
313            DO ji = i1+1, i2-1   ! vector opt.
314
315               IF (.NOT. tabspongedone_u(ji,jj)) THEN
316                  DO jk = 1, jpkm1                                 ! Horizontal slab
317                     ze2u = rotdiff (ji,jj,jk)
318                     ze1v = hdivdiff(ji,jj,jk)
319                     ! horizontal diffusive trends
320                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
321                           + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
322
323                     ! add it to the general momentum trends
324                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
325
326                  END DO
327               ENDIF
328
329            END DO
330         END DO
331
332         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
333
334         jmax = j2-1
335         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
336
337         DO jj = j1+1, jmax
338            DO ji = i1+1, i2   ! vector opt.
339
340               IF (.NOT. tabspongedone_v(ji,jj)) THEN
341                  DO jk = 1, jpkm1                                 ! Horizontal slab
342                     ze2u = rotdiff (ji,jj,jk)
343                     ze1v = hdivdiff(ji,jj,jk)
344
345                     ! horizontal diffusive trends
346                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
347                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
348
349                     ! add it to the general momentum trends
350                     va(ji,jj,jk) = va(ji,jj,jk) + zva
351                  END DO
352               ENDIF
353               !
354            END DO
355         END DO
356         !
357         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
358         !
359      ENDIF
360      !
361   END SUBROUTINE interpun_sponge
362
363
364   SUBROUTINE interpvn_sponge( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
365      !!----------------------------------------------------------------------
366      !!                 *** ROUTINE interpvn_sponge ***
367      !!----------------------------------------------------------------------
368      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
369      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   tabres
370      LOGICAL                               , INTENT(in   ) ::   before
371      INTEGER                               , INTENT(in   ) ::   nb , ndir
372      !
373      INTEGER ::   ji, jj, jk
374      INTEGER ::   imax
375      REAL(wp)::   ze2u, ze1v, zua, zva, zbtr
376      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) ::   vbdiff, rotdiff, hdivdiff
377      !!----------------------------------------------------------------------
378
379      IF( before ) THEN
380         tabres = vn(i1:i2,j1:j2,:)
381      ELSE
382         !
383         vbdiff(i1:i2,j1:j2,:) = ( vb(i1:i2,j1:j2,:) - tabres(:,:,:) ) * vmask(i1:i2,j1:j2,:)
384         !
385         DO jk = 1, jpkm1                                 ! Horizontal slab
386            !                                             ! ===============
387
388            !                                             ! --------
389            ! Horizontal divergence                       !   div
390            !                                             ! --------
391            DO jj = j1+1,j2
392               DO ji = i1,i2   ! vector opt.
393                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
394                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
395                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
396               END DO
397            END DO
398            DO jj = j1,j2
399               DO ji = i1,i2-1   ! vector opt.
400                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
401                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
402                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
403               END DO
404            END DO
405         END DO
406
407         !                                                ! ===============
408         !                                               
409
410         imax = i2 - 1
411         IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-3)
412
413         DO jj = j1+1, j2
414            DO ji = i1+1, imax   ! vector opt.
415               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
416                  DO jk = 1, jpkm1
417                     ua(ji,jj,jk) = ua(ji,jj,jk)                                                               &
418                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  &
419                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
420                  END DO
421               ENDIF
422            END DO
423         END DO
424         !
425         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
426         !
427         DO jj = j1+1, j2-1
428            DO ji = i1+1, i2-1   ! vector opt.
429               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
430                  DO jk = 1, jpkm1
431                     va(ji,jj,jk) = va(ji,jj,jk)                                                                  &
432                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
433                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
434                  END DO
435               ENDIF
436            END DO
437         END DO
438         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
439      ENDIF
440      !
441   END SUBROUTINE interpvn_sponge
442
443#else
444   !!----------------------------------------------------------------------
445   !!   Empty module                                          no AGRIF zoom
446   !!----------------------------------------------------------------------
447CONTAINS
448   SUBROUTINE agrif_opa_sponge_empty
449      !!----------------------------------------------------------------------
450      !!                 *** ROUTINE agrif_OPA_sponge_empty ***
451      !!----------------------------------------------------------------------
452      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
453   END SUBROUTINE agrif_opa_sponge_empty
454#endif
455
456   !!======================================================================
457END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.