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

Last change on this file since 9046 was 9046, checked in by timgraham, 6 years ago

Remove missed ">>>>right" line from merge
  • 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         !
403         DO jk = 1, jpkm1                                 ! Horizontal slab
404            !                                             ! ===============
405
406            !                                             ! --------
407            ! Horizontal divergence                       !   div
408            !                                             ! --------
409            DO jj = j1,j2
410               DO ji = i1+1,i2   ! vector opt.
411                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
412                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
413                                     &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
414               END DO
415            END DO
416
417            DO jj = j1,j2-1
418               DO ji = i1,i2   ! vector opt.
419                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
420                  rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
421                                    &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
422               END DO
423            END DO
424         END DO
425         !
426         DO jj = j1+1, j2-1
427            DO ji = i1+1, i2-1   ! vector opt.
428
429               IF (.NOT. tabspongedone_u(ji,jj)) THEN
430                  DO jk = 1, jpkm1                                 ! Horizontal slab
431                     ze2u = rotdiff (ji,jj,jk)
432                     ze1v = hdivdiff(ji,jj,jk)
433                     ! horizontal diffusive trends
434                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
435                           + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
436
437                     ! add it to the general momentum trends
438                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
439
440                  END DO
441               ENDIF
442
443            END DO
444         END DO
445
446         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
447
448         jmax = j2-1
449         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2)   ! North
450
451         DO jj = j1+1, jmax
452            DO ji = i1+1, i2   ! vector opt.
453
454               IF (.NOT. tabspongedone_v(ji,jj)) THEN
455                  DO jk = 1, jpkm1                                 ! Horizontal slab
456                     ze2u = rotdiff (ji,jj,jk)
457                     ze1v = hdivdiff(ji,jj,jk)
458
459                     ! horizontal diffusive trends
460                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
461                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
462
463                     ! add it to the general momentum trends
464                     va(ji,jj,jk) = va(ji,jj,jk) + zva
465                  END DO
466               ENDIF
467               !
468            END DO
469         END DO
470         !
471         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
472         !
473      ENDIF
474      !
475   END SUBROUTINE interpun_sponge
476
477   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
478      !!---------------------------------------------
479      !!   *** ROUTINE interpvn_sponge ***
480      !!---------------------------------------------
481      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
482      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
483      LOGICAL, INTENT(in) :: before
484      INTEGER, INTENT(in) :: nb , ndir
485      !
486      INTEGER ::   ji, jj, jk
487      INTEGER ::   imax
488      REAL(wp)::   ze2u, ze1v, zua, zva, zbtr
489      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) ::   vbdiff, rotdiff, hdivdiff
490      !!----------------------------------------------------------------------
491      INTEGER  ::   ji, jj, jk, imax
492      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff
493      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
494      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
495      ! vertical interpolation:
496      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
497      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
498      REAL(wp), DIMENSION(1:jpk) :: h_out
499      INTEGER :: N_in, N_out
500      !!---------------------------------------------
501>>>>>>> .merge-right.r9019
502     
503      IF( before ) THEN
504         DO jk=1,jpkm1
505            DO jj=j1,j2
506               DO ji=i1,i2
507                  tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
508# if defined key_vertical
509                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
510# endif
511               END DO
512            END DO
513         END DO
514      ELSE
515
516# if defined key_vertical
517         tabres_child(:,:,:) = 0._wp
518         DO jj=j1,j2
519            DO ji=i1,i2
520               N_in = 0
521               DO jk=k1,k2
522                  IF (tabres(ji,jj,jk,m2) == 0) EXIT
523                  N_in = N_in + 1
524                  tabin(jk) = tabres(ji,jj,jk,m1)
525                  h_in(N_in) = tabres(ji,jj,jk,m2)
526              ENDDO
527         
528              IF (N_in == 0) THEN
529                 tabres_child(ji,jj,:) = 0.
530                 CYCLE
531              ENDIF
532         
533              N_out = 0
534              DO jk=1,jpk
535                 if (vmask(ji,jj,jk) == 0) EXIT
536                 N_out = N_out + 1
537                 h_out(N_out) = e3v_n(ji,jj,jk)
538              ENDDO
539         
540              IF (N_in * N_out > 0) THEN
541                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
542                 if (h_diff < -1.e4) then
543                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
544                 endif
545              ENDIF
546              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)
547            ENDDO
548         ENDDO
549
550         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 
551# else
552         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
553# endif
554         !
555         DO jk = 1, jpkm1                                 ! Horizontal slab
556            !                                             ! ===============
557
558            !                                             ! --------
559            ! Horizontal divergence                       !   div
560            !                                             ! --------
561            DO jj = j1+1,j2
562               DO ji = i1,i2   ! vector opt.
563                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
564                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
565                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
566               END DO
567            END DO
568            DO jj = j1,j2
569               DO ji = i1,i2-1   ! vector opt.
570                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
571                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
572                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
573               END DO
574            END DO
575         END DO
576
577         !                                                ! ===============
578         !                                               
579
580         imax = i2 - 1
581         IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-nbghostcells-2)   ! East
582
583         DO jj = j1+1, j2
584            DO ji = i1+1, imax   ! vector opt.
585               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
586                  DO jk = 1, jpkm1
587                     ua(ji,jj,jk) = ua(ji,jj,jk)                                                               &
588                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  &
589                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
590                  END DO
591               ENDIF
592            END DO
593         END DO
594         !
595         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
596         !
597         DO jj = j1+1, j2-1
598            DO ji = i1+1, i2-1   ! vector opt.
599               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
600                  DO jk = 1, jpkm1
601                     va(ji,jj,jk) = va(ji,jj,jk)                                                                  &
602                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
603                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
604                  END DO
605               ENDIF
606            END DO
607         END DO
608         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
609      ENDIF
610      !
611   END SUBROUTINE interpvn_sponge
612
613#else
614   !!----------------------------------------------------------------------
615   !!   Empty module                                          no AGRIF zoom
616   !!----------------------------------------------------------------------
617CONTAINS
618   SUBROUTINE agrif_opa_sponge_empty
619      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
620   END SUBROUTINE agrif_opa_sponge_empty
621#endif
622
623   !!======================================================================
624END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.