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 @ 11741

Last change on this file since 11741 was 11741, checked in by jchanut, 4 years ago

#2222: correct definition of parent vertical grid on the child domain to perform vertical interpolation at boundaries. Use additionnal parent depths and number of levels arrays interpolated on the child grid domain to do so.
Correction of vertical interpolation of viscosity remains to be done as well as duplication of changes for passive tracers.

  • Property svn:keywords set to Id
File size: 49.3 KB
Line 
1MODULE agrif_oce_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_oce_interp  ***
4   !! AGRIF: interpolation package for the ocean dynamics (OPA)
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (L. Debreu)  Original cade
7   !!            3.2  !  2009-04  (R. Benshila)
8   !!            3.6  !  2014-09  (R. Benshila)
9   !!----------------------------------------------------------------------
10#if defined key_agrif
11   !!----------------------------------------------------------------------
12   !!   'key_agrif'                                              AGRIF zoom
13   !!----------------------------------------------------------------------
14   !!   Agrif_tra     :
15   !!   Agrif_dyn     :
16   !!   Agrif_ssh     :
17   !!   Agrif_dyn_ts  :
18   !!   Agrif_dta_ts  :
19   !!   Agrif_ssh_ts  :
20   !!   Agrif_avm     :
21   !!   interpu       :
22   !!   interpv       :
23   !!----------------------------------------------------------------------
24   USE par_oce
25   USE oce
26   USE dom_oce     
27   USE zdf_oce
28   USE agrif_oce
29   USE phycst
30   USE dynspg_ts, ONLY: un_adv, vn_adv
31   !
32   USE in_out_manager
33   USE agrif_oce_sponge
34   USE lib_mpp
35   USE vremap
36 
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-1
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               zhtot = 0._wp
722               DO jk=1,N_in !k2 = jpk of parent grid
723                  IF (jk==N_in) THEN
724                     h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
725                  ELSE
726                     h_in(jk) = ptab(ji,jj,jk,n2)
727                  ENDIF
728                  zhtot = zhtot + h_in(jk)
729                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
730               END DO
731               N_out = 0
732               DO jk=1,jpk ! jpk of child grid
733                  IF (tmask(ji,jj,jk) == 0) EXIT
734                  N_out = N_out + 1
735                  h_out(jk) = e3t_a(ji,jj,jk)
736               ENDDO
737               IF (N_in*N_out > 0) THEN
738                  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)
739               ENDIF
740            ENDDO
741         ENDDO
742# else
743         !
744         DO jn=1, jpts
745            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
746         END DO
747# endif
748
749      ENDIF
750      !
751   END SUBROUTINE interptsn
752
753   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
754      !!----------------------------------------------------------------------
755      !!                  ***  ROUTINE interpsshn  ***
756      !!---------------------------------------------------------------------- 
757      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
758      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
759      LOGICAL                         , INTENT(in   ) ::   before
760      !
761      !!---------------------------------------------------------------------- 
762      !
763      IF( before) THEN
764         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
765      ELSE
766         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
767      ENDIF
768      !
769   END SUBROUTINE interpsshn
770
771   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
772      !!----------------------------------------------------------------------
773      !!                  *** ROUTINE interpun ***
774      !!---------------------------------------------   
775      !!
776      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
777      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
778      LOGICAL, INTENT(in) :: before
779      !!
780      INTEGER :: ji,jj,jk
781      REAL(wp) :: zrhoy, zhtot
782      ! vertical interpolation:
783      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
784      REAL(wp), DIMENSION(1:jpk) :: h_out
785      INTEGER  :: N_in, N_out
786      REAL(wp) :: h_diff
787      !!---------------------------------------------   
788      !
789      IF (before) THEN
790         DO jk=1,jpk
791            DO jj=j1,j2
792               DO ji=i1,i2
793                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
794# if defined key_vertical
795                  ! Interpolate thicknesses (masked for subsequent extrapolation)
796                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
797# endif
798               END DO
799            END DO
800         END DO
801# if defined key_vertical
802         ! Extrapolate thicknesses in partial bottom cells:
803         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
804         IF (ln_zps) THEN
805            DO jj=j1,j2
806               DO ji=i1,i2
807                  jk = mbku(ji,jj)
808                  ptab(ji,jj,jk,2) = 0._wp
809               END DO
810            END DO           
811         END IF
812        ! Save ssh at last level:
813        ptab(i1:i2,j1:j2,k2,2) = 0._wp
814        IF (.NOT.ln_linssh) THEN
815           ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
816           DO jk=1,jpk
817              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)
818           END DO
819           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
820        END IF 
821# endif
822         !
823      ELSE
824         zrhoy = Agrif_rhoy()
825# if defined key_vertical
826! VERTICAL REFINEMENT BEGIN
827
828         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
829
830         DO ji=i1,i2
831            DO jj=j1,j2
832               ua(ji,jj,:) = 0._wp
833               N_in = mbku_parent(ji,jj)
834               zhtot = 0._wp
835               DO jk=1,N_in
836                  IF (jk==N_in) THEN
837                     h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
838                  ELSE
839                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
840                  ENDIF
841                  zhtot = zhtot + h_in(jk)
842                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
843              ENDDO
844                 
845              N_out = 0
846              DO jk=1,jpk
847                 if (umask(ji,jj,jk) == 0) EXIT
848                 N_out = N_out + 1
849                 h_out(N_out) = e3u_a(ji,jj,jk)
850              ENDDO
851              IF (N_in*N_out > 0) THEN
852                 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)
853              ENDIF
854            ENDDO
855         ENDDO
856
857# else
858         DO jk = 1, jpkm1
859            DO jj=j1,j2
860               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
861            END DO
862         END DO
863# endif
864
865      ENDIF
866      !
867   END SUBROUTINE interpun
868
869   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
870      !!----------------------------------------------------------------------
871      !!                  *** ROUTINE interpvn ***
872      !!----------------------------------------------------------------------
873      !
874      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
875      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
876      LOGICAL, INTENT(in) :: before
877      !
878      INTEGER :: ji,jj,jk
879      REAL(wp) :: zrhox
880      ! vertical interpolation:
881      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
882      REAL(wp), DIMENSION(1:jpk) :: h_out
883      INTEGER  :: N_in, N_out
884      REAL(wp) :: h_diff, zhtot
885      !!---------------------------------------------   
886      !     
887      IF (before) THEN         
888         DO jk=k1,k2
889            DO jj=j1,j2
890               DO ji=i1,i2
891                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
892# if defined key_vertical
893                  ! Interpolate thicknesses (masked for subsequent extrapolation)
894                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
895# endif
896               END DO
897            END DO
898         END DO
899# if defined key_vertical
900         ! Extrapolate thicknesses in partial bottom cells:
901         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
902         IF (ln_zps) THEN
903            DO jj=j1,j2
904               DO ji=i1,i2
905                  jk = mbkv(ji,jj)
906                  ptab(ji,jj,jk,2) = 0._wp
907               END DO
908            END DO           
909         END IF
910        ! Save ssh at last level:
911        ptab(i1:i2,j1:j2,k2,2) = 0._wp
912        IF (.NOT.ln_linssh) THEN
913           ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
914           DO jk=1,jpk
915              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)
916           END DO
917           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
918        END IF 
919# endif
920      ELSE       
921         zrhox = Agrif_rhox()
922# if defined key_vertical
923
924         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
925
926         DO jj=j1,j2
927            DO ji=i1,i2
928               va(ji,jj,:) = 0._wp
929               N_in = mbkv_parent(ji,jj)
930               zhtot = 0._wp
931               DO jk=1,N_in
932                  IF (jk==N_in) THEN
933                     h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
934                  ELSE
935                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
936                  ENDIF
937                  zhtot = zhtot + h_in(jk)
938                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
939              ENDDO
940         
941               N_out = 0
942               DO jk=1,jpk
943                  if (vmask(ji,jj,jk) == 0) EXIT
944                  N_out = N_out + 1
945                  h_out(N_out) = e3v_a(ji,jj,jk)
946               END DO
947
948               IF (N_in*N_out > 0) THEN
949                  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)
950               ENDIF
951            END DO
952         END DO
953# else
954         DO jk = 1, jpkm1
955            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) )
956         END DO
957# endif
958      ENDIF
959      !       
960   END SUBROUTINE interpvn
961
962   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
963      !!----------------------------------------------------------------------
964      !!                  ***  ROUTINE interpunb  ***
965      !!---------------------------------------------------------------------- 
966      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
967      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
968      LOGICAL                         , INTENT(in   ) ::   before
969      !
970      INTEGER  ::   ji, jj
971      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
972      !!---------------------------------------------------------------------- 
973      !
974      IF( before ) THEN
975         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
976      ELSE
977         zrhoy = Agrif_Rhoy()
978         zrhot = Agrif_rhot()
979         ! Time indexes bounds for integration
980         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
981         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
982         !
983         DO ji = i1, i2
984            DO jj = j1, j2
985               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
986                  IF    ( utint_stage(ji,jj) == 1  ) THEN
987                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
988                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
989                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
990                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
991                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
992                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
993                     ztcoeff = 1._wp
994                  ELSE
995                     ztcoeff = 0._wp
996                  ENDIF
997                  !   
998                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
999                  !           
1000                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
1001                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
1002                  ENDIF
1003                  !
1004                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
1005               ENDIF
1006            END DO
1007         END DO
1008      END IF
1009      !
1010   END SUBROUTINE interpunb
1011
1012
1013   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
1014      !!----------------------------------------------------------------------
1015      !!                  ***  ROUTINE interpvnb  ***
1016      !!---------------------------------------------------------------------- 
1017      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1018      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1019      LOGICAL                         , INTENT(in   ) ::   before
1020      !
1021      INTEGER  ::   ji, jj
1022      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1023      !!---------------------------------------------------------------------- 
1024      !
1025      IF( before ) THEN
1026         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1027      ELSE
1028         zrhox = Agrif_Rhox()
1029         zrhot = Agrif_rhot()
1030         ! Time indexes bounds for integration
1031         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1032         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
1033         !     
1034         DO ji = i1, i2
1035            DO jj = j1, j2
1036               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1037                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
1038                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1039                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1040                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
1041                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1042                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1043                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
1044                     ztcoeff = 1._wp
1045                  ELSE
1046                     ztcoeff = 0._wp
1047                  ENDIF
1048                  !   
1049                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
1050                  !           
1051                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
1052                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
1053                  ENDIF
1054                  !
1055                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
1056               ENDIF
1057            END DO
1058         END DO         
1059      ENDIF
1060      !
1061   END SUBROUTINE interpvnb
1062
1063
1064   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
1065      !!----------------------------------------------------------------------
1066      !!                  ***  ROUTINE interpub2b  ***
1067      !!---------------------------------------------------------------------- 
1068      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1069      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1070      LOGICAL                         , INTENT(in   ) ::   before
1071      !
1072      INTEGER  ::   ji,jj
1073      REAL(wp) ::   zrhot, zt0, zt1, zat
1074      !!---------------------------------------------------------------------- 
1075      IF( before ) THEN
1076         IF ( ln_bt_fw ) THEN
1077            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1078         ELSE
1079            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1080         ENDIF
1081      ELSE
1082         zrhot = Agrif_rhot()
1083         ! Time indexes bounds for integration
1084         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1085         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1086         ! Polynomial interpolation coefficients:
1087         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1088            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1089         !
1090         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1091         !
1092         ! Update interpolation stage:
1093         utint_stage(i1:i2,j1:j2) = 1
1094      ENDIF
1095      !
1096   END SUBROUTINE interpub2b
1097   
1098
1099   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
1100      !!----------------------------------------------------------------------
1101      !!                  ***  ROUTINE interpvb2b  ***
1102      !!---------------------------------------------------------------------- 
1103      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1104      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1105      LOGICAL                         , INTENT(in   ) ::   before
1106      !
1107      INTEGER ::   ji,jj
1108      REAL(wp) ::   zrhot, zt0, zt1, zat
1109      !!---------------------------------------------------------------------- 
1110      !
1111      IF( before ) THEN
1112         IF ( ln_bt_fw ) THEN
1113            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1114         ELSE
1115            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1116         ENDIF
1117      ELSE     
1118         zrhot = Agrif_rhot()
1119         ! Time indexes bounds for integration
1120         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1121         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1122         ! Polynomial interpolation coefficients:
1123         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1124            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1125         !
1126         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1127         !
1128         ! update interpolation stage:
1129         vtint_stage(i1:i2,j1:j2) = 1
1130      ENDIF
1131      !     
1132   END SUBROUTINE interpvb2b
1133
1134
1135   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1136      !!----------------------------------------------------------------------
1137      !!                  ***  ROUTINE interpe3t  ***
1138      !!---------------------------------------------------------------------- 
1139      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1140      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1141      LOGICAL                              , INTENT(in   ) :: before
1142      INTEGER                              , INTENT(in   ) :: nb , ndir
1143      !
1144      INTEGER :: ji, jj, jk
1145      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1146      !!---------------------------------------------------------------------- 
1147      !   
1148      IF( before ) THEN
1149         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1150      ELSE
1151         western_side  = (nb == 1).AND.(ndir == 1)
1152         eastern_side  = (nb == 1).AND.(ndir == 2)
1153         southern_side = (nb == 2).AND.(ndir == 1)
1154         northern_side = (nb == 2).AND.(ndir == 2)
1155         !
1156         DO jk = k1, k2
1157            DO jj = j1, j2
1158               DO ji = i1, i2
1159                  !
1160                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1161                     IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN
1162                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1163                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk) 
1164                        kindic_agr = kindic_agr + 1
1165                     ELSEIF (eastern_side.AND.(ptab(i2-nbghostcells+1,jj,jk)>0._wp)) THEN
1166                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1167                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1168                        kindic_agr = kindic_agr + 1
1169                     ELSEIF (southern_side.AND.(ptab(ji,j1+nbghostcells-1,jk)>0._wp)) THEN
1170                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1171                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1172                        kindic_agr = kindic_agr + 1
1173                     ELSEIF (northern_side.AND.(ptab(ji,j2-nbghostcells+1,jk)>0._wp)) THEN
1174                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1175                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1176                        kindic_agr = kindic_agr + 1
1177                     ENDIF
1178                  ENDIF
1179               END DO
1180            END DO
1181         END DO
1182         !
1183      ENDIF
1184      !
1185   END SUBROUTINE interpe3t
1186
1187
1188   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1189      !!----------------------------------------------------------------------
1190      !!                  ***  ROUTINE interpumsk  ***
1191      !!---------------------------------------------------------------------- 
1192      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1193      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1194      LOGICAL                              , INTENT(in   ) ::   before
1195      INTEGER                              , INTENT(in   ) ::   nb , ndir
1196      !
1197      INTEGER ::   ji, jj, jk
1198      LOGICAL ::   western_side, eastern_side   
1199      !!---------------------------------------------------------------------- 
1200      !   
1201      IF( before ) THEN
1202         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1203      ELSE
1204         western_side = (nb == 1).AND.(ndir == 1)
1205         eastern_side = (nb == 1).AND.(ndir == 2)
1206         DO jk = k1, k2
1207            DO jj = j1, j2
1208               DO ji = i1, i2
1209                   ! Velocity mask at boundary edge points:
1210                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1211                     IF (western_side) THEN
1212                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1213                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1214                        kindic_agr = kindic_agr + 1
1215                     ELSEIF (eastern_side) THEN
1216                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1217                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1218                        kindic_agr = kindic_agr + 1
1219                     ENDIF
1220                  ENDIF
1221               END DO
1222            END DO
1223         END DO
1224         !
1225      ENDIF
1226      !
1227   END SUBROUTINE interpumsk
1228
1229
1230   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1231      !!----------------------------------------------------------------------
1232      !!                  ***  ROUTINE interpvmsk  ***
1233      !!---------------------------------------------------------------------- 
1234      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1235      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1236      LOGICAL                              , INTENT(in   ) ::   before
1237      INTEGER                              , INTENT(in   ) :: nb , ndir
1238      !
1239      INTEGER ::   ji, jj, jk
1240      LOGICAL ::   northern_side, southern_side     
1241      !!---------------------------------------------------------------------- 
1242      !   
1243      IF( before ) THEN
1244         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1245      ELSE
1246         southern_side = (nb == 2).AND.(ndir == 1)
1247         northern_side = (nb == 2).AND.(ndir == 2)
1248         DO jk = k1, k2
1249            DO jj = j1, j2
1250               DO ji = i1, i2
1251                   ! Velocity mask at boundary edge points:
1252                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1253                     IF (southern_side) THEN
1254                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1255                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1256                        kindic_agr = kindic_agr + 1
1257                     ELSEIF (northern_side) THEN
1258                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1259                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1260                        kindic_agr = kindic_agr + 1
1261                     ENDIF
1262                  ENDIF
1263               END DO
1264            END DO
1265         END DO
1266         !
1267      ENDIF
1268      !
1269   END SUBROUTINE interpvmsk
1270
1271
1272   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1273      !!----------------------------------------------------------------------
1274      !!                  ***  ROUTINE interavm  ***
1275      !!---------------------------------------------------------------------- 
1276      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1277      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1278      LOGICAL                                    , INTENT(in   ) ::   before
1279      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
1280      REAL(wp), DIMENSION(1:jpk) :: h_out
1281      INTEGER  :: N_in, N_out, ji, jj, jk
1282      !!---------------------------------------------------------------------- 
1283      !     
1284      IF (before) THEN         
1285         DO jk=k1,k2
1286            DO jj=j1,j2
1287              DO ji=i1,i2
1288                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1289              END DO
1290           END DO
1291        END DO
1292#ifdef key_vertical         
1293        DO jk=k1,k2
1294           DO jj=j1,j2
1295              DO ji=i1,i2
1296                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk) 
1297              END DO
1298           END DO
1299        END DO
1300#endif
1301      ELSE 
1302#ifdef key_vertical         
1303         avm_k(i1:i2,j1:j2,1:jpk) = 0.
1304         DO jj=j1,j2
1305            DO ji=i1,i2
1306               N_in = 0
1307               DO jk=k1,k2 !k2 = jpk of parent grid
1308                  IF (ptab(ji,jj,jk,2) == 0) EXIT
1309                  N_in = N_in + 1
1310                  tabin(jk) = ptab(ji,jj,jk,1)
1311                  h_in(N_in) = ptab(ji,jj,jk,2)
1312               END DO
1313               N_out = 0
1314               DO jk=1,jpk ! jpk of child grid
1315                  IF (wmask(ji,jj,jk) == 0) EXIT
1316                  N_out = N_out + 1
1317                  h_out(jk) = e3t_n(ji,jj,jk)
1318               ENDDO
1319               IF (N_in > 0) THEN
1320                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out,1)
1321               ENDIF
1322            ENDDO
1323         ENDDO
1324#else
1325         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1326#endif
1327      ENDIF
1328      !
1329   END SUBROUTINE interpavm
1330
1331# if defined key_vertical
1332   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
1333      !!----------------------------------------------------------------------
1334      !!                  ***  ROUTINE interpsshn  ***
1335      !!---------------------------------------------------------------------- 
1336      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1337      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1338      LOGICAL                         , INTENT(in   ) ::   before
1339      !
1340      !!---------------------------------------------------------------------- 
1341      !
1342      IF( before) THEN
1343         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
1344      ELSE
1345         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
1346      ENDIF
1347      !
1348   END SUBROUTINE interpmbkt
1349
1350   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
1351      !!----------------------------------------------------------------------
1352      !!                  ***  ROUTINE interpsshn  ***
1353      !!---------------------------------------------------------------------- 
1354      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1355      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1356      LOGICAL                         , INTENT(in   ) ::   before
1357      !
1358      !!---------------------------------------------------------------------- 
1359      !
1360      IF( before) THEN
1361         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2)
1362      ELSE
1363         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
1364      ENDIF
1365      !
1366   END SUBROUTINE interpht0
1367#endif
1368
1369#else
1370   !!----------------------------------------------------------------------
1371   !!   Empty module                                          no AGRIF zoom
1372   !!----------------------------------------------------------------------
1373CONTAINS
1374   SUBROUTINE Agrif_OCE_Interp_empty
1375      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1376   END SUBROUTINE Agrif_OCE_Interp_empty
1377#endif
1378
1379   !!======================================================================
1380END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.