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_oce_sponge.F90 in NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90 @ 14048

Last change on this file since 14048 was 13937, checked in by jchanut, 4 years ago

#2222, added:

  • Divergence conserving interpolation for transport (Balsara's scheme)
  • Possibility to shift barotropic interface inside the domain (currently disabled - needed to match volume within sponge)
  • PPR library as default vertical interpolation scheme
  • Trilinear interpolation over partial bottom cells
  • Property svn:keywords set to Id
File size: 44.7 KB
Line 
1#define SPONGE && define SPONGE_TOP
2
3MODULE agrif_oce_sponge
4   !!======================================================================
5   !!                   ***  MODULE  agrif_oce_interp  ***
6   !! AGRIF: sponge package for the ocean dynamics (OPA)
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 lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE iom
25   USE vremap
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC Agrif_Sponge, Agrif_Sponge_2d, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
31   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
32   PUBLIC interpunb_sponge, interpvnb_sponge
33
34   !! * Substitutions
35#  include "do_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/NST 4.0 , NEMO Consortium (2018)
38   !! $Id$
39   !! Software governed by the CeCILL license (see ./LICENSE)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE Agrif_Sponge_Tra
44      !!----------------------------------------------------------------------
45      !!                 *** ROUTINE Agrif_Sponge_Tra ***
46      !!----------------------------------------------------------------------
47      REAL(wp) ::   zcoef   ! local scalar
48      INTEGER  :: istart, iend, jstart, jend 
49      !!----------------------------------------------------------------------
50      !
51#if defined SPONGE
52      !! Assume persistence:
53      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
54
55      Agrif_SpecialValue    = 0._wp
56      Agrif_UseSpecialValue = .TRUE.
57      l_vremap              = ln_vert_remap
58      tabspongedone_tsn     = .FALSE.
59      !
60      CALL Agrif_Bc_Variable( ts_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
61      !
62      Agrif_UseSpecialValue = .FALSE.
63      l_vremap              = .FALSE.
64#endif
65      !
66      CALL iom_put( 'agrif_spu', fspu(:,:))
67      CALL iom_put( 'agrif_spv', fspv(:,:))
68      !
69   END SUBROUTINE Agrif_Sponge_Tra
70
71
72   SUBROUTINE Agrif_Sponge_dyn
73      !!----------------------------------------------------------------------
74      !!                 *** ROUTINE Agrif_Sponge_dyn ***
75      !!----------------------------------------------------------------------
76      REAL(wp) ::   zcoef   ! local scalar
77      !!----------------------------------------------------------------------
78      !
79#if defined SPONGE
80      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
81
82      Agrif_SpecialValue    = 0._wp
83      Agrif_UseSpecialValue = ln_spc_dyn
84      l_vremap              = ln_vert_remap
85      use_sign_north        = .TRUE.
86      sign_north            = -1._wp
87      !
88      tabspongedone_u = .FALSE.
89      tabspongedone_v = .FALSE.         
90      CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
91      !
92      tabspongedone_u = .FALSE.
93      tabspongedone_v = .FALSE.
94      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
95     
96      IF ( nn_shift_bar>0 ) THEN ! then split sponge between 2d and 3d
97         zcoef = REAL(Agrif_NbStepint(),wp)/REAL(Agrif_rhot()) ! forward tsplit
98         tabspongedone_u = .FALSE.
99         tabspongedone_v = .FALSE.         
100         CALL Agrif_Bc_Variable( unb_sponge_id, calledweight=zcoef, procname=interpunb_sponge )
101         !
102         tabspongedone_u = .FALSE.
103         tabspongedone_v = .FALSE.
104         CALL Agrif_Bc_Variable( vnb_sponge_id, calledweight=zcoef, procname=interpvnb_sponge )
105      ENDIF
106      !
107      Agrif_UseSpecialValue = .FALSE.
108      use_sign_north        = .FALSE.
109      l_vremap              = .FALSE.
110      !
111#endif
112      !
113      CALL iom_put( 'agrif_spt', fspt(:,:))
114      CALL iom_put( 'agrif_spf', fspf(:,:))
115      !
116   END SUBROUTINE Agrif_Sponge_dyn
117
118
119   SUBROUTINE Agrif_Sponge
120      !!----------------------------------------------------------------------
121      !!                 *** ROUTINE  Agrif_Sponge ***
122      !!----------------------------------------------------------------------
123      INTEGER  ::   ji, jj, ind1, ind2
124      INTEGER  ::   ispongearea, jspongearea
125      REAL(wp) ::   z1_ispongearea, z1_jspongearea
126      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
127      !!----------------------------------------------------------------------
128      !
129      ! Sponge 1d example with:
130      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
131      !                       
132      !coarse :     U     T     U     T     U     T     U
133      !|            |           |           |           |
134      !fine :     t u t u t u t u t u t u t u t u t u t u t
135      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0
136      !           |   ghost     | <-- sponge area  -- > |
137      !           |   points    |                       |
138      !                         |--> dynamical interface
139
140#if defined SPONGE || defined SPONGE_TOP
141      ! Define ramp from boundaries towards domain interior at F-points
142      ! Store it in ztabramp
143
144      ispongearea    = nn_sponge_len * Agrif_irhox()
145      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp )
146      jspongearea    = nn_sponge_len * Agrif_irhoy()
147      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp )
148         
149      ztabramp(:,:) = 0._wp
150
151      IF( lk_west ) THEN                             ! --- West --- !
152         ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
153         ind2 = nn_hls + 1 + nbghostcells + ispongearea 
154         DO ji = mi0(ind1), mi1(ind2)   
155            DO jj = 1, jpj               
156               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea
157            END DO
158         END DO
159         ! ghost cells:
160         ind1 = 1
161         ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
162         DO ji = mi0(ind1), mi1(ind2)   
163            DO jj = 1, jpj               
164               ztabramp(ji,jj) = 1._wp
165            END DO
166         END DO
167      ENDIF
168      IF( lk_east ) THEN                             ! --- East --- !
169         ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1
170         ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1
171         DO ji = mi0(ind1), mi1(ind2)
172            DO jj = 1, jpj
173               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) 
174            END DO
175         END DO
176         ! ghost cells:
177         ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1
178         ind2 = jpiglo - 1
179         DO ji = mi0(ind1), mi1(ind2)
180            DO jj = 1, jpj
181               ztabramp(ji,jj) = 1._wp
182            END DO
183         END DO
184      ENDIF     
185      IF( lk_south ) THEN                            ! --- South --- !
186         ind1 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells
187         ind2 = nn_hls + 1 + nbghostcells + jspongearea 
188         DO jj = mj0(ind1), mj1(ind2) 
189            DO ji = 1, jpi
190               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
191            END DO
192         END DO
193         ! ghost cells:
194         ind1 = 1
195         ind2 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells
196         DO jj = mj0(ind1), mj1(ind2) 
197            DO ji = 1, jpi
198               ztabramp(ji,jj) = 1._wp
199            END DO
200         END DO
201      ENDIF
202      IF( lk_north ) THEN                            ! --- North --- !
203         ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1
204         ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1
205         DO jj = mj0(ind1), mj1(ind2)
206            DO ji = 1, jpi
207               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) 
208            END DO
209         END DO
210         ! ghost cells:
211         ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
212         ind2 = jpjglo
213         DO jj = mj0(ind1), mj1(ind2)
214            DO ji = 1, jpi
215               ztabramp(ji,jj) = 1._wp
216            END DO
217         END DO
218      ENDIF     
219      !
220      ! Tracers
221      fspu(:,:) = 0._wp
222      fspv(:,:) = 0._wp
223      DO_2D( 0, 0, 0, 0 )
224         fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj)
225         fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj)
226      END_2D
227
228      ! Dynamics
229      fspt(:,:) = 0._wp
230      fspf(:,:) = 0._wp
231      DO_2D( 0, 0, 0, 0 )
232         fspt(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) &
233                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
234         fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
235      END_2D
236     
237      CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp )
238      !
239      ! Remove vertical interpolation where not needed:
240      ! (A null value in mbkx arrays does the job)
241      WHERE (fspu(:,:) == 0._wp) mbku_parent(:,:) = 0
242      WHERE (fspv(:,:) == 0._wp) mbkv_parent(:,:) = 0
243      WHERE (fspt(:,:) == 0._wp) mbkt_parent(:,:) = 0
244      !
245#endif
246      !
247   END SUBROUTINE Agrif_Sponge
248
249
250   SUBROUTINE Agrif_Sponge_2d
251      !!----------------------------------------------------------------------
252      !!                 *** ROUTINE  Agrif_Sponge_2d ***
253      !!----------------------------------------------------------------------
254      INTEGER  ::   ji, jj, ind1, ind2, ishift, jshift
255      INTEGER  ::   ispongearea, jspongearea
256      REAL(wp) ::   z1_ispongearea, z1_jspongearea
257      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
258      !!----------------------------------------------------------------------
259      !
260      ! Sponge 1d example with:
261      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
262      !                       
263      !coarse :     U     T     U     T     U     T     U
264      !|            |           |           |           |
265      !fine :     t u t u t u t u t u t u t u t u t u t u t
266      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0
267      !           |   ghost     | <-- sponge area  -- > |
268      !           |   points    |                       |
269      !                         |--> dynamical interface
270
271#if defined SPONGE || defined SPONGE_TOP
272      ! Define ramp from boundaries towards domain interior at F-points
273      ! Store it in ztabramp
274
275      ispongearea    = nn_sponge_len * Agrif_irhox()
276      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp )
277      jspongearea    = nn_sponge_len * Agrif_irhoy()
278      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp )
279      ishift = nn_shift_bar * Agrif_irhox()
280      jshift = nn_shift_bar * Agrif_irhoy()
281       
282      ztabramp(:,:) = 0._wp
283
284      IF( lk_west ) THEN                             ! --- West --- !
285         ind1 = nn_hls + 1 + nbghostcells + ishift
286         ind2 = nn_hls + 1 + nbghostcells + ishift + ispongearea 
287         DO ji = mi0(ind1), mi1(ind2)   
288            DO jj = 1, jpj               
289               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea
290            END DO
291         END DO
292         ! ghost cells:
293         ind1 = 1
294         ind2 = nn_hls + 1 + nbghostcells + ishift               ! halo + land + nbghostcells
295         DO ji = mi0(ind1), mi1(ind2)   
296            DO jj = 1, jpj               
297               ztabramp(ji,jj) = 1._wp
298            END DO
299         END DO
300      ENDIF
301      IF( lk_east ) THEN                             ! --- East --- !
302         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - ispongearea - 1
303         ind2 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1
304         DO ji = mi0(ind1), mi1(ind2)
305            DO jj = 1, jpj
306               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) 
307            END DO
308         END DO
309         ! ghost cells:
310         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1
311         ind2 = jpiglo - 1
312         DO ji = mi0(ind1), mi1(ind2)
313            DO jj = 1, jpj
314               ztabramp(ji,jj) = 1._wp
315            END DO
316         END DO
317      ENDIF     
318      IF( lk_south ) THEN                            ! --- South --- !
319         ind1 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells
320         ind2 = nn_hls + 1 + nbghostcells + jshift + jspongearea 
321         DO jj = mj0(ind1), mj1(ind2) 
322            DO ji = 1, jpi
323               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
324            END DO
325         END DO
326         ! ghost cells:
327         ind1 = 1
328         ind2 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells
329         DO jj = mj0(ind1), mj1(ind2) 
330            DO ji = 1, jpi
331               ztabramp(ji,jj) = 1._wp
332            END DO
333         END DO
334      ENDIF
335      IF( lk_north ) THEN                            ! --- North --- !
336         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift) - jspongearea - 1
337         ind2 = jpjglo - ( nn_hls + nbghostcells + jshift) - 1    ! halo + land + nbghostcells - 1
338         DO jj = mj0(ind1), mj1(ind2)
339            DO ji = 1, jpi
340               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) 
341            END DO
342         END DO
343         ! ghost cells:
344         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift)      ! halo + land + nbghostcells - 1
345         ind2 = jpjglo
346         DO jj = mj0(ind1), mj1(ind2)
347            DO ji = 1, jpi
348               ztabramp(ji,jj) = 1._wp
349            END DO
350         END DO
351      ENDIF
352      !
353      ! Tracers
354      fspu_2d(:,:) = 0._wp
355      fspv_2d(:,:) = 0._wp
356      DO_2D( 0, 0, 0, 0 )
357         fspu_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj)
358         fspv_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj)
359      END_2D
360
361      ! Dynamics
362      fspt_2d(:,:) = 0._wp
363      fspf_2d(:,:) = 0._wp
364      DO_2D( 0, 0, 0, 0 )
365         fspt_2d(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) &
366                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
367         fspf_2d(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
368         END_2D
369      CALL lbc_lnk_multi( 'agrif_Sponge_2d', fspu_2d, 'U', 1._wp, fspv_2d, 'V', 1._wp, fspt_2d, 'T', 1._wp, fspf_2d, 'F', 1._wp )
370      !
371#endif
372      !
373   END SUBROUTINE Agrif_Sponge_2d
374
375
376   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
377      !!----------------------------------------------------------------------
378      !!                 *** ROUTINE interptsn_sponge ***
379      !!----------------------------------------------------------------------
380      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
381      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
382      LOGICAL                                     , INTENT(in   ) ::   before
383      !
384      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
385      INTEGER  ::   iku, ikv
386      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
387      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
388      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
389      ! vertical interpolation:
390      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
391      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i
392      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i
393      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
394      INTEGER :: N_in, N_out
395      !!----------------------------------------------------------------------
396      !
397      IF( before ) THEN
398         DO jn = 1, jpts
399            DO jk=k1,k2
400               DO jj=j1,j2
401                  DO ji=i1,i2
402                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
403                  END DO
404               END DO
405            END DO
406         END DO
407
408         IF ( l_vremap.OR.ln_zps ) THEN
409
410            ! Fill cell depths (i.e. gdept) to be interpolated
411            ! Warning: these are masked, hence extrapolated prior interpolation.
412            DO jj=j1,j2
413               DO ji=i1,i2
414                  tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a)
415                  DO jk=k1+1,k2
416                     tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
417                        & ( tabres(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) )
418                  END DO
419               END DO
420            END DO
421
422            ! Save ssh at last level:
423            IF ( .NOT.ln_linssh ) THEN
424               tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
425            END IF 
426   
427         END IF
428
429      ELSE 
430         !
431         IF ( l_vremap ) THEN
432
433            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
434
435            DO jj=j1,j2
436               DO ji=i1,i2
437
438                  tabres_child(ji,jj,:,:) = 0._wp 
439                  ! Build vertical grids:
440                  N_in = mbkt_parent(ji,jj)
441                  ! Input grid (account for partial cells if any):
442                  DO jk=1,N_in
443                     z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2)
444                     tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts)
445                  END DO
446                 
447                  ! Intermediate grid:
448                  DO jk = 1, N_in
449                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
450                       &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
451                  END DO
452                  z_in_i(1) = 0.5_wp * h_in_i(1)
453                  DO jk=2,N_in
454                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
455                  END DO
456                  z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2)
457
458                  ! Output (Child) grid:
459                  N_out = mbkt(ji,jj)
460                  DO jk=1,N_out
461                     h_out(jk) = e3t(ji,jj,jk,Kbb_a)
462                  END DO
463                  z_out(1) = 0.5_wp * h_out(1)
464                  DO jk=2,N_out
465                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
466                  END DO
467                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a)
468
469                  ! Account for small differences in the free-surface
470                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN
471                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) )
472                  ELSE
473                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) )
474                  END IF
475                  IF (N_in*N_out > 0) THEN
476                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),z_in_i(1:N_in),N_in,N_in,jpts)
477                     CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
478!                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),z_out(1:N_in),N_in,N_out,jpts) 
479                  ENDIF
480               END DO
481            END DO
482
483            DO jj=j1,j2
484               DO ji=i1,i2
485                  DO jk=1,jpkm1
486                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
487                  END DO
488               END DO
489            END DO
490
491         ELSE
492
493            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
494
495               DO jj=j1,j2
496                  DO ji=i1,i2
497                     !
498                     N_in  = mbkt(ji,jj) 
499                     N_out = mbkt(ji,jj) 
500                     z_in(1) = tabres(ji,jj,1,n2)
501                     tabin(1,1:jpts) = tabres(ji,jj,1,1:jpts)
502                     DO jk=2, N_in
503                        z_in(jk) = tabres(ji,jj,jk,n2)
504                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts)
505                     END DO
506                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2)
507
508                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a)
509                     DO jk=2, N_out
510                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a)) 
511                     END DO
512                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
513
514                     CALL remap_linear(tabin(1:N_in,1:jpts), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jpts), &
515                                         &   z_out(1:N_out), N_in, N_out, jpts)
516                  END DO
517               END DO
518            ENDIF
519
520            DO jj=j1,j2
521               DO ji=i1,i2
522                  DO jk=1,jpkm1
523                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
524                  END DO
525               END DO
526            END DO
527
528         END IF
529
530         DO jn = 1, jpts           
531            DO jk = 1, jpkm1
532               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
533               DO jj = j1,j2
534                  DO ji = i1,i2-1
535                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
536                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
537                  END DO
538               END DO
539               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
540               DO ji = i1,i2
541                  DO jj = j1,j2-1
542                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
543                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
544                  END DO
545               END DO
546               !
547               IF( ln_zps ) THEN      ! set gradient at partial step level
548                  DO jj = j1,j2
549                     DO ji = i1,i2
550                        ! last level
551                        iku = mbku(ji,jj)
552                        ikv = mbkv(ji,jj)
553                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
554                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
555                     END DO
556                  END DO
557               ENDIF
558            END DO
559            !
560! JC: there is something wrong with the Laplacian in corners
561            DO jk = 1, jpkm1
562               DO jj = j1,j2
563                  DO ji = i1,i2
564                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
565                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
566                        ! horizontal diffusive trends
567                        ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
568                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
569                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
570                        ! add it to the general tracer trends
571                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
572                     ENDIF
573                  END DO
574               END DO
575
576            END DO
577            !
578         END DO
579         !
580         tabspongedone_tsn(i1:i2,j1:j2) = .TRUE.
581         !
582      ENDIF
583      !
584   END SUBROUTINE interptsn_sponge
585
586   
587   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
588      !!---------------------------------------------
589      !!   *** ROUTINE interpun_sponge ***
590      !!---------------------------------------------   
591      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
592      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
593      LOGICAL, INTENT(in) :: before
594
595      INTEGER  :: ji,jj,jk,jmax
596      INTEGER  :: ind1
597      ! sponge parameters
598      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
599      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
600      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
601      ! vertical interpolation:
602      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
603      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
604      REAL(wp), DIMENSION(1:jpk) :: h_out
605      INTEGER ::N_in, N_out
606      !!---------------------------------------------   
607      !
608      IF( before ) THEN
609         DO jk=k1,k2
610            DO jj=j1,j2
611               DO ji=i1,i2
612                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
613               END DO
614            END DO
615         END DO
616
617         IF ( l_vremap ) THEN
618
619            DO jk=k1,k2
620               DO jj=j1,j2
621                  DO ji=i1,i2
622                     tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
623                  END DO
624               END DO
625            END DO
626
627            ! Extrapolate thicknesses in partial bottom cells:
628            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
629            IF (ln_zps) THEN
630               DO jj=j1,j2
631                  DO ji=i1,i2
632                     jk = mbku(ji,jj)
633                     tabres(ji,jj,jk,m2) = 0._wp
634                  END DO
635               END DO           
636            END IF
637            ! Save ssh at last level:
638            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
639            IF (.NOT.ln_linssh) THEN
640               ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
641               DO jk=1,jpk
642                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
643               END DO
644               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
645            END IF
646         END IF
647
648      ELSE
649
650         IF ( l_vremap ) THEN
651
652            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
653
654            DO jj=j1,j2
655               DO ji=i1,i2
656                  tabres_child(ji,jj,:) = 0._wp
657                  N_in = mbku_parent(ji,jj)
658                  zhtot = 0._wp
659                  DO jk=1,N_in
660                     !IF (jk==N_in) THEN
661                     !   h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
662                     !ELSE
663                     !   h_in(jk) = tabres(ji,jj,jk,m2)
664                     !ENDIF
665                     h_in(jk) = e3u0_parent(ji,jj,jk)
666                     zhtot = zhtot + h_in(jk)
667                     tabin(jk) = tabres(ji,jj,jk,m1)
668                  END DO
669                  !         
670                  N_out = 0
671                  DO jk=1,jpk
672                     IF (umask(ji,jj,jk) == 0) EXIT
673                     N_out = N_out + 1
674                     h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
675                  END DO
676
677                  ! Account for small differences in free-surface
678                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
679                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
680                  ELSE
681                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
682                  ENDIF
683                 
684                  IF (N_in * N_out > 0) THEN
685                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
686                  ENDIF
687               END DO
688            END DO
689
690            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
691         ELSE
692
693            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
694 
695         ENDIF
696         !
697         DO jk = 1, jpkm1                                 ! Horizontal slab
698            !                                             ! ===============
699
700            !                                             ! --------
701            ! Horizontal divergence                       !   div
702            !                                             ! --------
703            DO jj = j1,j2
704               DO ji = i1+1,i2   ! vector opt.
705                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
706                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kbb_a) * ubdiff(ji  ,jj,jk) &
707                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr
708               END DO
709            END DO
710
711            DO jj = j1,j2-1
712               DO ji = i1,i2   ! vector opt.
713                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
714                  rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
715                                    &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
716               END DO
717            END DO
718         END DO
719         !
720         DO jj = j1+1, j2-1
721            DO ji = i1+1, i2-1   ! vector opt.
722
723               IF (.NOT. tabspongedone_u(ji,jj)) THEN
724                  DO jk = 1, jpkm1                                 ! Horizontal slab
725                     ze2u = rotdiff (ji,jj,jk)
726                     ze1v = hdivdiff(ji,jj,jk)
727                     ! horizontal diffusive trends
728                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )   &
729                         & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 
730                         & - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk)
731
732                     ! add it to the general momentum trends
733                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua                                 
734                  END DO
735               ENDIF
736
737            END DO
738         END DO
739
740         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
741
742         jmax = j2-1
743         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
744         DO jj = mj0(ind1), mj1(ind1)                 
745            jmax = MIN(jmax,jj)
746         END DO
747
748         DO jj = j1+1, jmax
749            DO ji = i1+1, i2   ! vector opt.
750
751               IF (.NOT. tabspongedone_v(ji,jj)) THEN
752                  DO jk = 1, jpkm1                                 ! Horizontal slab
753                     ze2u = rotdiff (ji,jj,jk)
754                     ze1v = hdivdiff(ji,jj,jk)
755
756                     ! horizontal diffusive trends
757                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
758                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
759
760                     ! add it to the general momentum trends
761                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
762                  END DO
763               ENDIF
764               !
765            END DO
766         END DO
767         !
768         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
769         !
770      ENDIF
771      !
772   END SUBROUTINE interpun_sponge
773
774   
775   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
776      !!---------------------------------------------
777      !!   *** ROUTINE interpvn_sponge ***
778      !!---------------------------------------------
779      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
780      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
781      LOGICAL, INTENT(in) :: before
782      !
783      INTEGER  ::   ji, jj, jk, imax
784      INTEGER  :: ind1
785      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot
786      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
787      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
788      ! vertical interpolation:
789      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
790      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
791      REAL(wp), DIMENSION(1:jpk) :: h_out
792      INTEGER :: N_in, N_out
793      !!---------------------------------------------
794     
795      IF( before ) THEN
796         DO jk=k1,k2
797            DO jj=j1,j2
798               DO ji=i1,i2
799                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
800               END DO
801            END DO
802         END DO
803
804         IF ( l_vremap ) THEN
805
806            DO jk=k1,k2
807               DO jj=j1,j2
808                  DO ji=i1,i2
809                     tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
810                  END DO
811               END DO
812            END DO
813            ! Extrapolate thicknesses in partial bottom cells:
814            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
815            IF (ln_zps) THEN
816               DO jj=j1,j2
817                  DO ji=i1,i2
818                     jk = mbkv(ji,jj)
819                     tabres(ji,jj,jk,m2) = 0._wp
820                  END DO
821               END DO           
822            END IF
823            ! Save ssh at last level:
824            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
825            IF (.NOT.ln_linssh) THEN
826               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
827               DO jk=1,jpk
828                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
829               END DO
830               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
831            END IF
832
833         END IF
834
835      ELSE
836
837         IF ( l_vremap ) THEN
838            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
839            DO jj=j1,j2
840               DO ji=i1,i2
841                  tabres_child(ji,jj,:) = 0._wp
842                  N_in = mbkv_parent(ji,jj)
843                  zhtot = 0._wp
844                  DO jk=1,N_in
845                     !IF (jk==N_in) THEN
846                     !   h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
847                     !ELSE
848                     !   h_in(jk) = tabres(ji,jj,jk,m2)
849                     !ENDIF
850                     h_in(jk) = e3v0_parent(ji,jj,jk)
851                     zhtot = zhtot + h_in(jk)
852                     tabin(jk) = tabres(ji,jj,jk,m1)
853                  END DO
854                  !         
855                  N_out = 0
856                  DO jk=1,jpk
857                     IF (vmask(ji,jj,jk) == 0) EXIT
858                     N_out = N_out + 1
859                     h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
860                  END DO
861
862                  ! Account for small differences in free-surface
863                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
864                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
865                  ELSE
866                     h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
867                  ENDIF
868         
869                  IF (N_in * N_out > 0) THEN
870                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
871                  ENDIF
872               END DO
873            END DO
874
875            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 
876         ELSE
877
878            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
879
880         ENDIF
881         !
882         DO jk = 1, jpkm1                                 ! Horizontal slab
883            !                                             ! ===============
884
885            !                                             ! --------
886            ! Horizontal divergence                       !   div
887            !                                             ! --------
888            DO jj = j1+1,j2
889               DO ji = i1,i2   ! vector opt.
890                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
891                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kbb_a) * vbdiff(ji,jj  ,jk)  &
892                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk)  ) * zbtr
893               END DO
894            END DO
895            DO jj = j1,j2
896               DO ji = i1,i2-1   ! vector opt.
897                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
898                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
899                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
900               END DO
901            END DO
902         END DO
903
904         !                                                ! ===============
905         !                                               
906
907         imax = i2 - 1
908         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East
909         DO ji = mi0(ind1), mi1(ind1)               
910            imax = MIN(imax,ji)
911         END DO
912         
913         DO jj = j1+1, j2
914            DO ji = i1+1, imax   ! vector opt.
915               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
916                  DO jk = 1, jpkm1
917                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                     &
918                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
919                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
920                  END DO
921               ENDIF
922            END DO
923         END DO
924         !
925         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
926         !
927         DO jj = j1+1, j2-1
928            DO ji = i1+1, i2-1   ! vector opt.
929               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
930                  DO jk = 1, jpkm1
931                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                        &
932                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
933                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                          &
934                        &  - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
935                  END DO
936               ENDIF
937            END DO
938         END DO
939         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
940      ENDIF
941      !
942   END SUBROUTINE interpvn_sponge
943
944   SUBROUTINE interpunb_sponge(tabres,i1,i2,j1,j2, before)
945      !!---------------------------------------------
946      !!   *** ROUTINE interpunb_sponge ***
947      !!---------------------------------------------   
948      INTEGER, INTENT(in) :: i1,i2,j1,j2
949      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
950      LOGICAL, INTENT(in) :: before
951
952      INTEGER  :: ji, jj, ind1, jmax
953      ! sponge parameters
954      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
955      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ubdiff
956      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff
957      !!---------------------------------------------   
958      !
959      IF( before ) THEN
960         DO jj=j1,j2
961            DO ji=i1,i2
962               tabres(ji,jj) = uu_b(ji,jj,Kmm_a)
963            END DO
964         END DO
965
966      ELSE
967
968         ubdiff(i1:i2,j1:j2) = (uu_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*umask(i1:i2,j1:j2,1)
969         !
970         !                                             ! --------
971         ! Horizontal divergence                       !   div
972         !                                             ! --------
973         DO jj = j1,j2
974            DO ji = i1+1,i2   ! vector opt.
975               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj)
976               hdivdiff(ji,jj) = (  e2u(ji  ,jj)*hu(ji  ,jj,Kbb_a) * ubdiff(ji  ,jj) &
977                                  &-e2u(ji-1,jj)*hu(ji-1,jj,Kbb_a) * ubdiff(ji-1,jj) ) * zbtr
978            END DO
979         END DO
980
981         DO jj = j1,j2-1
982            DO ji = i1,i2   ! vector opt.
983               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)
984               rotdiff(ji,jj) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1)   &
985                              &   +e1u(ji,jj  ) * ubdiff(ji,jj  ) ) * fmask(ji,jj,1) * zbtr 
986            END DO
987         END DO
988         !
989         DO jj = j1+1, j2-1
990            DO ji = i1+1, i2-1   ! vector opt.
991               IF (.NOT. tabspongedone_u(ji,jj)) THEN
992                  ze2u = rotdiff (ji,jj)
993                  ze1v = hdivdiff(ji,jj)
994                  ! horizontal diffusive trends
995                  zua = - ( ze2u - rotdiff (ji,jj-1) ) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  &
996                      & + ( hdivdiff(ji+1,jj) - ze1v ) * r1_e1u(ji,jj)                       & 
997                      & - rn_trelax_dyn * r1_Dt * fspu_2d(ji,jj) * ubdiff(ji,jj)
998
999                  ! add it to the general momentum trends
1000                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua                                 
1001               ENDIF
1002            END DO
1003         END DO
1004
1005         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
1006
1007         jmax = j2-1
1008         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
1009         DO jj = mj0(ind1), mj1(ind1)                 
1010            jmax = MIN(jmax,jj)
1011         END DO
1012
1013         DO jj = j1+1, jmax
1014            DO ji = i1+1, i2   ! vector opt.
1015               IF (.NOT. tabspongedone_v(ji,jj)) THEN
1016                     ze2u = rotdiff (ji,jj)
1017                     ze1v = hdivdiff(ji,jj)
1018                     zva = + ( ze2u - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) * r1_hv(ji,jj,Kmm_a) &
1019                           + ( hdivdiff(ji,jj+1) - ze1v ) * r1_e2v(ji,jj)
1020                     vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva
1021               ENDIF
1022            END DO
1023         END DO
1024         !
1025         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
1026         !
1027      ENDIF
1028      !
1029   END SUBROUTINE interpunb_sponge
1030
1031   
1032   SUBROUTINE interpvnb_sponge(tabres,i1,i2,j1,j2, before)
1033      !!---------------------------------------------
1034      !!   *** ROUTINE interpvnb_sponge ***
1035      !!---------------------------------------------
1036      INTEGER, INTENT(in) :: i1,i2,j1,j2
1037      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1038      LOGICAL, INTENT(in) :: before
1039      !
1040      INTEGER  ::   ji, jj, ind1, imax
1041      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
1042      REAL(wp), DIMENSION(i1:i2,j1:j2) :: vbdiff
1043      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff
1044      !!---------------------------------------------
1045     
1046      IF( before ) THEN
1047         DO jj=j1,j2
1048            DO ji=i1,i2
1049               tabres(ji,jj) = vv_b(ji,jj,Kmm_a)
1050            END DO
1051         END DO
1052      ELSE
1053         vbdiff(i1:i2,j1:j2) = (vv_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*vmask(i1:i2,j1:j2,1)
1054         !                                             ! --------
1055         ! Horizontal divergence                       !   div
1056         !                                             ! --------
1057         DO jj = j1+1,j2
1058            DO ji = i1,i2   ! vector opt.
1059               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj)
1060               hdivdiff(ji,jj) = ( e1v(ji,jj  ) * hv(ji,jj  ,Kbb_a) * vbdiff(ji,jj  )  &
1061                               &  -e1v(ji,jj-1) * hv(ji,jj-1,Kbb_a) * vbdiff(ji,jj-1)  ) * zbtr
1062            END DO
1063         END DO
1064         DO jj = j1,j2
1065            DO ji = i1,i2-1   ! vector opt.
1066               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj) 
1067               rotdiff(ji,jj) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj) & 
1068                              &  -e2v(ji  ,jj) * vbdiff(ji  ,jj)  ) * fmask(ji,jj,1) * zbtr
1069            END DO
1070         END DO
1071         !                                                ! ===============
1072         !                                               
1073
1074         imax = i2 - 1
1075         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East
1076         DO ji = mi0(ind1), mi1(ind1)               
1077            imax = MIN(imax,ji)
1078         END DO
1079         
1080         DO jj = j1+1, j2
1081            DO ji = i1+1, imax   ! vector opt.
1082               IF( .NOT. tabspongedone_u(ji,jj) ) THEN                                                     
1083                  zua = - ( rotdiff (ji  ,jj) - rotdiff (ji,jj-1)) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  &
1084                      & + ( hdivdiff(ji+1,jj) - hdivdiff(ji,jj  )) * r1_e1u(ji,jj)
1085                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua
1086               ENDIF
1087            END DO
1088         END DO
1089         !
1090         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
1091         !
1092         DO jj = j1+1, j2-1
1093            DO ji = i1+1, i2-1   ! vector opt.
1094               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
1095                  zva  =  ( rotdiff (ji,jj  ) - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) *r1_hv(ji,jj,Kmm_a) &
1096                     &  + ( hdivdiff(ji,jj+1) - hdivdiff(ji  ,jj) ) * r1_e2v(ji,jj)                     &
1097                     &  - rn_trelax_dyn * r1_Dt * fspv_2d(ji,jj) * vbdiff(ji,jj)
1098                  vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva
1099               ENDIF
1100            END DO
1101         END DO
1102         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
1103      ENDIF
1104      !
1105   END SUBROUTINE interpvnb_sponge
1106
1107
1108#else
1109   !!----------------------------------------------------------------------
1110   !!   Empty module                                          no AGRIF zoom
1111   !!----------------------------------------------------------------------
1112CONTAINS
1113   SUBROUTINE agrif_oce_sponge_empty
1114      WRITE(*,*)  'agrif_oce_sponge : You should not have seen this print! error?'
1115   END SUBROUTINE agrif_oce_sponge_empty
1116#endif
1117
1118   !!======================================================================
1119END MODULE agrif_oce_sponge
Note: See TracBrowser for help on using the repository browser.