source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90 @ 12152

Last change on this file since 12152 was 12152, checked in by jchanut, 7 months ago

#2222: fixes linear vertical interpolation of viscosities

  • Property svn:keywords set to Id
File size: 44.8 KB
Line 
1MODULE agrif_oce_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_oce_interp  ***
4   !! AGRIF: interpolation package for the ocean dynamics (OPA)
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (L. Debreu)  Original cade
7   !!            3.2  !  2009-04  (R. Benshila)
8   !!            3.6  !  2014-09  (R. Benshila)
9   !!----------------------------------------------------------------------
10#if defined key_agrif
11   !!----------------------------------------------------------------------
12   !!   'key_agrif'                                              AGRIF zoom
13   !!----------------------------------------------------------------------
14   !!   Agrif_tra     :
15   !!   Agrif_dyn     :
16   !!   Agrif_ssh     :
17   !!   Agrif_dyn_ts  :
18   !!   Agrif_dta_ts  :
19   !!   Agrif_ssh_ts  :
20   !!   Agrif_avm     :
21   !!   interpu       :
22   !!   interpv       :
23   !!----------------------------------------------------------------------
24   USE par_oce
25   USE oce
26   USE dom_oce     
27   USE zdf_oce
28   USE agrif_oce
29   USE phycst
30   USE dynspg_ts, ONLY: un_adv, vn_adv
31   !
32   USE in_out_manager
33   USE agrif_oce_sponge
34   USE lib_mpp
35   USE vremap
36 
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts
41   PUBLIC   Agrif_tra, Agrif_avm
42   PUBLIC   interpun , interpvn
43   PUBLIC   interptsn, interpsshn, interpavm
44   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
45   PUBLIC   interpe3t
46#if defined key_vertical
47   PUBLIC   interpht0, interpmbkt
48# endif
49   INTEGER ::   bdy_tinterp = 0
50
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/NST 4.0 , NEMO Consortium (2018)
54   !! $Id$
55   !! Software governed by the CeCILL license (see ./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE Agrif_tra
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE Agrif_tra  ***
62      !!----------------------------------------------------------------------
63      !
64      IF( Agrif_Root() )   RETURN
65      !
66      Agrif_SpecialValue    = 0._wp
67      Agrif_UseSpecialValue = .TRUE.
68      !
69      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
70      !
71      Agrif_UseSpecialValue = .FALSE.
72      !
73   END SUBROUTINE Agrif_tra
74
75
76   SUBROUTINE Agrif_dyn( kt )
77      !!----------------------------------------------------------------------
78      !!                  ***  ROUTINE Agrif_DYN  ***
79      !!---------------------------------------------------------------------- 
80      INTEGER, INTENT(in) ::   kt
81      !
82      INTEGER ::   ji, jj, jk       ! dummy loop indices
83      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
84      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
85      !!---------------------------------------------------------------------- 
86      !
87      IF( Agrif_Root() )   RETURN
88      !
89      Agrif_SpecialValue    = 0._wp
90      Agrif_UseSpecialValue = ln_spc_dyn
91      !
92      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
93      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
94      !
95      Agrif_UseSpecialValue = .FALSE.
96      !
97      ! --- West --- !
98      ibdy1 = 2
99      ibdy2 = 1+nbghostcells 
100      !
101      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
102         DO ji = mi0(ibdy1), mi1(ibdy2)
103            ua_b(ji,:) = 0._wp
104
105            DO jk = 1, jpkm1
106               DO jj = 1, jpj
107                  ua_b(ji,jj) = ua_b(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
108               END DO
109            END DO
110
111            DO jj = 1, jpj
112               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
113            END DO
114         END DO
115      ENDIF
116      !
117      DO ji = mi0(ibdy1), mi1(ibdy2)
118         zub(ji,:) = 0._wp    ! Correct transport
119         DO jk = 1, jpkm1
120            DO jj = 1, jpj
121               zub(ji,jj) = zub(ji,jj) & 
122                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk)*umask(ji,jj,jk)
123            END DO
124         END DO
125         DO jj=1,jpj
126            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
127         END DO
128           
129         DO jk = 1, jpkm1
130            DO jj = 1, jpj
131               ua(ji,jj,jk) = ( ua(ji,jj,jk) + ua_b(ji,jj)-zub(ji,jj)) * umask(ji,jj,jk)
132            END DO
133         END DO
134      END DO
135           
136      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
137         DO ji = mi0(ibdy1), mi1(ibdy2)
138            zvb(ji,:) = 0._wp
139            DO jk = 1, jpkm1
140               DO jj = 1, jpj
141                  zvb(ji,jj) = zvb(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
142               END DO
143            END DO
144            DO jj = 1, jpj
145               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
146            END DO
147            DO jk = 1, jpkm1
148               DO jj = 1, jpj
149                  va(ji,jj,jk) = ( va(ji,jj,jk) + va_b(ji,jj)-zvb(ji,jj))*vmask(ji,jj,jk)
150               END DO
151            END DO
152         END DO
153      ENDIF
154
155      ! --- East --- !
156      ibdy1 = jpiglo-1-nbghostcells
157      ibdy2 = jpiglo-2 
158      !
159      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
160         DO ji = mi0(ibdy1), mi1(ibdy2)
161            ua_b(ji,:) = 0._wp
162            DO jk = 1, jpkm1
163               DO jj = 1, jpj
164                  ua_b(ji,jj) = ua_b(ji,jj) & 
165                      & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
166               END DO
167            END DO
168            DO jj = 1, jpj
169               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
170            END DO
171         END DO
172      ENDIF
173      !
174      DO ji = mi0(ibdy1), mi1(ibdy2)
175         zub(ji,:) = 0._wp    ! Correct transport
176         DO jk = 1, jpkm1
177            DO jj = 1, jpj
178               zub(ji,jj) = zub(ji,jj) & 
179                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk) * umask(ji,jj,jk)
180            END DO
181         END DO
182         DO jj=1,jpj
183            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
184         END DO
185           
186         DO jk = 1, jpkm1
187            DO jj = 1, jpj
188               ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
189                 & + ua_b(ji,jj)-zub(ji,jj))*umask(ji,jj,jk)
190            END DO
191         END DO
192      END DO
193           
194      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
195         ibdy1 = jpiglo-nbghostcells
196         ibdy2 = jpiglo-1 
197         DO ji = mi0(ibdy1), mi1(ibdy2)
198            zvb(ji,:) = 0._wp
199            DO jk = 1, jpkm1
200               DO jj = 1, jpj
201                  zvb(ji,jj) = zvb(ji,jj) &
202                     & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
203               END DO
204            END DO
205            DO jj = 1, jpj
206               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
207            END DO
208            DO jk = 1, jpkm1
209               DO jj = 1, jpj
210                  va(ji,jj,jk) = ( va(ji,jj,jk) & 
211                      & + va_b(ji,jj)-zvb(ji,jj)) * vmask(ji,jj,jk)
212               END DO
213            END DO
214         END DO
215      ENDIF
216
217      ! --- South --- !
218      jbdy1 = 2
219      jbdy2 = 1+nbghostcells 
220      !
221      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
222         DO jj = mj0(jbdy1), mj1(jbdy2)
223            va_b(:,jj) = 0._wp
224            DO jk = 1, jpkm1
225               DO ji = 1, jpi
226                  va_b(ji,jj) = va_b(ji,jj) & 
227                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
228               END DO
229            END DO
230            DO ji=1,jpi
231               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)     
232            END DO
233         END DO
234      ENDIF
235      !
236      DO jj = mj0(jbdy1), mj1(jbdy2)
237         zvb(:,jj) = 0._wp    ! Correct transport
238         DO jk=1,jpkm1
239            DO ji=1,jpi
240               zvb(ji,jj) = zvb(ji,jj) & 
241                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
242            END DO
243         END DO
244         DO ji = 1, jpi
245            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
246         END DO
247
248         DO jk = 1, jpkm1
249            DO ji = 1, jpi
250               va(ji,jj,jk) = ( va(ji,jj,jk) & 
251                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
252            END DO
253         END DO
254      END DO
255           
256      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
257         DO jj = mj0(jbdy1), mj1(jbdy2)
258            zub(:,jj) = 0._wp
259            DO jk = 1, jpkm1
260               DO ji = 1, jpi
261                  zub(ji,jj) = zub(ji,jj) & 
262                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
263               END DO
264            END DO
265            DO ji = 1, jpi
266               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
267            END DO
268               
269            DO jk = 1, jpkm1
270               DO ji = 1, jpi
271                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
272                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
273               END DO
274            END DO
275         END DO
276      ENDIF
277
278      ! --- North --- !
279      jbdy1 = jpjglo-1-nbghostcells
280      jbdy2 = jpjglo-2 
281      !
282      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
283         DO jj = mj0(jbdy1), mj1(jbdy2)
284            va_b(:,jj) = 0._wp
285            DO jk = 1, jpkm1
286               DO ji = 1, jpi
287                  va_b(ji,jj) = va_b(ji,jj) & 
288                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
289               END DO
290            END DO
291            DO ji=1,jpi
292               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)
293            END DO
294         END DO
295      ENDIF
296      !
297      DO jj = mj0(jbdy1), mj1(jbdy2)
298         zvb(:,jj) = 0._wp    ! Correct transport
299         DO jk=1,jpkm1
300            DO ji=1,jpi
301               zvb(ji,jj) = zvb(ji,jj) & 
302                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
303            END DO
304         END DO
305         DO ji = 1, jpi
306            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
307         END DO
308
309         DO jk = 1, jpkm1
310            DO ji = 1, jpi
311               va(ji,jj,jk) = ( va(ji,jj,jk) & 
312                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
313            END DO
314         END DO
315      END DO
316           
317      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
318         jbdy1 = jpjglo-nbghostcells
319         jbdy2 = jpjglo-1
320         DO jj = mj0(jbdy1), mj1(jbdy2)
321            zub(:,jj) = 0._wp
322            DO jk = 1, jpkm1
323               DO ji = 1, jpi
324                  zub(ji,jj) = zub(ji,jj) & 
325                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
326               END DO
327            END DO
328            DO ji = 1, jpi
329               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
330            END DO
331               
332            DO jk = 1, jpkm1
333               DO ji = 1, jpi
334                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
335                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
336               END DO
337            END DO
338         END DO
339      ENDIF
340      !
341   END SUBROUTINE Agrif_dyn
342
343
344   SUBROUTINE Agrif_dyn_ts( jn )
345      !!----------------------------------------------------------------------
346      !!                  ***  ROUTINE Agrif_dyn_ts  ***
347      !!---------------------------------------------------------------------- 
348      INTEGER, INTENT(in) ::   jn
349      !!
350      INTEGER :: ji, jj
351      INTEGER :: istart, iend, jstart, jend
352      !!---------------------------------------------------------------------- 
353      !
354      IF( Agrif_Root() )   RETURN
355      !
356      !--- West ---!
357      istart = 2
358      iend   = nbghostcells+1
359      DO ji = mi0(istart), mi1(iend)
360         DO jj=1,jpj
361            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
362            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
363         END DO
364      END DO
365      !
366      !--- East ---!
367      istart = jpiglo-nbghostcells
368      iend   = jpiglo-1
369      DO ji = mi0(istart), mi1(iend)
370         DO jj=1,jpj
371            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
372         END DO
373      END DO
374      istart = jpiglo-nbghostcells-1
375      iend   = jpiglo-2
376      DO ji = mi0(istart), mi1(iend)
377         DO jj=1,jpj
378            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
379         END DO
380      END DO
381      !
382      !--- South ---!
383      jstart = 2
384      jend   = nbghostcells+1
385      DO jj = mj0(jstart), mj1(jend)
386         DO ji=1,jpi
387            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
388            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
389         END DO
390      END DO
391      !
392      !--- North ---!
393      jstart = jpjglo-nbghostcells
394      jend   = jpjglo-1
395      DO jj = mj0(jstart), mj1(jend)
396         DO ji=1,jpi
397            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
398         END DO
399      END DO
400      jstart = jpjglo-nbghostcells-1
401      jend   = jpjglo-2
402      DO jj = mj0(jstart), mj1(jend)
403         DO ji=1,jpi
404            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
405         END DO
406      END DO
407      !
408   END SUBROUTINE Agrif_dyn_ts
409
410   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
411      !!----------------------------------------------------------------------
412      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
413      !!---------------------------------------------------------------------- 
414      INTEGER, INTENT(in) ::   jn
415      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
416      !!
417      INTEGER :: ji, jj
418      INTEGER :: istart, iend, jstart, jend
419      !!---------------------------------------------------------------------- 
420      !
421      IF( Agrif_Root() )   RETURN
422      !
423      !--- West ---!
424      istart = 2
425      iend   = nbghostcells+1
426      DO ji = mi0(istart), mi1(iend)
427         DO jj=1,jpj
428            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
429            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
430         END DO
431      END DO
432      !
433      !--- East ---!
434      istart = jpiglo-nbghostcells
435      iend   = jpiglo-1
436      DO ji = mi0(istart), mi1(iend)
437         DO jj=1,jpj
438            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
439         END DO
440      END DO
441      istart = jpiglo-nbghostcells-1
442      iend   = jpiglo-2
443      DO ji = mi0(istart), mi1(iend)
444         DO jj=1,jpj
445            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
446         END DO
447      END DO
448      !
449      !--- South ---!
450      jstart = 2
451      jend   = nbghostcells+1
452      DO jj = mj0(jstart), mj1(jend)
453         DO ji=1,jpi
454            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
455            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
456         END DO
457      END DO
458      !
459      !--- North ---!
460      jstart = jpjglo-nbghostcells
461      jend   = jpjglo-1
462      DO jj = mj0(jstart), mj1(jend)
463         DO ji=1,jpi
464            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
465         END DO
466      END DO
467      jstart = jpjglo-nbghostcells-1
468      jend   = jpjglo-2
469      DO jj = mj0(jstart), mj1(jend)
470         DO ji=1,jpi
471            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
472         END DO
473      END DO
474      !
475   END SUBROUTINE Agrif_dyn_ts_flux
476
477   SUBROUTINE Agrif_dta_ts( kt )
478      !!----------------------------------------------------------------------
479      !!                  ***  ROUTINE Agrif_dta_ts  ***
480      !!---------------------------------------------------------------------- 
481      INTEGER, INTENT(in) ::   kt
482      !!
483      INTEGER :: ji, jj
484      LOGICAL :: ll_int_cons
485      !!---------------------------------------------------------------------- 
486      !
487      IF( Agrif_Root() )   RETURN
488      !
489      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
490      !
491      ! Enforce volume conservation if no time refinement: 
492      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
493      !
494      ! Interpolate barotropic fluxes
495      Agrif_SpecialValue = 0._wp
496      Agrif_UseSpecialValue = ln_spc_dyn
497      !
498      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
499      utint_stage(:,:) = 0
500      vtint_stage(:,:) = 0
501      !
502      IF( ll_int_cons ) THEN  ! Conservative interpolation
503         ! order matters here !!!!!!
504         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
505         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
506         !
507         bdy_tinterp = 1
508         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
509         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  ) 
510         !
511         bdy_tinterp = 2
512         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
513         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )   
514      ELSE ! Linear interpolation
515         !
516         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
517         CALL Agrif_Bc_variable( unb_id, procname=interpunb )
518         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
519      ENDIF
520      Agrif_UseSpecialValue = .FALSE.
521      !
522   END SUBROUTINE Agrif_dta_ts
523
524
525   SUBROUTINE Agrif_ssh( kt )
526      !!----------------------------------------------------------------------
527      !!                  ***  ROUTINE Agrif_ssh  ***
528      !!---------------------------------------------------------------------- 
529      INTEGER, INTENT(in) ::   kt
530      !
531      INTEGER  :: ji, jj
532      INTEGER  :: istart, iend, jstart, jend
533      !!---------------------------------------------------------------------- 
534      !
535      IF( Agrif_Root() )   RETURN
536      !     
537      ! Linear time interpolation of sea level
538      !
539      Agrif_SpecialValue    = 0._wp
540      Agrif_UseSpecialValue = .TRUE.
541      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
542      Agrif_UseSpecialValue = .FALSE.
543      !
544      ! --- West --- !
545      istart = 2
546      iend   = 1 + nbghostcells
547      DO ji = mi0(istart), mi1(iend)
548         DO jj = 1, jpj
549            ssha(ji,jj) = hbdy(ji,jj)
550         ENDDO
551      ENDDO
552      !
553      ! --- East --- !
554      istart = jpiglo - nbghostcells
555      iend   = jpiglo - 1
556      DO ji = mi0(istart), mi1(iend)
557         DO jj = 1, jpj
558            ssha(ji,jj) = hbdy(ji,jj)
559         ENDDO
560      ENDDO
561      !
562      ! --- South --- !
563      jstart = 2
564      jend   = 1 + nbghostcells
565      DO jj = mj0(jstart), mj1(jend)
566         DO ji = 1, jpi
567            ssha(ji,jj) = hbdy(ji,jj)
568         ENDDO
569      ENDDO
570      !
571      ! --- North --- !
572      jstart = jpjglo - nbghostcells
573      jend   = jpjglo - 1
574      DO jj = mj0(jstart), mj1(jend)
575         DO ji = 1, jpi
576            ssha(ji,jj) = hbdy(ji,jj)
577         ENDDO
578      ENDDO
579      !
580   END SUBROUTINE Agrif_ssh
581
582
583   SUBROUTINE Agrif_ssh_ts( jn )
584      !!----------------------------------------------------------------------
585      !!                  ***  ROUTINE Agrif_ssh_ts  ***
586      !!---------------------------------------------------------------------- 
587      INTEGER, INTENT(in) ::   jn
588      !!
589      INTEGER :: ji, jj
590      INTEGER  :: istart, iend, jstart, jend
591      !!---------------------------------------------------------------------- 
592      !
593      IF( Agrif_Root() )   RETURN
594      !
595      ! --- West --- !
596      istart = 2
597      iend   = 1+nbghostcells
598      DO ji = mi0(istart), mi1(iend)
599         DO jj = 1, jpj
600            ssha_e(ji,jj) = hbdy(ji,jj)
601         ENDDO
602      ENDDO
603      !
604      ! --- East --- !
605      istart = jpiglo - nbghostcells
606      iend   = jpiglo - 1
607      DO ji = mi0(istart), mi1(iend)
608         DO jj = 1, jpj
609            ssha_e(ji,jj) = hbdy(ji,jj)
610         ENDDO
611      ENDDO
612      !
613      ! --- South --- !
614      jstart = 2
615      jend   = 1+nbghostcells
616      DO jj = mj0(jstart), mj1(jend)
617         DO ji = 1, jpi
618            ssha_e(ji,jj) = hbdy(ji,jj)
619         ENDDO
620      ENDDO
621      !
622      ! --- North --- !
623      jstart = jpjglo - nbghostcells
624      jend   = jpjglo - 1
625      DO jj = mj0(jstart), mj1(jend)
626         DO ji = 1, jpi
627            ssha_e(ji,jj) = hbdy(ji,jj)
628         ENDDO
629      ENDDO
630      !
631   END SUBROUTINE Agrif_ssh_ts
632
633   SUBROUTINE Agrif_avm
634      !!----------------------------------------------------------------------
635      !!                  ***  ROUTINE Agrif_avm  ***
636      !!---------------------------------------------------------------------- 
637      REAL(wp) ::   zalpha
638      !!---------------------------------------------------------------------- 
639      !
640      IF( Agrif_Root() )   RETURN
641      !
642      zalpha = 1._wp ! JC: proper time interpolation impossible 
643                     ! => use last available value from parent
644      !
645      Agrif_SpecialValue    = 0.e0
646      Agrif_UseSpecialValue = .TRUE.
647      !
648      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
649      !
650      Agrif_UseSpecialValue = .FALSE.
651      !
652   END SUBROUTINE Agrif_avm
653   
654
655   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
656      !!----------------------------------------------------------------------
657      !!                  *** ROUTINE interptsn ***
658      !!----------------------------------------------------------------------
659      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
660      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
661      LOGICAL                                     , INTENT(in   ) ::   before
662      !
663      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
664      INTEGER  ::   N_in, N_out
665      ! vertical interpolation:
666      REAL(wp) :: zhtot
667      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin
668      REAL(wp), DIMENSION(k1:k2) :: h_in
669      REAL(wp), DIMENSION(1:jpk) :: h_out
670      !!----------------------------------------------------------------------
671
672      IF( before ) THEN         
673         DO jn = 1,jpts
674            DO jk=k1,k2
675               DO jj=j1,j2
676                 DO ji=i1,i2
677                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
678                 END DO
679              END DO
680           END DO
681        END DO
682
683# if defined key_vertical
684        ! Interpolate thicknesses
685        ! Warning: these are masked, hence extrapolated prior interpolation.
686        DO jk=k1,k2
687           DO jj=j1,j2
688              DO ji=i1,i2
689                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
690              END DO
691           END DO
692        END DO
693
694        ! Extrapolate thicknesses in partial bottom cells:
695        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
696        IF (ln_zps) THEN
697           DO jj=j1,j2
698              DO ji=i1,i2
699                  jk = mbkt(ji,jj)
700                  ptab(ji,jj,jk,jpts+1) = 0._wp
701              END DO
702           END DO           
703        END IF
704     
705        ! Save ssh at last level:
706        IF (.NOT.ln_linssh) THEN
707           ptab(i1:i2,j1:j2,k2,jpts+1) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
708        ELSE
709           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
710        END IF     
711# endif
712      ELSE 
713
714# if defined key_vertical 
715         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
716           
717         DO jj=j1,j2
718            DO ji=i1,i2
719               tsa(ji,jj,:,:) = 0._wp
720               N_in = mbkt_parent(ji,jj)
721               zhtot = 0._wp
722               DO jk=1,N_in !k2 = jpk of parent grid
723                  IF (jk==N_in) THEN
724                     h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
725                  ELSE
726                     h_in(jk) = ptab(ji,jj,jk,n2)
727                  ENDIF
728                  zhtot = zhtot + h_in(jk)
729                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
730               END DO
731               N_out = 0
732               DO jk=1,jpk ! jpk of child grid
733                  IF (tmask(ji,jj,jk) == 0._wp) EXIT
734                  N_out = N_out + 1
735                  h_out(jk) = e3t_a(ji,jj,jk)
736               ENDDO
737               IF (N_in*N_out > 0) THEN
738                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tsa(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
739               ENDIF
740            ENDDO
741         ENDDO
742# else
743         !
744         DO jn=1, jpts
745            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
746         END DO
747# endif
748
749      ENDIF
750      !
751   END SUBROUTINE interptsn
752
753   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
754      !!----------------------------------------------------------------------
755      !!                  ***  ROUTINE interpsshn  ***
756      !!---------------------------------------------------------------------- 
757      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
758      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
759      LOGICAL                         , INTENT(in   ) ::   before
760      !
761      !!---------------------------------------------------------------------- 
762      !
763      IF( before) THEN
764         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
765      ELSE
766         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
767      ENDIF
768      !
769   END SUBROUTINE interpsshn
770
771   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
772      !!----------------------------------------------------------------------
773      !!                  *** ROUTINE interpun ***
774      !!---------------------------------------------   
775      !!
776      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
777      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
778      LOGICAL, INTENT(in) :: before
779      !!
780      INTEGER :: ji,jj,jk
781      REAL(wp) :: zrhoy, zhtot
782      ! vertical interpolation:
783      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
784      REAL(wp), DIMENSION(1:jpk) :: h_out
785      INTEGER  :: N_in, N_out
786      REAL(wp) :: h_diff
787      !!---------------------------------------------   
788      !
789      IF (before) THEN
790         DO jk=1,jpk
791            DO jj=j1,j2
792               DO ji=i1,i2
793                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
794# if defined key_vertical
795                  ! Interpolate thicknesses (masked for subsequent extrapolation)
796                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
797# endif
798               END DO
799            END DO
800         END DO
801# if defined key_vertical
802         ! Extrapolate thicknesses in partial bottom cells:
803         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
804         IF (ln_zps) THEN
805            DO jj=j1,j2
806               DO ji=i1,i2
807                  jk = mbku(ji,jj)
808                  ptab(ji,jj,jk,2) = 0._wp
809               END DO
810            END DO           
811         END IF
812        ! Save ssh at last level:
813        ptab(i1:i2,j1:j2,k2,2) = 0._wp
814        IF (.NOT.ln_linssh) THEN
815           ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
816           DO jk=1,jpk
817              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u_n(i1:i2,j1:j2,jk) * umask(i1:i2,j1:j2,jk)
818           END DO
819           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
820        END IF 
821# endif
822         !
823      ELSE
824         zrhoy = Agrif_rhoy()
825# if defined key_vertical
826! VERTICAL REFINEMENT BEGIN
827
828         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
829
830         DO ji=i1,i2
831            DO jj=j1,j2
832               ua(ji,jj,:) = 0._wp
833               N_in = mbku_parent(ji,jj)
834               zhtot = 0._wp
835               DO jk=1,N_in
836                  IF (jk==N_in) THEN
837                     h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
838                  ELSE
839                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
840                  ENDIF
841                  zhtot = zhtot + h_in(jk)
842                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
843              ENDDO
844                 
845              N_out = 0
846              DO jk=1,jpk
847                 if (umask(ji,jj,jk) == 0) EXIT
848                 N_out = N_out + 1
849                 h_out(N_out) = e3u_a(ji,jj,jk)
850              ENDDO
851              IF (N_in*N_out > 0) THEN
852                 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
853              ENDIF
854            ENDDO
855         ENDDO
856
857# else
858         DO jk = 1, jpkm1
859            DO jj=j1,j2
860               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
861            END DO
862         END DO
863# endif
864
865      ENDIF
866      !
867   END SUBROUTINE interpun
868
869   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
870      !!----------------------------------------------------------------------
871      !!                  *** ROUTINE interpvn ***
872      !!----------------------------------------------------------------------
873      !
874      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
875      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
876      LOGICAL, INTENT(in) :: before
877      !
878      INTEGER :: ji,jj,jk
879      REAL(wp) :: zrhox
880      ! vertical interpolation:
881      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
882      REAL(wp), DIMENSION(1:jpk) :: h_out
883      INTEGER  :: N_in, N_out
884      REAL(wp) :: h_diff, zhtot
885      !!---------------------------------------------   
886      !     
887      IF (before) THEN         
888         DO jk=k1,k2
889            DO jj=j1,j2
890               DO ji=i1,i2
891                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
892# if defined key_vertical
893                  ! Interpolate thicknesses (masked for subsequent extrapolation)
894                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
895# endif
896               END DO
897            END DO
898         END DO
899# if defined key_vertical
900         ! Extrapolate thicknesses in partial bottom cells:
901         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
902         IF (ln_zps) THEN
903            DO jj=j1,j2
904               DO ji=i1,i2
905                  jk = mbkv(ji,jj)
906                  ptab(ji,jj,jk,2) = 0._wp
907               END DO
908            END DO           
909         END IF
910        ! Save ssh at last level:
911        ptab(i1:i2,j1:j2,k2,2) = 0._wp
912        IF (.NOT.ln_linssh) THEN
913           ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
914           DO jk=1,jpk
915              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v_n(i1:i2,j1:j2,jk) * vmask(i1:i2,j1:j2,jk)
916           END DO
917           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
918        END IF 
919# endif
920      ELSE       
921         zrhox = Agrif_rhox()
922# if defined key_vertical
923
924         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
925
926         DO jj=j1,j2
927            DO ji=i1,i2
928               va(ji,jj,:) = 0._wp
929               N_in = mbkv_parent(ji,jj)
930               zhtot = 0._wp
931               DO jk=1,N_in
932                  IF (jk==N_in) THEN
933                     h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
934                  ELSE
935                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
936                  ENDIF
937                  zhtot = zhtot + h_in(jk)
938                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
939              ENDDO
940         
941               N_out = 0
942               DO jk=1,jpk
943                  if (vmask(ji,jj,jk) == 0) EXIT
944                  N_out = N_out + 1
945                  h_out(N_out) = e3v_a(ji,jj,jk)
946               END DO
947               IF (N_in*N_out > 0) THEN
948                  call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
949               ENDIF
950            END DO
951         END DO
952# else
953         DO jk = 1, jpkm1
954            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
955         END DO
956# endif
957      ENDIF
958      !       
959   END SUBROUTINE interpvn
960
961   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
962      !!----------------------------------------------------------------------
963      !!                  ***  ROUTINE interpunb  ***
964      !!---------------------------------------------------------------------- 
965      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
966      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
967      LOGICAL                         , INTENT(in   ) ::   before
968      !
969      INTEGER  ::   ji, jj
970      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
971      !!---------------------------------------------------------------------- 
972      !
973      IF( before ) THEN
974         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
975      ELSE
976         zrhoy = Agrif_Rhoy()
977         zrhot = Agrif_rhot()
978         ! Time indexes bounds for integration
979         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
980         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
981         !
982         DO ji = i1, i2
983            DO jj = j1, j2
984               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
985                  IF    ( utint_stage(ji,jj) == 1  ) THEN
986                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
987                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
988                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
989                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
990                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
991                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
992                     ztcoeff = 1._wp
993                  ELSE
994                     ztcoeff = 0._wp
995                  ENDIF
996                  !   
997                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
998                  !           
999                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
1000                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
1001                  ENDIF
1002                  !
1003                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
1004               ENDIF
1005            END DO
1006         END DO
1007      END IF
1008      !
1009   END SUBROUTINE interpunb
1010
1011
1012   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
1013      !!----------------------------------------------------------------------
1014      !!                  ***  ROUTINE interpvnb  ***
1015      !!---------------------------------------------------------------------- 
1016      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1017      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1018      LOGICAL                         , INTENT(in   ) ::   before
1019      !
1020      INTEGER  ::   ji, jj
1021      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1022      !!---------------------------------------------------------------------- 
1023      !
1024      IF( before ) THEN
1025         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1026      ELSE
1027         zrhox = Agrif_Rhox()
1028         zrhot = Agrif_rhot()
1029         ! Time indexes bounds for integration
1030         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1031         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
1032         !     
1033         DO ji = i1, i2
1034            DO jj = j1, j2
1035               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1036                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
1037                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1038                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1039                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
1040                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1041                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1042                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
1043                     ztcoeff = 1._wp
1044                  ELSE
1045                     ztcoeff = 0._wp
1046                  ENDIF
1047                  !   
1048                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
1049                  !           
1050                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
1051                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
1052                  ENDIF
1053                  !
1054                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
1055               ENDIF
1056            END DO
1057         END DO         
1058      ENDIF
1059      !
1060   END SUBROUTINE interpvnb
1061
1062
1063   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
1064      !!----------------------------------------------------------------------
1065      !!                  ***  ROUTINE interpub2b  ***
1066      !!---------------------------------------------------------------------- 
1067      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1068      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1069      LOGICAL                         , INTENT(in   ) ::   before
1070      !
1071      INTEGER  ::   ji,jj
1072      REAL(wp) ::   zrhot, zt0, zt1, zat
1073      !!---------------------------------------------------------------------- 
1074      IF( before ) THEN
1075         IF ( ln_bt_fw ) THEN
1076            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1077         ELSE
1078            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1079         ENDIF
1080      ELSE
1081         zrhot = Agrif_rhot()
1082         ! Time indexes bounds for integration
1083         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1084         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1085         ! Polynomial interpolation coefficients:
1086         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1087            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1088         !
1089         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1090         !
1091         ! Update interpolation stage:
1092         utint_stage(i1:i2,j1:j2) = 1
1093      ENDIF
1094      !
1095   END SUBROUTINE interpub2b
1096   
1097
1098   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
1099      !!----------------------------------------------------------------------
1100      !!                  ***  ROUTINE interpvb2b  ***
1101      !!---------------------------------------------------------------------- 
1102      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1103      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1104      LOGICAL                         , INTENT(in   ) ::   before
1105      !
1106      INTEGER ::   ji,jj
1107      REAL(wp) ::   zrhot, zt0, zt1, zat
1108      !!---------------------------------------------------------------------- 
1109      !
1110      IF( before ) THEN
1111         IF ( ln_bt_fw ) THEN
1112            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1113         ELSE
1114            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1115         ENDIF
1116      ELSE     
1117         zrhot = Agrif_rhot()
1118         ! Time indexes bounds for integration
1119         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1120         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1121         ! Polynomial interpolation coefficients:
1122         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1123            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1124         !
1125         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1126         !
1127         ! update interpolation stage:
1128         vtint_stage(i1:i2,j1:j2) = 1
1129      ENDIF
1130      !     
1131   END SUBROUTINE interpvb2b
1132
1133
1134   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before )
1135      !!----------------------------------------------------------------------
1136      !!                  ***  ROUTINE interpe3t  ***
1137      !!---------------------------------------------------------------------- 
1138      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1139      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1140      LOGICAL                              , INTENT(in   ) :: before
1141      !
1142      INTEGER :: ji, jj, jk
1143      !!---------------------------------------------------------------------- 
1144      !   
1145      IF( before ) THEN
1146         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1147      ELSE
1148         !
1149         DO jk = k1, k2
1150            DO jj = j1, j2
1151               DO ji = i1, i2
1152                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1153                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  & 
1154                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), &
1155                     &                 ji+nimpp-1, jj+njmpp-1, jk
1156                     kindic_agr = kindic_agr + 1
1157                  ENDIF
1158               END DO
1159            END DO
1160         END DO
1161         !
1162      ENDIF
1163      !
1164   END SUBROUTINE interpe3t
1165
1166
1167   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1168      !!----------------------------------------------------------------------
1169      !!                  ***  ROUTINE interavm  ***
1170      !!---------------------------------------------------------------------- 
1171      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1172      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1173      LOGICAL                                    , INTENT(in   ) ::   before
1174      !
1175      INTEGER  :: ji, jj, jk
1176      INTEGER  :: N_in, N_out
1177      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in
1178      REAL(wp), DIMENSION(1:jpk) :: z_out
1179      !!---------------------------------------------------------------------- 
1180      !     
1181      IF (before) THEN         
1182         DO jk=k1,k2
1183            DO jj=j1,j2
1184              DO ji=i1,i2
1185                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1186              END DO
1187           END DO
1188        END DO
1189
1190# if defined key_vertical
1191        ! Interpolate thicknesses
1192        ! Warning: these are masked, hence extrapolated prior interpolation.
1193        DO jk=k1,k2
1194           DO jj=j1,j2
1195              DO ji=i1,i2
1196                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
1197              END DO
1198           END DO
1199        END DO
1200
1201        ! Extrapolate thicknesses in partial bottom cells:
1202        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1203        IF (ln_zps) THEN
1204           DO jj=j1,j2
1205              DO ji=i1,i2
1206                  jk = mbkt(ji,jj)
1207                  ptab(ji,jj,jk,2) = 0._wp
1208              END DO
1209           END DO           
1210        END IF
1211     
1212        ! Save ssh at last level:
1213        IF (.NOT.ln_linssh) THEN
1214           ptab(i1:i2,j1:j2,k2,2) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
1215        ELSE
1216           ptab(i1:i2,j1:j2,k2,2) = 0._wp
1217        END IF     
1218# endif
1219      ELSE 
1220#ifdef key_vertical         
1221         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1222         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
1223           
1224         DO jj = j1, j2
1225            DO ji =i1, i2
1226               N_in = mbkt_parent(ji,jj)
1227               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
1228               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
1229               DO jk = N_in, 1, -1  ! Parent vertical grid               
1230                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
1231                    tabin(jk) = ptab(ji,jj,jk,1)
1232               END DO
1233               N_out = mbkt(ji,jj) 
1234               DO jk = 1, N_out        ! Child vertical grid
1235                  z_out(jk) = gdepw_n(ji,jj,jk)
1236               ENDDO
1237               IF (N_in*N_out > 0) THEN
1238                  CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
1239               ENDIF
1240            ENDDO
1241         ENDDO
1242#else
1243         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1244#endif
1245      ENDIF
1246      !
1247   END SUBROUTINE interpavm
1248
1249# if defined key_vertical
1250   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
1251      !!----------------------------------------------------------------------
1252      !!                  ***  ROUTINE interpsshn  ***
1253      !!---------------------------------------------------------------------- 
1254      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1255      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1256      LOGICAL                         , INTENT(in   ) ::   before
1257      !
1258      !!---------------------------------------------------------------------- 
1259      !
1260      IF( before) THEN
1261         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
1262      ELSE
1263         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
1264      ENDIF
1265      !
1266   END SUBROUTINE interpmbkt
1267
1268   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
1269      !!----------------------------------------------------------------------
1270      !!                  ***  ROUTINE interpsshn  ***
1271      !!---------------------------------------------------------------------- 
1272      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1273      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1274      LOGICAL                         , INTENT(in   ) ::   before
1275      !
1276      !!---------------------------------------------------------------------- 
1277      !
1278      IF( before) THEN
1279         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2)
1280      ELSE
1281         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
1282      ENDIF
1283      !
1284   END SUBROUTINE interpht0
1285#endif
1286
1287#else
1288   !!----------------------------------------------------------------------
1289   !!   Empty module                                          no AGRIF zoom
1290   !!----------------------------------------------------------------------
1291CONTAINS
1292   SUBROUTINE Agrif_OCE_Interp_empty
1293      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1294   END SUBROUTINE Agrif_OCE_Interp_empty
1295#endif
1296
1297   !!======================================================================
1298END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.