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/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_oce_sponge.F90 @ 15437

Last change on this file since 15437 was 15437, checked in by jchanut, 12 months ago

agrif fixes: i) Add masking prior interpolation on "before" tracers: these seem to be strangly unmasked which induces extrapolation issues. ii) Add flux limiter on tracer sponge (recall sponge is applied on anomalies) to prevent from generating spurious extrema.

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