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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_opa_interp.F90 in branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 6401

Last change on this file since 6401 was 6401, checked in by timgraham, 8 years ago

Reinstate svn keywords before update to head of trunk

  • Property svn:keywords set to Id
File size: 56.0 KB
Line 
1MODULE agrif_opa_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_opa_interp  ***
4   !! AGRIF: interpolation package
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (XXX)  Original cade
7   !!             -   !  2005-11  (XXX)
8   !!            3.2  !  2009-04  (R. Benshila)
9   !!            3.6  !  2014-09  (R. Benshila)
10   !!----------------------------------------------------------------------
11#if defined key_agrif && ! defined key_offline
12   !!----------------------------------------------------------------------
13   !!   'key_agrif'                                              AGRIF zoom
14   !!   NOT 'key_offline'                               NO off-line tracers
15   !!----------------------------------------------------------------------
16   !!   Agrif_tra     :
17   !!   Agrif_dyn     :
18   !!   interpu       :
19   !!   interpv       :
20   !!----------------------------------------------------------------------
21   USE par_oce
22   USE oce
23   USE dom_oce     
24   USE sol_oce
25   USE agrif_oce
26   USE phycst
27   USE in_out_manager
28   USE agrif_opa_sponge
29   USE lib_mpp
30   USE wrk_nemo
31   USE dynspg_oce
32   USE zdf_oce
33 
34   IMPLICIT NONE
35   PRIVATE
36
37   INTEGER :: bdy_tinterp = 0
38
39   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
40! VERTICAL REFINEMENT BEGIN   
41   PUBLIC   Agrif_Init_InterpScales
42! VERTICAL REFINEMENT END
43   PUBLIC   interpun, interpvn, interpun2d, interpvn2d 
44   PUBLIC   interptsn,  interpsshn
45   PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b
46   PUBLIC   interpe3t, interpumsk, interpvmsk
47# if defined key_zdftke
48   PUBLIC   Agrif_tke, interpavm
49# endif
50
51#  include "domzgr_substitute.h90" 
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/NST 3.6 , NEMO Consortium (2010)
55   !! $Id$
56   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58
59! VERTICAL REFINEMENT BEGIN
60   REAL, DIMENSION(:,:,:), ALLOCATABLE :: interp_scales_t, interp_scales_u, interp_scales_v
61!$AGRIF_DO_NOT_TREAT
62   LOGICAL :: scaleT, scaleU, scaleV = .FALSE.
63!$AGRIF_END_DO_NOT_TREAT
64! VERTICAL REFINEMENT END
65
66CONTAINS
67
68! VERTICAL REFINEMENT BEGIN
69
70   SUBROUTINE Agrif_Init_InterpScales
71
72    scaleT = .TRUE.
73    Call Agrif_Bc_Variable(scales_t_id,calledweight=1.,procname=interpscales)
74    scaleT = .FALSE.
75   
76    scaleU = .TRUE.
77    Call Agrif_Bc_Variable(scales_u_id,calledweight=1.,procname=interpscales)
78    scaleU = .FALSE.
79
80    scaleV = .TRUE.
81    Call Agrif_Bc_Variable(scales_v_id,calledweight=1.,procname=interpscales)
82    scaleV = .FALSE.
83
84   END SUBROUTINE Agrif_Init_InterpScales
85   
86   SUBROUTINE interpscales(ptab,i1,i2,j1,j2,k1,k2,before)
87      !!---------------------------------------------
88      !!   *** ROUTINE interpscales ***
89      !!---------------------------------------------
90     
91      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
92      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
93
94      INTEGER :: ji, jj, jk
95      LOGICAL :: before
96
97      IF (before) THEN
98      IF (scaleT ) THEN
99      DO jk=k1,k2
100         DO jj=j1,j2
101            DO ji=i1,i2
102!               ptab(ji,jj,jk) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
103               ptab(ji,jj,jk) = fse3t_n(ji,jj,jk)
104            END DO
105         END DO
106      END DO
107      ELSEIF (scaleU) THEN
108      DO jk=k1,k2
109         DO jj=j1,j2
110            DO ji=i1,i2
111!               ptab(ji,jj,jk) = fse3u_n(ji,jj,jk) * umask(ji,jj,jk)
112!               ptab(ji,jj,jk) = fse3u_n(ji,jj,jk)
113                ptab(ji,jj,jk) = umask(ji,jj,jk)
114            END DO
115         END DO
116      END DO
117      ELSEIF (scaleV) THEN
118      DO jk=k1,k2
119         DO jj=j1,j2
120            DO ji=i1,i2
121!               ptab(ji,jj,jk) = fse3v_n(ji,jj,jk) * vmask(ji,jj,jk)
122!               ptab(ji,jj,jk) = fse3v_n(ji,jj,jk)
123               ptab(ji,jj,jk) = vmask(ji,jj,jk)
124            END DO
125         END DO
126      END DO
127      ENDIF
128      ELSE
129      IF (scaleT ) THEN
130      IF (.not.allocated(interp_scales_t)) allocate(interp_scales_t(jpi,jpj,k1:k2))
131      DO jk=k1,k2
132         DO jj=j1,j2
133            DO ji=i1,i2
134               interp_scales_t(ji,jj,jk) = ptab(ji,jj,jk)
135            END DO
136         END DO
137      END DO
138      ELSEIF (scaleU) THEN
139      IF (.not.allocated(interp_scales_u)) allocate(interp_scales_u(jpi,jpj,k1:k2))
140      DO jk=k1,k2
141         DO jj=j1,j2
142            DO ji=i1,i2
143               interp_scales_u(ji,jj,jk) = ptab(ji,jj,jk)
144            END DO
145         END DO
146      END DO
147      ELSEIF (scaleV) THEN
148      IF (.not.allocated(interp_scales_v)) allocate(interp_scales_v(jpi,jpj,k1:k2))
149      DO jk=k1,k2
150         DO jj=j1,j2
151            DO ji=i1,i2
152               interp_scales_v(ji,jj,jk) = ptab(ji,jj,jk)
153            END DO
154         END DO
155      END DO
156      ENDIF
157      ENDIF
158
159   END SUBROUTINE interpscales
160
161! VERTICAL REFINEMENT END
162
163   SUBROUTINE Agrif_tra
164      !!----------------------------------------------------------------------
165      !!                  ***  ROUTINE Agrif_tra  ***
166      !!----------------------------------------------------------------------
167      !
168      IF( Agrif_Root() )   RETURN
169
170      Agrif_SpecialValue    = 0.e0
171      Agrif_UseSpecialValue = .TRUE.
172
173      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
174      Agrif_UseSpecialValue = .FALSE.
175      !
176   END SUBROUTINE Agrif_tra
177
178
179   SUBROUTINE Agrif_dyn( kt )
180      !!----------------------------------------------------------------------
181      !!                  ***  ROUTINE Agrif_DYN  ***
182      !!---------------------------------------------------------------------- 
183      !!
184      INTEGER, INTENT(in) ::   kt
185      !!
186      INTEGER :: ji,jj,jk, j1,j2, i1,i2
187      REAL(wp) :: timeref
188      REAL(wp) :: z2dt, znugdt
189      REAL(wp) :: zrhox, zrhoy
190      REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1
191      !!---------------------------------------------------------------------- 
192
193      IF( Agrif_Root() )   RETURN
194
195      CALL wrk_alloc( jpi, jpj, spgv1, spgu1 )
196
197      Agrif_SpecialValue=0.
198      Agrif_UseSpecialValue = ln_spc_dyn
199
200      CALL Agrif_Bc_variable(un_interp_id,procname=interpun)
201      CALL Agrif_Bc_variable(vn_interp_id,procname=interpvn)
202
203#if defined key_dynspg_flt
204      CALL Agrif_Bc_variable(e1u_id,calledweight=1., procname=interpun2d)
205      CALL Agrif_Bc_variable(e2v_id,calledweight=1., procname=interpvn2d)
206#endif
207
208      Agrif_UseSpecialValue = .FALSE.
209
210      zrhox = Agrif_Rhox()
211      zrhoy = Agrif_Rhoy()
212
213      timeref = 1.
214      ! time step: leap-frog
215      z2dt = 2. * rdt
216      ! time step: Euler if restart from rest
217      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
218      ! coefficients
219      znugdt =  grav * z2dt   
220
221      ! prevent smoothing in ghost cells
222      i1=1
223      i2=jpi
224      j1=1
225      j2=jpj
226      IF((nbondj == -1).OR.(nbondj == 2)) j1 = 3
227      IF((nbondj == +1).OR.(nbondj == 2)) j2 = nlcj-2
228      IF((nbondi == -1).OR.(nbondi == 2)) i1 = 3
229      IF((nbondi == +1).OR.(nbondi == 2)) i2 = nlci-2
230
231
232      IF((nbondi == -1).OR.(nbondi == 2)) THEN
233#if defined key_dynspg_flt
234         DO jk=1,jpkm1
235            DO jj=j1,j2
236               ua(2,jj,jk) = (ua(2,jj,jk) - z2dt * znugdt * laplacu(2,jj))*umask(2,jj,jk)
237            END DO
238         END DO
239
240         spgu(2,:)=0.
241
242         DO jk=1,jpkm1
243            DO jj=1,jpj
244               spgu(2,jj)=spgu(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
245            END DO
246         END DO
247
248         DO jj=1,jpj
249            IF (umask(2,jj,1).NE.0.) THEN
250               spgu(2,jj)=spgu(2,jj)/hu(2,jj)
251            ENDIF
252         END DO
253#else
254         spgu(2,:) = ua_b(2,:)
255#endif
256
257         DO jk=1,jpkm1
258            DO jj=j1,j2
259               ua(2,jj,jk) = 0.25*(ua(1,jj,jk)+2.*ua(2,jj,jk)+ua(3,jj,jk))
260               ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk)
261            END DO
262         END DO
263
264         spgu1(2,:)=0.
265
266         DO jk=1,jpkm1
267            DO jj=1,jpj
268               spgu1(2,jj)=spgu1(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
269            END DO
270         END DO
271
272         DO jj=1,jpj
273            IF (umask(2,jj,1).NE.0.) THEN
274               spgu1(2,jj)=spgu1(2,jj)/hu(2,jj)
275            ENDIF
276         END DO
277
278         DO jk=1,jpkm1
279            DO jj=j1,j2
280               ua(2,jj,jk) = (ua(2,jj,jk)+spgu(2,jj)-spgu1(2,jj))*umask(2,jj,jk)
281            END DO
282         END DO
283
284#if defined key_dynspg_ts
285         ! Set tangential velocities to time splitting estimate
286         spgv1(2,:)=0.
287         DO jk=1,jpkm1
288            DO jj=1,jpj
289               spgv1(2,jj)=spgv1(2,jj)+fse3v_a(2,jj,jk)*va(2,jj,jk)
290            END DO
291         END DO
292         DO jj=1,jpj
293            spgv1(2,jj)=spgv1(2,jj)*hvr_a(2,jj)
294         END DO
295         DO jk=1,jpkm1
296            DO jj=1,jpj
297               va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-spgv1(2,jj))*vmask(2,jj,jk)
298            END DO
299         END DO
300#endif
301
302      ENDIF
303
304      IF((nbondi == 1).OR.(nbondi == 2)) THEN
305#if defined key_dynspg_flt
306         DO jk=1,jpkm1
307            DO jj=j1,j2
308               ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)- z2dt * znugdt * laplacu(nlci-2,jj))*umask(nlci-2,jj,jk)
309            END DO
310         END DO
311         spgu(nlci-2,:)=0.
312         DO jk=1,jpkm1
313            DO jj=1,jpj
314               spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)
315            ENDDO
316         ENDDO
317         DO jj=1,jpj
318            IF (umask(nlci-2,jj,1).NE.0.) THEN
319               spgu(nlci-2,jj)=spgu(nlci-2,jj)/hu(nlci-2,jj)
320            ENDIF
321         END DO
322#else
323         spgu(nlci-2,:) = ua_b(nlci-2,:)
324#endif
325         DO jk=1,jpkm1
326            DO jj=j1,j2
327               ua(nlci-2,jj,jk) = 0.25*(ua(nlci-3,jj,jk)+2.*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk))
328
329               ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk)
330
331            END DO
332         END DO
333         spgu1(nlci-2,:)=0.
334         DO jk=1,jpkm1
335            DO jj=1,jpj
336               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk)
337            END DO
338         END DO
339         DO jj=1,jpj
340            IF (umask(nlci-2,jj,1).NE.0.) THEN
341               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)/hu(nlci-2,jj)
342            ENDIF
343         END DO
344         DO jk=1,jpkm1
345            DO jj=j1,j2
346               ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+spgu(nlci-2,jj)-spgu1(nlci-2,jj))*umask(nlci-2,jj,jk)
347            END DO
348         END DO
349
350#if defined key_dynspg_ts
351         ! Set tangential velocities to time splitting estimate
352         spgv1(nlci-1,:)=0._wp
353         DO jk=1,jpkm1
354            DO jj=1,jpj
355               spgv1(nlci-1,jj)=spgv1(nlci-1,jj)+fse3v_a(nlci-1,jj,jk)*va(nlci-1,jj,jk)*vmask(nlci-1,jj,jk)
356            END DO
357         END DO
358
359         DO jj=1,jpj
360            spgv1(nlci-1,jj)=spgv1(nlci-1,jj)*hvr_a(nlci-1,jj)
361         END DO
362
363         DO jk=1,jpkm1
364            DO jj=1,jpj
365               va(nlci-1,jj,jk) = (va(nlci-1,jj,jk)+va_b(nlci-1,jj)-spgv1(nlci-1,jj))*vmask(nlci-1,jj,jk)
366            END DO
367         END DO
368#endif
369
370      ENDIF
371
372      IF((nbondj == -1).OR.(nbondj == 2)) THEN
373
374#if defined key_dynspg_flt
375         DO jk=1,jpkm1
376            DO ji=1,jpi
377               va(ji,2,jk) = (va(ji,2,jk) - z2dt * znugdt * laplacv(ji,2))*vmask(ji,2,jk)
378            END DO
379         END DO
380
381         spgv(:,2)=0.
382
383         DO jk=1,jpkm1
384            DO ji=1,jpi
385               spgv(ji,2)=spgv(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)
386            END DO
387         END DO
388
389         DO ji=1,jpi
390            IF (vmask(ji,2,1).NE.0.) THEN
391               spgv(ji,2)=spgv(ji,2)/hv(ji,2)
392            ENDIF
393         END DO
394#else
395         spgv(:,2)=va_b(:,2)
396#endif
397
398         DO jk=1,jpkm1
399            DO ji=i1,i2
400               va(ji,2,jk)=0.25*(va(ji,1,jk)+2.*va(ji,2,jk)+va(ji,3,jk))
401               va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk)
402            END DO
403         END DO
404
405         spgv1(:,2)=0.
406
407         DO jk=1,jpkm1
408            DO ji=1,jpi
409               spgv1(ji,2)=spgv1(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk)
410            END DO
411         END DO
412
413         DO ji=1,jpi
414            IF (vmask(ji,2,1).NE.0.) THEN
415               spgv1(ji,2)=spgv1(ji,2)/hv(ji,2)
416            ENDIF
417         END DO
418
419         DO jk=1,jpkm1
420            DO ji=1,jpi
421               va(ji,2,jk) = (va(ji,2,jk)+spgv(ji,2)-spgv1(ji,2))*vmask(ji,2,jk)
422            END DO
423         END DO
424
425#if defined key_dynspg_ts
426         ! Set tangential velocities to time splitting estimate
427         spgu1(:,2)=0._wp
428         DO jk=1,jpkm1
429            DO ji=1,jpi
430               spgu1(ji,2)=spgu1(ji,2)+fse3u_a(ji,2,jk)*ua(ji,2,jk)*umask(ji,2,jk)
431            END DO
432         END DO
433
434         DO ji=1,jpi
435            spgu1(ji,2)=spgu1(ji,2)*hur_a(ji,2)
436         END DO
437
438         DO jk=1,jpkm1
439            DO ji=1,jpi
440               ua(ji,2,jk) = (ua(ji,2,jk)+ua_b(ji,2)-spgu1(ji,2))*umask(ji,2,jk)
441            END DO
442         END DO
443#endif
444      ENDIF
445
446      IF((nbondj == 1).OR.(nbondj == 2)) THEN
447
448#if defined key_dynspg_flt
449         DO jk=1,jpkm1
450            DO ji=1,jpi
451               va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)-z2dt * znugdt * laplacv(ji,nlcj-2))*vmask(ji,nlcj-2,jk)
452            END DO
453         END DO
454
455
456         spgv(:,nlcj-2)=0.
457
458         DO jk=1,jpkm1
459            DO ji=1,jpi
460               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
461            END DO
462         END DO
463
464         DO ji=1,jpi
465            IF (vmask(ji,nlcj-2,1).NE.0.) THEN
466               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)/hv(ji,nlcj-2)
467            ENDIF
468         END DO
469
470#else
471         spgv(:,nlcj-2)=va_b(:,nlcj-2)
472#endif
473
474         DO jk=1,jpkm1
475            DO ji=i1,i2
476               va(ji,nlcj-2,jk)=0.25*(va(ji,nlcj-3,jk)+2.*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk))
477               va(ji,nlcj-2,jk) = va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
478            END DO
479         END DO
480
481         spgv1(:,nlcj-2)=0.
482
483         DO jk=1,jpkm1
484            DO ji=1,jpi
485               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
486            END DO
487         END DO
488
489         DO ji=1,jpi
490            IF (vmask(ji,nlcj-2,1).NE.0.) THEN
491               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)/hv(ji,nlcj-2)
492            ENDIF
493         END DO
494
495         DO jk=1,jpkm1
496            DO ji=1,jpi
497               va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+spgv(ji,nlcj-2)-spgv1(ji,nlcj-2))*vmask(ji,nlcj-2,jk)
498            END DO
499         END DO
500
501#if defined key_dynspg_ts
502         ! Set tangential velocities to time splitting estimate
503         spgu1(:,nlcj-1)=0._wp
504         DO jk=1,jpkm1
505            DO ji=1,jpi
506               spgu1(ji,nlcj-1)=spgu1(ji,nlcj-1)+fse3u_a(ji,nlcj-1,jk)*ua(ji,nlcj-1,jk)
507            END DO
508         END DO
509
510         DO ji=1,jpi
511            spgu1(ji,nlcj-1)=spgu1(ji,nlcj-1)*hur_a(ji,nlcj-1)
512         END DO
513
514         DO jk=1,jpkm1
515            DO ji=1,jpi
516               ua(ji,nlcj-1,jk) = (ua(ji,nlcj-1,jk)+ua_b(ji,nlcj-1)-spgu1(ji,nlcj-1))*umask(ji,nlcj-1,jk)
517            END DO
518         END DO
519#endif
520
521      ENDIF
522      !
523      CALL wrk_dealloc( jpi, jpj, spgv1, spgu1 )
524      !
525   END SUBROUTINE Agrif_dyn
526
527   SUBROUTINE Agrif_dyn_ts( jn )
528      !!----------------------------------------------------------------------
529      !!                  ***  ROUTINE Agrif_dyn_ts  ***
530      !!---------------------------------------------------------------------- 
531      !!
532      INTEGER, INTENT(in) ::   jn
533      !!
534      INTEGER :: ji, jj
535      !!---------------------------------------------------------------------- 
536
537      IF( Agrif_Root() )   RETURN
538
539      IF((nbondi == -1).OR.(nbondi == 2)) THEN
540         DO jj=1,jpj
541            va_e(2,jj) = vbdy_w(jj) * hvr_e(2,jj)
542            ! Specified fluxes:
543            ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj)
544            ! Characteristics method:
545            !alt            ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
546            !alt                       &           - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
547         END DO
548      ENDIF
549
550      IF((nbondi == 1).OR.(nbondi == 2)) THEN
551         DO jj=1,jpj
552            va_e(nlci-1,jj) = vbdy_e(jj) * hvr_e(nlci-1,jj)
553            ! Specified fluxes:
554            ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj)
555            ! Characteristics method:
556            !alt            ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
557            !alt                            &           + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
558         END DO
559      ENDIF
560
561      IF((nbondj == -1).OR.(nbondj == 2)) THEN
562         DO ji=1,jpi
563            ua_e(ji,2) = ubdy_s(ji) * hur_e(ji,2)
564            ! Specified fluxes:
565            va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2)
566            ! Characteristics method:
567            !alt            va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
568            !alt                       &           - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
569         END DO
570      ENDIF
571
572      IF((nbondj == 1).OR.(nbondj == 2)) THEN
573         DO ji=1,jpi
574            ua_e(ji,nlcj-1) = ubdy_n(ji) * hur_e(ji,nlcj-1)
575            ! Specified fluxes:
576            va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2)
577            ! Characteristics method:
578            !alt            va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2)  + va_e(ji,nlcj-3) &
579            !alt                            &           + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
580         END DO
581      ENDIF
582      !
583   END SUBROUTINE Agrif_dyn_ts
584
585   SUBROUTINE Agrif_dta_ts( kt )
586      !!----------------------------------------------------------------------
587      !!                  ***  ROUTINE Agrif_dta_ts  ***
588      !!---------------------------------------------------------------------- 
589      !!
590      INTEGER, INTENT(in) ::   kt
591      !!
592      INTEGER :: ji, jj
593      LOGICAL :: ll_int_cons
594      REAL(wp) :: zrhot, zt
595      !!---------------------------------------------------------------------- 
596
597      IF( Agrif_Root() )   RETURN
598
599      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in
600      ! the forward case only
601
602      zrhot = Agrif_rhot()
603
604      ! "Central" time index for interpolation:
605      IF (ln_bt_fw) THEN
606         zt = REAL(Agrif_NbStepint()+0.5_wp,wp) / zrhot
607      ELSE
608         zt = REAL(Agrif_NbStepint(),wp) / zrhot
609      ENDIF
610
611      ! Linear interpolation of sea level
612      Agrif_SpecialValue    = 0.e0
613      Agrif_UseSpecialValue = .TRUE.
614      CALL Agrif_Bc_variable(sshn_id,calledweight=zt, procname=interpsshn )
615      Agrif_UseSpecialValue = .FALSE.
616
617      ! Interpolate barotropic fluxes
618      Agrif_SpecialValue=0.
619      Agrif_UseSpecialValue = ln_spc_dyn
620
621      IF (ll_int_cons) THEN ! Conservative interpolation
622         ! orders matters here !!!!!!
623         CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1._wp, procname=interpub2b) ! Time integrated
624         CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1._wp, procname=interpvb2b)
625         bdy_tinterp = 1
626         CALL Agrif_Bc_variable(unb_id ,calledweight=1._wp, procname=interpunb) ! After
627         CALL Agrif_Bc_variable(vnb_id ,calledweight=1._wp, procname=interpvnb)
628         bdy_tinterp = 2
629         CALL Agrif_Bc_variable(unb_id ,calledweight=0._wp, procname=interpunb) ! Before
630         CALL Agrif_Bc_variable(vnb_id ,calledweight=0._wp, procname=interpvnb)         
631      ELSE ! Linear interpolation
632         bdy_tinterp = 0
633         ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 
634         ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 
635         ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 
636         ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 
637         CALL Agrif_Bc_variable(unb_id,calledweight=zt, procname=interpunb)
638         CALL Agrif_Bc_variable(vnb_id,calledweight=zt, procname=interpvnb)
639      ENDIF
640      Agrif_UseSpecialValue = .FALSE.
641      !
642   END SUBROUTINE Agrif_dta_ts
643
644   SUBROUTINE Agrif_ssh( kt )
645      !!----------------------------------------------------------------------
646      !!                  ***  ROUTINE Agrif_DYN  ***
647      !!---------------------------------------------------------------------- 
648      INTEGER, INTENT(in) ::   kt
649      !!
650      !!---------------------------------------------------------------------- 
651
652      IF( Agrif_Root() )   RETURN
653
654      IF((nbondi == -1).OR.(nbondi == 2)) THEN
655         ssha(2,:)=ssha(3,:)
656         sshn(2,:)=sshn(3,:)
657      ENDIF
658
659      IF((nbondi == 1).OR.(nbondi == 2)) THEN
660         ssha(nlci-1,:)=ssha(nlci-2,:)
661         sshn(nlci-1,:)=sshn(nlci-2,:)
662      ENDIF
663
664      IF((nbondj == -1).OR.(nbondj == 2)) THEN
665         ssha(:,2)=ssha(:,3)
666         sshn(:,2)=sshn(:,3)
667      ENDIF
668
669      IF((nbondj == 1).OR.(nbondj == 2)) THEN
670         ssha(:,nlcj-1)=ssha(:,nlcj-2)
671         sshn(:,nlcj-1)=sshn(:,nlcj-2)
672      ENDIF
673
674   END SUBROUTINE Agrif_ssh
675
676   SUBROUTINE Agrif_ssh_ts( jn )
677      !!----------------------------------------------------------------------
678      !!                  ***  ROUTINE Agrif_ssh_ts  ***
679      !!---------------------------------------------------------------------- 
680      INTEGER, INTENT(in) ::   jn
681      !!
682      INTEGER :: ji,jj
683      !!---------------------------------------------------------------------- 
684
685      IF((nbondi == -1).OR.(nbondi == 2)) THEN
686         DO jj=1,jpj
687            ssha_e(2,jj) = hbdy_w(jj)
688         END DO
689      ENDIF
690
691      IF((nbondi == 1).OR.(nbondi == 2)) THEN
692         DO jj=1,jpj
693            ssha_e(nlci-1,jj) = hbdy_e(jj)
694         END DO
695      ENDIF
696
697      IF((nbondj == -1).OR.(nbondj == 2)) THEN
698         DO ji=1,jpi
699            ssha_e(ji,2) = hbdy_s(ji)
700         END DO
701      ENDIF
702
703      IF((nbondj == 1).OR.(nbondj == 2)) THEN
704         DO ji=1,jpi
705            ssha_e(ji,nlcj-1) = hbdy_n(ji)
706         END DO
707      ENDIF
708
709   END SUBROUTINE Agrif_ssh_ts
710
711# if defined key_zdftke
712   SUBROUTINE Agrif_tke
713      !!----------------------------------------------------------------------
714      !!                  ***  ROUTINE Agrif_tke  ***
715      !!---------------------------------------------------------------------- 
716      REAL(wp) ::   zalpha
717      !
718      return
719     
720      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
721      IF( zalpha > 1. )   zalpha = 1.
722     
723      Agrif_SpecialValue    = 0.e0
724      Agrif_UseSpecialValue = .TRUE.
725     
726      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
727             
728      Agrif_UseSpecialValue = .FALSE.
729      !
730   END SUBROUTINE Agrif_tke
731# endif
732
733   SUBROUTINE interptsn(ptab,i1,i2,j1,j2,k1,k2,n1,n2,before,nb,ndir)
734      !!---------------------------------------------
735      !!   *** ROUTINE interptsn ***
736      !!---------------------------------------------
737      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
738      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
739      LOGICAL, INTENT(in) :: before
740      INTEGER, INTENT(in) :: nb , ndir
741      !
742      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
743      INTEGER :: imin, imax, jmin, jmax
744      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
745      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
746      LOGICAL :: western_side, eastern_side,northern_side,southern_side
747! VERTICAL REFINEMENT BEGIN
748      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
749      REAL(wp) :: h_in(k1:k2)
750      REAL(wp) :: h_out(1:jpk)
751      INTEGER :: N_in, N_out
752      REAL(wp) :: h_diff
753! VERTICAL REFINEMENT END
754
755      IF (before) THEN         
756         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
757      ELSE 
758! VERTICAL REFINEMENT BEGIN
759
760         ptab_child(:,:,:,:) = 0.
761         do jj=j1,j2
762         do ji=i1,i2
763           h_in(k1:k2) = interp_scales_t(ji,jj,k1:k2)
764           h_out(1:jpk) = fse3t(ji,jj,1:jpk)
765           h_diff = sum(h_out(1:jpk-1))-sum(h_in(k1:k2-1))
766           N_in = k2-1
767           N_out = jpk-1
768           if (h_diff > 0) then
769             h_in(N_in+1) = h_diff
770             N_in = N_in + 1
771           else
772             h_out(N_out+1) = -h_diff
773             N_out = N_out + 1
774           endif
775           ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:)
776           do jn=1,jpts
777             call reconstructandremap(ptab(ji,jj,1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
778           enddo
779!           if (abs(h_diff) > 1000.) then
780!           do jn=1,jpts
781!             do jk=1,N_out
782!             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn)
783!             enddo
784!           enddo
785!         endif
786         enddo
787         enddo
788
789! VERTICAL REFINEMENT END
790
791         !
792         western_side  = (nb == 1).AND.(ndir == 1)
793         eastern_side  = (nb == 1).AND.(ndir == 2)
794         southern_side = (nb == 2).AND.(ndir == 1)
795         northern_side = (nb == 2).AND.(ndir == 2)
796         !
797         zrhox = Agrif_Rhox()
798         !
799         zalpha1 = ( zrhox - 1. ) * 0.5
800         zalpha2 = 1. - zalpha1
801         !
802         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
803         zalpha4 = 1. - zalpha3
804         !
805         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
806         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
807         zalpha5 = 1. - zalpha6 - zalpha7
808         !
809         imin = i1
810         imax = i2
811         jmin = j1
812         jmax = j2
813         !
814         ! Remove CORNERS
815         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
816         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
817         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
818         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2       
819         !
820! VERTICAL REFINEMENT BEGIN
821
822! WARNING :
823! ptab replaced by ptab_child in the following
824! k1 k2 replaced by 1 jpk
825! VERTICAL REFINEMENT END
826         IF( eastern_side) THEN
827            DO jn = 1, jpts
828               tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn)
829               DO jk = 1, jpkm1
830                  DO jj = jmin,jmax
831                     IF( umask(nlci-2,jj,jk) == 0.e0 ) THEN
832                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
833                     ELSE
834                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
835                        IF( un(nlci-2,jj,jk) > 0.e0 ) THEN
836                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
837                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
838                        ENDIF
839                     ENDIF
840                  END DO
841               END DO
842            ENDDO
843         ENDIF
844         !
845         IF( northern_side ) THEN           
846            DO jn = 1, jpts
847               tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn)
848               DO jk = 1, jpkm1
849                  DO ji = imin,imax
850                     IF( vmask(ji,nlcj-2,jk) == 0.e0 ) THEN
851                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
852                     ELSE
853                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)       
854                        IF (vn(ji,nlcj-2,jk) > 0.e0 ) THEN
855                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
856                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
857                        ENDIF
858                     ENDIF
859                  END DO
860               END DO
861            ENDDO
862         ENDIF
863         !
864         IF( western_side) THEN           
865            DO jn = 1, jpts
866               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn)
867               DO jk = 1, jpkm1
868                  DO jj = jmin,jmax
869                     IF( umask(2,jj,jk) == 0.e0 ) THEN
870                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
871                     ELSE
872                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)       
873                        IF( un(2,jj,jk) < 0.e0 ) THEN
874                           tsa(2,jj,jk,jn)=(zalpha6*tsa(3,jj,jk,jn)+zalpha5*tsa(1,jj,jk,jn)+zalpha7*tsa(4,jj,jk,jn))*tmask(2,jj,jk)
875                        ENDIF
876                     ENDIF
877                  END DO
878               END DO
879            END DO
880         ENDIF
881         !
882         IF( southern_side ) THEN           
883            DO jn = 1, jpts
884               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn)
885               DO jk=1,jpk     
886                  DO ji=imin,imax
887                     IF( vmask(ji,2,jk) == 0.e0 ) THEN
888                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
889                     ELSE
890                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
891                        IF( vn(ji,2,jk) < 0.e0 ) THEN
892                           tsa(ji,2,jk,jn)=(zalpha6*tsa(ji,3,jk,jn)+zalpha5*tsa(ji,1,jk,jn)+zalpha7*tsa(ji,4,jk,jn))*tmask(ji,2,jk)
893                        ENDIF
894                     ENDIF
895                  END DO
896               END DO
897            ENDDO
898         ENDIF
899         !
900         ! Treatment of corners
901         !
902         ! East south
903         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
904            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:)
905         ENDIF
906         ! East north
907         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
908            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:)
909         ENDIF
910         ! West south
911         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
912            tsa(2,2,:,:) = ptab_child(2,2,:,:)
913         ENDIF
914         ! West north
915         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
916            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:)
917         ENDIF
918         !
919      ENDIF
920      !
921   END SUBROUTINE interptsn
922
923   SUBROUTINE interpsshn(ptab,i1,i2,j1,j2,before,nb,ndir)
924      !!----------------------------------------------------------------------
925      !!                  ***  ROUTINE interpsshn  ***
926      !!---------------------------------------------------------------------- 
927      INTEGER, INTENT(in) :: i1,i2,j1,j2
928      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
929      LOGICAL, INTENT(in) :: before
930      INTEGER, INTENT(in) :: nb , ndir
931      LOGICAL :: western_side, eastern_side,northern_side,southern_side
932      !!---------------------------------------------------------------------- 
933      !
934      IF( before) THEN
935         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
936      ELSE
937         western_side  = (nb == 1).AND.(ndir == 1)
938         eastern_side  = (nb == 1).AND.(ndir == 2)
939         southern_side = (nb == 2).AND.(ndir == 1)
940         northern_side = (nb == 2).AND.(ndir == 2)
941         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
942         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
943         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
944         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
945      ENDIF
946      !
947   END SUBROUTINE interpsshn
948
949   SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir)
950      !!---------------------------------------------
951      !!   *** ROUTINE interpun ***
952      !!---------------------------------------------   
953      !!
954      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
955      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
956      LOGICAL, INTENT(in) :: before
957      INTEGER, INTENT(in) :: nb , ndir
958      !!
959      INTEGER :: ji,jj,jk
960      REAL(wp) :: zrhoy
961! VERTICAL REFINEMENT BEGIN
962      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
963      REAL(wp), DIMENSION(k1:k2) :: tabin
964      REAL(wp) :: h_in(k1:k2)
965      REAL(wp) :: h_out(1:jpk)
966      INTEGER :: N_in, N_out
967      REAL(wp) :: h_diff
968      LOGICAL :: western_side, eastern_side
969      INTEGER :: iref
970
971! VERTICAL REFINEMENT END
972      !!---------------------------------------------   
973      !
974      IF (before) THEN
975         DO jk=1,jpk
976            DO jj=j1,j2
977               DO ji=i1,i2
978                  ptab(ji,jj,jk,1) = e2u(ji,jj) * un(ji,jj,jk)
979                  ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * fse3u(ji,jj,jk)
980                  ptab(ji,jj,jk,2) = fse3u(ji,jj,jk)
981               END DO
982            END DO
983         END DO
984      ELSE
985! VERTICAL REFINEMENT BEGIN
986         western_side  = (nb == 1).AND.(ndir == 1)
987         eastern_side  = (nb == 1).AND.(ndir == 2)
988         
989         ptab_child(:,:,:) = 0.
990         do jj=j1,j2
991         do ji=i1,i2
992         iref = ji
993         IF (western_side) iref = 2
994         IF (eastern_side) iref = nlci-2
995
996         N_in = 0
997         do jk=k1,k2
998           if (interp_scales_u(ji,jj,jk) == 0) EXIT
999             N_in = N_in + 1
1000             tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1001             h_in(N_in) = ptab(ji,jj,jk,2)
1002         enddo
1003         
1004         IF (N_in == 0) THEN
1005           ptab_child(ji,jj,:) = 0.
1006           CYCLE
1007         ENDIF
1008         
1009         N_out = 0
1010         do jk=1,jpk
1011           if (umask(iref,jj,jk) == 0) EXIT
1012           N_out = N_out + 1
1013           h_out(N_out) = fse3u(ji,jj,jk)
1014         enddo
1015         
1016         IF (N_out == 0) THEN
1017           ptab_child(ji,jj,:) = 0.
1018           CYCLE
1019         ENDIF
1020         
1021         h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
1022         IF (h_diff > 0.) THEN
1023           N_in = N_in + 1
1024           h_in(N_in) = h_diff
1025           tabin(N_in) = 0.
1026         ELSE
1027           h_out(N_out) = h_out(N_out) - h_diff
1028         ENDIF
1029         
1030         call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
1031         
1032         ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / fse3u(ji,jj,N_out)
1033
1034         ENDDO
1035         ENDDO
1036
1037! in the following
1038! remove division of ua by fs e3u (already done)
1039! VERTICAL REFINEMENT END
1040
1041         zrhoy = Agrif_Rhoy()
1042         DO jk=1,jpkm1
1043            DO jj=j1,j2
1044               ua(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj)))
1045            END DO
1046         END DO
1047      ENDIF
1048      !
1049   END SUBROUTINE interpun
1050
1051
1052   SUBROUTINE interpun2d(ptab,i1,i2,j1,j2,before)
1053      !!---------------------------------------------
1054      !!   *** ROUTINE interpun ***
1055      !!---------------------------------------------   
1056      !
1057      INTEGER, INTENT(in) :: i1,i2,j1,j2
1058      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1059      LOGICAL, INTENT(in) :: before
1060      !
1061      INTEGER :: ji,jj
1062      REAL(wp) :: ztref
1063      REAL(wp) :: zrhoy 
1064      !!---------------------------------------------   
1065      !
1066      ztref = 1.
1067
1068      IF (before) THEN
1069         DO jj=j1,j2
1070            DO ji=i1,MIN(i2,nlci-1)
1071               ptab(ji,jj) = e2u(ji,jj) * ((gcx(ji+1,jj) - gcx(ji,jj))/e1u(ji,jj)) 
1072            END DO
1073         END DO
1074      ELSE
1075         zrhoy = Agrif_Rhoy()
1076         DO jj=j1,j2
1077            laplacu(i1:i2,jj) = ztref * (ptab(i1:i2,jj)/(zrhoy*e2u(i1:i2,jj))) !*umask(i1:i2,jj,1)
1078         END DO
1079      ENDIF
1080      !
1081   END SUBROUTINE interpun2d
1082
1083
1084   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir)
1085      !!---------------------------------------------
1086      !!   *** ROUTINE interpvn ***
1087      !!---------------------------------------------   
1088      !
1089      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
1090      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1091      LOGICAL, INTENT(in) :: before
1092      INTEGER, INTENT(in) :: nb , ndir
1093      !
1094      INTEGER :: ji,jj,jk
1095      REAL(wp) :: zrhox
1096! VERTICAL REFINEMENT BEGIN
1097      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
1098      REAL(wp), DIMENSION(k1:k2) :: tabin
1099      REAL(wp) :: h_in(k1:k2)
1100      REAL(wp) :: h_out(1:jpk)
1101      INTEGER :: N_in, N_out
1102      REAL(wp) :: h_diff
1103      LOGICAL :: northern_side,southern_side
1104      INTEGER :: jref
1105
1106! VERTICAL REFINEMENT END
1107      !!---------------------------------------------   
1108      !     
1109      IF (before) THEN         
1110         !interpv entre 1 et k2 et interpv2d en jpkp1
1111         DO jk=k1,jpk
1112            DO jj=j1,j2
1113               DO ji=i1,i2
1114                  ptab(ji,jj,jk,1) = e1v(ji,jj) * vn(ji,jj,jk)
1115                  ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * fse3v(ji,jj,jk)
1116                  ptab(ji,jj,jk,2) = fse3v(ji,jj,jk)
1117               END DO
1118            END DO
1119         END DO
1120      ELSE       
1121! VERTICAL REFINEMENT BEGIN
1122         ptab_child(:,:,:) = 0.
1123         southern_side = (nb == 2).AND.(ndir == 1)
1124         northern_side = (nb == 2).AND.(ndir == 2)
1125         do jj=j1,j2
1126         jref = jj
1127         IF (southern_side) jref = 2
1128         IF (northern_side) jref = nlcj-2
1129         do ji=i1,i2
1130
1131         N_in = 0
1132         do jk=k1,k2
1133           if (interp_scales_v(ji,jj,jk) == 0) EXIT
1134             N_in = N_in + 1
1135             tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1136             h_in(N_in) = ptab(ji,jj,jk,2)
1137         enddo
1138         IF (N_in == 0) THEN
1139           ptab_child(ji,jj,:) = 0.
1140           CYCLE
1141         ENDIF
1142         
1143         N_out = 0
1144         do jk=1,jpk
1145           if (vmask(ji,jref,jk) == 0) EXIT
1146           N_out = N_out + 1
1147           h_out(N_out) = fse3v(ji,jj,jk)
1148         enddo
1149         IF (N_out == 0) THEN
1150           ptab_child(ji,jj,:) = 0.
1151           CYCLE
1152         ENDIF
1153         
1154         h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
1155         IF (h_diff > 0.) THEN
1156           N_in = N_in + 1
1157           h_in(N_in) = h_diff
1158           tabin(N_in) = 0.
1159         ELSE
1160           h_out(N_out) = h_out(N_out) - h_diff
1161         ENDIF
1162         
1163         call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
1164         
1165         ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / fse3v(ji,jj,N_out)
1166
1167         enddo
1168         enddo
1169! in the following
1170! remove division of va by fs e3v (already done)
1171! VERTICAL REFINEMENT END
1172         zrhox= Agrif_Rhox()
1173         DO jk=1,jpkm1
1174            DO jj=j1,j2
1175               va(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj)))
1176            END DO
1177         END DO
1178      ENDIF
1179      !       
1180   END SUBROUTINE interpvn
1181
1182   SUBROUTINE interpvn2d(ptab,i1,i2,j1,j2,before)
1183      !!---------------------------------------------
1184      !!   *** ROUTINE interpvn ***
1185      !!---------------------------------------------   
1186      !
1187      INTEGER, INTENT(in) :: i1,i2,j1,j2
1188      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1189      LOGICAL, INTENT(in) :: before
1190      !
1191      INTEGER :: ji,jj
1192      REAL(wp) :: zrhox 
1193      REAL(wp) :: ztref
1194      !!---------------------------------------------   
1195      !
1196      ztref = 1.   
1197      IF (before) THEN 
1198         !interpv entre 1 et k2 et interpv2d en jpkp1
1199         DO jj=j1,MIN(j2,nlcj-1)
1200            DO ji=i1,i2
1201               ptab(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) * vmask(ji,jj,1)
1202            END DO
1203         END DO
1204      ELSE           
1205         zrhox = Agrif_Rhox()
1206         DO ji=i1,i2
1207            laplacv(ji,j1:j2) = ztref * (ptab(ji,j1:j2)/(zrhox*e1v(ji,j1:j2)))
1208         END DO
1209      ENDIF
1210      !     
1211   END SUBROUTINE interpvn2d
1212
1213   SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir)
1214      !!----------------------------------------------------------------------
1215      !!                  ***  ROUTINE interpunb  ***
1216      !!---------------------------------------------------------------------- 
1217      INTEGER, INTENT(in) :: i1,i2,j1,j2
1218      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1219      LOGICAL, INTENT(in) :: before
1220      INTEGER, INTENT(in) :: nb , ndir
1221      !!
1222      INTEGER :: ji,jj
1223      REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
1224      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1225      !!---------------------------------------------------------------------- 
1226      !
1227      IF (before) THEN
1228         DO jj=j1,j2
1229            DO ji=i1,i2
1230               ptab(ji,jj) = un_b(ji,jj) * e2u(ji,jj) * hu(ji,jj) 
1231            END DO
1232         END DO
1233      ELSE
1234         western_side  = (nb == 1).AND.(ndir == 1)
1235         eastern_side  = (nb == 1).AND.(ndir == 2)
1236         southern_side = (nb == 2).AND.(ndir == 1)
1237         northern_side = (nb == 2).AND.(ndir == 2)
1238         zrhoy = Agrif_Rhoy()
1239         zrhot = Agrif_rhot()
1240         ! Time indexes bounds for integration
1241         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1242         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1243         ! Polynomial interpolation coefficients:
1244         IF( bdy_tinterp == 1 ) THEN
1245            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1246                  &      - zt0**2._wp * (       zt0 - 1._wp)        )
1247         ELSEIF( bdy_tinterp == 2 ) THEN
1248            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1249                  &      - zt0        * (       zt0 - 1._wp)**2._wp ) 
1250
1251         ELSE
1252            ztcoeff = 1
1253         ENDIF
1254         !   
1255         IF(western_side) THEN
1256            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1257         ENDIF
1258         IF(eastern_side) THEN
1259            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1260         ENDIF
1261         IF(southern_side) THEN
1262            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1263         ENDIF
1264         IF(northern_side) THEN
1265            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1266         ENDIF
1267         !           
1268         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1269            IF(western_side) THEN
1270               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
1271                     &                                  * umask(i1,j1:j2,1)
1272            ENDIF
1273            IF(eastern_side) THEN
1274               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
1275                     &                                  * umask(i1,j1:j2,1)
1276            ENDIF
1277            IF(southern_side) THEN
1278               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
1279                     &                                  * umask(i1:i2,j1,1)
1280            ENDIF
1281            IF(northern_side) THEN
1282               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
1283                     &                                  * umask(i1:i2,j1,1)
1284            ENDIF
1285         ENDIF
1286      ENDIF
1287      !
1288   END SUBROUTINE interpunb
1289
1290   SUBROUTINE interpvnb(ptab,i1,i2,j1,j2,before,nb,ndir)
1291      !!----------------------------------------------------------------------
1292      !!                  ***  ROUTINE interpvnb  ***
1293      !!---------------------------------------------------------------------- 
1294      INTEGER, INTENT(in) :: i1,i2,j1,j2
1295      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1296      LOGICAL, INTENT(in) :: before
1297      INTEGER, INTENT(in) :: nb , ndir
1298      !!
1299      INTEGER :: ji,jj
1300      REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff   
1301      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1302      !!---------------------------------------------------------------------- 
1303      !
1304      IF (before) THEN
1305         DO jj=j1,j2
1306            DO ji=i1,i2
1307               ptab(ji,jj) = vn_b(ji,jj) * e1v(ji,jj) * hv(ji,jj) 
1308            END DO
1309         END DO
1310      ELSE
1311         western_side  = (nb == 1).AND.(ndir == 1)
1312         eastern_side  = (nb == 1).AND.(ndir == 2)
1313         southern_side = (nb == 2).AND.(ndir == 1)
1314         northern_side = (nb == 2).AND.(ndir == 2)
1315         zrhox = Agrif_Rhox()
1316         zrhot = Agrif_rhot()
1317         ! Time indexes bounds for integration
1318         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1319         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1320         IF( bdy_tinterp == 1 ) THEN
1321            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1322                  &      - zt0**2._wp * (       zt0 - 1._wp)        )
1323         ELSEIF( bdy_tinterp == 2 ) THEN
1324            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1325                  &      - zt0        * (       zt0 - 1._wp)**2._wp ) 
1326
1327         ELSE
1328            ztcoeff = 1
1329         ENDIF
1330         !
1331         IF(western_side) THEN
1332            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1333         ENDIF
1334         IF(eastern_side) THEN
1335            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1336         ENDIF
1337         IF(southern_side) THEN
1338            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1339         ENDIF
1340         IF(northern_side) THEN
1341            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1342         ENDIF
1343         !           
1344         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1345            IF(western_side) THEN
1346               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1347                     &                                  * vmask(i1,j1:j2,1)
1348            ENDIF
1349            IF(eastern_side) THEN
1350               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1351                     &                                  * vmask(i1,j1:j2,1)
1352            ENDIF
1353            IF(southern_side) THEN
1354               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1355                     &                                  * vmask(i1:i2,j1,1)
1356            ENDIF
1357            IF(northern_side) THEN
1358               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1359                     &                                  * vmask(i1:i2,j1,1)
1360            ENDIF
1361         ENDIF
1362      ENDIF
1363      !
1364   END SUBROUTINE interpvnb
1365
1366   SUBROUTINE interpub2b(ptab,i1,i2,j1,j2,before,nb,ndir)
1367      !!----------------------------------------------------------------------
1368      !!                  ***  ROUTINE interpub2b  ***
1369      !!---------------------------------------------------------------------- 
1370      INTEGER, INTENT(in) :: i1,i2,j1,j2
1371      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1372      LOGICAL, INTENT(in) :: before
1373      INTEGER, INTENT(in) :: nb , ndir
1374      !!
1375      INTEGER :: ji,jj
1376      REAL(wp) :: zrhot, zt0, zt1,zat
1377      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1378      !!---------------------------------------------------------------------- 
1379      IF( before ) THEN
1380         DO jj=j1,j2
1381            DO ji=i1,i2
1382               ptab(ji,jj) = ub2_b(ji,jj) * e2u(ji,jj)
1383            END DO
1384         END DO
1385      ELSE
1386         western_side  = (nb == 1).AND.(ndir == 1)
1387         eastern_side  = (nb == 1).AND.(ndir == 2)
1388         southern_side = (nb == 2).AND.(ndir == 1)
1389         northern_side = (nb == 2).AND.(ndir == 2)
1390         zrhot = Agrif_rhot()
1391         ! Time indexes bounds for integration
1392         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1393         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1394         ! Polynomial interpolation coefficients:
1395         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
1396               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
1397         !
1398         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1399         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1400         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1401         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1402      ENDIF
1403      !
1404   END SUBROUTINE interpub2b
1405
1406   SUBROUTINE interpvb2b(ptab,i1,i2,j1,j2,before,nb,ndir)
1407      !!----------------------------------------------------------------------
1408      !!                  ***  ROUTINE interpvb2b  ***
1409      !!---------------------------------------------------------------------- 
1410      INTEGER, INTENT(in) :: i1,i2,j1,j2
1411      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1412      LOGICAL, INTENT(in) :: before
1413      INTEGER, INTENT(in) :: nb , ndir
1414      !!
1415      INTEGER :: ji,jj
1416      REAL(wp) :: zrhot, zt0, zt1,zat
1417      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1418      !!---------------------------------------------------------------------- 
1419      !
1420      IF( before ) THEN
1421         DO jj=j1,j2
1422            DO ji=i1,i2
1423               ptab(ji,jj) = vb2_b(ji,jj) * e1v(ji,jj)
1424            END DO
1425         END DO
1426      ELSE     
1427         western_side  = (nb == 1).AND.(ndir == 1)
1428         eastern_side  = (nb == 1).AND.(ndir == 2)
1429         southern_side = (nb == 2).AND.(ndir == 1)
1430         northern_side = (nb == 2).AND.(ndir == 2)
1431         zrhot = Agrif_rhot()
1432         ! Time indexes bounds for integration
1433         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1434         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1435         ! Polynomial interpolation coefficients:
1436         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
1437               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
1438         !
1439         IF(western_side ) vbdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1440         IF(eastern_side ) vbdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1441         IF(southern_side) vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1442         IF(northern_side) vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1443      ENDIF
1444      !     
1445   END SUBROUTINE interpvb2b
1446
1447   SUBROUTINE interpe3t(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1448      !!----------------------------------------------------------------------
1449      !!                  ***  ROUTINE interpe3t  ***
1450      !!---------------------------------------------------------------------- 
1451      !
1452      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1453      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1454      LOGICAL :: before
1455      INTEGER, INTENT(in) :: nb , ndir
1456      !
1457      INTEGER :: ji, jj, jk
1458      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1459      REAL(wp) :: ztmpmsk     
1460      !!---------------------------------------------------------------------- 
1461      !   
1462      IF (before) THEN
1463         DO jk=k1,k2
1464            DO jj=j1,j2
1465               DO ji=i1,i2
1466                  ptab(ji,jj,jk) = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
1467               END DO
1468            END DO
1469         END DO
1470      ELSE
1471         western_side  = (nb == 1).AND.(ndir == 1)
1472         eastern_side  = (nb == 1).AND.(ndir == 2)
1473         southern_side = (nb == 2).AND.(ndir == 1)
1474         northern_side = (nb == 2).AND.(ndir == 2)
1475
1476         DO jk=k1,k2
1477            DO jj=j1,j2
1478               DO ji=i1,i2
1479                  ! Get velocity mask at boundary edge points:
1480                  IF (western_side)  ztmpmsk = umask(ji    ,jj    ,1)
1481                  IF (eastern_side)  ztmpmsk = umask(nlci-2,jj    ,1)
1482                  IF (northern_side) ztmpmsk = vmask(ji    ,nlcj-2,1)
1483                  IF (southern_side) ztmpmsk = vmask(ji    ,2     ,1)
1484
1485                  IF (ABS(ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk))*ztmpmsk > 1.D-2) THEN
1486                     IF (western_side) THEN
1487                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1488                     ELSEIF (eastern_side) THEN
1489                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1490                     ELSEIF (southern_side) THEN
1491                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1492                     ELSEIF (northern_side) THEN
1493                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1494                     ENDIF
1495                     WRITE(numout,*) '      ptab(ji,jj,jk), fse3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1496                     kindic_agr = kindic_agr + 1
1497                  ENDIF
1498               END DO
1499            END DO
1500         END DO
1501
1502      ENDIF
1503      !
1504   END SUBROUTINE interpe3t
1505
1506   SUBROUTINE interpumsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1507      !!----------------------------------------------------------------------
1508      !!                  ***  ROUTINE interpumsk  ***
1509      !!---------------------------------------------------------------------- 
1510      !
1511      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1512      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1513      LOGICAL :: before
1514      INTEGER, INTENT(in) :: nb , ndir
1515      !
1516      INTEGER :: ji, jj, jk
1517      LOGICAL :: western_side, eastern_side   
1518      !!---------------------------------------------------------------------- 
1519      !   
1520      IF (before) THEN
1521         DO jk=k1,k2
1522            DO jj=j1,j2
1523               DO ji=i1,i2
1524                  ptab(ji,jj,jk) = umask(ji,jj,jk)
1525               END DO
1526            END DO
1527         END DO
1528      ELSE
1529
1530         western_side  = (nb == 1).AND.(ndir == 1)
1531         eastern_side  = (nb == 1).AND.(ndir == 2)
1532         DO jk=k1,k2
1533            DO jj=j1,j2
1534               DO ji=i1,i2
1535                   ! Velocity mask at boundary edge points:
1536                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1537                     IF (western_side) THEN
1538                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1539                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1540                        kindic_agr = kindic_agr + 1
1541                     ELSEIF (eastern_side) THEN
1542                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1543                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1544                        kindic_agr = kindic_agr + 1
1545                     ENDIF
1546                  ENDIF
1547               END DO
1548            END DO
1549         END DO
1550
1551      ENDIF
1552      !
1553   END SUBROUTINE interpumsk
1554
1555   SUBROUTINE interpvmsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1556      !!----------------------------------------------------------------------
1557      !!                  ***  ROUTINE interpvmsk  ***
1558      !!---------------------------------------------------------------------- 
1559      !
1560      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1561      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1562      LOGICAL :: before
1563      INTEGER, INTENT(in) :: nb , ndir
1564      !
1565      INTEGER :: ji, jj, jk
1566      LOGICAL :: northern_side, southern_side     
1567      !!---------------------------------------------------------------------- 
1568      !   
1569      IF (before) THEN
1570         DO jk=k1,k2
1571            DO jj=j1,j2
1572               DO ji=i1,i2
1573                  ptab(ji,jj,jk) = vmask(ji,jj,jk)
1574               END DO
1575            END DO
1576         END DO
1577      ELSE
1578
1579         southern_side = (nb == 2).AND.(ndir == 1)
1580         northern_side = (nb == 2).AND.(ndir == 2)
1581         DO jk=k1,k2
1582            DO jj=j1,j2
1583               DO ji=i1,i2
1584                   ! Velocity mask at boundary edge points:
1585                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1586                     IF (southern_side) THEN
1587                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1588                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1589                        kindic_agr = kindic_agr + 1
1590                     ELSEIF (northern_side) THEN
1591                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1592                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1593                        kindic_agr = kindic_agr + 1
1594                     ENDIF
1595                  ENDIF
1596               END DO
1597            END DO
1598         END DO
1599
1600      ENDIF
1601      !
1602   END SUBROUTINE interpvmsk
1603
1604# if defined key_zdftke
1605
1606   SUBROUTINE interpavm(ptab,i1,i2,j1,j2,k1,k2,before)
1607      !!----------------------------------------------------------------------
1608      !!                  ***  ROUTINE interavm  ***
1609      !!---------------------------------------------------------------------- 
1610      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1611      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1612      LOGICAL, INTENT(in) :: before
1613      !!---------------------------------------------------------------------- 
1614      !     
1615      IF( before) THEN
1616         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1617      ELSE
1618         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1619      ENDIF
1620      !
1621   END SUBROUTINE interpavm
1622
1623# endif /* key_zdftke */
1624
1625#else
1626   !!----------------------------------------------------------------------
1627   !!   Empty module                                          no AGRIF zoom
1628   !!----------------------------------------------------------------------
1629CONTAINS
1630   SUBROUTINE Agrif_OPA_Interp_empty
1631      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1632   END SUBROUTINE Agrif_OPA_Interp_empty
1633#endif
1634
1635   !!======================================================================
1636END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.