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

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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_oce_interp.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh
  • Property svn:keywords set to Id
File size: 59.4 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 licence (./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            ua_b(ibdy1:ibdy2,:) = 0._wp
110            DO jk = 1, jpkm1
111               DO jj = 1, jpj
112                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) & 
113                      & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
114               END DO
115            END DO
116            DO jj = 1, jpj
117               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * 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                  ua(ibdy2,jj,jk) = 0.25_wp*(ua(ibdy2-1,jj,jk)+2._wp*ua(ibdy2,jj,jk)+ua(ibdy2+1,jj,jk))
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_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk)*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               ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) &
143                 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
144            END DO
145         END DO
146           
147         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
148            zvb(ibdy1:ibdy2,:) = 0._wp
149            DO jk = 1, jpkm1
150               DO jj = 1, jpj
151                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 
152                     & + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * 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                  va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 
161                    & + va_b(ibdy1:ibdy2,jj)-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               ua(1,jj,jk) = 0._wp
169               va(1,jj,jk) = 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            ua_b(ibdy1:ibdy2,:) = 0._wp
181            DO jk = 1, jpkm1
182               DO jj = 1, jpj
183                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) & 
184                      & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
185               END DO
186            END DO
187            DO jj = 1, jpj
188               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * 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                  ua(ibdy1,jj,jk) = 0.25_wp*(ua(ibdy1-1,jj,jk)+2._wp*ua(ibdy1,jj,jk)+ua(ibdy1+1,jj,jk))
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_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk) * 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               ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) & 
214                 & + ua_b(ibdy1:ibdy2,jj)-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_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * 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                  va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 
234                      & + va_b(ibdy1:ibdy2,jj)-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               ua(nlci-1,jj,jk) = 0._wp
242               va(nlci  ,jj,jk) = 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            va_b(:,jbdy1:jbdy2) = 0._wp
254            DO jk = 1, jpkm1
255               DO ji = 1, jpi
256                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) & 
257                      & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
258               END DO
259            END DO
260            DO ji=1,jpi
261               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * 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                  va(ji,jbdy2,jk) = 0.25_wp*(va(ji,jbdy2-1,jk)+2._wp*va(ji,jbdy2,jk)+va(ji,jbdy2+1,jk))
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_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * 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               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 
287                 & + va_b(ji,jbdy1:jbdy2) - 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_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * 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                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 
306                    & + ua_b(ji,jbdy1:jbdy2) - 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               ua(ji,1,jk) = 0._wp
314               va(ji,1,jk) = 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            va_b(:,jbdy1:jbdy2) = 0._wp
326            DO jk = 1, jpkm1
327               DO ji = 1, jpi
328                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) & 
329                      & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
330               END DO
331            END DO
332            DO ji=1,jpi
333               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * 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                  va(ji,jbdy1,jk) = 0.25_wp*(va(ji,jbdy1-1,jk)+2._wp*va(ji,jbdy1,jk)+va(ji,jbdy1+1,jk))
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_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * 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               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 
359                 & + va_b(ji,jbdy1:jbdy2) - 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_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * 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                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 
380                    & + ua_b(ji,jbdy1:jbdy2) - 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               ua(ji,nlcj  ,jk) = 0._wp
388               va(ji,nlcj-1,jk) = 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,Naa) = hbdy_w(ji-1,jj)
523            END DO
524         END DO
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,Naa) = hbdy_e(ji-indx+1,jj)
533            END DO
534         END DO
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,Naa) = hbdy_s(ji,jj-1)
543            END DO
544         END DO
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,Naa) = hbdy_n(ji,jj-indy+1)
553            END DO
554         END DO
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            END DO
579         END DO
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            END DO
589         END DO
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            END DO
599         END DO
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            END DO
609         END DO
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) = tsn(ji,jj,jk,jn)
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_n(ji,jj,jk) 
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_n(iref,jref,jk)
702               END DO
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                  END DO
707               ENDIF
708            END DO
709         END DO
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            tsa(i1:i2,j1:j2,1:jpk,jn)=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                  tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = 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                           tsa(ibdy,jj,jk,jn) = tsa(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)
744                        ELSE
745                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy+1,jj,jk,jn)+z3*tsa(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)
746                           IF( un(ibdy-1,jj,jk) > 0._wp ) THEN
747                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy-1,jj,jk,jn)+z5*tsa(ibdy+1,jj,jk,jn) & 
748                                                 + z7*tsa(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
749                           ENDIF
750                        ENDIF
751                     END DO
752                  END DO
753                  ! Restore ghost points:
754                  tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = 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                  tsa(imin:imax,jbdy+1,1:jpkm1,jn) = 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                           tsa(ji,jbdy,jk,jn) = tsa(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)
773                        ELSE
774                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy+1,jk,jn)+z3*tsa(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)       
775                           IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN
776                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy-1,jk,jn)+z5*tsa(ji,jbdy+1,jk,jn)  &
777                                                 + z7*tsa(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)
778                           ENDIF
779                        ENDIF
780                     END DO
781                  END DO
782                  ! Restore ghost points:
783                  tsa(imin:imax,jbdy+1,1:jpkm1,jn) = 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                  tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = 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                           tsa(ibdy,jj,jk,jn) = tsa(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)
802                        ELSE
803                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy-1,jj,jk,jn)+z3*tsa(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)       
804                           IF( un(ibdy,jj,jk) < 0._wp ) THEN
805                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy+1,jj,jk,jn)+z5*tsa(ibdy-1,jj,jk,jn) &
806                                                 + z7*tsa(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
807                           ENDIF
808                        ENDIF
809                     END DO
810                  END DO
811                  ! Restore ghost points:
812                  tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = 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                  tsa(imin:imax,jbdy-1,1:jpkm1,jn) = 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                           tsa(ji,jbdy,jk,jn)=tsa(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)
831                        ELSE
832                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy-1,jk,jn)+z3*tsa(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)
833                           IF( vn(ji,jbdy,jk) < 0._wp ) THEN
834                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy+1,jk,jn)+z5*tsa(ji,jbdy-1,jk,jn) & 
835                                                 + z7*tsa(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk)
836                           ENDIF
837                        ENDIF
838                     END DO
839                  END DO
840                  ! Restore ghost points:
841                  tsa(imin:imax,jbdy-1,1:jpkm1,jn) = 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
851   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
852      !!----------------------------------------------------------------------
853      !!                  ***  ROUTINE interpsshn  ***
854      !!---------------------------------------------------------------------- 
855      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
856      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
857      LOGICAL                         , INTENT(in   ) ::   before
858      INTEGER                         , INTENT(in   ) ::   nb , ndir
859      !
860      LOGICAL :: western_side, eastern_side,northern_side,southern_side
861      !!---------------------------------------------------------------------- 
862      !
863      IF( before) THEN
864         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Nnn)
865      ELSE
866         western_side  = (nb == 1).AND.(ndir == 1)
867         eastern_side  = (nb == 1).AND.(ndir == 2)
868         southern_side = (nb == 2).AND.(ndir == 1)
869         northern_side = (nb == 2).AND.(ndir == 2)
870         !! clem ghost
871         IF(western_side )   hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
872         IF(eastern_side )   hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
873         IF(southern_side)   hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
874         IF(northern_side)   hbdy_n(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
875      ENDIF
876      !
877   END SUBROUTINE interpsshn
878
879
880   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
881      !!----------------------------------------------------------------------
882      !!                  *** ROUTINE interpun ***
883      !!---------------------------------------------   
884      !!
885      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
886      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
887      LOGICAL, INTENT(in) :: before
888      INTEGER, INTENT(in) :: nb , ndir
889      !!
890      INTEGER :: ji,jj,jk
891      REAL(wp) :: zrhoy
892      ! vertical interpolation:
893      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
894      REAL(wp), DIMENSION(1:jpk) :: h_out
895      INTEGER  :: N_in, N_out, iref
896      REAL(wp) :: h_diff
897      LOGICAL  :: western_side, eastern_side
898      !!---------------------------------------------   
899      !
900      IF (before) THEN
901         DO jk=1,jpk
902            DO jj=j1,j2
903               DO ji=i1,i2
904                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
905# if defined key_vertical
906                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
907# endif
908               END DO
909            END DO
910         END DO
911      ELSE
912         zrhoy = Agrif_rhoy()
913# if defined key_vertical
914! VERTICAL REFINEMENT BEGIN
915         western_side  = (nb == 1).AND.(ndir == 1)
916         eastern_side  = (nb == 1).AND.(ndir == 2)
917
918         DO ji=i1,i2
919            iref = ji
920            IF (western_side) iref = MAX(2,ji)
921            IF (eastern_side) iref = MIN(nlci-2,ji)
922            DO jj=j1,j2
923               N_in = 0
924               DO jk=k1,k2
925                  IF (ptab(ji,jj,jk,2) == 0) EXIT
926                  N_in = N_in + 1
927                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
928                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
929              END DO
930         
931              IF (N_in == 0) THEN
932                 ua(ji,jj,:) = 0._wp
933                 CYCLE
934              ENDIF
935         
936              N_out = 0
937              DO jk=1,jpk
938                 if (umask(iref,jj,jk) == 0) EXIT
939                 N_out = N_out + 1
940                 h_out(N_out) = e3u_a(iref,jj,jk)
941              END DO
942         
943              IF (N_out == 0) THEN
944                 ua(ji,jj,:) = 0._wp
945                 CYCLE
946              ENDIF
947         
948              IF (N_in * N_out > 0) THEN
949                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
950! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
951                 if (h_diff < -1.e4) then
952                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
953!                    stop
954                 endif
955              ENDIF
956              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
957            END DO
958         END DO
959
960# else
961         DO jk = 1, jpkm1
962            DO jj=j1,j2
963               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
964            END DO
965         END DO
966# endif
967
968      ENDIF
969      !
970   END SUBROUTINE interpun
971
972   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
973      !!----------------------------------------------------------------------
974      !!                  *** ROUTINE interpvn ***
975      !!----------------------------------------------------------------------
976      !
977      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
978      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
979      LOGICAL, INTENT(in) :: before
980      INTEGER, INTENT(in) :: nb , ndir
981      !
982      INTEGER :: ji,jj,jk
983      REAL(wp) :: zrhox
984      ! vertical interpolation:
985      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
986      REAL(wp), DIMENSION(1:jpk) :: h_out
987      INTEGER  :: N_in, N_out, jref
988      REAL(wp) :: h_diff
989      LOGICAL  :: northern_side,southern_side
990      !!---------------------------------------------   
991      !     
992      IF (before) THEN         
993         DO jk=k1,k2
994            DO jj=j1,j2
995               DO ji=i1,i2
996                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
997# if defined key_vertical
998                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
999# endif
1000               END DO
1001            END DO
1002         END DO
1003      ELSE       
1004         zrhox = Agrif_rhox()
1005# if defined key_vertical
1006
1007         southern_side = (nb == 2).AND.(ndir == 1)
1008         northern_side = (nb == 2).AND.(ndir == 2)
1009
1010         DO jj=j1,j2
1011            jref = jj
1012            IF (southern_side) jref = MAX(2,jj)
1013            IF (northern_side) jref = MIN(nlcj-2,jj)
1014            DO ji=i1,i2
1015               N_in = 0
1016               DO jk=k1,k2
1017                  if (ptab(ji,jj,jk,2) == 0) EXIT
1018                  N_in = N_in + 1
1019                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1020                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1021               END DO
1022               IF (N_in == 0) THEN
1023                  va(ji,jj,:) = 0._wp
1024                  CYCLE
1025               ENDIF
1026         
1027               N_out = 0
1028               DO jk=1,jpk
1029                  if (vmask(ji,jref,jk) == 0) EXIT
1030                  N_out = N_out + 1
1031                  h_out(N_out) = e3v_a(ji,jref,jk)
1032               END DO
1033               IF (N_out == 0) THEN
1034                 va(ji,jj,:) = 0._wp
1035                 CYCLE
1036               ENDIF
1037               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
1038            END DO
1039         END DO
1040# else
1041         DO jk = 1, jpkm1
1042            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
1043         END DO
1044# endif
1045      ENDIF
1046      !       
1047   END SUBROUTINE interpvn
1048
1049   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
1050      !!----------------------------------------------------------------------
1051      !!                  ***  ROUTINE interpunb  ***
1052      !!---------------------------------------------------------------------- 
1053      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1054      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1055      LOGICAL                         , INTENT(in   ) ::   before
1056      INTEGER                         , INTENT(in   ) ::   nb , ndir
1057      !
1058      INTEGER  ::   ji, jj
1059      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1060      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1061      !!---------------------------------------------------------------------- 
1062      !
1063      IF( before ) THEN
1064         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
1065      ELSE
1066         western_side  = (nb == 1).AND.(ndir == 1)
1067         eastern_side  = (nb == 1).AND.(ndir == 2)
1068         southern_side = (nb == 2).AND.(ndir == 1)
1069         northern_side = (nb == 2).AND.(ndir == 2)
1070         zrhoy = Agrif_Rhoy()
1071         zrhot = Agrif_rhot()
1072         ! Time indexes bounds for integration
1073         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1074         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1075         ! Polynomial interpolation coefficients:
1076         IF( bdy_tinterp == 1 ) THEN
1077            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1078               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1079         ELSEIF( bdy_tinterp == 2 ) THEN
1080            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1081               &               - zt0        * (       zt0 - 1._wp)**2._wp )
1082         ELSE
1083            ztcoeff = 1
1084         ENDIF
1085         !   
1086         IF(western_side)   ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 
1087         IF(eastern_side)   ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 
1088         IF(southern_side)  ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1089         IF(northern_side)  ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 
1090         !           
1091         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1092            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)
1093            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)
1094            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)
1095            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)
1096         ENDIF
1097      ENDIF
1098      !
1099   END SUBROUTINE interpunb
1100
1101
1102   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1103      !!----------------------------------------------------------------------
1104      !!                  ***  ROUTINE interpvnb  ***
1105      !!---------------------------------------------------------------------- 
1106      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1107      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1108      LOGICAL                         , INTENT(in   ) ::   before
1109      INTEGER                         , INTENT(in   ) ::   nb , ndir
1110      !
1111      INTEGER  ::   ji,jj
1112      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1113      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1114      !!---------------------------------------------------------------------- 
1115      !
1116      IF( before ) THEN
1117         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1118      ELSE
1119         western_side  = (nb == 1).AND.(ndir == 1)
1120         eastern_side  = (nb == 1).AND.(ndir == 2)
1121         southern_side = (nb == 2).AND.(ndir == 1)
1122         northern_side = (nb == 2).AND.(ndir == 2)
1123         zrhox = Agrif_Rhox()
1124         zrhot = Agrif_rhot()
1125         ! Time indexes bounds for integration
1126         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1127         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1128         IF( bdy_tinterp == 1 ) THEN
1129            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1130               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1131         ELSEIF( bdy_tinterp == 2 ) THEN
1132            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1133               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1134         ELSE
1135            ztcoeff = 1
1136         ENDIF
1137         !! clem ghost
1138         IF(western_side)   vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 
1139         IF(eastern_side)   vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)   
1140         IF(southern_side)  vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1141         IF(northern_side)  vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 
1142         !           
1143         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1144            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)
1145            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)
1146            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)
1147            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)
1148         ENDIF
1149      ENDIF
1150      !
1151   END SUBROUTINE interpvnb
1152
1153
1154   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1155      !!----------------------------------------------------------------------
1156      !!                  ***  ROUTINE interpub2b  ***
1157      !!---------------------------------------------------------------------- 
1158      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1159      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1160      LOGICAL                         , INTENT(in   ) ::   before
1161      INTEGER                         , INTENT(in   ) ::   nb , ndir
1162      !
1163      INTEGER  ::   ji,jj
1164      REAL(wp) ::   zrhot, zt0, zt1,zat
1165      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1166      !!---------------------------------------------------------------------- 
1167      IF( before ) THEN
1168         IF ( ln_bt_fw ) THEN
1169            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1170         ELSE
1171            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1172         ENDIF
1173      ELSE
1174         western_side  = (nb == 1).AND.(ndir == 1)
1175         eastern_side  = (nb == 1).AND.(ndir == 2)
1176         southern_side = (nb == 2).AND.(ndir == 1)
1177         northern_side = (nb == 2).AND.(ndir == 2)
1178         zrhot = Agrif_rhot()
1179         ! Time indexes bounds for integration
1180         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1181         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1182         ! Polynomial interpolation coefficients:
1183         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1184            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1185         !! clem ghost
1186         IF(western_side ) ubdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1187         IF(eastern_side ) ubdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1188         IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1189         IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 
1190      ENDIF
1191      !
1192   END SUBROUTINE interpub2b
1193   
1194
1195   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1196      !!----------------------------------------------------------------------
1197      !!                  ***  ROUTINE interpvb2b  ***
1198      !!---------------------------------------------------------------------- 
1199      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1200      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1201      LOGICAL                         , INTENT(in   ) ::   before
1202      INTEGER                         , INTENT(in   ) ::   nb , ndir
1203      !
1204      INTEGER ::   ji,jj
1205      REAL(wp) ::   zrhot, zt0, zt1,zat
1206      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
1207      !!---------------------------------------------------------------------- 
1208      !
1209      IF( before ) THEN
1210         IF ( ln_bt_fw ) THEN
1211            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1212         ELSE
1213            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1214         ENDIF
1215      ELSE     
1216         western_side  = (nb == 1).AND.(ndir == 1)
1217         eastern_side  = (nb == 1).AND.(ndir == 2)
1218         southern_side = (nb == 2).AND.(ndir == 1)
1219         northern_side = (nb == 2).AND.(ndir == 2)
1220         zrhot = Agrif_rhot()
1221         ! Time indexes bounds for integration
1222         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1223         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1224         ! Polynomial interpolation coefficients:
1225         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1226            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1227         !
1228         IF(western_side )   vbdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1229         IF(eastern_side )   vbdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1230         IF(southern_side)   vbdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1231         IF(northern_side)   vbdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 
1232      ENDIF
1233      !     
1234   END SUBROUTINE interpvb2b
1235
1236
1237   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1238      !!----------------------------------------------------------------------
1239      !!                  ***  ROUTINE interpe3t  ***
1240      !!---------------------------------------------------------------------- 
1241      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1242      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1243      LOGICAL                              , INTENT(in   ) :: before
1244      INTEGER                              , INTENT(in   ) :: nb , ndir
1245      !
1246      INTEGER :: ji, jj, jk
1247      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1248      !!---------------------------------------------------------------------- 
1249      !   
1250      IF( before ) THEN
1251         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1252      ELSE
1253         western_side  = (nb == 1).AND.(ndir == 1)
1254         eastern_side  = (nb == 1).AND.(ndir == 2)
1255         southern_side = (nb == 2).AND.(ndir == 1)
1256         northern_side = (nb == 2).AND.(ndir == 2)
1257         !
1258         DO jk = k1, k2
1259            DO jj = j1, j2
1260               DO ji = i1, i2
1261                  !
1262                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1263                     IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN
1264                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1265                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk) 
1266                        kindic_agr = kindic_agr + 1
1267                     ELSEIF (eastern_side.AND.(ptab(i2-nbghostcells+1,jj,jk)>0._wp)) THEN
1268                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1269                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1270                        kindic_agr = kindic_agr + 1
1271                     ELSEIF (southern_side.AND.(ptab(ji,j1+nbghostcells-1,jk)>0._wp)) THEN
1272                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1273                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1274                        kindic_agr = kindic_agr + 1
1275                     ELSEIF (northern_side.AND.(ptab(ji,j2-nbghostcells+1,jk)>0._wp)) THEN
1276                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1277                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1278                        kindic_agr = kindic_agr + 1
1279                     ENDIF
1280                  ENDIF
1281               END DO
1282            END DO
1283         END DO
1284         !
1285      ENDIF
1286      !
1287   END SUBROUTINE interpe3t
1288
1289
1290   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1291      !!----------------------------------------------------------------------
1292      !!                  ***  ROUTINE interpumsk  ***
1293      !!---------------------------------------------------------------------- 
1294      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1295      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1296      LOGICAL                              , INTENT(in   ) ::   before
1297      INTEGER                              , INTENT(in   ) ::   nb , ndir
1298      !
1299      INTEGER ::   ji, jj, jk
1300      LOGICAL ::   western_side, eastern_side   
1301      !!---------------------------------------------------------------------- 
1302      !   
1303      IF( before ) THEN
1304         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1305      ELSE
1306         western_side = (nb == 1).AND.(ndir == 1)
1307         eastern_side = (nb == 1).AND.(ndir == 2)
1308         DO jk = k1, k2
1309            DO jj = j1, j2
1310               DO ji = i1, i2
1311                   ! Velocity mask at boundary edge points:
1312                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1313                     IF (western_side) THEN
1314                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1315                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1316                        kindic_agr = kindic_agr + 1
1317                     ELSEIF (eastern_side) THEN
1318                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1319                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1320                        kindic_agr = kindic_agr + 1
1321                     ENDIF
1322                  ENDIF
1323               END DO
1324            END DO
1325         END DO
1326         !
1327      ENDIF
1328      !
1329   END SUBROUTINE interpumsk
1330
1331
1332   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1333      !!----------------------------------------------------------------------
1334      !!                  ***  ROUTINE interpvmsk  ***
1335      !!---------------------------------------------------------------------- 
1336      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1337      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1338      LOGICAL                              , INTENT(in   ) ::   before
1339      INTEGER                              , INTENT(in   ) :: nb , ndir
1340      !
1341      INTEGER ::   ji, jj, jk
1342      LOGICAL ::   northern_side, southern_side     
1343      !!---------------------------------------------------------------------- 
1344      !   
1345      IF( before ) THEN
1346         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1347      ELSE
1348         southern_side = (nb == 2).AND.(ndir == 1)
1349         northern_side = (nb == 2).AND.(ndir == 2)
1350         DO jk = k1, k2
1351            DO jj = j1, j2
1352               DO ji = i1, i2
1353                   ! Velocity mask at boundary edge points:
1354                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1355                     IF (southern_side) THEN
1356                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1357                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1358                        kindic_agr = kindic_agr + 1
1359                     ELSEIF (northern_side) THEN
1360                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1361                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1362                        kindic_agr = kindic_agr + 1
1363                     ENDIF
1364                  ENDIF
1365               END DO
1366            END DO
1367         END DO
1368         !
1369      ENDIF
1370      !
1371   END SUBROUTINE interpvmsk
1372
1373
1374   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1375      !!----------------------------------------------------------------------
1376      !!                  ***  ROUTINE interavm  ***
1377      !!---------------------------------------------------------------------- 
1378      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1379      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1380      LOGICAL                                    , INTENT(in   ) ::   before
1381      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
1382      REAL(wp), DIMENSION(1:jpk) :: h_out
1383      INTEGER  :: N_in, N_out, ji, jj, jk
1384      !!---------------------------------------------------------------------- 
1385      !     
1386      IF (before) THEN         
1387         DO jk=k1,k2
1388            DO jj=j1,j2
1389              DO ji=i1,i2
1390                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1391              END DO
1392           END DO
1393        END DO
1394#ifdef key_vertical         
1395        DO jk=k1,k2
1396           DO jj=j1,j2
1397              DO ji=i1,i2
1398                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk) 
1399              END DO
1400           END DO
1401        END DO
1402#endif
1403      ELSE 
1404#ifdef key_vertical         
1405         avm_k(i1:i2,j1:j2,1:jpk) = 0.
1406         DO jj=j1,j2
1407            DO ji=i1,i2
1408               N_in = 0
1409               DO jk=k1,k2 !k2 = jpk of parent grid
1410                  IF (ptab(ji,jj,jk,2) == 0) EXIT
1411                  N_in = N_in + 1
1412                  tabin(jk) = ptab(ji,jj,jk,1)
1413                  h_in(N_in) = ptab(ji,jj,jk,2)
1414               END DO
1415               N_out = 0
1416               DO jk=1,jpk ! jpk of child grid
1417                  IF (wmask(ji,jj,jk) == 0) EXIT
1418                  N_out = N_out + 1
1419                  h_out(jk) = e3t_n(ji,jj,jk)
1420               END DO
1421               IF (N_in > 0) THEN
1422                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
1423               ENDIF
1424            END DO
1425         END DO
1426#else
1427         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1428#endif
1429      ENDIF
1430      !
1431   END SUBROUTINE interpavm
1432
1433#else
1434   !!----------------------------------------------------------------------
1435   !!   Empty module                                          no AGRIF zoom
1436   !!----------------------------------------------------------------------
1437CONTAINS
1438   SUBROUTINE Agrif_OCE_Interp_empty
1439      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1440   END SUBROUTINE Agrif_OCE_Interp_empty
1441#endif
1442
1443   !!======================================================================
1444END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.