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/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST – NEMO

source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce_sponge.F90 @ 11542

Last change on this file since 11542 was 11542, checked in by jchanut, 5 years ago

#2199: 1) Replace TWO_WAY cpp key by namelist flag ln_agrif_2way. 2) Clean agrif_user module. At this stage, code is reproducible in one way mode down to jpi/jpj=3. Issues for small subdomain sizes (e.g. aprox < 7x7) remain in update (i.e. in 2 way mode). Masking of sponge in unconnected areas has to be implemented.
-Cette ligne, et les suivantes ci-dessous, seront ignorées--

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