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/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_oce_sponge.F90 @ 9748

Last change on this file since 9748 was 9598, checked in by nicolasmartin, 6 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

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