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

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

#2222, remove useless mask checking (and restrict scale factor check at the boundary only until nesting tools are updated in sponge areas). Take into account special values in tracer updates, again, till nesting tools are updated.

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