source: branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 7988

Last change on this file since 7988 was 7988, checked in by jchanut, 3 years ago

Add AGRIF proper AGRIF bcs to GLS and TKE + vvl update

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