New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_oce_interp.F90 in NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90 @ 11590

Last change on this file since 11590 was 11590, checked in by jchanut, 5 years ago

#2222: 1) create remapping module (vremap) and integration of D. Engwirda piecewise polynomial recontruction package (PPR_LIB cpp key). 2) Various bug corrections with key_vertical activated.

  • Property svn:keywords set to Id
File size: 45.0 KB
Line 
1MODULE agrif_oce_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_oce_interp  ***
4   !! AGRIF: interpolation package for the ocean dynamics (OPA)
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (L. Debreu)  Original cade
7   !!            3.2  !  2009-04  (R. Benshila)
8   !!            3.6  !  2014-09  (R. Benshila)
9   !!----------------------------------------------------------------------
10#if defined key_agrif
11   !!----------------------------------------------------------------------
12   !!   'key_agrif'                                              AGRIF zoom
13   !!----------------------------------------------------------------------
14   !!   Agrif_tra     :
15   !!   Agrif_dyn     :
16   !!   Agrif_ssh     :
17   !!   Agrif_dyn_ts  :
18   !!   Agrif_dta_ts  :
19   !!   Agrif_ssh_ts  :
20   !!   Agrif_avm     :
21   !!   interpu       :
22   !!   interpv       :
23   !!----------------------------------------------------------------------
24   USE par_oce
25   USE oce
26   USE dom_oce     
27   USE zdf_oce
28   USE agrif_oce
29   USE phycst
30   USE dynspg_ts, ONLY: un_adv, vn_adv
31   !
32   USE in_out_manager
33   USE agrif_oce_sponge
34   USE lib_mpp
35   USE vremap
36 
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts
41   PUBLIC   Agrif_tra, Agrif_avm
42   PUBLIC   interpun , interpvn
43   PUBLIC   interptsn, interpsshn, interpavm
44   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
45   PUBLIC   interpe3t, interpumsk, interpvmsk
46
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/NST 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE Agrif_tra
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE Agrif_tra  ***
58      !!----------------------------------------------------------------------
59      !
60      IF( Agrif_Root() )   RETURN
61      !
62      Agrif_SpecialValue    = 0._wp
63      Agrif_UseSpecialValue = .TRUE.
64      !
65      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
66      !
67      Agrif_UseSpecialValue = .FALSE.
68      !
69   END SUBROUTINE Agrif_tra
70
71
72   SUBROUTINE Agrif_dyn( kt )
73      !!----------------------------------------------------------------------
74      !!                  ***  ROUTINE Agrif_DYN  ***
75      !!---------------------------------------------------------------------- 
76      INTEGER, INTENT(in) ::   kt
77      !
78      INTEGER ::   ji, jj, jk       ! dummy loop indices
79      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
80      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
81      !!---------------------------------------------------------------------- 
82      !
83      IF( Agrif_Root() )   RETURN
84      !
85      Agrif_SpecialValue    = 0._wp
86      Agrif_UseSpecialValue = ln_spc_dyn
87      !
88      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
89      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
90      !
91      Agrif_UseSpecialValue = .FALSE.
92      !
93      ! --- West --- !
94      ibdy1 = 2
95      ibdy2 = 1+nbghostcells 
96      !
97      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
98         DO ji = mi0(ibdy1), mi1(ibdy2)
99            ua_b(ji,:) = 0._wp
100
101            DO jk = 1, jpkm1
102               DO jj = 1, jpj
103                  ua_b(ji,jj) = ua_b(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
104               END DO
105            END DO
106
107            DO jj = 1, jpj
108               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
109            END DO
110         END DO
111      ENDIF
112      !
113      DO ji = mi0(ibdy1), mi1(ibdy2)
114         zub(ji,:) = 0._wp    ! Correct transport
115         DO jk = 1, jpkm1
116            DO jj = 1, jpj
117               zub(ji,jj) = zub(ji,jj) & 
118                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk)*umask(ji,jj,jk)
119            END DO
120         END DO
121         DO jj=1,jpj
122            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
123         END DO
124           
125         DO jk = 1, jpkm1
126            DO jj = 1, jpj
127               ua(ji,jj,jk) = ( ua(ji,jj,jk) + ua_b(ji,jj)-zub(ji,jj)) * umask(ji,jj,jk)
128            END DO
129         END DO
130      END DO
131           
132      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
133         DO ji = mi0(ibdy1), mi1(ibdy2)
134            zvb(ji,:) = 0._wp
135            DO jk = 1, jpkm1
136               DO jj = 1, jpj
137                  zvb(ji,jj) = zvb(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
138               END DO
139            END DO
140            DO jj = 1, jpj
141               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
142            END DO
143            DO jk = 1, jpkm1
144               DO jj = 1, jpj
145                  va(ji,jj,jk) = ( va(ji,jj,jk) + va_b(ji,jj)-zvb(ji,jj))*vmask(ji,jj,jk)
146               END DO
147            END DO
148         END DO
149      ENDIF
150
151      ! --- East --- !
152      ibdy1 = jpiglo-1-nbghostcells
153      ibdy2 = jpiglo-2 
154      !
155      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
156         DO ji = mi0(ibdy1), mi1(ibdy2)
157            ua_b(ji,:) = 0._wp
158            DO jk = 1, jpkm1
159               DO jj = 1, jpj
160                  ua_b(ji,jj) = ua_b(ji,jj) & 
161                      & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
162               END DO
163            END DO
164            DO jj = 1, jpj
165               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
166            END DO
167         END DO
168      ENDIF
169      !
170      DO ji = mi0(ibdy1), mi1(ibdy2)
171         zub(ji,:) = 0._wp    ! Correct transport
172         DO jk = 1, jpkm1
173            DO jj = 1, jpj
174               zub(ji,jj) = zub(ji,jj) & 
175                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk) * umask(ji,jj,jk)
176            END DO
177         END DO
178         DO jj=1,jpj
179            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
180         END DO
181           
182         DO jk = 1, jpkm1
183            DO jj = 1, jpj
184               ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
185                 & + ua_b(ji,jj)-zub(ji,jj))*umask(ji,jj,jk)
186            END DO
187         END DO
188      END DO
189           
190      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
191         ibdy1 = jpiglo-nbghostcells
192         ibdy2 = jpiglo-1 
193         DO ji = mi0(ibdy1), mi1(ibdy2)
194            zvb(ji,:) = 0._wp
195            DO jk = 1, jpkm1
196               DO jj = 1, jpj
197                  zvb(ji,jj) = zvb(ji,jj) &
198                     & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
199               END DO
200            END DO
201            DO jj = 1, jpj
202               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
203            END DO
204            DO jk = 1, jpkm1
205               DO jj = 1, jpj
206                  va(ji,jj,jk) = ( va(ji,jj,jk) & 
207                      & + va_b(ji,jj)-zvb(ji,jj)) * vmask(ji,jj,jk)
208               END DO
209            END DO
210         END DO
211      ENDIF
212
213      ! --- South --- !
214      jbdy1 = 2
215      jbdy2 = 1+nbghostcells 
216      !
217      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
218         DO jj = mj0(jbdy1), mj1(jbdy2)
219            va_b(:,jj) = 0._wp
220            DO jk = 1, jpkm1
221               DO ji = 1, jpi
222                  va_b(ji,jj) = va_b(ji,jj) & 
223                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
224               END DO
225            END DO
226            DO ji=1,jpi
227               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)     
228            END DO
229         END DO
230      ENDIF
231      !
232      DO jj = mj0(jbdy1), mj1(jbdy2)
233         zvb(:,jj) = 0._wp    ! Correct transport
234         DO jk=1,jpkm1
235            DO ji=1,jpi
236               zvb(ji,jj) = zvb(ji,jj) & 
237                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
238            END DO
239         END DO
240         DO ji = 1, jpi
241            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
242         END DO
243
244         DO jk = 1, jpkm1
245            DO ji = 1, jpi
246               va(ji,jj,jk) = ( va(ji,jj,jk) & 
247                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
248            END DO
249         END DO
250      END DO
251           
252      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
253         DO jj = mj0(jbdy1), mj1(jbdy2)
254            zub(:,jj) = 0._wp
255            DO jk = 1, jpkm1
256               DO ji = 1, jpi
257                  zub(ji,jj) = zub(ji,jj) & 
258                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
259               END DO
260            END DO
261            DO ji = 1, jpi
262               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
263            END DO
264               
265            DO jk = 1, jpkm1
266               DO ji = 1, jpi
267                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
268                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
269               END DO
270            END DO
271         END DO
272      ENDIF
273
274      ! --- North --- !
275      jbdy1 = jpjglo-1-nbghostcells
276      jbdy2 = jpjglo-2 
277      !
278      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
279         DO jj = mj0(jbdy1), mj1(jbdy2)
280            va_b(:,jj) = 0._wp
281            DO jk = 1, jpkm1
282               DO ji = 1, jpi
283                  va_b(ji,jj) = va_b(ji,jj) & 
284                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
285               END DO
286            END DO
287            DO ji=1,jpi
288               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)
289            END DO
290         END DO
291      ENDIF
292      !
293      DO jj = mj0(jbdy1), mj1(jbdy2)
294         zvb(:,jj) = 0._wp    ! Correct transport
295         DO jk=1,jpkm1
296            DO ji=1,jpi
297               zvb(ji,jj) = zvb(ji,jj) & 
298                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
299            END DO
300         END DO
301         DO ji = 1, jpi
302            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
303         END DO
304
305         DO jk = 1, jpkm1
306            DO ji = 1, jpi
307               va(ji,jj,jk) = ( va(ji,jj,jk) & 
308                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
309            END DO
310         END DO
311      END DO
312           
313      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
314         jbdy1 = jpjglo-nbghostcells
315         jbdy2 = jpjglo-1
316         DO jj = mj0(jbdy1), mj1(jbdy2)
317            zub(:,jj) = 0._wp
318            DO jk = 1, jpkm1
319               DO ji = 1, jpi
320                  zub(ji,jj) = zub(ji,jj) & 
321                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
322               END DO
323            END DO
324            DO ji = 1, jpi
325               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
326            END DO
327               
328            DO jk = 1, jpkm1
329               DO ji = 1, jpi
330                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
331                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
332               END DO
333            END DO
334         END DO
335      ENDIF
336      !
337   END SUBROUTINE Agrif_dyn
338
339
340   SUBROUTINE Agrif_dyn_ts( jn )
341      !!----------------------------------------------------------------------
342      !!                  ***  ROUTINE Agrif_dyn_ts  ***
343      !!---------------------------------------------------------------------- 
344      INTEGER, INTENT(in) ::   jn
345      !!
346      INTEGER :: ji, jj
347      INTEGER :: istart, iend, jstart, jend
348      !!---------------------------------------------------------------------- 
349      !
350      IF( Agrif_Root() )   RETURN
351      !
352      !--- West ---!
353      istart = 2
354      iend   = nbghostcells+1
355      DO ji = mi0(istart), mi1(iend)
356         DO jj=1,jpj
357            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
358            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
359         END DO
360      END DO
361      !
362      !--- East ---!
363      istart = jpiglo-nbghostcells
364      iend   = jpiglo-1
365      DO ji = mi0(istart), mi1(iend)
366         DO jj=1,jpj
367            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
368         END DO
369      END DO
370      istart = jpiglo-nbghostcells-1
371      iend   = jpiglo-2
372      DO ji = mi0(istart), mi1(iend)
373         DO jj=1,jpj
374            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
375         END DO
376      END DO
377      !
378      !--- South ---!
379      jstart = 2
380      jend   = nbghostcells+1
381      DO jj = mj0(jstart), mj1(jend)
382         DO ji=1,jpi
383            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
384            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
385         END DO
386      END DO
387      !
388      !--- North ---!
389      jstart = jpjglo-nbghostcells
390      jend   = jpjglo-1
391      DO jj = mj0(jstart), mj1(jend)
392         DO ji=1,jpi
393            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
394         END DO
395      END DO
396      jstart = jpjglo-nbghostcells-1
397      jend   = jpjglo-2
398      DO jj = mj0(jstart), mj1(jend)
399         DO ji=1,jpi
400            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
401         END DO
402      END DO
403      !
404   END SUBROUTINE Agrif_dyn_ts
405
406   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
407      !!----------------------------------------------------------------------
408      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
409      !!---------------------------------------------------------------------- 
410      INTEGER, INTENT(in) ::   jn
411      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
412      !!
413      INTEGER :: ji, jj
414      INTEGER :: istart, iend, jstart, jend
415      !!---------------------------------------------------------------------- 
416      !
417      IF( Agrif_Root() )   RETURN
418      !
419      !--- West ---!
420      istart = 2
421      iend   = nbghostcells+1
422      DO ji = mi0(istart), mi1(iend)
423         DO jj=1,jpj
424            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
425            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
426         END DO
427      END DO
428      !
429      !--- East ---!
430      istart = jpiglo-nbghostcells
431      iend   = jpiglo-1
432      DO ji = mi0(istart), mi1(iend)
433         DO jj=1,jpj
434            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
435         END DO
436      END DO
437      istart = jpiglo-nbghostcells-1
438      iend   = jpiglo-2
439      DO ji = mi0(istart), mi1(iend)
440         DO jj=1,jpj
441            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
442         END DO
443      END DO
444      !
445      !--- South ---!
446      jstart = 2
447      jend   = nbghostcells+1
448      DO jj = mj0(jstart), mj1(jend)
449         DO ji=1,jpi
450            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
451            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
452         END DO
453      END DO
454      !
455      !--- North ---!
456      jstart = jpjglo-nbghostcells
457      jend   = jpjglo-1
458      DO jj = mj0(jstart), mj1(jend)
459         DO ji=1,jpi
460            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
461         END DO
462      END DO
463      jstart = jpjglo-nbghostcells-1
464      jend   = jpjglo-2
465      DO jj = mj0(jstart), mj1(jend)
466         DO ji=1,jpi
467            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
468         END DO
469      END DO
470      !
471   END SUBROUTINE Agrif_dyn_ts_flux
472
473   SUBROUTINE Agrif_dta_ts( kt )
474      !!----------------------------------------------------------------------
475      !!                  ***  ROUTINE Agrif_dta_ts  ***
476      !!---------------------------------------------------------------------- 
477      INTEGER, INTENT(in) ::   kt
478      !!
479      INTEGER :: ji, jj
480      LOGICAL :: ll_int_cons
481      !!---------------------------------------------------------------------- 
482      !
483      IF( Agrif_Root() )   RETURN
484      !
485      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
486      !
487      ! Enforce volume conservation if no time refinement: 
488      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
489      !
490      ! Interpolate barotropic fluxes
491      Agrif_SpecialValue = 0._wp
492      Agrif_UseSpecialValue = ln_spc_dyn
493      !
494      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
495      utint_stage(:,:) = 0
496      vtint_stage(:,:) = 0
497      !
498      IF( ll_int_cons ) THEN  ! Conservative interpolation
499         ! order matters here !!!!!!
500         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
501         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
502         !
503         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
504         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )
505         !
506         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
507         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )         
508      ELSE ! Linear interpolation
509         !
510         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
511         CALL Agrif_Bc_variable( unb_id, procname=interpunb )
512         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
513      ENDIF
514      Agrif_UseSpecialValue = .FALSE.
515      !
516   END SUBROUTINE Agrif_dta_ts
517
518
519   SUBROUTINE Agrif_ssh( kt )
520      !!----------------------------------------------------------------------
521      !!                  ***  ROUTINE Agrif_ssh  ***
522      !!---------------------------------------------------------------------- 
523      INTEGER, INTENT(in) ::   kt
524      !
525      INTEGER  :: ji, jj
526      INTEGER  :: istart, iend, jstart, jend
527      !!---------------------------------------------------------------------- 
528      !
529      IF( Agrif_Root() )   RETURN
530      !     
531      ! Linear time interpolation of sea level
532      !
533      Agrif_SpecialValue    = 0._wp
534      Agrif_UseSpecialValue = .TRUE.
535      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
536      Agrif_UseSpecialValue = .FALSE.
537      !
538      ! --- West --- !
539      istart = 2
540      iend   = 1 + nbghostcells
541      DO ji = mi0(istart), mi1(iend)
542         DO jj = 1, jpj
543            ssha(ji,jj) = hbdy(ji,jj)
544         ENDDO
545      ENDDO
546      !
547      ! --- East --- !
548      istart = jpiglo - nbghostcells
549      iend   = jpiglo - 1
550      DO ji = mi0(istart), mi1(iend)
551         DO jj = 1, jpj
552            ssha(ji,jj) = hbdy(ji,jj)
553         ENDDO
554      ENDDO
555      !
556      ! --- South --- !
557      jstart = 2
558      jend   = 1 + nbghostcells
559      DO jj = mj0(jstart), mj1(jend)
560         DO ji = 1, jpi
561            ssha(ji,jj) = hbdy(ji,jj)
562         ENDDO
563      ENDDO
564      !
565      ! --- North --- !
566      jstart = jpjglo - nbghostcells
567      jend   = jpjglo - 1
568      DO jj = mj0(jstart), mj1(jend)
569         DO ji = 1, jpi
570            ssha(ji,jj) = hbdy(ji,jj)
571         ENDDO
572      ENDDO
573      !
574   END SUBROUTINE Agrif_ssh
575
576
577   SUBROUTINE Agrif_ssh_ts( jn )
578      !!----------------------------------------------------------------------
579      !!                  ***  ROUTINE Agrif_ssh_ts  ***
580      !!---------------------------------------------------------------------- 
581      INTEGER, INTENT(in) ::   jn
582      !!
583      INTEGER :: ji, jj
584      INTEGER  :: istart, iend, jstart, jend
585      !!---------------------------------------------------------------------- 
586      !
587      IF( Agrif_Root() )   RETURN
588      !
589      ! --- West --- !
590      istart = 2
591      iend   = 1+nbghostcells
592      DO ji = mi0(istart), mi1(iend)
593         DO jj = 1, jpj
594            ssha_e(ji,jj) = hbdy(ji,jj)
595         ENDDO
596      ENDDO
597      !
598      ! --- East --- !
599      istart = jpiglo - nbghostcells
600      iend   = jpiglo - 1
601      DO ji = mi0(istart), mi1(iend)
602         DO jj = 1, jpj
603            ssha_e(ji,jj) = hbdy(ji,jj)
604         ENDDO
605      ENDDO
606      !
607      ! --- South --- !
608      jstart = 2
609      jend   = 1+nbghostcells
610      DO jj = mj0(jstart), mj1(jend)
611         DO ji = 1, jpi
612            ssha_e(ji,jj) = hbdy(ji,jj)
613         ENDDO
614      ENDDO
615      !
616      ! --- North --- !
617      jstart = jpjglo - nbghostcells
618      jend   = jpjglo - 1
619      DO jj = mj0(jstart), mj1(jend)
620         DO ji = 1, jpi
621            ssha_e(ji,jj) = hbdy(ji,jj)
622         ENDDO
623      ENDDO
624      !
625   END SUBROUTINE Agrif_ssh_ts
626
627   SUBROUTINE Agrif_avm
628      !!----------------------------------------------------------------------
629      !!                  ***  ROUTINE Agrif_avm  ***
630      !!---------------------------------------------------------------------- 
631      REAL(wp) ::   zalpha
632      !!---------------------------------------------------------------------- 
633      !
634      IF( Agrif_Root() )   RETURN
635      !
636      zalpha = 1._wp ! JC: proper time interpolation impossible 
637                     ! => use last available value from parent
638      !
639      Agrif_SpecialValue    = 0.e0
640      Agrif_UseSpecialValue = .TRUE.
641      !
642      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
643      !
644      Agrif_UseSpecialValue = .FALSE.
645      !
646   END SUBROUTINE Agrif_avm
647   
648
649   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
650      !!----------------------------------------------------------------------
651      !!                  *** ROUTINE interptsn ***
652      !!----------------------------------------------------------------------
653      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
654      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
655      LOGICAL                                     , INTENT(in   ) ::   before
656      !
657      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
658      INTEGER  ::   N_in, N_out
659      ! vertical interpolation:
660      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
661      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
662      REAL(wp), DIMENSION(k1:k2) :: h_in
663      REAL(wp), DIMENSION(1:jpk) :: h_out
664      !!----------------------------------------------------------------------
665
666      IF( before ) THEN         
667         DO jn = 1,jpts
668            DO jk=k1,k2
669               DO jj=j1,j2
670                 DO ji=i1,i2
671                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
672                 END DO
673              END DO
674           END DO
675        END DO
676
677# if defined key_vertical
678        DO jk=k1,k2
679           DO jj=j1,j2
680              DO ji=i1,i2
681                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
682              END DO
683           END DO
684        END DO
685# endif
686      ELSE 
687
688# if defined key_vertical             
689         DO jj=j1,j2
690            DO ji=i1,i2
691               N_in = 0
692               DO jk=k1,k2 !k2 = jpk of parent grid
693                  IF (ptab(ji,jj,jk,n2) == 0) EXIT
694                  N_in = N_in + 1
695                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
696                  h_in(N_in) = ptab(ji,jj,jk,n2)
697               END DO
698               N_out = 0
699               DO jk=1,jpk ! jpk of child grid
700                  IF (tmask(ji,jj,jk) == 0) EXIT
701                  N_out = N_out + 1
702                  h_out(jk) = e3t_a(ji,jj,jk)
703               ENDDO
704               IF (N_in > 0) THEN
705                  DO jn=1,jpts
706                     call reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),ptab_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
707                  ENDDO
708               ENDIF
709            ENDDO
710         ENDDO
711# else
712         ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)
713# endif
714         !
715         DO jn=1, jpts
716            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
717         END DO
718
719      ENDIF
720      !
721   END SUBROUTINE interptsn
722
723   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
724      !!----------------------------------------------------------------------
725      !!                  ***  ROUTINE interpsshn  ***
726      !!---------------------------------------------------------------------- 
727      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
728      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
729      LOGICAL                         , INTENT(in   ) ::   before
730      !
731      !!---------------------------------------------------------------------- 
732      !
733      IF( before) THEN
734         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
735      ELSE
736         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
737      ENDIF
738      !
739   END SUBROUTINE interpsshn
740
741   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
742      !!----------------------------------------------------------------------
743      !!                  *** ROUTINE interpun ***
744      !!---------------------------------------------   
745      !!
746      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
747      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
748      LOGICAL, INTENT(in) :: before
749      !!
750      INTEGER :: ji,jj,jk
751      REAL(wp) :: zrhoy
752      ! vertical interpolation:
753      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
754      REAL(wp), DIMENSION(1:jpk) :: h_out
755      INTEGER  :: N_in, N_out
756      REAL(wp) :: h_diff
757      !!---------------------------------------------   
758      !
759      IF (before) THEN
760         DO jk=1,jpk
761            DO jj=j1,j2
762               DO ji=i1,i2
763                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
764# if defined key_vertical
765                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
766# endif
767               END DO
768            END DO
769         END DO
770      ELSE
771         zrhoy = Agrif_rhoy()
772# if defined key_vertical
773! VERTICAL REFINEMENT BEGIN
774
775         DO ji=i1,i2
776            DO jj=j1,j2
777               N_in = 0
778               DO jk=k1,k2
779                  IF (ptab(ji,jj,jk,2) == 0) EXIT
780                  N_in = N_in + 1
781                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
782                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
783              ENDDO
784         
785              IF (N_in == 0) THEN
786                 ua(ji,jj,:) = 0._wp
787                 CYCLE
788              ENDIF
789         
790              N_out = 0
791              DO jk=1,jpk
792                 if (umask(ji,jj,jk) == 0) EXIT
793                 N_out = N_out + 1
794                 h_out(N_out) = e3u_a(ji,jj,jk)
795              ENDDO
796         
797              IF (N_out == 0) THEN
798                 ua(ji,jj,:) = 0._wp
799                 CYCLE
800              ENDIF
801         
802              IF (N_in * N_out > 0) THEN
803                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
804! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
805                 if (h_diff < -1.e4) then
806                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
807!                    stop
808                 endif
809              ENDIF
810              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
811            ENDDO
812         ENDDO
813
814# else
815         DO jk = 1, jpkm1
816            DO jj=j1,j2
817               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
818            END DO
819         END DO
820# endif
821
822      ENDIF
823      !
824   END SUBROUTINE interpun
825
826   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
827      !!----------------------------------------------------------------------
828      !!                  *** ROUTINE interpvn ***
829      !!----------------------------------------------------------------------
830      !
831      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
832      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
833      LOGICAL, INTENT(in) :: before
834      !
835      INTEGER :: ji,jj,jk
836      REAL(wp) :: zrhox
837      ! vertical interpolation:
838      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
839      REAL(wp), DIMENSION(1:jpk) :: h_out
840      INTEGER  :: N_in, N_out
841      REAL(wp) :: h_diff
842      !!---------------------------------------------   
843      !     
844      IF (before) THEN         
845         DO jk=k1,k2
846            DO jj=j1,j2
847               DO ji=i1,i2
848                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
849# if defined key_vertical
850                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
851# endif
852               END DO
853            END DO
854         END DO
855      ELSE       
856         zrhox = Agrif_rhox()
857# if defined key_vertical
858
859         DO jj=j1,j2
860            DO ji=i1,i2
861               N_in = 0
862               DO jk=k1,k2
863                  if (ptab(ji,jj,jk,2) == 0) EXIT
864                  N_in = N_in + 1
865                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
866                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
867               END DO
868               IF (N_in == 0) THEN
869                  va(ji,jj,:) = 0._wp
870                  CYCLE
871               ENDIF
872         
873               N_out = 0
874               DO jk=1,jpk
875                  if (vmask(ji,jj,jk) == 0) EXIT
876                  N_out = N_out + 1
877                  h_out(N_out) = e3v_a(ji,jj,jk)
878               END DO
879               IF (N_out == 0) THEN
880                 va(ji,jj,:) = 0._wp
881                 CYCLE
882               ENDIF
883               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
884            END DO
885         END DO
886# else
887         DO jk = 1, jpkm1
888            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
889         END DO
890# endif
891      ENDIF
892      !       
893   END SUBROUTINE interpvn
894
895   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
896      !!----------------------------------------------------------------------
897      !!                  ***  ROUTINE interpunb  ***
898      !!---------------------------------------------------------------------- 
899      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
900      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
901      LOGICAL                         , INTENT(in   ) ::   before
902      !
903      INTEGER  ::   ji, jj
904      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
905      !!---------------------------------------------------------------------- 
906      !
907      IF( before ) THEN
908         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
909      ELSE
910         zrhoy = Agrif_Rhoy()
911         zrhot = Agrif_rhot()
912         ! Time indexes bounds for integration
913         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
914         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
915         !
916         DO ji = i1, i2
917            DO jj = j1, j2
918               IF    ( utint_stage(ji,jj) == 1  ) THEN
919                  ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
920                     &               - zt0**2._wp * (       zt0 - 1._wp)        )
921               ELSEIF( utint_stage(ji,jj) == 2  ) THEN
922                  ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
923                     &               - zt0        * (       zt0 - 1._wp)**2._wp )
924               ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
925                  ztcoeff = 1._wp
926               ELSE
927                  ztcoeff = 0._wp
928               ENDIF
929               !   
930               ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
931               !           
932               IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
933                  ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
934                  utint_stage(ji,jj) = 3
935               ELSE
936                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
937               ENDIF
938            END DO
939         END DO
940      END IF
941      !
942   END SUBROUTINE interpunb
943
944
945   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
946      !!----------------------------------------------------------------------
947      !!                  ***  ROUTINE interpvnb  ***
948      !!---------------------------------------------------------------------- 
949      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
950      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
951      LOGICAL                         , INTENT(in   ) ::   before
952      !
953      INTEGER  ::   ji, jj
954      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
955      !!---------------------------------------------------------------------- 
956      !
957      IF( before ) THEN
958         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
959      ELSE
960         zrhox = Agrif_Rhox()
961         zrhot = Agrif_rhot()
962         ! Time indexes bounds for integration
963         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
964         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
965         !     
966         DO ji = i1, i2
967            DO jj = j1, j2
968               IF    ( vtint_stage(ji,jj) == 1  ) THEN
969                  ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
970                     &               - zt0**2._wp * (       zt0 - 1._wp)        )
971               ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
972                  ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
973                     &               - zt0        * (       zt0 - 1._wp)**2._wp )
974               ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
975                  ztcoeff = 1._wp
976               ELSE
977                  ztcoeff = 0._wp
978               ENDIF
979               !   
980               vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
981               !           
982               IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
983                  vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
984                  vtint_stage(ji,jj) = 3
985               ELSE
986                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
987               ENDIF
988            END DO
989         END DO         
990      ENDIF
991      !
992   END SUBROUTINE interpvnb
993
994
995   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
996      !!----------------------------------------------------------------------
997      !!                  ***  ROUTINE interpub2b  ***
998      !!---------------------------------------------------------------------- 
999      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1000      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1001      LOGICAL                         , INTENT(in   ) ::   before
1002      !
1003      INTEGER  ::   ji,jj
1004      REAL(wp) ::   zrhot, zt0, zt1, zat
1005      !!---------------------------------------------------------------------- 
1006      IF( before ) THEN
1007         IF ( ln_bt_fw ) THEN
1008            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1009         ELSE
1010            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1011         ENDIF
1012      ELSE
1013         zrhot = Agrif_rhot()
1014         ! Time indexes bounds for integration
1015         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1016         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1017         ! Polynomial interpolation coefficients:
1018         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1019            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1020         !
1021         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1022         !
1023         ! Update interpolation stage:
1024         utint_stage(i1:i2,j1:j2) = 1
1025      ENDIF
1026      !
1027   END SUBROUTINE interpub2b
1028   
1029
1030   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
1031      !!----------------------------------------------------------------------
1032      !!                  ***  ROUTINE interpvb2b  ***
1033      !!---------------------------------------------------------------------- 
1034      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1035      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1036      LOGICAL                         , INTENT(in   ) ::   before
1037      !
1038      INTEGER ::   ji,jj
1039      REAL(wp) ::   zrhot, zt0, zt1, zat
1040      !!---------------------------------------------------------------------- 
1041      !
1042      IF( before ) THEN
1043         IF ( ln_bt_fw ) THEN
1044            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1045         ELSE
1046            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1047         ENDIF
1048      ELSE     
1049         zrhot = Agrif_rhot()
1050         ! Time indexes bounds for integration
1051         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1052         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1053         ! Polynomial interpolation coefficients:
1054         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1055            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1056         !
1057         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1058         !
1059         ! update interpolation stage:
1060         vtint_stage(i1:i2,j1:j2) = 1
1061      ENDIF
1062      !     
1063   END SUBROUTINE interpvb2b
1064
1065
1066   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1067      !!----------------------------------------------------------------------
1068      !!                  ***  ROUTINE interpe3t  ***
1069      !!---------------------------------------------------------------------- 
1070      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1071      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1072      LOGICAL                              , INTENT(in   ) :: before
1073      INTEGER                              , INTENT(in   ) :: nb , ndir
1074      !
1075      INTEGER :: ji, jj, jk
1076      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1077      !!---------------------------------------------------------------------- 
1078      !   
1079      IF( before ) THEN
1080         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1081      ELSE
1082         western_side  = (nb == 1).AND.(ndir == 1)
1083         eastern_side  = (nb == 1).AND.(ndir == 2)
1084         southern_side = (nb == 2).AND.(ndir == 1)
1085         northern_side = (nb == 2).AND.(ndir == 2)
1086         !
1087         DO jk = k1, k2
1088            DO jj = j1, j2
1089               DO ji = i1, i2
1090                  !
1091                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1092                     IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN
1093                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1094                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk) 
1095                        kindic_agr = kindic_agr + 1
1096                     ELSEIF (eastern_side.AND.(ptab(i2-nbghostcells+1,jj,jk)>0._wp)) THEN
1097                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1098                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1099                        kindic_agr = kindic_agr + 1
1100                     ELSEIF (southern_side.AND.(ptab(ji,j1+nbghostcells-1,jk)>0._wp)) THEN
1101                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1102                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1103                        kindic_agr = kindic_agr + 1
1104                     ELSEIF (northern_side.AND.(ptab(ji,j2-nbghostcells+1,jk)>0._wp)) THEN
1105                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1106                        WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1107                        kindic_agr = kindic_agr + 1
1108                     ENDIF
1109                  ENDIF
1110               END DO
1111            END DO
1112         END DO
1113         !
1114      ENDIF
1115      !
1116   END SUBROUTINE interpe3t
1117
1118
1119   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1120      !!----------------------------------------------------------------------
1121      !!                  ***  ROUTINE interpumsk  ***
1122      !!---------------------------------------------------------------------- 
1123      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1124      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1125      LOGICAL                              , INTENT(in   ) ::   before
1126      INTEGER                              , INTENT(in   ) ::   nb , ndir
1127      !
1128      INTEGER ::   ji, jj, jk
1129      LOGICAL ::   western_side, eastern_side   
1130      !!---------------------------------------------------------------------- 
1131      !   
1132      IF( before ) THEN
1133         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1134      ELSE
1135         western_side = (nb == 1).AND.(ndir == 1)
1136         eastern_side = (nb == 1).AND.(ndir == 2)
1137         DO jk = k1, k2
1138            DO jj = j1, j2
1139               DO ji = i1, i2
1140                   ! Velocity mask at boundary edge points:
1141                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1142                     IF (western_side) THEN
1143                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1144                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1145                        kindic_agr = kindic_agr + 1
1146                     ELSEIF (eastern_side) THEN
1147                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1148                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1149                        kindic_agr = kindic_agr + 1
1150                     ENDIF
1151                  ENDIF
1152               END DO
1153            END DO
1154         END DO
1155         !
1156      ENDIF
1157      !
1158   END SUBROUTINE interpumsk
1159
1160
1161   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1162      !!----------------------------------------------------------------------
1163      !!                  ***  ROUTINE interpvmsk  ***
1164      !!---------------------------------------------------------------------- 
1165      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1166      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1167      LOGICAL                              , INTENT(in   ) ::   before
1168      INTEGER                              , INTENT(in   ) :: nb , ndir
1169      !
1170      INTEGER ::   ji, jj, jk
1171      LOGICAL ::   northern_side, southern_side     
1172      !!---------------------------------------------------------------------- 
1173      !   
1174      IF( before ) THEN
1175         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1176      ELSE
1177         southern_side = (nb == 2).AND.(ndir == 1)
1178         northern_side = (nb == 2).AND.(ndir == 2)
1179         DO jk = k1, k2
1180            DO jj = j1, j2
1181               DO ji = i1, i2
1182                   ! Velocity mask at boundary edge points:
1183                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1184                     IF (southern_side) THEN
1185                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1186                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1187                        kindic_agr = kindic_agr + 1
1188                     ELSEIF (northern_side) THEN
1189                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1190                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1191                        kindic_agr = kindic_agr + 1
1192                     ENDIF
1193                  ENDIF
1194               END DO
1195            END DO
1196         END DO
1197         !
1198      ENDIF
1199      !
1200   END SUBROUTINE interpvmsk
1201
1202
1203   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1204      !!----------------------------------------------------------------------
1205      !!                  ***  ROUTINE interavm  ***
1206      !!---------------------------------------------------------------------- 
1207      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1208      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1209      LOGICAL                                    , INTENT(in   ) ::   before
1210      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
1211      REAL(wp), DIMENSION(1:jpk) :: h_out
1212      INTEGER  :: N_in, N_out, ji, jj, jk
1213      !!---------------------------------------------------------------------- 
1214      !     
1215      IF (before) THEN         
1216         DO jk=k1,k2
1217            DO jj=j1,j2
1218              DO ji=i1,i2
1219                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1220              END DO
1221           END DO
1222        END DO
1223#ifdef key_vertical         
1224        DO jk=k1,k2
1225           DO jj=j1,j2
1226              DO ji=i1,i2
1227                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk) 
1228              END DO
1229           END DO
1230        END DO
1231#endif
1232      ELSE 
1233#ifdef key_vertical         
1234         avm_k(i1:i2,j1:j2,1:jpk) = 0.
1235         DO jj=j1,j2
1236            DO ji=i1,i2
1237               N_in = 0
1238               DO jk=k1,k2 !k2 = jpk of parent grid
1239                  IF (ptab(ji,jj,jk,2) == 0) EXIT
1240                  N_in = N_in + 1
1241                  tabin(jk) = ptab(ji,jj,jk,1)
1242                  h_in(N_in) = ptab(ji,jj,jk,2)
1243               END DO
1244               N_out = 0
1245               DO jk=1,jpk ! jpk of child grid
1246                  IF (wmask(ji,jj,jk) == 0) EXIT
1247                  N_out = N_out + 1
1248                  h_out(jk) = e3t_n(ji,jj,jk)
1249               ENDDO
1250               IF (N_in > 0) THEN
1251                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
1252               ENDIF
1253            ENDDO
1254         ENDDO
1255#else
1256         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1257#endif
1258      ENDIF
1259      !
1260   END SUBROUTINE interpavm
1261
1262#else
1263   !!----------------------------------------------------------------------
1264   !!   Empty module                                          no AGRIF zoom
1265   !!----------------------------------------------------------------------
1266CONTAINS
1267   SUBROUTINE Agrif_OCE_Interp_empty
1268      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1269   END SUBROUTINE Agrif_OCE_Interp_empty
1270#endif
1271
1272   !!======================================================================
1273END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.