source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce_interp.F90 @ 11219

Last change on this file since 11219 was 11219, checked in by jchanut, 15 months ago

#2199
1) Define aditionnal arrays to correct the time interpolation of barotropic arrays in corners. Since multiple stages in the time interpolation are necessary, overlapping segments in corners give wrong results otherwise (corrects stage 2 in previous commit)..
2) Added subroutine to correct time extrapolated fluxes at bdy in time splitting routine (updates stage 3 in previous commit).
3) Completly remove non-specified open boundary case. Boundares are now exactly set from parent (no more filtering nor extrapolation for outgoing flows).
At this stage, use of nbondi, nbondj variables has been completly removed.

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