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_r12973_AGRIF_CMEMS/src/NST – NEMO

source: NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/NST/agrif_oce_interp.F90 @ 13162

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

#2129, corrections/add ons to initial state interpolation with AGRIF
1) add namelist flag for child grid initial state interpolation - ice not considered yet
2) provide depths and not thicknesses as inputs to vertical linear interpolation
3) extend initial state interpolation to a restart scenario for parent grid (warning should be added in that case in order to prevent users doing this at each model restart...)
The online interpolation seems to work fine in the VORTEX case (provided 0. is not considered as a special value in the initial velocity field, i.e. ln_spc_dyn=F)

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