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/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 9006

Last change on this file since 9006 was 9005, checked in by timgraham, 6 years ago

Add TOP update code

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