New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_oce_sponge.F90 in NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90 @ 11574

Last change on this file since 11574 was 11574, checked in by jchanut, 5 years ago

#2222, import changes from dev_r10973_AGRIF-01_jchanut_small_jpi_jpj (i.e. #2199)

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