source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90 @ 11463

Last change on this file since 11463 was 11463, checked in by acc, 14 months ago

Branch: dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap. Minor bugfix in step.F90 to enable AGRIF SETTE tests to run. Also merged prettification changes to NST routines from the dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch. Neither AGRIF_DEMO nor VORTEX restart perfectly (drifting after 8 and 121 timesteps, respectively).

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