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_r12973_AGRIF_CMEMS/src/NST – NEMO

source: NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/NST/agrif_oce_interp.F90 @ 13026

Last change on this file since 13026 was 13026, checked in by rblod, 4 years ago

AGRIF with northfold and perio, see ticket #2129
  • Property svn:keywords set to Id
File size: 51.2 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
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 = 2
103         ibdy2 = 1+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-1-nbghostcells
163         ibdy2 = jpiglo-2 
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 = 2
227         jbdy2 = 1+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-1-nbghostcells
290         jbdy2 = jpjglo-2 
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 = 2
370         iend   = nbghostcells+1
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-nbghostcells
382         iend   = jpiglo-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-nbghostcells-1
390         iend   = jpiglo-2
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 = 2
401         jend   = nbghostcells+1
402         DO jj = mj0(jstart), mj1(jend)
403
404            DO ji=1,jpi
405               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
406               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
407            END DO
408         END DO
409      ENDIF       
410      !
411      !--- North ---!
412      IF( lk_north ) THEN
413         jstart = jpjglo-nbghostcells
414         jend   = jpjglo-1
415         DO jj = mj0(jstart), mj1(jend)
416            DO ji=1,jpi
417               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
418            END DO
419         END DO
420         jstart = jpjglo-nbghostcells-1
421         jend   = jpjglo-2
422         DO jj = mj0(jstart), mj1(jend)
423            DO ji=1,jpi
424               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
425            END DO
426         END DO
427      ENDIF 
428      !
429   END SUBROUTINE Agrif_dyn_ts
430
431   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
432      !!----------------------------------------------------------------------
433      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
434      !!---------------------------------------------------------------------- 
435      INTEGER, INTENT(in) ::   jn
436      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
437      !!
438      INTEGER :: ji, jj
439      INTEGER :: istart, iend, jstart, jend
440      !!---------------------------------------------------------------------- 
441      !
442      IF( Agrif_Root() )   RETURN
443      !
444      !--- West ---!
445      IF( lk_west ) THEN
446         istart = 2
447         iend   = nbghostcells+1
448         DO ji = mi0(istart), mi1(iend)
449            DO jj=1,jpj
450               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
451               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
452            END DO
453         END DO
454      ENDIF
455      !
456      !--- East ---!
457      IF( lk_east ) THEN
458         istart = jpiglo-nbghostcells
459         iend   = jpiglo-1
460         DO ji = mi0(istart), mi1(iend)
461            DO jj=1,jpj
462               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
463            END DO
464         END DO
465         istart = jpiglo-nbghostcells-1
466         iend   = jpiglo-2
467         DO ji = mi0(istart), mi1(iend)
468            DO jj=1,jpj
469               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
470            END DO
471         END DO
472      ENDIF
473      !
474      !--- South ---!
475      IF( lk_south ) THEN
476         jstart = 2
477         jend   = nbghostcells+1
478         DO jj = mj0(jstart), mj1(jend)
479            DO ji=1,jpi
480               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
481               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
482            END DO
483         END DO
484      ENDIF
485      !
486      !--- North ---!
487      IF( lk_north ) THEN
488         jstart = jpjglo-nbghostcells
489         jend   = jpjglo-1
490         DO jj = mj0(jstart), mj1(jend)
491            DO ji=1,jpi
492               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
493            END DO
494         END DO
495         jstart = jpjglo-nbghostcells-1
496         jend   = jpjglo-2
497         DO jj = mj0(jstart), mj1(jend)
498            DO ji=1,jpi
499               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
500            END DO
501         END DO
502      ENDIF
503      !
504   END SUBROUTINE Agrif_dyn_ts_flux
505
506   SUBROUTINE Agrif_dta_ts( kt )
507      !!----------------------------------------------------------------------
508      !!                  ***  ROUTINE Agrif_dta_ts  ***
509      !!---------------------------------------------------------------------- 
510      INTEGER, INTENT(in) ::   kt
511      !!
512      INTEGER :: ji, jj
513      LOGICAL :: ll_int_cons
514      !!---------------------------------------------------------------------- 
515      !
516      IF( Agrif_Root() )   RETURN
517      !
518      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
519      !
520      ! Enforce volume conservation if no time refinement: 
521      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
522      !
523      ! Interpolate barotropic fluxes
524      Agrif_SpecialValue = 0._wp
525      Agrif_UseSpecialValue = ln_spc_dyn
526
527      use_sign_north = .TRUE.
528      sign_north = -1.
529
530      !
531      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
532      utint_stage(:,:) = 0
533      vtint_stage(:,:) = 0
534      !
535      IF( ll_int_cons ) THEN  ! Conservative interpolation
536         ! order matters here !!!!!!
537         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
538         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
539         !
540         bdy_tinterp = 1
541         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
542         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  ) 
543         !
544         bdy_tinterp = 2
545         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
546         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )   
547      ELSE ! Linear interpolation
548         !
549         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
550         CALL Agrif_Bc_variable( unb_id, procname=interpunb )
551         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
552      ENDIF
553      Agrif_UseSpecialValue = .FALSE.
554      use_sign_north = .FALSE.
555      !
556   END SUBROUTINE Agrif_dta_ts
557
558
559   SUBROUTINE Agrif_ssh( kt )
560      !!----------------------------------------------------------------------
561      !!                  ***  ROUTINE Agrif_ssh  ***
562      !!---------------------------------------------------------------------- 
563      INTEGER, INTENT(in) ::   kt
564      !
565      INTEGER  :: ji, jj
566      INTEGER  :: istart, iend, jstart, jend
567      !!---------------------------------------------------------------------- 
568      !
569      IF( Agrif_Root() )   RETURN
570      !     
571      ! Linear time interpolation of sea level
572      !
573      Agrif_SpecialValue    = 0._wp
574      Agrif_UseSpecialValue = .TRUE.
575      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
576      Agrif_UseSpecialValue = .FALSE.
577      !
578      ! --- West --- !
579      IF(lk_west) THEN
580         istart = 2
581         iend   = 1 + nbghostcells
582         DO ji = mi0(istart), mi1(iend)
583            DO jj = 1, jpj
584               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
585            ENDDO
586         ENDDO
587      ENDIF
588      !
589      ! --- East --- !
590      IF(lk_east) THEN
591         istart = jpiglo - nbghostcells
592         iend   = jpiglo - 1
593         DO ji = mi0(istart), mi1(iend)
594            DO jj = 1, jpj
595               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
596            ENDDO
597         ENDDO
598      ENDIF
599      !
600      ! --- South --- !
601      IF(lk_south) THEN
602         jstart = 2
603         jend   = 1 + nbghostcells
604         DO jj = mj0(jstart), mj1(jend)
605            DO ji = 1, jpi
606               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
607            ENDDO
608         ENDDO
609      ENDIF
610      !
611      ! --- North --- !
612      IF(lk_north) THEN
613         jstart = jpjglo - nbghostcells
614         jend   = jpjglo - 1
615         DO jj = mj0(jstart), mj1(jend)
616            DO ji = 1, jpi
617               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
618            ENDDO
619         ENDDO
620      ENDIF
621      !
622   END SUBROUTINE Agrif_ssh
623
624
625   SUBROUTINE Agrif_ssh_ts( jn )
626      !!----------------------------------------------------------------------
627      !!                  ***  ROUTINE Agrif_ssh_ts  ***
628      !!---------------------------------------------------------------------- 
629      INTEGER, INTENT(in) ::   jn
630      !!
631      INTEGER :: ji, jj
632      INTEGER  :: istart, iend, jstart, jend
633      !!---------------------------------------------------------------------- 
634      !
635      IF( Agrif_Root() )   RETURN
636      !
637      ! --- West --- !
638      IF(lk_west) THEN
639         istart = 2
640         iend   = 1+nbghostcells
641         DO ji = mi0(istart), mi1(iend)
642            DO jj = 1, jpj
643               ssha_e(ji,jj) = hbdy(ji,jj)
644            ENDDO
645         ENDDO
646      ENDIF
647      !
648      ! --- East --- !
649      IF(lk_east) THEN
650         istart = jpiglo - nbghostcells
651         iend   = jpiglo - 1
652         DO ji = mi0(istart), mi1(iend)
653            DO jj = 1, jpj
654               ssha_e(ji,jj) = hbdy(ji,jj)
655            ENDDO
656         ENDDO
657      ENDIF
658      !
659      ! --- South --- !
660      IF(lk_south) THEN
661         jstart = 2
662         jend   = 1+nbghostcells
663         DO jj = mj0(jstart), mj1(jend)
664            DO ji = 1, jpi
665               ssha_e(ji,jj) = hbdy(ji,jj)
666            ENDDO
667         ENDDO
668      ENDIF
669      !
670      ! --- North --- !
671      IF(lk_north) THEN
672         jstart = jpjglo - nbghostcells
673         jend   = jpjglo - 1
674         DO jj = mj0(jstart), mj1(jend)
675            DO ji = 1, jpi
676               ssha_e(ji,jj) = hbdy(ji,jj)
677            ENDDO
678         ENDDO
679      ENDIF
680      !
681   END SUBROUTINE Agrif_ssh_ts
682
683   SUBROUTINE Agrif_avm
684      !!----------------------------------------------------------------------
685      !!                  ***  ROUTINE Agrif_avm  ***
686      !!---------------------------------------------------------------------- 
687      REAL(wp) ::   zalpha
688      !!---------------------------------------------------------------------- 
689      !
690      IF( Agrif_Root() )   RETURN
691      !
692      zalpha = 1._wp ! JC: proper time interpolation impossible 
693                     ! => use last available value from parent
694      !
695      Agrif_SpecialValue    = 0.e0
696      Agrif_UseSpecialValue = .TRUE.
697      !
698      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
699      !
700      Agrif_UseSpecialValue = .FALSE.
701      !
702   END SUBROUTINE Agrif_avm
703   
704
705   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
706      !!----------------------------------------------------------------------
707      !!                  *** ROUTINE interptsn ***
708      !!----------------------------------------------------------------------
709      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
710      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
711      LOGICAL                                     , INTENT(in   ) ::   before
712      !
713      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
714      INTEGER  ::   N_in, N_out
715      INTEGER  :: item
716      ! vertical interpolation:
717      REAL(wp) :: zhtot
718      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin
719      REAL(wp), DIMENSION(k1:k2) :: h_in
720      REAL(wp), DIMENSION(1:jpk) :: h_out
721      !!----------------------------------------------------------------------
722
723      IF( before ) THEN
724
725         item = Kmm_a
726         IF( l_ini_child )   Kmm_a = Kbb_a 
727
728         DO jn = 1,jpts
729            DO jk=k1,k2
730               DO jj=j1,j2
731                 DO ji=i1,i2
732                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)
733                 END DO
734              END DO
735           END DO
736         END DO
737
738         IF( l_vremap .OR. l_ini_child) THEN
739            ! Interpolate thicknesses
740            ! Warning: these are masked, hence extrapolated prior interpolation.
741            DO jk=k1,k2
742               DO jj=j1,j2
743                  DO ji=i1,i2
744                      ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
745
746                  END DO
747               END DO
748            END DO
749
750            ! Extrapolate thicknesses in partial bottom cells:
751            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
752            IF (ln_zps) THEN
753               DO jj=j1,j2
754                  DO ji=i1,i2
755                      jk = mbkt(ji,jj)
756                      ptab(ji,jj,jk,jpts+1) = 0._wp
757                  END DO
758               END DO           
759            END IF
760       
761            ! Save ssh at last level:
762            IF (.NOT.ln_linssh) THEN
763               ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
764            ELSE
765               ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
766            END IF     
767         ENDIF
768         Kmm_a = item
769
770      ELSE
771         item = Krhs_a
772         IF( l_ini_child )   Krhs_a = Kbb_a 
773
774         IF( l_vremap .OR. l_ini_child ) THEN
775            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
776               
777            DO jj=j1,j2
778               DO ji=i1,i2
779                  ts(ji,jj,:,:,Krhs_a) = 0.                 
780               !   IF( l_ini_child) ts(ji,jj,:,:,Krhs_a) = ptab(ji,jj,:,1:jpts)
781                  N_in = mbkt_parent(ji,jj)
782                  zhtot = 0._wp
783                  DO jk=1,N_in !k2 = jpk of parent grid
784                     IF (jk==N_in) THEN
785                        h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
786                     ELSE
787                        h_in(jk) = ptab(ji,jj,jk,n2)
788                     ENDIF
789                     zhtot = zhtot + h_in(jk)
790                     tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
791                  END DO
792                  N_out = 0
793                  DO jk=1,jpk ! jpk of child grid
794                     IF (tmask(ji,jj,jk) == 0._wp) EXIT
795                     N_out = N_out + 1
796                     h_out(jk) = e3t(ji,jj,jk,Krhs_a)
797                  ENDDO
798                  IF (N_in*N_out > 0) THEN
799                     IF( l_ini_child ) THEN
800                        CALL remap_linear(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),          &
801                                      &   h_out(1:N_out),N_in,N_out,jpts) 
802                     ELSE
803                        CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   &
804                                      &   h_out(1:N_out),N_in,N_out,jpts) 
805                     ENDIF
806                  ENDIF
807               ENDDO
808            ENDDO
809            Krhs_a = item
810 
811         ELSE
812         
813            DO jn=1, jpts
814                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) 
815            END DO
816         ENDIF
817
818      ENDIF
819      !
820   END SUBROUTINE interptsn
821
822   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
823      !!----------------------------------------------------------------------
824      !!                  ***  ROUTINE interpsshn  ***
825      !!---------------------------------------------------------------------- 
826      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
827      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
828      LOGICAL                         , INTENT(in   ) ::   before
829      !
830      !!---------------------------------------------------------------------- 
831      !
832      IF( before) THEN
833         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
834      ELSE
835         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
836      ENDIF
837      !
838   END SUBROUTINE interpsshn
839
840   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
841      !!----------------------------------------------------------------------
842      !!                  *** ROUTINE interpun ***
843      !!---------------------------------------------   
844      !!
845      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
846      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
847      LOGICAL, INTENT(in) :: before
848      !!
849      INTEGER :: ji,jj,jk
850      REAL(wp) :: zrhoy, zhtot
851      ! vertical interpolation:
852      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
853      REAL(wp), DIMENSION(1:jpk) :: h_out
854      INTEGER  :: N_in, N_out,item
855      REAL(wp) :: h_diff
856      !!---------------------------------------------   
857      !
858      IF (before) THEN
859
860         item = Kmm_a
861         IF( l_ini_child )   Kmm_a = Kbb_a     
862
863         DO jk=1,jpk
864            DO jj=j1,j2
865               DO ji=i1,i2
866                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
867                  IF( l_vremap .OR. l_ini_child) THEN
868                     ! Interpolate thicknesses (masked for subsequent extrapolation)
869                     ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
870                  ENDIF
871               END DO
872            END DO
873         END DO
874
875        IF( l_vremap .OR. l_ini_child) THEN
876         ! Extrapolate thicknesses in partial bottom cells:
877         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
878            IF (ln_zps) THEN
879               DO jj=j1,j2
880                  DO ji=i1,i2
881                     jk = mbku(ji,jj)
882                     ptab(ji,jj,jk,2) = 0._wp
883                  END DO
884               END DO           
885            END IF
886
887           ! Save ssh at last level:
888           ptab(i1:i2,j1:j2,k2,2) = 0._wp
889           IF (.NOT.ln_linssh) THEN
890              ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
891              DO jk=1,jpk
892                 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)
893              END DO
894              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
895           END IF
896        ENDIF
897
898         Kmm_a = item
899         !
900      ELSE
901         zrhoy = Agrif_rhoy()
902
903        IF( l_vremap .OR. l_ini_child) THEN
904! VERTICAL REFINEMENT BEGIN
905
906            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
907
908            DO ji=i1,i2
909               DO jj=j1,j2
910                  uu(ji,jj,:,Krhs_a) = 0._wp
911                  N_in = mbku_parent(ji,jj)
912                  zhtot = 0._wp
913                  DO jk=1,N_in
914                     IF (jk==N_in) THEN
915                        h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
916                     ELSE
917                        h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
918                     ENDIF
919                     zhtot = zhtot + h_in(jk)
920                     IF( h_in(jk) .GT. 0. ) THEN
921                     tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
922                     ELSE
923                     tabin(jk) = 0.
924                     ENDIF
925                 ENDDO
926                     
927                 N_out = 0
928                 DO jk=1,jpk
929                    IF (umask(ji,jj,jk) == 0) EXIT
930                    N_out = N_out + 1
931                    h_out(N_out) = e3u(ji,jj,jk,Krhs_a)
932                 ENDDO
933                 IF (N_in*N_out > 0) THEN
934                     IF( l_ini_child ) THEN
935                        CALL remap_linear       (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)
936                     ELSE
937                        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)
938                     ENDIF   
939                 ENDIF
940               ENDDO
941            ENDDO
942         ELSE
943            DO jk = 1, jpkm1
944               DO jj=j1,j2
945                  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) )
946               END DO
947            END DO
948         ENDIF
949
950      ENDIF
951      !
952   END SUBROUTINE interpun
953
954   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
955      !!----------------------------------------------------------------------
956      !!                  *** ROUTINE interpvn ***
957      !!----------------------------------------------------------------------
958      !
959      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
960      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
961      LOGICAL, INTENT(in) :: before
962      !
963      INTEGER :: ji,jj,jk
964      REAL(wp) :: zrhox
965      ! vertical interpolation:
966      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
967      REAL(wp), DIMENSION(1:jpk) :: h_out
968      INTEGER  :: N_in, N_out, item
969      REAL(wp) :: h_diff, zhtot
970      !!---------------------------------------------   
971      !     
972      IF (before) THEN   
973
974         item = Kmm_a
975         IF( l_ini_child )   Kmm_a = Kbb_a     
976       
977         DO jk=k1,k2
978            DO jj=j1,j2
979               DO ji=i1,i2
980                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
981                  IF( l_vremap .OR. l_ini_child) THEN
982                     ! Interpolate thicknesses (masked for subsequent extrapolation)
983                     ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
984                  ENDIF
985               END DO
986            END DO
987         END DO
988
989         IF( l_vremap .OR. l_ini_child) THEN
990         ! Extrapolate thicknesses in partial bottom cells:
991         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
992            IF (ln_zps) THEN
993               DO jj=j1,j2
994                  DO ji=i1,i2
995                     jk = mbkv(ji,jj)
996                     ptab(ji,jj,jk,2) = 0._wp
997                  END DO
998               END DO           
999            END IF
1000            ! Save ssh at last level:
1001            ptab(i1:i2,j1:j2,k2,2) = 0._wp
1002            IF (.NOT.ln_linssh) THEN
1003               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
1004               DO jk=1,jpk
1005                  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)
1006               END DO
1007               ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
1008            END IF
1009         ENDIF
1010         item = Kmm_a
1011
1012      ELSE       
1013         zrhox = Agrif_rhox()
1014
1015         IF( l_vremap .OR. l_ini_child ) THEN
1016
1017            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1018
1019            DO jj=j1,j2
1020               DO ji=i1,i2
1021                  vv(ji,jj,:,Krhs_a) = 0._wp
1022                  N_in = mbkv_parent(ji,jj)
1023                  zhtot = 0._wp
1024                  DO jk=1,N_in
1025                     IF (jk==N_in) THEN
1026                        h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1027                     ELSE
1028                        h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
1029                     ENDIF
1030                     zhtot = zhtot + h_in(jk)
1031                    IF( h_in(jk) .GT. 0. ) THEN
1032                     tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
1033                    ELSE
1034                     tabin(jk)  = 0.
1035                    ENDIF
1036                 ENDDO
1037           
1038                  N_out = 0
1039                  DO jk=1,jpk
1040                     if (vmask(ji,jj,jk) == 0) EXIT
1041                     N_out = N_out + 1
1042                     h_out(N_out) = e3v(ji,jj,jk,Krhs_a)
1043                  END DO
1044                  IF (N_in*N_out > 0) THEN
1045                     IF( l_ini_child ) THEN
1046                        CALL remap_linear       (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)
1047                     ELSE
1048                        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)
1049                     ENDIF   
1050                  ENDIF
1051               END DO
1052            END DO
1053         ELSE
1054            DO jk = 1, jpkm1
1055               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) )
1056            END DO
1057         ENDIF
1058      ENDIF
1059      !       
1060   END SUBROUTINE interpvn
1061
1062   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
1063      !!----------------------------------------------------------------------
1064      !!                  ***  ROUTINE interpunb  ***
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) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1072      !!---------------------------------------------------------------------- 
1073      !
1074      IF( before ) THEN
1075         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)
1076      ELSE
1077         zrhoy = Agrif_Rhoy()
1078         zrhot = Agrif_rhot()
1079         ! Time indexes bounds for integration
1080         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1081         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1082         !
1083         DO ji = i1, i2
1084            DO jj = j1, j2
1085               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1086                  IF    ( utint_stage(ji,jj) == 1  ) THEN
1087                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1088                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1089                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
1090                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1091                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1092                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
1093                     ztcoeff = 1._wp
1094                  ELSE
1095                     ztcoeff = 0._wp
1096                  ENDIF
1097                  !   
1098                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
1099                  !           
1100                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
1101                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
1102                  ENDIF
1103                  !
1104                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
1105               ENDIF
1106            END DO
1107         END DO
1108      END IF
1109      !
1110   END SUBROUTINE interpunb
1111
1112
1113   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
1114      !!----------------------------------------------------------------------
1115      !!                  ***  ROUTINE interpvnb  ***
1116      !!---------------------------------------------------------------------- 
1117      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1118      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1119      LOGICAL                         , INTENT(in   ) ::   before
1120      !
1121      INTEGER  ::   ji, jj
1122      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1123      !!---------------------------------------------------------------------- 
1124      !
1125      IF( before ) THEN
1126         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)
1127      ELSE
1128         zrhox = Agrif_Rhox()
1129         zrhot = Agrif_rhot()
1130         ! Time indexes bounds for integration
1131         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1132         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
1133         !     
1134         DO ji = i1, i2
1135            DO jj = j1, j2
1136               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1137                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
1138                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1139                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1140                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
1141                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1142                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1143                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
1144                     ztcoeff = 1._wp
1145                  ELSE
1146                     ztcoeff = 0._wp
1147                  ENDIF
1148                  !   
1149                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
1150                  !           
1151                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
1152                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
1153                  ENDIF
1154                  !
1155                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
1156               ENDIF
1157            END DO
1158         END DO         
1159      ENDIF
1160      !
1161   END SUBROUTINE interpvnb
1162
1163
1164   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
1165      !!----------------------------------------------------------------------
1166      !!                  ***  ROUTINE interpub2b  ***
1167      !!---------------------------------------------------------------------- 
1168      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1169      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1170      LOGICAL                         , INTENT(in   ) ::   before
1171      !
1172      INTEGER  ::   ji,jj
1173      REAL(wp) ::   zrhot, zt0, zt1, zat
1174      !!---------------------------------------------------------------------- 
1175      IF( before ) THEN
1176         IF ( ln_bt_fw ) THEN
1177            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1178         ELSE
1179            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1180         ENDIF
1181      ELSE
1182         zrhot = Agrif_rhot()
1183         ! Time indexes bounds for integration
1184         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1185         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1186         ! Polynomial interpolation coefficients:
1187         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1188            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1189         !
1190         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1191         !
1192         ! Update interpolation stage:
1193         utint_stage(i1:i2,j1:j2) = 1
1194      ENDIF
1195      !
1196   END SUBROUTINE interpub2b
1197   
1198
1199   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
1200      !!----------------------------------------------------------------------
1201      !!                  ***  ROUTINE interpvb2b  ***
1202      !!---------------------------------------------------------------------- 
1203      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1204      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1205      LOGICAL                         , INTENT(in   ) ::   before
1206      !
1207      INTEGER ::   ji,jj
1208      REAL(wp) ::   zrhot, zt0, zt1, zat
1209      !!---------------------------------------------------------------------- 
1210      !
1211      IF( before ) THEN
1212         IF ( ln_bt_fw ) THEN
1213            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1214         ELSE
1215            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1216         ENDIF
1217      ELSE     
1218         zrhot = Agrif_rhot()
1219         ! Time indexes bounds for integration
1220         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1221         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1222         ! Polynomial interpolation coefficients:
1223         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1224            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1225         !
1226         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1227         !
1228         ! update interpolation stage:
1229         vtint_stage(i1:i2,j1:j2) = 1
1230      ENDIF
1231      !     
1232   END SUBROUTINE interpvb2b
1233
1234
1235   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before )
1236      !!----------------------------------------------------------------------
1237      !!                  ***  ROUTINE interpe3t  ***
1238      !!---------------------------------------------------------------------- 
1239      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1240      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1241      LOGICAL                              , INTENT(in   ) :: before
1242      !
1243      INTEGER :: ji, jj, jk
1244      !!---------------------------------------------------------------------- 
1245      !   
1246      IF( before ) THEN
1247         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1248      ELSE
1249         !
1250         DO jk = k1, k2
1251            DO jj = j1, j2
1252               DO ji = i1, i2
1253                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1254                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  & 
1255                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), &
1256                     &                 ji+nimpp-1, jj+njmpp-1, jk
1257                     kindic_agr = kindic_agr + 1
1258                  ENDIF
1259               END DO
1260            END DO
1261         END DO
1262         !
1263      ENDIF
1264      !
1265   END SUBROUTINE interpe3t
1266
1267   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1268      !!----------------------------------------------------------------------
1269      !!                  ***  ROUTINE interavm  ***
1270      !!---------------------------------------------------------------------- 
1271      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1272      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1273      LOGICAL                                    , INTENT(in   ) ::   before
1274      !
1275      INTEGER  :: ji, jj, jk
1276      INTEGER  :: N_in, N_out
1277      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in
1278      REAL(wp), DIMENSION(1:jpk) :: z_out
1279      !!---------------------------------------------------------------------- 
1280      !     
1281      IF (before) THEN         
1282         DO jk=k1,k2
1283            DO jj=j1,j2
1284              DO ji=i1,i2
1285                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1286              END DO
1287           END DO
1288         END DO
1289
1290         IF( l_vremap ) THEN
1291            ! Interpolate thicknesses
1292            ! Warning: these are masked, hence extrapolated prior interpolation.
1293            DO jk=k1,k2
1294               DO jj=j1,j2
1295                  DO ji=i1,i2
1296                      ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
1297                  END DO
1298               END DO
1299            END DO
1300
1301            ! Extrapolate thicknesses in partial bottom cells:
1302            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1303            IF (ln_zps) THEN
1304               DO jj=j1,j2
1305                  DO ji=i1,i2
1306                      jk = mbkt(ji,jj)
1307                      ptab(ji,jj,jk,2) = 0._wp
1308                  END DO
1309               END DO           
1310            END IF
1311       
1312           ! Save ssh at last level:
1313            IF (.NOT.ln_linssh) THEN
1314               ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
1315            ELSE
1316               ptab(i1:i2,j1:j2,k2,2) = 0._wp
1317            END IF     
1318          ENDIF
1319
1320      ELSE
1321
1322         IF( l_vremap ) THEN
1323            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1324            avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
1325               
1326            DO jj = j1, j2
1327               DO ji =i1, i2
1328                  N_in = mbkt_parent(ji,jj)
1329                  IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
1330                  z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
1331                  DO jk = N_in, 1, -1  ! Parent vertical grid               
1332                        z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
1333                       tabin(jk) = ptab(ji,jj,jk,1)
1334                  END DO
1335                  N_out = mbkt(ji,jj) 
1336                  DO jk = 1, N_out        ! Child vertical grid
1337                     z_out(jk) = gdepw(ji,jj,jk,Kmm_a)
1338                  ENDDO
1339                  IF (N_in*N_out > 0) THEN
1340                     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)
1341                  ENDIF
1342               ENDDO
1343            ENDDO
1344         ELSE
1345            avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1346         ENDIF
1347      ENDIF
1348      !
1349   END SUBROUTINE interpavm
1350
1351   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
1352      !!----------------------------------------------------------------------
1353      !!                  ***  ROUTINE interpsshn  ***
1354      !!---------------------------------------------------------------------- 
1355      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1356      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1357      LOGICAL                         , INTENT(in   ) ::   before
1358      !
1359      !!---------------------------------------------------------------------- 
1360      !
1361      IF( before) THEN
1362         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
1363      ELSE
1364         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
1365      ENDIF
1366      !
1367   END SUBROUTINE interpmbkt
1368
1369   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
1370      !!----------------------------------------------------------------------
1371      !!                  ***  ROUTINE interpsshn  ***
1372      !!---------------------------------------------------------------------- 
1373      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1374      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1375      LOGICAL                         , INTENT(in   ) ::   before
1376      !
1377      !!---------------------------------------------------------------------- 
1378      !
1379      IF( before) THEN
1380         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2)
1381      ELSE
1382         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
1383      ENDIF
1384      !
1385   END SUBROUTINE interpht0
1386
1387   SUBROUTINE agrif_initts(tabres,i1,i2,j1,j2,k1,k2,m1,m2,before)
1388       INTEGER :: i1, i2, j1, j2, k1, k2, m1, m2
1389       REAL(wp):: tabres(i1:i2,j1:j2,k1:k2,m1:m2)
1390       LOGICAL :: before
1391
1392       INTEGER :: jm
1393
1394       IF (before) THEN
1395         DO jm=1,jpts
1396             tabres(i1:i2,j1:j2,k1:k2,jm) = ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)
1397         END DO
1398       ELSE
1399         DO jm=1,jpts
1400             ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)=tabres(i1:i2,j1:j2,k1:k2,jm)
1401         END DO
1402       ENDIF
1403   END SUBROUTINE agrif_initts 
1404
1405   SUBROUTINE agrif_initssh( ptab, i1, i2, j1, j2, before )
1406      !!----------------------------------------------------------------------
1407      !!                  ***  ROUTINE interpsshn  ***
1408      !!---------------------------------------------------------------------- 
1409      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1410      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1411      LOGICAL                         , INTENT(in   ) ::   before
1412      !
1413      !!---------------------------------------------------------------------- 
1414      !
1415      IF( before) THEN
1416         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kbb_a)
1417      ELSE
1418         ssh(i1:i2,j1:j2,Kbb_a) = ptab(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1)
1419      ENDIF
1420      !
1421   END SUBROUTINE agrif_initssh
1422   
1423#else
1424   !!----------------------------------------------------------------------
1425   !!   Empty module                                          no AGRIF zoom
1426   !!----------------------------------------------------------------------
1427CONTAINS
1428   SUBROUTINE Agrif_OCE_Interp_empty
1429      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1430   END SUBROUTINE Agrif_OCE_Interp_empty
1431#endif
1432
1433   !!======================================================================
1434END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.