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/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90 @ 11802

Last change on this file since 11802 was 11802, checked in by jchanut, 5 years ago

#2222, 1) add linear interpolation in vremap module.
2) Switch remapping of viscosity from polynomial to linear.
3) Move to truly volume weighted averages for parent to child update.

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