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/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_oce_interp.F90 @ 14170

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

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