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

Last change on this file since 13286 was 13286, checked in by smasson, 2 months ago

trunk: merge extra halos branch in trunk, see #2366

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