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

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

Last change on this file since 13237 was 13216, checked in by rblod, 4 years ago

Merge dev_r12973_AGRIF_CMEMS

  • 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.