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 @ 14088

Last change on this file since 14088 was 14086, checked in by cetlod, 3 years ago

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