source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90 @ 13498

Last change on this file since 13498 was 13498, checked in by jchanut, 6 months ago

#2222 1) Define agrif sponge at f-points 2) Correct bathymetry check with vertical refinement just over needed points

  • Property svn:keywords set to Id
File size: 31.3 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   !! * Substitutions
34#  include "do_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/NST 4.0 , NEMO Consortium (2018)
37   !! $Id$
38   !! Software governed by the CeCILL license (see ./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE Agrif_Sponge_Tra
43      !!----------------------------------------------------------------------
44      !!                 *** ROUTINE Agrif_Sponge_Tra ***
45      !!----------------------------------------------------------------------
46      REAL(wp) ::   zcoef   ! local scalar
47     
48      !!----------------------------------------------------------------------
49      !
50#if defined SPONGE
51      !! Assume persistence:
52      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
53
54      CALL Agrif_Sponge
55      Agrif_SpecialValue    = 0._wp
56      Agrif_UseSpecialValue = .TRUE.
57      l_vremap              = ln_vert_remap
58      tabspongedone_tsn     = .FALSE.
59      !
60      CALL Agrif_Bc_Variable( ts_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
61      !
62      Agrif_UseSpecialValue = .FALSE.
63      l_vremap              = .FALSE.
64#endif
65      !
66      CALL iom_put( 'agrif_spu', fspu(:,:))
67      CALL iom_put( 'agrif_spv', fspv(:,:))
68      !
69   END SUBROUTINE Agrif_Sponge_Tra
70
71
72   SUBROUTINE Agrif_Sponge_dyn
73      !!----------------------------------------------------------------------
74      !!                 *** ROUTINE Agrif_Sponge_dyn ***
75      !!----------------------------------------------------------------------
76      REAL(wp) ::   zcoef   ! local scalar
77      !!----------------------------------------------------------------------
78      !
79#if defined SPONGE
80      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
81
82      Agrif_SpecialValue    = 0._wp
83      Agrif_UseSpecialValue = ln_spc_dyn
84      l_vremap              = ln_vert_remap
85      use_sign_north        = .TRUE.
86      sign_north            = -1._wp
87      !
88      tabspongedone_u = .FALSE.
89      tabspongedone_v = .FALSE.         
90      CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
91      !
92      tabspongedone_u = .FALSE.
93      tabspongedone_v = .FALSE.
94      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
95      !
96      Agrif_UseSpecialValue = .FALSE.
97      use_sign_north        = .FALSE.
98      l_vremap              = .FALSE.
99#endif
100      !
101      CALL iom_put( 'agrif_spt', fspt(:,:))
102      CALL iom_put( 'agrif_spf', fspf(:,:))
103      !
104   END SUBROUTINE Agrif_Sponge_dyn
105
106
107   SUBROUTINE Agrif_Sponge
108      !!----------------------------------------------------------------------
109      !!                 *** ROUTINE  Agrif_Sponge ***
110      !!----------------------------------------------------------------------
111      INTEGER  ::   ji, jj, ind1, ind2
112      INTEGER  ::   ispongearea, jspongearea
113      REAL(wp) ::   z1_ispongearea, z1_jspongearea
114      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
115      !!----------------------------------------------------------------------
116      !
117      ! Sponge 1d example with:
118      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
119      !                       
120      !coarse :     U     T     U     T     U     T     U
121      !|            |           |           |           |
122      !fine :     t u t u t u t u t u t u t u t u t u t u t
123      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0
124      !           |   ghost     | <-- sponge area  -- > |
125      !           |   points    |                       |
126      !                         |--> dynamical interface
127
128#if defined SPONGE || defined SPONGE_TOP
129      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
130         ! Define ramp from boundaries towards domain interior at F-points
131         ! Store it in ztabramp
132
133         ispongearea    = nn_sponge_len * Agrif_irhox()
134         z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp )
135         jspongearea    = nn_sponge_len * Agrif_irhoy()
136         z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp )
137         
138         ztabramp(:,:) = 0._wp
139
140         ! Trick to remove sponge in 2DV domains:
141         IF ( nbcellsx <= 3 ) ispongearea = -1
142         IF ( nbcellsy <= 3 ) jspongearea = -1
143
144         IF( lk_west ) THEN                             ! --- West --- !
145            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
146            ind2 = nn_hls + 1 + nbghostcells + ispongearea 
147            DO ji = mi0(ind1), mi1(ind2)   
148               DO jj = 1, jpj               
149                  ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea
150               END DO
151            END DO
152            ! ghost cells:
153            ind1 = 1
154            ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
155            DO ji = mi0(ind1), mi1(ind2)   
156               DO jj = 1, jpj               
157                  ztabramp(ji,jj) = 1._wp
158               END DO
159            END DO
160         ENDIF
161         IF( lk_east ) THEN                             ! --- East --- !
162            ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1
163            ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1
164            DO ji = mi0(ind1), mi1(ind2)
165               DO jj = 1, jpj
166                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) 
167               END DO
168            END DO
169            ! ghost cells:
170            ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1
171            ind2 = jpiglo - 1
172            DO ji = mi0(ind1), mi1(ind2)
173               DO jj = 1, jpj
174                  ztabramp(ji,jj) = 1._wp
175               END DO
176            END DO
177         ENDIF     
178         IF( lk_south ) THEN                            ! --- South --- !
179            ind1 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells
180            ind2 = nn_hls + 1 + nbghostcells + jspongearea 
181            DO jj = mj0(ind1), mj1(ind2) 
182               DO ji = 1, jpi
183                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
184               END DO
185            END DO
186            ! ghost cells:
187            ind1 = 1
188            ind2 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells
189            DO jj = mj0(ind1), mj1(ind2) 
190               DO ji = 1, jpi
191                  ztabramp(ji,jj) = 1._wp
192               END DO
193            END DO
194         ENDIF
195         IF( lk_north ) THEN                            ! --- North --- !
196            ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1
197            ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1
198            DO jj = mj0(ind1), mj1(ind2)
199               DO ji = 1, jpi
200                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) 
201               END DO
202            END DO
203            ! ghost cells:
204            ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
205            ind2 = jpjglo
206            DO jj = mj0(ind1), mj1(ind2)
207               DO ji = 1, jpi
208                  ztabramp(ji,jj) = 1._wp
209               END DO
210            END DO
211         ENDIF
212         !
213      ENDIF
214
215      ! Tracers
216      IF( .NOT. spongedoneT ) THEN
217         fspu(:,:) = 0._wp
218         fspv(:,:) = 0._wp
219         DO_2D( 0, 0, 0, 0 )
220            fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj)
221            fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj)
222         END_2D
223      ENDIF
224
225      ! Dynamics
226      IF( .NOT. spongedoneU ) THEN
227         fspt(:,:) = 0._wp
228         fspf(:,:) = 0._wp
229         DO_2D( 0, 0, 0, 0 )
230            fspt(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) &
231                                  &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
232            fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
233         END_2D
234      ENDIF
235     
236      IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
237         CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp )
238         spongedoneT = .TRUE.
239         spongedoneU = .TRUE.
240      ENDIF
241      IF( .NOT. spongedoneT ) THEN
242         CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp )
243         spongedoneT = .TRUE.
244      ENDIF
245      IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
246         CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1._wp, fspf, 'F', 1._wp )
247         spongedoneU = .TRUE.
248      ENDIF
249      !
250      ! Remove vertical interpolation where not needed:
251      ! (A null value in mbkx arrays does the job)
252      WHERE (fspu(:,:) == 0._wp) mbku_parent(:,:) = 0
253      WHERE (fspv(:,:) == 0._wp) mbkv_parent(:,:) = 0
254      WHERE (fspt(:,:) == 0._wp) mbkt_parent(:,:) = 0
255      !
256#endif
257      !
258   END SUBROUTINE Agrif_Sponge
259
260   
261   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
262      !!----------------------------------------------------------------------
263      !!                 *** ROUTINE interptsn_sponge ***
264      !!----------------------------------------------------------------------
265      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
266      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
267      LOGICAL                                     , INTENT(in   ) ::   before
268      !
269      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
270      INTEGER  ::   iku, ikv
271      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
272      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
273      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
274      ! vertical interpolation:
275      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
276      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
277      REAL(wp), DIMENSION(k1:k2) :: h_in
278      REAL(wp), DIMENSION(1:jpk) :: h_out
279      INTEGER :: N_in, N_out
280      !!----------------------------------------------------------------------
281      !
282      IF( before ) THEN
283         DO jn = 1, jpts
284            DO jk=k1,k2
285               DO jj=j1,j2
286                  DO ji=i1,i2
287                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
288                  END DO
289               END DO
290            END DO
291         END DO
292
293        IF ( l_vremap ) THEN
294
295           ! Interpolate thicknesses
296           ! Warning: these are masked, hence extrapolated prior interpolation.
297           DO jk=k1,k2
298              DO jj=j1,j2
299                 DO ji=i1,i2
300                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)
301                 END DO
302              END DO
303           END DO
304
305           ! Extrapolate thicknesses in partial bottom cells:
306           ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
307           IF (ln_zps) THEN
308              DO jj=j1,j2
309                 DO ji=i1,i2
310                    jk = mbkt(ji,jj)
311                    tabres(ji,jj,jk,jpts+1) = 0._wp
312                 END DO
313              END DO           
314           END IF
315     
316           ! Save ssh at last level:
317           IF (.NOT.ln_linssh) THEN
318              tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
319           ELSE
320              tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
321           END IF     
322        END IF
323
324      ELSE   
325         !
326         IF ( l_vremap ) THEN
327
328            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
329
330            DO jj=j1,j2
331               DO ji=i1,i2
332                  tabres_child(ji,jj,:,:) = 0._wp 
333                  N_in = mbkt_parent(ji,jj)
334                  zhtot = 0._wp
335                  DO jk=1,N_in !k2 = jpk of parent grid
336                     IF (jk==N_in) THEN
337                        h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot
338                     ELSE
339                        h_in(jk) = tabres(ji,jj,jk,n2)
340                     END IF
341                     zhtot = zhtot + h_in(jk)
342                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
343                  END DO
344                  N_out = 0
345                  DO jk=1,jpk ! jpk of child grid
346                     IF (tmask(ji,jj,jk) == 0) EXIT
347                     N_out = N_out + 1
348                     h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
349                  END DO
350
351                  ! Account for small differences in the free-surface
352                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
353                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
354                  ELSE
355                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
356                  END IF
357                  IF (N_in*N_out > 0) THEN
358                     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)
359                  ENDIF
360               END DO
361            END DO
362
363            DO jj=j1,j2
364               DO ji=i1,i2
365                  DO jk=1,jpkm1
366                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
367                  END DO
368               END DO
369            END DO
370
371         ELSE
372
373            DO jj=j1,j2
374               DO ji=i1,i2
375                  DO jk=1,jpkm1
376                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
377                  END DO
378               END DO
379            END DO
380         END IF
381
382         DO jn = 1, jpts           
383            DO jk = 1, jpkm1
384               ztu(i1:i2,j1:j2,jk) = 0._wp
385               DO jj = j1,j2
386                  DO ji = i1,i2-1
387                     zabe1 = fspu(ji,jj) * rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
388                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
389                  END DO
390               END DO
391               ztv(i1:i2,j1:j2,jk) = 0._wp
392               DO ji = i1,i2
393                  DO jj = j1,j2-1
394                     zabe2 = fspv(ji,jj) * rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
395                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
396                  END DO
397               END DO
398               !
399               IF( ln_zps ) THEN      ! set gradient at partial step level
400                  DO jj = j1,j2
401                     DO ji = i1,i2
402                        ! last level
403                        iku = mbku(ji,jj)
404                        ikv = mbkv(ji,jj)
405                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
406                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
407                     END DO
408                  END DO
409               ENDIF
410            END DO
411            !
412            DO jk = 1, jpkm1
413               DO jj = j1+1,j2-1
414                  DO ji = i1+1,i2-1
415                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
416                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
417                        ! horizontal diffusive trends
418                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) &
419                             &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
420                        ! add it to the general tracer trends
421                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
422                     ENDIF
423                  END DO
424               END DO
425            END DO
426            !
427         END DO
428         !
429         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
430         !
431      ENDIF
432      !
433   END SUBROUTINE interptsn_sponge
434
435   
436   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
437      !!---------------------------------------------
438      !!   *** ROUTINE interpun_sponge ***
439      !!---------------------------------------------   
440      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
441      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
442      LOGICAL, INTENT(in) :: before
443
444      INTEGER  :: ji,jj,jk,jmax
445      INTEGER  :: ind1
446      ! sponge parameters
447      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
448      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
449      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
450      ! vertical interpolation:
451      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
452      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
453      REAL(wp), DIMENSION(1:jpk) :: h_out
454      INTEGER ::N_in, N_out
455      !!---------------------------------------------   
456      !
457      IF( before ) THEN
458         DO jk=k1,k2
459            DO jj=j1,j2
460               DO ji=i1,i2
461                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
462               END DO
463            END DO
464         END DO
465
466         IF ( l_vremap ) THEN
467
468            DO jk=k1,k2
469               DO jj=j1,j2
470                  DO ji=i1,i2
471                     tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
472                  END DO
473               END DO
474            END DO
475
476            ! Extrapolate thicknesses in partial bottom cells:
477            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
478            IF (ln_zps) THEN
479               DO jj=j1,j2
480                  DO ji=i1,i2
481                     jk = mbku(ji,jj)
482                     tabres(ji,jj,jk,m2) = 0._wp
483                  END DO
484               END DO           
485            END IF
486            ! Save ssh at last level:
487            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
488            IF (.NOT.ln_linssh) THEN
489               ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
490               DO jk=1,jpk
491                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
492               END DO
493               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
494            END IF
495         END IF
496
497      ELSE
498
499         IF ( l_vremap ) THEN
500
501            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
502
503            DO jj=j1,j2
504               DO ji=i1,i2
505                  tabres_child(ji,jj,:) = 0._wp
506                  N_in = mbku_parent(ji,jj)
507                  zhtot = 0._wp
508                  DO jk=1,N_in
509                     IF (jk==N_in) THEN
510                        h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
511                     ELSE
512                        h_in(jk) = tabres(ji,jj,jk,m2)
513                     ENDIF
514                     zhtot = zhtot + h_in(jk)
515                     tabin(jk) = tabres(ji,jj,jk,m1)
516                  END DO
517                  !         
518                  N_out = 0
519                  DO jk=1,jpk
520                     IF (umask(ji,jj,jk) == 0) EXIT
521                     N_out = N_out + 1
522                     h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
523                  END DO
524
525                  ! Account for small differences in free-surface
526                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
527                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
528                  ELSE
529                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
530                  ENDIF
531                 
532                  IF (N_in * N_out > 0) THEN
533                     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)
534                  ENDIF
535               END DO
536            END DO
537
538            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
539
540         ELSE
541
542            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
543
544         ENDIF
545         !
546         DO jk = 1, jpkm1                                 ! Horizontal slab
547            !                                             ! ===============
548
549            !                                             ! --------
550            ! Horizontal divergence                       !   div
551            !                                             ! --------
552            DO jj = j1,j2
553               DO ji = i1+1,i2   ! vector opt.
554                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
555                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kbb_a) * ubdiff(ji  ,jj,jk) &
556                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr
557               END DO
558            END DO
559
560            DO jj = j1,j2-1
561               DO ji = i1,i2   ! vector opt.
562                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
563                  rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
564                                    &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
565               END DO
566            END DO
567         END DO
568         !
569         DO jj = j1+1, j2-1
570            DO ji = i1+1, i2-1   ! vector opt.
571
572               IF (.NOT. tabspongedone_u(ji,jj)) THEN
573                  DO jk = 1, jpkm1                                 ! Horizontal slab
574                     ze2u = rotdiff (ji,jj,jk)
575                     ze1v = hdivdiff(ji,jj,jk)
576                     ! horizontal diffusive trends
577                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )   &
578                         & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 
579                         & - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk)
580
581                     ! add it to the general momentum trends
582                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua                                 
583                  END DO
584               ENDIF
585
586            END DO
587         END DO
588
589         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
590
591         jmax = j2-1
592         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
593         DO jj = mj0(ind1), mj1(ind1)                 
594            jmax = MIN(jmax,jj)
595         END DO
596
597         DO jj = j1+1, jmax
598            DO ji = i1+1, i2   ! vector opt.
599
600               IF (.NOT. tabspongedone_v(ji,jj)) THEN
601                  DO jk = 1, jpkm1                                 ! Horizontal slab
602                     ze2u = rotdiff (ji,jj,jk)
603                     ze1v = hdivdiff(ji,jj,jk)
604
605                     ! horizontal diffusive trends
606                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
607                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
608
609                     ! add it to the general momentum trends
610                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
611                  END DO
612               ENDIF
613               !
614            END DO
615         END DO
616         !
617         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
618         !
619      ENDIF
620      !
621   END SUBROUTINE interpun_sponge
622
623   
624   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
625      !!---------------------------------------------
626      !!   *** ROUTINE interpvn_sponge ***
627      !!---------------------------------------------
628      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
629      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
630      LOGICAL, INTENT(in) :: before
631      !
632      INTEGER  ::   ji, jj, jk, imax
633      INTEGER  :: ind1
634      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot
635      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
636      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
637      ! vertical interpolation:
638      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
639      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
640      REAL(wp), DIMENSION(1:jpk) :: h_out
641      INTEGER :: N_in, N_out
642      !!---------------------------------------------
643     
644      IF( before ) THEN
645         DO jk=k1,k2
646            DO jj=j1,j2
647               DO ji=i1,i2
648                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
649               END DO
650            END DO
651         END DO
652
653         IF ( l_vremap ) THEN
654
655            DO jk=k1,k2
656               DO jj=j1,j2
657                  DO ji=i1,i2
658                     tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
659                  END DO
660               END DO
661            END DO
662            ! Extrapolate thicknesses in partial bottom cells:
663            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
664            IF (ln_zps) THEN
665               DO jj=j1,j2
666                  DO ji=i1,i2
667                     jk = mbkv(ji,jj)
668                     tabres(ji,jj,jk,m2) = 0._wp
669                  END DO
670               END DO           
671            END IF
672            ! Save ssh at last level:
673            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
674            IF (.NOT.ln_linssh) THEN
675               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
676               DO jk=1,jpk
677                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
678               END DO
679               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
680            END IF
681
682         END IF
683
684      ELSE
685
686         IF ( l_vremap ) THEN
687            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
688            DO jj=j1,j2
689               DO ji=i1,i2
690                  tabres_child(ji,jj,:) = 0._wp
691                  N_in = mbkv_parent(ji,jj)
692                  zhtot = 0._wp
693                  DO jk=1,N_in
694                     IF (jk==N_in) THEN
695                        h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
696                     ELSE
697                        h_in(jk) = tabres(ji,jj,jk,m2)
698                     ENDIF
699                     zhtot = zhtot + h_in(jk)
700                     tabin(jk) = tabres(ji,jj,jk,m1)
701                  END DO
702                  !         
703                  N_out = 0
704                  DO jk=1,jpk
705                     IF (vmask(ji,jj,jk) == 0) EXIT
706                     N_out = N_out + 1
707                     h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
708                  END DO
709
710                  ! Account for small differences in free-surface
711                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
712                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
713                  ELSE
714                     h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
715                  ENDIF
716         
717                  IF (N_in * N_out > 0) THEN
718                     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)
719                  ENDIF
720               END DO
721            END DO
722
723            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 
724
725         ELSE
726
727            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
728
729         ENDIF
730         !
731         DO jk = 1, jpkm1                                 ! Horizontal slab
732            !                                             ! ===============
733
734            !                                             ! --------
735            ! Horizontal divergence                       !   div
736            !                                             ! --------
737            DO jj = j1+1,j2
738               DO ji = i1,i2   ! vector opt.
739                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
740                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kbb_a) * vbdiff(ji,jj  ,jk)  &
741                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk)  ) * zbtr
742               END DO
743            END DO
744            DO jj = j1,j2
745               DO ji = i1,i2-1   ! vector opt.
746                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
747                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
748                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
749               END DO
750            END DO
751         END DO
752
753         !                                                ! ===============
754         !                                               
755
756         imax = i2 - 1
757         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East
758         DO ji = mi0(ind1), mi1(ind1)               
759            imax = MIN(imax,ji)
760         END DO
761         
762         DO jj = j1+1, j2
763            DO ji = i1+1, imax   ! vector opt.
764               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
765                  DO jk = 1, jpkm1
766                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                     &
767                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
768                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
769                  END DO
770               ENDIF
771            END DO
772         END DO
773         !
774         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
775         !
776         DO jj = j1+1, j2-1
777            DO ji = i1+1, i2-1   ! vector opt.
778               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
779                  DO jk = 1, jpkm1
780                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                        &
781                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
782                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                          &
783                        &  - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
784                  END DO
785               ENDIF
786            END DO
787         END DO
788         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
789      ENDIF
790      !
791   END SUBROUTINE interpvn_sponge
792
793#else
794   !!----------------------------------------------------------------------
795   !!   Empty module                                          no AGRIF zoom
796   !!----------------------------------------------------------------------
797CONTAINS
798   SUBROUTINE agrif_oce_sponge_empty
799      WRITE(*,*)  'agrif_oce_sponge : You should not have seen this print! error?'
800   END SUBROUTINE agrif_oce_sponge_empty
801#endif
802
803   !!======================================================================
804END MODULE agrif_oce_sponge
Note: See TracBrowser for help on using the repository browser.