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

source: branches/2015/dev_r5847_MERCATOR9_solveur_simplification/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 5868

Last change on this file since 5868 was 5868, checked in by jchanut, 8 years ago

Free surface simplification #1620. Step 1: suppress filtered free surface

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