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/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

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

Last change on this file since 13565 was 13565, checked in by jchanut, 4 years ago

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