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

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

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

Last change on this file since 4785 was 4785, checked in by rblod, 10 years ago

dev_r4765_CNRS_agrif: First update of AGRIF for dynamic only (_flt and _ts), see ticket #1380 and associated wiki page

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