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_opa_interp.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 9116

Last change on this file since 9116 was 9116, checked in by jchanut, 6 years ago

Finalize AGRIF ghost cells implementation: ensure compatibility with bdy smoothing or extrapolation
  • Property svn:keywords set to Id
File size: 56.5 KB
Line 
1MODULE agrif_opa_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_opa_interp  ***
4   !! AGRIF: interpolation package for the ocean dynamics (OPA)
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (L. Debreu)  Original cade
7   !!            3.2  !  2009-04  (R. Benshila)
8   !!            3.6  !  2014-09  (R. Benshila)
9   !!----------------------------------------------------------------------
10#if defined key_agrif
11   !!----------------------------------------------------------------------
12   !!   'key_agrif'                                              AGRIF zoom
13   !!----------------------------------------------------------------------
14   !!   Agrif_tra     :
15   !!   Agrif_dyn     :
16   !!   Agrif_ssh     :
17   !!   Agrif_dyn_ts  :
18   !!   Agrif_dta_ts  :
19   !!   Agrif_ssh_ts  :
20   !!   Agrif_avm     :
21   !!   interpu       :
22   !!   interpv       :
23   !!----------------------------------------------------------------------
24   USE par_oce
25   USE oce
26   USE dom_oce     
27   USE zdf_oce
28   USE agrif_oce
29   USE phycst
30   USE dynspg_ts, ONLY: un_adv, vn_adv
31   !
32   USE in_out_manager
33   USE agrif_opa_sponge
34   USE lib_mpp
35 
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
40   PUBLIC   Agrif_tra, Agrif_avm
41   PUBLIC   interpun , interpvn
42   PUBLIC   interptsn, interpsshn, interpavm
43   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
44   PUBLIC   interpe3t, interpumsk, interpvmsk
45
46   INTEGER ::   bdy_tinterp = 0
47
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/NST 4.0 , NEMO Consortium (2017)
51   !! $Id$
52   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE Agrif_tra
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE Agrif_tra  ***
59      !!----------------------------------------------------------------------
60      !
61      IF( Agrif_Root() )   RETURN
62      !
63      Agrif_SpecialValue    = 0._wp
64      Agrif_UseSpecialValue = .TRUE.
65      !
66      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
67      !
68      Agrif_UseSpecialValue = .FALSE.
69      !
70   END SUBROUTINE Agrif_tra
71
72
73   SUBROUTINE Agrif_dyn( kt )
74      !!----------------------------------------------------------------------
75      !!                  ***  ROUTINE Agrif_DYN  ***
76      !!---------------------------------------------------------------------- 
77      INTEGER, INTENT(in) ::   kt
78      !
79      INTEGER ::   ji, jj, jk       ! dummy loop indices
80      INTEGER ::   j1, j2, i1, i2
81      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
82      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
83      !!---------------------------------------------------------------------- 
84      !
85      IF( Agrif_Root() )   RETURN
86      !
87      Agrif_SpecialValue    = 0._wp
88      Agrif_UseSpecialValue = ln_spc_dyn
89      !
90      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
91      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
92      !
93      Agrif_UseSpecialValue = .FALSE.
94      !
95      ! prevent smoothing in ghost cells
96      i1 =  1   ;   i2 = jpi
97      j1 =  1   ;   j2 = jpj
98      IF( nbondj == -1 .OR. nbondj == 2 )   j1 = 2 + nbghostcells
99      IF( nbondj == +1 .OR. nbondj == 2 )   j2 = nlcj - nbghostcells - 1
100      IF( nbondi == -1 .OR. nbondi == 2 )   i1 = 2 + nbghostcells 
101      IF( nbondi == +1 .OR. nbondi == 2 )   i2 = nlci - nbghostcells - 1
102
103      ! --- West --- !
104      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
105         ibdy1 = 2
106         ibdy2 = 1+nbghostcells 
107         !
108         ! Smoothing
109         ! ---------
110         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
111            ua_b(ibdy1:ibdy2,:) = 0._wp
112            DO jk = 1, jpkm1
113               DO jj = 1, jpj
114                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk)
115               END DO
116            END DO
117            DO jj = 1, jpj
118               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
119            END DO
120         ENDIF
121         !
122         IF( .NOT.lk_agrif_clp ) THEN
123            DO jk=1,jpkm1                 ! Smooth
124               DO jj=j1,j2
125                  ua(ibdy2,jj,jk) = 0.25_wp*(ua(ibdy2-1,jj,jk)+2._wp*ua(ibdy2,jj,jk)+ua(ibdy2+1,jj,jk))
126                  ua(ibdy2,jj,jk) = ua(ibdy2,jj,jk) * umask(ibdy2,jj,jk)
127               END DO
128            END DO
129         ENDIF
130         !
131         zub(ibdy1:ibdy2,:) = 0._wp              ! Correct transport
132         DO jk = 1, jpkm1
133            DO jj = 1, jpj
134               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) + e3u_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk)
135            END DO
136         END DO
137         DO jj=1,jpj
138            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
139         END DO
140           
141         DO jk = 1, jpkm1
142            DO jj = 1, jpj
143               ua(ibdy1:ibdy2,jj,jk) = (ua(ibdy1:ibdy2,jj,jk)+ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
144            END DO
145         END DO
146           
147         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
148            zvb(ibdy1:ibdy2,:) = 0._wp
149            DO jk = 1, jpkm1
150               DO jj = 1, jpj
151                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk)
152               END DO
153            END DO
154            DO jj = 1, jpj
155               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
156            END DO
157            DO jk = 1, jpkm1
158               DO jj = 1, jpj
159                  va(ibdy1:ibdy2,jj,jk) = (va(ibdy1:ibdy2,jj,jk)+va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
160               END DO
161            END DO
162         ENDIF
163         !
164      ENDIF
165
166      ! --- East --- !
167      IF( nbondi ==  1 .OR. nbondi == 2 ) THEN
168         ibdy1 = nlci-1-nbghostcells
169         ibdy2 = nlci-2 
170         !
171         ! Smoothing
172         ! ---------
173         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
174            ua_b(ibdy1:ibdy2,:) = 0._wp
175            DO jk = 1, jpkm1
176               DO jj = 1, jpj
177                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk)
178               END DO
179            END DO
180            DO jj = 1, jpj
181               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
182            END DO
183         ENDIF
184         !
185         IF( .NOT.lk_agrif_clp ) THEN
186            DO jk=1,jpkm1                 ! Smooth
187               DO jj=j1,j2
188                  ua(ibdy1,jj,jk) = 0.25_wp*(ua(ibdy1-1,jj,jk)+2._wp*ua(ibdy1,jj,jk)+ua(ibdy1+1,jj,jk))
189                  ua(ibdy1,jj,jk) = ua(ibdy1,jj,jk) * umask(ibdy1,jj,jk)
190               END DO
191            END DO
192         ENDIF
193         !
194         zub(ibdy1:ibdy2,:) = 0._wp              ! Correct transport
195         DO jk = 1, jpkm1
196            DO jj = 1, jpj
197               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) + e3u_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk)
198            END DO
199         END DO
200         DO jj=1,jpj
201            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
202         END DO
203           
204         DO jk = 1, jpkm1
205            DO jj = 1, jpj
206               ua(ibdy1:ibdy2,jj,jk) = (ua(ibdy1:ibdy2,jj,jk)+ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
207            END DO
208         END DO
209           
210         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
211            ibdy1 = ibdy1 + 1
212            ibdy2 = ibdy2 + 1 
213            zvb(ibdy1:ibdy2,:) = 0._wp
214            DO jk = 1, jpkm1
215               DO jj = 1, jpj
216                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk)
217               END DO
218            END DO
219            DO jj = 1, jpj
220               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
221            END DO
222            DO jk = 1, jpkm1
223               DO jj = 1, jpj
224                  va(ibdy1:ibdy2,jj,jk) = (va(ibdy1:ibdy2,jj,jk)+va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
225               END DO
226            END DO
227         ENDIF
228         !
229      ENDIF
230
231      ! --- South --- !
232      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
233         jbdy1 = 2
234         jbdy2 = 1+nbghostcells 
235         !
236         ! Smoothing
237         ! ---------
238         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
239            va_b(:,jbdy1:jbdy2) = 0._wp
240            DO jk = 1, jpkm1
241               DO ji = 1, jpi
242                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk)
243               END DO
244            END DO
245            DO ji=1,jpi
246               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
247            END DO
248         ENDIF
249         !
250         IF ( .NOT.lk_agrif_clp ) THEN
251            DO jk = 1, jpkm1              ! Smooth
252               DO ji = i1, i2
253                  va(ji,jbdy2,jk) = 0.25_wp * vmask(ji,jbdy2,jk)    &
254                     &        * ( va(ji,jbdy2-1,jk) + 2._wp*va(ji,jbdy2,jk) + va(ji,jbdy2+1,jk) )
255               END DO
256            END DO
257         ENDIF
258         !
259         zvb(:,jbdy1:jbdy2) = 0._wp              ! Correct transport
260         DO jk=1,jpkm1
261            DO ji=1,jpi
262               zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
263            END DO
264         END DO
265         DO ji = 1, jpi
266            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
267         END DO
268         DO jk = 1, jpkm1
269            DO ji = 1, jpi
270               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
271            END DO
272         END DO
273           
274         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
275            zub(:,2) = 0._wp
276            DO jk = 1, jpkm1
277               DO ji = 1, jpi
278                  zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) + e3u_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)
279               END DO
280            END DO
281            DO ji = 1, jpi
282               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
283            END DO
284               
285            DO jk = 1, jpkm1
286               DO ji = 1, jpi
287                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
288               END DO
289            END DO
290         ENDIF
291         !
292      ENDIF
293
294      ! --- North --- !
295      IF( nbondj ==  1 .OR. nbondj == 2 ) THEN
296         jbdy1 = nlcj-1-nbghostcells
297         jbdy2 = nlcj-2 
298         !
299         ! Smoothing
300         ! ---------
301         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
302            va_b(:,jbdy1:jbdy2) = 0._wp
303            DO jk = 1, jpkm1
304               DO ji = 1, jpi
305                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk)
306               END DO
307            END DO
308            DO ji=1,jpi
309               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
310            END DO
311         ENDIF
312         !
313         IF ( .NOT.lk_agrif_clp ) THEN
314            DO jk = 1, jpkm1              ! Smooth
315               DO ji = i1, i2
316                  va(ji,jbdy1,jk) = 0.25_wp * vmask(ji,jbdy1,jk)    &
317                     &        * ( va(ji,jbdy1-1,jk) + 2._wp*va(ji,jbdy1,jk) + va(ji,jbdy1+1,jk) )
318               END DO
319            END DO
320         ENDIF
321         !
322         zvb(:,jbdy1:jbdy2) = 0._wp              ! Correct transport
323         DO jk=1,jpkm1
324            DO ji=1,jpi
325               zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
326            END DO
327         END DO
328         DO ji = 1, jpi
329            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
330         END DO
331         DO jk = 1, jpkm1
332            DO ji = 1, jpi
333               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
334            END DO
335         END DO
336           
337         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
338            jbdy1 = jbdy1 + 1
339            jbdy2 = jbdy2 + 1 
340            zub(:,2) = 0._wp
341            DO jk = 1, jpkm1
342               DO ji = 1, jpi
343                  zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) + e3u_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)
344               END DO
345            END DO
346            DO ji = 1, jpi
347               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
348            END DO
349               
350            DO jk = 1, jpkm1
351               DO ji = 1, jpi
352                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
353               END DO
354            END DO
355         ENDIF
356         !
357      ENDIF
358      !
359   END SUBROUTINE Agrif_dyn
360
361
362   SUBROUTINE Agrif_dyn_ts( jn )
363      !!----------------------------------------------------------------------
364      !!                  ***  ROUTINE Agrif_dyn_ts  ***
365      !!---------------------------------------------------------------------- 
366      INTEGER, INTENT(in) ::   jn
367      !!
368      INTEGER :: ji, jj
369      !!---------------------------------------------------------------------- 
370      !
371      IF( Agrif_Root() )   RETURN
372      !
373      IF((nbondi == -1).OR.(nbondi == 2)) THEN
374         DO jj=1,jpj
375            va_e(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * hvr_e(2:nbghostcells+1,jj)
376            ! Specified fluxes:
377            ua_e(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * hur_e(2:nbghostcells+1,jj)
378            ! Characteristics method (only if ghostcells=1):
379            !alt            ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
380            !alt                       &           - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
381         END DO
382      ENDIF
383      !
384      IF((nbondi == 1).OR.(nbondi == 2)) THEN
385         DO jj=1,jpj
386            va_e(nlci-nbghostcells:nlci-1,jj)   = vbdy_e(1:nbghostcells,jj) * hvr_e(nlci-nbghostcells:nlci-1,jj)
387            ! Specified fluxes:
388            ua_e(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * hur_e(nlci-nbghostcells-1:nlci-2,jj)
389            ! Characteristics method (only if ghostcells=1):
390            !alt            ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
391            !alt                            &           + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
392         END DO
393      ENDIF
394      !
395      IF((nbondj == -1).OR.(nbondj == 2)) THEN
396         DO ji=1,jpi
397            ua_e(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * hur_e(ji,2:nbghostcells+1)
398            ! Specified fluxes:
399            va_e(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * hvr_e(ji,2:nbghostcells+1)
400            ! Characteristics method (only if ghostcells=1):
401            !alt            va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
402            !alt                       &           - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
403         END DO
404      ENDIF
405      !
406      IF((nbondj == 1).OR.(nbondj == 2)) THEN
407         DO ji=1,jpi
408            ua_e(ji,nlcj-nbghostcells:nlcj-1)   = ubdy_n(ji,1:nbghostcells) * hur_e(ji,nlcj-nbghostcells:nlcj-1)
409            ! Specified fluxes:
410            va_e(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * hvr_e(ji,nlcj-nbghostcells-1:nlcj-2)
411            ! Characteristics method (only if ghostcells=1):
412            !alt            va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2)  + va_e(ji,nlcj-3) &
413            !alt                            &           + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
414         END DO
415      ENDIF
416      !
417   END SUBROUTINE Agrif_dyn_ts
418
419
420   SUBROUTINE Agrif_dta_ts( kt )
421      !!----------------------------------------------------------------------
422      !!                  ***  ROUTINE Agrif_dta_ts  ***
423      !!---------------------------------------------------------------------- 
424      INTEGER, INTENT(in) ::   kt
425      !!
426      INTEGER :: ji, jj
427      LOGICAL :: ll_int_cons
428      !!---------------------------------------------------------------------- 
429      !
430      IF( Agrif_Root() )   RETURN
431      !
432      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
433      !
434      ! Enforce volume conservation if no time refinement: 
435      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
436      !
437      ! Interpolate barotropic fluxes
438      Agrif_SpecialValue=0._wp
439      Agrif_UseSpecialValue = ln_spc_dyn
440      !
441      IF( ll_int_cons ) THEN  ! Conservative interpolation
442         ! order matters here !!!!!!
443         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
444         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )
445         bdy_tinterp = 1
446         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
447         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )
448         bdy_tinterp = 2
449         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
450         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )         
451      ELSE ! Linear interpolation
452         bdy_tinterp = 0
453         ubdy_w(:,:) = 0._wp   ;   vbdy_w(:,:) = 0._wp 
454         ubdy_e(:,:) = 0._wp   ;   vbdy_e(:,:) = 0._wp 
455         ubdy_n(:,:) = 0._wp   ;   vbdy_n(:,:) = 0._wp 
456         ubdy_s(:,:) = 0._wp   ;   vbdy_s(:,:) = 0._wp
457         CALL Agrif_Bc_variable( unb_id, procname=interpunb )
458         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
459      ENDIF
460      Agrif_UseSpecialValue = .FALSE.
461      !
462   END SUBROUTINE Agrif_dta_ts
463
464
465   SUBROUTINE Agrif_ssh( kt )
466      !!----------------------------------------------------------------------
467      !!                  ***  ROUTINE Agrif_ssh  ***
468      !!---------------------------------------------------------------------- 
469      INTEGER, INTENT(in) ::   kt
470      !
471      INTEGER  :: ji, jj, indx, indy
472      !!---------------------------------------------------------------------- 
473      !
474      IF( Agrif_Root() )   RETURN
475      !     
476      ! Linear time interpolation of sea level
477      !
478      Agrif_SpecialValue    = 0._wp
479      Agrif_UseSpecialValue = .TRUE.
480      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
481      Agrif_UseSpecialValue = .FALSE.
482      !
483      ! --- West --- !
484      IF((nbondi == -1).OR.(nbondi == 2)) THEN
485         indx = 1+nbghostcells
486         DO jj = 1, jpj
487            DO ji = 2, indx
488               ssha(ji,jj) = hbdy_w(ji-1,jj)
489            ENDDO
490         ENDDO
491      ENDIF
492      !
493      ! --- East --- !
494      IF((nbondi == 1).OR.(nbondi == 2)) THEN
495         indx = nlci-nbghostcells
496         DO jj = 1, jpj
497            DO ji = indx, nlci-1
498               ssha(ji,jj) = hbdy_e(ji-indx+1,jj)
499            ENDDO
500         ENDDO
501      ENDIF
502      !
503      ! --- South --- !
504      IF((nbondj == -1).OR.(nbondj == 2)) THEN
505         indy = 1+nbghostcells
506         DO jj = 2, indy
507            DO ji = 1, jpi
508               ssha(ji,jj) = hbdy_s(ji,jj-1)
509            ENDDO
510         ENDDO
511      ENDIF
512      !
513      ! --- North --- !
514      IF((nbondj == 1).OR.(nbondj == 2)) THEN
515         indy = nlcj-nbghostcells
516         DO jj = indy, nlcj-1
517            DO ji = 1, jpi
518               ssha(ji,jj) = hbdy_n(ji,jj-indy+1)
519            ENDDO
520         ENDDO
521      ENDIF
522      !
523   END SUBROUTINE Agrif_ssh
524
525
526   SUBROUTINE Agrif_ssh_ts( jn )
527      !!----------------------------------------------------------------------
528      !!                  ***  ROUTINE Agrif_ssh_ts  ***
529      !!---------------------------------------------------------------------- 
530      INTEGER, INTENT(in) ::   jn
531      !!
532      INTEGER :: ji, jj, indx, indy
533      !!---------------------------------------------------------------------- 
534      !! clem ghost (starting at i,j=1 is important I think otherwise you introduce a grad(ssh)/=0 at point 2)
535      !
536      IF( Agrif_Root() )   RETURN
537      !
538      ! --- West --- !
539      IF((nbondi == -1).OR.(nbondi == 2)) THEN
540         indx = 1+nbghostcells
541         DO jj = 1, jpj
542            DO ji = 2, indx
543               ssha_e(ji,jj) = hbdy_w(ji-1,jj)
544            ENDDO
545         ENDDO
546      ENDIF
547      !
548      ! --- East --- !
549      IF((nbondi == 1).OR.(nbondi == 2)) THEN
550         indx = nlci-nbghostcells
551         DO jj = 1, jpj
552            DO ji = indx, nlci-1
553               ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj)
554            ENDDO
555         ENDDO
556      ENDIF
557      !
558      ! --- South --- !
559      IF((nbondj == -1).OR.(nbondj == 2)) THEN
560         indy = 1+nbghostcells
561         DO jj = 2, indy
562            DO ji = 1, jpi
563               ssha_e(ji,jj) = hbdy_s(ji,jj-1)
564            ENDDO
565         ENDDO
566      ENDIF
567      !
568      ! --- North --- !
569      IF((nbondj == 1).OR.(nbondj == 2)) THEN
570         indy = nlcj-nbghostcells
571         DO jj = indy, nlcj-1
572            DO ji = 1, jpi
573               ssha_e(ji,jj) = hbdy_n(ji,jj-indy+1)
574            ENDDO
575         ENDDO
576      ENDIF
577      !
578   END SUBROUTINE Agrif_ssh_ts
579
580   SUBROUTINE Agrif_avm
581      !!----------------------------------------------------------------------
582      !!                  ***  ROUTINE Agrif_avm  ***
583      !!---------------------------------------------------------------------- 
584      REAL(wp) ::   zalpha
585      !!---------------------------------------------------------------------- 
586      !
587      IF( Agrif_Root() )   RETURN
588      !
589      zalpha = 1._wp ! JC: proper time interpolation impossible 
590                     ! => use last available value from parent
591      !
592      Agrif_SpecialValue    = 0.e0
593      Agrif_UseSpecialValue = .TRUE.
594      !
595      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
596      !
597      Agrif_UseSpecialValue = .FALSE.
598      !
599   END SUBROUTINE Agrif_avm
600   
601
602   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
603      !!----------------------------------------------------------------------
604      !!                  *** ROUTINE interptsn ***
605      !!----------------------------------------------------------------------
606      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
607      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
608      LOGICAL                                     , INTENT(in   ) ::   before
609      INTEGER                                     , INTENT(in   ) ::   nb , ndir
610      !
611      INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices
612      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
613      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
614      LOGICAL :: western_side, eastern_side,northern_side,southern_side
615      ! vertical interpolation:
616      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
617      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
618      REAL(wp), DIMENSION(k1:k2) :: h_in
619      REAL(wp), DIMENSION(1:jpk) :: h_out
620      REAL(wp) :: h_diff
621
622      IF( before ) THEN         
623         DO jn = 1,jpts
624            DO jk=k1,k2
625               DO jj=j1,j2
626                 DO ji=i1,i2
627                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
628                 END DO
629              END DO
630           END DO
631        END DO
632
633# if defined key_vertical
634        DO jk=k1,k2
635           DO jj=j1,j2
636              DO ji=i1,i2
637                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
638              END DO
639           END DO
640        END DO
641# endif
642      ELSE
643
644         western_side  = (nb == 1).AND.(ndir == 1)   ;   eastern_side  = (nb == 1).AND.(ndir == 2)
645         southern_side = (nb == 2).AND.(ndir == 1)   ;   northern_side = (nb == 2).AND.(ndir == 2)
646
647# if defined key_vertical             
648         DO jj=j1,j2
649            DO ji=i1,i2
650               iref = ji
651               jref = jj
652               if(western_side) iref=MAX(2,ji)
653               if(eastern_side) iref=MIN(nlci-1,ji)
654               if(southern_side) jref=MAX(2,jj)
655               if(northern_side) jref=MIN(nlcj-1,jj)
656               N_in = 0
657               DO jk=k1,k2 !k2 = jpk of parent grid
658                  IF (ptab(ji,jj,jk,n2) == 0) EXIT
659                  N_in = N_in + 1
660                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
661                  h_in(N_in) = ptab(ji,jj,jk,n2)
662               END DO
663               N_out = 0
664               DO jk=1,jpk ! jpk of child grid
665                  IF (tmask(iref,jref,jk) == 0) EXIT
666                  N_out = N_out + 1
667                  h_out(jk) = e3t_n(iref,jref,jk)
668               ENDDO
669               IF (N_in > 0) THEN
670                  DO jn=1,jpts
671                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
672                  ENDDO
673               ENDIF
674            ENDDO
675         ENDDO
676# else
677         ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)
678# endif
679         !
680         tsa(i1:i2,j1:j2,1:jpk,1:jpts) = ptab_child(i1:i2,j1:j2,1:jpk,1:jpts)
681
682         IF ( .NOT.lk_agrif_clp ) THEN 
683            !
684            zrhox = Agrif_Rhox()
685            z1 = ( zrhox - 1. ) * 0.5
686            z3 = ( zrhox - 1. ) / ( zrhox + 1. )
687            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
688            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
689            !
690            z2 = 1. - z1
691            z4 = 1. - z3
692            z5 = 1. - z6 - z7
693            !
694            imin = i1 ; imax = i2
695            jmin = j1 ; jmax = j2
696            !
697            ! Remove CORNERS
698            IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells
699            IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1
700            IF((nbondi == -1).OR.(nbondi == 2)) imin = 1 + nbghostcells
701            IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1     
702            !
703            IF( eastern_side ) THEN
704               ibdy = nlci-nbghostcells
705               DO jn = 1, jpts
706                  tsa(ibdy+1,jmin:jmax,k1:k2,jn) = z1 * ptab_child(ibdy+1,jmin:jmax,k1:k2,jn) + z2 * ptab_child(ibdy,jmin:jmax,k1:k2,jn)
707                  DO jk = 1, jpkm1
708                     DO jj = jmin,jmax
709                        IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
710                           tsa(ibdy,jj,jk,jn) = tsa(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)
711                        ELSE
712                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy+1,jj,jk,jn)+z3*tsa(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)
713                           IF( un(ibdy-1,jj,jk) > 0._wp ) THEN
714                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy-1,jj,jk,jn)+z5*tsa(ibdy+1,jj,jk,jn) & 
715                                                 + z7*tsa(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
716                           ENDIF
717                        ENDIF
718                     END DO
719                  END DO
720                  ! Restore ghost points:
721                  tsa(ibdy+1,jmin:jmax,k1:k2,jn) = ptab_child(ibdy+1,jmin:jmax,k1:k2,jn) * tmask(ibdy+1,jmin:jmax,k1:k2)
722               END DO
723            ENDIF
724            !
725            IF( northern_side ) THEN
726               jbdy = nlcj-nbghostcells         
727               DO jn = 1, jpts
728                  tsa(imin:imax,jbdy+1,k1:k2,jn) = z1 * ptab_child(imin:imax,jbdy+1,k1:k2,jn) + z2 * ptab_child(imin:imax,jbdy,k1:k2,jn)
729                  DO jk = 1, jpkm1
730                     DO ji = imin,imax
731                        IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
732                           tsa(ji,jbdy,jk,jn) = tsa(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)
733                        ELSE
734                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy+1,jk,jn)+z3*tsa(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)       
735                           IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN
736                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy-1,jk,jn)+z5*tsa(ji,jbdy+1,jk,jn)  &
737                                                 + z7*tsa(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)
738                           ENDIF
739                        ENDIF
740                     END DO
741                  END DO
742                  ! Restore ghost points:
743                  tsa(imin:imax,jbdy+1,k1:k2,jn) = ptab_child(imin:imax,jbdy+1,k1:k2,jn) * tmask(imin:imax,jbdy+1,k1:k2)
744               END DO
745            ENDIF
746            !
747            IF( western_side ) THEN   
748               ibdy = 1+nbghostcells       
749               DO jn = 1, jpts
750                  tsa(ibdy-1,jmin:jmax,k1:k2,jn) = z1 * ptab_child(ibdy-1,jmin:jmax,k1:k2,jn) + z2 * ptab_child(ibdy,jmin:jmax,k1:k2,jn)
751                  DO jk = 1, jpkm1
752                     DO jj = jmin,jmax
753                        IF( umask(ibdy,jj,jk) == 0._wp ) THEN
754                           tsa(ibdy,jj,jk,jn) = tsa(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)
755                        ELSE
756                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy-1,jj,jk,jn)+z3*tsa(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)       
757                           IF( un(ibdy,jj,jk) < 0._wp ) THEN
758                              tsa(ibdy,jj,jk,jn)=(z6*tsa(ibdy+1,jj,jk,jn)+z5*tsa(ibdy-1,jj,jk,jn)+z7*tsa(ibdy+2,jj,jk,jn))*tmask(ibdy,jj,jk)
759                           ENDIF
760                        ENDIF
761                     END DO
762                  END DO
763                  ! Restore ghost points:
764                  tsa(ibdy-1,jmin:jmax,k1:k2,jn) = ptab_child(ibdy-1,jmin:jmax,k1:k2,jn) * tmask(ibdy-1,jmin:jmax,k1:k2)
765               END DO
766            ENDIF
767            !
768            IF( southern_side ) THEN 
769               jbdy=1+nbghostcells       
770               DO jn = 1, jpts
771                  tsa(imin:imax,jbdy-1,k1:k2,jn) = z1 * ptab_child(imin:imax,jbdy-1,k1:k2,jn) + z2 * ptab_child(imin:imax,jbdy,k1:k2,jn)
772                  DO jk = 1, jpk     
773                     DO ji=imin,imax
774                        IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
775                           tsa(ji,jbdy,jk,jn)=tsa(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)
776                        ELSE
777                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy-1,jk,jn)+z3*tsa(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)
778                           IF( vn(ji,jbdy,jk) < 0._wp ) THEN
779                              tsa(ji,jbdy,jk,jn)=(z6*tsa(ji,jbdy+1,jk,jn)+z5*tsa(ji,jbdy-1,jk,jn)+z7*tsa(ji,jbdy+2,jk,jn))*tmask(ji,jbdy,jk)
780                           ENDIF
781                        ENDIF
782                     END DO
783                  END DO
784                  ! Restore ghost points:
785                  tsa(imin:imax,jbdy-1,k1:k2,jn) = tsa(imin:imax,jbdy-1,k1:k2,jn) * tmask(imin:imax,jbdy-1,k1:k2)
786               END DO
787            ENDIF
788            !
789         ENDIF
790      ENDIF
791      !
792   END SUBROUTINE interptsn
793
794   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
795      !!----------------------------------------------------------------------
796      !!                  ***  ROUTINE interpsshn  ***
797      !!---------------------------------------------------------------------- 
798      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
799      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
800      LOGICAL                         , INTENT(in   ) ::   before
801      INTEGER                         , INTENT(in   ) ::   nb , ndir
802      !
803      LOGICAL :: western_side, eastern_side,northern_side,southern_side
804      !!---------------------------------------------------------------------- 
805      !
806      IF( before) THEN
807         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
808      ELSE
809         western_side  = (nb == 1).AND.(ndir == 1)
810         eastern_side  = (nb == 1).AND.(ndir == 2)
811         southern_side = (nb == 2).AND.(ndir == 1)
812         northern_side = (nb == 2).AND.(ndir == 2)
813         !! clem ghost
814         IF(western_side)  hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
815         IF(eastern_side)  hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
816         IF(southern_side) hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
817         IF(northern_side) hbdy_n(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
818      ENDIF
819      !
820   END SUBROUTINE interpsshn
821
822   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
823      !!----------------------------------------------------------------------
824      !!                  *** ROUTINE interpun ***
825      !!---------------------------------------------   
826      !!
827      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
828      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
829      LOGICAL, INTENT(in) :: before
830      INTEGER, INTENT(in) :: nb , ndir
831      !!
832      INTEGER :: ji,jj,jk
833      REAL(wp) :: zrhoy
834      ! vertical interpolation:
835      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
836      REAL(wp), DIMENSION(1:jpk) :: h_out
837      INTEGER  :: N_in, N_out, iref
838      REAL(wp) :: h_diff
839      LOGICAL  :: western_side, eastern_side
840      !!---------------------------------------------   
841      !
842      IF (before) THEN
843         DO jk=1,jpk
844            DO jj=j1,j2
845               DO ji=i1,i2
846                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
847# if defined key_vertical
848                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
849# endif
850               END DO
851            END DO
852         END DO
853      ELSE
854         zrhoy = Agrif_rhoy()
855# if defined key_vertical
856! VERTICAL REFINEMENT BEGIN
857         western_side  = (nb == 1).AND.(ndir == 1)
858         eastern_side  = (nb == 1).AND.(ndir == 2)
859
860         DO ji=i1,i2
861            iref = ji
862            IF (western_side) iref = MAX(2,ji)
863            IF (eastern_side) iref = MIN(nlci-2,ji)
864            DO jj=j1,j2
865               N_in = 0
866               DO jk=k1,k2
867                  IF (ptab(ji,jj,jk,2) == 0) EXIT
868                  N_in = N_in + 1
869                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
870                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
871              ENDDO
872         
873              IF (N_in == 0) THEN
874                 ua(ji,jj,:) = 0._wp
875                 CYCLE
876              ENDIF
877         
878              N_out = 0
879              DO jk=1,jpk
880                 if (umask(iref,jj,jk) == 0) EXIT
881                 N_out = N_out + 1
882                 h_out(N_out) = e3u_a(iref,jj,jk)
883              ENDDO
884         
885              IF (N_out == 0) THEN
886                 ua(ji,jj,:) = 0._wp
887                 CYCLE
888              ENDIF
889         
890              IF (N_in * N_out > 0) THEN
891                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
892! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
893                 if (h_diff < -1.e4) then
894                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
895!                    stop
896                 endif
897              ENDIF
898              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
899            ENDDO
900         ENDDO
901
902# else
903         DO jk = 1, jpkm1
904            DO jj=j1,j2
905               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
906            END DO
907         END DO
908# endif
909
910      ENDIF
911      !
912   END SUBROUTINE interpun
913
914   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
915      !!----------------------------------------------------------------------
916      !!                  *** ROUTINE interpvn ***
917      !!----------------------------------------------------------------------
918      !
919      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
920      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
921      LOGICAL, INTENT(in) :: before
922      INTEGER, INTENT(in) :: nb , ndir
923      !
924      INTEGER :: ji,jj,jk
925      REAL(wp) :: zrhox
926      ! vertical interpolation:
927      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
928      REAL(wp), DIMENSION(1:jpk) :: h_out
929      INTEGER  :: N_in, N_out, jref
930      REAL(wp) :: h_diff
931      LOGICAL  :: northern_side,southern_side
932      !!---------------------------------------------   
933      !     
934      IF (before) THEN         
935         DO jk=k1,k2
936            DO jj=j1,j2
937               DO ji=i1,i2
938                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
939# if defined key_vertical
940                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
941# endif
942               END DO
943            END DO
944         END DO
945      ELSE       
946         zrhox = Agrif_rhox()
947# if defined key_vertical
948
949         southern_side = (nb == 2).AND.(ndir == 1)
950         northern_side = (nb == 2).AND.(ndir == 2)
951
952         DO jj=j1,j2
953            jref = jj
954            IF (southern_side) jref = MAX(2,jj)
955            IF (northern_side) jref = MIN(nlcj-2,jj)
956            DO ji=i1,i2
957               N_in = 0
958               DO jk=k1,k2
959                  if (ptab(ji,jj,jk,2) == 0) EXIT
960                  N_in = N_in + 1
961                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
962                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
963               END DO
964               IF (N_in == 0) THEN
965                  va(ji,jj,:) = 0._wp
966                  CYCLE
967               ENDIF
968         
969               N_out = 0
970               DO jk=1,jpk
971                  if (vmask(ji,jref,jk) == 0) EXIT
972                  N_out = N_out + 1
973                  h_out(N_out) = e3v_a(ji,jref,jk)
974               END DO
975               IF (N_out == 0) THEN
976                 va(ji,jj,:) = 0._wp
977                 CYCLE
978               ENDIF
979               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
980            END DO
981         END DO
982# else
983         DO jk = 1, jpkm1
984            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
985         END DO
986# endif
987      ENDIF
988      !       
989   END SUBROUTINE interpvn
990
991   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
992      !!----------------------------------------------------------------------
993      !!                  ***  ROUTINE interpunb  ***
994      !!---------------------------------------------------------------------- 
995      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
996      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
997      LOGICAL                         , INTENT(in   ) ::   before
998      INTEGER                         , INTENT(in   ) ::   nb , ndir
999      !
1000      INTEGER  ::   ji, jj
1001      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1002      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1003      !!---------------------------------------------------------------------- 
1004      !
1005      IF( before ) THEN
1006         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
1007      ELSE
1008         western_side  = (nb == 1).AND.(ndir == 1)
1009         eastern_side  = (nb == 1).AND.(ndir == 2)
1010         southern_side = (nb == 2).AND.(ndir == 1)
1011         northern_side = (nb == 2).AND.(ndir == 2)
1012         zrhoy = Agrif_Rhoy()
1013         zrhot = Agrif_rhot()
1014         ! Time indexes bounds for integration
1015         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1016         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1017         ! Polynomial interpolation coefficients:
1018         IF( bdy_tinterp == 1 ) THEN
1019            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1020               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1021         ELSEIF( bdy_tinterp == 2 ) THEN
1022            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1023               &               - zt0        * (       zt0 - 1._wp)**2._wp )
1024         ELSE
1025            ztcoeff = 1
1026         ENDIF
1027         !   
1028         IF(western_side)   ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 
1029         IF(eastern_side)   ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 
1030         IF(southern_side)  ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1031         IF(northern_side)  ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 
1032         !           
1033         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1034            IF(western_side)   ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
1035            IF(eastern_side)   ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
1036            IF(southern_side)  ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
1037            IF(northern_side)  ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
1038         ENDIF
1039      ENDIF
1040      !
1041   END SUBROUTINE interpunb
1042
1043
1044   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1045      !!----------------------------------------------------------------------
1046      !!                  ***  ROUTINE interpvnb  ***
1047      !!---------------------------------------------------------------------- 
1048      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1049      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1050      LOGICAL                         , INTENT(in   ) ::   before
1051      INTEGER                         , INTENT(in   ) ::   nb , ndir
1052      !
1053      INTEGER  ::   ji,jj
1054      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1055      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1056      !!---------------------------------------------------------------------- 
1057      !
1058      IF( before ) THEN
1059         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1060      ELSE
1061         western_side  = (nb == 1).AND.(ndir == 1)
1062         eastern_side  = (nb == 1).AND.(ndir == 2)
1063         southern_side = (nb == 2).AND.(ndir == 1)
1064         northern_side = (nb == 2).AND.(ndir == 2)
1065         zrhox = Agrif_Rhox()
1066         zrhot = Agrif_rhot()
1067         ! Time indexes bounds for integration
1068         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1069         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1070         IF( bdy_tinterp == 1 ) THEN
1071            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1072               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1073         ELSEIF( bdy_tinterp == 2 ) THEN
1074            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1075               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1076         ELSE
1077            ztcoeff = 1
1078         ENDIF
1079         !! clem ghost
1080         IF(western_side)   vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 
1081         IF(eastern_side)   vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)   
1082         IF(southern_side)  vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1083         IF(northern_side)  vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 
1084         !           
1085         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1086            IF(western_side)   vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
1087            IF(eastern_side)   vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
1088            IF(southern_side)  vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
1089            IF(northern_side)  vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
1090         ENDIF
1091      ENDIF
1092      !
1093   END SUBROUTINE interpvnb
1094
1095
1096   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1097      !!----------------------------------------------------------------------
1098      !!                  ***  ROUTINE interpub2b  ***
1099      !!---------------------------------------------------------------------- 
1100      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1101      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1102      LOGICAL                         , INTENT(in   ) ::   before
1103      INTEGER                         , INTENT(in   ) ::   nb , ndir
1104      !
1105      INTEGER  ::   ji,jj
1106      REAL(wp) ::   zrhot, zt0, zt1,zat
1107      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1108      !!---------------------------------------------------------------------- 
1109      IF( before ) THEN
1110         IF ( ln_bt_fw ) THEN
1111            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1112         ELSE
1113            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1114         ENDIF
1115      ELSE
1116         western_side  = (nb == 1).AND.(ndir == 1)
1117         eastern_side  = (nb == 1).AND.(ndir == 2)
1118         southern_side = (nb == 2).AND.(ndir == 1)
1119         northern_side = (nb == 2).AND.(ndir == 2)
1120         zrhot = Agrif_rhot()
1121         ! Time indexes bounds for integration
1122         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1123         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1124         ! Polynomial interpolation coefficients:
1125         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1126            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1127         !! clem ghost
1128         IF(western_side ) ubdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1129         IF(eastern_side ) ubdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1130         IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1131         IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 
1132      ENDIF
1133      !
1134   END SUBROUTINE interpub2b
1135   
1136
1137   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1138      !!----------------------------------------------------------------------
1139      !!                  ***  ROUTINE interpvb2b  ***
1140      !!---------------------------------------------------------------------- 
1141      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1142      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1143      LOGICAL                         , INTENT(in   ) ::   before
1144      INTEGER                         , INTENT(in   ) ::   nb , ndir
1145      !
1146      INTEGER ::   ji,jj
1147      REAL(wp) ::   zrhot, zt0, zt1,zat
1148      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
1149      !!---------------------------------------------------------------------- 
1150      !
1151      IF( before ) THEN
1152         IF ( ln_bt_fw ) THEN
1153            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1154         ELSE
1155            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1156         ENDIF
1157      ELSE     
1158         western_side  = (nb == 1).AND.(ndir == 1)
1159         eastern_side  = (nb == 1).AND.(ndir == 2)
1160         southern_side = (nb == 2).AND.(ndir == 1)
1161         northern_side = (nb == 2).AND.(ndir == 2)
1162         zrhot = Agrif_rhot()
1163         ! Time indexes bounds for integration
1164         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1165         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1166         ! Polynomial interpolation coefficients:
1167         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1168            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1169         !
1170         IF(western_side )   vbdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1171         IF(eastern_side )   vbdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1172         IF(southern_side)   vbdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1173         IF(northern_side)   vbdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 
1174      ENDIF
1175      !     
1176   END SUBROUTINE interpvb2b
1177
1178
1179   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1180      !!----------------------------------------------------------------------
1181      !!                  ***  ROUTINE interpe3t  ***
1182      !!---------------------------------------------------------------------- 
1183      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1184      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1185      LOGICAL                              , INTENT(in   ) :: before
1186      INTEGER                              , INTENT(in   ) :: nb , ndir
1187      !
1188      INTEGER :: ji, jj, jk
1189      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1190      !!---------------------------------------------------------------------- 
1191      !   
1192      IF( before ) THEN
1193         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1194      ELSE
1195         western_side  = (nb == 1).AND.(ndir == 1)
1196         eastern_side  = (nb == 1).AND.(ndir == 2)
1197         southern_side = (nb == 2).AND.(ndir == 1)
1198         northern_side = (nb == 2).AND.(ndir == 2)
1199         !
1200         DO jk = k1, k2
1201            DO jj = j1, j2
1202               DO ji = i1, i2
1203                  !
1204                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1205                     IF (western_side) THEN
1206                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1207                     ELSEIF (eastern_side) THEN
1208                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1209                     ELSEIF (southern_side) THEN
1210                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1211                     ELSEIF (northern_side) THEN
1212                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1213                     ENDIF
1214                     WRITE(numout,*) '      ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1215                     kindic_agr = kindic_agr + 1
1216                  ENDIF
1217               END DO
1218            END DO
1219         END DO
1220         !
1221      ENDIF
1222      !
1223   END SUBROUTINE interpe3t
1224
1225
1226   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1227      !!----------------------------------------------------------------------
1228      !!                  ***  ROUTINE interpumsk  ***
1229      !!---------------------------------------------------------------------- 
1230      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1231      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1232      LOGICAL                              , INTENT(in   ) ::   before
1233      INTEGER                              , INTENT(in   ) ::   nb , ndir
1234      !
1235      INTEGER ::   ji, jj, jk
1236      LOGICAL ::   western_side, eastern_side   
1237      !!---------------------------------------------------------------------- 
1238      !   
1239      IF( before ) THEN
1240         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1241      ELSE
1242         western_side = (nb == 1).AND.(ndir == 1)
1243         eastern_side = (nb == 1).AND.(ndir == 2)
1244         DO jk = k1, k2
1245            DO jj = j1, j2
1246               DO ji = i1, i2
1247                   ! Velocity mask at boundary edge points:
1248                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1249                     IF (western_side) THEN
1250                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1251                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1252                        kindic_agr = kindic_agr + 1
1253                     ELSEIF (eastern_side) THEN
1254                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1255                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1256                        kindic_agr = kindic_agr + 1
1257                     ENDIF
1258                  ENDIF
1259               END DO
1260            END DO
1261         END DO
1262         !
1263      ENDIF
1264      !
1265   END SUBROUTINE interpumsk
1266
1267
1268   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1269      !!----------------------------------------------------------------------
1270      !!                  ***  ROUTINE interpvmsk  ***
1271      !!---------------------------------------------------------------------- 
1272      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1273      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1274      LOGICAL                              , INTENT(in   ) ::   before
1275      INTEGER                              , INTENT(in   ) :: nb , ndir
1276      !
1277      INTEGER ::   ji, jj, jk
1278      LOGICAL ::   northern_side, southern_side     
1279      !!---------------------------------------------------------------------- 
1280      !   
1281      IF( before ) THEN
1282         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1283      ELSE
1284         southern_side = (nb == 2).AND.(ndir == 1)
1285         northern_side = (nb == 2).AND.(ndir == 2)
1286         DO jk = k1, k2
1287            DO jj = j1, j2
1288               DO ji = i1, i2
1289                   ! Velocity mask at boundary edge points:
1290                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1291                     IF (southern_side) THEN
1292                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1293                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1294                        kindic_agr = kindic_agr + 1
1295                     ELSEIF (northern_side) THEN
1296                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1297                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1298                        kindic_agr = kindic_agr + 1
1299                     ENDIF
1300                  ENDIF
1301               END DO
1302            END DO
1303         END DO
1304         !
1305      ENDIF
1306      !
1307   END SUBROUTINE interpvmsk
1308
1309
1310   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1311      !!----------------------------------------------------------------------
1312      !!                  ***  ROUTINE interavm  ***
1313      !!---------------------------------------------------------------------- 
1314      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1315      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1316      LOGICAL                                    , INTENT(in   ) ::   before
1317      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
1318      REAL(wp), DIMENSION(1:jpk) :: h_out
1319      INTEGER  :: N_in, N_out, ji, jj, jk
1320      !!---------------------------------------------------------------------- 
1321      !     
1322      IF (before) THEN         
1323         DO jk=k1,k2
1324            DO jj=j1,j2
1325              DO ji=i1,i2
1326                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1327              END DO
1328           END DO
1329        END DO
1330#ifdef key_vertical         
1331        DO jk=k1,k2
1332           DO jj=j1,j2
1333              DO ji=i1,i2
1334                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk) 
1335              END DO
1336           END DO
1337        END DO
1338#endif
1339      ELSE 
1340#ifdef key_vertical         
1341         avm_k(i1:i2,j1:j2,1:jpk) = 0.
1342         DO jj=j1,j2
1343            DO ji=i1,i2
1344               N_in = 0
1345               DO jk=k1,k2 !k2 = jpk of parent grid
1346                  IF (ptab(ji,jj,jk,2) == 0) EXIT
1347                  N_in = N_in + 1
1348                  tabin(jk) = ptab(ji,jj,jk,1)
1349                  h_in(N_in) = ptab(ji,jj,jk,2)
1350               END DO
1351               N_out = 0
1352               DO jk=1,jpk ! jpk of child grid
1353                  IF (wmask(ji,jj,jk) == 0) EXIT
1354                  N_out = N_out + 1
1355                  h_out(jk) = e3t_n(ji,jj,jk)
1356               ENDDO
1357               IF (N_in > 0) THEN
1358                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
1359               ENDIF
1360            ENDDO
1361         ENDDO
1362#else
1363         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1364#endif
1365      ENDIF
1366      !
1367   END SUBROUTINE interpavm
1368
1369#else
1370   !!----------------------------------------------------------------------
1371   !!   Empty module                                          no AGRIF zoom
1372   !!----------------------------------------------------------------------
1373CONTAINS
1374   SUBROUTINE Agrif_OPA_Interp_empty
1375      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1376   END SUBROUTINE Agrif_OPA_Interp_empty
1377#endif
1378
1379   !!======================================================================
1380END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.