source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_interp.F90 @ 13130

Last change on this file since 13130 was 13130, checked in by smasson, 3 months ago

Extra_Halo: supress halos from outputs and coupling, see #2366

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