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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90 @ 10989

Last change on this file since 10989 was 10989, checked in by acc, 5 years ago

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