source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90 @ 10989

Last change on this file since 10989 was 10989, checked in by acc, 19 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert NST routines in preparation for getting AGRIF back up and running. AGRIF conv stage now works but requires some renaming of recently changes DIU modules (included in this commit). AGRIF compile and link stage not yet working (agrif routines need to be passed the time-level indices) but non-AGRIF SETTE tests are all OK

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