source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/NST/agrif_oce_sponge.F90 @ 12460

Last change on this file since 12460 was 12377, checked in by acc, 16 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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