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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
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) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
213                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_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_e1e2t(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_e1e2t(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_e1e2f(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_e1e2t(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_e1e2f(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.