New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_oce_interp.F90 in NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST – NEMO

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

Last change on this file since 13230 was 13230, checked in by smasson, 4 years ago

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