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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90 @ 11603

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

#2222, 1) correct time interpolation of barotropic velocities in corners. 2) Clean remapping module and enable remapping several variables at the same time. At this stage, vertical remapping doesn't change VORTEX results with an identical vertical grid ONLY in one way mode and a linearized free surface (within truncature errors).

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