source: NEMO/trunk/src/NST/agrif_oce_interp.F90

Last change on this file was 15119, checked in by jchanut, 3 months ago

#2638, changes to accomodate nn_hls=2 and AGRIF zooms crossing cyclic boundaries. E-W case ok, update for North-Fold still needed.

  • Property svn:keywords set to Id
File size: 74.1 KB
Line 
1MODULE agrif_oce_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_oce_interp  ***
4   !! AGRIF: interpolation package for the ocean dynamics (OCE)
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   interpglamt, interpgphit
47   PUBLIC   interpht0, interpmbkt, interpe3t0_vremap
48   PUBLIC   agrif_istate_oce, agrif_istate_ssh   ! called by icestate.F90 and domvvl.F90
49   PUBLIC   agrif_check_bat
50
51   INTEGER ::   bdy_tinterp = 0
52
53   !! * Substitutions
54#  include "domzgr_substitute.h90"
55   !! NEMO/NST 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE Agrif_istate_oce( Kbb, Kmm, Kaa )
62      !!----------------------------------------------------------------------
63      !!                 *** ROUTINE agrif_istate_oce ***
64      !!
65      !!                 set initial t, s, u, v, ssh from parent
66      !!----------------------------------------------------------------------
67      !
68      IMPLICIT NONE
69      !
70      INTEGER, INTENT(in)  :: Kbb, Kmm, Kaa
71      INTEGER :: jn
72      !!----------------------------------------------------------------------
73      IF(lwp) WRITE(numout,*) ' '
74      IF(lwp) WRITE(numout,*) 'Agrif_istate_oce : interp child initial state from parent'
75      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~'
76      IF(lwp) WRITE(numout,*) ' '
77
78      IF ( .NOT.Agrif_Parent(l_1st_euler) ) & 
79         & CALL ctl_stop('AGRIF hot start requires to force Euler first step on parent')
80
81      l_ini_child           = .TRUE.
82      Agrif_SpecialValue    = 0.0_wp
83      Agrif_UseSpecialValue = .TRUE.
84
85      ts(:,:,:,:,Kbb) = 0.0_wp
86      uu(:,:,:,Kbb)   = 0.0_wp
87      vv(:,:,:,Kbb)   = 0.0_wp 
88       
89      Krhs_a = Kbb   ;   Kmm_a = Kbb
90
91      CALL Agrif_Init_Variable(tsini_id, procname=interptsn)
92
93      Agrif_UseSpecialValue = ln_spc_dyn
94      use_sign_north = .TRUE.
95      sign_north = -1._wp
96      CALL Agrif_Init_Variable(uini_id , procname=interpun )
97      CALL Agrif_Init_Variable(vini_id , procname=interpvn )
98      use_sign_north = .FALSE.
99
100      Agrif_UseSpecialValue = .FALSE.
101      l_ini_child           = .FALSE.
102
103      Krhs_a = Kaa   ;   Kmm_a = Kmm
104
105      DO jn = 1, jpts
106         ts(:,:,:,jn,Kbb) = ts(:,:,:,jn,Kbb) * tmask(:,:,:)
107      END DO
108      uu(:,:,:,Kbb) = uu(:,:,:,Kbb) * umask(:,:,:)     
109      vv(:,:,:,Kbb) = vv(:,:,:,Kbb) * vmask(:,:,:) 
110
111      CALL lbc_lnk( 'agrif_istate_oce', uu(:,:,:  ,Kbb), 'U', -1.0_wp , vv(:,:,:,Kbb), 'V', -1.0_wp )
112      CALL lbc_lnk( 'agrif_istate_oce', ts(:,:,:,:,Kbb), 'T',  1.0_wp )
113
114   END SUBROUTINE Agrif_istate_oce
115
116
117   SUBROUTINE Agrif_istate_ssh( Kbb, Kmm, Kaa )
118      !!----------------------------------------------------------------------
119      !!                 *** ROUTINE agrif_istate_ssh ***
120      !!
121      !!                    set initial ssh from parent
122      !!----------------------------------------------------------------------
123      !
124      IMPLICIT NONE
125      !
126      INTEGER, INTENT(in)  :: Kbb, Kmm, Kaa 
127      !!----------------------------------------------------------------------
128      IF(lwp) WRITE(numout,*) ' '
129      IF(lwp) WRITE(numout,*) 'Agrif_istate_ssh : interp child ssh from parent'
130      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~'
131      IF(lwp) WRITE(numout,*) ' '
132
133      IF ( .NOT.Agrif_Parent(l_1st_euler) ) & 
134         & CALL ctl_stop('AGRIF hot start requires to force Euler first step on parent')
135
136      Krhs_a = Kbb   ;   Kmm_a = Kbb
137      !
138      Agrif_SpecialValue    = 0._wp
139      Agrif_UseSpecialValue = .TRUE.
140      l_ini_child           = .TRUE.
141      !
142      ssh(:,:,Kbb) = 0._wp
143      CALL Agrif_Init_Variable(sshini_id, procname=interpsshn)
144      !
145      Agrif_UseSpecialValue = .FALSE.
146      l_ini_child           = .FALSE.
147      !
148      Krhs_a = Kaa   ;   Kmm_a = Kmm
149      !
150      CALL lbc_lnk( 'Agrif_istate_ssh', ssh(:,:,Kbb), 'T', 1._wp )
151      !
152      ssh(:,:,Kmm) = ssh(:,:,Kbb)
153      ssh(:,:,Kaa) = 0._wp
154
155   END SUBROUTINE Agrif_istate_ssh
156
157
158   SUBROUTINE Agrif_tra
159      !!----------------------------------------------------------------------
160      !!                  ***  ROUTINE Agrif_tra  ***
161      !!----------------------------------------------------------------------
162      !
163      IF( Agrif_Root() )   RETURN
164      !
165      Agrif_SpecialValue    = 0._wp
166      Agrif_UseSpecialValue = .TRUE.
167      l_vremap           = ln_vert_remap
168      !
169      CALL Agrif_Bc_variable( ts_interp_id, procname=interptsn )
170      !
171      Agrif_UseSpecialValue = .FALSE.
172      l_vremap              = .FALSE.
173      !
174   END SUBROUTINE Agrif_tra
175
176
177   SUBROUTINE Agrif_dyn( kt )
178      !!----------------------------------------------------------------------
179      !!                  ***  ROUTINE Agrif_DYN  ***
180      !!---------------------------------------------------------------------- 
181      INTEGER, INTENT(in) ::   kt
182      !
183      INTEGER ::   ji, jj, jk       ! dummy loop indices
184      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
185      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
186      !!---------------------------------------------------------------------- 
187      !
188      IF( Agrif_Root() )   RETURN
189      !
190      Agrif_SpecialValue    = 0.0_wp
191      Agrif_UseSpecialValue = ln_spc_dyn
192      l_vremap              = ln_vert_remap
193      !
194      use_sign_north = .TRUE.
195      sign_north = -1.0_wp
196      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
197      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
198
199      IF( .NOT.ln_dynspg_ts ) THEN ! Get transports
200         ubdy(:,:) = 0._wp    ;  vbdy(:,:) = 0._wp
201         utint_stage(:,:) = 0 ;  vtint_stage(:,:) = 0
202         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb )
203         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb )
204      ENDIF
205
206      use_sign_north = .FALSE.
207      !
208      Agrif_UseSpecialValue = .FALSE.
209      l_vremap              = .FALSE.
210      !
211      ! Ensure below that vertically integrated transports match
212      ! either transports out of time splitting procedure (ln_dynspg_ts=.TRUE.)
213      ! or parent grid transports (ln_dynspg_ts=.FALSE.)
214      !
215      ! --- West --- !
216      IF( lk_west ) THEN
217         ibdy1 = nn_hls + 2                  ! halo + land + 1
218         ibdy2 = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhox()   ! halo + land + nbghostcells
219         !
220         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
221            DO ji = mi0(ibdy1), mi1(ibdy2)
222               DO jj = 1, jpj
223                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
224                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
225               END DO
226            END DO
227         ENDIF
228         !
229         DO ji = mi0(ibdy1), mi1(ibdy2)
230            zub(ji,:) = 0._wp 
231            DO jk = 1, jpkm1
232               DO jj = 1, jpj
233                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
234               END DO
235            END DO
236            DO jj=1,jpj
237               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
238            END DO
239            DO jk = 1, jpkm1
240               DO jj = 1, jpj
241                  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)
242               END DO
243            END DO
244         END DO
245         !   
246         DO ji = mi0(ibdy1), mi1(ibdy2)
247            zvb(ji,:) = 0._wp
248            DO jk = 1, jpkm1
249               DO jj = 1, jpj
250                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
251               END DO
252            END DO
253            DO jj = 1, jpj
254               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
255            END DO
256            DO jk = 1, jpkm1
257               DO jj = 1, jpj
258                  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)
259               END DO
260            END DO
261         END DO
262         !
263      ENDIF
264
265      ! --- East --- !
266      IF( lk_east) THEN
267         ibdy1 = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox()   
268         ibdy2 = jpiglo - ( nn_hls + 2 )                 
269         !
270         IF( .NOT.ln_dynspg_ts ) THEN
271            DO ji = mi0(ibdy1), mi1(ibdy2)
272               DO jj = 1, jpj
273                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
274               END DO
275            END DO
276         ENDIF
277         !
278         DO ji = mi0(ibdy1), mi1(ibdy2)
279            zub(ji,:) = 0._wp   
280            DO jk = 1, jpkm1
281               DO jj = 1, jpj
282                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
283               END DO
284            END DO
285            DO jj=1,jpj
286               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
287            END DO
288            DO jk = 1, jpkm1
289               DO jj = 1, jpj
290                  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)
291               END DO
292            END DO
293         END DO
294         !
295         ibdy1 = jpiglo - ( nn_hls + nbghostcells - 1 ) - nn_shift_bar*Agrif_Rhox() 
296         ibdy2 = jpiglo - ( nn_hls + 1 )     
297         !
298         IF( .NOT.ln_dynspg_ts ) THEN
299            DO ji = mi0(ibdy1), mi1(ibdy2)
300               DO jj = 1, jpj
301                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
302               END DO
303            END DO
304         ENDIF
305         !
306         DO ji = mi0(ibdy1), mi1(ibdy2)
307            zvb(ji,:) = 0._wp
308            DO jk = 1, jpkm1
309               DO jj = 1, jpj
310                     zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
311               END DO
312            END DO
313            DO jj = 1, jpj
314               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
315            END DO
316            DO jk = 1, jpkm1
317               DO jj = 1, jpj
318                     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)
319               END DO
320            END DO
321         END DO
322         !
323      ENDIF
324
325      ! --- South --- !
326      IF( lk_south ) THEN
327         jbdy1 = nn_hls + 2                 
328         jbdy2 = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhoy()   
329         !
330         IF( .NOT.ln_dynspg_ts ) THEN
331            DO jj = mj0(jbdy1), mj1(jbdy2)
332               DO ji = 1, jpi
333                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
334                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
335               END DO
336            END DO
337         ENDIF
338         !
339         DO jj = mj0(jbdy1), mj1(jbdy2)
340            zvb(:,jj) = 0._wp
341            DO jk=1,jpkm1
342               DO ji=1,jpi
343                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
344               END DO
345            END DO
346            DO ji = 1, jpi
347               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
348            END DO
349            DO jk = 1, jpkm1
350               DO ji = 1, jpi
351                  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)
352               END DO
353            END DO
354         END DO
355         !
356         DO jj = mj0(jbdy1), mj1(jbdy2)
357            zub(:,jj) = 0._wp
358            DO jk = 1, jpkm1
359               DO ji = 1, jpi
360                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
361               END DO
362            END DO
363            DO ji = 1, jpi
364               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
365            END DO
366            DO jk = 1, jpkm1
367               DO ji = 1, jpi
368                  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)
369               END DO
370            END DO
371         END DO
372         !
373      ENDIF
374
375      ! --- North --- !
376      IF( lk_north ) THEN
377         jbdy1 = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy() 
378         jbdy2 = jpjglo - ( nn_hls + 2 )
379         !
380         IF( .NOT.ln_dynspg_ts ) THEN
381            DO jj = mj0(jbdy1), mj1(jbdy2)
382               DO ji = 1, jpi
383                  vv_b(ji,jj,Krhs_a) = vbdy(ji,jj) * r1_hv(ji,jj,Krhs_a)
384               END DO
385            END DO
386         ENDIF
387         !
388         DO jj = mj0(jbdy1), mj1(jbdy2)
389            zvb(:,jj) = 0._wp 
390            DO jk=1,jpkm1
391               DO ji=1,jpi
392                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk)
393               END DO
394            END DO
395            DO ji = 1, jpi
396               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a)
397            END DO
398            DO jk = 1, jpkm1
399               DO ji = 1, jpi
400                  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)
401               END DO
402            END DO
403         END DO
404         !
405         jbdy1 = jpjglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhoy() 
406         jbdy2 = jpjglo - ( nn_hls + 1 )
407         !
408         IF( .NOT.ln_dynspg_ts ) THEN
409            DO jj = mj0(jbdy1), mj1(jbdy2)
410               DO ji = 1, jpi
411                  uu_b(ji,jj,Krhs_a) = ubdy(ji,jj) * r1_hu(ji,jj,Krhs_a)
412               END DO
413            END DO
414         ENDIF
415         !
416         DO jj = mj0(jbdy1), mj1(jbdy2)
417            zub(:,jj) = 0._wp
418            DO jk = 1, jpkm1
419               DO ji = 1, jpi
420                  zub(ji,jj) = zub(ji,jj) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk)
421               END DO
422            END DO
423            DO ji = 1, jpi
424               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
425            END DO
426            DO jk = 1, jpkm1
427               DO ji = 1, jpi
428                     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)
429               END DO
430            END DO
431         END DO
432         !
433      ENDIF
434      !
435   END SUBROUTINE Agrif_dyn
436
437
438   SUBROUTINE Agrif_dyn_ts( jn )
439      !!----------------------------------------------------------------------
440      !!                  ***  ROUTINE Agrif_dyn_ts  ***
441      !!---------------------------------------------------------------------- 
442      INTEGER, INTENT(in) ::   jn
443      !!
444      INTEGER :: ji, jj
445      INTEGER :: istart, iend, jstart, jend
446      !!---------------------------------------------------------------------- 
447      !
448      IF( Agrif_Root() )   RETURN
449      !
450      !--- West ---!
451      IF( lk_west ) THEN
452         istart = nn_hls + 2                              ! halo + land + 1
453         iend   = nn_hls + nbghostcells  + nn_shift_bar*Agrif_Rhox()              ! halo + land + nbghostcells
454         DO ji = mi0(istart), mi1(iend)
455            DO jj=1,jpj
456               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
457               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
458            END DO
459         END DO
460      ENDIF
461      !
462      !--- East ---!
463      IF( lk_east ) THEN
464         istart = jpiglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhox() 
465         iend   = jpiglo - ( nn_hls + 1 )               
466         DO ji = mi0(istart), mi1(iend)
467
468            DO jj=1,jpj
469               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
470            END DO
471         END DO
472         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox() 
473         iend   = jpiglo - ( nn_hls + 2 )               
474         DO ji = mi0(istart), mi1(iend)
475            DO jj=1,jpj
476               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
477            END DO
478         END DO
479      ENDIF 
480      !
481      !--- South ---!
482      IF( lk_south ) THEN
483         jstart = nn_hls + 2                             
484         jend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhoy()           
485         DO jj = mj0(jstart), mj1(jend)
486
487            DO ji=1,jpi
488               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
489               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
490            END DO
491         END DO
492      ENDIF       
493      !
494      !--- North ---!
495      IF( lk_north ) THEN
496         jstart = jpjglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhoy()     
497         jend   = jpjglo - ( nn_hls + 1 )               
498         DO jj = mj0(jstart), mj1(jend)
499            DO ji=1,jpi
500               ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
501            END DO
502         END DO
503         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy() 
504         jend   = jpjglo - ( nn_hls + 2 )               
505         DO jj = mj0(jstart), mj1(jend)
506            DO ji=1,jpi
507               va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
508            END DO
509         END DO
510      ENDIF 
511      !
512   END SUBROUTINE Agrif_dyn_ts
513
514   
515   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
516      !!----------------------------------------------------------------------
517      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
518      !!---------------------------------------------------------------------- 
519      INTEGER, INTENT(in) ::   jn
520      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
521      !!
522      INTEGER :: ji, jj
523      INTEGER :: istart, iend, jstart, jend
524      !!---------------------------------------------------------------------- 
525      !
526      IF( Agrif_Root() )   RETURN
527      !
528      !--- West ---!
529      IF( lk_west ) THEN
530         istart = nn_hls + 2                             
531         iend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhox() 
532         DO ji = mi0(istart), mi1(iend)
533            DO jj=1,jpj
534               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
535               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
536            END DO
537         END DO
538      ENDIF
539      !
540      !--- East ---!
541      IF( lk_east ) THEN
542         istart = jpiglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhox()
543         iend   = jpiglo - ( nn_hls + 1 )                 
544         DO ji = mi0(istart), mi1(iend)
545            DO jj=1,jpj
546               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
547            END DO
548         END DO
549         istart = jpiglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhox() 
550         iend   = jpiglo - ( nn_hls + 2 )                 
551         DO ji = mi0(istart), mi1(iend)
552            DO jj=1,jpj
553               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
554            END DO
555         END DO
556      ENDIF
557      !
558      !--- South ---!
559      IF( lk_south ) THEN
560         jstart = nn_hls + 2                             
561         jend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhoy() 
562         DO jj = mj0(jstart), mj1(jend)
563            DO ji=1,jpi
564               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
565               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
566            END DO
567         END DO
568      ENDIF
569      !
570      !--- North ---!
571      IF( lk_north ) THEN
572         jstart = jpjglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhoy() 
573         jend   = jpjglo - ( nn_hls + 1 )               
574         DO jj = mj0(jstart), mj1(jend)
575            DO ji=1,jpi
576               zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
577            END DO
578         END DO
579         jstart = jpjglo - ( nn_hls + nbghostcells ) - nn_shift_bar*Agrif_Rhoy() 
580         jend   = jpjglo - ( nn_hls + 2 )               
581         DO jj = mj0(jstart), mj1(jend)
582            DO ji=1,jpi
583               zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
584            END DO
585         END DO
586      ENDIF
587      !
588   END SUBROUTINE Agrif_dyn_ts_flux
589
590   
591   SUBROUTINE Agrif_dta_ts( kt )
592      !!----------------------------------------------------------------------
593      !!                  ***  ROUTINE Agrif_dta_ts  ***
594      !!---------------------------------------------------------------------- 
595      INTEGER, INTENT(in) ::   kt
596      !!
597      LOGICAL :: ll_int_cons
598      !!---------------------------------------------------------------------- 
599      !
600      IF( Agrif_Root() )   RETURN
601      !
602      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
603      !
604      ! Enforce volume conservation if no time refinement: 
605      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
606      !
607      ! Interpolate barotropic fluxes
608      Agrif_SpecialValue = 0._wp
609      Agrif_UseSpecialValue = ln_spc_dyn
610
611      use_sign_north = .TRUE.
612      sign_north = -1.
613
614      !
615      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
616      utint_stage(:,:) = 0
617      vtint_stage(:,:) = 0
618      !
619      IF( ll_int_cons ) THEN  ! Conservative interpolation
620         IF ( lk_tint2d_notinterp ) THEN
621            Agrif_UseSpecialValue = .FALSE.
622            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b_const )
623            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b_const ) 
624            ! Divergence conserving correction terms:
625            IF ( Agrif_Rhox()>1 ) CALL Agrif_Bc_variable(    ub2b_cor_id, calledweight=1._wp, procname=ub2b_cor )
626            IF ( Agrif_Rhoy()>1 ) CALL Agrif_Bc_variable(    vb2b_cor_id, calledweight=1._wp, procname=vb2b_cor )
627         ELSE
628            ! order matters here !!!!!!
629            CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
630            CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
631            !
632            bdy_tinterp = 1
633            CALL Agrif_Bc_variable( unb_interp_id , calledweight=1._wp, procname=interpunb  ) ! After
634            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=1._wp, procname=interpvnb  ) 
635            !
636            bdy_tinterp = 2
637            CALL Agrif_Bc_variable( unb_interp_id , calledweight=0._wp, procname=interpunb  ) ! Before
638            CALL Agrif_Bc_variable( vnb_interp_id , calledweight=0._wp, procname=interpvnb  )   
639         ENDIF
640      ELSE ! Linear interpolation
641         !
642         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
643         CALL Agrif_Bc_variable( unb_interp_id, procname=interpunb )
644         CALL Agrif_Bc_variable( vnb_interp_id, procname=interpvnb )
645      ENDIF
646      Agrif_UseSpecialValue = .FALSE.
647      use_sign_north = .FALSE.
648      !
649   END SUBROUTINE Agrif_dta_ts
650
651
652   SUBROUTINE Agrif_ssh( kt )
653      !!----------------------------------------------------------------------
654      !!                  ***  ROUTINE Agrif_ssh  ***
655      !!---------------------------------------------------------------------- 
656      INTEGER, INTENT(in) ::   kt
657      !
658      INTEGER  :: ji, jj
659      INTEGER  :: istart, iend, jstart, jend
660      !!---------------------------------------------------------------------- 
661      !
662      IF( Agrif_Root() )   RETURN
663      !     
664      ! Linear time interpolation of sea level
665      !
666      Agrif_SpecialValue    = 0._wp
667      Agrif_UseSpecialValue = .TRUE.
668      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
669      Agrif_UseSpecialValue = .FALSE.
670      !
671      ! --- West --- !
672      IF(lk_west) THEN
673         istart = nn_hls + 2                                                          ! halo + land + 1
674         iend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhox()               ! halo + land + nbghostcells
675         DO ji = mi0(istart), mi1(iend)
676            DO jj = 1, jpj
677               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
678            END DO
679         END DO
680      ENDIF
681      !
682      ! --- East --- !
683      IF(lk_east) THEN
684         istart = jpiglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhox()       ! halo + land + nbghostcells - 1
685         iend   = jpiglo - ( nn_hls + 1 )                                              ! halo + land + 1            - 1
686         DO ji = mi0(istart), mi1(iend)
687            DO jj = 1, jpj
688               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
689            END DO
690         END DO
691      ENDIF
692      !
693      ! --- South --- !
694      IF(lk_south) THEN
695         jstart = nn_hls + 2                                                          ! halo + land + 1
696         jend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhoy()               ! halo + land + nbghostcells
697         DO jj = mj0(jstart), mj1(jend)
698            DO ji = 1, jpi
699               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
700            END DO
701         END DO
702      ENDIF
703      !
704      ! --- North --- !
705      IF(lk_north) THEN
706         jstart = jpjglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhoy()     ! halo + land + nbghostcells - 1
707         jend   = jpjglo - ( nn_hls + 1 )                                            ! halo + land + 1            - 1
708         DO jj = mj0(jstart), mj1(jend)
709            DO ji = 1, jpi
710               ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
711            END DO
712         END DO
713      ENDIF
714      !
715   END SUBROUTINE Agrif_ssh
716
717
718   SUBROUTINE Agrif_ssh_ts( jn )
719      !!----------------------------------------------------------------------
720      !!                  ***  ROUTINE Agrif_ssh_ts  ***
721      !!---------------------------------------------------------------------- 
722      INTEGER, INTENT(in) ::   jn
723      !!
724      INTEGER :: ji, jj
725      INTEGER  :: istart, iend, jstart, jend
726      !!---------------------------------------------------------------------- 
727      !
728      IF( Agrif_Root() )   RETURN
729      !
730      ! --- West --- !
731      IF(lk_west) THEN
732         istart = nn_hls + 2                                                        ! halo + land + 1
733         iend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhox()             ! halo + land + nbghostcells
734         DO ji = mi0(istart), mi1(iend)
735            DO jj = 1, jpj
736               ssha_e(ji,jj) = hbdy(ji,jj)
737            END DO
738         END DO
739      ENDIF
740      !
741      ! --- East --- !
742      IF(lk_east) THEN
743         istart = jpiglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhox()    ! halo + land + nbghostcells - 1
744         iend   = jpiglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1
745         DO ji = mi0(istart), mi1(iend)
746            DO jj = 1, jpj
747               ssha_e(ji,jj) = hbdy(ji,jj)
748            END DO
749         END DO
750      ENDIF
751      !
752      ! --- South --- !
753      IF(lk_south) THEN
754         jstart = nn_hls + 2                                                        ! halo + land + 1
755         jend   = nn_hls + nbghostcells + nn_shift_bar*Agrif_Rhoy()             ! halo + land + nbghostcells
756         DO jj = mj0(jstart), mj1(jend)
757            DO ji = 1, jpi
758               ssha_e(ji,jj) = hbdy(ji,jj)
759            END DO
760         END DO
761      ENDIF
762      !
763      ! --- North --- !
764      IF(lk_north) THEN
765         jstart = jpjglo - ( nn_hls + nbghostcells -1 ) - nn_shift_bar*Agrif_Rhoy()    ! halo + land + nbghostcells - 1
766         jend   = jpjglo - ( nn_hls + 1 )                                           ! halo + land + 1            - 1
767         DO jj = mj0(jstart), mj1(jend)
768            DO ji = 1, jpi
769               ssha_e(ji,jj) = hbdy(ji,jj)
770            END DO
771         END DO
772      ENDIF
773      !
774   END SUBROUTINE Agrif_ssh_ts
775
776   
777   SUBROUTINE Agrif_avm
778      !!----------------------------------------------------------------------
779      !!                  ***  ROUTINE Agrif_avm  ***
780      !!---------------------------------------------------------------------- 
781      REAL(wp) ::   zalpha
782      !!---------------------------------------------------------------------- 
783      !
784      IF( Agrif_Root() )   RETURN
785      !
786      zalpha = 1._wp ! JC: proper time interpolation impossible 
787                     ! => use last available value from parent
788      !
789      Agrif_SpecialValue    = 0.e0
790      Agrif_UseSpecialValue = .TRUE.
791      l_vremap              = ln_vert_remap
792      !
793      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
794      !
795      Agrif_UseSpecialValue = .FALSE.
796      l_vremap              = .FALSE.
797      !
798   END SUBROUTINE Agrif_avm
799
800
801   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
802      !!----------------------------------------------------------------------
803      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
804      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
805      LOGICAL                                     , INTENT(in   ) ::   before
806      !
807      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
808      INTEGER  ::   N_in, N_out
809      INTEGER  :: item
810      ! vertical interpolation:
811      REAL(wp) :: zhtot, zwgt
812      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin, tabin_i
813      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i
814      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
815      !!----------------------------------------------------------------------
816
817      IF( before ) THEN
818
819         item = Kmm_a
820         IF( l_ini_child )   Kmm_a = Kbb_a 
821
822         DO jn = 1,jpts
823            DO jk=k1,k2
824               DO jj=j1,j2
825                 DO ji=i1,i2
826                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)
827                 END DO
828              END DO
829           END DO
830         END DO
831
832         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN
833
834            ! Fill cell depths (i.e. gdept) to be interpolated
835            ! Warning: these are masked, hence extrapolated prior interpolation.
836            DO jj=j1,j2
837               DO ji=i1,i2
838                  ptab(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a)
839                  DO jk=k1+1,k2
840                     ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
841                        & ( ptab(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) )
842                  END DO
843               END DO
844            END DO
845         
846            ! Save ssh at last level:
847            IF (.NOT.ln_linssh) THEN
848               ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
849            END IF     
850         ENDIF
851         Kmm_a = item
852
853      ELSE
854         item = Krhs_a
855         IF( l_ini_child )   Krhs_a = Kbb_a 
856
857         IF( l_vremap .OR. l_ini_child ) THEN
858            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
859            DO jj=j1,j2
860               DO ji=i1,i2
861                  ts(ji,jj,:,:,Krhs_a) = 0. 
862                  !
863                  ! Build vertical grids:
864                  N_in = mbkt_parent(ji,jj)
865                  N_out = mbkt(ji,jj)
866                  IF (N_in*N_out > 0) THEN
867                     ! Input grid (account for partial cells if any):
868                     DO jk=1,N_in
869                        z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2)
870                        tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
871                     END DO
872                 
873                     ! Intermediate grid:
874                     IF ( l_vremap ) THEN
875                        DO jk = 1, N_in
876                           h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
877                             &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
878                        END DO
879                        z_in_i(1) = 0.5_wp * h_in_i(1)
880                        DO jk=2,N_in
881                           z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
882                        END DO
883                        z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2)
884                     ENDIF                             
885
886                     ! Output (Child) grid:
887                     DO jk=1,N_out
888                        h_out(jk) = e3t(ji,jj,jk,Krhs_a)
889                     END DO
890                     z_out(1) = 0.5_wp * h_out(1)
891                     DO jk=2,N_out
892                        z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
893                     END DO
894                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)
895
896                     IF( l_ini_child ) THEN
897                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),              &
898                                      &   z_out(1:N_out),N_in,N_out,jpts) 
899                     ELSE
900                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),                       &
901                                     &   z_in_i(1:N_in),N_in,N_in,jpts)
902                        CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   &
903                                      &   h_out(1:N_out),N_in,N_out,jpts) 
904                     ENDIF
905                  ENDIF
906               END DO
907            END DO
908            Krhs_a = item
909 
910         ELSE
911         
912            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
913                                             ! linear vertical interpolation
914               DO jj=j1,j2
915                  DO ji=i1,i2
916                     !
917                     N_in  = mbkt(ji,jj)
918                     N_out = mbkt(ji,jj)
919                     z_in(1) = ptab(ji,jj,1,n2)
920                     tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts)
921                     DO jk=2, N_in
922                        z_in(jk) = ptab(ji,jj,jk,n2)
923                        tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
924                     END DO
925                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
926                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a)
927                     DO jk=2, N_out
928                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a))
929                     END DO
930                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
931                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), &
932                                   &   z_out(1:N_out),N_in,N_out,jpts) 
933                  END DO
934               END DO
935            ENDIF
936
937            DO jn =1, jpts
938               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)
939            END DO
940         ENDIF
941
942      ENDIF
943      !
944   END SUBROUTINE interptsn
945
946   
947   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
948      !!----------------------------------------------------------------------
949      !!                  ***  ROUTINE interpsshn  ***
950      !!---------------------------------------------------------------------- 
951      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
952      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
953      LOGICAL                         , INTENT(in   ) ::   before
954      !
955      !!---------------------------------------------------------------------- 
956      !
957      IF( before) THEN
958         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
959      ELSE
960         IF( l_ini_child ) THEN
961            ssh(i1:i2,j1:j2,Krhs_a) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
962         ELSE
963            hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
964         ENDIF
965      ENDIF
966      !
967   END SUBROUTINE interpsshn
968
969   
970   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
971      !!----------------------------------------------------------------------
972      !!                  *** ROUTINE interpun ***
973      !!---------------------------------------------   
974      !!
975      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
976      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
977      LOGICAL, INTENT(in) :: before
978      !!
979      INTEGER :: ji,jj,jk
980      REAL(wp) :: zrhoy, zhtot
981      ! vertical interpolation:
982      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
983      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
984      INTEGER  :: N_in, N_out,item
985      REAL(wp) :: h_diff
986      !!---------------------------------------------   
987      !
988      IF (before) THEN
989
990         item = Kmm_a
991         IF( l_ini_child )   Kmm_a = Kbb_a     
992
993         DO jk=1,jpk
994            DO jj=j1,j2
995               DO ji=i1,i2
996                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
997                  IF( l_vremap .OR. l_ini_child) THEN
998                     ! Interpolate thicknesses (masked for subsequent extrapolation)
999                     ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
1000                  ENDIF
1001               END DO
1002            END DO
1003         END DO
1004
1005        IF( l_vremap .OR. l_ini_child ) THEN
1006         ! Extrapolate thicknesses in partial bottom cells:
1007         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1008            IF (ln_zps) THEN
1009               DO jj=j1,j2
1010                  DO ji=i1,i2
1011                     jk = mbku(ji,jj)
1012                     ptab(ji,jj,jk,2) = 0._wp
1013                  END DO
1014               END DO           
1015            END IF
1016
1017           ! Save ssh at last level:
1018           ptab(i1:i2,j1:j2,k2,2) = 0._wp
1019           IF (.NOT.ln_linssh) THEN
1020              ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
1021              DO jk=1,jpk
1022                 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)
1023              END DO
1024              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
1025           END IF
1026        ENDIF
1027
1028         Kmm_a = item
1029         !
1030      ELSE
1031         zrhoy = Agrif_rhoy()
1032
1033        IF( l_vremap .OR. l_ini_child) THEN
1034! VERTICAL REFINEMENT BEGIN
1035
1036            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1037
1038            DO ji=i1,i2
1039               DO jj=j1,j2
1040                  uu(ji,jj,:,Krhs_a) = 0._wp
1041                  N_in = mbku_parent(ji,jj)
1042                  N_out = mbku(ji,jj)
1043                  IF (N_in*N_out > 0) THEN
1044                     zhtot = 0._wp
1045                     DO jk=1,N_in
1046                        !IF (jk==N_in) THEN
1047                        !   h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1048                        !ELSE
1049                        !   h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)
1050                        !ENDIF
1051                        IF ( l_vremap ) THEN
1052                           h_in(jk) = e3u0_parent(ji,jj,jk) 
1053                        ELSE
1054                           IF (jk==N_in) THEN
1055                              h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1056                           ELSE
1057                              h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
1058                           ENDIF
1059                        ENDIF
1060                        zhtot = zhtot + h_in(jk)
1061                        IF( h_in(jk) .GT. 0. ) THEN
1062                           tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
1063                        ELSE
1064                           tabin(jk) = 0.
1065                        ENDIF
1066                    END DO
1067                    z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj) 
1068                    DO jk=2,N_in
1069                       z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk)+h_in(jk-1))
1070                    END DO
1071                     
1072                    DO jk=1, N_out
1073                       h_out(jk) = e3u(ji,jj,jk,Krhs_a)
1074                    END DO
1075
1076                    z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj)
1077                    DO jk=2,N_out
1078                       z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1) + h_out(jk)) 
1079                    END DO 
1080
1081                    IF( l_ini_child ) THEN
1082                       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)
1083                    ELSE
1084                       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)
1085                    ENDIF   
1086                 ENDIF
1087               END DO
1088            END DO
1089         ELSE
1090            DO jk = 1, jpkm1
1091               DO jj=j1,j2
1092                  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) )
1093               END DO
1094            END DO
1095         ENDIF
1096
1097      ENDIF
1098      !
1099   END SUBROUTINE interpun
1100
1101   
1102   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1103      !!----------------------------------------------------------------------
1104      !!                  *** ROUTINE interpvn ***
1105      !!----------------------------------------------------------------------
1106      !
1107      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
1108      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1109      LOGICAL, INTENT(in) :: before
1110      !
1111      INTEGER :: ji,jj,jk
1112      REAL(wp) :: zrhox
1113      ! vertical interpolation:
1114      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
1115      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
1116      INTEGER  :: N_in, N_out, item
1117      REAL(wp) :: h_diff, zhtot
1118      !!---------------------------------------------   
1119      !     
1120      IF (before) THEN   
1121
1122         item = Kmm_a
1123         IF( l_ini_child )   Kmm_a = Kbb_a     
1124       
1125         DO jk=k1,k2
1126            DO jj=j1,j2
1127               DO ji=i1,i2
1128                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
1129                  IF( l_vremap .OR. l_ini_child) THEN
1130                     ! Interpolate thicknesses (masked for subsequent extrapolation)
1131                     ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
1132                  ENDIF
1133               END DO
1134            END DO
1135         END DO
1136
1137         IF( l_vremap .OR. l_ini_child) THEN
1138         ! Extrapolate thicknesses in partial bottom cells:
1139         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1140            IF (ln_zps) THEN
1141               DO jj=j1,j2
1142                  DO ji=i1,i2
1143                     jk = mbkv(ji,jj)
1144                     ptab(ji,jj,jk,2) = 0._wp
1145                  END DO
1146               END DO           
1147            END IF
1148            ! Save ssh at last level:
1149            ptab(i1:i2,j1:j2,k2,2) = 0._wp
1150            IF (.NOT.ln_linssh) THEN
1151               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
1152               DO jk=1,jpk
1153                  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)
1154               END DO
1155               ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
1156            END IF
1157         ENDIF
1158         item = Kmm_a
1159
1160      ELSE       
1161         zrhox = Agrif_rhox()
1162
1163         IF( l_vremap .OR. l_ini_child ) THEN
1164
1165            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1166
1167            DO jj=j1,j2
1168               DO ji=i1,i2
1169                  vv(ji,jj,:,Krhs_a) = 0._wp
1170                  N_in = mbkv_parent(ji,jj)
1171                  N_out = mbkv(ji,jj)
1172
1173                  IF (N_in*N_out > 0) THEN
1174                     zhtot = 0._wp
1175                     DO jk=1,N_in
1176                        !IF (jk==N_in) THEN
1177                        !   h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1178                        !ELSE
1179                        !   h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1180                        !ENDIF
1181                        IF (l_vremap) THEN
1182                           h_in(jk) = e3v0_parent(ji,jj,jk)
1183                        ELSE
1184                           IF (jk==N_in) THEN
1185                              h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
1186                           ELSE
1187                              h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
1188                           ENDIF
1189                        ENDIF
1190                        zhtot = zhtot + h_in(jk)
1191                        IF( h_in(jk) .GT. 0. ) THEN
1192                          tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
1193                        ELSE
1194                          tabin(jk)  = 0.
1195                        ENDIF
1196                     END DO
1197
1198                     z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj)
1199                     DO jk=2,N_in
1200                        z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk-1)+h_in(jk))
1201                     END DO
1202
1203                     DO jk=1,N_out
1204                        h_out(jk) = e3v(ji,jj,jk,Krhs_a)
1205                     END DO
1206
1207                     z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj)
1208                     DO jk=2,N_out
1209                        z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1)+h_out(jk))
1210                     END DO
1211 
1212                     IF( l_ini_child ) THEN
1213                        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)
1214                     ELSE
1215                        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)
1216                     ENDIF   
1217                  ENDIF
1218               END DO
1219            END DO
1220         ELSE
1221            DO jk = 1, jpkm1
1222               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) )
1223            END DO
1224         ENDIF
1225      ENDIF
1226      !       
1227   END SUBROUTINE interpvn
1228
1229   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
1230      !!----------------------------------------------------------------------
1231      !!                  ***  ROUTINE interpunb  ***
1232      !!---------------------------------------------------------------------- 
1233      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1234      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1235      LOGICAL                         , INTENT(in   ) ::   before
1236      !
1237      INTEGER  ::   ji, jj
1238      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1239      !!---------------------------------------------------------------------- 
1240      !
1241      IF( before ) THEN
1242         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)
1243      ELSE
1244         zrhoy = Agrif_Rhoy()
1245         zrhot = Agrif_rhot()
1246         ! Time indexes bounds for integration
1247         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1248         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1249         !
1250         DO ji = i1, i2
1251            DO jj = j1, j2
1252               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1253                  IF    ( utint_stage(ji,jj) == 1  ) THEN
1254                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1255                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1256                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
1257                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1258                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1259                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
1260                     ztcoeff = 1._wp
1261                  ELSE
1262                     ztcoeff = 0._wp
1263                  ENDIF
1264                  !   
1265                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
1266                  !           
1267                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
1268                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
1269                  ENDIF
1270                  !
1271                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
1272               ENDIF
1273            END DO
1274         END DO
1275      END IF
1276      !
1277   END SUBROUTINE interpunb
1278
1279
1280   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
1281      !!----------------------------------------------------------------------
1282      !!                  ***  ROUTINE interpvnb  ***
1283      !!---------------------------------------------------------------------- 
1284      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1285      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1286      LOGICAL                         , INTENT(in   ) ::   before
1287      !
1288      INTEGER  ::   ji, jj
1289      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1290      !!---------------------------------------------------------------------- 
1291      !
1292      IF( before ) THEN
1293         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)
1294      ELSE
1295         zrhox = Agrif_Rhox()
1296         zrhot = Agrif_rhot()
1297         ! Time indexes bounds for integration
1298         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1299         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
1300         !     
1301         DO ji = i1, i2
1302            DO jj = j1, j2
1303               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1304                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
1305                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1306                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1307                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
1308                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1309                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1310                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
1311                     ztcoeff = 1._wp
1312                  ELSE
1313                     ztcoeff = 0._wp
1314                  ENDIF
1315                  !   
1316                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
1317                  !           
1318                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
1319                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
1320                  ENDIF
1321                  !
1322                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
1323               ENDIF
1324            END DO
1325         END DO         
1326      ENDIF
1327      !
1328   END SUBROUTINE interpvnb
1329
1330
1331   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
1332      !!----------------------------------------------------------------------
1333      !!                  ***  ROUTINE interpub2b  ***
1334      !!---------------------------------------------------------------------- 
1335      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1336      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1337      LOGICAL                         , INTENT(in   ) ::   before
1338      !
1339      INTEGER  ::   ji,jj
1340      REAL(wp) ::   zrhot, zt0, zt1, zat
1341      !!---------------------------------------------------------------------- 
1342      IF( before ) THEN
1343!         IF ( ln_bt_fw ) THEN
1344            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1345!         ELSE
1346!            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1347!         ENDIF
1348      ELSE
1349         zrhot = Agrif_rhot()
1350         ! Time indexes bounds for integration
1351         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1352         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1353         ! Polynomial interpolation coefficients:
1354         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1355            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1356         !
1357         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1358         !
1359         ! Update interpolation stage:
1360         utint_stage(i1:i2,j1:j2) = 1
1361      ENDIF
1362      !
1363   END SUBROUTINE interpub2b
1364   
1365   SUBROUTINE interpub2b_const( ptab, i1, i2, j1, j2, before )
1366      !!----------------------------------------------------------------------
1367      !!                  ***  ROUTINE interpub2b_const  ***
1368      !!---------------------------------------------------------------------- 
1369      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1370      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1371      LOGICAL                         , INTENT(in   ) ::   before
1372      !
1373      REAL(wp) :: zrhoy
1374      !!---------------------------------------------------------------------- 
1375      IF( before ) THEN
1376!         IF ( ln_bt_fw ) THEN
1377            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1378!         ELSE
1379!            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1380!         ENDIF
1381      ELSE
1382         zrhoy = Agrif_Rhoy()
1383         !
1384         ubdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) & 
1385                           & / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
1386         !
1387      ENDIF
1388      !
1389   END SUBROUTINE interpub2b_const
1390
1391
1392   SUBROUTINE ub2b_cor( ptab, i1, i2, j1, j2, before )
1393      !!----------------------------------------------------------------------
1394      !!                  ***  ROUTINE ub2b_cor  ***
1395      !!---------------------------------------------------------------------- 
1396      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1397      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1398      LOGICAL                         , INTENT(in   ) ::   before
1399      !
1400      INTEGER  :: ji, jj
1401      INTEGER  :: imin, imax, jmin, jmax
1402      REAL(wp) :: zrhox, zrhoy, zx
1403      !!---------------------------------------------------------------------- 
1404      IF( before ) THEN
1405         ptab(:,:) = 0._wp
1406         imin = MAX(i1, 2) ; imax = MIN(i2, jpi-1)
1407         jmin = MAX(j1, 2) ; jmax = MIN(j2, jpj-1)
1408         DO ji=imin,imax
1409            DO jj=jmin,jmax
1410               ptab(ji,jj) = 0.25_wp*( ( vb2_b(ji+1,jj  )*e1v(ji+1,jj  )   & 
1411                           &            -vb2_b(ji-1,jj  )*e1v(ji-1,jj  ) ) &
1412                           &          -( vb2_b(ji+1,jj-1)*e1v(ji+1,jj-1)   &
1413                           &            -vb2_b(ji-1,jj-1)*e1v(ji-1,jj-1) ) )
1414            END DO
1415         END DO
1416      ELSE
1417         !
1418         zrhox = Agrif_Rhox() 
1419         zrhoy = Agrif_Rhoy()
1420         DO ji=i1,i2
1421            DO jj=j1,j2
1422               IF (utint_stage(ji,jj)==0) THEN
1423                  zx = 2._wp*MOD(ABS(mig0(ji)-nbghostcells_x_w), INT(Agrif_Rhox()))/zrhox - 1._wp 
1424                  ubdy(ji,jj) = ubdy(ji,jj) + 0.25_wp*(1._wp-zx*zx) * ptab(ji,jj) & 
1425                              &         / zrhoy *r1_e2u(ji,jj) * umask(ji,jj,1) 
1426                  utint_stage(ji,jj) = 1 
1427               ENDIF
1428            END DO
1429         END DO 
1430         !
1431      ENDIF
1432      !
1433   END SUBROUTINE ub2b_cor
1434
1435
1436   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
1437      !!----------------------------------------------------------------------
1438      !!                  ***  ROUTINE interpvb2b  ***
1439      !!---------------------------------------------------------------------- 
1440      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1441      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1442      LOGICAL                         , INTENT(in   ) ::   before
1443      !
1444      INTEGER ::   ji,jj
1445      REAL(wp) ::   zrhot, zt0, zt1, zat
1446      !!---------------------------------------------------------------------- 
1447      !
1448      IF( before ) THEN
1449!         IF ( ln_bt_fw ) THEN
1450            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1451!         ELSE
1452!            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1453!         ENDIF
1454      ELSE     
1455         zrhot = Agrif_rhot()
1456         ! Time indexes bounds for integration
1457         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1458         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1459         ! Polynomial interpolation coefficients:
1460         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1461            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1462         !
1463         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1464         !
1465         ! update interpolation stage:
1466         vtint_stage(i1:i2,j1:j2) = 1
1467      ENDIF
1468      !     
1469   END SUBROUTINE interpvb2b
1470
1471
1472   SUBROUTINE interpvb2b_const( ptab, i1, i2, j1, j2, before )
1473      !!----------------------------------------------------------------------
1474      !!                  ***  ROUTINE interpub2b_const  ***
1475      !!---------------------------------------------------------------------- 
1476      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1477      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1478      LOGICAL                         , INTENT(in   ) ::   before
1479      !
1480      REAL(wp) :: zrhox
1481      !!---------------------------------------------------------------------- 
1482      IF( before ) THEN
1483!         IF ( ln_bt_fw ) THEN
1484            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1485!         ELSE
1486!            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1487!         ENDIF
1488      ELSE
1489         zrhox = Agrif_Rhox()
1490         !
1491         vbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) &
1492                           & / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
1493         !
1494      ENDIF
1495      !
1496   END SUBROUTINE interpvb2b_const
1497
1498 
1499   SUBROUTINE vb2b_cor( ptab, i1, i2, j1, j2, before )
1500      !!----------------------------------------------------------------------
1501      !!                  ***  ROUTINE vb2b_cor  ***
1502      !!---------------------------------------------------------------------- 
1503      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1504      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1505      LOGICAL                         , INTENT(in   ) ::   before
1506      !
1507      INTEGER  :: ji, jj
1508      INTEGER  :: imin, imax, jmin, jmax
1509      REAL(wp) :: zrhox, zrhoy, zy
1510      !!---------------------------------------------------------------------- 
1511      IF( before ) THEN
1512         ptab(:,:) = 0._wp
1513         imin = MAX(i1, 2) ; imax = MIN(i2, jpi-1)
1514         jmin = MAX(j1, 2) ; jmax = MIN(j2, jpj-1)
1515         DO ji=imin,imax
1516            DO jj=jmin,jmax
1517               ptab(ji,jj) = 0.25_wp*( ( ub2_b(ji  ,jj+1)*e2u(ji  ,jj+1)   & 
1518                           &            -ub2_b(ji  ,jj-1)*e2u(ji  ,jj-1) ) &
1519                           &          -( ub2_b(ji-1,jj+1)*e2u(ji-1,jj+1)   &
1520                           &            -ub2_b(ji-1,jj-1)*e2u(ji-1,jj-1) ) )
1521            END DO
1522         END DO
1523      ELSE
1524         !
1525         zrhox = Agrif_Rhox() 
1526         zrhoy = Agrif_Rhoy()
1527         DO ji=i1,i2
1528            DO jj=j1,j2
1529               IF (vtint_stage(ji,jj)==0) THEN
1530                  zy = 2._wp*MOD(ABS(mjg0(jj)-nbghostcells_y_s), INT(Agrif_Rhoy()))/zrhoy - 1._wp 
1531                  vbdy(ji,jj) = vbdy(ji,jj) + 0.25_wp*(1._wp-zy*zy) * ptab(ji,jj) & 
1532                              &         / zrhox * r1_e1v(ji,jj) * vmask(ji,jj,1) 
1533                  vtint_stage(ji,jj) = 1 
1534               ENDIF
1535            END DO
1536         END DO 
1537         !
1538      ENDIF
1539      !
1540   END SUBROUTINE vb2b_cor
1541
1542
1543   SUBROUTINE interpe3t0_vremap( ptab, i1, i2, j1, j2, k1, k2, before )
1544      !!----------------------------------------------------------------------
1545      !!                  ***  ROUTINE interpe3t0_vremap  ***
1546      !!---------------------------------------------------------------------- 
1547      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1548      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1549      LOGICAL                              , INTENT(in   ) :: before
1550      !
1551      INTEGER :: ji, jj, jk
1552      REAL(wp) :: zh
1553      !!---------------------------------------------------------------------- 
1554      !   
1555      IF( before ) THEN
1556         IF ( ln_zps ) THEN
1557            DO jk = k1, k2
1558               DO jj = j1, j2
1559                  DO ji = i1, i2
1560                     ptab(ji, jj, jk) = e3t_1d(jk)
1561                  END DO
1562               END DO
1563            END DO
1564         ELSE
1565            ptab(i1:i2,j1:j2,k1:k2) = e3t_0(i1:i2,j1:j2,k1:k2)
1566         ENDIF
1567      ELSE
1568         !
1569         DO jk = k1, k2
1570            DO jj = j1, j2
1571               DO ji = i1, i2
1572                  e3t0_parent(ji,jj,jk) = ptab(ji,jj,jk)
1573               END DO
1574            END DO
1575         END DO
1576
1577         ! Retrieve correct scale factor at the bottom:
1578         DO jj = j1, j2
1579            DO ji = i1, i2
1580               zh = 0._wp
1581               DO jk = 1, mbkt_parent(ji, jj)-1
1582                  zh = zh + e3t0_parent(ji,jj,jk)
1583               END DO
1584               e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh
1585            END DO
1586         END DO
1587         
1588      ENDIF
1589      !
1590   END SUBROUTINE interpe3t0_vremap
1591
1592
1593   SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before )
1594      !!----------------------------------------------------------------------
1595      !!                  ***  ROUTINE interpglamt  ***
1596      !!---------------------------------------------------------------------- 
1597      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
1598      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1599      LOGICAL                        , INTENT(in   ) :: before
1600      !
1601      INTEGER :: ji, jj, jk
1602      REAL(wp):: ztst
1603      !!---------------------------------------------------------------------- 
1604      !   
1605      IF( before ) THEN
1606         ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2)
1607      ELSE
1608         ztst = MAXVAL(ABS(glamt(i1:i2,j1:j2)))*1.e-4
1609         DO jj = j1, j2
1610            DO ji = i1, i2
1611               IF( ABS( ptab(ji,jj) - glamt(ji,jj) ) > ztst ) THEN
1612                  WRITE(numout,*) ' Agrif error for glamt: parent, child, i, j ', ptab(ji,jj), glamt(ji,jj), mig0(ji), mig0(jj)
1613!                  kindic_agr = kindic_agr + 1
1614               ENDIF
1615            END DO
1616         END DO
1617      ENDIF
1618      !
1619   END SUBROUTINE interpglamt
1620
1621
1622   SUBROUTINE interpgphit( ptab, i1, i2, j1, j2, before )
1623      !!----------------------------------------------------------------------
1624      !!                  ***  ROUTINE interpgphit  ***
1625      !!---------------------------------------------------------------------- 
1626      INTEGER                        , INTENT(in   ) :: i1, i2, j1, j2
1627      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1628      LOGICAL                        , INTENT(in   ) :: before
1629      !
1630      INTEGER :: ji, jj, jk
1631      REAL(wp):: ztst
1632      !!---------------------------------------------------------------------- 
1633      !   
1634      IF( before ) THEN
1635         ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2)
1636      ELSE
1637         ztst = MAXVAL(ABS(gphit(i1:i2,j1:j2)))*1.e-4
1638         DO jj = j1, j2
1639            DO ji = i1, i2
1640               IF( ABS( ptab(ji,jj) - gphit(ji,jj) ) > ztst ) THEN
1641                  WRITE(numout,*) ' Agrif error for gphit: parent, child, i, j ', ptab(ji,jj), gphit(ji,jj), mig0(ji), mig0(jj)
1642!                  kindic_agr = kindic_agr + 1
1643               ENDIF
1644            END DO
1645         END DO
1646      ENDIF
1647      !
1648   END SUBROUTINE interpgphit
1649
1650
1651   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1652      !!----------------------------------------------------------------------
1653      !!                  ***  ROUTINE interavm  ***
1654      !!---------------------------------------------------------------------- 
1655      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1656      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1657      LOGICAL                                    , INTENT(in   ) ::   before
1658      !
1659      INTEGER  :: ji, jj, jk
1660      INTEGER  :: N_in, N_out
1661      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in
1662      REAL(wp), DIMENSION(1:jpk) :: z_out
1663      !!---------------------------------------------------------------------- 
1664      !     
1665      IF (before) THEN         
1666         DO jk=k1,k2
1667            DO jj=j1,j2
1668              DO ji=i1,i2
1669                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1670              END DO
1671           END DO
1672         END DO
1673
1674         IF( l_vremap ) THEN
1675            ! Interpolate interfaces
1676            ! Warning: these are masked, hence extrapolated prior interpolation.
1677            DO jk=k1,k2
1678               DO jj=j1,j2
1679                  DO ji=i1,i2
1680                      ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * gdepw(ji,jj,jk,Kmm_a)
1681                  END DO
1682               END DO
1683            END DO
1684       
1685           ! Save ssh at last level:
1686            IF (.NOT.ln_linssh) THEN
1687               ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
1688            ELSE
1689               ptab(i1:i2,j1:j2,k2,2) = 0._wp
1690            END IF     
1691          ENDIF
1692
1693      ELSE
1694
1695         IF( l_vremap ) THEN
1696            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1697            avm_k(i1:i2,j1:j2,1:jpkm1) = 0._wp
1698               
1699            DO jj = j1, j2
1700               DO ji =i1, i2
1701                  N_in = mbkt_parent(ji,jj)
1702                  N_out = mbkt(ji,jj)
1703                  IF (N_in*N_out > 0) THEN
1704                     DO jk = 1, N_in  ! Parent vertical grid               
1705                        z_in(jk)  = ptab(ji,jj,jk,2) - ptab(ji,jj,k2,2)
1706                        tabin(jk) = ptab(ji,jj,jk,1)
1707                     END DO
1708                     DO jk = 1, N_out        ! Child vertical grid
1709                        z_out(jk) = gdepw(ji,jj,jk,Kmm_a) - ssh(ji,jj,Kmm_a)
1710                     END DO
1711                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kmm_a)
1712
1713                     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)
1714                  ENDIF
1715               END DO
1716            END DO
1717         ELSE
1718            avm_k(i1:i2,j1:j2,1:jpkm1) = ptab (i1:i2,j1:j2,1:jpkm1,1)
1719         ENDIF
1720      ENDIF
1721      !
1722   END SUBROUTINE interpavm
1723
1724   
1725   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
1726      !!----------------------------------------------------------------------
1727      !!                  ***  ROUTINE interpmbkt  ***
1728      !!---------------------------------------------------------------------- 
1729      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1730      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1731      LOGICAL                         , INTENT(in   ) ::   before
1732      !
1733      !!---------------------------------------------------------------------- 
1734      !
1735      IF( before) THEN
1736         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
1737      ELSE
1738         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
1739      ENDIF
1740      !
1741   END SUBROUTINE interpmbkt
1742
1743   
1744   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
1745      !!----------------------------------------------------------------------
1746      !!                  ***  ROUTINE interpht0  ***
1747      !!---------------------------------------------------------------------- 
1748      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1749      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1750      LOGICAL                         , INTENT(in   ) ::   before
1751      !
1752      !!---------------------------------------------------------------------- 
1753      !
1754      IF( before) THEN
1755         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2)
1756      ELSE
1757         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
1758      ENDIF
1759      !
1760   END SUBROUTINE interpht0
1761
1762   SUBROUTINE Agrif_check_bat( iindic )
1763      !!----------------------------------------------------------------------
1764      !!                  ***  ROUTINE Agrif_check_bat  ***
1765      !!---------------------------------------------------------------------- 
1766      INTEGER, INTENT(inout) ::   iindic
1767      !!
1768      INTEGER :: ji, jj, jk
1769      INTEGER  :: istart, iend, jstart, jend, ispon
1770      !!---------------------------------------------------------------------- 
1771      !
1772      !
1773      ! --- West --- !
1774      IF(lk_west) THEN
1775         ispon  = nn_sponge_len * Agrif_irhox()
1776         istart = nn_hls + 2                                  ! halo + land + 1
1777         iend   = nn_hls + nbghostcells + ispon           ! halo + land + nbghostcells + sponge
1778         jstart = nn_hls + 2 
1779         jend   = jpjglo - nn_hls - 1 
1780         DO ji = mi0(istart), mi1(iend)
1781            DO jj = mj0(jstart), mj1(jend)
1782               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1783               IF ( .NOT.ln_vert_remap) THEN
1784                  DO jk = 1, jpkm1
1785                     IF ( ABS(e3t0_parent(ji,jj,jk)-e3t_0(ji,jj,jk))*tmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1786                  END DO
1787               ENDIF
1788            END DO
1789            DO jj = mj0(jstart), mj1(jend-1)
1790               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1791               IF ( .NOT.ln_vert_remap) THEN
1792                  DO jk = 1, jpkm1
1793                     IF ( ABS(e3v0_parent(ji,jj,jk)-e3v_0(ji,jj,jk))*vmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1794                  END DO
1795               ENDIF
1796            END DO
1797         END DO
1798         DO ji = mi0(istart), mi1(iend-1)
1799            DO jj = mj0(jstart), mj1(jend)
1800               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1801               IF ( .NOT.ln_vert_remap) THEN
1802                  DO jk = 1, jpkm1
1803                     IF ( ABS(e3u0_parent(ji,jj,jk)-e3u_0(ji,jj,jk))*umask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1804                  END DO
1805               ENDIF
1806            END DO
1807         END DO
1808      ENDIF
1809      !
1810      ! --- East --- !
1811      IF(lk_east) THEN
1812         ispon  = nn_sponge_len * Agrif_irhox() 
1813         istart = jpiglo - ( nn_hls + nbghostcells + ispon -1 )  ! halo + land + nbghostcells + sponge - 1
1814         iend   = jpiglo - nn_hls - 1                            ! halo + land + 1                     - 1
1815         jstart = nn_hls + 2 
1816         jend   = jpjglo - nn_hls - 1
1817         DO ji = mi0(istart), mi1(iend)
1818            DO jj = mj0(jstart), mj1(jend)
1819               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1820               IF ( .NOT.ln_vert_remap) THEN
1821                  DO jk = 1, jpkm1
1822                     IF ( ABS(e3t0_parent(ji,jj,jk)-e3t_0(ji,jj,jk))*tmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1823                  END DO
1824               ENDIF
1825            END DO
1826            DO jj = mj0(jstart), mj1(jend-1)
1827               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1828               IF ( .NOT.ln_vert_remap) THEN
1829                  DO jk = 1, jpkm1
1830                     IF ( ABS(e3v0_parent(ji,jj,jk)-e3v_0(ji,jj,jk))*vmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1831                  END DO
1832               ENDIF
1833            END DO
1834         END DO
1835         DO ji = mi0(istart), mi1(iend-1)
1836            DO jj = mj0(jstart), mj1(jend)
1837               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1838               IF ( .NOT.ln_vert_remap) THEN
1839                  DO jk = 1, jpkm1
1840                     IF ( ABS(e3u0_parent(ji,jj,jk)-e3u_0(ji,jj,jk))*umask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1841                  END DO
1842               ENDIF
1843            END DO
1844         END DO
1845      ENDIF
1846      !
1847      ! --- South --- !
1848      IF(lk_south) THEN
1849         ispon  = nn_sponge_len * Agrif_irhoy() 
1850         jstart = nn_hls + 2                                 ! halo + land + 1
1851         jend   = nn_hls + nbghostcells + ispon          ! halo + land + nbghostcells + sponge
1852         istart = nn_hls + 2 
1853         iend   = jpiglo - nn_hls - 1 
1854         DO jj = mj0(jstart), mj1(jend)
1855            DO ji = mi0(istart), mi1(iend)
1856               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1857               IF ( .NOT.ln_vert_remap) THEN
1858                  DO jk = 1, jpkm1
1859                     IF ( ABS(e3t0_parent(ji,jj,jk)-e3t_0(ji,jj,jk))*tmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1860                  END DO
1861               ENDIF
1862            END DO
1863            DO ji = mi0(istart), mi1(iend-1)
1864               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1865               IF ( .NOT.ln_vert_remap) THEN
1866                  DO jk = 1, jpkm1
1867                     IF ( ABS(e3u0_parent(ji,jj,jk)-e3u_0(ji,jj,jk))*umask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1868                  END DO
1869               ENDIF
1870            END DO
1871         END DO
1872         DO jj = mj0(jstart), mj1(jend-1)
1873            DO ji = mi0(istart), mi1(iend)
1874               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1875               IF ( .NOT.ln_vert_remap) THEN
1876                  DO jk = 1, jpkm1
1877                     IF ( ABS(e3v0_parent(ji,jj,jk)-e3v_0(ji,jj,jk))*vmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1878                  END DO
1879               ENDIF
1880            END DO
1881         END DO
1882      ENDIF
1883      !
1884      ! --- North --- !
1885      IF(lk_north) THEN
1886         ispon  = nn_sponge_len * Agrif_irhoy() 
1887         jstart = jpjglo - ( nn_hls + nbghostcells + ispon - 1)  ! halo + land + nbghostcells +sponge - 1
1888         jend   = jpjglo - nn_hls - 1                            ! halo + land + 1            - 1
1889         istart = nn_hls + 2 
1890         iend   = jpiglo - nn_hls - 1 
1891         DO jj = mj0(jstart), mj1(jend)
1892            DO ji = mi0(istart), mi1(iend)
1893               IF ( ABS(ht0_parent(ji,jj)-ht_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1894               IF ( .NOT.ln_vert_remap) THEN
1895                  DO jk = 1, jpkm1
1896                     IF ( ABS(e3t0_parent(ji,jj,jk)-e3t_0(ji,jj,jk))*tmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1897                  END DO
1898               ENDIF
1899            END DO
1900            DO ji = mi0(istart), mi1(iend-1)
1901               IF ( ABS(hu0_parent(ji,jj)-hu_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1902               IF ( .NOT.ln_vert_remap) THEN
1903                  DO jk = 1, jpkm1
1904                     IF ( ABS(e3u0_parent(ji,jj,jk)-e3u_0(ji,jj,jk))*umask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1905                  END DO
1906               ENDIF
1907            END DO
1908         END DO
1909         DO jj = mj0(jstart), mj1(jend-1)
1910            DO ji = mi0(istart), mi1(iend)
1911               IF ( ABS(hv0_parent(ji,jj)-hv_0(ji,jj)) > 1.e-3 ) iindic = iindic + 1
1912               IF ( .NOT.ln_vert_remap) THEN
1913                  DO jk = 1, jpkm1
1914                     IF ( ABS(e3v0_parent(ji,jj,jk)-e3v_0(ji,jj,jk))*vmask(ji,jj,jk) > 1.e-3 ) iindic = iindic + 1
1915                  END DO
1916               ENDIF
1917            END DO
1918         END DO
1919      ENDIF
1920      !
1921   END SUBROUTINE Agrif_check_bat
1922   
1923#else
1924   !!----------------------------------------------------------------------
1925   !!   Empty module                                          no AGRIF zoom
1926   !!----------------------------------------------------------------------
1927CONTAINS
1928   SUBROUTINE Agrif_OCE_Interp_empty
1929      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1930   END SUBROUTINE Agrif_OCE_Interp_empty
1931#endif
1932
1933   !!======================================================================
1934END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.