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/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 4984

Last change on this file since 4984 was 4984, checked in by jchanut, 9 years ago

AGRIF: Improve bathymetry checks at child boundaries

  • Property svn:keywords set to Id
File size: 50.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   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
45   PUBLIC   Agrif_tke, interpavt, interpavm, interpavmu, interpavmv
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(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)/hu(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(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)/hu(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(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)/hu(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(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)/hu(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(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)/hv(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(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)/hv(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(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)/hv(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(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)/hv(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      CALL wrk_dealloc( jpi, jpj, spgv1, spgu1 )
419      !
420   END SUBROUTINE Agrif_dyn
421
422   SUBROUTINE Agrif_dyn_ts( jn )
423      !!----------------------------------------------------------------------
424      !!                  ***  ROUTINE Agrif_dyn_ts  ***
425      !!---------------------------------------------------------------------- 
426      !!
427      INTEGER, INTENT(in) ::   jn
428      !!
429      INTEGER :: ji, jj
430      !!---------------------------------------------------------------------- 
431
432      IF( Agrif_Root() )   RETURN
433
434      IF((nbondi == -1).OR.(nbondi == 2)) THEN
435         DO jj=1,jpj
436            va_e(2,jj) = vbdy_w(jj) * hvr_e(2,jj)
437            ! Specified fluxes:
438            ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj)
439            ! Characteristics method:
440            !alt            ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
441            !alt                       &           - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
442         END DO
443      ENDIF
444
445      IF((nbondi == 1).OR.(nbondi == 2)) THEN
446         DO jj=1,jpj
447            va_e(nlci-1,jj) = vbdy_e(jj) * hvr_e(nlci-1,jj)
448            ! Specified fluxes:
449            ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj)
450            ! Characteristics method:
451            !alt            ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
452            !alt                            &           + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
453         END DO
454      ENDIF
455
456      IF((nbondj == -1).OR.(nbondj == 2)) THEN
457         DO ji=1,jpi
458            ua_e(ji,2) = ubdy_s(ji) * hur_e(ji,2)
459            ! Specified fluxes:
460            va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2)
461            ! Characteristics method:
462            !alt            va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
463            !alt                       &           - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
464         END DO
465      ENDIF
466
467      IF((nbondj == 1).OR.(nbondj == 2)) THEN
468         DO ji=1,jpi
469            ua_e(ji,nlcj-1) = ubdy_n(ji) * hur_e(ji,nlcj-1)
470            ! Specified fluxes:
471            va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2)
472            ! Characteristics method:
473            !alt            va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2)  + va_e(ji,nlcj-3) &
474            !alt                            &           + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
475         END DO
476      ENDIF
477      !
478   END SUBROUTINE Agrif_dyn_ts
479
480   SUBROUTINE Agrif_dta_ts( kt )
481      !!----------------------------------------------------------------------
482      !!                  ***  ROUTINE Agrif_dta_ts  ***
483      !!---------------------------------------------------------------------- 
484      !!
485      INTEGER, INTENT(in) ::   kt
486      !!
487      INTEGER :: ji, jj
488      LOGICAL :: ll_int_cons
489      REAL(wp) :: zrhot, zt
490      !!---------------------------------------------------------------------- 
491
492      IF( Agrif_Root() )   RETURN
493
494      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in
495      ! the forward case only
496
497      zrhot = Agrif_rhot()
498
499      ! "Central" time index for interpolation:
500      IF (ln_bt_fw) THEN
501         zt = REAL(Agrif_NbStepint()+0.5_wp,wp) / zrhot
502      ELSE
503         zt = REAL(Agrif_NbStepint(),wp) / zrhot
504      ENDIF
505
506      ! Linear interpolation of sea level
507      Agrif_SpecialValue    = 0.e0
508      Agrif_UseSpecialValue = .TRUE.
509      CALL Agrif_Bc_variable(sshn_id,calledweight=zt, procname=interpsshn )
510      Agrif_UseSpecialValue = .FALSE.
511
512      ! Interpolate barotropic fluxes
513      Agrif_SpecialValue=0.
514      Agrif_UseSpecialValue = ln_spc_dyn
515
516      IF (ll_int_cons) THEN ! Conservative interpolation
517         ! orders matters here !!!!!!
518         CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1._wp, procname=interpub2b) ! Time integrated
519         CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1._wp, procname=interpvb2b)
520         bdy_tinterp = 1
521         CALL Agrif_Bc_variable(unb_id ,calledweight=1._wp, procname=interpunb) ! After
522         CALL Agrif_Bc_variable(vnb_id ,calledweight=1._wp, procname=interpvnb)
523         bdy_tinterp = 2
524         CALL Agrif_Bc_variable(unb_id ,calledweight=0._wp, procname=interpunb) ! Before
525         CALL Agrif_Bc_variable(vnb_id ,calledweight=0._wp, procname=interpvnb)         
526      ELSE ! Linear interpolation
527         bdy_tinterp = 0
528         ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 
529         ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 
530         ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 
531         ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 
532         CALL Agrif_Bc_variable(unb_id,calledweight=zt, procname=interpunb)
533         CALL Agrif_Bc_variable(vnb_id,calledweight=zt, procname=interpvnb)
534      ENDIF
535      Agrif_UseSpecialValue = .FALSE.
536      !
537   END SUBROUTINE Agrif_dta_ts
538
539   SUBROUTINE Agrif_ssh( kt )
540      !!----------------------------------------------------------------------
541      !!                  ***  ROUTINE Agrif_DYN  ***
542      !!---------------------------------------------------------------------- 
543      INTEGER, INTENT(in) ::   kt
544      !!
545      !!---------------------------------------------------------------------- 
546
547      IF( Agrif_Root() )   RETURN
548
549      IF((nbondi == -1).OR.(nbondi == 2)) THEN
550         ssha(2,:)=ssha(3,:)
551         sshn(2,:)=sshn(3,:)
552      ENDIF
553
554      IF((nbondi == 1).OR.(nbondi == 2)) THEN
555         ssha(nlci-1,:)=ssha(nlci-2,:)
556         sshn(nlci-1,:)=sshn(nlci-2,:)
557      ENDIF
558
559      IF((nbondj == -1).OR.(nbondj == 2)) THEN
560         ssha(:,2)=ssha(:,3)
561         sshn(:,2)=sshn(:,3)
562      ENDIF
563
564      IF((nbondj == 1).OR.(nbondj == 2)) THEN
565         ssha(:,nlcj-1)=ssha(:,nlcj-2)
566         sshn(:,nlcj-1)=sshn(:,nlcj-2)
567      ENDIF
568
569   END SUBROUTINE Agrif_ssh
570
571   SUBROUTINE Agrif_ssh_ts( jn )
572      !!----------------------------------------------------------------------
573      !!                  ***  ROUTINE Agrif_ssh_ts  ***
574      !!---------------------------------------------------------------------- 
575      INTEGER, INTENT(in) ::   jn
576      !!
577      INTEGER :: ji,jj
578      !!---------------------------------------------------------------------- 
579
580      IF((nbondi == -1).OR.(nbondi == 2)) THEN
581         DO jj=1,jpj
582            ssha_e(2,jj) = hbdy_w(jj)
583         END DO
584      ENDIF
585
586      IF((nbondi == 1).OR.(nbondi == 2)) THEN
587         DO jj=1,jpj
588            ssha_e(nlci-1,jj) = hbdy_e(jj)
589         END DO
590      ENDIF
591
592      IF((nbondj == -1).OR.(nbondj == 2)) THEN
593         DO ji=1,jpi
594            ssha_e(ji,2) = hbdy_s(ji)
595         END DO
596      ENDIF
597
598      IF((nbondj == 1).OR.(nbondj == 2)) THEN
599         DO ji=1,jpi
600            ssha_e(ji,nlcj-1) = hbdy_n(ji)
601         END DO
602      ENDIF
603
604   END SUBROUTINE Agrif_ssh_ts
605
606# if defined key_zdftke
607   SUBROUTINE Agrif_tke
608      !!----------------------------------------------------------------------
609      !!                  ***  ROUTINE Agrif_tke  ***
610      !!---------------------------------------------------------------------- 
611      !
612      IF( Agrif_Root() )   RETURN
613
614     
615      Agrif_SpecialValue    = 0.e0
616      Agrif_UseSpecialValue = .TRUE.
617     
618      CALL Agrif_Bc_variable(avt_id , procname=interpavt)       
619      CALL Agrif_Bc_variable(avm_id , procname=interpavm)       
620      CALL Agrif_Bc_variable(avmu_id, procname=interpavmu)
621      CALL Agrif_Bc_variable(avmv_id, procname=interpavmv)
622             
623      Agrif_UseSpecialValue = .FALSE.
624      !
625   END SUBROUTINE Agrif_tke
626# endif
627
628   SUBROUTINE interptsn(ptab,i1,i2,j1,j2,k1,k2,n1,n2,before,nb,ndir)
629      !!---------------------------------------------
630      !!   *** ROUTINE interptsn ***
631      !!---------------------------------------------
632      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
633      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
634      LOGICAL, INTENT(in) :: before
635      INTEGER, INTENT(in) :: nb , ndir
636      !
637      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
638      INTEGER :: imin, imax, jmin, jmax
639      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
640      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
641      LOGICAL :: western_side, eastern_side,northern_side,southern_side
642
643      IF (before) THEN         
644         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
645      ELSE
646         !
647         western_side  = (nb == 1).AND.(ndir == 1)
648         eastern_side  = (nb == 1).AND.(ndir == 2)
649         southern_side = (nb == 2).AND.(ndir == 1)
650         northern_side = (nb == 2).AND.(ndir == 2)
651         !
652         zrhox = Agrif_Rhox()
653         !
654         zalpha1 = ( zrhox - 1. ) * 0.5
655         zalpha2 = 1. - zalpha1
656         !
657         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
658         zalpha4 = 1. - zalpha3
659         !
660         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
661         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
662         zalpha5 = 1. - zalpha6 - zalpha7
663         !
664         imin = i1
665         imax = i2
666         jmin = j1
667         jmax = j2
668         !
669         ! Remove CORNERS
670         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
671         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
672         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
673         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2       
674         !
675         IF( eastern_side) THEN
676            DO jn = 1, jpts
677               tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
678               DO jk = 1, jpkm1
679                  DO jj = jmin,jmax
680                     IF( umask(nlci-2,jj,jk) == 0.e0 ) THEN
681                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
682                     ELSE
683                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
684                        IF( un(nlci-2,jj,jk) > 0.e0 ) THEN
685                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
686                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
687                        ENDIF
688                     ENDIF
689                  END DO
690               END DO
691            ENDDO
692         ENDIF
693         !
694         IF( northern_side ) THEN           
695            DO jn = 1, jpts
696               tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)
697               DO jk = 1, jpkm1
698                  DO ji = imin,imax
699                     IF( vmask(ji,nlcj-2,jk) == 0.e0 ) THEN
700                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
701                     ELSE
702                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)       
703                        IF (vn(ji,nlcj-2,jk) > 0.e0 ) THEN
704                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
705                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
706                        ENDIF
707                     ENDIF
708                  END DO
709               END DO
710            ENDDO
711         ENDIF
712         !
713         IF( western_side) THEN           
714            DO jn = 1, jpts
715               tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
716               DO jk = 1, jpkm1
717                  DO jj = jmin,jmax
718                     IF( umask(2,jj,jk) == 0.e0 ) THEN
719                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
720                     ELSE
721                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)       
722                        IF( un(2,jj,jk) < 0.e0 ) THEN
723                           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)
724                        ENDIF
725                     ENDIF
726                  END DO
727               END DO
728            END DO
729         ENDIF
730         !
731         IF( southern_side ) THEN           
732            DO jn = 1, jpts
733               tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
734               DO jk=1,jpk     
735                  DO ji=imin,imax
736                     IF( vmask(ji,2,jk) == 0.e0 ) THEN
737                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
738                     ELSE
739                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
740                        IF( vn(ji,2,jk) < 0.e0 ) THEN
741                           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)
742                        ENDIF
743                     ENDIF
744                  END DO
745               END DO
746            ENDDO
747         ENDIF
748         !
749         ! Treatment of corners
750         !
751         ! East south
752         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
753            tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)
754         ENDIF
755         ! East north
756         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
757            tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)
758         ENDIF
759         ! West south
760         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
761            tsa(2,2,:,:) = ptab(2,2,:,:)
762         ENDIF
763         ! West north
764         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
765            tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)
766         ENDIF
767         !
768      ENDIF
769      !
770   END SUBROUTINE interptsn
771
772   SUBROUTINE interpsshn(ptab,i1,i2,j1,j2,before,nb,ndir)
773      !!----------------------------------------------------------------------
774      !!                  ***  ROUTINE interpsshn  ***
775      !!---------------------------------------------------------------------- 
776      INTEGER, INTENT(in) :: i1,i2,j1,j2
777      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
778      LOGICAL, INTENT(in) :: before
779      INTEGER, INTENT(in) :: nb , ndir
780      LOGICAL :: western_side, eastern_side,northern_side,southern_side
781      !!---------------------------------------------------------------------- 
782      !
783      IF( before) THEN
784         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
785      ELSE
786         western_side  = (nb == 1).AND.(ndir == 1)
787         eastern_side  = (nb == 1).AND.(ndir == 2)
788         southern_side = (nb == 2).AND.(ndir == 1)
789         northern_side = (nb == 2).AND.(ndir == 2)
790         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
791         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
792         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
793         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
794      ENDIF
795      !
796   END SUBROUTINE interpsshn
797
798   SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2, before)
799      !!---------------------------------------------
800      !!   *** ROUTINE interpun ***
801      !!---------------------------------------------   
802      !!
803      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
804      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
805      LOGICAL, INTENT(in) :: before
806      !!
807      INTEGER :: ji,jj,jk
808      REAL(wp) :: zrhoy 
809      !!---------------------------------------------   
810      !
811      IF (before) THEN
812         DO jk=1,jpk
813            DO jj=j1,j2
814               DO ji=i1,i2
815                  ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
816                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u(ji,jj,jk)
817               END DO
818            END DO
819         END DO
820      ELSE
821         zrhoy = Agrif_Rhoy()
822         DO jk=1,jpkm1
823            DO jj=j1,j2
824               ua(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj)))
825               ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u(i1:i2,jj,jk)
826            END DO
827         END DO
828      ENDIF
829      !
830   END SUBROUTINE interpun
831
832
833   SUBROUTINE interpun2d(ptab,i1,i2,j1,j2,before)
834      !!---------------------------------------------
835      !!   *** ROUTINE interpun ***
836      !!---------------------------------------------   
837      !
838      INTEGER, INTENT(in) :: i1,i2,j1,j2
839      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
840      LOGICAL, INTENT(in) :: before
841      !
842      INTEGER :: ji,jj
843      REAL(wp) :: ztref
844      REAL(wp) :: zrhoy 
845      !!---------------------------------------------   
846      !
847      ztref = 1.
848
849      IF (before) THEN
850         DO jj=j1,j2
851            DO ji=i1,MIN(i2,nlci-1)
852               ptab(ji,jj) = e2u(ji,jj) * ((gcx(ji+1,jj) - gcx(ji,jj))/e1u(ji,jj)) 
853            END DO
854         END DO
855      ELSE
856         zrhoy = Agrif_Rhoy()
857         DO jj=j1,j2
858            laplacu(i1:i2,jj) = ztref * (ptab(i1:i2,jj)/(zrhoy*e2u(i1:i2,jj))) !*umask(i1:i2,jj,1)
859         END DO
860      ENDIF
861      !
862   END SUBROUTINE interpun2d
863
864
865   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before)
866      !!---------------------------------------------
867      !!   *** ROUTINE interpvn ***
868      !!---------------------------------------------   
869      !
870      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
871      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
872      LOGICAL, INTENT(in) :: before
873      !
874      INTEGER :: ji,jj,jk
875      REAL(wp) :: zrhox 
876      !!---------------------------------------------   
877      !     
878      IF (before) THEN         
879         !interpv entre 1 et k2 et interpv2d en jpkp1
880         DO jk=k1,jpk
881            DO jj=j1,j2
882               DO ji=i1,i2
883                  ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
884                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v(ji,jj,jk)
885               END DO
886            END DO
887         END DO
888      ELSE         
889         zrhox= Agrif_Rhox()
890         DO jk=1,jpkm1
891            DO jj=j1,j2
892               va(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj)))
893               va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v(i1:i2,jj,jk)
894            END DO
895         END DO
896      ENDIF
897      !       
898   END SUBROUTINE interpvn
899
900   SUBROUTINE interpvn2d(ptab,i1,i2,j1,j2,before)
901      !!---------------------------------------------
902      !!   *** ROUTINE interpvn ***
903      !!---------------------------------------------   
904      !
905      INTEGER, INTENT(in) :: i1,i2,j1,j2
906      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
907      LOGICAL, INTENT(in) :: before
908      !
909      INTEGER :: ji,jj
910      REAL(wp) :: zrhox 
911      REAL(wp) :: ztref
912      !!---------------------------------------------   
913      !
914      ztref = 1.   
915      IF (before) THEN 
916         !interpv entre 1 et k2 et interpv2d en jpkp1
917         DO jj=j1,MIN(j2,nlcj-1)
918            DO ji=i1,i2
919               ptab(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) * vmask(ji,jj,1)
920            END DO
921         END DO
922      ELSE           
923         zrhox = Agrif_Rhox()
924         DO ji=i1,i2
925            laplacv(ji,j1:j2) = ztref * (ptab(ji,j1:j2)/(zrhox*e1v(ji,j1:j2)))
926         END DO
927      ENDIF
928      !     
929   END SUBROUTINE interpvn2d
930
931   SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir)
932      !!----------------------------------------------------------------------
933      !!                  ***  ROUTINE interpunb  ***
934      !!---------------------------------------------------------------------- 
935      INTEGER, INTENT(in) :: i1,i2,j1,j2
936      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
937      LOGICAL, INTENT(in) :: before
938      INTEGER, INTENT(in) :: nb , ndir
939      !!
940      INTEGER :: ji,jj
941      REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
942      LOGICAL :: western_side, eastern_side,northern_side,southern_side
943      !!---------------------------------------------------------------------- 
944      !
945      IF (before) THEN
946         DO jj=j1,j2
947            DO ji=i1,i2
948               ptab(ji,jj) = un_b(ji,jj) * e2u(ji,jj) * hu(ji,jj) 
949            END DO
950         END DO
951      ELSE
952         western_side  = (nb == 1).AND.(ndir == 1)
953         eastern_side  = (nb == 1).AND.(ndir == 2)
954         southern_side = (nb == 2).AND.(ndir == 1)
955         northern_side = (nb == 2).AND.(ndir == 2)
956         zrhoy = Agrif_Rhoy()
957         zrhot = Agrif_rhot()
958         ! Time indexes bounds for integration
959         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
960         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
961         ! Polynomial interpolation coefficients:
962         IF( bdy_tinterp == 1 ) THEN
963            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
964                  &      - zt0**2._wp * (       zt0 - 1._wp)        )
965         ELSEIF( bdy_tinterp == 2 ) THEN
966            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
967                  &      - zt0        * (       zt0 - 1._wp)**2._wp ) 
968
969         ELSE
970            ztcoeff = 1
971         ENDIF
972         !   
973         IF(western_side) THEN
974            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
975         ENDIF
976         IF(eastern_side) THEN
977            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
978         ENDIF
979         IF(southern_side) THEN
980            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
981         ENDIF
982         IF(northern_side) THEN
983            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
984         ENDIF
985         !           
986         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
987            IF(western_side) THEN
988               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
989                     &                                  * umask(i1,j1:j2,1)
990            ENDIF
991            IF(eastern_side) THEN
992               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
993                     &                                  * umask(i1,j1:j2,1)
994            ENDIF
995            IF(southern_side) THEN
996               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
997                     &                                  * umask(i1:i2,j1,1)
998            ENDIF
999            IF(northern_side) THEN
1000               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
1001                     &                                  * umask(i1:i2,j1,1)
1002            ENDIF
1003         ENDIF
1004      ENDIF
1005      !
1006   END SUBROUTINE interpunb
1007
1008   SUBROUTINE interpvnb(ptab,i1,i2,j1,j2,before,nb,ndir)
1009      !!----------------------------------------------------------------------
1010      !!                  ***  ROUTINE interpvnb  ***
1011      !!---------------------------------------------------------------------- 
1012      INTEGER, INTENT(in) :: i1,i2,j1,j2
1013      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1014      LOGICAL, INTENT(in) :: before
1015      INTEGER, INTENT(in) :: nb , ndir
1016      !!
1017      INTEGER :: ji,jj
1018      REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff   
1019      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1020      !!---------------------------------------------------------------------- 
1021      !
1022      IF (before) THEN
1023         DO jj=j1,j2
1024            DO ji=i1,i2
1025               ptab(ji,jj) = vn_b(ji,jj) * e1v(ji,jj) * hv(ji,jj) 
1026            END DO
1027         END DO
1028      ELSE
1029         western_side  = (nb == 1).AND.(ndir == 1)
1030         eastern_side  = (nb == 1).AND.(ndir == 2)
1031         southern_side = (nb == 2).AND.(ndir == 1)
1032         northern_side = (nb == 2).AND.(ndir == 2)
1033         zrhox = Agrif_Rhox()
1034         zrhot = Agrif_rhot()
1035         ! Time indexes bounds for integration
1036         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1037         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1038         IF( bdy_tinterp == 1 ) THEN
1039            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1040                  &      - zt0**2._wp * (       zt0 - 1._wp)        )
1041         ELSEIF( bdy_tinterp == 2 ) THEN
1042            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1043                  &      - zt0        * (       zt0 - 1._wp)**2._wp ) 
1044
1045         ELSE
1046            ztcoeff = 1
1047         ENDIF
1048         !
1049         IF(western_side) THEN
1050            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1051         ENDIF
1052         IF(eastern_side) THEN
1053            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1054         ENDIF
1055         IF(southern_side) THEN
1056            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1057         ENDIF
1058         IF(northern_side) THEN
1059            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1060         ENDIF
1061         !           
1062         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1063            IF(western_side) THEN
1064               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1065                     &                                  * vmask(i1,j1:j2,1)
1066            ENDIF
1067            IF(eastern_side) THEN
1068               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1069                     &                                  * vmask(i1,j1:j2,1)
1070            ENDIF
1071            IF(southern_side) THEN
1072               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1073                     &                                  * vmask(i1:i2,j1,1)
1074            ENDIF
1075            IF(northern_side) THEN
1076               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1077                     &                                  * vmask(i1:i2,j1,1)
1078            ENDIF
1079         ENDIF
1080      ENDIF
1081      !
1082   END SUBROUTINE interpvnb
1083
1084   SUBROUTINE interpub2b(ptab,i1,i2,j1,j2,before,nb,ndir)
1085      !!----------------------------------------------------------------------
1086      !!                  ***  ROUTINE interpub2b  ***
1087      !!---------------------------------------------------------------------- 
1088      INTEGER, INTENT(in) :: i1,i2,j1,j2
1089      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1090      LOGICAL, INTENT(in) :: before
1091      INTEGER, INTENT(in) :: nb , ndir
1092      !!
1093      INTEGER :: ji,jj
1094      REAL(wp) :: zrhot, zt0, zt1,zat
1095      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1096      !!---------------------------------------------------------------------- 
1097      IF( before ) THEN
1098         DO jj=j1,j2
1099            DO ji=i1,i2
1100               ptab(ji,jj) = ub2_b(ji,jj) * e2u(ji,jj)
1101            END DO
1102         END DO
1103      ELSE
1104         western_side  = (nb == 1).AND.(ndir == 1)
1105         eastern_side  = (nb == 1).AND.(ndir == 2)
1106         southern_side = (nb == 2).AND.(ndir == 1)
1107         northern_side = (nb == 2).AND.(ndir == 2)
1108         zrhot = Agrif_rhot()
1109         ! Time indexes bounds for integration
1110         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1111         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1112         ! Polynomial interpolation coefficients:
1113         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
1114               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
1115         !
1116         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1117         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1118         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1119         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1120      ENDIF
1121      !
1122   END SUBROUTINE interpub2b
1123
1124   SUBROUTINE interpvb2b(ptab,i1,i2,j1,j2,before,nb,ndir)
1125      !!----------------------------------------------------------------------
1126      !!                  ***  ROUTINE interpvb2b  ***
1127      !!---------------------------------------------------------------------- 
1128      INTEGER, INTENT(in) :: i1,i2,j1,j2
1129      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1130      LOGICAL, INTENT(in) :: before
1131      INTEGER, INTENT(in) :: nb , ndir
1132      !!
1133      INTEGER :: ji,jj
1134      REAL(wp) :: zrhot, zt0, zt1,zat
1135      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1136      !!---------------------------------------------------------------------- 
1137      !
1138      IF( before ) THEN
1139         DO jj=j1,j2
1140            DO ji=i1,i2
1141               ptab(ji,jj) = vb2_b(ji,jj) * e1v(ji,jj)
1142            END DO
1143         END DO
1144      ELSE     
1145         western_side  = (nb == 1).AND.(ndir == 1)
1146         eastern_side  = (nb == 1).AND.(ndir == 2)
1147         southern_side = (nb == 2).AND.(ndir == 1)
1148         northern_side = (nb == 2).AND.(ndir == 2)
1149         zrhot = Agrif_rhot()
1150         ! Time indexes bounds for integration
1151         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1152         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1153         ! Polynomial interpolation coefficients:
1154         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
1155               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
1156         !
1157         IF(western_side ) vbdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1158         IF(eastern_side ) vbdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1159         IF(southern_side) vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1160         IF(northern_side) vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1161      ENDIF
1162      !     
1163   END SUBROUTINE interpvb2b
1164
1165   SUBROUTINE interpe3t(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1166      !!----------------------------------------------------------------------
1167      !!                  ***  ROUTINE interpe3t  ***
1168      !!---------------------------------------------------------------------- 
1169      !
1170      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1171      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1172      LOGICAL :: before
1173      INTEGER, INTENT(in) :: nb , ndir
1174      !
1175      INTEGER :: ji, jj, jk
1176      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1177      REAL(wp) :: ztmpmsk     
1178      !!---------------------------------------------------------------------- 
1179      !   
1180      IF (before) THEN
1181         DO jk=k1,k2
1182            DO jj=j1,j2
1183               DO ji=i1,i2
1184                  ptab(ji,jj,jk) = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
1185               END DO
1186            END DO
1187         END DO
1188      ELSE
1189         western_side  = (nb == 1).AND.(ndir == 1)
1190         eastern_side  = (nb == 1).AND.(ndir == 2)
1191         southern_side = (nb == 2).AND.(ndir == 1)
1192         northern_side = (nb == 2).AND.(ndir == 2)
1193
1194         DO jk=k1,k2
1195            DO jj=j1,j2
1196               DO ji=i1,i2
1197                  ! Get velocity mask at boundary edge points:
1198                  IF (western_side)  ztmpmsk = umask(ji    ,jj    ,1)
1199                  IF (eastern_side)  ztmpmsk = umask(nlci-2,jj    ,1)
1200                  IF (northern_side) ztmpmsk = vmask(ji    ,nlcj-2,1)
1201                  IF (southern_side) ztmpmsk = vmask(ji    ,2     ,1)
1202
1203                  IF (ABS(ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk))*ztmpmsk > 1.D-2) THEN
1204                     IF (western_side) THEN
1205                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1206                     ELSEIF (eastern_side) THEN
1207                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1208                     ELSEIF (southern_side) THEN
1209                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1210                     ELSEIF (northern_side) THEN
1211                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1212                     ENDIF
1213                     WRITE(numout,*) '      ptab(ji,jj,jk), fse3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1214                     kindic_agr = kindic_agr + 1
1215                  ENDIF
1216               END DO
1217            END DO
1218         END DO
1219
1220      ENDIF
1221      !
1222   END SUBROUTINE interpe3t
1223
1224   SUBROUTINE interpumsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1225      !!----------------------------------------------------------------------
1226      !!                  ***  ROUTINE interpumsk  ***
1227      !!---------------------------------------------------------------------- 
1228      !
1229      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1230      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1231      LOGICAL :: before
1232      INTEGER, INTENT(in) :: nb , ndir
1233      !
1234      INTEGER :: ji, jj, jk
1235      LOGICAL :: western_side, eastern_side   
1236      !!---------------------------------------------------------------------- 
1237      !   
1238      IF (before) THEN
1239         DO jk=k1,k2
1240            DO jj=j1,j2
1241               DO ji=i1,i2
1242                  ptab(ji,jj,jk) = umask(ji,jj,jk)
1243               END DO
1244            END DO
1245         END DO
1246      ELSE
1247
1248         western_side  = (nb == 1).AND.(ndir == 1)
1249         eastern_side  = (nb == 1).AND.(ndir == 2)
1250         DO jk=k1,k2
1251            DO jj=j1,j2
1252               DO ji=i1,i2
1253                   ! Velocity mask at boundary edge points:
1254                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1255                     IF (western_side) THEN
1256                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1257                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1258                        kindic_agr = kindic_agr + 1
1259                     ELSEIF (eastern_side) THEN
1260                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1261                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1262                        kindic_agr = kindic_agr + 1
1263                     ENDIF
1264                  ENDIF
1265               END DO
1266            END DO
1267         END DO
1268
1269      ENDIF
1270      !
1271   END SUBROUTINE interpumsk
1272
1273   SUBROUTINE interpvmsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1274      !!----------------------------------------------------------------------
1275      !!                  ***  ROUTINE interpvmsk  ***
1276      !!---------------------------------------------------------------------- 
1277      !
1278      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1279      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1280      LOGICAL :: before
1281      INTEGER, INTENT(in) :: nb , ndir
1282      !
1283      INTEGER :: ji, jj, jk
1284      LOGICAL :: northern_side, southern_side     
1285      !!---------------------------------------------------------------------- 
1286      !   
1287      IF (before) THEN
1288         DO jk=k1,k2
1289            DO jj=j1,j2
1290               DO ji=i1,i2
1291                  ptab(ji,jj,jk) = vmask(ji,jj,jk)
1292               END DO
1293            END DO
1294         END DO
1295      ELSE
1296
1297         southern_side = (nb == 2).AND.(ndir == 1)
1298         northern_side = (nb == 2).AND.(ndir == 2)
1299         DO jk=k1,k2
1300            DO jj=j1,j2
1301               DO ji=i1,i2
1302                   ! Velocity mask at boundary edge points:
1303                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1304                     IF (southern_side) THEN
1305                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1306                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1307                        kindic_agr = kindic_agr + 1
1308                     ELSEIF (northern_side) THEN
1309                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1310                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1311                        kindic_agr = kindic_agr + 1
1312                     ENDIF
1313                  ENDIF
1314               END DO
1315            END DO
1316         END DO
1317
1318      ENDIF
1319      !
1320   END SUBROUTINE interpvmsk
1321
1322# if defined key_zdftke
1323   SUBROUTINE interpavt(ptab,i1,i2,j1,j2,k1,k2,before)
1324      !!----------------------------------------------------------------------
1325      !!                  ***  ROUTINE interavt  ***
1326      !!---------------------------------------------------------------------- 
1327      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1328      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab     
1329      LOGICAL, INTENT(in) :: before
1330      !!---------------------------------------------------------------------- 
1331      !     
1332      IF( before) THEN
1333         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1334      ELSE
1335         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1336      ENDIF
1337      !
1338     
1339   END SUBROUTINE interpavt
1340
1341
1342   SUBROUTINE interpavm(ptab,i1,i2,j1,j2,k1,k2,before)
1343      !!----------------------------------------------------------------------
1344      !!                  ***  ROUTINE interavm  ***
1345      !!---------------------------------------------------------------------- 
1346      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1347      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1348      LOGICAL, INTENT(in) :: before
1349      !!---------------------------------------------------------------------- 
1350      !     
1351      IF( before) THEN
1352         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1353      ELSE
1354         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1355      ENDIF
1356      !
1357   END SUBROUTINE interpavm
1358
1359
1360   SUBROUTINE interpavmu(ptab,i1,i2,j1,j2,k1,k2,before)
1361      !!----------------------------------------------------------------------
1362      !!                  ***  ROUTINE interavmu  ***
1363      !!---------------------------------------------------------------------- 
1364      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1365      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1366      LOGICAL, INTENT(in) :: before
1367      !!---------------------------------------------------------------------- 
1368      !
1369      IF( before) THEN
1370         ptab  (i1:i2,j1:j2,k1:k2) = avmu_k(i1:i2,j1:j2,k1:k2)
1371      ELSE
1372         avmu_k(i1:i2,j1:j2,k1:k2) = ptab   (i1:i2,j1:j2,k1:k2)
1373      ENDIF
1374      !     
1375   END SUBROUTINE interpavmu
1376
1377
1378   SUBROUTINE interpavmv(ptab,i1,i2,j1,j2,k1,k2,before)
1379      !!----------------------------------------------------------------------
1380      !!                  ***  ROUTINE interavmv  ***
1381      !!---------------------------------------------------------------------- 
1382      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1383      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1384      LOGICAL, INTENT(in) :: before
1385      !!---------------------------------------------------------------------- 
1386      !
1387      IF( before) THEN
1388         ptab  (i1:i2,j1:j2,k1:k2) = avmv_k(i1:i2,j1:j2,k1:k2)
1389      ELSE
1390         avmv_k(i1:i2,j1:j2,k1:k2) = ptab  (i1:i2,j1:j2,k1:k2)
1391      ENDIF
1392      !
1393   END SUBROUTINE interpavmv
1394# endif /* key_zdftke */
1395
1396#else
1397   !!----------------------------------------------------------------------
1398   !!   Empty module                                          no AGRIF zoom
1399   !!----------------------------------------------------------------------
1400CONTAINS
1401   SUBROUTINE Agrif_OPA_Interp_empty
1402      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1403   END SUBROUTINE Agrif_OPA_Interp_empty
1404#endif
1405
1406   !!======================================================================
1407END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.