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 @ 14201

Last change on this file since 14201 was 14170, checked in by jchanut, 3 years ago

#2222, 2129: 1) Corrected ssh initialization from parent in line with what has been introduced by Sibylle 2) Fixed bug in dyn interp with expliciit free surface 3) Added check on number of levels in child grid without vertical remapping (must be < jpk_parent) 4) Removed the constrain on initialization from parent only when starting from climatology (requires Euler first step though).
  • Property svn:keywords set to Id
File size: 44.9 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 "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 + 1 + nbghostcells               ! halo + land + nbghostcells
154         ind2 = nn_hls + 1 + 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 + 1 + nbghostcells               ! halo + land + 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 ) - ispongearea - 1
171         ind2 = jpiglo - ( nn_hls + nbghostcells ) - 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    ! 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 + 1 + nbghostcells                 ! halo + land + nbghostcells
188         ind2 = nn_hls + 1 + 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 + 1 + nbghostcells                 ! halo + land + 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 ) - jspongearea - 1
205         ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + 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 )      ! 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_multi( '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 + 1 + nbghostcells + ishift
287         ind2 = nn_hls + 1 + 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 + 1 + nbghostcells + ishift               ! halo + land + 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 + ishift) - ispongearea - 1
304         ind2 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + 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 + ishift) - 1    ! halo + land + 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 + 1 + nbghostcells + jshift                ! halo + land + nbghostcells
321         ind2 = nn_hls + 1 + 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 + 1 + 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 + jshift) - jspongearea - 1
338         ind2 = jpjglo - ( nn_hls + nbghostcells + 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 + 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_multi( '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), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
389      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
390      ! vertical interpolation:
391      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
392      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i
393      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i
394      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
395      INTEGER :: N_in, N_out
396      !!----------------------------------------------------------------------
397      !
398      IF( before ) THEN
399         DO jn = 1, jpts
400            DO jk=k1,k2
401               DO jj=j1,j2
402                  DO ji=i1,i2
403                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
404                  END DO
405               END DO
406            END DO
407         END DO
408
409         IF ( l_vremap.OR.ln_zps ) THEN
410
411            ! Fill cell depths (i.e. gdept) to be interpolated
412            ! Warning: these are masked, hence extrapolated prior interpolation.
413            DO jj=j1,j2
414               DO ji=i1,i2
415                  tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a)
416                  DO jk=k1+1,k2
417                     tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
418                        & ( tabres(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) )
419                  END DO
420               END DO
421            END DO
422
423            ! Save ssh at last level:
424            IF ( .NOT.ln_linssh ) THEN
425               tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
426            END IF 
427   
428         END IF
429
430      ELSE 
431         !
432         IF ( l_vremap ) THEN
433
434            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
435
436            DO jj=j1,j2
437               DO ji=i1,i2
438
439                  tabres_child(ji,jj,:,:) = 0._wp 
440                  ! Build vertical grids:
441                  N_in = mbkt_parent(ji,jj)
442                  ! Input grid (account for partial cells if any):
443                  IF ( N_in > 0 ) THEN
444                     DO jk=1,N_in
445                        z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2)
446                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts)
447                     END DO
448                 
449                     ! Intermediate grid:
450                     DO jk = 1, N_in
451                        h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
452                          &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
453                     END DO
454                     z_in_i(1) = 0.5_wp * h_in_i(1)
455                     DO jk=2,N_in
456                        z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
457                     END DO
458                     z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2)
459                  END IF
460                  ! Output (Child) grid:
461                  N_out = mbkt(ji,jj)
462                  DO jk=1,N_out
463                     h_out(jk) = e3t(ji,jj,jk,Kbb_a)
464                  END DO
465                  z_out(1) = 0.5_wp * h_out(1)
466                  DO jk=2,N_out
467                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
468                  END DO
469                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a)
470
471                  ! Account for small differences in the free-surface
472                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN
473                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) )
474                  ELSE
475                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) )
476                  END IF
477                  IF (N_in*N_out > 0) THEN
478                     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)
479                     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)
480!                     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) 
481                  ENDIF
482               END DO
483            END DO
484
485            DO jj=j1,j2
486               DO ji=i1,i2
487                  DO jk=1,jpkm1
488                     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)
489                  END DO
490               END DO
491            END DO
492
493         ELSE
494
495            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
496
497               DO jj=j1,j2
498                  DO ji=i1,i2
499                     !
500                     N_in  = mbkt(ji,jj) 
501                     N_out = mbkt(ji,jj) 
502                     z_in(1) = tabres(ji,jj,1,n2)
503                     tabin(1,1:jpts) = tabres(ji,jj,1,1:jpts)
504                     DO jk=2, N_in
505                        z_in(jk) = tabres(ji,jj,jk,n2)
506                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts)
507                     END DO
508                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2)
509
510                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a)
511                     DO jk=2, N_out
512                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a)) 
513                     END DO
514                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
515
516                     CALL remap_linear(tabin(1:N_in,1:jpts), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jpts), &
517                                         &   z_out(1:N_out), N_in, N_out, jpts)
518                  END DO
519               END DO
520            ENDIF
521
522            DO jj=j1,j2
523               DO ji=i1,i2
524                  DO jk=1,jpkm1
525                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
526                  END DO
527               END DO
528            END DO
529
530         END IF
531
532         DO jn = 1, jpts           
533            DO jk = 1, jpkm1
534               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
535               DO jj = j1,j2
536                  DO ji = i1,i2-1
537                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
538                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
539                  END DO
540               END DO
541               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
542               DO ji = i1,i2
543                  DO jj = j1,j2-1
544                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
545                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
546                  END DO
547               END DO
548               !
549               IF( ln_zps ) THEN      ! set gradient at partial step level
550                  DO jj = j1,j2
551                     DO ji = i1,i2
552                        ! last level
553                        iku = mbku(ji,jj)
554                        ikv = mbkv(ji,jj)
555                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
556                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
557                     END DO
558                  END DO
559               ENDIF
560            END DO
561            !
562! JC: there is something wrong with the Laplacian in corners
563            DO jk = 1, jpkm1
564               DO jj = j1,j2
565                  DO ji = i1,i2
566                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
567                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
568                        ! horizontal diffusive trends
569                        ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
570                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
571                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
572                        ! add it to the general tracer trends
573                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
574                     ENDIF
575                  END DO
576               END DO
577
578            END DO
579            !
580         END DO
581         !
582         tabspongedone_tsn(i1:i2,j1:j2) = .TRUE.
583         !
584      ENDIF
585      !
586   END SUBROUTINE interptsn_sponge
587
588   
589   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
590      !!---------------------------------------------
591      !!   *** ROUTINE interpun_sponge ***
592      !!---------------------------------------------   
593      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
594      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
595      LOGICAL, INTENT(in) :: before
596
597      INTEGER  :: ji,jj,jk,jmax
598      INTEGER  :: ind1
599      ! sponge parameters
600      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
601      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
602      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
603      ! vertical interpolation:
604      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
605      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
606      REAL(wp), DIMENSION(1:jpk) :: h_out
607      INTEGER ::N_in, N_out
608      !!---------------------------------------------   
609      !
610      IF( before ) THEN
611         DO jk=k1,k2
612            DO jj=j1,j2
613               DO ji=i1,i2
614                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
615               END DO
616            END DO
617         END DO
618
619         IF ( l_vremap ) THEN
620
621            DO jk=k1,k2
622               DO jj=j1,j2
623                  DO ji=i1,i2
624                     tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
625                  END DO
626               END DO
627            END DO
628
629            ! Extrapolate thicknesses in partial bottom cells:
630            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
631            IF (ln_zps) THEN
632               DO jj=j1,j2
633                  DO ji=i1,i2
634                     jk = mbku(ji,jj)
635                     tabres(ji,jj,jk,m2) = 0._wp
636                  END DO
637               END DO           
638            END IF
639            ! Save ssh at last level:
640            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
641            IF (.NOT.ln_linssh) THEN
642               ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
643               DO jk=1,jpk
644                  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)
645               END DO
646               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
647            END IF
648         END IF
649
650      ELSE
651
652         IF ( l_vremap ) THEN
653
654            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
655
656            DO jj=j1,j2
657               DO ji=i1,i2
658                  tabres_child(ji,jj,:) = 0._wp
659                  N_in = mbku_parent(ji,jj)
660                  zhtot = 0._wp
661                  DO jk=1,N_in
662                     !IF (jk==N_in) THEN
663                     !   h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
664                     !ELSE
665                     !   h_in(jk) = tabres(ji,jj,jk,m2)
666                     !ENDIF
667                     h_in(jk) = e3u0_parent(ji,jj,jk)
668                     zhtot = zhtot + h_in(jk)
669                     tabin(jk) = tabres(ji,jj,jk,m1)
670                  END DO
671                  !         
672                  N_out = 0
673                  DO jk=1,jpk
674                     IF (umask(ji,jj,jk) == 0) EXIT
675                     N_out = N_out + 1
676                     h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
677                  END DO
678
679                  ! Account for small differences in free-surface
680                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
681                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
682                  ELSE
683                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
684                  ENDIF
685                 
686                  IF (N_in * N_out > 0) THEN
687                     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)
688                  ENDIF
689               END DO
690            END DO
691
692            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
693         ELSE
694
695            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
696 
697         ENDIF
698         !
699         DO jk = 1, jpkm1                                 ! Horizontal slab
700            !                                             ! ===============
701
702            !                                             ! --------
703            ! Horizontal divergence                       !   div
704            !                                             ! --------
705            DO jj = j1,j2
706               DO ji = i1+1,i2   ! vector opt.
707                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
708                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kbb_a) * ubdiff(ji  ,jj,jk) &
709                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr
710               END DO
711            END DO
712
713            DO jj = j1,j2-1
714               DO ji = i1,i2   ! vector opt.
715                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
716                  rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
717                                    &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
718               END DO
719            END DO
720         END DO
721         !
722         DO jj = j1+1, j2-1
723            DO ji = i1+1, i2-1   ! vector opt.
724
725               IF (.NOT. tabspongedone_u(ji,jj)) THEN
726                  DO jk = 1, jpkm1                                 ! Horizontal slab
727                     ze2u = rotdiff (ji,jj,jk)
728                     ze1v = hdivdiff(ji,jj,jk)
729                     ! horizontal diffusive trends
730                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )   &
731                         & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 
732                         & - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk)
733
734                     ! add it to the general momentum trends
735                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua                                 
736                  END DO
737               ENDIF
738
739            END DO
740         END DO
741
742         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
743
744         jmax = j2-1
745         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
746         DO jj = mj0(ind1), mj1(ind1)                 
747            jmax = MIN(jmax,jj)
748         END DO
749
750         DO jj = j1+1, jmax
751            DO ji = i1+1, i2   ! vector opt.
752
753               IF (.NOT. tabspongedone_v(ji,jj)) THEN
754                  DO jk = 1, jpkm1                                 ! Horizontal slab
755                     ze2u = rotdiff (ji,jj,jk)
756                     ze1v = hdivdiff(ji,jj,jk)
757
758                     ! horizontal diffusive trends
759                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
760                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
761
762                     ! add it to the general momentum trends
763                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
764                  END DO
765               ENDIF
766               !
767            END DO
768         END DO
769         !
770         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
771         !
772      ENDIF
773      !
774   END SUBROUTINE interpun_sponge
775
776   
777   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
778      !!---------------------------------------------
779      !!   *** ROUTINE interpvn_sponge ***
780      !!---------------------------------------------
781      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
782      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
783      LOGICAL, INTENT(in) :: before
784      !
785      INTEGER  ::   ji, jj, jk, imax
786      INTEGER  :: ind1
787      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot
788      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
789      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
790      ! vertical interpolation:
791      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
792      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
793      REAL(wp), DIMENSION(1:jpk) :: h_out
794      INTEGER :: N_in, N_out
795      !!---------------------------------------------
796     
797      IF( before ) THEN
798         DO jk=k1,k2
799            DO jj=j1,j2
800               DO ji=i1,i2
801                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
802               END DO
803            END DO
804         END DO
805
806         IF ( l_vremap ) THEN
807
808            DO jk=k1,k2
809               DO jj=j1,j2
810                  DO ji=i1,i2
811                     tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
812                  END DO
813               END DO
814            END DO
815            ! Extrapolate thicknesses in partial bottom cells:
816            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
817            IF (ln_zps) THEN
818               DO jj=j1,j2
819                  DO ji=i1,i2
820                     jk = mbkv(ji,jj)
821                     tabres(ji,jj,jk,m2) = 0._wp
822                  END DO
823               END DO           
824            END IF
825            ! Save ssh at last level:
826            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
827            IF (.NOT.ln_linssh) THEN
828               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
829               DO jk=1,jpk
830                  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)
831               END DO
832               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
833            END IF
834
835         END IF
836
837      ELSE
838
839         IF ( l_vremap ) THEN
840            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
841            DO jj=j1,j2
842               DO ji=i1,i2
843                  tabres_child(ji,jj,:) = 0._wp
844                  N_in = mbkv_parent(ji,jj)
845                  zhtot = 0._wp
846                  DO jk=1,N_in
847                     !IF (jk==N_in) THEN
848                     !   h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
849                     !ELSE
850                     !   h_in(jk) = tabres(ji,jj,jk,m2)
851                     !ENDIF
852                     h_in(jk) = e3v0_parent(ji,jj,jk)
853                     zhtot = zhtot + h_in(jk)
854                     tabin(jk) = tabres(ji,jj,jk,m1)
855                  END DO
856                  !         
857                  N_out = 0
858                  DO jk=1,jpk
859                     IF (vmask(ji,jj,jk) == 0) EXIT
860                     N_out = N_out + 1
861                     h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
862                  END DO
863
864                  ! Account for small differences in free-surface
865                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
866                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
867                  ELSE
868                     h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
869                  ENDIF
870         
871                  IF (N_in * N_out > 0) THEN
872                     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)
873                  ENDIF
874               END DO
875            END DO
876
877            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 
878         ELSE
879
880            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
881
882         ENDIF
883         !
884         DO jk = 1, jpkm1                                 ! Horizontal slab
885            !                                             ! ===============
886
887            !                                             ! --------
888            ! Horizontal divergence                       !   div
889            !                                             ! --------
890            DO jj = j1+1,j2
891               DO ji = i1,i2   ! vector opt.
892                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
893                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kbb_a) * vbdiff(ji,jj  ,jk)  &
894                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk)  ) * zbtr
895               END DO
896            END DO
897            DO jj = j1,j2
898               DO ji = i1,i2-1   ! vector opt.
899                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
900                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
901                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
902               END DO
903            END DO
904         END DO
905
906         !                                                ! ===============
907         !                                               
908
909         imax = i2 - 1
910         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East
911         DO ji = mi0(ind1), mi1(ind1)               
912            imax = MIN(imax,ji)
913         END DO
914         
915         DO jj = j1+1, j2
916            DO ji = i1+1, imax   ! vector opt.
917               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
918                  DO jk = 1, jpkm1
919                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                     &
920                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
921                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
922                  END DO
923               ENDIF
924            END DO
925         END DO
926         !
927         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
928         !
929         DO jj = j1+1, j2-1
930            DO ji = i1+1, i2-1   ! vector opt.
931               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
932                  DO jk = 1, jpkm1
933                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                        &
934                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
935                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                          &
936                        &  - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
937                  END DO
938               ENDIF
939            END DO
940         END DO
941         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
942      ENDIF
943      !
944   END SUBROUTINE interpvn_sponge
945
946   SUBROUTINE interpunb_sponge(tabres,i1,i2,j1,j2, before)
947      !!---------------------------------------------
948      !!   *** ROUTINE interpunb_sponge ***
949      !!---------------------------------------------   
950      INTEGER, INTENT(in) :: i1,i2,j1,j2
951      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
952      LOGICAL, INTENT(in) :: before
953
954      INTEGER  :: ji, jj, ind1, jmax
955      ! sponge parameters
956      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
957      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ubdiff
958      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff
959      !!---------------------------------------------   
960      !
961      IF( before ) THEN
962         DO jj=j1,j2
963            DO ji=i1,i2
964               tabres(ji,jj) = uu_b(ji,jj,Kmm_a)
965            END DO
966         END DO
967
968      ELSE
969
970         ubdiff(i1:i2,j1:j2) = (uu_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*umask(i1:i2,j1:j2,1)
971         !
972         !                                             ! --------
973         ! Horizontal divergence                       !   div
974         !                                             ! --------
975         DO jj = j1,j2
976            DO ji = i1+1,i2   ! vector opt.
977               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj)
978               hdivdiff(ji,jj) = (  e2u(ji  ,jj)*hu(ji  ,jj,Kbb_a) * ubdiff(ji  ,jj) &
979                                  &-e2u(ji-1,jj)*hu(ji-1,jj,Kbb_a) * ubdiff(ji-1,jj) ) * zbtr
980            END DO
981         END DO
982
983         DO jj = j1,j2-1
984            DO ji = i1,i2   ! vector opt.
985               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)
986               rotdiff(ji,jj) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1)   &
987                              &   +e1u(ji,jj  ) * ubdiff(ji,jj  ) ) * fmask(ji,jj,1) * zbtr 
988            END DO
989         END DO
990         !
991         DO jj = j1+1, j2-1
992            DO ji = i1+1, i2-1   ! vector opt.
993               IF (.NOT. tabspongedone_u(ji,jj)) THEN
994                  ze2u = rotdiff (ji,jj)
995                  ze1v = hdivdiff(ji,jj)
996                  ! horizontal diffusive trends
997                  zua = - ( ze2u - rotdiff (ji,jj-1) ) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  &
998                      & + ( hdivdiff(ji+1,jj) - ze1v ) * r1_e1u(ji,jj)                       & 
999                      & - rn_trelax_dyn * r1_Dt * fspu_2d(ji,jj) * ubdiff(ji,jj)
1000
1001                  ! add it to the general momentum trends
1002                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua                                 
1003               ENDIF
1004            END DO
1005         END DO
1006
1007         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
1008
1009         jmax = j2-1
1010         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
1011         DO jj = mj0(ind1), mj1(ind1)                 
1012            jmax = MIN(jmax,jj)
1013         END DO
1014
1015         DO jj = j1+1, jmax
1016            DO ji = i1+1, i2   ! vector opt.
1017               IF (.NOT. tabspongedone_v(ji,jj)) THEN
1018                     ze2u = rotdiff (ji,jj)
1019                     ze1v = hdivdiff(ji,jj)
1020                     zva = + ( ze2u - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) * r1_hv(ji,jj,Kmm_a) &
1021                           + ( hdivdiff(ji,jj+1) - ze1v ) * r1_e2v(ji,jj)
1022                     vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva
1023               ENDIF
1024            END DO
1025         END DO
1026         !
1027         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
1028         !
1029      ENDIF
1030      !
1031   END SUBROUTINE interpunb_sponge
1032
1033   
1034   SUBROUTINE interpvnb_sponge(tabres,i1,i2,j1,j2, before)
1035      !!---------------------------------------------
1036      !!   *** ROUTINE interpvnb_sponge ***
1037      !!---------------------------------------------
1038      INTEGER, INTENT(in) :: i1,i2,j1,j2
1039      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1040      LOGICAL, INTENT(in) :: before
1041      !
1042      INTEGER  ::   ji, jj, ind1, imax
1043      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
1044      REAL(wp), DIMENSION(i1:i2,j1:j2) :: vbdiff
1045      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff
1046      !!---------------------------------------------
1047     
1048      IF( before ) THEN
1049         DO jj=j1,j2
1050            DO ji=i1,i2
1051               tabres(ji,jj) = vv_b(ji,jj,Kmm_a)
1052            END DO
1053         END DO
1054      ELSE
1055         vbdiff(i1:i2,j1:j2) = (vv_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*vmask(i1:i2,j1:j2,1)
1056         !                                             ! --------
1057         ! Horizontal divergence                       !   div
1058         !                                             ! --------
1059         DO jj = j1+1,j2
1060            DO ji = i1,i2   ! vector opt.
1061               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj)
1062               hdivdiff(ji,jj) = ( e1v(ji,jj  ) * hv(ji,jj  ,Kbb_a) * vbdiff(ji,jj  )  &
1063                               &  -e1v(ji,jj-1) * hv(ji,jj-1,Kbb_a) * vbdiff(ji,jj-1)  ) * zbtr
1064            END DO
1065         END DO
1066         DO jj = j1,j2
1067            DO ji = i1,i2-1   ! vector opt.
1068               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj) 
1069               rotdiff(ji,jj) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj) & 
1070                              &  -e2v(ji  ,jj) * vbdiff(ji  ,jj)  ) * fmask(ji,jj,1) * zbtr
1071            END DO
1072         END DO
1073         !                                                ! ===============
1074         !                                               
1075
1076         imax = i2 - 1
1077         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East
1078         DO ji = mi0(ind1), mi1(ind1)               
1079            imax = MIN(imax,ji)
1080         END DO
1081         
1082         DO jj = j1+1, j2
1083            DO ji = i1+1, imax   ! vector opt.
1084               IF( .NOT. tabspongedone_u(ji,jj) ) THEN                                                     
1085                  zua = - ( rotdiff (ji  ,jj) - rotdiff (ji,jj-1)) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  &
1086                      & + ( hdivdiff(ji+1,jj) - hdivdiff(ji,jj  )) * r1_e1u(ji,jj)
1087                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua
1088               ENDIF
1089            END DO
1090         END DO
1091         !
1092         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
1093         !
1094         DO jj = j1+1, j2-1
1095            DO ji = i1+1, i2-1   ! vector opt.
1096               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
1097                  zva  =  ( rotdiff (ji,jj  ) - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) *r1_hv(ji,jj,Kmm_a) &
1098                     &  + ( hdivdiff(ji,jj+1) - hdivdiff(ji  ,jj) ) * r1_e2v(ji,jj)                     &
1099                     &  - rn_trelax_dyn * r1_Dt * fspv_2d(ji,jj) * vbdiff(ji,jj)
1100                  vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva
1101               ENDIF
1102            END DO
1103         END DO
1104         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
1105      ENDIF
1106      !
1107   END SUBROUTINE interpvnb_sponge
1108
1109
1110#else
1111   !!----------------------------------------------------------------------
1112   !!   Empty module                                          no AGRIF zoom
1113   !!----------------------------------------------------------------------
1114CONTAINS
1115   SUBROUTINE agrif_oce_sponge_empty
1116      WRITE(*,*)  'agrif_oce_sponge : You should not have seen this print! error?'
1117   END SUBROUTINE agrif_oce_sponge_empty
1118#endif
1119
1120   !!======================================================================
1121END MODULE agrif_oce_sponge
Note: See TracBrowser for help on using the repository browser.