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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/NST – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/NST/agrif_oce_interp.F90 @ 13678

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

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