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_opa_interp.F90 in branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 7824

Last change on this file since 7824 was 7824, checked in by timgraham, 7 years ago

Commit changes to interp & update routines for U, V, T&S.
This version now works correctly for simple GYRE test case with 41 levels on the child grid.
  • Property svn:keywords set to Id
File size: 57.2 KB
Line 
1MODULE agrif_opa_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_opa_interp  ***
4   !! AGRIF: interpolation package
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (XXX)  Original cade
7   !!             -   !  2005-11  (XXX)
8   !!            3.2  !  2009-04  (R. Benshila)
9   !!            3.6  !  2014-09  (R. Benshila)
10   !!----------------------------------------------------------------------
11#if defined key_agrif && ! defined key_offline
12   !!----------------------------------------------------------------------
13   !!   'key_agrif'                                              AGRIF zoom
14   !!   NOT 'key_offline'                               NO off-line tracers
15   !!----------------------------------------------------------------------
16   !!   Agrif_tra     :
17   !!   Agrif_dyn     :
18   !!   interpu       :
19   !!   interpv       :
20   !!----------------------------------------------------------------------
21   USE par_oce
22   USE oce
23   USE dom_oce     
24   USE zdf_oce
25   USE agrif_oce
26   USE phycst
27   !
28   USE in_out_manager
29   USE agrif_opa_sponge
30   USE lib_mpp
31   USE wrk_nemo
32 
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
37! VERTICAL REFINEMENT BEGIN   
38   PUBLIC   Agrif_Init_InterpScales
39! VERTICAL REFINEMENT END
40   PUBLIC   interpun, interpvn
41   PUBLIC   interptsn,  interpsshn
42   PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b
43   PUBLIC   interpe3t, interpumsk, interpvmsk
44# if defined key_zdftke
45   PUBLIC   Agrif_tke, interpavm
46# endif
47
48   INTEGER ::   bdy_tinterp = 0
49
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/NST 3.7 , NEMO Consortium (2015)
53   !! $Id$
54   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56! VERTICAL REFINEMENT BEGIN
57   REAL, DIMENSION(:,:,:), ALLOCATABLE :: interp_scales_t, interp_scales_u, interp_scales_v
58!$AGRIF_DO_NOT_TREAT
59   LOGICAL :: scaleT, scaleU, scaleV = .FALSE.
60!$AGRIF_END_DO_NOT_TREAT
61! VERTICAL REFINEMENT END
62
63CONTAINS
64
65! VERTICAL REFINEMENT BEGIN
66
67   SUBROUTINE Agrif_Init_InterpScales
68
69    scaleT = .TRUE.
70    Call Agrif_Bc_Variable(scales_t_id,calledweight=1.,procname=interpscales)
71    scaleT = .FALSE.
72   
73    scaleU = .TRUE.
74    Call Agrif_Bc_Variable(scales_u_id,calledweight=1.,procname=interpscales)
75    scaleU = .FALSE.
76
77    scaleV = .TRUE.
78    Call Agrif_Bc_Variable(scales_v_id,calledweight=1.,procname=interpscales)
79    scaleV = .FALSE.
80
81   END SUBROUTINE Agrif_Init_InterpScales
82   
83   SUBROUTINE interpscales(ptab,i1,i2,j1,j2,k1,k2,before)
84      !!---------------------------------------------
85      !!   *** ROUTINE interpscales ***
86      !!---------------------------------------------
87     
88      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
89      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
90
91      INTEGER :: ji, jj, jk
92      LOGICAL :: before
93
94      IF (before) THEN
95      IF (scaleT ) THEN
96      DO jk=k1,k2
97         DO jj=j1,j2
98            DO ji=i1,i2
99!               ptab(ji,jj,jk) = e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
100               ptab(ji,jj,jk) = e3t_n(ji,jj,jk)
101            END DO
102         END DO
103      END DO
104      ELSEIF (scaleU) THEN
105      DO jk=k1,k2
106         DO jj=j1,j2
107            DO ji=i1,i2
108!               ptab(ji,jj,jk) = e3u_n(ji,jj,jk) * umask(ji,jj,jk)
109!               ptab(ji,jj,jk) = e3u_n(ji,jj,jk)
110                ptab(ji,jj,jk) = umask(ji,jj,jk)
111            END DO
112         END DO
113      END DO
114      ELSEIF (scaleV) THEN
115      DO jk=k1,k2
116         DO jj=j1,j2
117            DO ji=i1,i2
118!               ptab(ji,jj,jk) = e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
119!               ptab(ji,jj,jk) = e3v_n(ji,jj,jk)
120               ptab(ji,jj,jk) = vmask(ji,jj,jk)
121            END DO
122         END DO
123      END DO
124      ENDIF
125      ELSE
126      IF (scaleT ) THEN
127      IF (.not.allocated(interp_scales_t)) allocate(interp_scales_t(jpi,jpj,k1:k2))
128      DO jk=k1,k2
129         DO jj=j1,j2
130            DO ji=i1,i2
131               interp_scales_t(ji,jj,jk) = ptab(ji,jj,jk)
132            END DO
133         END DO
134      END DO
135      ELSEIF (scaleU) THEN
136      IF (.not.allocated(interp_scales_u)) allocate(interp_scales_u(jpi,jpj,k1:k2))
137      DO jk=k1,k2
138         DO jj=j1,j2
139            DO ji=i1,i2
140               interp_scales_u(ji,jj,jk) = ptab(ji,jj,jk)
141            END DO
142         END DO
143      END DO
144      ELSEIF (scaleV) THEN
145      IF (.not.allocated(interp_scales_v)) allocate(interp_scales_v(jpi,jpj,k1:k2))
146      DO jk=k1,k2
147         DO jj=j1,j2
148            DO ji=i1,i2
149               interp_scales_v(ji,jj,jk) = ptab(ji,jj,jk)
150            END DO
151         END DO
152      END DO
153      ENDIF
154      ENDIF
155
156   END SUBROUTINE interpscales
157
158! VERTICAL REFINEMENT END
159
160   SUBROUTINE Agrif_tra
161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE Agrif_tra  ***
163      !!----------------------------------------------------------------------
164      !
165      IF( Agrif_Root() )   RETURN
166      !
167      Agrif_SpecialValue    = 0._wp
168      Agrif_UseSpecialValue = .TRUE.
169      !
170      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
171      !
172      Agrif_UseSpecialValue = .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 ::   j1, j2, i1, i2
185      REAL(wp), POINTER, DIMENSION(:,:) ::   zub, zvb
186      !!---------------------------------------------------------------------- 
187      !
188      IF( Agrif_Root() )   RETURN
189      !
190      CALL wrk_alloc( jpi,jpj,   zub, zvb )
191      !
192      Agrif_SpecialValue    = 0._wp
193      Agrif_UseSpecialValue = ln_spc_dyn
194      !
195      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
196      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
197      !
198      Agrif_UseSpecialValue = .FALSE.
199      !
200      ! prevent smoothing in ghost cells
201      i1 =  1   ;   i2 = jpi
202      j1 =  1   ;   j2 = jpj
203      IF( nbondj == -1 .OR. nbondj == 2 )   j1 = 3
204      IF( nbondj == +1 .OR. nbondj == 2 )   j2 = nlcj-2
205      IF( nbondi == -1 .OR. nbondi == 2 )   i1 = 3
206      IF( nbondi == +1 .OR. nbondi == 2 )   i2 = nlci-2
207
208      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
209         !
210         ! Smoothing
211         ! ---------
212         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
213            ua_b(2,:) = 0._wp
214            DO jk = 1, jpkm1
215               DO jj = 1, jpj
216                  ua_b(2,jj) = ua_b(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
217               END DO
218            END DO
219            DO jj = 1, jpj
220               ua_b(2,jj) = ua_b(2,jj) * r1_hu_a(2,jj)           
221            END DO
222         ENDIF
223         !
224         DO jk=1,jpkm1                 ! Smooth
225            DO jj=j1,j2
226               ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk))
227               ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk)
228            END DO
229         END DO
230         !
231         zub(2,:) = 0._wp              ! Correct transport
232         DO jk = 1, jpkm1
233            DO jj = 1, jpj
234               zub(2,jj) = zub(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
235            END DO
236         END DO
237         DO jj=1,jpj
238            zub(2,jj) = zub(2,jj) * r1_hu_a(2,jj)
239         END DO
240
241         DO jk=1,jpkm1
242            DO jj=1,jpj
243               ua(2,jj,jk) = (ua(2,jj,jk)+ua_b(2,jj)-zub(2,jj))*umask(2,jj,jk)
244            END DO
245         END DO
246
247         ! Set tangential velocities to time splitting estimate
248         !-----------------------------------------------------
249         IF( ln_dynspg_ts ) THEN
250            zvb(2,:) = 0._wp
251            DO jk = 1, jpkm1
252               DO jj = 1, jpj
253                  zvb(2,jj) = zvb(2,jj) + e3v_a(2,jj,jk) * va(2,jj,jk)
254               END DO
255            END DO
256            DO jj = 1, jpj
257               zvb(2,jj) = zvb(2,jj) * r1_hv_a(2,jj)
258            END DO
259            DO jk = 1, jpkm1
260               DO jj = 1, jpj
261                  va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj)) * vmask(2,jj,jk)
262               END DO
263            END DO
264         ENDIF
265         !
266         ! Mask domain edges:
267         !-------------------
268         DO jk = 1, jpkm1
269            DO jj = 1, jpj
270               ua(1,jj,jk) = 0._wp
271               va(1,jj,jk) = 0._wp
272            END DO
273         END DO         
274         !
275      ENDIF
276
277      IF( nbondi == 1 .OR. nbondi == 2 ) THEN
278
279         ! Smoothing
280         ! ---------
281         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
282            ua_b(nlci-2,:) = 0._wp
283            DO jk=1,jpkm1
284               DO jj=1,jpj
285                  ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
286               END DO
287            END DO
288            DO jj=1,jpj
289               ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * r1_hu_a(nlci-2,jj)           
290            END DO
291         ENDIF
292
293         DO jk = 1, jpkm1              ! Smooth
294            DO jj = j1, j2
295               ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk)      &
296                  &             * ( ua(nlci-3,jj,jk) + 2._wp*ua(nlci-2,jj,jk) + ua(nlci-1,jj,jk) )
297            END DO
298         END DO
299
300         zub(nlci-2,:) = 0._wp        ! Correct transport
301         DO jk = 1, jpkm1
302            DO jj = 1, jpj
303               zub(nlci-2,jj) = zub(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
304            END DO
305         END DO
306         DO jj = 1, jpj
307            zub(nlci-2,jj) = zub(nlci-2,jj) * r1_hu_a(nlci-2,jj)
308         END DO
309
310         DO jk = 1, jpkm1
311            DO jj = 1, jpj
312               ua(nlci-2,jj,jk) = ( ua(nlci-2,jj,jk) + ua_b(nlci-2,jj) - zub(nlci-2,jj) ) * umask(nlci-2,jj,jk)
313            END DO
314         END DO
315         !
316         ! Set tangential velocities to time splitting estimate
317         !-----------------------------------------------------
318         IF( ln_dynspg_ts ) THEN
319            zvb(nlci-1,:) = 0._wp
320            DO jk = 1, jpkm1
321               DO jj = 1, jpj
322                  zvb(nlci-1,jj) = zvb(nlci-1,jj) + e3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk)
323               END DO
324            END DO
325            DO jj=1,jpj
326               zvb(nlci-1,jj) = zvb(nlci-1,jj) * r1_hv_a(nlci-1,jj)
327            END DO
328            DO jk = 1, jpkm1
329               DO jj = 1, jpj
330                  va(nlci-1,jj,jk) = ( va(nlci-1,jj,jk) + va_b(nlci-1,jj) - zvb(nlci-1,jj) ) * vmask(nlci-1,jj,jk)
331               END DO
332            END DO
333         ENDIF
334         !
335         ! Mask domain edges:
336         !-------------------
337         DO jk = 1, jpkm1
338            DO jj = 1, jpj
339               ua(nlci-1,jj,jk) = 0._wp
340               va(nlci  ,jj,jk) = 0._wp
341            END DO
342         END DO 
343         !
344      ENDIF
345
346      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
347
348         ! Smoothing
349         ! ---------
350         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
351            va_b(:,2) = 0._wp
352            DO jk = 1, jpkm1
353               DO ji = 1, jpi
354                  va_b(ji,2) = va_b(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk)
355               END DO
356            END DO
357            DO ji=1,jpi
358               va_b(ji,2) = va_b(ji,2) * r1_hv_a(ji,2)           
359            END DO
360         ENDIF
361         !
362         DO jk = 1, jpkm1              ! Smooth
363            DO ji = i1, i2
364               va(ji,2,jk) = 0.25_wp * vmask(ji,2,jk)    &
365                  &        * ( va(ji,1,jk) + 2._wp*va(ji,2,jk) + va(ji,3,jk) )
366            END DO
367         END DO
368         !
369         zvb(:,2) = 0._wp              ! Correct transport
370         DO jk=1,jpkm1
371            DO ji=1,jpi
372               zvb(ji,2) = zvb(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk)
373            END DO
374         END DO
375         DO ji = 1, jpi
376            zvb(ji,2) = zvb(ji,2) * r1_hv_a(ji,2)
377         END DO
378         DO jk = 1, jpkm1
379            DO ji = 1, jpi
380               va(ji,2,jk) = ( va(ji,2,jk) + va_b(ji,2) - zvb(ji,2) ) * vmask(ji,2,jk)
381            END DO
382         END DO
383
384         ! Set tangential velocities to time splitting estimate
385         !-----------------------------------------------------
386         IF( ln_dynspg_ts ) THEN
387            zub(:,2) = 0._wp
388            DO jk = 1, jpkm1
389               DO ji = 1, jpi
390                  zub(ji,2) = zub(ji,2) + e3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk)
391               END DO
392            END DO
393            DO ji = 1, jpi
394               zub(ji,2) = zub(ji,2) * r1_hu_a(ji,2)
395            END DO
396
397            DO jk = 1, jpkm1
398               DO ji = 1, jpi
399                  ua(ji,2,jk) = ( ua(ji,2,jk) + ua_b(ji,2) - zub(ji,2) ) * umask(ji,2,jk)
400               END DO
401            END DO
402         ENDIF
403
404         ! Mask domain edges:
405         !-------------------
406         DO jk = 1, jpkm1
407            DO ji = 1, jpi
408               ua(ji,1,jk) = 0._wp
409               va(ji,1,jk) = 0._wp
410            END DO
411         END DO
412
413      ENDIF
414
415      IF( nbondj == 1 .OR. nbondj == 2 ) THEN
416         !
417         ! Smoothing
418         ! ---------
419         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
420            va_b(:,nlcj-2) = 0._wp
421            DO jk = 1, jpkm1
422               DO ji = 1, jpi
423                  va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk)
424               END DO
425            END DO
426            DO ji = 1, jpi
427               va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)           
428            END DO
429         ENDIF
430         !
431         DO jk = 1, jpkm1              ! Smooth
432            DO ji = i1, i2
433               va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk)   &
434                  &             * ( va(ji,nlcj-3,jk) + 2._wp * va(ji,nlcj-2,jk) + va(ji,nlcj-1,jk) )
435            END DO
436         END DO
437         !
438         zvb(:,nlcj-2) = 0._wp         ! Correct transport
439         DO jk = 1, jpkm1
440            DO ji = 1, jpi
441               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
442            END DO
443         END DO
444         DO ji = 1, jpi
445            zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)
446         END DO
447         DO jk = 1, jpkm1
448            DO ji = 1, jpi
449               va(ji,nlcj-2,jk) = ( va(ji,nlcj-2,jk) + va_b(ji,nlcj-2) - zvb(ji,nlcj-2) ) * vmask(ji,nlcj-2,jk)
450            END DO
451         END DO
452         !
453         ! Set tangential velocities to time splitting estimate
454         !-----------------------------------------------------
455         IF( ln_dynspg_ts ) THEN
456            zub(:,nlcj-1) = 0._wp
457            DO jk = 1, jpkm1
458               DO ji = 1, jpi
459                  zub(ji,nlcj-1) = zub(ji,nlcj-1) + e3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk)
460               END DO
461            END DO
462            DO ji = 1, jpi
463               zub(ji,nlcj-1) = zub(ji,nlcj-1) * r1_hu_a(ji,nlcj-1)
464            END DO
465            !
466            DO jk = 1, jpkm1
467               DO ji = 1, jpi
468                  ua(ji,nlcj-1,jk) = ( ua(ji,nlcj-1,jk) + ua_b(ji,nlcj-1) - zub(ji,nlcj-1) ) * umask(ji,nlcj-1,jk)
469               END DO
470            END DO
471         ENDIF
472         !
473         ! Mask domain edges:
474         !-------------------
475         DO jk = 1, jpkm1
476            DO ji = 1, jpi
477               ua(ji,nlcj  ,jk) = 0._wp
478               va(ji,nlcj-1,jk) = 0._wp
479            END DO
480         END DO 
481         !
482      ENDIF
483      !
484      CALL wrk_dealloc( jpi,jpj,   zub, zvb )
485      !
486   END SUBROUTINE Agrif_dyn
487
488
489   SUBROUTINE Agrif_dyn_ts( jn )
490      !!----------------------------------------------------------------------
491      !!                  ***  ROUTINE Agrif_dyn_ts  ***
492      !!---------------------------------------------------------------------- 
493      !!
494      INTEGER, INTENT(in) ::   jn
495      !!
496      INTEGER :: ji, jj
497      !!---------------------------------------------------------------------- 
498      !
499      IF( Agrif_Root() )   RETURN
500      !
501      IF((nbondi == -1).OR.(nbondi == 2)) THEN
502         DO jj=1,jpj
503            va_e(2,jj) = vbdy_w(jj) * hvr_e(2,jj)
504            ! Specified fluxes:
505            ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj)
506            ! Characteristics method:
507            !alt            ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
508            !alt                       &           - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
509         END DO
510      ENDIF
511      !
512      IF((nbondi == 1).OR.(nbondi == 2)) THEN
513         DO jj=1,jpj
514            va_e(nlci-1,jj) = vbdy_e(jj) * hvr_e(nlci-1,jj)
515            ! Specified fluxes:
516            ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj)
517            ! Characteristics method:
518            !alt            ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
519            !alt                            &           + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
520         END DO
521      ENDIF
522      !
523      IF((nbondj == -1).OR.(nbondj == 2)) THEN
524         DO ji=1,jpi
525            ua_e(ji,2) = ubdy_s(ji) * hur_e(ji,2)
526            ! Specified fluxes:
527            va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2)
528            ! Characteristics method:
529            !alt            va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
530            !alt                       &           - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
531         END DO
532      ENDIF
533      !
534      IF((nbondj == 1).OR.(nbondj == 2)) THEN
535         DO ji=1,jpi
536            ua_e(ji,nlcj-1) = ubdy_n(ji) * hur_e(ji,nlcj-1)
537            ! Specified fluxes:
538            va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2)
539            ! Characteristics method:
540            !alt            va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2)  + va_e(ji,nlcj-3) &
541            !alt                            &           + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
542         END DO
543      ENDIF
544      !
545   END SUBROUTINE Agrif_dyn_ts
546
547
548   SUBROUTINE Agrif_dta_ts( kt )
549      !!----------------------------------------------------------------------
550      !!                  ***  ROUTINE Agrif_dta_ts  ***
551      !!---------------------------------------------------------------------- 
552      !!
553      INTEGER, INTENT(in) ::   kt
554      !!
555      INTEGER :: ji, jj
556      LOGICAL :: ll_int_cons
557      REAL(wp) :: zrhot, zt
558      !!---------------------------------------------------------------------- 
559      !
560      IF( Agrif_Root() )   RETURN
561      !
562      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
563      !
564      zrhot = Agrif_rhot()
565      !
566      ! "Central" time index for interpolation:
567      IF( ln_bt_fw ) THEN
568         zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot
569      ELSE
570         zt = REAL( Agrif_NbStepint()       , wp ) / zrhot
571      ENDIF
572      !
573      ! Linear interpolation of sea level
574      Agrif_SpecialValue    = 0._wp
575      Agrif_UseSpecialValue = .TRUE.
576      CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn )
577      Agrif_UseSpecialValue = .FALSE.
578      !
579      ! Interpolate barotropic fluxes
580      Agrif_SpecialValue=0.
581      Agrif_UseSpecialValue = ln_spc_dyn
582      !
583      IF( ll_int_cons ) THEN  ! Conservative interpolation
584         ! orders matters here !!!!!!
585         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
586         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )
587         bdy_tinterp = 1
588         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
589         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )
590         bdy_tinterp = 2
591         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
592         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )         
593      ELSE ! Linear interpolation
594         bdy_tinterp = 0
595         ubdy_w(:) = 0._wp   ;   vbdy_w(:) = 0._wp 
596         ubdy_e(:) = 0._wp   ;   vbdy_e(:) = 0._wp 
597         ubdy_n(:) = 0._wp   ;   vbdy_n(:) = 0._wp 
598         ubdy_s(:) = 0._wp   ;   vbdy_s(:) = 0._wp
599         CALL Agrif_Bc_variable( unb_id, calledweight=zt, procname=interpunb )
600         CALL Agrif_Bc_variable( vnb_id, calledweight=zt, procname=interpvnb )
601      ENDIF
602      Agrif_UseSpecialValue = .FALSE.
603      !
604   END SUBROUTINE Agrif_dta_ts
605
606
607   SUBROUTINE Agrif_ssh( kt )
608      !!----------------------------------------------------------------------
609      !!                  ***  ROUTINE Agrif_DYN  ***
610      !!---------------------------------------------------------------------- 
611      INTEGER, INTENT(in) ::   kt
612      !!
613      !!---------------------------------------------------------------------- 
614      !
615      IF( Agrif_Root() )   RETURN
616      !
617      IF((nbondi == -1).OR.(nbondi == 2)) THEN
618         ssha(2,:)=ssha(3,:)
619         sshn(2,:)=sshn(3,:)
620      ENDIF
621      !
622      IF((nbondi == 1).OR.(nbondi == 2)) THEN
623         ssha(nlci-1,:)=ssha(nlci-2,:)
624         sshn(nlci-1,:)=sshn(nlci-2,:)
625      ENDIF
626      !
627      IF((nbondj == -1).OR.(nbondj == 2)) THEN
628         ssha(:,2)=ssha(:,3)
629         sshn(:,2)=sshn(:,3)
630      ENDIF
631      !
632      IF((nbondj == 1).OR.(nbondj == 2)) THEN
633         ssha(:,nlcj-1)=ssha(:,nlcj-2)
634         sshn(:,nlcj-1)=sshn(:,nlcj-2)
635      ENDIF
636      !
637   END SUBROUTINE Agrif_ssh
638
639
640   SUBROUTINE Agrif_ssh_ts( jn )
641      !!----------------------------------------------------------------------
642      !!                  ***  ROUTINE Agrif_ssh_ts  ***
643      !!---------------------------------------------------------------------- 
644      INTEGER, INTENT(in) ::   jn
645      !!
646      INTEGER :: ji,jj
647      !!---------------------------------------------------------------------- 
648      !
649      IF((nbondi == -1).OR.(nbondi == 2)) THEN
650         DO jj = 1, jpj
651            ssha_e(2,jj) = hbdy_w(jj)
652         END DO
653      ENDIF
654      !
655      IF((nbondi == 1).OR.(nbondi == 2)) THEN
656         DO jj = 1, jpj
657            ssha_e(nlci-1,jj) = hbdy_e(jj)
658         END DO
659      ENDIF
660      !
661      IF((nbondj == -1).OR.(nbondj == 2)) THEN
662         DO ji = 1, jpi
663            ssha_e(ji,2) = hbdy_s(ji)
664         END DO
665      ENDIF
666      !
667      IF((nbondj == 1).OR.(nbondj == 2)) THEN
668         DO ji = 1, jpi
669            ssha_e(ji,nlcj-1) = hbdy_n(ji)
670         END DO
671      ENDIF
672      !
673   END SUBROUTINE Agrif_ssh_ts
674
675# if defined key_zdftke
676
677   SUBROUTINE Agrif_tke
678      !!----------------------------------------------------------------------
679      !!                  ***  ROUTINE Agrif_tke  ***
680      !!---------------------------------------------------------------------- 
681      REAL(wp) ::   zalpha
682      !!---------------------------------------------------------------------- 
683      !
684      return
685     
686      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
687      IF( zalpha > 1. )   zalpha = 1.
688      !
689      Agrif_SpecialValue    = 0.e0
690      Agrif_UseSpecialValue = .TRUE.
691      !
692      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
693      !
694      Agrif_UseSpecialValue = .FALSE.
695      !
696   END SUBROUTINE Agrif_tke
697   
698# endif
699
700   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
701      !!----------------------------------------------------------------------
702      !!   *** ROUTINE interptsn ***
703      !!----------------------------------------------------------------------
704      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
705      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
706      LOGICAL                                     , INTENT(in   ) ::   before
707      INTEGER                                     , INTENT(in   ) ::   nb , ndir
708      !
709      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
710      INTEGER  ::   imin, imax, jmin, jmax
711      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
712      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
713      LOGICAL :: western_side, eastern_side,northern_side,southern_side
714! VERTICAL REFINEMENT BEGIN
715      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
716      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
717      REAL(wp) :: h_in(k1:k2)
718      REAL(wp) :: h_out(1:jpk)
719      INTEGER :: N_in, N_out
720      REAL(wp) :: h_diff
721      REAL(wp) :: zrhoxy
722! VERTICAL REFINEMENT END
723
724      zrhoxy = Agrif_rhox()*Agrif_rhoy()
725      IF (before) THEN         
726         DO jn = n1,n2-1
727            DO jk=k1,k2
728               DO jj=j1,j2
729                  DO ji=i1,i2
730                     ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
731                  END DO
732               END DO
733            END DO
734         END DO
735         DO jk=k1,k2
736            DO jj=j1,j2
737               DO ji=i1,i2
738                  ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
739               END DO
740            END DO
741         END DO
742
743      ELSE 
744! VERTICAL REFINEMENT BEGIN
745
746         ptab_child(:,:,:,:) = 0.
747         do jj=j1,j2
748         do ji=i1,i2
749           N_in = 0
750           DO jk=k1,k2 !k2 = jpk of parent grid
751             IF (ptab(ji,jj,jk,n2) == 0) EXIT
752             N_in = N_in + 1
753             tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2)
754             h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy)
755           END DO
756           N_out = jpk-1
757           DO jk=1,jpk ! jpk of child grid
758!             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment. Is it just a GYRE issue???????
759!             N_out = N_out + 1
760             h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
761           ENDDO
762           IF (N_in > 0) THEN
763             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
764!             if (h_diff > 0) then
765!               h_in(N_in+1) = h_diff
766!               N_in = N_in + 1
767!             else
768!               h_out(N_out+1) = -h_diff
769!               N_out = N_out + 1
770!             endif
771           ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for?????
772           do jn=1,jpts
773             call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
774           enddo
775!           if (abs(h_diff) > 1000.) then
776!           do jn=1,jpts
777!             do jk=1,N_out
778!             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn)
779!             enddo
780!           enddo
781!         endif
782           ENDIF
783         enddo
784         enddo
785   
786
787! VERTICAL REFINEMENT END
788
789         !
790         western_side  = (nb == 1).AND.(ndir == 1)
791         eastern_side  = (nb == 1).AND.(ndir == 2)
792         southern_side = (nb == 2).AND.(ndir == 1)
793         northern_side = (nb == 2).AND.(ndir == 2)
794         !
795         zrhox = Agrif_Rhox()
796         !
797         zalpha1 = ( zrhox - 1. ) * 0.5
798         zalpha2 = 1. - zalpha1
799         !
800         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
801         zalpha4 = 1. - zalpha3
802         !
803         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
804         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
805         zalpha5 = 1. - zalpha6 - zalpha7
806         !
807         imin = i1
808         imax = i2
809         jmin = j1
810         jmax = j2
811         !
812         ! Remove CORNERS
813         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
814         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
815         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
816         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2       
817         !
818! VERTICAL REFINEMENT BEGIN
819
820! WARNING :
821! ptab replaced by ptab_child in the following
822! k1 k2 replaced by 1 jpk
823! VERTICAL REFINEMENT END
824         IF( eastern_side) THEN
825            DO jn = 1, jpts
826               tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn)
827               DO jk = 1, jpkm1
828                  DO jj = jmin,jmax
829                     IF( umask(nlci-2,jj,jk) == 0._wp ) THEN
830                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
831                     ELSE
832                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
833                        IF( un(nlci-2,jj,jk) > 0._wp ) THEN
834                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
835                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
836                        ENDIF
837                     ENDIF
838                  END DO
839               END DO
840               tsa(nlci,j1:j2,k1:k2,jn) = 0._wp
841            END DO
842         ENDIF
843         !
844         IF( northern_side ) THEN           
845            DO jn = 1, jpts
846               tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn)
847               DO jk = 1, jpkm1
848                  DO ji = imin,imax
849                     IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN
850                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
851                     ELSE
852                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)       
853                        IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN
854                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
855                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
856                        ENDIF
857                     ENDIF
858                  END DO
859               END DO
860               tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp
861            END DO
862         ENDIF
863         !
864         IF( western_side ) THEN           
865            DO jn = 1, jpts
866               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn)
867               DO jk = 1, jpkm1
868                  DO jj = jmin,jmax
869                     IF( umask(2,jj,jk) == 0._wp ) THEN
870                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
871                     ELSE
872                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)       
873                        IF( un(2,jj,jk) < 0._wp ) THEN
874                           tsa(2,jj,jk,jn)=(zalpha6*tsa(3,jj,jk,jn)+zalpha5*tsa(1,jj,jk,jn)+zalpha7*tsa(4,jj,jk,jn))*tmask(2,jj,jk)
875                        ENDIF
876                     ENDIF
877                  END DO
878               END DO
879               tsa(1,j1:j2,k1:k2,jn) = 0._wp
880            END DO
881         ENDIF
882         !
883         IF( southern_side ) THEN           
884            DO jn = 1, jpts
885               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn)
886               DO jk = 1, jpk     
887                  DO ji=imin,imax
888                     IF( vmask(ji,2,jk) == 0._wp ) THEN
889                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
890                     ELSE
891                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
892                        IF( vn(ji,2,jk) < 0._wp ) THEN
893                           tsa(ji,2,jk,jn)=(zalpha6*tsa(ji,3,jk,jn)+zalpha5*tsa(ji,1,jk,jn)+zalpha7*tsa(ji,4,jk,jn))*tmask(ji,2,jk)
894                        ENDIF
895                     ENDIF
896                  END DO
897               END DO
898               tsa(i1:i2,1,k1:k2,jn) = 0._wp
899            END DO
900         ENDIF
901         !
902         ! Treatment of corners
903         !
904         ! East south
905         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
906            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:)
907         ENDIF
908         ! East north
909         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
910            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:)
911         ENDIF
912         ! West south
913         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
914            tsa(2,2,:,:) = ptab_child(2,2,:,:)
915         ENDIF
916         ! West north
917         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
918            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:)
919         ENDIF
920         !
921      ENDIF
922      !
923   END SUBROUTINE interptsn
924
925
926   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
927      !!----------------------------------------------------------------------
928      !!                  ***  ROUTINE interpsshn  ***
929      !!---------------------------------------------------------------------- 
930      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
931      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
932      LOGICAL                         , INTENT(in   ) ::   before
933      INTEGER                         , INTENT(in   ) ::   nb , ndir
934      !
935      LOGICAL :: western_side, eastern_side,northern_side,southern_side
936      !!---------------------------------------------------------------------- 
937      !
938      IF( before) THEN
939         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
940      ELSE
941         western_side  = (nb == 1).AND.(ndir == 1)
942         eastern_side  = (nb == 1).AND.(ndir == 2)
943         southern_side = (nb == 2).AND.(ndir == 1)
944         northern_side = (nb == 2).AND.(ndir == 2)
945         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
946         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
947         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
948         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
949      ENDIF
950      !
951   END SUBROUTINE interpsshn
952
953   SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir)
954      !!---------------------------------------------
955      !!   *** ROUTINE interpun ***
956      !!---------------------------------------------   
957      !!
958      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
959      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
960      LOGICAL, INTENT(in) :: before
961      INTEGER, INTENT(in) :: nb , ndir
962      !!
963      INTEGER :: ji,jj,jk
964      REAL(wp) :: zrhoy
965! VERTICAL REFINEMENT BEGIN
966      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
967      REAL(wp), DIMENSION(k1:k2) :: tabin
968      REAL(wp) :: h_in(k1:k2)
969      REAL(wp) :: h_out(1:jpk)
970      INTEGER :: N_in, N_out
971      REAL(wp) :: h_diff
972      LOGICAL :: western_side, eastern_side
973      INTEGER :: iref
974
975! VERTICAL REFINEMENT END
976      !!---------------------------------------------   
977      !
978      zrhoy = Agrif_rhoy()
979      IF (before) THEN
980         DO jk=1,jpk
981            DO jj=j1,j2
982               DO ji=i1,i2
983                  ptab(ji,jj,jk,1) = e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)
984                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
985               END DO
986            END DO
987         END DO
988      ELSE
989! VERTICAL REFINEMENT BEGIN
990         western_side  = (nb == 1).AND.(ndir == 1)
991         eastern_side  = (nb == 1).AND.(ndir == 2)
992         
993         ptab_child(:,:,:) = 0.
994         DO jj=j1,j2
995            DO ji=i1,i2
996               iref = ji
997               IF (western_side) iref = 2
998               IF (eastern_side) iref = nlci-2
999
1000               N_in = 0
1001               DO jk=k1,k2
1002                  IF (ptab(ji,jj,jk,2) == 0) EXIT
1003                  N_in = N_in + 1
1004                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1005                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj) * zrhoy) 
1006              ENDDO
1007         
1008              IF (N_in == 0) THEN
1009                 ptab_child(ji,jj,:) = 0.
1010                 CYCLE
1011              ENDIF
1012         
1013              N_out = 0
1014              DO jk=1,jpk
1015                 if (umask(iref,jj,jk) == 0) EXIT
1016                 N_out = N_out + 1
1017                 h_out(N_out) = e3u_n(ji,jj,jk)
1018              ENDDO
1019         
1020              IF (N_out == 0) THEN
1021                 ptab_child(ji,jj,:) = 0.
1022                 CYCLE
1023              ENDIF
1024         
1025!              IF (N_in * N_out > 0) THEN
1026!                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
1027! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
1028!                 if (h_diff < 0.) then
1029!                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
1030!                    stop
1031!                 endif
1032!              ENDIF
1033              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
1034         
1035            ENDDO
1036         ENDDO
1037
1038! in the following
1039! remove division of ua by fs e3u (already done) and also zrhoy and e2u
1040! VERTICAL REFINEMENT END
1041
1042         DO jk = 1, jpkm1
1043            DO jj=j1,j2
1044               ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
1045!/(zrhoy*e2u(i1:i2,jj)))
1046            END DO
1047         END DO
1048      ENDIF
1049      !
1050   END SUBROUTINE interpun
1051
1052
1053
1054   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir)
1055      !!---------------------------------------------
1056      !!   *** ROUTINE interpvn ***
1057      !!----------------------------------------------------------------------
1058      !
1059      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
1060      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1061      LOGICAL, INTENT(in) :: before
1062      INTEGER, INTENT(in) :: nb , ndir
1063      !
1064      INTEGER :: ji,jj,jk
1065      REAL(wp) :: zrhox
1066! VERTICAL REFINEMENT BEGIN
1067      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
1068      REAL(wp), DIMENSION(k1:k2) :: tabin
1069      REAL(wp) :: h_in(k1:k2)
1070      REAL(wp) :: h_out(1:jpk)
1071      INTEGER :: N_in, N_out
1072      REAL(wp) :: h_diff
1073      LOGICAL :: northern_side,southern_side
1074      INTEGER :: jref
1075
1076! VERTICAL REFINEMENT END
1077      !!---------------------------------------------   
1078      !     
1079      zrhox = Agrif_rhox()
1080      IF (before) THEN         
1081         DO jk=k1,k2
1082            DO jj=j1,j2
1083               DO ji=i1,i2
1084                  ptab(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
1085                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
1086               END DO
1087            END DO
1088         END DO
1089      ELSE       
1090! VERTICAL REFINEMENT BEGIN
1091         ptab_child(:,:,:) = 0.
1092         southern_side = (nb == 2).AND.(ndir == 1)
1093         northern_side = (nb == 2).AND.(ndir == 2)
1094         do jj=j1,j2
1095            jref = jj
1096            IF (southern_side) jref = 2
1097            IF (northern_side) jref = nlcj-2
1098            do ji=i1,i2
1099
1100               N_in = 0
1101               do jk=k1,k2
1102                  if (ptab(ji,jj,jk,2) == 0) EXIT
1103                  N_in = N_in + 1
1104                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1105                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1106               enddo
1107               IF (N_in == 0) THEN
1108                  ptab_child(ji,jj,:) = 0.
1109                  CYCLE
1110               ENDIF
1111         
1112               N_out = 0
1113               do jk=1,jpk
1114                  if (vmask(ji,jref,jk) == 0) EXIT
1115                  N_out = N_out + 1
1116                  h_out(N_out) = e3v_n(ji,jj,jk)
1117               enddo
1118               IF (N_out == 0) THEN
1119                 ptab_child(ji,jj,:) = 0.
1120                 CYCLE
1121               ENDIF
1122
1123!               IF (N_in * N_out > 0) THEN
1124!                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
1125! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
1126!                  if (h_diff < 0.) then
1127!                     print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
1128!                     stop
1129!                  endif
1130!               ENDIF
1131               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
1132         
1133            enddo
1134         enddo
1135! in the following
1136! remove division of va by fs e3v, zrhox and e1v (already done)
1137! VERTICAL REFINEMENT END
1138         DO jk=1,jpkm1
1139            DO jj=j1,j2
1140               va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
1141            END DO
1142         END DO
1143      ENDIF
1144      !       
1145   END SUBROUTINE interpvn
1146   
1147
1148   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
1149      !!----------------------------------------------------------------------
1150      !!                  ***  ROUTINE interpunb  ***
1151      !!---------------------------------------------------------------------- 
1152      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1153      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1154      LOGICAL                         , INTENT(in   ) ::   before
1155      INTEGER                         , INTENT(in   ) ::   nb , ndir
1156      !
1157      INTEGER  ::   ji, jj
1158      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1159      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1160      !!---------------------------------------------------------------------- 
1161      !
1162      IF( before ) THEN
1163         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
1164      ELSE
1165         western_side  = (nb == 1).AND.(ndir == 1)
1166         eastern_side  = (nb == 1).AND.(ndir == 2)
1167         southern_side = (nb == 2).AND.(ndir == 1)
1168         northern_side = (nb == 2).AND.(ndir == 2)
1169         zrhoy = Agrif_Rhoy()
1170         zrhot = Agrif_rhot()
1171         ! Time indexes bounds for integration
1172         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1173         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1174         ! Polynomial interpolation coefficients:
1175         IF( bdy_tinterp == 1 ) THEN
1176            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1177               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1178         ELSEIF( bdy_tinterp == 2 ) THEN
1179            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1180               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1181
1182         ELSE
1183            ztcoeff = 1
1184         ENDIF
1185         !   
1186         IF(western_side) THEN
1187            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1188         ENDIF
1189         IF(eastern_side) THEN
1190            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1191         ENDIF
1192         IF(southern_side) THEN
1193            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1194         ENDIF
1195         IF(northern_side) THEN
1196            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1197         ENDIF
1198         !           
1199         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1200            IF(western_side) THEN
1201               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
1202            ENDIF
1203            IF(eastern_side) THEN
1204               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
1205            ENDIF
1206            IF(southern_side) THEN
1207               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
1208            ENDIF
1209            IF(northern_side) THEN
1210               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
1211            ENDIF
1212         ENDIF
1213      ENDIF
1214      !
1215   END SUBROUTINE interpunb
1216
1217
1218   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1219      !!----------------------------------------------------------------------
1220      !!                  ***  ROUTINE interpvnb  ***
1221      !!---------------------------------------------------------------------- 
1222      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1223      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1224      LOGICAL                         , INTENT(in   ) ::   before
1225      INTEGER                         , INTENT(in   ) ::   nb , ndir
1226      !
1227      INTEGER  ::   ji,jj
1228      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1229      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1230      !!---------------------------------------------------------------------- 
1231      !
1232      IF( before ) THEN
1233         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1234      ELSE
1235         western_side  = (nb == 1).AND.(ndir == 1)
1236         eastern_side  = (nb == 1).AND.(ndir == 2)
1237         southern_side = (nb == 2).AND.(ndir == 1)
1238         northern_side = (nb == 2).AND.(ndir == 2)
1239         zrhox = Agrif_Rhox()
1240         zrhot = Agrif_rhot()
1241         ! Time indexes bounds for integration
1242         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1243         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1244         IF( bdy_tinterp == 1 ) THEN
1245            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1246               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1247         ELSEIF( bdy_tinterp == 2 ) THEN
1248            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1249               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1250         ELSE
1251            ztcoeff = 1
1252         ENDIF
1253         !
1254         IF(western_side) THEN
1255            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1256         ENDIF
1257         IF(eastern_side) THEN
1258            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1259         ENDIF
1260         IF(southern_side) THEN
1261            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1262         ENDIF
1263         IF(northern_side) THEN
1264            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1265         ENDIF
1266         !           
1267         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1268            IF(western_side) THEN
1269               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1270                     &                                  * vmask(i1,j1:j2,1)
1271            ENDIF
1272            IF(eastern_side) THEN
1273               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1274                     &                                  * vmask(i1,j1:j2,1)
1275            ENDIF
1276            IF(southern_side) THEN
1277               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1278                     &                                  * vmask(i1:i2,j1,1)
1279            ENDIF
1280            IF(northern_side) THEN
1281               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1282                     &                                  * vmask(i1:i2,j1,1)
1283            ENDIF
1284         ENDIF
1285      ENDIF
1286      !
1287   END SUBROUTINE interpvnb
1288
1289
1290   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1291      !!----------------------------------------------------------------------
1292      !!                  ***  ROUTINE interpub2b  ***
1293      !!---------------------------------------------------------------------- 
1294      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1295      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1296      LOGICAL                         , INTENT(in   ) ::   before
1297      INTEGER                         , INTENT(in   ) ::   nb , ndir
1298      !
1299      INTEGER  ::   ji,jj
1300      REAL(wp) ::   zrhot, zt0, zt1,zat
1301      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1302      !!---------------------------------------------------------------------- 
1303      IF( before ) THEN
1304         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1305      ELSE
1306         western_side  = (nb == 1).AND.(ndir == 1)
1307         eastern_side  = (nb == 1).AND.(ndir == 2)
1308         southern_side = (nb == 2).AND.(ndir == 1)
1309         northern_side = (nb == 2).AND.(ndir == 2)
1310         zrhot = Agrif_rhot()
1311         ! Time indexes bounds for integration
1312         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1313         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1314         ! Polynomial interpolation coefficients:
1315         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1316            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1317         !
1318         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1319         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1320         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1321         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1322      ENDIF
1323      !
1324   END SUBROUTINE interpub2b
1325   
1326
1327   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1328      !!----------------------------------------------------------------------
1329      !!                  ***  ROUTINE interpvb2b  ***
1330      !!---------------------------------------------------------------------- 
1331      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1332      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1333      LOGICAL                         , INTENT(in   ) ::   before
1334      INTEGER                         , INTENT(in   ) ::   nb , ndir
1335      !
1336      INTEGER ::   ji,jj
1337      REAL(wp) ::   zrhot, zt0, zt1,zat
1338      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
1339      !!---------------------------------------------------------------------- 
1340      !
1341      IF( before ) THEN
1342         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1343      ELSE     
1344         western_side  = (nb == 1).AND.(ndir == 1)
1345         eastern_side  = (nb == 1).AND.(ndir == 2)
1346         southern_side = (nb == 2).AND.(ndir == 1)
1347         northern_side = (nb == 2).AND.(ndir == 2)
1348         zrhot = Agrif_rhot()
1349         ! Time indexes bounds for integration
1350         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1351         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1352         ! Polynomial interpolation coefficients:
1353         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1354            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1355         !
1356         IF(western_side )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1357         IF(eastern_side )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1358         IF(southern_side)   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1359         IF(northern_side)   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1360      ENDIF
1361      !     
1362   END SUBROUTINE interpvb2b
1363
1364
1365   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1366      !!----------------------------------------------------------------------
1367      !!                  ***  ROUTINE interpe3t  ***
1368      !!---------------------------------------------------------------------- 
1369      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1370      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1371      LOGICAL                              , INTENT(in   ) :: before
1372      INTEGER                              , INTENT(in   ) :: nb , ndir
1373      !
1374      INTEGER :: ji, jj, jk
1375      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1376      REAL(wp) :: ztmpmsk     
1377      !!---------------------------------------------------------------------- 
1378      !   
1379      IF( before ) THEN
1380         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1381      ELSE
1382         western_side  = (nb == 1).AND.(ndir == 1)
1383         eastern_side  = (nb == 1).AND.(ndir == 2)
1384         southern_side = (nb == 2).AND.(ndir == 1)
1385         northern_side = (nb == 2).AND.(ndir == 2)
1386
1387         DO jk = k1, k2
1388            DO jj = j1, j2
1389               DO ji = i1, i2
1390                  ! Get velocity mask at boundary edge points:
1391                  IF( western_side )   ztmpmsk = umask(ji    ,jj    ,1)
1392                  IF( eastern_side )   ztmpmsk = umask(nlci-2,jj    ,1)
1393                  IF( northern_side)   ztmpmsk = vmask(ji    ,nlcj-2,1)
1394                  IF( southern_side)   ztmpmsk = vmask(ji    ,2     ,1)
1395                  !
1396                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
1397                     IF (western_side) THEN
1398                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1399                     ELSEIF (eastern_side) THEN
1400                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1401                     ELSEIF (southern_side) THEN
1402                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1403                     ELSEIF (northern_side) THEN
1404                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1405                     ENDIF
1406                     WRITE(numout,*) '      ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1407                     kindic_agr = kindic_agr + 1
1408                  ENDIF
1409               END DO
1410            END DO
1411         END DO
1412         !
1413      ENDIF
1414      !
1415   END SUBROUTINE interpe3t
1416
1417
1418   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1419      !!----------------------------------------------------------------------
1420      !!                  ***  ROUTINE interpumsk  ***
1421      !!---------------------------------------------------------------------- 
1422      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1423      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1424      LOGICAL                              , INTENT(in   ) ::   before
1425      INTEGER                              , INTENT(in   ) ::   nb , ndir
1426      !
1427      INTEGER ::   ji, jj, jk
1428      LOGICAL ::   western_side, eastern_side   
1429      !!---------------------------------------------------------------------- 
1430      !   
1431      IF( before ) THEN
1432         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1433      ELSE
1434         western_side = (nb == 1).AND.(ndir == 1)
1435         eastern_side = (nb == 1).AND.(ndir == 2)
1436         DO jk = k1, k2
1437            DO jj = j1, j2
1438               DO ji = i1, i2
1439                   ! Velocity mask at boundary edge points:
1440                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1441                     IF (western_side) THEN
1442                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1443                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1444                        kindic_agr = kindic_agr + 1
1445                     ELSEIF (eastern_side) THEN
1446                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1447                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1448                        kindic_agr = kindic_agr + 1
1449                     ENDIF
1450                  ENDIF
1451               END DO
1452            END DO
1453         END DO
1454         !
1455      ENDIF
1456      !
1457   END SUBROUTINE interpumsk
1458
1459
1460   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1461      !!----------------------------------------------------------------------
1462      !!                  ***  ROUTINE interpvmsk  ***
1463      !!---------------------------------------------------------------------- 
1464      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1465      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1466      LOGICAL                              , INTENT(in   ) ::   before
1467      INTEGER                              , INTENT(in   ) :: nb , ndir
1468      !
1469      INTEGER ::   ji, jj, jk
1470      LOGICAL ::   northern_side, southern_side     
1471      !!---------------------------------------------------------------------- 
1472      !   
1473      IF( before ) THEN
1474         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1475      ELSE
1476         southern_side = (nb == 2).AND.(ndir == 1)
1477         northern_side = (nb == 2).AND.(ndir == 2)
1478         DO jk = k1, k2
1479            DO jj = j1, j2
1480               DO ji = i1, i2
1481                   ! Velocity mask at boundary edge points:
1482                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1483                     IF (southern_side) THEN
1484                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1485                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1486                        kindic_agr = kindic_agr + 1
1487                     ELSEIF (northern_side) THEN
1488                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1489                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1490                        kindic_agr = kindic_agr + 1
1491                     ENDIF
1492                  ENDIF
1493               END DO
1494            END DO
1495         END DO
1496         !
1497      ENDIF
1498      !
1499   END SUBROUTINE interpvmsk
1500
1501# if defined key_zdftke
1502
1503   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
1504      !!----------------------------------------------------------------------
1505      !!                  ***  ROUTINE interavm  ***
1506      !!---------------------------------------------------------------------- 
1507      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1508      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1509      LOGICAL                              , INTENT(in   ) ::   before
1510      !!---------------------------------------------------------------------- 
1511      !     
1512      IF( before ) THEN
1513         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1514      ELSE
1515         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1516      ENDIF
1517      !
1518   END SUBROUTINE interpavm
1519
1520# endif /* key_zdftke */
1521
1522#else
1523   !!----------------------------------------------------------------------
1524   !!   Empty module                                          no AGRIF zoom
1525   !!----------------------------------------------------------------------
1526CONTAINS
1527   SUBROUTINE Agrif_OPA_Interp_empty
1528      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1529   END SUBROUTINE Agrif_OPA_Interp_empty
1530#endif
1531
1532   !!======================================================================
1533END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.