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/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 10395

Last change on this file since 10395 was 8058, checked in by jgraham, 7 years ago

Clear keywords

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