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/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 6079

Last change on this file since 6079 was 6079, checked in by jamesharle, 8 years ago

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