source: NEMO/releases/r4.0/r4.0-HEAD/src/NST/agrif_oce_sponge.F90 @ 12737

Last change on this file since 12737 was 12737, checked in by jchanut, 8 months ago

Fixes AGRIF reproductibility with land processors removal, i.e. #2240. Trunk is not concerned by this problem since nbondi/nbondj variables are not used anymore.

  • Property svn:keywords set to Id
File size: 23.2 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
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
29   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
30
31   !!----------------------------------------------------------------------
32   !! NEMO/NST 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE Agrif_Sponge_Tra
39      !!----------------------------------------------------------------------
40      !!                 *** ROUTINE Agrif_Sponge_Tra ***
41      !!----------------------------------------------------------------------
42      REAL(wp) ::   zcoef   ! local scalar
43     
44      !!----------------------------------------------------------------------
45      !
46#if defined SPONGE
47      !! Assume persistence:
48      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
49
50      CALL Agrif_Sponge
51      Agrif_SpecialValue    = 0._wp
52      Agrif_UseSpecialValue = .TRUE.
53      tabspongedone_tsn     = .FALSE.
54      !
55      CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
56      !
57      Agrif_UseSpecialValue = .FALSE.
58#endif
59      !
60   END SUBROUTINE Agrif_Sponge_Tra
61
62
63   SUBROUTINE Agrif_Sponge_dyn
64      !!----------------------------------------------------------------------
65      !!                 *** ROUTINE Agrif_Sponge_dyn ***
66      !!----------------------------------------------------------------------
67      REAL(wp) ::   zcoef   ! local scalar
68      !!----------------------------------------------------------------------
69      !
70#if defined SPONGE
71      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
72
73      Agrif_SpecialValue=0.
74      Agrif_UseSpecialValue = ln_spc_dyn
75      !
76      tabspongedone_u = .FALSE.
77      tabspongedone_v = .FALSE.         
78      CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
79      !
80      tabspongedone_u = .FALSE.
81      tabspongedone_v = .FALSE.
82      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
83      !
84      Agrif_UseSpecialValue = .FALSE.
85#endif
86      !
87   END SUBROUTINE Agrif_Sponge_dyn
88
89
90   SUBROUTINE Agrif_Sponge
91      !!----------------------------------------------------------------------
92      !!                 *** ROUTINE  Agrif_Sponge ***
93      !!----------------------------------------------------------------------
94      INTEGER  ::   ji, jj, ind1, ind2
95      INTEGER  ::   ispongearea
96      REAL(wp) ::   z1_spongearea
97      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
98      !!----------------------------------------------------------------------
99      !
100#if defined SPONGE || defined SPONGE_TOP
101      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
102         ! Define ramp from boundaries towards domain interior at T-points
103         ! Store it in ztabramp
104
105         ispongearea  = 1 + nn_sponge_len * Agrif_irhox()
106         z1_spongearea = 1._wp / REAL( ispongearea )
107         
108         ztabramp(:,:) = 0._wp
109
110         ! --- West --- !
111         IF( l_Westedge ) THEN
112            ind1 = 1+nbghostcells
113            ind2 = 1+nbghostcells + ispongearea 
114            DO jj = 1, jpj
115               DO ji = ind1, ind2               
116                  ztabramp(ji,jj) = REAL( ind2 - ji ) * z1_spongearea * umask(ind1,jj,1)
117               END DO
118            ENDDO
119         ENDIF
120
121         ! --- East --- !
122         IF( l_Eastedge ) THEN
123            ind1 = nlci - nbghostcells - ispongearea
124            ind2 = nlci - nbghostcells
125            DO jj = 1, jpj
126               DO ji = ind1, ind2
127                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ji - ind1 ) * z1_spongearea * umask(ind2-1,jj,1) )
128               ENDDO
129            ENDDO
130         ENDIF
131
132         ! --- South --- !
133         IF( l_Southedge ) THEN
134            ind1 = 1+nbghostcells
135            ind2 = 1+nbghostcells + ispongearea
136            DO jj = ind1, ind2 
137               DO ji = 1, jpi
138                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - jj ) * z1_spongearea * vmask(ji,ind1,1) )
139               END DO
140            ENDDO
141         ENDIF
142
143         ! --- North --- !
144         IF( l_Northedge ) THEN
145            ind1 = nlcj - nbghostcells - ispongearea
146            ind2 = nlcj - nbghostcells
147            DO jj = ind1, ind2
148               DO ji = 1, jpi
149                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( jj - ind1 ) * z1_spongearea * vmask(ji,ind2-1,1) )
150               END DO
151            ENDDO
152         ENDIF
153
154      ENDIF
155
156      ! Tracers
157      IF( .NOT. spongedoneT ) THEN
158         fsaht_spu(:,:) = 0._wp
159         fsaht_spv(:,:) = 0._wp
160         DO jj = 2, jpjm1
161            DO ji = 2, jpim1   ! vector opt.
162               fsaht_spu(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) )
163               fsaht_spv(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) )
164            END DO
165         END DO
166         CALL lbc_lnk( 'agrif_oce_sponge', fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
167         CALL lbc_lnk( 'agrif_oce_sponge', fsaht_spv, 'V', 1. )
168         
169         spongedoneT = .TRUE.
170      ENDIF
171
172      ! Dynamics
173      IF( .NOT. spongedoneU ) THEN
174         fsahm_spt(:,:) = 0._wp
175         fsahm_spf(:,:) = 0._wp
176         DO jj = 2, jpjm1
177            DO ji = 2, jpim1   ! vector opt.
178               fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
179               fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1) &
180                                                     &  +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj  ) )
181            END DO
182         END DO
183         CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
184         CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spf, 'F', 1. )
185         
186         spongedoneU = .TRUE.
187      ENDIF
188      !
189#endif
190      !
191   END SUBROUTINE Agrif_Sponge
192
193   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
194      !!----------------------------------------------------------------------
195      !!                 *** ROUTINE interptsn_sponge ***
196      !!----------------------------------------------------------------------
197      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
198      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
199      LOGICAL                                     , INTENT(in   ) ::   before
200      !
201      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
202      INTEGER  ::   iku, ikv
203      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
204      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
205      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
206      ! vertical interpolation:
207      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
208      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
209      REAL(wp), DIMENSION(k1:k2) :: h_in
210      REAL(wp), DIMENSION(1:jpk) :: h_out
211      INTEGER :: N_in, N_out
212      REAL(wp) :: h_diff
213      !!----------------------------------------------------------------------
214      !
215      IF( before ) THEN
216         DO jn = 1, jpts
217            DO jk=k1,k2
218               DO jj=j1,j2
219                  DO ji=i1,i2
220                     tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)
221                  END DO
222               END DO
223            END DO
224         END DO
225
226# if defined key_vertical
227         DO jk=k1,k2
228            DO jj=j1,j2
229               DO ji=i1,i2
230                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
231               END DO
232            END DO
233         END DO
234# endif
235
236      ELSE   
237         !
238# if defined key_vertical
239         tabres_child(:,:,:,:) = 0.
240         DO jj=j1,j2
241            DO ji=i1,i2
242               N_in = 0
243               DO jk=k1,k2 !k2 = jpk of parent grid
244                  IF (tabres(ji,jj,jk,n2) == 0) EXIT
245                  N_in = N_in + 1
246                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
247                  h_in(N_in) = tabres(ji,jj,jk,n2)
248               END DO
249               N_out = 0
250               DO jk=1,jpk ! jpk of child grid
251                  IF (tmask(ji,jj,jk) == 0) EXIT
252                  N_out = N_out + 1
253                  h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
254               ENDDO
255               IF (N_in > 0) THEN
256                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
257                  tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
258                  DO jn=1,jpts
259                     call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
260                  ENDDO
261               ENDIF
262            ENDDO
263         ENDDO
264# endif
265
266         DO jj=j1,j2
267            DO ji=i1,i2
268               DO jk=1,jpkm1
269# if defined key_vertical
270                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)
271# else
272                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
273# endif
274               ENDDO
275            ENDDO
276         ENDDO
277
278         DO jn = 1, jpts           
279            DO jk = 1, jpkm1
280               ztu(i1:i2,j1:j2,jk) = 0._wp
281               DO jj = j1,j2
282                  DO ji = i1,i2-1
283                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
284                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
285                  END DO
286               END DO
287               ztv(i1:i2,j1:j2,jk) = 0._wp
288               DO ji = i1,i2
289                  DO jj = j1,j2-1
290                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
291                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
292                  END DO
293               END DO
294               !
295               IF( ln_zps ) THEN      ! set gradient at partial step level
296                  DO jj = j1,j2
297                     DO ji = i1,i2
298                        ! last level
299                        iku = mbku(ji,jj)
300                        ikv = mbkv(ji,jj)
301                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
302                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
303                     END DO
304                  END DO
305               ENDIF
306            END DO
307            !
308            DO jk = 1, jpkm1
309               DO jj = j1+1,j2-1
310                  DO ji = i1+1,i2-1
311                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
312                        zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
313                        ! horizontal diffusive trends
314                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
315                        ! add it to the general tracer trends
316                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
317                     ENDIF
318                  END DO
319               END DO
320            END DO
321            !
322         END DO
323         !
324         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
325         !
326      ENDIF
327      !
328   END SUBROUTINE interptsn_sponge
329
330   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
331      !!---------------------------------------------
332      !!   *** ROUTINE interpun_sponge ***
333      !!---------------------------------------------   
334      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
335      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
336      LOGICAL, INTENT(in) :: before
337
338      INTEGER :: ji,jj,jk,jmax
339
340      ! sponge parameters
341      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
342      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
343      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
344      ! vertical interpolation:
345      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
346      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
347      REAL(wp), DIMENSION(1:jpk) :: h_out
348      INTEGER ::N_in,N_out
349      !!---------------------------------------------   
350      !
351      IF( before ) THEN
352         DO jk=1,jpkm1
353            DO jj=j1,j2
354               DO ji=i1,i2
355                  tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
356# if defined key_vertical
357                  tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
358# endif
359               END DO
360            END DO
361         END DO
362
363      ELSE
364
365# if defined key_vertical
366         tabres_child(:,:,:) = 0._wp
367         DO jj=j1,j2
368            DO ji=i1,i2
369               N_in = 0
370               DO jk=k1,k2
371                  IF (tabres(ji,jj,jk,m2) == 0) EXIT
372                  N_in = N_in + 1
373                  tabin(jk) = tabres(ji,jj,jk,m1)
374                  h_in(N_in) = tabres(ji,jj,jk,m2)
375              ENDDO
376              !
377              IF (N_in == 0) THEN
378                 tabres_child(ji,jj,:) = 0.
379                 CYCLE
380              ENDIF
381         
382              N_out = 0
383              DO jk=1,jpk
384                 if (umask(ji,jj,jk) == 0) EXIT
385                 N_out = N_out + 1
386                 h_out(N_out) = e3u_n(ji,jj,jk)
387              ENDDO
388         
389              IF (N_out == 0) THEN
390                 tabres_child(ji,jj,:) = 0.
391                 CYCLE
392              ENDIF
393         
394              IF (N_in * N_out > 0) THEN
395                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
396                 if (h_diff < -1.e4) then
397                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
398                 endif
399              ENDIF
400              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)
401         
402            ENDDO
403         ENDDO
404
405         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
406#else
407         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
408#endif
409         !
410         DO jk = 1, jpkm1                                 ! Horizontal slab
411            !                                             ! ===============
412
413            !                                             ! --------
414            ! Horizontal divergence                       !   div
415            !                                             ! --------
416            DO jj = j1,j2
417               DO ji = i1+1,i2   ! vector opt.
418                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
419                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
420                                     &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
421               END DO
422            END DO
423
424            DO jj = j1,j2-1
425               DO ji = i1,i2   ! vector opt.
426                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
427                  rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
428                                    &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
429               END DO
430            END DO
431         END DO
432         !
433         DO jj = j1+1, j2-1
434            DO ji = i1+1, i2-1   ! vector opt.
435
436               IF (.NOT. tabspongedone_u(ji,jj)) THEN
437                  DO jk = 1, jpkm1                                 ! Horizontal slab
438                     ze2u = rotdiff (ji,jj,jk)
439                     ze1v = hdivdiff(ji,jj,jk)
440                     ! horizontal diffusive trends
441                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
442                           + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
443
444                     ! add it to the general momentum trends
445                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
446
447                  END DO
448               ENDIF
449
450            END DO
451         END DO
452
453         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
454
455         jmax = j2-1
456         IF ( l_Northedge ) jmax = MIN(jmax,nlcj-nbghostcells-2)   ! North
457
458         DO jj = j1+1, jmax
459            DO ji = i1+1, i2   ! vector opt.
460
461               IF (.NOT. tabspongedone_v(ji,jj)) THEN
462                  DO jk = 1, jpkm1                                 ! Horizontal slab
463                     ze2u = rotdiff (ji,jj,jk)
464                     ze1v = hdivdiff(ji,jj,jk)
465
466                     ! horizontal diffusive trends
467                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
468                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
469
470                     ! add it to the general momentum trends
471                     va(ji,jj,jk) = va(ji,jj,jk) + zva
472                  END DO
473               ENDIF
474               !
475            END DO
476         END DO
477         !
478         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
479         !
480      ENDIF
481      !
482   END SUBROUTINE interpun_sponge
483
484   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
485      !!---------------------------------------------
486      !!   *** ROUTINE interpvn_sponge ***
487      !!---------------------------------------------
488      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
489      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
490      LOGICAL, INTENT(in) :: before
491      INTEGER, INTENT(in) :: nb , ndir
492      !
493      INTEGER  ::   ji, jj, jk, imax
494      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff
495      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
496      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
497      ! vertical interpolation:
498      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
499      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
500      REAL(wp), DIMENSION(1:jpk) :: h_out
501      INTEGER :: N_in, N_out
502      !!---------------------------------------------
503     
504      IF( before ) THEN
505         DO jk=1,jpkm1
506            DO jj=j1,j2
507               DO ji=i1,i2
508                  tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
509# if defined key_vertical
510                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
511# endif
512               END DO
513            END DO
514         END DO
515      ELSE
516
517# if defined key_vertical
518         tabres_child(:,:,:) = 0._wp
519         DO jj=j1,j2
520            DO ji=i1,i2
521               N_in = 0
522               DO jk=k1,k2
523                  IF (tabres(ji,jj,jk,m2) == 0) EXIT
524                  N_in = N_in + 1
525                  tabin(jk) = tabres(ji,jj,jk,m1)
526                  h_in(N_in) = tabres(ji,jj,jk,m2)
527              ENDDO
528         
529              IF (N_in == 0) THEN
530                 tabres_child(ji,jj,:) = 0.
531                 CYCLE
532              ENDIF
533         
534              N_out = 0
535              DO jk=1,jpk
536                 if (vmask(ji,jj,jk) == 0) EXIT
537                 N_out = N_out + 1
538                 h_out(N_out) = e3v_n(ji,jj,jk)
539              ENDDO
540         
541              IF (N_in * N_out > 0) THEN
542                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
543                 if (h_diff < -1.e4) then
544                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
545                 endif
546              ENDIF
547              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)
548            ENDDO
549         ENDDO
550
551         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 
552# else
553         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
554# endif
555         !
556         DO jk = 1, jpkm1                                 ! Horizontal slab
557            !                                             ! ===============
558
559            !                                             ! --------
560            ! Horizontal divergence                       !   div
561            !                                             ! --------
562            DO jj = j1+1,j2
563               DO ji = i1,i2   ! vector opt.
564                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
565                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
566                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
567               END DO
568            END DO
569            DO jj = j1,j2
570               DO ji = i1,i2-1   ! vector opt.
571                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
572                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
573                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
574               END DO
575            END DO
576         END DO
577
578         !                                                ! ===============
579         !                                               
580
581         imax = i2 - 1
582         IF ( l_Eastedge )   imax = MIN(imax,nlci-nbghostcells-2)   ! East
583
584         DO jj = j1+1, j2
585            DO ji = i1+1, imax   ! vector opt.
586               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
587                  DO jk = 1, jpkm1
588                     ua(ji,jj,jk) = ua(ji,jj,jk)                                                               &
589                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  &
590                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
591                  END DO
592               ENDIF
593            END DO
594         END DO
595         !
596         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
597         !
598         DO jj = j1+1, j2-1
599            DO ji = i1+1, i2-1   ! vector opt.
600               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
601                  DO jk = 1, jpkm1
602                     va(ji,jj,jk) = va(ji,jj,jk)                                                                  &
603                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
604                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
605                  END DO
606               ENDIF
607            END DO
608         END DO
609         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
610      ENDIF
611      !
612   END SUBROUTINE interpvn_sponge
613
614#else
615   !!----------------------------------------------------------------------
616   !!   Empty module                                          no AGRIF zoom
617   !!----------------------------------------------------------------------
618CONTAINS
619   SUBROUTINE agrif_oce_sponge_empty
620      WRITE(*,*)  'agrif_oce_sponge : You should not have seen this print! error?'
621   END SUBROUTINE agrif_oce_sponge_empty
622#endif
623
624   !!======================================================================
625END MODULE agrif_oce_sponge
Note: See TracBrowser for help on using the repository browser.