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

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 5819

Last change on this file since 5819 was 5819, checked in by timgraham, 8 years ago

Deleted fcm keywords

File size: 16.2 KB
Line 
1#define SPONGE && define SPONGE_TOP
2
3MODULE agrif_opa_sponge
4#if defined key_agrif  && ! defined key_offline
5   USE par_oce
6   USE oce
7   USE dom_oce
8   USE in_out_manager
9   USE agrif_oce
10   USE wrk_nemo 
11   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
12
13   IMPLICIT NONE
14   PRIVATE
15
16   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
17   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
18
19   !! * Substitutions
20#  include "domzgr_substitute.h90"
21   !!----------------------------------------------------------------------
22   !! NEMO/NST 3.3 , NEMO Consortium (2010)
23   !! $Id$
24   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
25   !!----------------------------------------------------------------------
26
27CONTAINS
28
29   SUBROUTINE Agrif_Sponge_Tra
30      !!---------------------------------------------
31      !!   *** ROUTINE Agrif_Sponge_Tra ***
32      !!---------------------------------------------
33      !!
34      REAL(wp) :: timecoeff
35
36#if defined SPONGE
37      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
38
39      CALL Agrif_Sponge
40      Agrif_SpecialValue=0.
41      Agrif_UseSpecialValue = .TRUE.
42      tabspongedone_tsn = .FALSE.
43
44      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
45
46      Agrif_UseSpecialValue = .FALSE.
47#endif
48
49   END SUBROUTINE Agrif_Sponge_Tra
50
51   SUBROUTINE Agrif_Sponge_dyn
52      !!---------------------------------------------
53      !!   *** ROUTINE Agrif_Sponge_dyn ***
54      !!---------------------------------------------
55      !!
56      REAL(wp) :: timecoeff
57
58#if defined SPONGE
59      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
60
61      Agrif_SpecialValue=0.
62      Agrif_UseSpecialValue = ln_spc_dyn
63
64      tabspongedone_u = .FALSE.
65      tabspongedone_v = .FALSE.         
66      CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
67
68      tabspongedone_u = .FALSE.
69      tabspongedone_v = .FALSE.
70      CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
71
72      Agrif_UseSpecialValue = .FALSE.
73#endif
74
75   END SUBROUTINE Agrif_Sponge_dyn
76
77   SUBROUTINE Agrif_Sponge
78      !!---------------------------------------------
79      !!   *** ROUTINE  Agrif_Sponge ***
80      !!---------------------------------------------
81      INTEGER  :: ji,jj,jk
82      INTEGER  :: ispongearea, ilci, ilcj
83      LOGICAL  :: ll_spdone
84      REAL(wp) :: z1spongearea, zramp
85      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
86
87#if defined SPONGE || defined SPONGE_TOP
88      ll_spdone=.TRUE.
89      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
90         ! Define ramp from boundaries towards domain interior
91         ! at T-points
92         ! Store it in ztabramp
93         ll_spdone=.FALSE.
94
95         CALL wrk_alloc( jpi, jpj, ztabramp )
96
97         ispongearea  = 2 + nn_sponge_len * Agrif_irhox()
98         ilci = nlci - ispongearea
99         ilcj = nlcj - ispongearea 
100         z1spongearea = 1._wp / REAL( ispongearea - 2 )
101
102         ztabramp(:,:) = 0._wp
103
104         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
105            DO jj = 1, jpj
106               IF ( umask(2,jj,1) == 1._wp ) THEN
107                 DO ji = 2, ispongearea                 
108                    ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
109                 END DO
110               ENDIF
111            ENDDO
112         ENDIF
113
114         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
115            DO jj = 1, jpj
116               IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
117                  DO ji = ilci+1,nlci-1
118                     zramp = (ji - (ilci+1) ) * z1spongearea
119                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
120                  ENDDO
121               ENDIF
122            ENDDO
123         ENDIF
124
125         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
126            DO ji = 1, jpi
127               IF ( vmask(ji,2,1) == 1._wp ) THEN
128                  DO jj = 2, ispongearea
129                     zramp = ( ispongearea-jj ) * z1spongearea
130                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
131                  END DO
132               ENDIF
133            ENDDO
134         ENDIF
135
136         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
137            DO ji = 1, jpi
138               IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
139                  DO jj = ilcj+1,nlcj-1
140                     zramp = (jj - (ilcj+1) ) * z1spongearea
141                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
142                  END DO
143               ENDIF
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
160         CALL lbc_lnk( fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
161         CALL lbc_lnk( fsaht_spv, 'V', 1. )
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
177         CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
178         CALL lbc_lnk( fsahm_spf, 'F', 1. )
179         spongedoneU = .TRUE.
180      ENDIF
181      !
182      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
183      !
184#endif
185
186   END SUBROUTINE Agrif_Sponge
187
188   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
189      !!---------------------------------------------
190      !!   *** ROUTINE interptsn_sponge ***
191      !!---------------------------------------------
192      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
193      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
194      LOGICAL, INTENT(in) :: before
195
196
197      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
198      INTEGER  ::   iku, ikv
199      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
200      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
201      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
202      !
203      IF (before) THEN
204         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
205      ELSE   
206   
207         tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)   
208         DO jn = 1, jpts           
209            DO jk = 1, jpkm1
210               DO jj = j1,j2-1
211                  DO ji = i1,i2-1
212                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
213                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
214                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
215                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
216                  ENDDO
217               ENDDO
218
219               IF( ln_zps ) THEN      ! set gradient at partial step level
220                  DO jj = j1,j2-1
221                     DO ji = i1,i2-1
222                        ! last level
223                        iku = mbku(ji,jj)
224                        ikv = mbkv(ji,jj)
225                        IF( iku == jk ) THEN
226                           ztu(ji,jj,jk) = 0._wp
227                        ENDIF
228                        IF( ikv == jk ) THEN
229                           ztv(ji,jj,jk) = 0._wp
230                        ENDIF
231                     END DO
232                  END DO
233               ENDIF
234            ENDDO
235
236            DO jk = 1, jpkm1
237               DO jj = j1+1,j2-1
238                  DO ji = i1+1,i2-1
239
240                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
241                        zbtr = r1_e12t(ji,jj) / fse3t_n(ji,jj,jk)
242                        ! horizontal diffusive trends
243                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  )
244                        ! add it to the general tracer trends
245                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
246                     ENDIF
247
248                  ENDDO
249               ENDDO
250
251            ENDDO
252         ENDDO
253
254         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
255
256      ENDIF
257
258   END SUBROUTINE interptsn_sponge
259
260   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
261      !!---------------------------------------------
262      !!   *** ROUTINE interpun_sponge ***
263      !!---------------------------------------------   
264      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
265      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
266      LOGICAL, INTENT(in) :: before
267
268      INTEGER :: ji,jj,jk
269
270      ! sponge parameters
271      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
272      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
273      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
274      INTEGER :: jmax
275      !
276
277
278      IF (before) THEN
279         tabres = un(i1:i2,j1:j2,:)
280      ELSE
281
282         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
283
284         DO jk = 1, jpkm1                                 ! Horizontal slab
285            !                                             ! ===============
286
287            !                                             ! --------
288            ! Horizontal divergence                       !   div
289            !                                             ! --------
290            DO jj = j1,j2
291               DO ji = i1+1,i2   ! vector opt.
292                  zbtr = r1_e12t(ji,jj) / fse3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
293                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*fse3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
294                                     &   -e2u(ji-1,jj)*fse3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
295               END DO
296            END DO
297
298            DO jj = j1,j2-1
299               DO ji = i1,i2   ! vector opt.
300                  zbtr = r1_e12f(ji,jj) * fse3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
301                  rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
302                                       +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) & 
303                                    & ) * fmask(ji,jj,jk) * zbtr 
304               END DO
305            END DO
306         ENDDO
307
308         !
309
310
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) * fse3u_n(ji,jj,jk) )   &
321                           + ( hdivdiff(ji+1,jj,jk) - ze1v  ) / 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) * fse3v_n(ji,jj,jk) )   &
347                           + ( hdivdiff(ji,jj+1,jk) - ze1v  ) / 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
358         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
359
360      ENDIF
361
362
363   END SUBROUTINE interpun_sponge
364
365
366   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
367      !!---------------------------------------------
368      !!   *** ROUTINE interpvn_sponge ***
369      !!---------------------------------------------
370      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
371      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
372      LOGICAL, INTENT(in) :: before
373      INTEGER, INTENT(in) :: nb , ndir
374
375      INTEGER :: ji,jj,jk
376
377      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
378
379      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
380      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
381      INTEGER :: imax
382      !
383
384      IF (before) THEN
385         tabres = vn(i1:i2,j1:j2,:)
386      ELSE
387
388         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
389
390         DO jk = 1, jpkm1                                 ! Horizontal slab
391            !                                             ! ===============
392
393            !                                             ! --------
394            ! Horizontal divergence                       !   div
395            !                                             ! --------
396            DO jj = j1+1,j2
397               DO ji = i1,i2   ! vector opt.
398                  zbtr = r1_e12t(ji,jj) / fse3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
399                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
400                                     &  -e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
401               END DO
402            END DO
403            DO jj = j1,j2
404               DO ji = i1,i2-1   ! vector opt.
405                  zbtr = r1_e12f(ji,jj) * fse3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
406                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
407                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk) &
408                                    & ) * fmask(ji,jj,jk) * zbtr
409               END DO
410            END DO
411         ENDDO
412
413         !                                                ! ===============
414         !                                               
415
416         imax = i2-1
417         IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3)
418
419         DO jj = j1+1, j2
420            DO ji = i1+1, imax   ! vector opt.
421               IF (.NOT. tabspongedone_u(ji,jj)) THEN
422                  DO jk = 1, jpkm1                                 ! Horizontal slab
423                     ze2u = rotdiff (ji,jj,jk)
424                     ze1v = hdivdiff(ji,jj,jk)
425                     ! horizontal diffusive trends
426                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u_n(ji,jj,jk) ) + ( hdivdiff(ji+1,jj,jk) - ze1v) &
427                           / e1u(ji,jj)
428
429
430                     ! add it to the general momentum trends
431                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
432                  END DO
433
434               ENDIF
435            END DO
436         END DO
437
438         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
439
440         DO jj = j1+1, j2-1
441            DO ji = i1+1, i2-1   ! vector opt.
442               IF (.NOT. tabspongedone_v(ji,jj)) THEN
443                  DO jk = 1, jpkm1                                 ! Horizontal slab
444                     ze2u = rotdiff (ji,jj,jk)
445                     ze1v = hdivdiff(ji,jj,jk)
446                     ! horizontal diffusive trends
447
448                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v_n(ji,jj,jk) ) + ( hdivdiff(ji,jj+1,jk) - ze1v) &
449                           / e2v(ji,jj)
450
451                     ! add it to the general momentum trends
452                     va(ji,jj,jk) = va(ji,jj,jk) + zva
453                  END DO
454               ENDIF
455            END DO
456         END DO
457         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
458      ENDIF
459
460   END SUBROUTINE interpvn_sponge
461
462#else
463CONTAINS
464
465   SUBROUTINE agrif_opa_sponge_empty
466      !!---------------------------------------------
467      !!   *** ROUTINE agrif_OPA_sponge_empty ***
468      !!---------------------------------------------
469      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
470   END SUBROUTINE agrif_opa_sponge_empty
471#endif
472
473END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.