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

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

#2222, 1) correct time interpolation of barotropic velocities in corners. 2) Clean remapping module and enable remapping several variables at the same time. At this stage, vertical remapping doesn't change VORTEX results with an identical vertical grid ONLY in one way mode and a linearized free surface (within truncature errors).

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