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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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