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

Last change on this file since 12377 was 12377, checked in by acc, 8 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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