New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_opa_sponge.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 9031

Last change on this file since 9031 was 9031, checked in by timgraham, 7 years ago

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