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_oce_interp.F90 in NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_interp.F90 @ 13937

Last change on this file since 13937 was 13937, checked in by jchanut, 4 years ago

#2222, added:

  • Divergence conserving interpolation for transport (Balsara's scheme)
  • Possibility to shift barotropic interface inside the domain (currently disabled - needed to match volume within sponge)
  • PPR library as default vertical interpolation scheme
  • Trilinear interpolation over partial bottom cells
  • Property svn:keywords set to Id
File size: 73.5 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   USE lbclnk
37 
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts
42   PUBLIC   Agrif_tra, Agrif_avm
43   PUBLIC   interpun , interpvn
44   PUBLIC   interptsn, interpsshn, interpavm
45   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
46   PUBLIC   interpe3t, interpglamt, interpgphit
47   PUBLIC   interpht0, interpmbkt, interpe3t0_vremap
48   PUBLIC   agrif_istate_oce, agrif_istate_ssh   ! called by icestate.F90 and domvvl.F90
49   PUBLIC   agrif_check_bat
50
51   INTEGER ::   bdy_tinterp = 0
52
53   !!----------------------------------------------------------------------
54   !! NEMO/NST 4.0 , NEMO Consortium (2018)
55   !! $Id$
56   !! Software governed by the CeCILL license (see ./LICENSE)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE Agrif_istate_oce( Kbb, Kmm, Kaa )
61      !!----------------------------------------------------------------------
62      !!                 *** ROUTINE agrif_istate_oce ***
63      !!
64      !!                 set initial t, s, u, v, ssh from parent
65      !!----------------------------------------------------------------------
66      !
67      IMPLICIT NONE
68      !
69      INTEGER, INTENT(in)  :: Kbb, Kmm, Kaa
70      INTEGER :: jn
71      !!----------------------------------------------------------------------
72      IF(lwp) WRITE(numout,*) ' '
73      IF(lwp) WRITE(numout,*) 'Agrif_istate_oce : interp child initial state from parent'
74      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~'
75      IF(lwp) WRITE(numout,*) ' '
76
77      IF ( ln_rstart ) & 
78         & CALL ctl_stop('AGRIF hot start should be desactivated in restarting mode')
79
80      IF ( .NOT.Agrif_Parent(l_1st_euler) ) & 
81         & CALL ctl_stop('AGRIF hot start requires to force Euler first step on parent')
82
83      l_ini_child           = .TRUE.
84      Agrif_SpecialValue    = 0.0_wp
85      Agrif_UseSpecialValue = .TRUE.
86
87      ts(:,:,:,:,:) = 0.0_wp
88      uu(:,:,:,:)   = 0.0_wp
89      vv(:,:,:,:)   = 0.0_wp 
90      ssh(:,:,:)    = 0._wp
91       
92      Krhs_a = Kbb   ;   Kmm_a = Kbb
93
94      CALL Agrif_Init_Variable(tsini_id, procname=interptsn)
95      CALL Agrif_Init_Variable(sshini_id, procname=interpsshn)
96
97      Agrif_UseSpecialValue = ln_spc_dyn
98      use_sign_north = .TRUE.
99      sign_north = -1._wp
100      CALL Agrif_Init_Variable(uini_id , procname=interpun )
101      CALL Agrif_Init_Variable(vini_id , procname=interpvn )
102      use_sign_north = .FALSE.
103
104      Agrif_UseSpecialValue = .FALSE.
105      l_ini_child           = .FALSE.
106
107      Krhs_a = Kaa   ;   Kmm_a = Kmm
108
109      ssh(:,:,Kbb) = ssh(:,:,Kbb) * tmask(:,:,1)
110
111      DO jn = 1, jpts
112         ts(:,:,:,jn,Kbb) = ts(:,:,:,jn,Kbb) * tmask(:,:,:)
113      END DO
114      uu(:,:,:,Kbb) = uu(:,:,:,Kbb) * umask(:,:,:)     
115      vv(:,:,:,Kbb) = vv(:,:,:,Kbb) * vmask(:,:,:) 
116
117      CALL lbc_lnk_multi( 'agrif_istate_oce', uu(:,:,:  ,Kbb), 'U', -1.0_wp , vv(:,:,:,Kbb), 'V', -1.0_wp )
118      CALL lbc_lnk( 'agrif_istate_oce', ts(:,:,:,:,Kbb), 'T',  1.0_wp )
119      CALL lbc_lnk( 'agrif_istate_oce', ssh(:,:,Kbb), 'T', 1.0_wp )
120
121   END SUBROUTINE Agrif_istate_oce
122
123
124   SUBROUTINE Agrif_istate_ssh( Kbb, Kmm )
125      !!----------------------------------------------------------------------
126      !!                 *** ROUTINE agrif_istate_ssh ***
127      !!
128      !!                    set initial ssh from parent
129      !!----------------------------------------------------------------------
130      !
131      IMPLICIT NONE
132      !
133      INTEGER, INTENT(in)  :: Kbb, Kmm 
134      !!----------------------------------------------------------------------
135      IF(lwp) WRITE(numout,*) ' '
136      IF(lwp) WRITE(numout,*) 'Agrif_istate_ssh : interp child ssh from parent'
137      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~'
138      IF(lwp) WRITE(numout,*) ' '
139
140      IF ( ln_rstart ) & 
141         & CALL ctl_stop('AGRIF hot start should be desactivated in restarting mode')
142
143      IF ( .NOT.Agrif_Parent(l_1st_euler) ) & 
144         & CALL ctl_stop('AGRIF hot start requires to force Euler first step on parent')
145
146      Kmm_a = Kmm
147      ssh(:,:,Kmm) = 0._wp
148
149      Agrif_SpecialValue    = 0._wp
150      Agrif_UseSpecialValue = .TRUE.
151      l_ini_child           = .TRUE.
152      !
153      CALL Agrif_Init_Variable(sshini_id, procname=interpsshn)
154      !
155      Agrif_UseSpecialValue = .FALSE.
156      l_ini_child           = .FALSE.
157      CALL lbc_lnk( 'dom_vvl_rst', ssh(:,:,Kmm), 'T', 1._wp )
158
159   END SUBROUTINE Agrif_istate_ssh
160
161
162   SUBROUTINE Agrif_tra
163      !!----------------------------------------------------------------------
164      !!                  ***  ROUTINE Agrif_tra  ***
165      !!----------------------------------------------------------------------
166      !
167      IF( Agrif_Root() )   RETURN
168      !
169      Agrif_SpecialValue    = 0._wp
170      Agrif_UseSpecialValue = .TRUE.
171      l_vremap           = ln_vert_remap
172      !
173      CALL Agrif_Bc_variable( ts_interp_id, procname=interptsn )
174      !
175      Agrif_UseSpecialValue = .FALSE.
176      l_vremap              = .FALSE.
177      !
178   END SUBROUTINE Agrif_tra
179
180
181   SUBROUTINE Agrif_dyn( kt )
182      !!----------------------------------------------------------------------
183      !!                  ***  ROUTINE Agrif_DYN  ***
184      !!---------------------------------------------------------------------- 
185      INTEGER, INTENT(in) ::   kt
186      !
187      INTEGER ::   ji, jj, jk       ! dummy loop indices
188      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
189      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
190      !!---------------------------------------------------------------------- 
191      !
192      IF( Agrif_Root() )   RETURN
193      !
194      Agrif_SpecialValue    = 0.0_wp
195      Agrif_UseSpecialValue = ln_spc_dyn
196      l_vremap              = ln_vert_remap
197      !
198      use_sign_north = .TRUE.
199      sign_north = -1.0_wp
200      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
201      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
202
203      IF( .NOT.ln_dynspg_ts ) THEN ! Get transports
204         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
205         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb )
206         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb )
207      ENDIF
208
209      use_sign_north = .FALSE.
210      !
211      Agrif_UseSpecialValue = .FALSE.
212      l_vremap              = .FALSE.
213      !
214      ! Ensure below that vertically integrated transports match
215      ! either transports out of time splitting procedure (ln_dynspg_ts=.TRUE.)
216      ! or parent grid transports (ln_dynspg_ts=.FALSE.)
217      !
218      ! --- West --- !
219      IF( lk_west ) THEN
220         ibdy1 = nn_hls + 2                  ! halo + land + 1
221         ibdy2 = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()   ! halo + land + nbghostcells
222         !
223         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
224            DO ji = mi0(ibdy1), mi1(ibdy2)
225               DO jj = 1, jpj
226                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
227                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
228               END DO
229            END DO
230         ENDIF
231         !
232         DO ji = mi0(ibdy1), mi1(ibdy2)
233            zub(ji,:) = 0._wp 
234            DO jk = 1, jpkm1
235               DO jj = 1, jpj
236                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
237               END DO
238            END DO
239            DO jj=1,jpj
240               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
241            END DO
242            DO jk = 1, jpkm1
243               DO jj = 1, jpj
244                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
245               END DO
246            END DO
247         END DO
248         !   
249         DO ji = mi0(ibdy1), mi1(ibdy2)
250            zvb(ji,:) = 0._wp
251            DO jk = 1, jpkm1
252               DO jj = 1, jpj
253                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
254               END DO
255            END DO
256            DO jj = 1, jpj
257               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
258            END DO
259            DO jk = 1, jpkm1
260               DO jj = 1, jpj
261                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) )*vmask(ji,jj,jk)
262               END DO
263            END DO
264         END DO
265         !
266      ENDIF
267
268      ! --- East --- !
269      IF( lk_east) THEN
270         ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox()   
271         ibdy2 = jpiglo - ( nn_hls + 2 )                 
272         !
273         IF( .NOT.ln_dynspg_ts ) THEN
274            DO ji = mi0(ibdy1), mi1(ibdy2)
275               uu_b(ji,:,Krhs_a) = 0._wp
276               DO jk = 1, jpkm1
277                  DO jj = 1, jpj
278                     uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
279                  END DO
280               END DO
281               DO jj = 1, jpj
282                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
283               END DO
284            END DO
285         ENDIF
286         !
287         DO ji = mi0(ibdy1), mi1(ibdy2)
288            zub(ji,:) = 0._wp   
289            DO jk = 1, jpkm1
290               DO jj = 1, jpj
291                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
292               END DO
293            END DO
294            DO jj=1,jpj
295               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
296            END DO
297            DO jk = 1, jpkm1
298               DO jj = 1, jpj
299                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
300               END DO
301            END DO
302         END DO
303         !
304         ibdy1 = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox() 
305         ibdy2 = jpiglo - ( nn_hls + 1 )              !
306         IF( .NOT.ln_dynspg_ts ) THEN
307            DO ji = mi0(ibdy1), mi1(ibdy2)
308               vv_b(ji,:,Krhs_a) = 0._wp
309               DO jk = 1, jpkm1
310                  DO jj = 1, jpj
311                     vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
312                  END DO
313               END DO
314               DO jj = 1, jpj
315                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
316               END DO
317            END DO
318         ENDIF
319
320         DO ji = mi0(ibdy1), mi1(ibdy2)
321            zvb(ji,:) = 0._wp
322            DO jk = 1, jpkm1
323               DO jj = 1, jpj
324                     zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
325               END DO
326            END DO
327            DO jj = 1, jpj
328               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
329            END DO
330            DO jk = 1, jpkm1
331               DO jj = 1, jpj
332                     vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
333               END DO
334            END DO
335         END DO
336         !
337      ENDIF
338
339      ! --- South --- !
340      IF( lk_south ) THEN
341         jbdy1 = nn_hls + 2                 
342         jbdy2 = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()   
343         !
344         IF( .NOT.ln_dynspg_ts ) THEN
345            DO jj = mj0(jbdy1), mj1(jbdy2)
346               vv_b(:,jj,Krhs_a) = 0._wp
347               DO jk = 1, jpkm1
348                  DO ji = 1, jpi
349                     vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
350                  END DO
351               END DO
352               DO ji=1,jpi
353                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 
354                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a)
355               END DO
356            END DO
357         ENDIF
358         !
359         DO jj = mj0(jbdy1), mj1(jbdy2)
360            zvb(:,jj) = 0._wp
361            DO jk=1,jpkm1
362               DO ji=1,jpi
363                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
364               END DO
365            END DO
366            DO ji = 1, jpi
367               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
368            END DO
369            DO jk = 1, jpkm1
370               DO ji = 1, jpi
371                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
372               END DO
373            END DO
374         END DO
375         !
376         DO jj = mj0(jbdy1), mj1(jbdy2)
377            zub(:,jj) = 0._wp
378            DO jk = 1, jpkm1
379               DO ji = 1, jpi
380                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
381               END DO
382            END DO
383            DO ji = 1, jpi
384               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
385            END DO
386            DO jk = 1, jpkm1
387               DO ji = 1, jpi
388                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
389               END DO
390            END DO
391         END DO
392         !
393      ENDIF
394
395      ! --- North --- !
396      IF( lk_north ) THEN
397         jbdy1 = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy() 
398         jbdy2 = jpjglo - ( nn_hls + 2 )
399         !
400         IF( .NOT.ln_dynspg_ts ) THEN
401            DO jj = mj0(jbdy1), mj1(jbdy2)
402               vv_b(:,jj,Krhs_a) = 0._wp
403               DO jk = 1, jpkm1
404                  DO ji = 1, jpi
405                     vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
406                  END DO
407               END DO
408               DO ji=1,jpi
409                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
410               END DO
411            END DO
412         ENDIF
413         !
414         DO jj = mj0(jbdy1), mj1(jbdy2)
415            zvb(:,jj) = 0._wp 
416            DO jk=1,jpkm1
417               DO ji=1,jpi
418                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
419               END DO
420            END DO
421            DO ji = 1, jpi
422               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
423            END DO
424            DO jk = 1, jpkm1
425               DO ji = 1, jpi
426                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
427               END DO
428            END DO
429         END DO
430         !
431         jbdy1 = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy() 
432         jbdy2 = jpjglo - ( nn_hls + 1 )
433         IF( .NOT.ln_dynspg_ts ) THEN
434            DO jj = mj0(jbdy1), mj1(jbdy2)
435               uu_b(:,jj,Krhs_a) = 0._wp
436               DO jk = 1, jpkm1
437                  DO ji = 1, jpi
438                     uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
439                  END DO
440               END DO
441               DO ji=1,jpi
442                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
443               END DO
444            END DO
445         ENDIF
446
447         DO jj = mj0(jbdy1), mj1(jbdy2)
448            zub(:,jj) = 0._wp
449            DO jk = 1, jpkm1
450               DO ji = 1, jpi
451                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
452               END DO
453            END DO
454            DO ji = 1, jpi
455               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
456            END DO
457            DO jk = 1, jpkm1
458               DO ji = 1, jpi
459                     uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
460               END DO
461            END DO
462         END DO
463         !
464      ENDIF
465      !
466   END SUBROUTINE Agrif_dyn
467
468
469   SUBROUTINE Agrif_dyn_ts( jn )
470      !!----------------------------------------------------------------------
471      !!                  ***  ROUTINE Agrif_dyn_ts  ***
472      !!---------------------------------------------------------------------- 
473      INTEGER, INTENT(in) ::   jn
474      !!
475      INTEGER :: ji, jj
476      INTEGER :: istart, iend, jstart, jend
477      !!---------------------------------------------------------------------- 
478      !
479      IF( Agrif_Root() )   RETURN
480      !
481      !--- West ---!
482      IF( lk_west ) THEN
483         istart = nn_hls + 2                              ! halo + land + 1
484         iend   = nn_hls + 1 + nbghostcells  + nn_shift_bar*Agrif_Rhox()              ! halo + land + nbghostcells
485         DO ji = mi0(istart), mi1(iend)
486            DO jj=1,jpj
487               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
488               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
489            END DO
490         END DO
491      ENDIF
492      !
493      !--- East ---!
494      IF( lk_east ) THEN
495         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox() 
496         iend   = jpiglo - ( nn_hls + 1 )               
497         DO ji = mi0(istart), mi1(iend)
498
499            DO jj=1,jpj
500               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
501            END DO
502         END DO
503         istart = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox() 
504         iend   = jpiglo - ( nn_hls + 2 )               
505         DO ji = mi0(istart), mi1(iend)
506            DO jj=1,jpj
507               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
508            END DO
509         END DO
510      ENDIF 
511      !
512      !--- South ---!
513      IF( lk_south ) THEN
514         jstart = nn_hls + 2                             
515         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()           
516         DO jj = mj0(jstart), mj1(jend)
517
518            DO ji=1,jpi
519               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
520               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
521            END DO
522         END DO
523      ENDIF       
524      !
525      !--- North ---!
526      IF( lk_north ) THEN
527         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()     
528         jend   = jpjglo - ( nn_hls + 1 )               
529         DO jj = mj0(jstart), mj1(jend)
530            DO ji=1,jpi
531               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
532            END DO
533         END DO
534         jstart = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy() 
535         jend   = jpjglo - ( nn_hls + 2 )               
536         DO jj = mj0(jstart), mj1(jend)
537            DO ji=1,jpi
538               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
539            END DO
540         END DO
541      ENDIF 
542      !
543   END SUBROUTINE Agrif_dyn_ts
544
545   
546   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
547      !!----------------------------------------------------------------------
548      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
549      !!---------------------------------------------------------------------- 
550      INTEGER, INTENT(in) ::   jn
551      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
552      !!
553      INTEGER :: ji, jj
554      INTEGER :: istart, iend, jstart, jend
555      !!---------------------------------------------------------------------- 
556      !
557      IF( Agrif_Root() )   RETURN
558      !
559      !--- West ---!
560      IF( lk_west ) THEN
561         istart = nn_hls + 2                             
562         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox() 
563         DO ji = mi0(istart), mi1(iend)
564            DO jj=1,jpj
565               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
566               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
567            END DO
568         END DO
569      ENDIF
570      !
571      !--- East ---!
572      IF( lk_east ) THEN
573         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()
574         iend   = jpiglo - ( nn_hls + 1 )                 
575         DO ji = mi0(istart), mi1(iend)
576            DO jj=1,jpj
577               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
578            END DO
579         END DO
580         istart = jpiglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhox() 
581         iend   = jpiglo - ( nn_hls + 2 )                 
582         DO ji = mi0(istart), mi1(iend)
583            DO jj=1,jpj
584               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
585            END DO
586         END DO
587      ENDIF
588      !
589      !--- South ---!
590      IF( lk_south ) THEN
591         jstart = nn_hls + 2                             
592         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy() 
593         DO jj = mj0(jstart), mj1(jend)
594            DO ji=1,jpi
595               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
596               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
597            END DO
598         END DO
599      ENDIF
600      !
601      !--- North ---!
602      IF( lk_north ) THEN
603         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy() 
604         jend   = jpjglo - ( nn_hls + 1 )               
605         DO jj = mj0(jstart), mj1(jend)
606            DO ji=1,jpi
607               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
608            END DO
609         END DO
610         jstart = jpjglo - ( nn_hls + nbghostcells + 1) - nn_shift_bar*Agrif_Rhoy() 
611         jend   = jpjglo - ( nn_hls + 2 )               
612         DO jj = mj0(jstart), mj1(jend)
613            DO ji=1,jpi
614               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
615            END DO
616         END DO
617      ENDIF
618      !
619   END SUBROUTINE Agrif_dyn_ts_flux
620
621   
622   SUBROUTINE Agrif_dta_ts( kt )
623      !!----------------------------------------------------------------------
624      !!                  ***  ROUTINE Agrif_dta_ts  ***
625      !!---------------------------------------------------------------------- 
626      INTEGER, INTENT(in) ::   kt
627      !!
628      LOGICAL :: ll_int_cons
629      !!---------------------------------------------------------------------- 
630      !
631      IF( Agrif_Root() )   RETURN
632      !
633      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
634      !
635      ! Enforce volume conservation if no time refinement: 
636      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
637      !
638      ! Interpolate barotropic fluxes
639      Agrif_SpecialValue = 0._wp
640      Agrif_UseSpecialValue = ln_spc_dyn
641
642      use_sign_north = .TRUE.
643      sign_north = -1.
644
645      !
646      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
647      utint_stage(:,:) = 0
648      vtint_stage(:,:) = 0
649      !
650      IF( ll_int_cons ) THEN  ! Conservative interpolation
651         IF ( lk_tint2d_notinterp ) THEN
652            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b_const )
653            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b_const ) 
654            ! Divergence conserving correction terms:
655            CALL Agrif_Bc_variable(    ub2b_cor_id, calledweight=1._wp, procname=        ub2b_cor )
656            CALL Agrif_Bc_variable(    vb2b_cor_id, calledweight=1._wp, procname=        vb2b_cor )
657         ELSE
658            ! order matters here !!!!!!
659            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
660            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
661            !
662            bdy_tinterp = 1
663            CALL Agrif_Bc_variable( unb_interp_id , calledweight=1._wp, procname=interpunb  ) ! After
664            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=1._wp, procname=interpvnb  ) 
665            !
666            bdy_tinterp = 2
667            CALL Agrif_Bc_variable( unb_interp_id , calledweight=0._wp, procname=interpunb  ) ! Before
668            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=0._wp, procname=interpvnb  )   
669         ENDIF
670      ELSE ! Linear interpolation
671         !
672         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
673         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb )
674         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb )
675      ENDIF
676      Agrif_UseSpecialValue = .FALSE.
677      use_sign_north = .FALSE.
678      !
679   END SUBROUTINE Agrif_dta_ts
680
681
682   SUBROUTINE Agrif_ssh( kt )
683      !!----------------------------------------------------------------------
684      !!                  ***  ROUTINE Agrif_ssh  ***
685      !!---------------------------------------------------------------------- 
686      INTEGER, INTENT(in) ::   kt
687      !
688      INTEGER  :: ji, jj
689      INTEGER  :: istart, iend, jstart, jend
690      !!---------------------------------------------------------------------- 
691      !
692      IF( Agrif_Root() )   RETURN
693      !     
694      ! Linear time interpolation of sea level
695      !
696      Agrif_SpecialValue    = 0._wp
697      Agrif_UseSpecialValue = .TRUE.
698      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
699      Agrif_UseSpecialValue = .FALSE.
700      !
701      ! --- West --- !
702      IF(lk_west) THEN
703         istart = nn_hls + 2                                                          ! halo + land + 1
704         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()               ! halo + land + nbghostcells
705         DO ji = mi0(istart), mi1(iend)
706            DO jj = 1, jpj
707               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
708            END DO
709         END DO
710      ENDIF
711      !
712      ! --- East --- !
713      IF(lk_east) THEN
714         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()       ! halo + land + nbghostcells - 1
715         iend   = jpiglo - ( nn_hls + 1 )                                              ! halo + land + 1            - 1
716         DO ji = mi0(istart), mi1(iend)
717            DO jj = 1, jpj
718               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
719            END DO
720         END DO
721      ENDIF
722      !
723      ! --- South --- !
724      IF(lk_south) THEN
725         jstart = nn_hls + 2                                                          ! halo + land + 1
726         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()               ! halo + land + nbghostcells
727         DO jj = mj0(jstart), mj1(jend)
728            DO ji = 1, jpi
729               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
730            END DO
731         END DO
732      ENDIF
733      !
734      ! --- North --- !
735      IF(lk_north) THEN
736         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()     ! halo + land + nbghostcells - 1
737         jend   = jpjglo - ( nn_hls + 1 )                                            ! halo + land + 1            - 1
738         DO jj = mj0(jstart), mj1(jend)
739            DO ji = 1, jpi
740               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
741            END DO
742         END DO
743      ENDIF
744      !
745   END SUBROUTINE Agrif_ssh
746
747
748   SUBROUTINE Agrif_ssh_ts( jn )
749      !!----------------------------------------------------------------------
750      !!                  ***  ROUTINE Agrif_ssh_ts  ***
751      !!---------------------------------------------------------------------- 
752      INTEGER, INTENT(in) ::   jn
753      !!
754      INTEGER :: ji, jj
755      INTEGER  :: istart, iend, jstart, jend
756      !!---------------------------------------------------------------------- 
757      !
758      IF( Agrif_Root() )   RETURN
759      !
760      ! --- West --- !
761      IF(lk_west) THEN
762         istart = nn_hls + 2                                                        ! halo + land + 1
763         iend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhox()             ! halo + land + nbghostcells
764         DO ji = mi0(istart), mi1(iend)
765            DO jj = 1, jpj
766               ssha_e(ji,jj) = hbdy(ji,jj)
767            END DO
768         END DO
769      ENDIF
770      !
771      ! --- East --- !
772      IF(lk_east) THEN
773         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()    ! halo + land + nbghostcells - 1
774         iend   = jpiglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1
775         DO ji = mi0(istart), mi1(iend)
776            DO jj = 1, jpj
777               ssha_e(ji,jj) = hbdy(ji,jj)
778            END DO
779         END DO
780      ENDIF
781      !
782      ! --- South --- !
783      IF(lk_south) THEN
784         jstart = nn_hls + 2                                                        ! halo + land + 1
785         jend   = nn_hls + 1 + nbghostcells + nn_shift_bar*Agrif_Rhoy()             ! halo + land + nbghostcells
786         DO jj = mj0(jstart), mj1(jend)
787            DO ji = 1, jpi
788               ssha_e(ji,jj) = hbdy(ji,jj)
789            END DO
790         END DO
791      ENDIF
792      !
793      ! --- North --- !
794      IF(lk_north) THEN
795         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy()    ! halo + land + nbghostcells - 1
796         jend   = jpjglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1
797         DO jj = mj0(jstart), mj1(jend)
798            DO ji = 1, jpi
799               ssha_e(ji,jj) = hbdy(ji,jj)
800            END DO
801         END DO
802      ENDIF
803      !
804   END SUBROUTINE Agrif_ssh_ts
805
806   
807   SUBROUTINE Agrif_avm
808      !!----------------------------------------------------------------------
809      !!                  ***  ROUTINE Agrif_avm  ***
810      !!---------------------------------------------------------------------- 
811      REAL(wp) ::   zalpha
812      !!---------------------------------------------------------------------- 
813      !
814      IF( Agrif_Root() )   RETURN
815      !
816      zalpha = 1._wp ! JC: proper time interpolation impossible 
817                     ! => use last available value from parent
818      !
819      Agrif_SpecialValue    = 0.e0
820      Agrif_UseSpecialValue = .TRUE.
821      l_vremap              = ln_vert_remap
822      !
823      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
824      !
825      Agrif_UseSpecialValue = .FALSE.
826      l_vremap              = .FALSE.
827      !
828   END SUBROUTINE Agrif_avm
829
830
831   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
832      !!----------------------------------------------------------------------
833      !!                  *** ROUTINE interptsn ***
834      !!----------------------------------------------------------------------
835      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
836      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
837      LOGICAL                                     , INTENT(in   ) ::   before
838      !
839      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
840      INTEGER  ::   N_in, N_out
841      INTEGER  :: item
842      ! vertical interpolation:
843      REAL(wp) :: zhtot, zwgt
844      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin, tabin_i
845      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i
846      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
847      !!----------------------------------------------------------------------
848
849      IF( before ) THEN
850
851         item = Kmm_a
852         IF( l_ini_child )   Kmm_a = Kbb_a 
853
854         DO jn = 1,jpts
855            DO jk=k1,k2
856               DO jj=j1,j2
857                 DO ji=i1,i2
858                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)
859                 END DO
860              END DO
861           END DO
862         END DO
863
864         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN
865
866            ! Fill cell depths (i.e. gdept) to be interpolated
867            ! Warning: these are masked, hence extrapolated prior interpolation.
868            DO jj=j1,j2
869               DO ji=i1,i2
870                  ptab(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a)
871                  DO jk=k1+1,k2
872                     ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
873                        & ( ptab(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) )
874                  END DO
875               END DO
876            END DO
877         
878            ! Save ssh at last level:
879            IF (.NOT.ln_linssh) THEN
880               ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
881            END IF     
882         ENDIF
883         Kmm_a = item
884
885      ELSE
886         item = Krhs_a
887         IF( l_ini_child )   Krhs_a = Kbb_a 
888
889         IF( l_vremap .OR. l_ini_child ) THEN
890            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
891            DO jj=j1,j2
892               DO ji=i1,i2
893                  ts(ji,jj,:,:,Krhs_a) = 0. 
894                  !
895                  ! Build vertical grids:
896                  N_in = mbkt_parent(ji,jj)
897                  ! Input grid (account for partial cells if any):
898                  DO jk=1,N_in
899                     z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2)
900                     tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
901                  END DO
902                 
903                  ! Intermediate grid:
904                  DO jk = 1, N_in
905                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
906                       &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
907                  END DO
908                  z_in_i(1) = 0.5_wp * h_in_i(1)
909                  DO jk=2,N_in
910                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
911                  END DO
912                  z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2)
913
914                  ! Output (Child) grid:
915                  N_out = mbkt(ji,jj)
916                  DO jk=1,N_out
917                     h_out(jk) = e3t(ji,jj,jk,Krhs_a)
918                  END DO
919                  z_out(1) = 0.5_wp * h_out(1)
920                  DO jk=2,N_out
921                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
922                  END DO
923                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)
924
925                  IF (N_in*N_out > 0) THEN
926                     IF( l_ini_child ) THEN
927                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),              &
928                                      &   z_out(1:N_out),N_in,N_out,jpts) 
929                     ELSE
930                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),                       &
931                                     &   z_in_i(1:N_in),N_in,N_in,jpts)
932                        CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   &
933                                      &   h_out(1:N_out),N_in,N_out,jpts) 
934                     ENDIF
935                  ENDIF
936               END DO
937            END DO
938            Krhs_a = item
939 
940         ELSE
941         
942            DO jn =1, jpts
943
944               IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
945                                                ! linear vertical interpolation
946                  DO jj=j1,j2
947                     DO ji=i1,i2
948                        !
949                        N_in  = mbkt(ji,jj)
950                        N_out = mbkt(ji,jj)
951                        z_in(1) = ptab(ji,jj,1,n2)
952                        tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts)
953                        DO jk=2, N_in
954                           z_in(jk) = ptab(ji,jj,jk,n2)
955                           tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
956                        END DO
957                        IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
958                        z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a)
959                        DO jk=2, N_out
960                           z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a))
961                        END DO
962                        IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
963                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), &
964                                      &   z_out(1:N_out),N_in,N_out,jpts) 
965                     END DO
966                  END DO
967
968               ENDIF
969
970               ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a) = ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)
971
972            END DO
973         ENDIF
974
975      ENDIF
976      !
977   END SUBROUTINE interptsn
978
979   
980   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
981      !!----------------------------------------------------------------------
982      !!                  ***  ROUTINE interpsshn  ***
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      !
988      !!---------------------------------------------------------------------- 
989      !
990      IF( before) THEN
991         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
992      ELSE
993         IF( l_ini_child ) THEN
994            ssh(i1:i2,j1:j2,Kmm_a) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
995         ELSE
996            hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
997         ENDIF
998      ENDIF
999      !
1000   END SUBROUTINE interpsshn
1001
1002   
1003   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1004      !!----------------------------------------------------------------------
1005      !!                  *** ROUTINE interpun ***
1006      !!---------------------------------------------   
1007      !!
1008      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
1009      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1010      LOGICAL, INTENT(in) :: before
1011      !!
1012      INTEGER :: ji,jj,jk
1013      REAL(wp) :: zrhoy, zhtot
1014      ! vertical interpolation:
1015      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
1016      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
1017      INTEGER  :: N_in, N_out,item
1018      REAL(wp) :: h_diff
1019      !!---------------------------------------------   
1020      !
1021      IF (before) THEN
1022
1023         item = Kmm_a
1024         IF( l_ini_child )   Kmm_a = Kbb_a     
1025
1026         DO jk=1,jpk
1027            DO jj=j1,j2
1028               DO ji=i1,i2
1029                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
1030                  IF( l_vremap .OR. l_ini_child) THEN
1031                     ! Interpolate thicknesses (masked for subsequent extrapolation)
1032                     ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
1033                  ENDIF
1034               END DO
1035            END DO
1036         END DO
1037
1038        IF( l_vremap .OR. l_ini_child ) THEN
1039         ! Extrapolate thicknesses in partial bottom cells:
1040         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1041            IF (ln_zps) THEN
1042               DO jj=j1,j2
1043                  DO ji=i1,i2
1044                     jk = mbku(ji,jj)
1045                     ptab(ji,jj,jk,2) = 0._wp
1046                  END DO
1047               END DO           
1048            END IF
1049
1050           ! Save ssh at last level:
1051           ptab(i1:i2,j1:j2,k2,2) = 0._wp
1052           IF (.NOT.ln_linssh) THEN
1053              ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
1054              DO jk=1,jpk
1055                 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk)
1056              END DO
1057              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
1058           END IF
1059        ENDIF
1060
1061         Kmm_a = item
1062         !
1063      ELSE
1064         zrhoy = Agrif_rhoy()
1065
1066        IF( l_vremap .OR. l_ini_child) THEN
1067! VERTICAL REFINEMENT BEGIN
1068
1069            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1070
1071            DO ji=i1,i2
1072               DO jj=j1,j2
1073                  uu(ji,jj,:,Krhs_a) = 0._wp
1074                  N_in = mbku_parent(ji,jj)
1075                  zhtot = 0._wp
1076                  DO jk=1,N_in
1077                     !IF (jk==N_in) THEN
1078                     !   h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1079                     !ELSE
1080                     !   h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)
1081                     !ENDIF
1082                     h_in(jk) = e3u0_parent(ji,jj,jk) 
1083                     zhtot = zhtot + h_in(jk)
1084                     IF( h_in(jk) .GT. 0. ) THEN
1085                     tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
1086                     ELSE
1087                     tabin(jk) = 0.
1088                     ENDIF
1089                 END DO
1090                 z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj) 
1091                 DO jk=2,N_in
1092                    z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk)+h_in(jk-1))
1093                 END DO
1094                     
1095                 N_out = 0
1096                 DO jk=1,jpk
1097                    IF (umask(ji,jj,jk) == 0) EXIT
1098                    N_out = N_out + 1
1099                    h_out(N_out) = e3u(ji,jj,jk,Krhs_a)
1100                 END DO
1101
1102                 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj)
1103                 DO jk=2,N_out
1104                    z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1) + h_out(jk)) 
1105                 END DO 
1106
1107                 IF (N_in*N_out > 0) THEN
1108                     IF( l_ini_child ) THEN
1109                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
1110                     ELSE
1111                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
1112                     ENDIF   
1113                 ENDIF
1114               END DO
1115            END DO
1116         ELSE
1117            DO jk = 1, jpkm1
1118               DO jj=j1,j2
1119                  uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
1120               END DO
1121            END DO
1122         ENDIF
1123
1124      ENDIF
1125      !
1126   END SUBROUTINE interpun
1127
1128   
1129   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1130      !!----------------------------------------------------------------------
1131      !!                  *** ROUTINE interpvn ***
1132      !!----------------------------------------------------------------------
1133      !
1134      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
1135      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1136      LOGICAL, INTENT(in) :: before
1137      !
1138      INTEGER :: ji,jj,jk
1139      REAL(wp) :: zrhox
1140      ! vertical interpolation:
1141      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
1142      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
1143      INTEGER  :: N_in, N_out, item
1144      REAL(wp) :: h_diff, zhtot
1145      !!---------------------------------------------   
1146      !     
1147      IF (before) THEN   
1148
1149         item = Kmm_a
1150         IF( l_ini_child )   Kmm_a = Kbb_a     
1151       
1152         DO jk=k1,k2
1153            DO jj=j1,j2
1154               DO ji=i1,i2
1155                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
1156                  IF( l_vremap .OR. l_ini_child) THEN
1157                     ! Interpolate thicknesses (masked for subsequent extrapolation)
1158                     ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
1159                  ENDIF
1160               END DO
1161            END DO
1162         END DO
1163
1164         IF( l_vremap .OR. l_ini_child) THEN
1165         ! Extrapolate thicknesses in partial bottom cells:
1166         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1167            IF (ln_zps) THEN
1168               DO jj=j1,j2
1169                  DO ji=i1,i2
1170                     jk = mbkv(ji,jj)
1171                     ptab(ji,jj,jk,2) = 0._wp
1172                  END DO
1173               END DO           
1174            END IF
1175            ! Save ssh at last level:
1176            ptab(i1:i2,j1:j2,k2,2) = 0._wp
1177            IF (.NOT.ln_linssh) THEN
1178               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
1179               DO jk=1,jpk
1180                  ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk)
1181               END DO
1182               ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
1183            END IF
1184         ENDIF
1185         item = Kmm_a
1186
1187      ELSE       
1188         zrhox = Agrif_rhox()
1189
1190         IF( l_vremap .OR. l_ini_child ) THEN
1191
1192            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1193
1194            DO jj=j1,j2
1195               DO ji=i1,i2
1196                  vv(ji,jj,:,Krhs_a) = 0._wp
1197                  N_in = mbkv_parent(ji,jj)
1198                  zhtot = 0._wp
1199                  DO jk=1,N_in
1200                     !IF (jk==N_in) THEN
1201                     !   h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1202                     !ELSE
1203                     !   h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1204                     !ENDIF
1205                     h_in(jk) = e3v0_parent(ji,jj,jk)
1206                     zhtot = zhtot + h_in(jk)
1207                     IF( h_in(jk) .GT. 0. ) THEN
1208                       tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
1209                     ELSE
1210                       tabin(jk)  = 0.
1211                     ENDIF
1212                  END DO
1213
1214                  z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj)
1215                  DO jk=2,N_in
1216                     z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk-1)+h_in(jk))
1217                  END DO
1218
1219                  N_out = 0
1220                  DO jk=1,jpk
1221                     IF (vmask(ji,jj,jk) == 0) EXIT
1222                     N_out = N_out + 1
1223                     h_out(N_out) = e3v(ji,jj,jk,Krhs_a)
1224                  END DO
1225
1226                  z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj)
1227                  DO jk=2,N_out
1228                     z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1)+h_out(jk))
1229                  END DO
1230 
1231                  IF (N_in*N_out > 0) THEN
1232                     IF( l_ini_child ) THEN
1233                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
1234                     ELSE
1235                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
1236                     ENDIF   
1237                  ENDIF
1238               END DO
1239            END DO
1240         ELSE
1241            DO jk = 1, jpkm1
1242               vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
1243            END DO
1244         ENDIF
1245      ENDIF
1246      !       
1247   END SUBROUTINE interpvn
1248
1249   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
1250      !!----------------------------------------------------------------------
1251      !!                  ***  ROUTINE interpunb  ***
1252      !!---------------------------------------------------------------------- 
1253      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1254      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1255      LOGICAL                         , INTENT(in   ) ::   before
1256      !
1257      INTEGER  ::   ji, jj
1258      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1259      !!---------------------------------------------------------------------- 
1260      !
1261      IF( before ) THEN
1262         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu(i1:i2,j1:j2,Kmm_a) * uu_b(i1:i2,j1:j2,Kmm_a)
1263      ELSE
1264         zrhoy = Agrif_Rhoy()
1265         zrhot = Agrif_rhot()
1266         ! Time indexes bounds for integration
1267         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1268         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1269         !
1270         DO ji = i1, i2
1271            DO jj = j1, j2
1272               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1273                  IF    ( utint_stage(ji,jj) == 1  ) THEN
1274                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1275                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1276                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
1277                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1278                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1279                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
1280                     ztcoeff = 1._wp
1281                  ELSE
1282                     ztcoeff = 0._wp
1283                  ENDIF
1284                  !   
1285                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
1286                  !           
1287                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
1288                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
1289                  ENDIF
1290                  !
1291                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
1292               ENDIF
1293            END DO
1294         END DO
1295      END IF
1296      !
1297   END SUBROUTINE interpunb
1298
1299
1300   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
1301      !!----------------------------------------------------------------------
1302      !!                  ***  ROUTINE interpvnb  ***
1303      !!---------------------------------------------------------------------- 
1304      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1305      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1306      LOGICAL                         , INTENT(in   ) ::   before
1307      !
1308      INTEGER  ::   ji, jj
1309      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1310      !!---------------------------------------------------------------------- 
1311      !
1312      IF( before ) THEN
1313         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv(i1:i2,j1:j2,Kmm_a) * vv_b(i1:i2,j1:j2,Kmm_a)
1314      ELSE
1315         zrhox = Agrif_Rhox()
1316         zrhot = Agrif_rhot()
1317         ! Time indexes bounds for integration
1318         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1319         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
1320         !     
1321         DO ji = i1, i2
1322            DO jj = j1, j2
1323               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1324                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
1325                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1326                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1327                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
1328                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1329                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1330                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
1331                     ztcoeff = 1._wp
1332                  ELSE
1333                     ztcoeff = 0._wp
1334                  ENDIF
1335                  !   
1336                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
1337                  !           
1338                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
1339                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
1340                  ENDIF
1341                  !
1342                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
1343               ENDIF
1344            END DO
1345         END DO         
1346      ENDIF
1347      !
1348   END SUBROUTINE interpvnb
1349
1350
1351   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
1352      !!----------------------------------------------------------------------
1353      !!                  ***  ROUTINE interpub2b  ***
1354      !!---------------------------------------------------------------------- 
1355      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1356      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1357      LOGICAL                         , INTENT(in   ) ::   before
1358      !
1359      INTEGER  ::   ji,jj
1360      REAL(wp) ::   zrhot, zt0, zt1, zat
1361      !!---------------------------------------------------------------------- 
1362      IF( before ) THEN
1363         IF ( ln_bt_fw ) THEN
1364            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1365         ELSE
1366            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1367         ENDIF
1368      ELSE
1369         zrhot = Agrif_rhot()
1370         ! Time indexes bounds for integration
1371         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1372         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1373         ! Polynomial interpolation coefficients:
1374         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1375            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1376         !
1377         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1378         !
1379         ! Update interpolation stage:
1380         utint_stage(i1:i2,j1:j2) = 1
1381      ENDIF
1382      !
1383   END SUBROUTINE interpub2b
1384   
1385   SUBROUTINE interpub2b_const( ptab, i1, i2, j1, j2, before )
1386      !!----------------------------------------------------------------------
1387      !!                  ***  ROUTINE interpub2b_const  ***
1388      !!---------------------------------------------------------------------- 
1389      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1390      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1391      LOGICAL                         , INTENT(in   ) ::   before
1392      !
1393      REAL(wp) :: zrhoy
1394      !!---------------------------------------------------------------------- 
1395      IF( before ) THEN
1396         IF ( ln_bt_fw ) THEN
1397            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1398         ELSE
1399            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1400         ENDIF
1401      ELSE
1402         zrhoy = Agrif_Rhoy()
1403         !
1404         ubdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) & 
1405                           & / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
1406         !
1407      ENDIF
1408      !
1409   END SUBROUTINE interpub2b_const
1410
1411
1412   SUBROUTINE ub2b_cor( ptab, i1, i2, j1, j2, before )
1413      !!----------------------------------------------------------------------
1414      !!                  ***  ROUTINE ub2b_cor  ***
1415      !!---------------------------------------------------------------------- 
1416      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1417      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1418      LOGICAL                         , INTENT(in   ) ::   before
1419      !
1420      INTEGER  :: ji, jj
1421      REAL(wp) :: zrhox, zrhoy, zx
1422      !!---------------------------------------------------------------------- 
1423      IF( before ) THEN
1424         ptab(:,:) = 0._wp
1425         DO ji=i1+1,i2-1
1426            DO jj=j1+1,j2
1427               ptab(ji,jj) = 0.25_wp*( ( vb2_b(ji+1,jj  )*e1v(ji+1,jj  )   & 
1428                           &            -vb2_b(ji-1,jj  )*e1v(ji-1,jj  ) ) &
1429                           &          -( vb2_b(ji+1,jj-1)*e1v(ji+1,jj-1)   &
1430                           &            -vb2_b(ji-1,jj-1)*e1v(ji-1,jj-1) ) )
1431            END DO
1432         END DO
1433      ELSE
1434         !
1435         zrhox = Agrif_Rhox() 
1436         zrhoy = Agrif_Rhoy()
1437         DO ji=i1,i2
1438            DO jj=j1,j2
1439               IF (utint_stage(ji,jj)==0) THEN
1440                  zx = 2._wp*MOD(ABS(mig0(ji)-nbghostcells-1), INT(Agrif_Rhox()))/zrhox - 1._wp 
1441                  ubdy(ji,jj) = ubdy(ji,jj) + 0.25_wp*(1._wp-zx*zx) * ptab(ji,jj) & 
1442                              &         / zrhoy *r1_e2u(ji,jj) * umask(ji,jj,1) 
1443                  utint_stage(ji,jj) = 1 
1444               ENDIF
1445            END DO
1446         END DO 
1447         !
1448      ENDIF
1449      !
1450   END SUBROUTINE ub2b_cor
1451
1452
1453   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
1454      !!----------------------------------------------------------------------
1455      !!                  ***  ROUTINE interpvb2b  ***
1456      !!---------------------------------------------------------------------- 
1457      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1458      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1459      LOGICAL                         , INTENT(in   ) ::   before
1460      !
1461      INTEGER ::   ji,jj
1462      REAL(wp) ::   zrhot, zt0, zt1, zat
1463      !!---------------------------------------------------------------------- 
1464      !
1465      IF( before ) THEN
1466         IF ( ln_bt_fw ) THEN
1467            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1468         ELSE
1469            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1470         ENDIF
1471      ELSE     
1472         zrhot = Agrif_rhot()
1473         ! Time indexes bounds for integration
1474         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1475         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1476         ! Polynomial interpolation coefficients:
1477         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1478            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1479         !
1480         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1481         !
1482         ! update interpolation stage:
1483         vtint_stage(i1:i2,j1:j2) = 1
1484      ENDIF
1485      !     
1486   END SUBROUTINE interpvb2b
1487
1488
1489   SUBROUTINE interpvb2b_const( ptab, i1, i2, j1, j2, before )
1490      !!----------------------------------------------------------------------
1491      !!                  ***  ROUTINE interpub2b_const  ***
1492      !!---------------------------------------------------------------------- 
1493      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1494      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1495      LOGICAL                         , INTENT(in   ) ::   before
1496      !
1497      REAL(wp) :: zrhox
1498      !!---------------------------------------------------------------------- 
1499      IF( before ) THEN
1500         IF ( ln_bt_fw ) THEN
1501            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1502         ELSE
1503            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1504         ENDIF
1505      ELSE
1506         zrhox = Agrif_Rhox()
1507         !
1508         vbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) &
1509                           & / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
1510         !
1511      ENDIF
1512      !
1513   END SUBROUTINE interpvb2b_const
1514
1515 
1516   SUBROUTINE vb2b_cor( ptab, i1, i2, j1, j2, before )
1517      !!----------------------------------------------------------------------
1518      !!                  ***  ROUTINE vb2b_cor  ***
1519      !!---------------------------------------------------------------------- 
1520      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1521      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1522      LOGICAL                         , INTENT(in   ) ::   before
1523      !
1524      INTEGER  :: ji, jj
1525      REAL(wp) :: zrhox, zrhoy, zy
1526      !!---------------------------------------------------------------------- 
1527      IF( before ) THEN
1528         ptab(:,:) = 0._wp
1529         DO ji=i1+1,i2-1
1530            DO jj=j1+1,j2
1531               ptab(ji,jj) = 0.25_wp*( ( ub2_b(ji  ,jj+1)*e2u(ji  ,jj+1)   & 
1532                           &            -ub2_b(ji  ,jj-1)*e2u(ji  ,jj-1) ) &
1533                           &          -( ub2_b(ji-1,jj+1)*e2u(ji-1,jj+1)   &
1534                           &            -ub2_b(ji-1,jj-1)*e2u(ji-1,jj-1) ) )
1535            END DO
1536         END DO
1537      ELSE
1538         !
1539         zrhox = Agrif_Rhox() 
1540         zrhoy = Agrif_Rhoy()
1541         DO ji=i1,i2
1542            DO jj=j1,j2
1543               IF (vtint_stage(ji,jj)==0) THEN
1544                  zy = 2._wp*MOD(ABS(mjg0(jj)-nbghostcells-1), INT(Agrif_Rhoy()))/zrhoy - 1._wp 
1545                  vbdy(ji,jj) = vbdy(ji,jj) + 0.25_wp*(1._wp-zy*zy) * ptab(ji,jj) & 
1546                              &         / zrhox * r1_e1v(ji,jj) * vmask(ji,jj,1) 
1547                  vtint_stage(ji,jj) = 1 
1548               ENDIF
1549            END DO
1550         END DO 
1551         !
1552      ENDIF
1553      !
1554   END SUBROUTINE vb2b_cor
1555
1556
1557   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before )
1558      !!----------------------------------------------------------------------
1559      !!                  ***  ROUTINE interpe3t  ***
1560      !!---------------------------------------------------------------------- 
1561      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1562      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1563      LOGICAL                              , INTENT(in   ) :: before
1564      !
1565      INTEGER :: ji, jj, jk
1566      !!---------------------------------------------------------------------- 
1567      !   
1568      IF( before ) THEN
1569         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1570      ELSE
1571         !
1572         DO jk = k1, k2
1573            DO jj = j1, j2
1574               DO ji = i1, i2
1575                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1576                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  & 
1577                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), &
1578                     &                 mig0(ji), mjg0(jj), jk
1579                !     kindic_agr = kindic_agr + 1
1580                  ENDIF
1581               END DO
1582            END DO
1583         END DO
1584         !
1585      ENDIF
1586      !
1587   END SUBROUTINE interpe3t
1588
1589
1590   SUBROUTINE interpe3t0_vremap( ptab, i1, i2, j1, j2, k1, k2, before )
1591      !!----------------------------------------------------------------------
1592      !!                  ***  ROUTINE interpe3t0_vremap  ***
1593      !!---------------------------------------------------------------------- 
1594      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1595      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1596      LOGICAL                              , INTENT(in   ) :: before
1597      !
1598      INTEGER :: ji, jj, jk
1599      REAL(wp) :: zh
1600      !!---------------------------------------------------------------------- 
1601      !   
1602      IF( before ) THEN
1603         IF ( ln_zps ) THEN
1604            DO jk = k1, k2
1605               DO jj = j1, j2
1606                  DO ji = i1, i2
1607                     ptab(ji, jj, jk) = e3t_1d(jk)
1608                  END DO
1609               END DO
1610            END DO
1611         ELSE
1612            ptab(i1:i2,j1:j2,k1:k2) = e3t_0(i1:i2,j1:j2,k1:k2)
1613         ENDIF
1614      ELSE
1615         !
1616         DO jk = k1, k2
1617            DO jj = j1, j2
1618               DO ji = i1, i2
1619                  e3t0_parent(ji,jj,jk) = ptab(ji,jj,jk)
1620               END DO
1621            END DO
1622         END DO
1623
1624         ! Retrieve correct scale factor at the bottom:
1625         DO jj = j1, j2
1626            DO ji = i1, i2
1627               zh = 0._wp
1628               DO jk = 1, mbkt_parent(ji, jj)-1
1629                  zh = zh + e3t0_parent(ji,jj,jk)
1630               END DO
1631               e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh
1632            END DO
1633         END DO
1634         
1635      ENDIF
1636      !
1637   END SUBROUTINE interpe3t0_vremap
1638
1639
1640   SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before )
1641      !!----------------------------------------------------------------------
1642      !!                  ***  ROUTINE interpglamt  ***
1643      !!---------------------------------------------------------------------- 
1644      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
1645      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1646      LOGICAL                        , INTENT(in   ) :: before
1647      !
1648      INTEGER :: ji, jj, jk
1649      REAL(wp):: ztst
1650      !!---------------------------------------------------------------------- 
1651      !   
1652      IF( before ) THEN
1653         ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2)
1654      ELSE
1655         ztst = MAXVAL(ABS(glamt(i1:i2,j1:j2)))*1.e-4
1656         DO jj = j1, j2
1657            DO ji = i1, i2
1658               IF( ABS( ptab(ji,jj) - glamt(ji,jj) ) > ztst ) THEN
1659                  WRITE(numout,*) ' Agrif error for glamt: parent, child, i, j ', ptab(ji,jj), glamt(ji,jj), mig0(ji), mig0(jj)
1660!                  kindic_agr = kindic_agr + 1
1661               ENDIF
1662            END DO
1663         END DO
1664      ENDIF
1665      !
1666   END SUBROUTINE interpglamt
1667
1668
1669   SUBROUTINE interpgphit( ptab, i1, i2, j1, j2, before )
1670      !!----------------------------------------------------------------------
1671      !!                  ***  ROUTINE interpgphit  ***
1672      !!---------------------------------------------------------------------- 
1673      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
1674      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1675      LOGICAL                        , INTENT(in   ) :: before
1676      !
1677      INTEGER :: ji, jj, jk
1678      REAL(wp):: ztst
1679      !!---------------------------------------------------------------------- 
1680      !   
1681      IF( before ) THEN
1682         ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2)
1683      ELSE
1684         ztst = MAXVAL(ABS(gphit(i1:i2,j1:j2)))*1.e-4
1685         DO jj = j1, j2
1686            DO ji = i1, i2
1687               IF( ABS( ptab(ji,jj) - gphit(ji,jj) ) > ztst ) THEN
1688                  WRITE(numout,*) ' Agrif error for gphit: parent, child, i, j ', ptab(ji,jj), gphit(ji,jj), mig0(ji), mig0(jj)
1689!                  kindic_agr = kindic_agr + 1
1690               ENDIF
1691            END DO
1692         END DO
1693      ENDIF
1694      !
1695   END SUBROUTINE interpgphit
1696
1697
1698   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1699      !!----------------------------------------------------------------------
1700      !!                  ***  ROUTINE interavm  ***
1701      !!---------------------------------------------------------------------- 
1702      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1703      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1704      LOGICAL                                    , INTENT(in   ) ::   before
1705      !
1706      INTEGER  :: ji, jj, jk
1707      INTEGER  :: N_in, N_out
1708      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in
1709      REAL(wp), DIMENSION(1:jpk) :: z_out
1710      !!---------------------------------------------------------------------- 
1711      !     
1712      IF (before) THEN         
1713         DO jk=k1,k2
1714            DO jj=j1,j2
1715              DO ji=i1,i2
1716                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1717              END DO
1718           END DO
1719         END DO
1720
1721         IF( l_vremap ) THEN
1722            ! Interpolate thicknesses
1723            ! Warning: these are masked, hence extrapolated prior interpolation.
1724            DO jk=k1,k2
1725               DO jj=j1,j2
1726                  DO ji=i1,i2
1727                      ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
1728                  END DO
1729               END DO
1730            END DO
1731
1732            ! Extrapolate thicknesses in partial bottom cells:
1733            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1734            IF (ln_zps) THEN
1735               DO jj=j1,j2
1736                  DO ji=i1,i2
1737                      jk = mbkt(ji,jj)
1738                      ptab(ji,jj,jk,2) = 0._wp
1739                  END DO
1740               END DO           
1741            END IF
1742       
1743           ! Save ssh at last level:
1744            IF (.NOT.ln_linssh) THEN
1745               ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
1746            ELSE
1747               ptab(i1:i2,j1:j2,k2,2) = 0._wp
1748            END IF     
1749          ENDIF
1750
1751      ELSE
1752
1753         IF( l_vremap ) THEN
1754            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1755            avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
1756               
1757            DO jj = j1, j2
1758               DO ji =i1, i2
1759                  N_in = mbkt_parent(ji,jj)
1760                  IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
1761                  z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
1762                  DO jk = N_in, 1, -1  ! Parent vertical grid               
1763                        z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
1764                       tabin(jk) = ptab(ji,jj,jk,1)
1765                  END DO
1766                  N_out = mbkt(ji,jj) 
1767                  DO jk = 1, N_out        ! Child vertical grid
1768                     z_out(jk) = gdepw(ji,jj,jk,Kmm_a)
1769                  END DO
1770                  IF (N_in*N_out > 0) THEN
1771                     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)
1772                  ENDIF
1773               END DO
1774            END DO
1775         ELSE
1776            avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1777         ENDIF
1778      ENDIF
1779      !
1780   END SUBROUTINE interpavm
1781
1782   
1783   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
1784      !!----------------------------------------------------------------------
1785      !!                  ***  ROUTINE interpmbkt  ***
1786      !!---------------------------------------------------------------------- 
1787      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1788      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1789      LOGICAL                         , INTENT(in   ) ::   before
1790      !
1791      !!---------------------------------------------------------------------- 
1792      !
1793      IF( before) THEN
1794         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
1795      ELSE
1796         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
1797      ENDIF
1798      !
1799   END SUBROUTINE interpmbkt
1800
1801   
1802   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
1803      !!----------------------------------------------------------------------
1804      !!                  ***  ROUTINE interpht0  ***
1805      !!---------------------------------------------------------------------- 
1806      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1807      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1808      LOGICAL                         , INTENT(in   ) ::   before
1809      !
1810      !!---------------------------------------------------------------------- 
1811      !
1812      IF( before) THEN
1813         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2)
1814      ELSE
1815         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
1816      ENDIF
1817      !
1818   END SUBROUTINE interpht0
1819
1820   SUBROUTINE Agrif_check_bat( iindic )
1821      !!----------------------------------------------------------------------
1822      !!                  ***  ROUTINE Agrif_check_bat  ***
1823      !!---------------------------------------------------------------------- 
1824      INTEGER, INTENT(inout) ::   iindic
1825      !!
1826      INTEGER :: ji, jj
1827      INTEGER  :: istart, iend, jstart, jend, ispon
1828      !!---------------------------------------------------------------------- 
1829      !
1830      !
1831      ! --- West --- !
1832      IF(lk_west) THEN
1833         ispon  = nn_sponge_len * Agrif_irhox()
1834         istart = nn_hls + 2                                  ! halo + land + 1
1835         iend   = nn_hls + 1 + nbghostcells + ispon           ! halo + land + nbghostcells + sponge
1836         jstart = nn_hls + 2
1837         jend   = jpjglo - nn_hls - 1
1838         DO ji = mi0(istart), mi1(iend)
1839            DO jj = mj0(jstart), mj1(jend)
1840               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1841            END DO
1842            DO jj = mj0(jstart), mj1(jend-1)
1843               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1844            END DO
1845         END DO
1846         DO ji = mi0(istart), mi1(iend-1)
1847            DO jj = mj0(jstart), mj1(jend)
1848               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1849            END DO
1850         END DO
1851      ENDIF
1852      !
1853      ! --- East --- !
1854      IF(lk_east) THEN
1855         ispon  = nn_sponge_len * Agrif_irhox() 
1856         istart = jpiglo - ( nn_hls + nbghostcells + ispon )  ! halo + land + nbghostcells + sponge - 1
1857         iend   = jpiglo - ( nn_hls + 1 )                     ! halo + land + 1                     - 1
1858         jstart = nn_hls + 2
1859         jend   = jpjglo - nn_hls - 1 
1860         DO ji = mi0(istart), mi1(iend)
1861            DO jj = mj0(jstart), mj1(jend)
1862               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1863            END DO
1864            DO jj = mj0(jstart), mj1(jend-1)
1865               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1866            END DO
1867         END DO
1868         DO ji = mi0(istart+1), mi1(iend-1)
1869            DO jj = mj0(jstart), mj1(jend)
1870               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1871            END DO
1872         END DO
1873      ENDIF
1874      !
1875      ! --- South --- !
1876      IF(lk_south) THEN
1877         ispon  = nn_sponge_len * Agrif_irhoy() 
1878         jstart = nn_hls + 2                                 ! halo + land + 1
1879         jend   = nn_hls + 1 + nbghostcells + ispon          ! halo + land + nbghostcells + sponge
1880         istart = nn_hls + 2
1881         iend   = jpiglo - nn_hls - 1
1882         DO jj = mj0(jstart), mj1(jend)
1883            DO ji = mi0(istart), mi1(iend)
1884               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1885            END DO
1886            DO ji = mi0(istart), mi1(iend-1)
1887               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1888            END DO
1889         END DO
1890         DO jj = mj0(jstart), mj1(jend-1)
1891            DO ji = mi0(istart), mi1(iend)
1892               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1893            END DO
1894         END DO
1895      ENDIF
1896      !
1897      ! --- North --- !
1898      IF(lk_north) THEN
1899         ispon  = nn_sponge_len * Agrif_irhoy() 
1900         jstart = jpjglo - ( nn_hls + nbghostcells + ispon)  ! halo + land + nbghostcells +sponge - 1
1901         jend   = jpjglo - ( nn_hls + 1 )                    ! halo + land + 1            - 1
1902         istart = nn_hls + 2
1903         iend   = jpiglo - nn_hls - 1
1904         DO jj = mj0(jstart), mj1(jend)
1905            DO ji = mi0(istart), mi1(iend)
1906               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1907            END DO
1908            DO ji = mi0(istart), mi1(iend-1)
1909               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1910            END DO
1911         END DO
1912         DO jj = mj0(jstart+1), mj1(jend-1)
1913            DO ji = mi0(istart), mi1(iend)
1914               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1915            END DO
1916         END DO
1917      ENDIF
1918      !
1919   END SUBROUTINE Agrif_check_bat
1920   
1921#else
1922   !!----------------------------------------------------------------------
1923   !!   Empty module                                          no AGRIF zoom
1924   !!----------------------------------------------------------------------
1925CONTAINS
1926   SUBROUTINE Agrif_OCE_Interp_empty
1927      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1928   END SUBROUTINE Agrif_OCE_Interp_empty
1929#endif
1930
1931   !!======================================================================
1932END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.