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

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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 9 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

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